Fourier Transform and Inverse Fourier Transform for Signal Processing in MATLAB

 

Contents
  • Software for this post
  • Introduction to FFT and IFFT
  • FFT of Time Series Signal using MATLAB
  • IFFT of Spectrum to Time Series Signal using MATLAB
  • Conclusion
Software for this post
  • MATLAB/Simulink (verified by version R2015b);
Introduction to FFT and IFFT

It is known that any waveform (signal) can be rewritten as the sum of sinusoidal functions. Fourier Transform is the tool for the realization. By analyzing a signal in the frequency domain, we may know the characteristics of a system.

Let’s start from Fourier Series written as

(1)    \begin{equation*} f(t) = a_0+\sum_{m=1}^{\infty} a_m \cos\left(\frac{2\pi m t}{T}\right) + \sum_{m=1}^{\infty} b_m \sin\left(\frac{2\pi m t}{T}\right), \end{equation*}

where T is the fundamental period, a_m and b_m are the coefficients of the Fourier Series. The coefficients are given by

(2)    \begin{equation*} a_0=\frac{1}{T} \int_0^T f(t) dt , \end{equation*}

(3)    \begin{equation*} a_m = \frac{2}{T} \int_0^T f(t) \cos \left(\frac{2\pi m t}{T}\right) dt , \end{equation*}

(4)    \begin{equation*} b_m = \frac{2}{T} \int_0^T f(t) \sin \left(\frac{2\pi m t}{T}\right) dt. \end{equation*}

Equation (1) can be written using complex exponential function that

(5)    \begin{equation*} f(t) = \sum_{n=-\infty}^{\infty} c_n exp\left(i\frac{2\pi n t}{T}\right), \end{equation*}

where

(6)    \begin{equation*} c_n=\frac{1}{T} \int_0^T f(t) exp\left(i\frac{2\pi n t}{T}\right)dt. \end{equation*}

Fourier transform and inverse Fourier transform are defined as

(7)    \begin{equation*} \mathcal{F}\{f(t)\} = F(f) = \int_{-\infty}^{\infty}f(t) exp(-2\pi i f t) dt , \end{equation*}

(8)    \begin{equation*} \mathcal{F}^{-1}\{F(f)\} = \int_{-\infty}^{\infty}F(f) exp(2\pi i f t) df = f(t). \end{equation*}

Note that Laplace transform has a similar form as Fourier transform that

(9)    \begin{equation*} \mathcal{L}\{f(t)\} = F(s) = \int_{t=0}^{\infty}f(t)exp(-st) dt, \end{equation*}

where s is a complex number frequency parameter s=\sigma+i\omega.

FFT of Time Series Signal using MATLAB

MATLAB function for FFT of time series signal is

function fft_freq(SAMPL,H)
  N = length(SAMPL);
  if(mod(N,2) ~= 0)
    SAMPL = SAMPL(1:N-1);
    N = N-1;
  end
  f = (0:N/2);
  f = f/N*1/H;   
  sampl_fft2 = fft(SAMPL,N);%
  sampl_fft1 = sampl_fft2(1:N/2+1);% nyquist frequency
  sampl_fft1(2:end-1) = abs(sampl_fft1(2:end-1)/N) * 2;

  fonts = 20;
  figure;
  h=semilogx(f,sampl_fft1,'-');
  grid;
  ylabel('Amplitude','FontSize',fonts);
  set(gca,'Fontsize',fonts);
  set(h,'Linewidth',2.5)

  xlabel('Frequency [Hz]','FontSize',fonts);
  title('FFT','FontSize',fonts);
end

The function can be downloaded from here. fft_freq
Here is an example of using the function fft_freq.

t=0:0.02:50;
y=sin(2*pi*t)+cos(2*pi*2.5*t);
fft_freq(y,0.02);

The following results can be obtained.

IFFT of Spectrum to Time Series Signal using MATLAB

An IFFT MATLAB function transforming spectrum to time series signal is given in the following.

function [y,t]=Freq2Time(w,H)
%This function calculate the time signal using ifft function for the
%signal's spectrum

dw=w(2)-w(1);
 if(w(1)>0)
     add_w=0:dw:w(1)-dw; 
     NH=length(add_w);
     [m,n]=size(w);
     [p,q]=size(H);
     if(m>1)
         w=[add_w';w];
         if(p>1)
             H=[zeros(NH,1);H];
         else
             H=[zeros(1,NH) H];
         end
     else
         w=[add_w w];
         
         if(p>1)
             H=[H(1)*ones(NH,1);H];
         else
             H=[H(1)*ones(1,NH) H]';
         end
     end
        
 end


 Ts=(2*pi/max(w));
 N=length(w);
 t=0:Ts:(N-1)*Ts;
 
% compute ifft
 IR_data=real((ifft(H, 'symmetric' )/Ts)); 
 y=IR_data(1:N);
 if(y(1)>0) 
    y(1)=2*y(1);
 end
end

The function can be downloaded from here. Freq2Time
Here, we also give a comparison of the result using the above function with the result using impulse function in MATLAB.

sys=tf([1,2,300],[1,2,3,1]);
w=0:0.01:62.83;
N=length(w);
[mag,phase]=bode(sys,w);

H=mag(:).*exp(1j*phase(:)*pi/180);

% setup time step for ifft
Ts=(2*pi/max(w));
t=0:Ts:(N-1)*Ts;
y2 = Freq2Time(w,H);
[y1,t1]=impulse(sys,t);
figure(1);
plot(t,y2,'r',t1,y1,'b--');
xlim([0 10]);

And the results are shown in the following figure:

Conclusion

Fourier Transform and Inverse Fourier Transform are powerful mathematical tools for signal and system analysis. For the advanced application, please read some more about Nyquist frequency and window functions.

Leave a Reply

Your email address will not be published. Required fields are marked *