Fast Fourier Transforms

Fast Fourier Transform Algorithms

10/23/08

1

By Akshaya Srivatsa and Ajith J Kocheri

Overview  Definition

and need for FFT  Principle behind the algorithms  Direct computation of FFT  Radix 2 FFT algorithm  Goertzel algorithm  Chirp Z transform  FFT in Matlab and Applications of FFT 10/23/08

2

By Akshaya Srivatsa and Ajith J Kocheri

Definition Given a sequence of N complex valued numbers x(n), the DFT X(k) is defined as

And the inverse DFT is given by the equation

10/23/08

3

By Akshaya Srivatsa and Ajith J Kocheri

Need for FFT 

10/23/08

The direct computation of FFT requires N2 complex multiplications and N(N-1) complex additions, which has polynomial complexity.

4

By Akshaya Srivatsa and Ajith J Kocheri

Some Examples Speech can be sampled “Adequately” with a sampling rate of 8kHz (N=8000).  Digital audio sampling rate is about 44.1 kHz (N=44100). This sampling rate is employed in Digital Audio CD’s  Video sampling rate is much higher than audio sampling rates. 

10/23/08

5

By Akshaya Srivatsa and Ajith J Kocheri

Principle behind the Algorithms 

10/23/08

Essentially all these algorithms use symmetry and periodicity property to their advantages to make sure that smaller number of steps are used to computer the DFT. Where WN= e-j2π/N

6

By Akshaya Srivatsa and Ajith J Kocheri

Radix 2 Approach 

V[k] = Σn even WNkn v[n] + Σ n odd WNkn v[n]



V[k] = Σr=0,1,…N-1 WNkr v[2r] + Σr=0,1,…N-1WNk(2r+1) v[2r+1]



V[k] = Σr=0,1,…N-1 WNkr v[2r] + WNkΣr=0,1,…N-1WNk(2r) v[2r+1]





V[k] =F1(k) + WNk F2(k) where F1(k) and F2(k) are N/2 point DFT’s for even and odd number sequences. But WNk+(N/2) = -WNk



Hence V[k] =F1(k) + WNk F2(k)

and

V[k+(N/2)] =F1(k) - WNk F2(k) k=0,1,2,…(N/2)-1 

10/23/08

We can further “Decimate” the even and odd sequences to give smaller and smaller divisions. 7

By Akshaya Srivatsa and Ajith J Kocheri

Butterfly Calculations Three stages in computation when N=8

V[0] = W20 v[0] + W20 v[1] = v[0] + W20 v[1] V[1] = W20 v[0] + W21 v[1] = v[0] + W21 v[1] Note that the V(k) is not in the same sequence as v(n). It can be observed that v(k) is in bit reversed order. 10/23/08

8

By Akshaya Srivatsa and Ajith J Kocheri

  

Decimation in Time

This is a special case of Divide and Conquer approach. Here N is represented be the equation N=2γ So we need to decimate X(n) γ= log2N times. This will give us decimation where each butterfly computation involves only 2 inputs at a time. -Total number of complex multiplications---(N/2)log2N -Total number of complex additions----------Nlog2N Here we see that we need only 2N computer registers as these 2N registers which initially stored x(n) will finally hold X(k). Hence we say that the Computation is in place.

10/23/08

9

By Akshaya Srivatsa and Ajith J Kocheri

Decimation in Frequency  

This is the reverse process of Decimation in time. Here we decimate X(k) successively to get x(n). Here we choose M=2 and L=N/2 Unlike decimation in time, here we split x(n) as first N/2 number and the next N/2 numbers.

10/23/08

10 By Akshaya Srivatsa and Ajith J Kocheri

How good is the Algorithm??

10/23/08

11 By Akshaya Srivatsa and Ajith J Kocheri

Efficient Computing DFT of Two Real sequences 



Fact: DFT computation is the same whether we handle real numbers of complex numbers. If the two real sequences are x1(n) and x2(n) then let x(n) = x1(n) + jx2(n)



So x1(n) = [x(n) + x*(n)]/2 and x2(n) =[x(n) – x*(n)]/2j



Also DFT of x*(n) = X*(N-k) Therefore X1(k) = [X(k) + X*(N-k)]/2 and



X1(k) = [X(k) - X*(N-k)]/2j

10/23/08

12 By Akshaya Srivatsa and Ajith J Kocheri

Goertzel Algorithm 



This algorithm uses linear filtering operation where the linear filter takes the form of N parallel resonators where each resonators selects one of the N frequencies in the DFT. We know that WN-kN=1, we can rewrite X(k) as X(k)=Σ0N-1 x(m) WN–k(N-m)



This is of the form of circular convolution where hk(n)=WN-knu(n)



Lets take the z transform of hk(n), we get Hk(z) = 1/[1-WN-kz-1] which has a single pole at ω = 2πk/N



10/23/08

From this, we get yk(n)=wN-kyk(n-1) + x(n). This makes evaluation of yk(n) an iterative process. Hence for values of M
13 By Akshaya Srivatsa and Ajith J Kocheri

Goertzel Algorithm (Continued) 

For the system having two pole filters, Hk(z) is of the form Hk(z) = [1-WN-k]/[1-2cos(2πk/N)z-1 +z-2] Given yk(-1) and yk(-2)

10/23/08

14 By Akshaya Srivatsa and Ajith J Kocheri

Chirp-z Transform Approach 



Z transform of x(n) is same as the DFT of x(n)r-n as DFT is obtained by setting r=1 in the Z transform of x(n) More generally speaking, in the z-plane, if we chose zk=roejθ(RoejØ)k where the first complex number is a point in the z plane and the second complex number is a circle or radius Ro, and if Ø=2π/N,θ=0,ro=0, we obtain the frequencies as points on the z plane.

nk=[1/2][n2+k2-(n-k)2]

10/23/08

15 By Akshaya Srivatsa and Ajith J Kocheri



So we have transformed Z transform into convolution which can be easily obtained using convolution theorem where 2



h(n)=ejπ/n(k-n) For the convolution, we got go consider an M point sequence of h(n) where M=L+N-1



Let n vary such that –(N-1)


Compute the M point DFT of h1(n) and M point DFT of an(After padding with zeros). Therefore Y1(k)=An(k)H1(k)



Take IDFT of Y1(k); We know that the first N-1 points are corrupted by aliasing. So discard the first N-1 point and consider the rest.

10/23/08

16 By Akshaya Srivatsa and Ajith J Kocheri

FFT in Matlab Y = fft(X) It returns FFT of the matrix X given  Y = fft(X,n) It returns FFT of X taking only n points 

10/23/08

17 By Akshaya Srivatsa and Ajith J Kocheri

Applications

x(n)=sin(2*50πt)+sin(2*150πt) + random normalized noise

10/23/08

18 By Akshaya Srivatsa and Ajith J Kocheri

 

10/23/08

Calculation of DTFS = [fft(x(n),N)]/N Frequency response of LTI filter is given by DTFT of its impulse response. This can be obtained by zero padding x(n) upto ∞ (Very large) and the DFT of it will give us the DTFT.

19 By Akshaya Srivatsa and Ajith J Kocheri