Assignment 3 Fast Fourier Transform

1、Fast Fourier Transform1 IntroductionA fast Fourier transform (FFT) algorithm computes the discrete Fourier transform (DFT) of a sequence, or its inversion (IFFT). Fourier analysis converts a signal from its original domain (often time or space) to a representation in the frequency domain and vice ve

2、rsa. A FFT rapidly computes such transformations by factorizing the DFT matrix into a product of sparse (mostly zero) factors.1 As a result, it manages to reduce the complexity of computing the DFT from O(n2), which arises if one simply applies the definition of DFT, to O(n log n), where n is the da

3、ta size.The development of fast algorithms for DFT can be traced to Gauss's unpublished work in 1805 when he needed it to interpolate the orbit of asteroids Pallas and Juno from sample observations.3His method was very similar to the one published in 1965 by Cooley and Tukey, who are generally c

4、redited for the invention of the modern generic FFT algorithm. While Gauss's work predated even Fourier's results in 1822, he did not analyze the computation time and eventually used other methods to achieve his goal.2 Definition2.1 DFT FormulaThere are many different FFT algorithms involvin

5、g a wide range of mathematics, from simple complex-number arithmetic to group theory and number theory. However there have a fundamental formula, let x0, x1, ,xn be complex numbers. The DFT is defined by the formula:Xk=n=0N-1xne-i2knN k=0,1,N-1For convenience in notation, these equations are often w

6、ritten in terms of the complex quantity,WN=e-j2NWith this notation, the DFS analysis-synthesis pair is expressed as follows:Analysis equation: Xk=n=0N-1xnWNkn.Synthesis equation: xn=1Nn=0N-1XkWN-kn. In both of these equations,Xk and xn are periodic sequences. We will sometimes find it convenient to

7、use the notationxnDFTXk2.2 Speed and core theory of FFTEvaluating above definition directly requires O(N2) operations: there are N outputs Xk, and each output requires a sum of N terms. An FFT is any method to compute the same results in O(N log N) operations. More precisely, all known FFT algorithm

8、s require O(N log N) operations (technically, O only denotes an upper bound), although there is no known proof that a lower complexity score is impossible.To illustrate the savings of an FFT, consider the count of complex multiplications and additions. Evaluating the DFT's sums directly involves

9、 N2 complex multiplications and N(N1) complex additions of which O(N) operations can be saved by eliminating trivial operations such as multiplications by 1. The well-known radix-2 CooleyTukey algorithm, for N a power of 2, can compute the same result with only (N/2)log2(N) complex multiplications (

10、again, ignoring simplifications of multiplications by 1 and similar) and Nlog2(N) complex additions. In practice, actual performance on modern computers is usually dominated by factors other than the speed of arithmetic operations and the analysis is a complicated subject (see, e.g., Frigo & Joh

11、nson, 2005), but the overall improvement from O(N2) to O(N log N) remains.The regular computation algorithms, both DIT-FFT and DIF-FFT, are established on above thoughts.3 DIF-FFT3.1 Decimation-In-Frequency AlgorithmsIn computing the DFT, dramatic efficiency results from decomposing the computation

12、into successively smaller DFT computation. In this process, we exploit both the symmetry and the periodicity of thee complex exponential WNkn=e-j2Nkn.Algorithms in which the decomposition is based on decomposing the sequence xn into successively smaller subsequences are called decimation-in-algorith

13、ms.To develop these FFT algorithms, let us restrict the discussion to N a power of 2 and consider computing separately the even-numbered frequency samples and odd- numbered frequency samples. Since the even-numbered frequency samples are X2r=n=0N-1xnWNn(2r), r=0,1,N/2-1,Which can be expressed as, Eq

14、1X2r=n=0N/2-1xnWN2rn+n=N/2N-1xnWN2rnWith a substitution of variables in the second summation in Eq, we obtainX2r=n=0N/2-1xnWN2rn+n=0N/2-1xn+(N/2)WN2r(n+(N/2)Finally, because of the periodicity of WN2rn,WN2rn=WN2rnWNrN=WN2r(n+N/2)Then, Eq1 can be expressed as,X2r=n=0N/2-1(xn+xn+N/2)WN/22rn, r=0,1,N/2

15、-1With the same deduction processes, we can obtain odd-numbered frequency points, given byX2r+1=n=0N-1xnWNn(2r+1), r=0,1,N/2-1Finally we can obtain Eq2X2r=n=0N2-1xn-xn+N/2WNnWN/2rn, r=0,1,N/2-1.3.2 Flow graph of decimation-in-frequencyAccording above equation (Eq1 and Eq2), we draw below Butterfly F

16、low graph,Figure 1 Flow graph of a typical 2-point DFT as required in the last stage of DIF-FFTWe note that since N is power of 2,N/2 is even; consequently, the N/2-point DFTs can be computed by computing the even-numbered and odd-numbered for those DFTs separately. We can obtain below graph of DIF

17、decomposition of an 8-point DFT four 2-point DFT computations,Figure 2 Flow graph of DIF decomposition of an 8-point DFT four 2-point DFT computations3.3 Programming3.3.1 Processing of FFT programmingFrom Flow graph of decimation-in-frequency decomposition of 8-point DFT, the basic component of FFT

18、computation is 2-point DFT Butterfly Flow Graph, it can be implemented by several comment codes. Then use loop code to compute every unit.For 8-point Butterfly computation, it has 3 class. The first class have 1 group and contains 4 rotation elements: W80，W81，W82，W83；the second class have 2 groups w

19、hich include 2 butterfly units , and contain 2 rotation elementsW40, W41; the third class have 4 groups, each group has 1 butterfly unit, and the rotation element W20. Actually, the rotation index follows the equation:r=J*2M-1, J=0,1,2,2M-1-1M is class, J is point in M class. Then we can set steps o

20、f loop programming.STARTSet input sequence xnConsider N=2m(if n<N, add 0 for rest sequence)Loop-for m=1 to 3Corresponding rotation elements WNM-1JLoop-for 1 to 2m-1 groupsLoop-for 1 to 23-m groupCompute Butterfly unit by Eq1&Eq2Reverse the sequence code and draw outputENDYYYNNNFigure 3 Progra

21、m flow graph3.3.2 Programming Code1) Use Matlab fuction fft()N=8;%8-point n=0:N-1;%axis of inputx=0 2 4 6 0 2 4 6 ;%initial sequence xny=fft(x,N)%call fft() directly and obtain y=X(k)mag=abs(y); %Amplitude subplot(2,1,1);stem(n,x);%Original sequence graphtitle('输入序列x(n)');subplot(2,1,2);stem

22、(n,mag);%plot X(k)title('8点调用FFT函数计算结果')2) Program DIF-FFT m-code :a. Function partFunction Xk=DIF_FFT(xn,N) %Compute the class M and add zero-squence for n<N situation %This DIF_FFT function is callback function N=log2(2nextpow2(length(xn) N=2n If length(xn)

23、<N Xn=xn,zeros(1,N-length(xn) EndM=log2(N) %Start loop for class layerFor m=0 :M-1Num_of_Group=N/2m Iterval_of_Unit=N/2(m+1) Cycle_Count=Interval_of_Unit-1 Wn=exp(-j*2*pi/Interval_of_Group) %start loop for group layerFor g=1 :Num_of_Group Interval_1=(g

24、-1)*Interval_of-Group Interval_2=(g-1)*Interval_of_Group+Interval_of_Unit %Start loop for counting butterfly unitsfor r=0 :Cycle_Count  k=r+1  xn(k+Interval_1)=xn(k+Interval_1)+xn(k+Interval_2); xn(k+Interval_2)=xn(k+Interval_1)-xn(k+Interval_2)*Wnr  end  endendn1=

25、fliplr(dec2bin(0 :N-1) n2=bin2dec(n1) for i=1 :N Xk(i)=xn(n2(i)+1) Endb. Test part(main function)xn=1 :8 m=1 :8 N=8 x1=DIF_FFT(xn,N) x2=abs(x1) x4=abs(x2) x5=angle(x1) x6=angle(x2) subplot(3,1,1)stem(m,xn),grid subplot(3

26、,1,2)stem(m,x3),grid title(Amplitude spectrum)subplot(3,1,3)stem(m,x5),grid title(Phase Spectrum)figure(2)subplot(311)stem(m,xn),grid title(input discrete sequence )subplot(313)stem(m,x4),grid title(The return Amplitude Spectrum of fft)subplot(313)stem(m,x6),grid title(

27、The return Phase Spectrum of fft)3.3.3 ResultsFigure 5 Output sequence by fft() functionFigure 6 Output sequence by DIF_FFT() function4 Application4.1 Image ProcessingAs I know, there have a research about rice phenotypic traits detection in Insitute of Mechatronics & Logistic Equipment, which m

28、ay use image processing. And there have some work accomplished by other scholar. “Rice phenotypic traits can be divided into growth traits and yield traits, this study focuses on obtaining rice yield traits acqusition, including panicle length measurement, rice threshing and grain traits extraction. However, traditional measurement method is manual by ruler and steelyard which is slow and inconsistent. Therefore the traditional method is far from meeting