next up previous
Next: Plotting Capabilities Up: Data Analysis Previous: Difference Equations & Filters


Fourier Transforms

Fourier transforms are used in numerous scientific applications. Let me give you a bit of background on what it is all about.

For any one-dimensional function (usually in time) f(x), the fourier transform is defined as:

\begin{displaymath}F(u) = \int_{-\infty}^{\infty} f(x) \exp{(-t2{\pi}ux)} dx \end{displaymath}

And the inverse fourier transform (IFT) is defined as:

\begin{displaymath}f(x) = \int_{-\infty}^{\infty} F(u) \exp{(t2{\pi}ux)} du \end{displaymath}

Note that the above shown transforms are based on a continuous function f(x). Hence they are not computationally feasible. Fast Fourier transform (FFT) is just a modified discretized version of the Fourier transform algorithm (DFT) which reduces the number of computations needed for N points from $2N^2$ to $2N \lg{N}$.

DFT for a sequence `x' of length `N' goes like this (`X' is a vector of length `N' again)...


\begin{displaymath}X(k) = \sum_{n=1}^N x(n) \exp{(-j2{\pi}(k-1)\left(\frac{n-1}{N}\right)), 1\le k \le N} \end{displaymath}

And the inverse discrete fourier transform (IDFT) is:

\begin{displaymath}x(n) = \frac{1}{N} \sum_{k=1}^N X(k) \exp{(j2{\pi}(k-1)\left(\frac{n-1}{N}\right)), 1\le n \le N} \end{displaymath}

Fourier transforms have a whole bunch of useful properties, which is the reason they are so widely used in the scientific community. For such information refer [5] and [6] as well as other references you might find!

Alright now with the theory of Fourier transforms out of the way, let me show how you can do DFT in Matlab. As usual a discussion follows the source code.

% Hint: use ; at end of line to prevent debug'ish output from being printed

x = [4 3 7 -9 1 0 0 0]' 
y = fft(x)               % Compute Discrete FFT
ifft(y)                  % Inverse

pause;                   % Pause before proceeding to next example

t = 0:1/100:10-1/100;
x = sin(2*pi*15*t) + sin(2*pi*40*t);

y = fft(x); 
m = abs(y);              % get magnitude and phase angle (next line)
p = unwrap(angle(y));    % Unwrap - removes multiples of Pi

f = (0:length(y)-1)'*100/length(y);
subplot(2,1,1), plot(f,m), 
ylabel('Abs. Magnitude'), grid on
subplot(2,1,2), plot(f,p*180/pi)
ylabel('Phase [Degrees]'), grid on
xlabel('Frequency [Hertz]')

I've shown two examples in the above code. The first one defines a simple vector `x' and performs FFT on it. I've also shown the computation of inverse transform using the ifft() function.

The second function defines a time span on the variable `t'. Then `x' is defined as sinuisoidal function of `t'. Then the fft() function is used.

Also note I've used a few other functions: abs(), angle() and unwrap(). abs() essentially gives the magnitude of the result while angle() gives the phase angle. unwrap() removes any redundant multiples of $\pi$ in the phase angle.

Following these computations is a plotting of the data. I will come back to it, if possible, after looking into plotting capabilities in Matlab.

x =
     4
     3
     7
    -9
     1
     0
     0
     0

y =
   6.0000          
  11.4853 - 2.7574i
  -2.0000 -12.0000i
  -5.4853 +11.2426i
  18.0000          
  -5.4853 -11.2426i
  -2.0000 +12.0000i
  11.4853 + 2.7574i

ans =
    4.0000
    3.0000
    7.0000
   -9.0000
    1.0000
   -0.0000
         0
    0.0000
Figure 7: Plot of Magnitude and Phase after doing a FFT
\begin{figure}\epsfig{file=src/fft1.eps, width=13cm}
\end{figure}


next up previous
Next: Plotting Capabilities Up: Data Analysis Previous: Difference Equations & Filters
Arvind Gopu 2006-03-24