2
$\begingroup$

I am trying to understand FFT, DFT through Matlab and I am fairly inexperienced in all of these. I have a time vector t and a corresponding simple sine wave x = 0.25*sin(t). Now I would like to view the 1-sided FFT of this

%MATLAB 
%define time and signal
delta_t = 2.6;
t = 0:delta_t:600;
x=0.25*sin(t);
L=length(x);

%define frequency vector
f_low  =1/t(end); %lowest frequency for this range
fs = 1/delta_t; %sampling frequency
f=[1:L]*f_low;

%calculate FFT
FFT = fft(x);

%double sided spectrum
P1=abs(FFT)/L;

%single sided spectrum
P2=P1(1:L/2+1);
P2(2:end-1)=2*P2(2:end-1);
f2=f(1:L/2+1);

%plot
figure()
plot(f2,P2,'LineWidth',1.25);
grid on
xlabel('f')
ylabel('abs(FFT)')

The result of this varies with delta_t (all the while respecting Nyquist theorem). The amplitude of the FFT is different for differnt delta_t and ideally shouldn't it be 0.25, equal to the amplitude of the wave?

Can someone please explain me why does the amplitude vary with delta_t?

Thank you! FFT Amplitude for deltat = 0.9FFT Amplitude for deltat = 2.6

$\endgroup$
  • $\begingroup$ the dft has bin centers. your sine frequency is off the center frequency $\endgroup$ – Stanley Pawlukiewicz Jul 10 at 10:55
  • $\begingroup$ Thanks! I looked about it and found about padding and windowing. $\endgroup$ – Chandramouli Santhanam Jul 10 at 11:47
  • $\begingroup$ those are important topics but all you really need for your experiment is to have an exact number of whole cycles of your sine wave. the last sample needs to continue the sine wave to the first sample. $\endgroup$ – Stanley Pawlukiewicz Jul 10 at 12:03
2
$\begingroup$

Your code has some troubles, and perhaps also your theoretical understanding. Let me put here a mini summary of observing a practical sine wave and its frequency spectrum using FFT function of MATLAB / Octave.

Assume that there's a continuous-time ideal infinitely long sinusoidal wave with frequency $\Omega_0$ given as:

$$ x(t) = A \cos(\Omega_0 t) , \tag{1} $$

Let's sample this infinetely long signal within a finite duration of $[0 , T_s]$ using a sampling rate of $F_s$ or the sampling period of $T_s = \frac{1}{F_s}$, so that $t = nT_s$, and we obtain $L$ samples of the signal :

$$ x_L[n] = x[n]w[n] = A \cos(\omega_0 n) ~~,~~n = 0,1,...,L-1 \tag{2}$$

where $\omega_0 = \Omega_0 T_s$ and $w[n]$ is a rectangular window, used to represent the finite duration of observation.

$$ w[n] = \begin{cases} { 1 ~~~,~~~~ 0 \leq n \leq L-1 \\ 0 ~~~,~~~ ow }\end{cases} \tag{3}$$

DTFT (Discrete-time Fourier transform) modulation property can be used to analyse the frequency domain representation of $x_L[n]$ :

$$ X_L(e^{j\omega}) = \frac{1}{2\pi} X(e^{j\omega}) \star W(e^{j\omega}) \tag{4} $$

where $W(e^{j\omega})$ is the DTFT of the rectangular window :

$$ W(e^{j\omega}) = e^{-j \omega (\frac{L-1}{2}) } \cdot \frac{ \sin(\omega L /2)}{\sin(\omega/2)} \tag{5} $$

A plot of this function is shown below :

enter image description here

Again we also we know that the DTFT of $x[n] = A\cos(\omega_0 n)$ is : $$ X(e^{j\omega}) = A\pi \delta(\omega + \omega_0) + A\pi \delta(\omega - \omega_0) \tag{6} $$

Hence we see that the DTFT of the finite length observation $x_L[n]$ is

$$ X_L(e^{j\omega}) = \frac{A}{2} [ W(e^{j(\omega+\omega_0)}) + W(e^{j(\omega-\omega_0)})] \tag{7} $$

Now, when you take $M$-point FFT of the $L$-point sequence $x_L[n]$, then you are effectively calculating the $M$-point DFT (discrete Fourier transform) which is equal to $M$ uniform samples of DTFT $X_L(e^{j\omega})$ as given by

$$ X_L^M[k] = X_L(e^{j\omega})|_{w = \frac{2\pi}{M} k} ~~~, ~~~k = 0,1,...,M-1 \tag{8} $$

which yields the complex frequency vector that you compute using FFT function:

$$ X_L^M[k] = \frac{A}{2} [ W(e^{j(\omega_k+\omega_0)}) + W(e^{j(\omega_k-\omega_0)})] ~~,~~k = 0,...,M-1 \tag{9} $$

So with an M-point FFT of $x_L[n]$ we calculate the samples of the shifted versions of $W(e^{j\omega})$.

An example plot of $|X_L^M[k]|$ for $L = 16$ , $M = 1024$, $A = 5$ and $\omega_0 = 0.2 \pi$ is shown.

enter image description here

Note that the peak values (which are $A \cdot L /2$ when $\omega_0$ meets some criterion, see below) are obtained for some $k=k_0$ and $k = M-k_0$. Consider the first peak (to the right of the plot here) and it can be shown that its value is

$$ |X_L^M[k_0]| = |0.5 A [ W( \frac{2\pi}{M} k_0 - w_0) + W( \frac{2\pi}{M} k_0 + w_0) ]| \tag{10} $$

Where the integer DFT index $k_0$ is found to be $k_0 = \text{round}( \frac{M \omega_0}{2\pi} )$. So the peak value of the spectrum is in general not equal to $0.5 A \cdot L$ which you are expecting, which would only be the case when

$$ \frac{2\pi}{M} k_0 = w_0 \tag{11}$$ for some integer $k_0$ which requires that $\frac{M \omega_0}{2\pi}$ is also an integer, or equivalently

$$ \omega_0 = \frac{ 2\pi }{M}k_0 \tag{12}$$ ,in that special case $|X_L^M[k_0]| = 0.5 A L$ will happen. Otherwise, you will have two complex numbers added and their magnitude is computed.

Note that by selecting $M = L$ when $\omega_0$ meets the criterion (i.e, $ \frac{M \omega_0}{2\pi}$ is an integer), then one can get what is an illusion of lollipop plot indicating only two nonzero samples at the frequencies of $\pm \omega_0$ with magnitude 0.5 $A \cdot L $ and all zeros. But this is a conseqeunce of frequency sampling and not to be taken unnecessarily into too far conclusions.

Here is a MATLAB / OCTAVE code to produce the results

clc; clear all; close all

% Define the sampling grid  
Fs = 1000;             % sampling frequency in Hz
Ts = 1/Fs;             % sampling period in seconds

ti = 0;                % initial time (s)
tf = .15;               % final time   (s)
t =  ti :Ts : tf;  % sampling time index t = n*Ts

L = length(t);          % number of samples in x[n] = A*sin(2*pi*f0*n*Ts)
M = L;                  % Number of DFT points to be computed

% Define the sampled signal
A = 0.5;               % amplitude of sine wave

% Signal Frequency Type-1 :
%f0 = 100;              % sine wave frequency in Hz.
%w0 = 2*pi*f0*Ts;        % set an arbitrary frequency
%x = A*cos(w0*[0:L-1]);   % Sine wave sampled at Fs.

% Signal Frequency Type-2:
w0 = 2*pi/M*floor(L/3); % set a special frequency so that DFT PEAK = A/2            
x = A*cos(w0*[0:L-1]);   % Sine wave sampled at Fs.


% Calculate DFT of x[n] using FFT function
k0 = ( M*w0/(2*pi));


figure,stem(linspace(-pi,pi,M),fftshift(abs(fft(x,M))/L))
title(['Magnitude of X_L^M[k] for L=',num2str(L),' , M=', num2str(M),...
      ' , A=',num2str(A),' , \omega_0 =', num2str(w0), ' , k_0 = ' num2str(k0)]);
grid on
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.