FFT快速傅里叶变换

1、快速傅里叶变换FFTyL2013-9我们关心的问题u快速解决多项式乘法问题u衍生问题高精度乘法问题的描述u记一个多项式次数界为n的多项式A(x)u则u其中a为每一项的系数u注意最高次系数为n-110)(njjjxaxA问题的描述u两个多项式相乘u我们记一般形式为uC的次数界为A与B次数界的和u普通的时间复杂度为O(n2)()()(xBxAxCPART1中心思想2 N次单位复数根3 FFT4 演示转换思路u普通的相乘方法u提出概念：点值，插值ijjijibac02 N次单位复数根3 FFT4 演示点值u一个次数界为n的多项式A(x)的点值表达就是一个由n个点值对所组成的集合：u其中每一个x都不相

2、同，且uE.G.u多项式 的一个合法点值表达是),( ,),(),(111100nnyxyxyx)(kkxAy 23)(2xxxA)12, 5(),2 , 0(2 N次单位复数根3 FFT4 演示插值u插值运算是点值运算的逆运算u假设我们得到了一个有n个点值对点值对的点点值表达值表达u那我们可以确定唯一的一个次数界为n的多项式2 N次单位复数根3 FFT4 演示多项式乘法u我们来探究一下如何用点值与插值来完成多项式乘法u我们确定一组x，求得A与B的点值表达u那我们可以得知C的点值表达u通过插值运算，我们可以知道多项式C的系数),( ,),(),(:111100nnyxyxyxA),( ,),(

3、),(:111100nnyxyxyxB),( ,),(),(:111111000nnnyyxyyxyyxC2 N次单位复数根3 FFT4 演示注意的地方u设A与B的次数界均为nu则C的次数界为2nu则我们要找出2n个x来求点值表达u否则不可以进行插值运算2 N次单位复数根3 FFT4 演示算法流程u对于次数界均为n的多项式A与Bu1点值运算u构造长度为2n的点值表达u2逐点相乘u计算出C的点值表达u3插值运算u通过C的点值表达求出多项式C的每项系数2 N次单位复数根3 FFT4 演示时间复杂度u可以证明，若选取n个xu计算点值与插值的时间复杂度均为O(N2)u本质上没有解决时间的问题u但我们可

4、以巧妙的选择这些数来优化时间复杂度。2 N次单位复数根3 FFT4 演示PART2N次单位复数根及其性质1 点值与插值3 FFT4 演示N次单位复数根un次单位复数根是满足 的复数w。un次单位复数根根恰好有n个，对于k=0, 1, . , n-1，这些根是 。为了解释这个表达式，我们利用复数的指数形式的定义：u下一页图说明n个单位复数根均匀地分布在以复平面的原点为圆心的单位半径的圆周上。1nwnike/2)sin()cos(uiueiu1 点值与插值3 FFT4 演示N次单位复数根28w18w08w78w68w58w48wii38w,N1210nnnnnwwww次单位复数根为我们记1 点值与

5、插值3 FFT4 演示性质u我们需要N次单位复数根u我们首先来探究这些根的性质1 点值与插值3 FFT4 演示性质1主n次单位根u我们称 为主n次单位根u同时注意到，n次单位复数根都是经过旋转而得到的u每次旋转都是一定角度un次单位复数根可视为公比为主n次单位根的等比数列1nnww nininwww11 点值与插值3 FFT4 演示性质2群的性质u因为u所以u推论10nnnwwnkjnkjnknjnwwwwmod)( nknknwwknknww22)(1 点值与插值3 FFT4 演示11nnnww性质3消去引理&折半引理.u消去引理：u推论：u折半引理：如果n0为偶数，那么n次单次单位复数根位

6、复数根的平方的集合就是n/2次单位复次单位复数根数根的集合。u证明：可以知道 的平方相同。kndkdnww122/ wwnn)2/(nknknww与222222/)()(knknnnknnknnknwwwwwwknknww2/2)(1 点值与插值3 FFT4 演示性质4求和引理u求和引理：对于任意整数n1和不能被n整除的非负整数k，有u等比数列求和u所以u注意k不能是n的倍数，否则分母为0100)(njjknw011) 1 (11)(11)()(10knkknknnknnknnjjknwwwwww1 点值与插值3 FFT4 演示PART3FFT及其关键算法1 点值与插值2 N次单位复数根4 演

7、示DFT离散傅里叶变换u我们希望计算次数界为n的多项式u在n次单位复数根处的值（共n个）u接下来定义结果yuy即为a的离散傅里叶变换（DFT）u我们也可记为10)(njjjxaxA10)(njkjnjknkwawAy)(aDFTyn1 点值与插值2 N次单位复数根4 演示FFT快速傅里叶变换u大前提：n为2的整数幂（方便计算）u利用复数单位根复数根的特殊性质u我们可以在 时间内计算出uFFT利用了分治策略)lg(nnO)(aDFTn1 点值与插值2 N次单位复数根4 演示PART3.1点值运算1 点值与插值2 N次单位复数根4 演示分治策略u如何求出单个数x的函数值A(x)？u我们定义两个新多

8、项式u观察两个多项式的特点u1分别拥有奇数下标的系数与偶数下标的分别拥有奇数下标的系数与偶数下标的系数系数u2次数界变为次数界变为n/2（缩小了一半）（缩小了一半）12/12531112/224200)()(nnnnxaxaxaaxAxaxaxaaxA1 点值与插值2 N次单位复数根4 演示分治策略u对于一个数x，求A(x)u则根据上两个多项式)()()(2120 xxAxAxAy1 点值与插值2 N次单位复数根4 演示分治策略u至此我们成功的转换了问题u原问题：求一个多项式A(x)在 的函数值。u现问题：求两个多项式 在 的函数值。110,nnnnwww)(A)(10 xxA与212120)

9、( ,)( ,)(nnnnwww1 点值与插值2 N次单位复数根4 演示分治策略u现问题：求两个多项式 在 的函数值。u根据折半引理， 并不是n个不同的值，而是由n/2个n/2次单位复数根组成，而每个根恰好出现2次u于是，我们递归地对n/2的多项式A0(x)与A1(x)在n/2个n/2次单位复数根进行求值)(A)(10 xxA与212120)( ,)( ,)(nnnnwww212120)( ,)( ,)(nnnnwww1 点值与插值2 N次单位复数根4 演示程序实现1 点值与插值2 N次单位复数根4 演示我们关心的程序实现问题u1/2：规定程序递归出口u3/4/12：定义主n次单位根，更新w值

10、u5/6/7/8：定义两个多项式并递归求解u13：返回DFTu9/10/11：递归结束后的处理工作1 点值与插值2 N次单位复数根4 演示递归结束后的处理工作u10：u11：)()()(212010knknknknkknkkwAwAwwAywyy)()()()()()2/(21)2/(2021)2/(201)2/(010)2/(nknnknnknnknknnknknknknkkknknkwAwAwwAwAwwAywyywyy1 点值与插值2 N次单位复数根4 演示FFT时间复杂度u对于运行时间有以下的递归式u所以采用FFT，我们可以在O(nlgn)时间内实现点值运算（求出次数界为n的多项式在n

11、次单位复数根处的值）。)lg()()2/(2)(nnOnOnTnT1 点值与插值2 N次单位复数根4 演示PART3.2插值运算1 点值与插值2 N次单位复数根4 演示矩阵乘积u我们可以把点值运算表示成一个矩阵方程u我们把DFT写成矩阵乘积y=Vna13210)1)(1()1(3)1(21)1(3963)1(2642132113210111111111nnnnnnnnnnnnnnnnnnnnnnnnnnaaaaawwwwwwwwwwwwwwwwyyyyy1 点值与插值2 N次单位复数根4 演示逆矩阵u到此为止我们把插值运算改写成y与 的逆矩阵 的乘积1nnVyaaVy1nVnV1 点值与插值2

12、 N次单位复数根4 演示逆矩阵u定理：u证明比较冗长。u我们证明 ，其中In为n*n的单位矩阵。考虑 中(j,j)处的元素：u如果j=j，则此和为1u否则根据求和引理，此和为0./),(1nwkjVkjnn处元素为的nnnIVV11nnVV10)(101/)(/(nkjjknnkj knkjnj jnnnwwnwVV1 点值与插值2 N次单位复数根4 演示逆DFTu给定矩阵 ，可以推导出 ：u通过比较DFT的运算u我们可以看到，对FFT作以下修改就可以计算出逆DFT：u把a与y互换，用 ，再把计算结果的每个元素除以n。u于是我们也可以在O(nlgn)时间内算出逆DFT。1nV)(1yDFTn1

13、01nkkjnkjwyna10njkjnjkwaynnww 代替11 点值与插值2 N次单位复数根4 演示PART4标程演示1 点值与插值2 N次单位复数根3 FFT程序实现优化u因为像伪代码中，赋值数组是不切实际的u但我们发现，a中的应用的系数是有规律的u所以根据位移与起始点的不同u来优化FFT的实现1 点值与插值2 N次单位复数根3 FFT程序实现优化(a0,a1,a2,a3,a4,a5,a6,a7)(a0,a2,a4,a6)(a1,a3,a5,a7)(a0,a4)(a2,a6)(a1,a5)(a3,a7)(a0)(a4)(a2) (a6)(a1) (a5)(a3)(a7)1 点值与插值2

14、 N次单位复数根3 FFTu typeuNode = recordux,y:double;uend;u arr = array0.200000 of Node;u operator + (a,b:Node) c:Node;u begin c.x := a.x + b.x;c.y := a.y + b.y;end;u operator - (a,b:Node) c:Node;u begin c.x := a.x - b.x;c.y := a.y - b.y;end;u operator * (a,b:Node) c:Node;u begin c.x := a.x * b.x - a.y * b.

15、y;c.y := a.x * b.y + a.y * b.x;end;u /定义复数类型，复数运算u procedure Dft(var a:arr;s,t:longint);/a答案数组，s起始量，t位移量u beginuif (n t) = 1 then exit;uDft(a,s,t + 1);Dft(a, s + 1 (t + 1) - 1 doubeginuj := s + i (t + 1);uwt := wi (t) * aj + 1 (t + 1) := aj - wt;uend;ufor i := 0 to n t - 1 do as + i 0) and(ansi = 0)

16、 do dec(i);ufor j := i downto 0 douwrite(ansj);uwriteln;uend.u/ 感谢感谢wwt同学提供的程序同学提供的程序1 点值与插值2 N次单位复数根3 FFTu#includeu#includeu#includeu#includeu#define N 50010/uusing namespace std;uconst double pi=acos(-1);uconst complex I(0,1);uconst double eps=1e-6;uchar aaN,bbN;uint ans4*N;/char ans4*N;!ucomplex

17、a4*N,b4*N;/an-1,an-2,a1,a0uint n;uvoid initial()uu int lenA=strlen(aa),lenB=strlen(bb);u n=max(lenA,lenB);u double t=log2(n);u int tt=int(t);u if(t-tt0)tt+;u n=1(tt+1);/double lengthu int i;u for(i=0; ilenA; i+)ai=aalenA-1-i-0;u while(in)ai+=0;u for(i=0; ilenB; i+)bi=bblenB-1-i-0;u while(in)bi+=0;u1

18、 点值与插值2 N次单位复数根3 FFTuvoid bitReverse(complex * a)uu for(int i=1; in-1; i+)u u int j=0;u for(int k=1,tmp=i; kn; j=(j1)|(tmp&1),k=1);u if(ji)swap(ai,aj);u uuvoid iterative_FFT(complex *a,int sig)uu bitReverse(a);u for(int m=2; m=n; m1;u for(int i=0; imh; i+)u u complex wi=exp(i*sig*pi/mh*I);u for(int

19、j=i; jn; j+=m)u u int k=j+mh;u complex u=aj,t=wi*ak;u aj=u+t;u ak=u-t;u u u u if(sig=-1)for(int i=0; in; i+)ai/=n;u1 点值与插值2 N次单位复数根3 FFTuvoid convolution()uu iterative_FFT(a,1);u iterative_FFT(b,1);u for(int i=0; in; i+)ai*=bi;/a*bu iterative_FFT(a,-1);uuvoid output()uu int k=0;u ans0=0;u for(int i=0; i=0)coutansk