# FFT of a Square wave, what are frequency of each Harmonic here

After getting suggestion to first study FFT on square wave to understand FFT of discrete signal, I'm trying to understand the FFT of square wave.

I've coded a program, here is the details, Frequency of Wave, 100KHz

f= 100e3;


Sampling Frequency, 1MHz

Fs = 1e6;
Ts = 1/Fs;   %Sampling Rate
% Ts = 0.000001

t = 1/f;       %Time period of 1 Oscillation = 1/f

% t = 0.00001
n = 0.000001:Ts:t;  %Generating Samples

% t/Ts = 10


i.e 10 Samples will be generated

% 1st  Sample at 0.000001s
% 2nd  Sample at 0.000002s
% 3rd  Sample at 0.000003s
% 4th  Sample at 0.000004s
% 5th  Sample at 0.000005s
% 6th  Sample at 0.000006s
% 7th  Sample at 0.000007s
% 8th  Sample at 0.000008s
% 9th  Sample at 0.000009s
% 10th Sample at 0.00001s


defining and plotting the square wave, Since Amplitude is not defined it will be 1 by default.

x=square(2*pi*f*n);
subplot(2,1 ,1)
stem(n,x)


Output wave,

Now plotting its FFT

subplot(2,1 ,2)
F=fft(x)
plot(2*abs(F)./(t/Ts));


Output wave,

For the First point,

Fs / # of Samples

1000,000/10 = 100KHz <- First point = 100KHz

Number of Harmonics should be = 10 but What are the Frequency of each Harmonic ? also, why the amplitude has changed to 1.2 which is supposed to be 1

-
Shouldn't the sample at time 10 be -1, not +1? – Jim Clay Oct 21 '12 at 19:52
Since this question has been resurrected, I will echo @JimClay's note that the signal that you showed above is not a square wave. That can help to explain why you aren't seeing the frequency-domain characteristics that you expect. – Jason R Feb 5 '13 at 18:26

From the Wikipedia article on the Discrete Fourier Transform:

The sequence of $N$ complex numbers $x_0, ..., x_{N−1}$ is transformed into an $N$-periodic sequence of complex numbers according to the DFT formula: $$X_k=\sum_{n=0}^{N-1} x_n e^{-2\pi ikn/N}.$$

What you have done is taken $N=10$ integers, $n=\{1,2,...,N\}$, and turned them into $10$ time samples $t=\{T_s,2T_s,...,NT_s\}$, or $t_n=nT_s$. (To clarify my notation change, this second set is what you called $n$.)

Then you created the signal $x_n=\text{square}(2\pi ft_n)$ with period $f=(NT_s)^{-1}$and fed it into MATLAB's FFT algorithm, attempting to take the DFT of it (as the FFT is just a fast DFT). However, this is NOT what the FFT expects to see!

Your signal starts with $x_1$ and ends with $x_N$. However, if you see the definition of the DFT I gave above, it expects the signal to start with $x_0$ and end with $x_{N-1}$. Ordinarily, this would not be much of a problem, because if your signal $x$ is actually $N$-periodic, then $x_0=x_N$. However, as noted in the comments by Jim Clay and Jason R above, the signal you start with is not actually a square wave. As you can see in your screenshot, there are six "1" values and only four "-1" values. The square wave should have an equal number of "1"s and "-1"s. I do not know why the values you put in are not a proper square wave, and I suspect there is some odd detail in how MATLAB implimented the function $\text{square}$. To create a square wave, you should change the line

n = 0.000001:Ts:t;  %Generating Samples


to

n = 0:Ts:t-Ts;  %Generating Samples


or, even better, to

N=t/Ts;
n=(0:N-1)*Ts;


which makes it clear that you are sampling at integer multiples of your sampling time. The signal $x$ you generate in this way will be equivalent to what Jim Clay has generated in his answer.

As to why your signal has magnitude $\approx 1.2$ instead of $1$, you need to remember how the square wave is defined. From the Wikipedia article on the square wave:

$$x_{\mathrm{square}}(t) =\frac{4}{\pi}\left (\sin(2\pi ft) + {1\over3}\sin(6\pi ft) + {1\over5}\sin(10\pi ft) + \cdots\right ).$$

The first term of this function has frequency $f$ and magnitude $\frac{4}{\pi}\approx 1.27$. If you look at Jim Clay's plot this is exactly the magnitude in bin 2 of the function he has plotted. Up to bin $N/2+1$, the value that will be plotted in bin $k$ is the coefficient of the term in the square wave with frequency $(k-1)f$. (The $-1$ comes from MATLAB indexing beginning with one rather than zero).

-

When I do everything exactly as you outline in your post, except change x to the following-

x = [ones(1,5), -ones(1,5)];


I get the following plot-

FFT bin 2 is the first harmonic and FFT bin 4 is the third harmonic (square waves only have odd harmonics). Likewise, FFT bin 10 is the negative first harmonic, and FFT bins 8 and 6 are the negative third and fifth harmonics, respectively. You can see that the harmonics gradually diminish in power.

-