Frequency spectrum of a square wave in MATLAB

Question

I need to plot the frequency spectrum for a square wave using MATLAB. The wave is HIGH (5mV) between 0 and -2 and LOW (omv) between 0 and 2. I have already obtained the fourier seires for this function and i have the first ten components of the series.

(5/2) + ((10/pi)*sin((pi*t)/2)) + ((10/(3*pi))*sin((3*pi*t)/2)) + ((10/(5*pi))*sin((5*pi*t)/2)) + ((10/(7*pi))*sin((7*pi*t)/2))+ ((10/(9*pi))*sin((9*pi*t)/2))+ ((10/(11*pi))*sin((11*pi*t)/2))+ ((10/(13*pi))*sin((13*pi*t)/2))+ ((10/(15*pi))*sin((15*pi*t)/2))+ ((10/(17*pi))*sin((17*pi*t)/2))+ ((10/(19*pi))*sin((19*pi*t)/2))

How do I plot the frequency spectrum for this wave using MATLAB? I have tried using FFTs, but I really don't know how it works to plot the graph.

Solution

This is adapted straight from the fft docs but using a square wave as the signal:

Fs = 1000;                    %// Sampling frequency
freq = 50;
T = 1/Fs;                     %// Sample time
L = 1000;                     %// Length of signal
t = (0:L-1)*T;                %// Time vector
A = 10;                       %// Amplitude

y = A*(sin(2*pi*freq*t) > 0); %// Make a square wave


plot(Fs*t,y);
xlabel('time (milliseconds)')

NFFT = 2^nextpow2(L);         %// Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

%// Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1))) 
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)') %// probably would be more meaningful if you convert this to radians per second, check out the xtickmark and xtickmarklabel properties of plot
ylabel('|Y(f)|')

Comments

[user2802349]

I tried this out with some variations to fit my square wave, but the graph isn't exactly what it's supposed to be. I get 11 peaks as required, but shouldn't i get a (5/2) peak at a frequency of 0 (from the fourier series)? I get a peak at a little over 5. code Fs = 1000; T = 1/Fs; L = length(a); t2 = (0:L-1)*T; NFFT = 2^nextpow2(L); Y = fft(a,NFFT)/L; f = Fs/2*linspace(0,1,NFFT/2+1); plot(f,2*abs(Y(1:NFFT/2+1))) axis([0, 5 0, 5]) code Where a is the waveform I posted in the original question

[Dan]

I don't understand... you have to create t before you create a, and you need L before t. You can't have a different t to calculate a to the one used for the fft plots...

[Dan]

@user2802349 hmmm... I'm not sure why there is a 2* in the fft plot, but that would explain why you get twice the amplitude for the 0Hz signal. I would investigate that if I were you