6
votes

I'm trying to visualize a signal and its frequency spectrum in Julia.

I found the FFTW package that provides the FFT and DSP for the frequencies.

Here is what I'm trying, with a sinusoidal signal:

using Plots
using FFTW
using DSP

# Number of points
N = 2^14 - 1
# Sample rate
fs = 1 / (1.1 * N)
# Start time
t0 = 0
tmax = t0 + N * fs

# time coordinate
t = [t0:fs:tmax;]

# signal
signal = sin.(2π * 60 * t)  # sin (2π f t)

# Fourier Transform of it
F = fft(signal)
freqs = fftfreq(length(t), fs)
freqs = fftshift(freqs)

# plots
time_domain = plot(t, signal, title = "Signal")
freq_domain = plot(freqs, abs.(F), title = "Spectrum")
plot(time_domain, freq_domain, layout = 2)
savefig("Wave.pdf")

I expected to see a nice plot with a peak in the 60 Hz, but all I got was a weird result:

Imgur

I'm ignoring the negative frequencies for now.

How should I do that in Julia?

1

1 Answers

10
votes

What you call fs in your code is not your sampling rate but the inverse of it: the sampling period.

The function fftfreq takes the sampling rate as its second argument. Since what you give as the second argument is the sampling period, the frequencies returned by the function are incorrectly scaled by (1/(Ts^2)).

I renamed fs to Ts and changed the second argument to fftfreq to the sampling rate 1.0/Ts. I think you also need to shift the result of fft.

# Number of points 
N = 2^14 - 1 
# Sample period
Ts = 1 / (1.1 * N) 
# Start time 
t0 = 0 
tmax = t0 + N * Ts
# time coordinate
t = t0:Ts:tmax

# signal 
signal = sin.(2π * 60 .* t) # sin (2π f t) 

# Fourier Transform of it 
F = fft(signal) |> fftshift
freqs = fftfreq(length(t), 1.0/Ts) |> fftshift

# plots 
time_domain = plot(t, signal, title = "Signal")
freq_domain = plot(freqs, abs.(F), title = "Spectrum", xlim=(-1000, +1000)) 
plot(time_domain, freq_domain, layout = 2)
savefig("Wave.pdf")