注册 登录  
 加关注
查看详情
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

n+e

NewWeb:http://trinkle.is-programmer.com/

 
 
 

日志

 
 

FFT小结  

2014-07-20 21:26:07|  分类: 算法学习 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

离散傅立叶变换在信号处理中扮演者重要的角色。利用傅立叶变换,可以实现信号在时域和频域之间的转换。
对于一个给定的长度为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。

由于所有的复数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)

注意到上式中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=w2ii
所以当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)的时间内计算傅立叶变换。

上模板
NTT:

#include<cstdio>
#include<cmath>
const int P=479<<21|1;
typedef long long ll;
ll power(ll t,int k,int mod){ll f=1;for(;k;k>>=1,t=t*t%mod)if(k&1)f=f*t%mod;return f;}
int a[1<<20],b[1<<20],n,m,k,w[2][1<<20],f;
void FFT(int x[],int k,int v)
{
int i,j,l,tmp;
for(i=j=0;i<k;i++)
{
if(i>j)tmp=x[i],x[i]=x[j],x[j]=tmp;
for(l=k>>1;(j^=l)<l;l>>=1);
}
for(i=2;i<=k;i<<=1)
for(j=0;j<k;j+=i)
for(l=0;l<i>>1;l++)
{
tmp=1LL*x[j+l+(i>>1)]*w[v][k/i*l]%P;
x[j+l+(i>>1)]=(1LL*x[j+l]-tmp+P)%P;
x[j+l]=(1LL*x[j+l]+tmp)%P;
}
}
int main(){
scanf("%d",&n);
for(i=0;i<n;i++)scanf("%d%d",&a[i],&b[i]);
for(k=1;k<n<<1;k<<=1);
w[0][0]=w[0][k]=1;j=power(3,(P-1)/k,P);
for(i=1;i<k;i++)w[0][i]=1LL*w[0][i-1]*j%P;
for(i=0;i<=k;i++)w[1][i]=w[0][k-i];
FFT(a,k,0);FFT(b,k,0);
for(i=0;i<k;i++)a[i]=1LL*a[i]*b[i]%P;
FFT(a,k,1);j=power(k,P-2,P);
for(i=0;i<2*n-1;i++)printf("%d\n",1LL*a[i]*j%P);
}

FFT:

#include<cstdio>
#include<cmath>
typedef long long ll;
typedef double ld;
const ld PI=2*asin(1);
struct P{ld x,y;};
P operator+(const P&a,const P&b){return (P){a.x+b.x,a.y+b.y};}
P operator-(const P&a,const P&b){return (P){a.x-b.x,a.y-b.y};}
P operator*(const P&a,const P&b){double d=a.x*b.x,e=a.y*b.y,f=(a.x+a.y)*(b.x+b.y);return (P){d-e,f-e-d};}

int a[1<<20],b[1<<20],n,m,k,f;ll c[1<<20];
P w[2][1<<20],x[1<<20],y[1<<20];
void FFT(P*x,int k,int v)
{
int i,j,l;P tmp;
for(i=j=0;i<k;i++)
{
if(i>j)tmp=x[i],x[i]=x[j],x[j]=tmp;
for(l=k>>1;(j^=l)<l;l>>=1);
}
for(i=2;i<=k;i<<=1)
for(j=0;j<k;j+=i)
for(l=0;l<i>>1;l++)
{
tmp=x[j+l+(i>>1)]*w[v][k/i*l];
x[j+l+(i>>1)]=x[j+l]-tmp;
x[j+l]=x[j+l]+tmp;
}
}
int main(){
scanf("%d",&n);
for(i=0;i<n;i++)scanf("%d%d",&a[i],&b[i]);
for(k=1;k<n<<1;k<<=1);
for(i=0;i<=k;i++)w[1][k-i]=w[0][i]=(P){cos(PI*2*i/k),sin(PI*2*i/k)};
for(i=0;i<k;i++)x[i]=(P){a[i],0};FFT(x,k,0);
for(i=0;i<k;i++)y[i]=(P){b[i],0};FFT(y,k,0);
for(i=0;i<k;i++)x[i]=x[i]*y[i];FFT(x,k,1);

for(i=0;i<2*n-1;i++)c[i]=(ll)(x[i].x/k+0.5);
for(i=0;i<2*n-1;i++)printf("%lld\n",c[i]);
}


注意几点:

1. 理论上有c=8,实际算了下,c大概在80左右,还是NTT,FFT就更高了。

2. NTT中注意乘爆的地方,一定要加1LL*,否则呵呵

3. FFT其实是可以撑过1048575的,只要你的PI精度足够高并且被乘数<32767。亲自测试,不服来辩。

说什么FFT精度炸翔的,应该是这样子的:

const double PI=3.14159265359;

不炸翔才怪。调了一个上午,发现跟std比对后,第i个数的误差正比于sin(2*PI*i/N),然后就在那边调常数,过了大点小点又Wa。其实= =不想say。

4. 这个模板的蝶形变换就是在模拟进位。


//rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1))

  评论这张
 
阅读(114)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2018