Hello all,
I’ve been working with the FFTW library over the past few days and am struggling to understand what exactly I have when I take the fft of some array. In particular, taking f to be some function defined over [-\pi, \pi], I want to compute the coefficients c_n such that:
f(x) \approx \sum_{n=-N+1}^N c_n \exp(i n x)
From my understanding of the FFTW library and what it returns (based on what the documentation says), if c is the array returned by a fft of said function over a grid of size N, then I should be able to say
f(x) \approx \frac{1}{N} \sum_{n=1}^N c[n] \exp(i (n-1) x)
For example, take this small code snippet.
using FFTW
using LaTeXStrings
using Plots
N = 101;
x = LinRange(-pi, pi, N);
sinx = sin.(x);
sinx_h = fft(sinx);
f(z) = real( sum(sinx_h .* exp.(im*(0:N-1)*z)) / N );
plot(xlabel=L"x \in [-\pi,\pi]", legend=:topleft)
plot!(x, sinx, lw=3, label="sin(x)")
plot!(x, f.(x), lw=3, ls=:dash, label="Fourier Representation");
savefig("fourier_test.png")
I’d expect both sinx and f.(x) to be roughly equal, but they’re flipped and off by about a factor of 2. Any ideas as to why this is happening / where my misunderstand of either the math or what FFTW is doing is?
Hm, In addition tested this with the fourier coefficients computed from a simple for loop as detailed in Numerical Methods by Dahlquist and Bjorck chapter 9.3.1. Something like:
c = complex.(zeros(size(x)));
for j=1:N
for i=1:N
c[j] += sinx[i]*exp(-im*2*pi*(i-1)*(j-1)/N);
end
end
and I see the same result as the (incorrect) Fourier Representation plot above. Perhaps this indicates that the problem isn’t in me misunderstanding the FFTW library, but rather misunderstanding the math here?
Ah, indeed that works. Indeed, I do (hazily) remember that the fft was sensitive to even/odd grid sizes and endpoints over the period but will need to look up exactly why. Also thank you for the sign convention correction.
@Abhijit_Chowdhary, sorry but I’ve deleted my previous post as the answer was not general. Meanwhile, I am linking this other post where the Fourier coefficients were written more carefully by following their definition.
The source of confusion might be that with FFT we need to reorder the positive and negative frequency bins with fftshift().
Please check the code below where complex coefficients are used.
using FFTW, LaTeXStrings, Plots; gr()
function FourierSeries(x, fx, T)
ck = 1/N * fftshift(fft(fx))
n = 1 + length(x)÷2
yr = zero(fx) .+ 0.0im
for (i, ci) in pairs(ck)
yr .+= ci * exp.(im*2π*(i-n)/T * x)
end
return real(fftshift(yr)), ck
end
N = 100 # use N even
x = LinRange(-π, π, N+1)[1:N]
T = 2π
fx = @. sin(2π*x/T) + cos(2π*4*x/T)
fsx, ck = FourierSeries(x, fx, T)
p1 = plot(xlabel=L"x \in [-\pi,\pi]", legend=:topleft)
plot!(x, fx, lw=3, label="f(x)")
plot!(x, fsx, lw=1, lc=:black, ls=:dash, label="Fourier Series")
xp = 2π*(-N÷2:N÷2-1)/T
p2 = plot(xp, real.(ck), st=:stem, lc=:blue, label="Real(ck)")
plot!(xp, imag.(ck), st=:stem, lc=:red, xlabel="k", label="Imag(ck)")
plot(p1,p2)
1 Like
