FFT Secrets Revealed

QUESTION: Can you explain the IDL FFT function to me?

ANSWER: The following FFT tutorial was first published on the IDL newsgroup by Peter Brooker (ra5589@email.sps.mot.com). It is published here with his permission.

Introduction

Below is a repeat of my orginal post to the group in which I tried to explain FFT use from the point of view of somebody who has little experience of FFT algorithms but does understand the concept of Fourier series as well as Fourier transforms. This update includes a large modification to part 1 that develops the complex expansion from the series expansion. In it I show how all the imaginary components of the complex FFT expression cancel to zero. Modifications were also made to part 2.

This reference is organized as follows:

PART 1: Relate complex expansion to real Fourier series
PART 2: IDL form of complex expansion
PART 3: Specific example: f(t) = sin ( 4*pi*t)
PART 4: Specific example: T(x,y) =sin(6pi*x) + cos(4pi*y)

PART 1: Relate complex expansion to real Fourier series.

Assume you have a function f(t) that is periodic in t with a period T. Then there exists coefficients a_n and b_n such that:

   f(t) = a_o + sum_n(a_n*cos(2*pi*n*t/T)+b_n*sin(2*pi*n*t/T)) , n=1,2,...

The coefficients are given by:

   a_o = 2/T*integral(f(t)*dt)
   a_n = 2/T*integral( f(t)*cos(2*pi*n*t/T)*dt )
   b_n = 2/T*integral( f(t)*sin(2*pi*n*t/T)*dt )

The limits of the integration are from -T/2 to T/2.

This is just the Fourier series of function with period T. Nothing new here. See eq #4, section 10.3, Advanced Engineering Mathematics, Kreyszig. This expansion though assumes -T/2 < t < T/2.

Let's now define the complex coefficients:

   c_o = a_o
   c_n = 1/2*(a_n + j*b_n)

Then for n = 1,2,...

   c_n =  1/2*(a_n + j*b_n)
       = 2/T*1/2*(integral( f(t)*cos(2*pi*n*t/T)*dt ) +
                            j*2/T*integral( f(t)*sin(2*pi*n*t/T)*dt ))
       = 1/T*integral( f(t)*[cos(2*pi*n*t/T) + j*sin(2*pi*n*t/T)]*dt)
       = 1/T*integral( f(t)*exp(j*2*pi*n*t/T)*dt )

Since exp(0)=1 we can write:

   c_n = 1/T*integral( f(t)*exp(j*2*pi*n*t/T)*dt ) , n= 0,1,2,...

Note also that:

  [(c_n)^*] = 1/T*integral( f(t)*exp(-j*2*pi*n*t/T)*dt )
            = 1/T*integral( f(t)*exp(j*2*pi*(-n)*t/T)*dt )
            = c_(-n)

Here [(c_n)^*] denotes the complex conjugate of the element inside ( ). (Sorry for the clumsy notoation!!!)

Now here is something completely different ...

Given the definition of c_n above:

   f(t) = sum_n( c_n*exp(-j*2*pi*n*t/T) ), n = ...-2,-1,0,1,2,...

Proof

Define AA and BB by:

   sum_n( c_n*exp(-j*2*pi*n*t/T) ) = AA + c_o + BB,

where A is:

   A  = sum_n( c_n*exp(-j*2*pi*n*t/T) ), n = ...,-2,-1
   BB = sum_n( c_n*exp(-j*2*pi*n*t/T) ), n = 1,2,....

Let's work on AA some more:

   AA = sum_n( c_n*exp(-j*2*pi*n*t/T) ), n = ...,-2,-1
      = sum_n( c_(-n)* exp(-j*2*(-n)*t/T)), n = 1,2,....

Now:

   c_(-n) = [(c_n)^*]

and:

   exp(-j*2*(-n)*t/T) = [( exp(-j*2*pi*n*t/T))^*]

therefore:

   c_(-n)* exp(-j*2*(-n)*t/T) = [(c_n)^*] * [( exp(-j*2*pi*n*t/T) )^*]
                              = [(c_n* exp(-j*2*pi*n*t/T))^*]

Using this we have:

   AA  = sum_n( [c_n*exp(-j*2*n*t/T)]^*), n = 1,2,...
       = [(   sum_n( c_n*exp(-j*2*n*t/T) )   )]^*

Comparing this expression to sum BB we see that:

  AA = [(BB)^*]

Let's write the original sum again with this information:

   sum_n( c_n*exp(-j*2*pi*n*t/T) ) = AA + c_o + BB = [(BB)^*] + c_o +BB
                                   = c_o + 2*Re(BB)

Now:

   BB= sum_n( c_n*exp(-j*u) ), n = 1,2,....

and:

   u = 2*pi*n*t/T
     = sum_n( 1/2*(a_n+j*b_n)*(cos(u) - j*sin(u) )
     = 1/2*sum_n( a_n*cos(u) + j*(-j)*b_n*sin(u) +j*(b_n*cos(u)-a_n*sin(u) )
     = 1/2*sum_n(a_n*cos(u) + b_n*sin(u) +j*( b_n*cos(u)-a_n*sin(u) )

From this we see that:

   2*Re(BB) = 2*1/2*sum_n(a_n*cos(u) + b_n*sin(u))
            = sum_n(a_n*cos(u) + b_n*sin(u))

Now we have:

   sum_n( c_n*exp(-j*2*pi*n*t/T) ) = c_o + 2*Re(BB)
       = a_o + sum_n(a_n*cos(u) + b_n*sin(u)), n = 1,2,...

This is exactly the definition for f(t). I have therefore proved that:

   f(t) = sum_n( c_n*exp(-j*2*pi*n*t/T) ) ,  n = ...,-2,-1,0,1,2,...

Part 2: IDL form of complex expansion

Let f(t) be a periodic function with period T defined on an interval [0,T]. Then there exist complex A_n such that:

   f(t)= sum_n(A_n*exp(j*2*pi*n*t/T)), n=...,-1,0,1,...   j*j = -1

Divide the interval into N sections. t~t_i = i*T/N Then:

   f( t_i ) = sum_n(A_n*exp(j*2*pi*n*t_i/T))
            = sum_n(A_n*exp(j*2*pi*n*i*(T/N)/T))
            = sum_n( A_n*exp(j*2*pi*n*i/N) ) , n=...,-1,0,1,...

What about the ns?

IDL truncates the sum from -N/2+1, -N/2 + 2, ...., -1,0,1,...,N/2.

This is how IDL performs the FFT. This is found in the IDL manual under the section for FFT. The only difference is that t has been replaced by u and A_n has been replaced by F(u). Note that the period T has dropped out. Also note that t has been replaced by t_i = i*T/N. In order for this to happen, the interval over which t is defined must be from [0,T]. This is different from the definition of t being defined over the interval [-T/2,T/2].

It gets more complicated. From the manual we have:

   F(u) = 1/N*sum_x(f(x)*exp(-j*2pi*ux/N)) , x=0,1,...N-1

First thing to realize is that F(u) is really F_n. Where n is an integer. This comes from the fact that f(x) is periodic in x. The manual also mentions that the "frequencies" are:

   Fo, 1/(NT), 2/(NT), ..., 1/2T,-(N-2)/(2NT), ..., -1/NT

After trial and error I have determined that the value of the ns range for -N/2 to N/2. Futhermore, the F_n are stored in the order associated with the following values of n:

   0,1,2,...,N/2,-(N/2-1),-(N/2-1),...,-1  (Note, this is bizarre!!)

Let N=8. Then N/2=4. The F_n would be stored in an array. The array of n values associated with this array would be:

   [0,1,2,3,4,-3,-2,-1]

Part 3: Specific Example

Consider the interval t = [0,1]. This choice of interval implies T=1. Let f(t) = sin ( 4*pi*t). Then:

   f(t_i)=sin(2pi*2*i/N), i=0,1,...N

   f(t_i)= sum_n(A_n*exp(-j*2pi*n*i/N)) , n=-N/2,...-1,0,1,...N/2
         = A_nN/2... + A_n2*(cos(2pi*(-2)*i/N)+j*sin(2pi*(-2)*i/N))+
                           + A_o+A_n1*exp()+A_1*exp()
                           + A_2*(cos(2pi*(2)*i/N)+j*sin(2pi*(2)*i/N))
                           + A_3*exp()+...
        = ... + A_n2*cos(2pi*2*i/N)+A_2*cos(2pi*2*i/N) +
              + A_n2*j*(-1)*sin(2pi*2*i/N)+A_2*j*sin(2pi*2*i/N)) + ....
        = ... + (A_n2+A_2)*cos(2pi*2*i/N)+j*( -A_n2 + A_2)*sin(2pi*2*i/N)

   + ...

where A_n2 stands for A_n where n= -2.

Equating the series to sin(2pi*2*i/n) we conclude A_n = 0 for all n except n = -2 or n = 2.

   A_n2+A_2=0
   j*(-A_n2 + A_2) = 1

Let A_n2=(a_n2+j*b_n2) and A_2=(a_2 + j*b_2). The above equations imply:

   (a_n2 + a_2)   + j*(  b_n2+b_2) = 0

and

   j*[( -a_n2 + a_2) + j*( -b_n2 + b_2)] = 1

which implies:

   a_n2 + a_2=0, b_n2+b_2 =0
   a_n2 = - a_n2, b_n2 = -b_n2
   2*a_n2=0
   a_n2 = a_2 = 0
   j*j*(2*b_2) = 1 
   2*b_2 = -1/2, b_2 = 1/2
   A_n2 = 0 + j*(1/2)
   A_2 = 0 + j*(-1/2)

We now have calculated the solutions.

The following code calculates this and displays the correct answers. It shows how to plot A_n vs n correctly.

   ;idl_program fft_sine.pro
   TT=1
   Npts=100
   t=findgen(Npts)/(Npts-1)*TT
   f_t=sin(4.*!pi*t/TT)
   !p.multi=[0,2,3]
   plot,t,f_t, title='f(t) vs t'
   A_n=fft(f_t,-1) ; complex fourier coefficients

   plot,float(A_n),yrange=[-.5,.5],title='float(A_n)'
   plot,imaginary(A_n),yrange=[-.5,.5],title='imaginary(A_n)'

   a=findgen(Npts/2+1)
   b=-reverse(findgen(Npts/2-1)+1)
   c=[a,b] ; c=[-N/2+1,-N/2+2, ...,-1,0,1,...,N/2]
   print,c
   sub=sort(c)
   plot,c(sub),float(A_n(sub)),yrange=[-.5,.5],title='float(A_n) vs n'
   plot,c(sub),imaginary(A_n(sub)),yrange=[-.5,.5], $
           title='imaginary(A_n) vs n'
   plot,c(sub),imaginary(A_n(sub)),xrange=[-5,5],$
           title='imaginary(An) vs n' ; finer x scale
   end

Part 4: Specific Example:2D

   T(x,y)= sin(6pi*x) + cos(4pi*y)

Hopefully the program below is somewhat self explantory. It calculates the C=FFT(T,-1).

   ;idl_program fft_2D.pro
   !p.multi=0
   Tx=1
   Nx=100
   x=findgen(Nx)/(Nx-1)*Tx
   Ty=1
   Ny=100
   y=findgen(Ny)/(Ny-1)*Ty
   T=fltarr(Nx,Ny)
   for i=0,Nx-1 do begin
      for j=0,Ny-1 do begin
         T(i,j)= sin(2.*!pi*3.*i/Nx) + cos(2.*!pi*2*j/Ny)
         endfor
      endfor
   ;
   !p.multi=[0,2,2]
   shade_surf,T,x,y,xtitle='x',ytitle='y',title='T(x,y)'
   ;
   ;
   C=fft(T,-1) ; complex fourier coefficients
   surface,float(C)
   aaa=where(float(c) gt .4)
   surface,imaginary(C)
   ;
   a=findgen(Nx/2+1)
   b=-reverse(findgen(Nx/2-1)+1)
   ns=[a,b] ; this is the array of n's associated with C(n). n goes with x
   ;print,ns ; n goes from 0,...,Nx/2, -(Nx/2-1),...,-1
   subn=sort(ns) ;  n goes with x
   n_sort=ns(subn)
   ;
   a=findgen(Ny/2+1)
   b=-reverse(findgen(Ny/2-1)+1)
   ms=[a,b] ; this is the array of m's associated with C(n,m)
   print,ms ; m goes from 0,...,Ny/2, -(Ny/2-1),...,-1
   subm=sort(ms) ; m goes with y
   m_sort=ms(subm)
   ;
   sub_n_p3=where(ns eq 3)
   sub_n_n3=where(ns eq -3)
   sub_n_0=where(ns eq 0)
   ;
   sub_m_p2=where(ms eq 2)
   sub_m_n2=where(ms eq -2)
   sub_m_0=where(ms eq 0)
   ;
   print,'C(3,0),c(-3,0)=',C(sub_n_p3,sub_m_0),C(sub_n_n3,sub_m_0)
   ;
   print,'C(0,2),c(0,-2)=',C(sub_n_0,sub_m_p2),C(sub_n_0,sub_m_n2)
   ;
   ; now we need to define CC(n,m) to have normal scaling in n & m.
   ;
   CC=C*0.
   ;
   for n=0,Nx-1 do begin
      for m=0,Ny-1 do begin
         CC(subn(n),subm(m))= C(n,m)
         endfor
      endfor
   ;
   surface,ABS(CC),n_sort,m_sort
   ;
   end