目录
开场白
最近研究数字滤波器的时候接触到了离散傅里叶变换(DFT),它的实现非常简单但是使用非常不便,时间复杂度很高,对于性能一般且不带FPU的辣鸡单片机直接求解难度很大,在点数较少( )时勉强可用,点数再高计算几乎登天。因此需要一种能够快速计算DFT的算法减少时间复杂度,此快速傅里叶变换(FFT)的初衷是也。虽然现在有相当多的FFT库可以使用,这些库性能都相对不错,敲起代码来得心应手,丝毫没有无知的羞耻(^_^)。但作为一只好奇心甚重的小白,还是希望亲自实现一把,于是乎趁晚上散步的几个小时找来书研究了一下,推了一遍公式,敲了一遍代码,也算是搞懂了FFT的实现。
基础数学
在涉及傅里叶变换的核心算法之前,先对其中的基础数学简要归纳
复数由实部和虚部构成,形如(
),这里
为单位虚数(
)
单位虚数的性质:
复数也有四则运算,这里我们用到其中两种:复乘、复加
1.复数运算
对于两个复数
、
,其乘积满足
对于两个复数
、
,其和
2.欧拉公式
对于复指数形式,欧拉大人为我们指了一条明路,这就是著名的欧拉公式
这条公式将奠定快速傅里叶变换的理论基础
基础信号论
奈奎斯特采样定理
奈奎斯特采样定理(又名香农采样定理)是对信号离散化的准则,为了获取信号更多的信息,其必须满足
采样速率 必须大于信号最高频率的两倍,这个频率称为奈奎斯特频率 ,即
从解方程到傅里叶变换
1.何方神圣
法国大数学家傅里叶提出,任何一个函数都可以化为由正弦函数构成的无穷级数,简而言之,表现任何一个复杂的信号(如声音),只需要足够多的正弦函数相加。例如信号
可以写为
这里每一个小的正弦函数都包含属于自己的幅值、频率、相位,这样一个庞大的组合体构成了任意的 ,也构成我们能看到、听到的信息,狭义的说,傅里叶变换正是将一个完整的 打散为正弦之和的手段,每一个正弦代表了信号在该频率的分量,于是我们说xx信号200KHz的分量是…
2.分量求解:从方程开始
在上面我们提到任意信号都可以由不同频率的分量合成,也就是
如果我们换一个角度看上式,看成一个方程,结果如何?为了简单起见,我们先抛弃令人不舒服的相位
这里,每一个正弦函数我们都抽象地将其看作一个种频率,即
,我们把这些频率看作已知,也就是统统都存在(大不了幅值为0嘛)
于是我们需要求解上述方程。经验告诉我们一个方程式不够嘀,有多少个未知数就有多少个方程组成方程组,否则方程组就是欠定,解可能不唯一。很不幸这里
的取值为无穷大,也就是说我们需要找到无穷多个方程联立才能求解每一个频率分量,但这无疑为我们解决问题找到了一个途径。
事实上我们真的需要知道每一个频率分量吗?
世界上漂亮的女孩子千千万万,难道….
在信号论基础里面俺特意放了一条奈奎斯特采样定律,可见我们能获得的信号细节(也就是频率的种数吧)是有条件限制的,至少受采样速率的限制,采样速率又受产品要求、器件工艺、人民币张数等等的限制….
假若以200Hz的速率采样,能准确看清的信号频率最高也就100Hz吧?200Hz总得有个分辨率,分辨率为1Hz,就需要201个正弦函数。对!201个正弦函数,到此问题可解:为了在采样频率为200Hz时得到分辨率为1Hz的分量信息,需要求解未知数为201的方程:
这个方程求解可以采用LUP法,前面博文中详细介绍过。
但是这仅仅时幅值信息,相位信息又在何处?
又如何求解?
3.得到DFT(可读)
以下为纯数学推导,仅仅阐释思路,若不了解亦无大碍,可直接掠过。传说每多增加一个公式将会少去一半人
为了求解
我们对包含相位的正弦用公式展开,这里我们单独考虑这一部分
到这里三角变形貌似没有特别好的办法了,我们看到欧拉公式
对两式做加减可以得到
带入
可以得到
到此似乎依然很复杂,没关系先留着
上面我们只考虑了n=1的情况,对于任意n
令 ,带入之可得
到这里我们得到了
与
的关系,注意我们需要利用这个关系从中求解出
到这里肯定小伙伴们都晕了,没关系缓一缓仔细看看前面的东西,好好回味一下,下方高能预警!!!
为求解
我i们不惜代价地在等式左右两边同时乘上
并积分,看起来有点复杂我们一步一步来
先乘:
到这里应该还是容易的,只是莫名其妙多了一个k而已
下面等号两边同时积分
积分号搬进求和符号内
此时奇迹出现了:一个周期内当 时积分项为零,否则积分项为T,现在非常清楚了
即
等价于
正式的对于点数为N的傅里叶变换,我们有
实际上从正弦到DFT我们拼拼凑凑得到了这样的一个式子,中间省略了太多的细节与积分变换,推导方式也不严谨。
快速傅里叶变换的原理与细节
1.旋转因子
在这一节开始之前我们先做约定:定义一个旋转因子
其中
为单位虚根 ,那么理所当然的旋转因子的k次方就是
注意这个小小的旋转因子跟上面那个角频率不是一个东西
为什么这个叫做旋转因子?旋转因子如何转起来?
首先我们需要一个复平面
在这个复平面上我们表示了n=8,n=4,n=2的旋转因子
在这里我们观察到向量每转过八分之一个圆周,n=8的旋转因子上标就加上1,可以想见随着上标的增加,向量周而复始的运动,例如
不仅仅是下标为8的旋转因子拥有这个性质,任何一个旋转因子都具有这个性质哦
我们再来观察一下从
和
好像是一个东西,再一般一点
和
,也是一个东西,好像旋转因子的上下角标是可以约分的(^_^),这一点非常重要!!我们对这个性质加以概括
当n为偶数时
这个定理称为折半定理,他向我们揭示了旋转因子的计算可以蜕化为两个一半规模的计算,这点奠定了快速计算DFT的基础。
2.DFT到FFT
上面我们定义了一个旋转因子,再回到DFT,我们可以这么表示
显然这个是个关于
的多项式,他的最高次数为N-1次。
但是上面那个表达方式实在难懂,我们把它写成接地气的形式
为了能利用旋转因子的折半定理,我们把其拆成奇数和偶数两部分
(注:even,偶数;odd,奇数)
显然
使用折半定理,我们又可以 成为
现在我们把求解分割成为两部分,两部分都是独立的DFT。对于 h和 这两个DFT,折半定理依然起作用,我们可以不断拆分下去直到成为两个N/2个独立的2点DFT,为了保证N能一直被2整除,需要N为2的整数次,即
这个就是基2的FFT,同理有基4甚至任意基,分裂基的^_^
显然迭代的次数也就时k次。要实现这样的一个算法递归显然是非常好的,实现起来也相当容易,但是递归的效率你懂的,一般我们会将其重写成循环的形式。
3.纸上谈兵
为了能用循环完成FFT,我们不得不把递归的过程扒开,这里我们以N=8为例
每一个框框代表了一次DFT,这里我们最终将N=8的DFT转化为4个N=2的DFT。由于我们每一次选择偶数项、奇数项来完成DFT,这里DFT系数的脚标也在变化,不是01234567,而是04261537!
为了能看清每一个阶段做的操作,我们引入蝶形图
这里每一个交叉就是一次蝴蝶操作,蝴蝶操作定义如下:
很形象吧,虽然长得不像蝴蝶…
在每一次蝴蝶操作中我们需要完成两次复乘、两次复加
既然鲽形图贴出来了,俺们就图论图
图中一共有4各阶段,其实第一阶段我们可以忽略,N=1的DFT就是自己。在第一阶段,最主要的任务就是换换位置。
在第二到第四阶段中,每一个阶段都在重复做蝴蝶操作,我们把每一阶段的蝴蝶操作分组:
第二阶段:4组,每组1次蝴蝶操作,每一组旋转因子只有
第三阶段:2组,每组2次蝴蝶操作每一组旋转因子为
第四阶段,1组,每组4次蝴蝶操作每组旋转因子从 到 递增
注:图中没有表示出旋转因子的N,这里我们自己总结
每一次蝴蝶操作有两个点,上面的那个点我们叫
,下面的叫
我们又可以观察到:
第二阶段
每组最后一次蝴蝶操作的 与下一组第一次的 间隔为1, 与 间隔为1
第三阶段
每组最后一次蝴蝶操作的 与下一组第一次的 间隔为2, 与 间隔为2
第四阶段
没有下一组蝴蝶操作 与 间隔为4
有了这些信息,我们就可以写出精确的伪代码
4.FFT伪代码
function:fft,
array a is a complex number array;
n is length of array a,which is a power of 2;
input:a;
Reverse(a);
for i=0 to log(n)
for j=0 to n/pow(2,i+1)
for k=0 to pow(2,i)
theta=2*Pi*k/pow(2,i+1);
Pre=k+j*pow(2,i+1);
Next=k+j*pow(2,i+1)+pow(2,i);
Butterfly_Operation(Pre,Next,theta);
end for
end for
end for
end function
这段伪代码描述了基2FFT运行的方式,我们来详细看一下
for i=0 to log(n)
这个循环定义了FFT的级数,若N=8,显然为3级
for j=0 to n/pow(2,i+1)
在这个循环中,定义了蝴蝶操作的组数。对照上面我们总结出来N=8的情况:
第一级有4组,第二级有2组,第三级只有一组,组数为N不断除以2得到,这个就是
的来处
for k=0 to pow(2,i)
这个循环定义了每一级每组蝴蝶操作的数量。对照N=8的情况:
第一级每组1个,第二级每组2个,第三集每组4个,很显然是按2的整数次递增的。
theta=2*Pi*k/pow(2,i+1);
Pre=k+j*pow(2,i+1);
Next=k+j*pow(2,i+1)+pow(2,i);
theta是旋转因子的用欧拉公式展开得到的,我们有
Pre表示蝴蝶操作中上面的点,这个点的位置非常复杂,每除了组内递增之外还会隔组递增。上面我们看到了k表示蝴蝶操作的数量,在每组每一次蝴蝶操作中,Pre向下递增1,完成一组中Pre跳到下一组的第一位。对照N=8的情况,以第二级为例:
第二级中每一组间隔为2,每一组两个蝴蝶操作,第二组的第一个蝴蝶操作Pre为4
表明每k次增加一个组间隔,间隔为
,这样一来,Next无非就是在Pre上增加一个Pre与Next之间的间隔,结合我们对N=8的观察,
第一组每个间隔为1,第二个每个间隔为2,第三组每个间隔为4,于是只需要再Pre的基础上增加
即可。
这段伪代码中出现了两个函数,Reverse(),Butterfly_Operation(),这两组分别是位置更换、蝴蝶操作的实现函数。
5.倒位序
观察N=8的情况
我们看到输入脚标为0,1,2,3,4,5,6,7
更换位置以后为0,4,2,6,1,5,3,7
这些数字有什么特征?我们注意到这些数字都是每一组选出来排在偶数、奇数位上的数字,我们不妨把他们转换成二进制
DEC | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
---|---|---|---|---|---|---|---|---|
-BIN | 000 | 001 | 010 | 011 | 100 | 101 | 110 | 111 |
而换过位置以后的脚标
DEC | 0 | 4 | 2 | 6 | 1 | 5 | 3 | 7 |
---|---|---|---|---|---|---|---|---|
-BIN | 000 | 100 | 010 | 110 | 001 | 101 | 011 | 111 |
额看出来了吗?前者是从低向高进位,后者是从高向低进位。无论如何,知道这点就有很多种方法来实现了,比如说查表^_^实际上我们一般使用雷德算法来完成倒位序,这里俺不打算写在这里,这里用一个更加简明易懂的算法,下面给出伪代码
function Reverse
o is a orginal number
r is a reversed number
o=0;
for i=0 to log(n)
if o & (1<<i)
p = p + 1<<(log(n)-i-1);
end if
end for
end function
这个算法的本质还是把从低到高进位转化为从高向低进位,非常易懂。实际上对于N较小的情况,在运算力相对较差的单片机中查表当然是不二之选,后文我们将看到这一说法不无道理。STM32F1的DSP Library中包含了一个汇编代码写成的FFT,利用查表最高实现了1024位的FFT,性能也很不错。
6.蝴蝶操作
蝴蝶操作这个真的没啥好说的,主要是复乘复加,计算旋转因子 时用欧拉公式先展开,以下是伪代码
function:Butterfly_Opreation
a is a complex number array
temp0,temp1,omega are complex numbers
omega = cos(theta)- i * sin(theta);
temp0 = a[Pre] + omega * a[Next];
temp1 = a[Pre] - omega * a[Next];
a[Pre] = temp0;
a[Next] = temp1;
end function
基2FFT的C语言实现
首先为了表示复数,我们定义了一个结构体
struct Signal
{
double Real = 0;
double Image = 0;
};
为了倒位序,我们将上面的伪代码翻译成c语言代码
void BitReverse()
{
int i, j, r;
for (j = 0; j < N; j++)
{
p = 0;
for (i = 0; i < LogN; i++)
{
if (j&(1 << i))
{
r += 1 << (LogN - i - 1);
}
}
Sig[j].Real = Input[r];
}
}
然后是蝴蝶操作的实现
void Butterfly_Opreation(int Pre, int Next, double theta)
{
double Sin_Omega = -sin(theta); //sin_omega is nagetive
double Cos_Omega = cos(theta);
Signal temp[2];
temp[0].Real = Sig[Pre].Real + Cos_Omega*Sig[Next].Real - Sin_Omega*Sig[Next].Image;
temp[0].Image = Sig[Pre].Image + Cos_Omega*Sig[Next].Image + Sin_Omega*Sig[Next].Real;
temp[1].Real = Sig[Pre].Real - Cos_Omega*Sig[Next].Real + Sin_Omega*Sig[Next].Image;
temp[1].Image = Sig[Pre].Image - Cos_Omega*Sig[Next].Image - Sin_Omega*Sig[Next].Real;
Sig[Pre].Real = temp[0].Real;
Sig[Pre].Image = temp[0].Image;
Sig[Next].Real = temp[1].Real;
Sig[Next].Image = temp[1].Image;
}
接着是FFT的主程序,非常简洁
void FFT2()
{
int i, j, k;
double theta;
for (i = 0; i < LogN; i++)
{
for (j = 0; j < (N >> (i+1)); j++)
{
for (k = 0; k < (1<<i); k++)
{
theta = 2 * Pi / (1 << (i + 1)) * k;
Butterfly_Opreation(k + j*(1 << (i + 1)), k + j*(1 << (i+1))+(1<<i), theta);
}
}
}
}
到这里看起来万事大吉了。我们输入数据测试一下
input:
8
1.0 2.0 3.0 4.0 5.0 6.0 7.0
output:
36.0000 + 0.0000i -4.0000 + 9.6569i -4.0000 + 4.0000i
-4.0000 + 1.6569i -4.0000 + 0.0000i -4.0000 - 1.6569i
-4.0000 - 4.0000i -4.0000 - 9.6569i
用MATLAB验证计算结果一致。
反思:慢在哪里
我们对
的情况进行测试。这里我们用到了@ liangbch大神制作的7种FFT,里面包含基4,基2,混合基算法,我们用其进行测试
我的电脑配置如下
intel Core i3-2310M @2.1GHz
RAM 4.0GB
Windows 10 专业版 32位
Visual Studio 2013
结果如下
程序 | CSDN fft | galois fft | Ibcfft | mixfft | ooura fft | my fft |
---|---|---|---|---|---|---|
用时 | 0.95410588 | 0.83164607 | 0.26548133 | 0.76520639 | 0.14121538 | 1.09136531 |
我们看到最快的是ooura fft,达到了0.14秒,最慢的是我的fft,用时在1.1秒,我们考虑慢在哪里?
首先我们考虑循环内每一次都需要执行的部分。在FFT中我们需要进行蝴蝶操作,其中包含的数学运算是我们优化的重点。
我们回到蝴蝶操作的函数,以下截取片断
void Butterfly_Opreation(int Pre, int Next, double theta)
{
double Sin_Omega = -sin(theta); //sin_omega is nagetive
double Cos_Omega = cos(theta);
......
}
程序中进行蝴蝶操作的次数是可以预测的,为
次,在这么多次运算中,我们看到我们每次都需要执行一次正弦余弦的求值。经验告诉我们,进行数学运算例如开方、三角反三角求职将会花去大量的时间。将这两行注释掉,我们再次运行程序,这次时间仅为0.33秒。而我们知道每一组每一次旋转因子
都将乘上自己,也就是说以幂次上升,这样耗时的三角计算就可以变成简单复乘复加,这将大大减小运算的时间,优化以后的FFT在N=2697152中用时仅0.36秒。
除此之外还有什么可以优化?有时候细枝末节的东西往往消耗了大量的时间,例如倒位序函数。我们单独测试倒位序,用时竟高达0.23秒,几乎是鲽形操作的用时两倍。单独进行蝶形运算,实际用时0.12秒,这里给京都大学的ooura FFT跪下了!!
针对倒位序算法的优化,我将在另外一篇文章中介绍。
最后也不得不提一句,在对程序进行计时的时候,务必使用Release编译。用Debug说话这不是无赖耍流氓么?
结束
从敲出文章的第一个字到结束,花了一下午的时间,中间浏览器意外崩溃,重写了伪代码那一章。实际写出FFT也不过花了十几分钟而已,值得深思的是算法背后的原理,程序的优化。总而言之程序只是个工具,数学才是硬功夫。
参考资料
[1] 《算法导论》,第四版,第30章 多项式与快速傅里叶变换,P532-P535
注:《算法导论》中默认的旋转因子
本文使用的为
,两者等价。实际在程序中,两者输出的结果虚部相反。