Compute the discrete Fourier transform of the following cosine signal: $x(t) = cos(2πft)$ , for 𝑓 = 50 Hz, sampled during 1 second at every ∆= 0.001 s.
The power spectrum does not look like a line spectrum, i.e. a 𝛿-function at 50 Hz. Why?
In [1]:
using PlotlyJS
In [5]:
# Define signal
f = 50
dt = 0.0001
Fs=1/dt; #sampling rate
fNyquist= Fs/2
t= collect(0:1*Fs-1)*dt
s = cos.(2*π*f*t);
In [20]:
"""
compute_periodogram(x; nfft=length(x), fs=1)
### Arguments:
x: signal
T_pad: time to pad
Note, this was done as an excercise. Could also use
DSP.periodogram for a more rocust function
"""
function compute_periodogram(x; T_pad=0, fs=1, T=1)
if T_pad > 0
T = T + 2*T_pad
pad_size = Int(T_pad*fs)
x = vcat(zeros(Float64, pad_size), x, zeros(Float64, pad_size))
end
N= T*Fs
N_half = Int(round(N/2+1))
# Periodogram (rectangular window): Discrete Fourier transform (via Fast Fourier Transform, FFT)
x = x - mean(x) #subtract mean: zero DC shift
xf = fft(x)
# One-sided power spectrum (periodogram)
Sp = dt^2 * 1/T * abs.(xf).^2 # dt^2 * 1/T * Xf.*conj(Xf);
Sp = 2 * Sp[1:N_half] #One-sided spectrum (i.e. only positive frequencies) for even N;
Sp_db = 10*log10.(Sp/maximum(Sp))
return Sp_db
end
Out[20]:
In [29]:
function plot_periodogram(f, Sp_db; pad=0)
trace = PlotlyJS.scatter(x = f, y=Sp_db)
layout = Layout(;title="Pad = $(pad)", showlegend=false,
xaxis=attr(title="Frequency (Hz)", zeroline=false, range=[40, 60]),
yaxis=attr(title="Power (dB)", zeroline=false, range=[-60, 0]))
p = plot([trace], layout)
end
Out[29]:
*Compute the discrete Fourier transform of the following cosine signal: $x(t) = cos(2πft)$ , for 𝑓 = 50 Hz, sampled during 1 second at every ∆= 0.001 s.
R: See computation above
*Plot the power spectrum using the periodogram (i.e. the one-sided power spectrum based on the DFT) in dB. Use three different zero padding conditions: 1, 2 and 6 second zero padding.
R: See plots below
In [32]:
pad = 0
df = 1/(pad*2 + 1); #frequency resolution
F = collect(0:df:fNyquist)#frequency grid
p_dB = compute_periodogram(s, T_pad = pad, fs = Fs)
p0 = plot_periodogram(F, p_dB)
pad = 1
df = 1/(pad*2 + 1); #frequency resolution
F = collect(0:df:fNyquist)#frequency grid
p1_dB = compute_periodogram(s, T_pad = pad, fs = Fs)
p1 = plot_periodogram(F, p1_dB, pad=pad)
pad = 2
df = 1/(pad*2 + 1); #frequency resolution
F = collect(0:df:fNyquist)#frequency grid
p2_dB = compute_periodogram(s, T_pad = pad, fs = Fs)
p2 = plot_periodogram(F, p2_dB, pad=pad)
pad = 6
df = 1/(pad*2 + 1); #frequency resolution
F = collect(0:df:fNyquist)#frequency grid
p6_dB = compute_periodogram(s, T_pad = pad, fs = Fs)
p6 = plot_periodogram(F, p6_dB, pad=pad)
[p0 p1 p2 p6]
Out[32]:
Consider the following signal x(t)= cos(2πf1t)+cos(2πf2t), for 𝑓1 = 57 Hz and 𝑓2 = 60 Hz.
What is the maximum frequency resolution for a Fourier transform applied to a 200 ms time window?
R: fs = 1/0.2 = 5Hz.
Would you be able to discriminate between the two frequencies?
R: No since the f2-f1 < fs
Would zero padding help the discrimination? What should one do in order to discriminate the two frequencies?
R: No, padding does not effectively increase sampling rate. We could increase the duration of the window
In [ ]: