当前位置：首页 > news >正文

FFT剖析

news 2025/8/10 19:04:57

快速傅里叶变换 (fast Fourier transform)

xn={x0,x1,…xn-1} (num:N)

旋转因子系数：

d=2pik/N

旋转因子

wk(n)=(cos(dn)+isin(dn)) n=[0,N-1]
y(k)= sum(x(n)wk(n),0,N-1)
y(k)={y(0),y(1),…y(N-1)} 傅里叶级数
x(n)wk(n)的级数是：
1.d=2pik/N 这个系数决定转动圈数，不懂看第二个再说
2.y(k)=x(0)wk(0)+x(1)wk(2)+…+x(n-1)wk(n-1)
旋转因子：cos(t)+isin(t)
t=2pikn/N
T=2pi就是一圈，那么可以得出步进量：2pi*k/N, 圈数:k
在一圈中：wk(n)的单位圆上t<=pi/2，是他的展开区 t<=pi是他的对称区

#include <bits/stdc++.h>
#define fu(x , y , z) for(int x = y ; x <= z ; x ++)
#define fd(x , y , z) for(int x = y ; x >= z ; x --)
#define LL long long
using namespace std;
const int N = 4e6 + 5;
const double pi = acos (-1.0);
struct node {double x , y;
} a[N] , b[N];
int n , m , len = 1 , r[N] , l;
node operator + (node a, node b) { return (node){a.x + b.x , a.y + b.y};}
node operator - (node a, node b) { return (node){a.x - b.x , a.y - b.y};}
node operator * (node a, node b) { return (node){a.x * b.x - a.y * b.y , a.x * b.y + a.y * b.x};}
int read () {int val = 0 , fu = 1;char ch = getchar ();while (ch < '0' || ch > '9') {if (ch == '-') fu = -1;ch = getchar ();}while (ch >= '0' && ch <= '9') {val = val * 10 + (ch - '0');ch = getchar ();}return val * fu;
}
void fft (node *A , int inv) {for (int i = 0 ; i < len ; i ++)if (i < r[i]) swap (A[i] , A[r[i]]);for (int mid = 1 ; mid < len ; mid <<= 1) {  //mid={1,2,4,8,16,32,64,...}node wn = (node){cos (1.0 * pi / mid) , inv * sin (1.0 * pi / mid)};for (int R = mid << 1 , j = 0 ; j < len ; j += R) {node w = (node){1 , 0};for (int k = 0 ; k < mid ; k ++ , w = w * wn) {node x = A[j + k] , y = w * A[j + mid + k];A[j + k] = x + y;A[j + mid + k] = x - y;}}}
}
int main () {n = read () , m = read ();fu (i , 0 , n) a[i].x = read ();fu (i , 0 , m) b[i].x = read ();while (len <= n + m) len <<= 1 , l ++;for (int i = 0 ; i < len ; i ++)r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));fft (a , 1);fft (b , 1);fu (i , 0 , len) a[i] = a[i] * b[i];fft (a , -1);// for (int i = 0 ; i <= n + m ; i ++) cout << a[i].x << " " << a[i].y << "\n";fu (i , 0 , n + m) printf ("%d " , (int)(a[i].x / len + 0.5));return 0;
}

计算2^n的逆排长度：2的6次方：64=0b100-0000；逆排长度0b11-1111

while (len <= n + m) len <<= 1 , l ++;
产生逆排：0,32，
for (int i = 0 ; i < len ; i ++)
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
逆排交换数据
for (int i = 0 ; i < len ; i ++)
if (i < r[i])
swap (A[i] , A[r[i]]);

#include <bits/stdc++.h>
#define fu(x , y , z) for(int x = y ; x <= z ; x ++)
using namespace std;
const int N = 4e6 + 5;
const double pi = acos (-1.0);
int n , m1 , m2 , rev[N];
complex<double> a[N] , b[N];
void fft (complex<double> *a , int type) {fu (i , 0 , n - 1) if (i < rev[i]) swap (a[i] , a[rev[i]]);for (int j = 1 ; j < n ; j <<= 1) {complex<double> W(cos (pi / j) , sin (pi / j) * type);for (int k = 0 ; k < n ; k += (j << 1)) {complex<double> w(1.0 , 0.0);fu (i , 0 , j - 1) {complex<double> ye , yo;ye = a[i + k] , yo = a[i + j + k] * w;a[i + k] = ye + yo;a[i + k] = ye + yo;a[i + j + k] = ye - yo;w *= W;}}}
}
int main () {scanf ("%d%d" , &m1 , &m2);fu (i , 0 , m1) cin >> a[i];fu (i , 0 , m2) cin >> b[i];n = m1 + m2;int t = 0;while (n >= (1 << t)) t ++;n = (1 << t);fu (i , 0 , n - 1) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (n >> 1) : 0);fft (a , 1) , fft (b , 1);fu (i , 0 , n) a[i] *= b[i];fft (a , -1);fu (i , 0 , m1 + m2) printf ("%d " , (int)(a[i].real() / (double)n + 0.5));return 0;
}

f0=(a0,a1,…an-2)
f1=(a1,a3,…an-1)
f(wnk)=f0(wnk/2)+wnkf1(wkn/2)
f(wn(k+n/2)=f0(wnk/2)-wnkf1(wkn/2)

#include<cstdio>2 #include<iostream>3 #include<cmath>4 #include<cstring>5 #include<algorithm>6 #include<cstdlib>7 using namespace std;8 const int mod=1e9+7;9 const double pi=acos(-1);
10 struct cn
11 {
12     double x,y;
13     cn (double x=0,double y=0):x(x),y(y) {}
14 }a[300005],b[300005],c[300005];
15 cn operator + (const cn &a,const cn &b) {return cn(a.x+b.x,a.y+b.y);}
16 cn operator - (const cn &a,const cn &b) {return cn(a.x-b.x,a.y-b.y);}
17 cn operator * (const cn &a,const cn &b) {return cn(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
18 void fft(cn a[],int n,int l,int f)                                   
19 {                                                                                
20     int rev[n+5];
21     rev[0]=0;
22     for (int i=1; i<n; i++){
23         rev[i]=(rev[i>>1]>>1)|((i&1)<<l-1);
24         if (i<rev[i]) swap(a[i],a[rev[i]]);
25     }
26     for (int i=1; i<n; i<<=1){
27         cn wi(cos(pi/i),f*sin(pi/i));
28         for (int j=0; j<n; j+=i*2){
29             cn w(1,0);
30             for (int k=0; k<i; k++){
31                 cn x=a[j+k],y=w*a[j+k+i];
32                 a[j+k]=x+y;
33                 a[j+k+i]=x-y;
34                 w=w*wi;
35             }
36         }
37     }
38     if (f==-1)
39         for (int i=0; i<n; i++){
40             a[i].x/=n; a[i].y/=n;
41         }
42 }
43 int main()
44 {
45     int n,m;
46     scanf("%d%d",&n,&m); n++; m++;
47     for (int i=0; i<n; i++) scanf("%lf",&a[i].x);
48     for (int i=0; i<m; i++) scanf("%lf",&b[i].x);
49     int l=0,N=1;
50     while (N<n+m-1) N<<=1,l++;
51     fft(a,N,l,1);
52     fft(b,N,l,1);
53     for (int i=0; i<N; i++) c[i]=a[i]*b[i];
54     fft(c,N,l,-1);
55     for (int i=0; i<n+m-1; i++) printf("%d ",(int)(c[i].x+0.5));
56     return 0;
57 }