Using fft to plot frequency spectrum of sum of rectangular pulses

42 views (last 30 days)
fs = 10000;
t = -2:1/fs:2;
g_a = rectpuls(t/4) + rectpuls(t/2);
fn = fs/2;
f = linspace(-fn, fn, numel(t));
G_a = fft(g_a); % Fourier transform of function g_a
G_a_mag = abs(G_a);
plot(t,g_a) % plot function wrt time
plot(f, G_a_mag) % plot freq spectrum
axis([0 2 0 2])
Here is my code, I would like to obtain the correct plot for the frequency spectrum, but what I obtain is weird and wrong. By computing the result analytically, the Fourier transform should be the following : 2*sinc(omega) + 4*sinc(2*omega).
I can't figure out what I have done wrong, maybe I am not using the fft function correctly ?
Thank you for your help.

Accepted Answer

Paul on 6 Mar 2022
The code below illustrates the analysis for one component of the signal. It can be adapted to include both components of the signal.
Note the two key differences between this code and the code in the question, both of which are noted below.
First, compute the exact CTFT
syms t w f
g(t) = rectangularPulse(-0.5,0.5,t/2);
G(w) = simplify(fourier(g(t),t,w))
G(w) = 
which is the expected result stated in the question.
Now compute the samples of the signal g(t) at the speciifed sampling period over the specified time duration.
fs = 10000; Ts = 1/fs;
tval = -2:1/fs:2; % might be bettter to use linspace
gval = rectpuls(tval/2);
Note that rectpuls() results in gval = 1 at tval = -1, but gval = 0 at tval = 1 (which I believe is a desired result).
[tval(10001) gval(10001)]
ans = 1×2
-1 1
[tval(30001) gval(30001)]
ans = 1×2
1 0
Consequently, gval has an even number of non-zero samples for parameters of the question
ans = 20000
This fact will be important next.
Now compute the DTFT of gval using freqz() over the range -30 to 30 rad/sec. freqz() assumes that the first point in gval corresponds to the discrete-time sample n = 0, so we have to shift the output of freqz() usnig the DTFT time-shifting property to account for the fact that the first element of gval actually corresponds to n = -20000 (-20000/fs = -2)
wdtft = linspace(-1,1,1001)*30; % -30 rad/s to 30 rad/sec
dtft = freqz(gval,1,wdtft*Ts).*exp(1j*wdtft*Ts*20000);
dtft has a small, non-zero imaginary part because nnz(gval) is even, i.e., gval is not an odd sequence.
ans = 1.0000
To compare to G(w), we'll just look at the real part.
dtft = real(dtft);
Now compute the DFT of gval, the elements of which are samples of the DTFT. Shift the DFT in the frequency domain for the same reason we shifted the DTFT.
dft = fft(gval).*exp(1j*(0:(numel(gval)-1))/numel(gval)*2*pi*20000);
wdft = (0:(numel(dft)-1))/numel(dft)*2*pi*fs;
Get rid of the small imaginary part.
ans = 1.0000
dft = real(dft);
Now compare the CTFT and DTFT over the frequency range -30 to 30 rad/sec, and the DFT over 0-30 rad/sec. Could use fftshift() if we wanted to get the image of the DFT onto negative frequencies. Note that the code in the question did not use fftsthift() even though f was defined over from -fn to fn. Also, the DTFT and DFT have to be scaled by the sample period to approximate the CTFT; this scaling also needs to be added ot the code in the question.
hold on
fplot(G(w),[-30 30])
Note that starting with the third point of the DFT every other point is zero, which would result in a triangular look if we plotted the abs() of the three curves and connected the samples of the DFT.

More Answers (1)

Scott MacKenzie
Scott MacKenzie on 5 Mar 2022
Edited: Scott MacKenzie on 6 Mar 2022
I think this is what you are after. Mostly, this code just follows the example in the fft documentation, particularly w.r.t. getting the single-sided frequency spectrum from the FFT. The top wavefrom is a square wave of period 4 seconds, or frequency 0.25 Hz. The bottom waveform is the frequency spectrum. The peaks are the harmonics of the square wave. These occur at 0.25 Hz, 0.75 Hz, 1.25 Hz, and 1.75 Hz, as expected. I'm not sure why the waveform is triangular. It might be because the fft function is only working with one cycle of the square wave.
fs = 10000;
t = -2:1/fs:2;
g_a = rectpuls(t/4) + rectpuls(t/2);
g_a(end) = 1; % fix last sample (was 0)
g_a = rescale(g_a, -1, 1); % remove dc offset in signal
G_a = fft(g_a); % Fourier transform of function g_a
% get frequency spectrum from G_a (see fft documenation)
n = numel(t)-1;
P2 = abs(G_a/n);
G_a_mag = P2(1:n/2+1);
G_a_mag(2:end-1) = 2*G_a_mag(2:end-1);
f = fs*(0:(n/2))/n;
plot(t,g_a) % plot function wrt time
axis([-2 2 -1.5 1.5])
plot(f, G_a_mag) % plot freq spectrum
axis([0 2 0 2])

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!