Dear Ferret users,
The functions
ffta and
fftp return the wrong frequency axis when the length of the data is an odd number. The error is minor as long as the length of data is large. See below.
The source of the error comes from this statement:
"the frequency axis runs from freq1 = yquist/ float(nfreq) to freqn = yquist"
from "Ch3 Sec2.3.26. FFTA" of the official manual.
What's wrong with the statement?
The fundamental formula given in the manual
f(t) = S(j=1 to N/2)[a(j) cos(jwt + F(j))]
where w=2 pi/T
is correct. According to this formula, the angular frequencies are
w, 2w, 3w, . . .
In terms of frequency,
(1/T), 2(1/T), 3(1/T), . . .
[Note that angular frequency = 2 pi frequency.] So, the FFT function should return this frequency axis.
On the other hand, yquist = 1/(2 Δt) and Δt = T/N. When N is even, nfreq = N/2. So,
yquist/nfreq = 1/T
When N is even, therefore, yquist/nfreq is the correct fundamental frequency.
But when N is odd, nfreq = (N−1)/2. So,
yquist/nfreq = [N/(N−1)](1/T)
So, when N is odd, yquist/nfreq is the wrong fundamental frequency. Regardless of whether N is odd or not, the fundamental frequency is 1/T.
I paste a script that illustrates the problem, comparing the frequency axis returned by ffta with the correct frequencies.
Best regards,
Ryo
==================
let ndaysyear = 365 ! or 360
!let ndaysyear = 360 ! or 365
define axis/t=0:`ndaysyear`/edges/npoints=`ndaysyear`\
/calendear=`ndaysyear`_day/modulo=`ndaysyear`/units=days mytax
let tvals = t[gt=mytax]
let pi = 4 * atan(1.0)
let fr1 = 1.0 / ndaysyear ! fundamental freq in cycles/day
let w1 = 2*pi*fr1 ! fundamenal freq in rad/day
let mode1 = cos(w1 * tvals)
let mode2 = cos(2*w1 * tvals)
let ft = mode1 + mode2
let amp = ffta(ft)
list/L=1:6/format=(G0,1X,G0) amp, t[gt=amp]
list/format=(G0) fr1
list/format=(G0) 2*fr1