Extracting Fourier Coefficients of an arbitrary periodic signal?

New to Julia
1

Hello!

Does a package exist in Julia which takes a signal and return the Fourier Coefficients needed to approximate this signal?

Kind regards

2

If by a signal you mean a vector of floating point numbers, then FFTW.jl will transform it using fft for you.

4 Likes
3

May I profits from this question to get better understanding of the fft. Following the wikipedia formula for fourier series (Fourier series - Wikipedia) I though that this simple code will reconstruct a real signal from its harmonics:

using FFTW, Plots
fs = 256
Δt = 1 / fs
P = 1.0
x = 0.0:Δt:P
y = [sin(2 * pi * 7 * t) + sin(2 * pi * 15 * t) + sin(2 * pi * 30 * t) for t in x] # mixture of simple wave signal
plot(x, y, legend = false)
Fy = rfft(y)
freq = 0.0:(fs / length(y)):(fs / 2)
N = length(Fy)
ak = real.(Fy ./ N)
bk = imag.(Fy ./ N)
yr = [ak[1] / 2 for _ in x]
for i in 2:N
    yr .+= ak[i] .* cos.(2 * pi * freq[i] .* x) .+ bk[i] .* sin.(2 * pi *freq[i] .* x)
end
plot!(x, yr)

But actually this give me a signal that is flip from the x axis and shift in time by Δt, to get the correct reconstruction I found that I have to do:

yr = [ak[1] / 2 for _ in x]
for i in 2:N
    yr .+= ak[i] .* cos.(2 * pi * freq[i] .* (x .+ Δt)) .+ bk[i] .* sin.(2 * pi *freq[i] .* (x .+ Δt))
end
yr .= -yr

May I ask why ?

4

You might need fftshift or ifftshift

1 Like
5

As I use rfft I already get only half the spectrum and from my understanding fftshift only perform a circular shit to put the 0 frequency in the middle of the spectrum. For example if I want the reconstruction using fft instead I will have to do this:

plot(x, y, legend = false)
Fy = fftshift(fft(y))
freq = 0.0:(fs / length(y)):(fs / 2)
freq = [freq[end:-1:2]; freq]
N = length(Fy)
ak = real.(Fy ./ N)
bk = imag.(Fy ./ N)
yr = [ak[div(N, 2) + 1] / 2 for _ in x]
for i in 1:div(N, 2) # negative frequency
    yr .+= ak[i] .* cos.(2 * pi * freq[i] .* (x .+ Δt)) .- bk[i] .* sin.(2 * pi * freq[i] .* (x .+ Δt))
end
for i in (div(N, 2) + 2):N # positive frequency
    yr .+= ak[i] .* cos.(2 * pi * freq[i] .* (x .+ Δt)) .+ bk[i] .* sin.(2 * pi * freq[i] .* (x .+ Δt))
end
yr .= -yr
plot!(x, yr)

But once again I have to add a time shift of Δt and to reverse the signal from the x axis… It seems I’m missing something there

6

I want the specific coefficients, i.e. a0, a1,b1,a2,b2 etc.

When I use that library and the function fft I get the Fourier transformed signal, but how do I get it to output the coefficients it used?

Kind regards

7

@Ahmed_Salih the ak and bk coefficient from my code should be what your are looking for… just need to figure out why the reconstruction need to be inverted and shifted !

1 Like
8

@Geoffrey, I have edited your code below to be closer to the definition of the Fourier coefficients, used the standard fft and adjusted the length of input signal x to contain N samples. The minus sign required in the imaginary part of the fft, is probably due to the fft sign convention in Julia.
Edit: defined a1 and had loop starting at 1

using FFTW, Plots
N = 256;
P = 1.0;
Δt = P / N
x = 0.0:Δt:(P-Δt)   # lenght(x) == N
y = [sin(2π*7*t) + sin(2π*15*t) + sin(2π*30*t) for t in x] # mixture of simple wave signal
plot(x, y, legend = false, linewidth=2)
Fy = fft(y)[1:N÷2]
ak =  2/N * real.(Fy)
bk = -2/N * imag.(Fy)  # fft sign convention
ak[1] = ak[1]/2
yr = zeros(N,1)
for i in 1:N÷2
    yr .+= ak[i] * cos.(2π*(i-1)/P * x) .+ bk[i] * sin.(2π*(i-1)/P * x)
end
plot!(x, yr, linestyle=:dash, linewidth=2)

fourier_coefficients

6 Likes
9

@rafael.guerra thank you very much for your answer, indeed my mistake was with the sign of the bk coefficients !

1 Like