前言:
在学习FFT过程中看了很多博客,但发现在看博客的时候博客上的内容大多都晦涩难懂,于是乎想自己写一篇博客来记录一下自己学习的心得体会。
知其源:
先来讲讲FFT的起源:
快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少,特别是被变换的抽样点数N越多,FFT算法计算量的节省就越显著。
它是根据离散傅氏变换的奇、偶、虚、实等特性,对离散傅立叶变换的算法进行改进获得的。
(转载自百度百科)
举个例子:
上面的式子的时间复杂度为0(n^2),FFT能够把上式运算的时间复杂度将为O(nlogn)。
既然说到时间复杂度,那就来看看这里的O(logn),从O(n)到O(logn),显然用到了折半的思想,这里只是简单的提一下。
相信很多没有接触过这个类型的算法的人来说,对于怎么实现是一头雾水的,不积跬步无以至千里,下面就来聊聊FFT的前置知识。
1.1:多项式的系数表示法和点值表示法
为什么要了解这一部分的内容呢?
上面是一个由系数为{a0,a1,a2.....an}组成的多项式。
在函数图像中,这个图像可以被(n+1)个点唯一表示,也就是说我们使用{(x0,y0).....(xn,yn))这n+1个点可以表示这个图像,如果不理解的话可以类比解方程,当一个方程里面有x个未知数的时候,就需要有x个方程,根据这个想法我们可以把上面的多项式使用矩阵乘积的形式进行表示:
上面的一部分就是多项式的系数表示法,下面就来谈谈多项式的点值表示法.....
把这过多项式放到平面直角坐标系里面,然后把n+1个不同的点带入到这个多项式中就会得到n+1个不同的点(位置),那么这n+1个点就唯一确定了这个表达式,也就是说A(x)还可以表示成A(x)={(x0,A(x0)),(x1,A(x1)),......,(xn,A(xn))}。
1.2:线性卷积
在高精度乘法中,我们手推一下计算过程:
很显然,使用系数表示法进行高精度乘法的时间复杂度为0(n^2),我对这种乘法的理解是——这种乘法是横向进行的,369计算完之后计算246,最后计算123,但是呢现在如果转换一个方向的话就能改变时间复杂度,也就是我们把9看成一个部分(常数项),把66看成一个部分(x),把343看成一个部分(x^2)........那么时间复杂度就变成O(n)了。
前面提到的多项式是从常数项开始到x^n的(低位到高位),而上面我举例的乘法运算是从高位到低位的,不过变化的只是方向由x^0到x^n变成了从x^n到x^0,想象一下,当我们把运算过程模拟上面的运算列出来后左下角有一个空白的三角形,长度分别为0、1、2、3.......如果你有编程基础的话,能够想象我们在打印图形的时候遇到这样的问题的解决方法,填空,或者说往后移,那么上面的公式就实现了后移操作。
你可别觉得前面提到的多项式的系数表示法和点值表示法没有用,这里可就派上大用场了~
你瞧,这里的把常数项,x,x^2......看做一个个整体不就是点系数表示法嘛。
这样我们就可以将多项式的系数表示法转换成点值表示法来做,那时间复杂度不就降下来了嘛!
不过多项式转换成点值表示的朴素算法是O(n2)的,很遗憾,我们白高兴了一场,现在既然不能把时间复杂度从O(n^2)降到O(n),那么我们就得想其它的办法把时间复杂度降下来,毕竟时间就是金钱呐!
于是乎,傅里叶"高呼":这个简单,我会。
1.3:离散傅里叶变换
1.3.1:复数
在R范围内sqrt(x)(x>=0),但是在复数范围内,就不是这样了,定义i^2=-1,一个复数z可以表示成z=a+bi(a、b∈R,其中a为实部,b为虚部,i为虚部单位)
sqrt(-1)=sqrt(1)*i(在复数集里面可以用i表示负数的平方根)
还可以把复述看成复平面直角坐标系上面的一个点:
这个点(2,6)表示的复数就是2+6i,或者它表达的向量为(2,6),一个复数的共轭复数为虚部前面的数字*-1。
复数的运算如下:
复数不像点或者向量,它和实数一样可以进行四则运算
设两个负数分别为z1=a+bi,z2=c+di,那么:
另外,复数相乘还有一个与极角有关的小性质(这里记为①):
C++里面的的STL提供了复数模板!
头文件:#include <complex>
定义: complex<double> x;
我的实现过程如下:
struct Complex{
double x,y;//分别表示实部和虚部
Complex(double _x=0.0,double _y=0.0){
x=_x;
y=_y;
}
Complex operator -(const Complex &b)const{
return Complex(x-b.x,y-b.y);
}
Complex operator +(const Complex &b)const{
return Complex(x+b.x,y+b.y);
}
Complex operator *(const Complex &b)const{
return Complex(x*b.x-y*b.y,x*b.y+y*b.x);
}
};
FFT就是将系数表示法转化成点值表示法相乘,再由点值表示法转化为系数表示法的过程,第一个过程叫做求值(DFT),第二个过程叫做插值(IDFT)
- Duan2baka
1.3.2:NFT
傅里叶发明了一种办法:规定点值表示中的n个x为n个模长为1的复数(~~~~所以补充了复数的知识~~~~)。
对于模长为1的复数我是类比单位圆理解的~~~~
但是,傅里叶要用到的n个模长为1的复数可不是任意选择的,而是从一个圆中均匀选择的。
我们记从(1,0)开始的n个点表示为,使用复数的形式表示为:,根据①我们可以知道就是的k次方。
把带入到多项式中会得到一种特殊的点值表示,称为离散傅里叶变换!简单的说就是求值。
下面是一些性质及其证明:
性质一:
性质二:
性质三:
回到前面提到的系数表示法:
我们要将一个系数表示法的函数转换成点值表示法,那么我们就需要n+1组(不同的)数据带入到多项式里面,对于每组数据带入后的时间复杂度为O(n),那么总的时间复杂度就是O(n^2)。
这里我假设f(x)分为两部分,一个是奇数部分一个是偶数部分,表示如下:
根据前面提到的性质2我们可以再次进行化简:
在最开始我有提到折半的思想,也就是二分,想一下,上面式子里的两个部分不就可以使用二分来解决吗。
至于括号里面的,我们只需求出来,然后k次方就行了。
1.3.3:IDFT
最开始提到函数的矩阵表示,注意观察等号的两边。
我假设上面矩阵的三部分分别为A V C,那么最后的C就是我们最后要求的结果,而C是通过A和B求出来的,因为A和B位于等号的两端,所以我们可以通过C=V的逆矩阵*A进行表示。
一个矩阵乘上它的逆矩阵等于单位矩阵,根据这个性质我们可以求出逆矩阵
矩阵乘法是前一个矩阵的行中的元素乘上后一个矩阵中的列中的元素。
参考资料:
百度百科:https://baike.baidu.com/item/%E5%BF%AB%E9%80%9F%E5%82%85%E9%87%8C%E5%8F%B6%E5%8F%98%E6%8D%A2/214957?fromtitle=FFT&fromid=4766072&fr=aladdin