快速傅里叶变换(FFT)算法原理及代码解析

回复 星标
更多

快速傅里叶变换(FFT)算法原理及代码解析

510702

  • FFT与DFT关系:

快速傅里叶变换(Fast Fourier Transform)是离散傅里叶(DFT)变换的一种快速算法,简称FFT,通过FFT可以将一个信号从时域变换到频域;FFT(快速傅里叶变换)其本质就是DFT,只不过可以快速的计算出DFT结果,它只是傅立叶变换算法实现过程的一种改进。

510702

要弄懂FFT,必须先弄懂DFT,DFT(DiscreteFourier Transform) 离散傅里叶变换的缩写,咱们先来简单讨论一下DFT。DFT(FFT)的作用:可以将信号从时域变换到频域,而且时域和频域都是离散的,通俗的说,可以求出一个信号由哪些正弦波叠加而成,求出的结果就是这些正弦波的幅度和相位。

510702

  • DFT的公式:

510702

其中X(k)表示DFT变换后的数据,x(n)为采样的模拟信号,公式中的x(n)可以为复信号,实际当中x(n)都是实信号,即虚部为0,此时公式可以展开为:

510702

那么,对于一个 的序列进行不断分解,就可以得出如下所谓的蝶形图:

510702

  • FFT处理蝶形运算

蝶形运算的规律:同一级中所有蝶形的输入点在同一竖直线上,意味着我们可以按级来运算,对于M级的蝶形,编个M次循环就好了;所有数据点在运算后不会窜位,即计算后可以将结果存入原来的地址空间。

每级N/2个蝶都需要用到系数WN,这里称它为旋转因子。我们来观察旋转因子WN的规律。以8点的蝶形图为例:

510702

可见,第L级的旋转因子为:

510702

可以看到,每个蝶的两个输入点下标跨度是不一样的。比如第一级中是相邻两个数据作蝶运算,第二级中是两个输入点隔开一个点,而第三级隔开三个点。不难找到规律:第L级中,第二个输入点的坐标是第一个点的坐标+space,space=Math.Pow(2, L)=num。

FFT的算法是写一个三重循环:

  • 第一重循环对每一级运算(每级含num=Math.Pow(2, L)个蝶形);

  • 第二重对每一个旋转因子对应的蝶运算, 那么有几个蝶呢?很简单,每级都应该有N/2个蝶,而每个因子对应N/2 / num个蝶;

  • 第三重循环对每个蝶进行计算,需要注意的一是循环下标起始点的位置,二是每次计算需要申明临时变量来保存输入数据点。

实现代码:

510702

FFT算法的原理是通过许多小的更加容易进行的变换去实现大规模的变换,降低了运算要求,提高了与运算速度。FFT不是DFT的近似运算,它们完全是等效的。

510702

  • 傅里叶变换的C语言编程

对于快速傅里叶变换FFT,第一个要解决的问题就是码位倒序。假设一个N点的输入序列,那么它的序号二进制数位数就是t=log2N。码位倒序要解决两个问题:

  1. 将t位二进制数倒序;

  2. 将倒序后的两个存储单元进行交换。如果输入序列的自然顺序号i用二进制数表示,例如若最大序号为15,即用4位就可表示n3n2n1n0,则其倒序后j对应的二进制数就是n0n1n2n3。

代码如下:

复数类型定义及其运算

#define N 64 //64点

#define log2N 6 //log2N=6

/*复数类型*/

typedef struct

{

float real;

float img;

}complex;

complex xdata x[N]; //输入序列

/*复数加法*/

complex add(complex a,complex b)

{

complex c;

c.real=a.real+b.real;

c.img=a.img+b.img;

return c;

}

/*复数减法*/

complex sub(complex a,complex b)

{

complex c;

c.real=a.real-b.real;

c.img=a.img-b.img;

return c;

}

/*复数乘法*/

complex mul(complex a,complex b)

{

complex c;

c.real=a.real*b.real - a.img*b.img;

c.img=a.real*b.img + a.img*b.real;

return c;

}

/***码位倒序函数***/

void Reverse(void)

{

unsigned int i,j,k;

unsigned int t;

complex temp;//临时交换变量

for(i=0;i<N;i++)//从第0个序号到第N-1个序号

{

k=i;//当前第i个序号

j=0;//存储倒序后的序号,先初始化为0

for(t=0;t<log2N;t++)//共移位t次,其中log2N是事先宏定义算好的

{

j<<=1;

j|=(k&1);//j左移一位然后加上k的最低位

k>>=1;//k右移一位,次低位变为最低位

}

if(j>i)//如果倒序后大于原序数,就将两个存储单元进行交换(判断j>i是为了防止重复交换)

{

temp=x;

x=x[j];

x[j]=temp;

}

}

}

510702

  • 总结:

FFT是离散傅立叶变换的快速算法,可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的,但是如果变换到频域之后,就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外,FFT可以将一个信号的频谱提取出来,这在频谱分析方面也是经常用的。

此帖已被锁定,无法回复
新窗口打开 关闭