基数为2的FFT算法

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

按时间分解的FFT算法

设N=2n，n为正整数。考虑输入序列x0(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个矩阵，只不过这里的x1(l)与x0(l)中的l是用二进制数表示而已。

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

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

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

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

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

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

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

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

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

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

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

根据(7-2-13)，将其中的W4用-W0代替，画出从x0(r)到x1(r)的流程图。这一迭代过程用符号#1表示再根据(7-2-15)，将其中的与W4、W6分别换成-W0与-W2，画出从x1(r)到x2(r)的流程图，这一迭代过程记为#2最后，根据(7-2-17)，将其中的W4、W6、W5、W7分别换成-W0、-W2、-W1、-W3，画出流程图7-2-3合并图7-2-1～7-2-3，就得到从x0(r)到x3(r)的完整流程图7-2-4。

图7-2-1 第一次递推

图7-2-2 第二次递推

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

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

(1)从x0到终值xr的最大迭代次数r，由r=log2N确定。

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

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

图7-2-3 第三次递推

例如，N=8，在第一次迭代中，要乘因子W0在第二次迭代中要乘因子W0，W1，W2，W3。

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

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

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

图7-2-4 x0(r)到x3(r)的完整流程图

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

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

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

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

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

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

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

按照程春蠢亏序框图来看是这样的

1.for循环输出了一个文件中的一维数组，这个档判数组是经过了特殊索引规则后的数组。

2.一维数组连接了一个拆分一维数组，前512个元素构成的数组是时域信号，后面的数据是单片机的频谱。

3.时域信号连接了Butterworth滤波扒神器，经滤波后的数据再连接一个汉宁窗，FFT（X）变换得到一个复数组，做复数至极坐标转换之后取r这列数组，除以数组大小，实际是得到了幅度谱，连接下面theta的话就是做的是相位谱。这个步骤完全是FFT（x）的算法，你拿来用就行。

4.做了两次平方这，我明白第一次平方有可能是为了做功率谱，但是第二次何解？有可能是编程写错了

5.告诉你一个更好的方法，直接用“编程——信号处理——波形测量——FFT频谱（幅度-相位）”这个函数

二维FFT相当于对行和列分别进行一维FFT运算。

先对各行逐一进行一维FFT，然后再对变换后的新矩阵的各列逐一进行一维FFT。相应的伪代码如下所示：for (int i=0i<Mi++)FFT_1D(ROW[i]，N)for (int j=0j<Nj++)FFT_1D(COL[j]，M)其中，ROW[i]表示矩阵的第i行。

例：

#include <stdio.h>

#include <math.h>举老液

#include <stdlib.h>

#define N 1000

/*定义复数类型*/

typedef struct{

double real

double img

}complex

complex x[N], *W/*输入序列,变换核*/

int size_x=0/*输入序列的大小，在本程序中仅限2的次幂*/

double PI/*圆周率*/

void fft()/*快速傅里叶变换*/

void initW()  /*初始化变换核*/

void change()/*变址*/

void add(complex ,complex ,complex *)/*复数加法*/

void mul(complex ,complex ,complex *)/*复数乘法*/

void sub(complex ,complex ,complex *)/*复数减法*/

void output()

int main(){

int i/*输出结果*/

system("cls")

PI=atan(1)*4

printf("Please input the size of x:\n")

scanf("%d",&size_x)

printf("Please input the data in x[N]:\n")

for(i=0i<size_xi++)

scanf("%lf%lf",&x[i].real,&x[i].img)

initW()

fft()

output()

return 0

}

/*快速傅里叶变换*/

void fft(){

int i=0,j=0,k=0,l=0

complex up,down,product

change()

for(i=0i<log(size_x)/log(2) i++){   /*一级蝶形运算*/

l=1<<i

扩展资料：

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为基底的混合基算法。

参考资料来源：百度百科-快速傅里叶变换

---