很魔法,总结写在最前面,因为这是我觉得最关键的内容。
本人专业是自动化,离不开信号与系统这门课,提到信号与系统,自然是离不开大名鼎鼎的傅立叶,傅立叶级数、傅立叶变换都是常用的计算工具。
同时,我也搞搞算法竞赛,众所周知,FFT(快速傅立叶变换)是一种计算多项式乘法的算法。
虽然我都接触过这两个东西,但是快速傅立叶变换和傅立叶变换,有什么关系吗?
当然有,在此之前,你得先知道有个东西叫离散傅立叶变换(DFT),FFT的最早提出其实是为了快速计算出DFT。也就是说,FFT是一种算法,它的命名是因为它最早被用于求解离散傅立叶变换,而同时它也可以解决一些其他问题。
然后我现在总结一下这些东西。
关于拉普拉斯变换和Z变换写在了另一篇博客中。
傅立叶变换是针对周期函数而言的。数学家经常这么想,能不能用一些简单的函数累加去组合出(逼近)一个复杂的函数,其中最有名的肯定是泰勒公式。
傅立叶级数也是类似的想法,周期函数能不能分解成很多简单的函数和。和出来的函数是周期性的,那自然想到组成函数也是周期性的,最简单的周期函数是什么呢?当然是三角函数(简单在容易积分和求导,这个简单还包括正交性,乘积的积分等于0)。
然后就有了傅立叶级数:
推导可以看这个博客,也包括了后面傅立叶变换啥的
然后再增加一点,非周期信号也可以计算傅立叶级数。就是把函数当成:存在定义域的周期函数,即定义域内为整个函数的一个周期。
公式就变成了:
傅立叶级数是针对周期信号的,那如果是定义在整个周期上的非周期信号,就不可以用周期信号去逼近了吗?答案是可以,但是是连续的,既然是连续的,那就不能写成求和的形式,而是积分的形式。公式如下,包括逆变换公式:
谈谈理解,F(w)是频谱,什么是频谱,可以理解为一个个不同频率的周期函数,它们可以组成整个函数,但是这些函数的个数是无穷多的,所以函数的横坐标实际表达的是在某个频率附近函数的密度。如果你学过概率论,知道概率密度函数,应该可以比较轻松的理解。
然后,我们再反过来思考,对于周期信号,傅立叶变换是什么样的?如果你去看书,会发现频谱图是一个个离散的冲激信号(不懂自己查),也就是说这个地方频率密度无限高,换言之就是这个频率。这其实非常符合我们用傅立叶变换去逼近函数,用于逼近的周期函数的频率正好是离散的。
PS:提一个很简单的应用,稍微给大家建立一点实际应用的感觉。两个同时说话,声调不同,我们一开始采集到的声轨是时间关于声音(音量?振动速度?原谅我物理没学好)的函数。我们把它傅立叶变换,由于两个人声带振动频率不同,变换后的函数会在两个不同的频率处有两个峰,我们消除其中一个并拟变换回来,就消除了其中一个人的声音。类似的应用还有很多,比如说消除图片中的高频噪声等等。
离散傅里叶级数是针对周期离散时间信号而言的。
离散序列傅里叶级数(变换)之后是周期性的,这很关键。因为下一个周期又会重复某个(4pi和2pi转起来是一样的),才能抵消掉多余的部分,不然一圈内的频率,肯定组成的原函数是连续的。(一种感觉罢了 )
推出来是这个样子的。
先用单位冲激序列对函数采样,就变成了离散时间序列。然后使用傅里叶变换可以推出来DTFT的公式。(交换积分次序,利用冲激函数的积分求)
对于计算机来说,无穷是无法求解的。就像电脑算积分,其实是把函数分成非常多小条矩形,加起来的面积就是积分值,离散傅立叶变换也是类似的东西。
对于频谱密度函数(频谱F(w)),可以感受一下,如果每隔一小段,用其中某个值代替频率密度,函数就变成离散的了,然后可以用这些频率和其函数值作为系数去组成(累加)原函数,这不就是所谓的离散傅立叶变换吗。
这是DFT的公式:
看傅立叶级数的指数形式:
有没有觉得很像,注意:
标准化的傅立叶变换T趋向于N(隔一个点采一下),傅立叶级数指数形式中括号内就是傅立叶变换,而用离散傅立叶变换代替,T换成N,就变成了傅立叶逆变换。我并没有推过公式,但是这证明的感受是对的。可以说,离散傅立叶变换将傅里叶级数推广到了所有函数。
然后这个公式还可以写成矩阵变换的形式,记住这个
W
n
t
W_n^t
Wnt,这是快速傅里叶变换的关键,快速傅里叶变换指的就是快速求出所有的
W
n
t
W_n^t
Wnt(nlogn复杂度)。
然后还可以看一下这篇链接,不从傅里叶变换得到离散傅里叶变换,更容易理解也更好记。
首先是多项式的存储方式,一是系数存储法,那么在做多项式乘法的时候,显然会产生
n
2
n^2
n2的系数,然后再进行合并,那太慢了。二是使用点值存储法,n个点可以确定一个n阶多项式,并且点与点的对应相乘是O(n)的,非常快。要解决的就是如何将系数多项式转变为点值,以及如何将点值多项式转回去。
首先,将系数多项式转为转为点值,就是带n个点进去算,但是计算一个点的n次方是很耗时的,我们需要选择合适的点,加速这个过程。
我们重新来看看
ω
n
k
\omega^k_n
ωnk,将复平面单位圆n等分,
ω
n
0
\omega^0_n
ωn0代表点(1,0),
ω
n
k
\omega^k_n
ωnk代表绕单位元逆时针k号点。
显然
ω
n
k
=
ω
2
n
2
k
\omega^k_n=\omega^{2k}_{2n}
ωnk=ω2n2k,
ω
n
k
=
−
ω
n
k
+
n
/
2
\omega^k_n=-\omega^{k+n/2}_{n}
ωnk=−ωnk+n/2
推导:
含义解读:一个长度为n的序列,A[k]代表带入单位圆k分点对应的函数值,这个A[k]可以根据奇偶分成两部分A1[n/2]和A2[n/2],并且已知A1和A2可以O(n)求出A。并且A1和A2保持了原来的形式,可以递归求解。
优化:递归常数太大,考虑迭代。假设我们知道最后一层每个点的位置,我们就可以一层层地往上迭代了。另外,迭代过程只需要一维数组,因为每次都是两个点刷两个点。
观察得出:
最终每个元素所在位置就是其二进制翻转。
code:
for(int i=0;i<limit;i++) r[i]= ( r[i>>1]>>1 )| ( (i&1)<<(l-1) ) ; for(int i=0;i<limit;i++) if(i<r[i]) swap(A[i],A[r[i]]);//求出要迭代的序列,i<r[i]保证只换一次
PS:FFT一定要n是2的整数次幂,如果不是,将n增大到2的幂次,高次项系数为0。
然后是逆变换(IFFT),经过一套推算之后,在分治的过程中乘上单位根的共轭复数(
ω
n
k
\omega^k_n
ωnk的共轭复数),分治完的每一项除以n即为原多项式的每一项系数。
fft和ifft共用一个函数,ifft处理完后所有项要/n,
注意要这么写c[i]=(int)(a[i].real()/n+0.5);因为存在精度问题,1.99999999这种
本质是一个东西,应该说是多项式乘法被转化为了和DFT相同的形式。
把上面DFT公式用欧拉公式展开,就变了上面的形式,不信你试试看。