用MATLAB编写FFT程序

一、基2FFT的原理

基2FFT算法分为两类：时域抽取法FFT（Decimation-In-Time FFT，简称DIT-FFT）；频域抽取法FFT（Decimation-In-Frequently FFT，简称DIF-FFT）。这里只做了DIT-FFT的实验。

设序列 的长度为N，且满足 ，M为自然数。按n的奇偶把 分解为两个 点的子序列。然后分别对其进行DFT，最终可得以下式子：

       

   

根据以上两个式子可以画出蝶形图如下：

根据以上蝶形图的画法，可以画出8点DIF-FFT的运算流图：

二、DIF-FFT的运算规律及编程思想

1、原位计算

由8点DIF-FFT的运算流图可以看出，点的FFT共进行了M级运算，每级由N/2个蝶形运算组成。同一级中，每个蝶形运算的两个输入数据只对计算本蝶形有用，而且每个蝶形的输入和输出又在同一平线上，所以计算完一个蝶形运算后，所得的输出数据可立即存入原输入数据所占有的内存空间。这叫做原位（址）计算。

2、旋转因子

N点DIF-FFT每级都有N/2个蝶形，每个蝶形都要乘以因子，称为旋转因子，p为旋转因子指数。L表示运算的级数。

由归纳可得，对于的情况有：

      

3、蝶形运算

设序列经时域抽选（倒序）后存入数组A中。如果蝶形运算的两个输入数据相距B个点，则蝶形运算可以表示成如下：

其中：      ；

下标L表示第L级运算。

4、序列的倒序

流程图如下所示：

5、程序框图

三、MATLAB程序

这是基2的FFT，8点。

N=8;

M=3;

LH=N/2;

J=LH+1;

N1=N-1;

A=[1 2 3 4 5 6 7 8];

WN=[exp(j*2*pi*0/N)  exp(j*2*pi*1/N)  exp(j*2*pi*2/N)  exp(j*2*pi*3/N)];

for I=2:N1

    if (I>=J)

        K=LH;

        while(J>=K)

            J=J-K;

            K=K/2;

        end

        J=J+K;

    else

        T=A(I);

        A(I)=A(J);

        A(J)=T;

         K=LH;

        while(J>=K)

            J=J-K;

            K=K/2;

        end

        J=J+K;

    end

end

for L=1:M

    B=2^(L-1);

    for J=0:B-1

        P=2^(M-L)*J;

        for k=J+1:2^L:N-1

            T=A(k)+A(k+B)*WN(P+1);

            A(k+B)=A(k)-A(k+B)*WN(P+1);

            A(k)=T;

        end

    end

end

四、结果分析

使用MATLAB库函数FFT（），将结果做比较。

N=8;

n=0:N-1;

xn=[1 2 3 4 5 6 7 8];

Xk=fft(xn);

结果比较：

这是自己写的FFT运行结果：

这是MATLAB的库函数运行出来的结果：