FFT,快速傅里叶变换学习笔记
来源:互联网 发布:qt 编程 编辑:程序博客网 时间:2024/06/06 03:06
FFT,快速傅里叶变换
描述问题
计算两个n阶多项式的加法,需要
多项式有两种表示法,系数表示法与点值表示法。系数表示法表示一个多项式A可以记为:
点值法表示多项式B可以记为:
可以证明,一个包含n个点的点值表示法可以唯一确定一个n-1次多项式(最高项次数为n-1)。
点值表示法得优势在于:对于点值表示的多项式,相加相乘都很方便,我们只需要将同一x对应的y相乘,就得到新多项式的点值表示。注意,因为两个n次多项式的乘积是2n次,所以我们在选点时要选2n个点。这样乘起来才会得到2n个点。所以我们快速做多项式乘法的思路就是先将系数表示的多项式转化成点值表示,相乘,在转化回系数表示。
但是一般而言,对于系数表示法的多项式A,我们要求出2n个点,需要
算法
DFT
在DFT变换中, 希望计算多项式A(x)在复数根
称向量
称向量
FFT
直接计算DFT的复杂度仍是
我们假定n为2的整次幂,若不是,在前面填0补足即可。每一步将当前的多项式A(x), 次数是2的倍数, 分成两个部分:
于是就有了:
那么我们如果能求出次数界是
处的值。
折半引理:
ωk+n/2n=ωkn∗ωn/2n=−ωn/2n
根据折半引理,这n个数其实只有
RECURSIVE_FFT(a[]){ n=a.lenth if(n==1) return a wn=e^(2*pi*i/n) w=1 a0[]=(a0,a2,a4...) a1[]=(a1,a3,a5...) y0[]=RECURSIVE_FFT(a0[]) y1[]=RECURSIVE_FFT(a1[]) for k=0 to n/2-1 y[k]=y0[k]+w*y1[k] y[k+n/2]=y0[k]-w*y1[k] w=w*wn return y}
FFT算法复杂度
因为我们将n扩展为了2的整次幂,所以每次递归子问题是原来问题规模的一半。
FFT算法的实现
递归版本的FFT实现:
typedef complex<double> Complex;const double pi=acos(-1.0);vector<Complex> recursive_fft(vector<Complex> a,int oper){ int n=a.size(); if(n==1) { return a; } Complex omgn; omgn=Complex(cos(2*pi/n*oper), sin(2*pi/n*oper)); Complex omg=Complex(1, 0); vector<Complex> a0, a1; for(int i=0;i<n;i++) { if(i%2) a1.push_back(a[i]); else a0.push_back(a[i]); } vector<Complex> y0=recursive_fft(a0, oper); vector<Complex> y1=recursive_fft(a1, oper); vector<Complex> y;y.resize(n); for(int k=0;k<n/2;k++) { Complex tmp=omg*y1[k]; y[k]=y0[k]+tmp; y[k+n/2]=y0[k]-tmp; omg=omg*omgn; } return y;}
调用之前保证vector a的size已经被置为2的整数次幂。
算法改进
时间复杂度基本不能优化,空间复杂度可以。不采用递归,而采用迭代实现,可以避免重复开数组。
首先,观察a数组的递归调用树,可以发现
可以事先将a数组排序成这种形式,然后从底层算起,每算完一层就提高一层,每层中计算
迭代版FFT
int bit_reverst(int n, int ma){ int res=0; ma--; while(ma) { res|=(n&1); res<<=1; n>>=1; ma>>=1; } res>>=1; return res;}void bit_reverse_copy(vector<Complex> &a, vector<Complex> &A){ int n=a.size(); A.resize(n); for(int k=0;k<n;k++) { int revk=bit_reverst(k, n); A[revk]=a[k]; }}vector<Complex> iterative_fft(vector<Complex> &a, int oper){ vector<Complex> A; bit_reverse_copy(a, A); int n=a.size(); for(int s=0; (1<<s)<n; s++) { int m=1<<s, m2=m*2; Complex omgm=Complex(cos(pi/m*oper), sin(pi/m*oper)); for(int k=0;k<n;k+=m2) { Complex omg=Complex(1, 0); for(int j=0;j<m;j++) { Complex t=omg*A[k+j+m]; Complex u=A[k+j]; A[k+j]=u+t; A[k+j+m]=u-t; omg=omg*omgm; } } } return A;}
bit_reverst
函数是将一个数进行二进制位的反转,比如6=110,rev(6)=011=3。
bit_reverse_copy
函数将系数数组拷贝并重排序成我们需要的顺序。
iterative_fft
实现了一个迭代版本的FFT算法。
迭代版FFT复杂度
s层for循环执行了logn次,k层for循环执行了n/m2=n/(2*m)次,j层for循环执行了m次。总复杂度:
完整代码
#include <cstdio>#include <cstdlib>#include <iostream>#include <algorithm>#include <string>#include <cstring>#include <vector>#include <cmath>#include <queue>#include <stack>#include <set>#include <map>#include <iomanip>#include <complex>using namespace std;const int MAXN=2000000+7;typedef complex<double> Complex;const double pi=acos(-1.0);vector<Complex> recursive_fft(vector<Complex> a,int oper){ int n=a.size(); if(n==1) { return a; } Complex omgn; omgn=Complex(cos(2*pi/n*oper), sin(2*pi/n*oper)); Complex omg=Complex(1, 0); vector<Complex> a0, a1; for(int i=0;i<n;i++) { if(i%2) a1.push_back(a[i]); else a0.push_back(a[i]); } vector<Complex> y0=recursive_fft(a0, oper); vector<Complex> y1=recursive_fft(a1, oper); vector<Complex> y;y.resize(n); for(int k=0;k<n/2;k++) { Complex tmp=omg*y1[k]; y[k]=y0[k]+tmp; y[k+n/2]=y0[k]-tmp; omg=omg*omgn; } return y;}int bit_reverst(int n, int ma){ int res=0; ma--; while(ma) { res|=(n&1); res<<=1; n>>=1; ma>>=1; } res>>=1; return res;}void bit_reverse_copy(vector<Complex> &a, vector<Complex> &A){ int n=a.size(); A.resize(n); for(int k=0;k<n;k++) { int revk=bit_reverst(k, n); A[revk]=a[k]; }}vector<Complex> iterative_fft(vector<Complex> &a, int oper){ vector<Complex> A; bit_reverse_copy(a, A); int n=a.size(); for(int s=0; (1<<s)<n; s++) { int m=1<<s, m2=m*2; Complex omgm=Complex(cos(pi/m*oper), sin(pi/m*oper)); for(int k=0;k<n;k+=m2) { Complex omg=Complex(1, 0); for(int j=0;j<m;j++) { Complex t=omg*A[k+j+m]; Complex u=A[k+j]; A[k+j]=u+t; A[k+j+m]=u-t; omg=omg*omgm; } } } return A;}int main(){ int n; vector<Complex> xs1; vector<Complex> xs2; while(cin>>n) { for(int i=0;i<n;i++) { int tmp;cin>>tmp; xs1.push_back(Complex(tmp, 0)); } for(int i=0;i<n;i++) { int tmp;cin>>tmp; xs2.push_back(Complex(tmp, 0)); } vector<Complex> res1; vector<Complex> res2; vector<Complex> res3; res1.resize(n*2);res2.resize(n*2);res3.resize(n*2); xs1.resize(n*2);xs2.resize(2*n); res1=iterative_fft(xs1, 1); res2=iterative_fft(xs2, 1); for(int i=0;i<res1.size();i++) { res3[i]=res1[i]*res2[i]; } vector<Complex> res;res.resize(n*2); res=iterative_fft(res3, -1); for(int i=0;i<res.size();i++) { cout<<res[i].real()/n/2<<endl; } res1=recursive_fft(xs1, 1); res2=recursive_fft(xs2, 1); for(int i=0;i<res1.size();i++) { res3[i]=res1[i]*res2[i]; } res=recursive_fft(res3, -1); for(int i=0;i<res.size();i++) { cout<<res[i].real()/n/2<<endl; } } return 0;}
- [FFT] 快速傅里叶变换学习笔记
- FFT,快速傅里叶变换学习笔记
- [学习笔记]快速傅里叶变换 FFT
- FFT(快速傅里叶变换)算法学习笔记
- 快速傅里叶变换 FFT 学习小记
- FFT—快速傅里叶变换学习小记
- 快速傅里叶变换(FFT)(ZZ)
- FFT快速傅里叶变换;
- 快速傅里叶变换(FFT)
- 关于FFT快速傅里叶变换
- FFT快速傅里叶变换
- FFT - 快速傅里叶变换
- 快速傅里叶变换(FFT)
- 快速傅里叶变换(FFT)
- GSL快速傅里叶变换FFT
- 快速傅里叶变换FFT
- 【数学】快速傅里叶变换(FFT)
- FFT 快速傅里叶变换 初探
- Django06模型(Model)
- 字符串替换
- Scala-基础语法
- QT面试题
- Echarts绘制横向柱状图(圆角+渐变)
- FFT,快速傅里叶变换学习笔记
- Are virtual interfaces supported on Quagga v0.98.3 (on Debian GNU/Linux 2.6.16)?
- Android View体系(九)自定义View
- Win7 64位系统打开 chm 文件右侧窗口显示空白的解决方法
- zabbix3.2监控Esxi5.5主机
- ycb去接水
- Halcon导出dxf文件转换至不同的图层
- “图解服务器端网络架构”小结
- CentOS7中配置NFS服务