Incorrect Frequency axis returned from FFTA and other FFT functions

[ACManke]

Ryo Furue reports this issue here,
https://www.pmel.noaa.gov/maillists/tmap/ferret_users/fu_2019/msg01191.html @ryofurue

The frequency axis defined by the function FFTA is half the length of the time axis, the number rounded down if NT is odd. But the fundamental frequency should be based on the full time-span, whether there is an even number or an odd number of timesteps.

A simple example:

yes? set list/prec=8

! Define a time axis with 100 timesteps

yes? let nt = 100
yes? define axis/t=1:`nt`:1/units=minutes tax
 !-> define axis/t=1:100:1/units=minutes tax
 
! Define a function of time on this axis

yes? let tvals = t[gt=tax]
 
yes? let pi = 4 * atan(1.0)
yes? let fr1 = 1.0 / nt ! fundamental freq in cycles/minute
yes? let w1  = 2*pi*fr1 ! fundamental freq in rad/minute
yes? let mode1 = cos(w1 * tvals)
yes? let mode2 = cos(4*w1 * tvals)
yes? let ft = mode1 + mode2

yes? let amp = ffta(ft)
 
! The first coordinate on the frequency axis returned by FFTA 
! should equal fr1, the fundamental cycle.  It does.
 
yes? list fr1
             VARIABLE : 1.0 / NT
          0.0100000000
yes? list 2*fr1
             VARIABLE : 2*FR1
          0.020000000
 
yes? list/L=1:6 amp, t[gt=amp]
             T (CYC/MINUTES): 0.005 to 0.065
 Column  1: AMP is FFTA(FT)
 Column  2: T is T (axis (AX002))
                 AMP        T
0.01  / 1:  1.0000000  0.010000000
0.02  / 2:  0.0000000  0.020000000
0.03  / 3:  0.0000000  0.030000000
0.04  / 4:  1.0000000  0.040000000
0.05  / 5:  0.0000000  0.050000000
0.06  / 6:  0.0000000  0.060000000
 


! Redefine the time axis with 101 timesteps. The frequency axis
! from FFTA should be based on the fundamental frequency of 1/101
 
yes? let nt = 101
yes? define axis/t=1:`nt`:1/units=minutes tax
 !-> define axis/t=1:101:1/units=minutes tax
Replacing definition of axis TAX
 
! The frequency axis coordinates are incorrect. They are the same 
! as those of the 100-point axis.

yes? list fr1
             VARIABLE : 1.0 / NT
          0.0099009901
yes? list 2*fr1
             VARIABLE : 2*FR1
          0.019801980
 
yes? list/L=1:6 amp, t[gt=amp]
             T (CYC/MINUTES): 0.005 to 0.065
 Column  1: AMP is FFTA(FT)
 Column  2: T is T (axis (AX001))
                  AMP        T
0.01  / 1:  1.00000000  0.010000000
0.02  / 2:  0.00000000  0.020000000
0.03  / 3:  0.00000000  0.030000000
0.04  / 4:  1.00000000  0.040000000
0.05  / 5:  0.00000000  0.050000000
0.06  / 6:  0.00000000  0.060000000

The routine ef_set_freq_axis should be removed, and the FFT* functions that call it changed to instead correctly call ef_set_custom_axis.

[ACManke]

This is fixed and the changes are committed to my fork. As Ryo says, the differences are minor.