基2fft算法

A. FFT算法分几种

FFT算法分析FFT算法的基本原理是把长序列的DFT逐次分解为较短序列的DFT。按照抽取方式的不同可分为DIT-FFT（按时间抽取）和DIF-FFT（按频率抽取）算法。按照蝶形运算的构成不同可分为基2、基4、基8以及任意因子（2n,n为大于1的整数），基2、基4算法较为常用。 网上有帮助文档: http://www.5doc.com/doc/123035(右上角有点击下载)

B. 基2—fft算法的软件实现(MATLAB代码)

参考网络： clc; clear all; close all; x=ones(1,128); %输入的信号，自己可以改变 %整体运用原位计算 m=nextpow2(x);N=2^m; % 求x的长度对应的2的最低幂次m if length(x)<N x=[x,zeros(1,N-length(x))]; % 若x的长度不是2的幂，补零到2的整数幂 end nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1; % 求1:2^m数列序号的倒序 y=x(nxd); % 将x倒序排列作为y的初始值 for mm=1:m % 将DFT作m次基2分解,从左到右，对每次分解作DFT运算,共做m级蝶形运算，每一级都有2^(mm-1)个蝶形结 Nz=2^mm;u=1; % 旋转因子u初始化为WN^0=1 WN=exp(-i*2*pi/Nz); % 本次分解的基本DFT因子WN=exp(-i*2*pi/Nz) for j=1:Nz/2 % 本次跨越间隔内的各次蝶形运算,在进行第mm级运算时需要2^(mm-1)个 蝶形 for k=j:Nz:N % 本次蝶形运算的跨越间隔为Nz=2^mm kp=k+Nz/2; % 蝶形运算的两个因子对应单元下标的关系 t=y(kp)*u; % 蝶形运算的乘积项 y(kp)=y(k)-t; % 蝶形运算 y(k)=y(k)+t; % 蝶形运算 end u=u*WN; % 修改旋转因子,多乘一个基本DFT因子WN end end y y1=fft(x) %自己编的FFT跟直接调用的函数运算以后的结果进行对比

C. 快速傅里叶变换的算法类型

FFT算法很多，根据实现运算过程是否有指数因子WN可分为有、无指数因子的两类算法。
有指数因子的算法
经典库利-图基算法 当输入序列的长度N不是素数(素数只能被1而它本身整除）而是可以高度分解的复合数，即N=N1N2N3…Nr时，若N1=N2=…=Nr=2，N=2则N点DFT的计算可分解为N=2×N/2，即两个N/2点DFT计算的组合，而N/2点DFT的计算又可分解为N/2=2×N/4，即两个N/4点DFT计算的组合。依此类推，使DFT的计算形成有规则的模式，故称之为以2为基底的FFT算法。同理，当N=4时，则称之为以4为基底的FFT算法。当N=N1·N2时，称为以N1和N2为基底的混合基算法。
在这些算法中，基2算法用得最普遍。通常按序列在时域或在频域分解过程的不同，又可分为两种：一种是时间抽取FFT算法（DIT），将N点DFT输入序列x(n)、在时域分解成2个N/2点序列而x1(n)和x2(n)。前者是从原序列中按偶数序号抽取而成，而后者则按奇数序号抽取而成。DIT就是这样有规律地按奇、偶次序逐次进行分解所构成的一种快速算法。
分裂基算法(RSFFT) 1984年由P．杜哈美尔和H．赫尔曼等导出的一种比库利图基算法更加有效的改进算法，其基本思想是在变换式的偶部采用基2算法，在变换式的奇部采用基4算法。优点是具有相对简单的结构，非常适用于实对称数据，对长度N=2能获得最少的运算量(乘法和加法），所以是选用固定基算法中的一种最佳折衷算法。

D. Matlab的时间抽取基2FFT算法

基于Matlab的时间抽取基2FFT算法
function y=myditfft(x)
%本程序对输入序列实现DIT-FFT基2算法，点数取大于等于长度的2的幂次
%------------------------------------
% Leo's fft program(改编网上的一个程序)
%------------------------------------
m=log2(2^nextpow2(length(x))); %求的x长度对应的2的最低幂次m
N=2^m;
if length(x)<N
x=[x,zeros(1,N-length(x))]; %若长度不是2的幂，补0到2的整数幂
end
x;
%--------------------------------------------------------------------------
%对输入序列进行倒序
%如果输入序列的自然顺序号I用二进制数（例如n2n1n0）表示
%则其倒位序J对应的二进制数就是（n0n1n2），这样，在原来自然顺序时应该放x(I)的
%单元，现在倒位序后应放x(J)。
%--------------------------------------------------------------------------
%以下程序相当于以下程序:
%nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1; %求1：2^m数列的倒序
%y=x(nxd); %将倒序排列作为初始值
%--------------------------------------------------------------------------
NV2=N/2;
NM1=N-1;
I=0;
J=0;
while I<NM1
if I<J
T=x(J+1);
x(J+1)=x(I+1);
x(I+1)=T;
end
K=NV2;

while K<=J
J=J-K;
K=K/2;
end
J=J+K;
I=I+1;
end
x;
%--------------------------------------------------------------------------
%以下程序解释:
%第一级从x(0)开始,跨接一阶蝶形,再取每条对称
%第二级从x(0)开始,跨接两阶蝶形,再取每条对称
%第m级从x(0)开始,跨接2^(m-1)阶蝶形,再取每条对称....
%--------------------------------------------------------------------------
for mm=1:m %将DFT做m次基2分解，从左到右，对每次分解作DFT运算
Nmr=2^mm;
u=1; %旋转因子u初始化
WN=exp(-j*2*pi/Nmr); %本次分解的基本DFT因子WN＝exp(-i*2*pi/Nmr)
for n=1:Nmr/2 %本次跨越间隔内的各次碟形运算
for k=n:Nmr:N %本次碟形运算的跨越间隔为Nmr=2^mm
kp=k+Nmr/2; %确定碟形运算的对应单元下标(对称性)
t=x(kp)*u; %碟形运算的乘积项
x(kp)=x(k)-t; %碟形运算的加法项
x(k)=x(k)+t;
end
u=u*WN; %修改旋转因子，多乘一个基本DFT因子WN
end
end
y=x; %输出

E. 基数为2的FFT算法

从上节所述，FFT算法快速的关键在于将原来傅氏矩阵分解为每一行仅含有两个非零项l与Wⁱ的矩阵的乘积。下面用基数为2，即N=2ⁿ的情形讨论矩阵的分解过程.并主要按时间分解的情况讨论。

按时间分解的FFT算法

设N=2ⁿ，n为正整数。考虑输入序列x₀(l)的DFT

物探数字信号分析与处理技术

将l与m用二进制表示

物探数字信号分析与处理技术

将(7-2-2)代入(7-2-1)中，得到

物探数字信号分析与处理技术

为了说明问题，我们取N=8，于是从(7-2-2)得到

物探数字信号分析与处理技术

从(7-2-4)和(7-2-3)得到

物探数字信号分析与处理技术

将(7-2-5)中的W的指数按时间l进行分解，有

物探数字信号分析与处理技术

因为

物探数字信号分析与处理技术

故从(7-2-6)得到

物探数字信号分析与处理技术

将上式代入(7-2-5)中得到

物探数字信号分析与处理技术

物探数字信号分析与处理技术

我们在公式(7-2-7)中由里往外求和，并置

物探数字信号分析与处理技术

于是得到

物探数字信号分析与处理技术

首先写出(7-2-8)的所有式子

物探数字信号分析与处理技术

将方程组(7-2-12)写成矩阵形式，得到

物探数字信号分析与处理技术

我们看到(7-2-13)中的方阵，正好是(7-1-13)中的方阵，也就是(7-1-12)中被分解出来的第3个矩阵，只不过这里的x₁(l)与x₀(l)中的l是用二进制数表示而已。

再写出(7-2-9)的所有式子，得到

物探数字信号分析与处理技术

将方程组写成矩阵形式，则有

物探数字信号分析与处理技术

显然(7-2-15)中的矩阵，正好是(7-1-14)中的方阵，也就是(7-1-12)中被分解出来的第二个矩阵，只不过这里的x₂(l)与x₁(l)是用二进制数表示而已，最后将(7-2-10)中的全部式子写出，得到

物探数字信号分析与处理技术

将方程组(7-2-16)写成矩阵形式，则有

物探数字信号分析与处理技术

显然，(7-2-17)中的方阵正是(7-1-15)中的方阵，也就是(7-1-12)中被分解出来的第1个矩阵，只不过这里的x₃(l)与x₂(l)中的l是用二进制数表示。

由此可见，(7-2-7)中由里往外的三个求和式(7-2-8)、(7-2-9)及(7-2-10)，完全确定了(7-1-12)中三个被分解的矩阵因子。求和得到的最终结果x₃(m₀，m₁，m₂)，与我们所要求的X(m₂，m₁，m₀)正好是逆序的。

到此为止，我们就看到(7-1-11)中的方阵是怎样被分解成三个方阵因子的。对于N=8，方程(7-2-8)～(7-2-11)就是计算DFT的FFT算法。为了对FFT算法有一个直观的了解并便于编制程序，我们以N=8为例，画出其流程图(图7-2-1)。

根据(7-2-13)，将其中的W⁴用-W⁰代替，画出从x₀(r)到x₁(r)的流程图。这一迭代过程用符号#1表示;再根据(7-2-15)，将其中的与W⁴、W⁶分别换成-W⁰与-W²，画出从x₁(r)到x₂(r)的流程图，这一迭代过程记为#2;最后，根据(7-2-17)，将其中的W⁴、W⁶、W⁵、W⁷分别换成-W⁰、-W²、-W¹、-W³，画出流程图7-2-3合并图7-2-1～7-2-3，就得到从x₀(r)到x₃(r)的完整流程图7-2-4。

图7-2-1 第一次递推

图7-2-2 第二次递推

在图7-2-5中，画出N=16=2⁴的FFT算法流程图:

根据从x₀到谱X的FFT算法流程图7-2-4与图7-2-5，我们总结出如下几点:

(1)从x₀到终值x_r的最大迭代次数r，由r=log₂N确定。

例如，N=8，最大迭代次数r=3;N=16，最大迭代次数r=4。

(2)在第r次迭代中，要乘的W因子为

图7-2-3 第三次递推

例如，N=8，在第一次迭代中，要乘因子W⁰;在第二次迭代中要乘因子W⁰，W¹，W²，W³。

(3)在第r次迭代中，包含2^r-1个组，每个组

包含 。例如N=8，第一次迭代r=1，有

一个组，每组包含8个x(s);在第二次迭代中包含2个组，每组包含4个x(s);第三次迭代中包含4个组，每组2个x(s)。

图7-2-4 x₀(r)到x₃(r)的完整流程图

(4)在第r次迭代中，各组包含的W因子各不相同，且每一组仅包含一种类型的因子 ，此因子在组中一半数为正，另一半数为负。例如N=8，第二次迭代中，第一组包含因子W⁰，且在该组中一半数为正，另一半数为负;第二组包含因子W²，在该组中也是一半数为正，另一半数为负。

(5)在第r次迭代中，各组包含的W因子除正负号外类型均相同。所以只须确定每组中第一个因子，之后按半数反号，即得到所求W因子。具体做法是，在每组第一个因子W^SN2r对应的x_r(k)中，将k表示成n位的二进制数，n=log₂N，并把这个二进制数右移n-r位，左边空出的地方添零补足n位，之后再将此n位二进制数逆位，即得到所求W因子的指数。例如，N=8，迭代#2包含两组，每组包含4个x₂(k)，第二组第一个因子W对应于x₂(4)。将4表示成3位的二进制数为100，把它右移1位成10，右边添零补成3位为010，逆位仍为010，故所求因子为W²，第2组全部W因子为W²，W²，-W²，-W²。又如，N=16，迭代#3中包含4个组，每组包含4个x₃(k)，第4组第1个因子W对应于x₃(12)。将12表示成4位的二进制数为1100。右移1位变成110，将左边空处添零补成4位为0110，逆位仍为0110，故所求因子为W⁶，从而第4组全部W因子为W⁶，W⁶，-W⁶与-W⁶。

图7-2-5 N=16=2⁴的FFT算法流程图

(6)如果已知N=2的FFT算法，容易从它求得n=2的FFT算法。具体作法是，在n=2^n-1的FFT算法中，将所有x_r(l)的个数加倍，所有W的个数及其乘幂加倍，就得N=2ⁿ中前n-1次迭代的全部结果。之后，将新得到的第n-1次迭代中乘幂相同的W个数减半，就是第n次迭代中前2ⁿ/2个W，将这些W的乘幂依次加1，就得到后2ⁿ/2个W。例如N=16的前3次迭代，都是N=8的三次迭代中所有x_r(l)的个数加倍，W的个数及其乘幂加倍的结果。再将N=16的第三次迭代中乘幂相同的W个数减半，就是第4次迭代中前8个W。

(7)在第r-1次迭代中的x_r-1(i)与x_r-1(j)仅用于计算r次迭代中的x_r(i)与x_r(j)，而不会用于计算任何其它的x_r(k)与x_r(l)。例如N=16的第二次迭代中的x₂(0)与x₂(2)，只用于计算第三次迭代中的x₃(0)与x₃(2);第三次迭代中的x₃(8)与x₃(9)也只用于计算第四次迭代中的x₄(8)与x₄(9)。因此，我们可以把第r次迭代中的x_r(i)与x_r(j)，存放到r-1次迭代x_r-1(i)与x_r-1(j)所占用的原存储单元中去，这样，所需要的计算机存储容量就只限于原数据序列占据的单元数。如果是复序列的话，存储单元要加倍。

(8)上述FFT算法也可用于计算逆离散傅氏变换(IDFT)(图7-2-6)，只不过在计算时要将上述FFT算法中的W因子用其共轭W^*代替，并将最后迭代计算的结果统统乘以1/N.

图7-2-6 N=8的逆离散富氏变换流程图