刷清橙OJ--A1084.快速傅里叶变换

问题：

A1084. 快速傅立叶变换

时间限制：1.0s   内存限制：512.0MB  

总提交次数：1414   AC次数：177   平均分：43.67

问题描述

　　离散傅立叶变换在信号处理中扮演者重要的角色。利用傅立叶变换，可以实现信号在时域和频域之间的转换。

　　对于一个给定的长度为n=2m (m为整数) 的复数序列X0, X1, …, Xn-1，离散傅立叶变换将得到另一个长度为n的复数序列Y0, Y1, …, Yn-1。其中

　　Yi=X0+X1wi+ X2w2i+ X3w3i+…+ Xn-1w(n-1)i

　　其中w=e2πI/n=cos(2π/n)+I sin(2π/n)，称为旋转因子，其中I为虚数单位，I2= –1。

　　给定输入序列X，请输出傅立叶变换后的输出序列。

　　由于所有的复数C都可以表示成C=a+Ib的形式，即由实部和虚部的和表示，所以C可以用一个二元组(a, b)来表示，用这种方法w可表示为(cos(2π/n), sin(2π/n))。复数的计算规则如下：

　　(a1, b1)+(a2, b2)=(a1+a2, b1+b2)

　　(a1, b1)(a2, b2)=(a1, b1)*(a2, b2)=(a1*a2-b1*b2, a1*b2+a2*b1)

　　对于本题，你可以按照上面的公式直接计算，也可以按下面的方法计算。使用上面的公式的复杂度为O(n2)，可以得到一半的分数，而下面的方法的复杂度为O(n log n)，可以得到全部的分数。如果要使用上面的公式直接计算，请跳过下面的算法描述部分。

算法描述

　　注意到上式中w=e2πI/n=cos(2π/n)+I sin(2π/n)，所以wn+k=wk，这个公式将在下面的计算用用到。

　　对上面的公式进行变形，得到：

　　Yi

　　=X0 + X1wi + X2w2i + X3w3i +…+ Xn-1w(n-1)i

　　=X0 + X2w2i + X4w4i +…+ Xn-2w(n-2)i + wi(X1 + X3w2i + X5w4i +…+ Xn-1w(n-2)i)

　　=(X0 + X2φi + X4φ2i +…+ Xn-2φ(n/2-1)i) + wi(X1 + X3φi + X5φ2i +…+ Xn-1φ(n/2-1)i)

　　其中φ=w2。利用这个公式可得

　　Yi+n/2=(X0 + X2φi+n/2 + X4φ2(i+n/2) +…+ Xn-2φ(n/2-1) (i+n/2)) + wi(X1 + X3φ(i+n/2) + X5φ2(i+n/2) +…+ Xn-1φ(n/2-1) (i+n/2))

　　其中φi+n/2=w2i+n=w2i=φi。

　　所以当i<n/2时，令pi=X0 + X2φi + X4φ2i +…+ Xn-2φ(n/2-1)i，qi= X1 + X3φi + X5φ2i +…+ Xn-1φ(n/2-1)i，则Yi=pi+wiqi，Yi+n/2= pi+wi+n/2qi。

　　利用上面的公式，我们可以得到一种快速计算旋转因子为w的傅立叶变换的方法如下（快速傅立叶变换）：

　　算法Y=Fourier(X, n, w)

　　参数：X={Xi}为待变换的序列，n为序列的长度（2的整数次幂），w为旋转因子

　　结果：X的傅立叶变换Y={Yi}

　　1. 计算{X0, X2, X4, …, Xn-2}在旋转因子为φ=w2时的傅立叶变换序列{pi}

　　2. 计算{ X1, X3, X5, …, Xn-1}在旋转因子为φ=w2时的傅立叶变换序列{qi}

　　3. 对于0<=i<n，计算Yi=pi+wiqi。其中w0=(1, 0)，wi=wi-1*w，你要设置一个临时变量进行累乘而不能每次重新计算wi。

　　使用这种方法，即可在O(n log n)的时间内计算傅立叶变换。

输入格式

　　输入的第一行包含一个整数n，表示待变换的序列的长度，n是2的整数次幂。n<=16384。

　　接下来n行，每行2个实数a, b表示序列中的一个元素为(a, b)。

输出格式

　　输出n行，每行两个实数，表示经过变换后的序列。为了防止运算时的误差影响结果的输出，请将结果中的每一个实数除以n后保留两位小数。

样例输入

4

1.0 0.0

1.0 0.0

0.0 0.0

0.0 0.0

样例输出

0.50 0.00

0.25 0.25

0.00 0.00

0.25 -0.25

代码：

#include <cstdio>#include <cstdlib>#include <math.h>using namespace std;#define pi 3.1415926535897932384626433832795int N;struct vec{double x,y;void operator = (const vec a){if(fabs(a.x) > 1e-20) x = a.x;else x = 0;if(fabs(a.y) > 1e-20) y = a.y;else y = 0;}vec operator +(const vec a){vec tmp;tmp.x = x + a.x;tmp.y = y + a.y;return tmp;}vec operator *(const vec a){vec tmp;tmp.x = x*a.x - y*a.y;tmp.y = x*a.y + a.x*y;return tmp;}};vec* fft(vec *X,int n,vec w){if(n == 2){vec *y = (vec*)malloc(n* sizeof(vec));y[0] = X[0] + X[1];y[1] = X[0] + w*X[1];return y;}vec *odd = (vec*)malloc(n/2*sizeof(vec));vec *even = (vec*)malloc(n/2*sizeof(vec));int k1 = 0,k2 = 0;for(int i = 0;i < n;i += 2){odd[k1++] = X[i + 1];even[k2++] = X[i];}vec *p = fft(even,n/2,w*w);vec *q = fft(odd,n/2,w*w);free(even);free(odd);vec *y = (vec*)malloc(n* sizeof(vec));vec wi;wi.x = 1;wi.y = 0;for(int i = 0;i < n/2;i++){y[i] = p[i] + wi*q[i];wi = wi*w;}for(int i = 0;i < n /2;i++){y[i + n / 2] = p[i] + wi*q[i];wi = wi*w;}free(p);free(q);return y;}int main(){vec X[16400];scanf("%d",&N);for(int i = 0;i < N;i++){scanf("%lf %lf",&X[i].x,&X[i].y);}vec w;w.x = cos(2*pi/N);w.y = sin(2*pi/N);vec *Y = fft(X,N,w);double m = N;for(int i = 0;i < N;i++)printf("%.2lf %.2lf\n",Y[i].x/m,Y[i].y/m);return 0;}

个人想法：代码来自试题讨论