Incorrect Frequency axis returned from FFTA and other FFT functions
Closed this issue · 1 comments
ACManke commented
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 commented
This is fixed and the changes are committed to my fork. As Ryo says, the differences are minor.