Using fft to plot frequency spectrum of sum of rectangular pulses

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); 

subplot(2,1,1)
plot(t,g_a) % plot function wrt time
subplot(2,1,2)
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.

Answer 1

[Paul]

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

nnz(gval)

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.

max(abs(imag(dtft)))

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.

max(abs(imag(dft)))

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.

figure
hold on
fplot(G(w),[-30 30])
plot(wdtft,dtft*Ts)
plot(wdft(wdft<30),Ts*dft(wdft<30),'ko'),grid
xlabel('rad/sec')
legend('CTFT','T*DTFT','T*DFT')

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.

[Japhten]

Thank you for your answer, it was very useful !

Answer 2

[Scott MacKenzie]

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;

subplot(2,1,1)
plot(t,g_a) % plot function wrt time
axis([-2 2 -1.5 1.5])
subplot(2,1,2)
plot(f, G_a_mag) % plot freq spectrum
axis([0 2 0 2])

[Japhten]

Thanks a lot for your answer !