# 快速傅里叶变换（FFT）随笔

01 介绍FFT

FFT的过程一般为：

02快速傅里叶变换

$A(x)=a_0+a_1x+a_2x^2+\dots+a_{n-1}x^{n-1}$

$A(x)=(a_0+a_2x^2+\dots+a_{n-2}x^{n-2})+(a_1x+a_3x^3+\dots+a_{n-1}x^{n-1})$

$A(x)=(a_0+a_2x^2+\dots+a_{n-2}x^{n-2})+x\cdot(a_1+a_3x^2+\dots+a_{n-1}x^{n-2})$

$A_1(x)=a_0+a_2x+\dots+a_{n-2}x^{\tfrac{n-2}{2}}$

$A_2(x)=a_1+a_3x+\dots+a_{n-1}x^{\tfrac{n-2}{2}}$

不难发现

$A(x)=A_1(x^2)+xA_2(x^2)$

$A(\omega_{n}^k)=A_1(\omega_{n}^{2k})+\omega_{n}^{k}A_2(\omega_{n}^{2k})$

$A(\omega_{n}^k)=A_1(\omega_{\tfrac{n}{2}}^{k})+\omega_{n}^{k}A_2(\omega_{\tfrac{n}{2}}^{k})$

$A(\omega_{n}^{k}+\tfrac{n}{2})=A_1(\omega_{n}^{2k+n})+\omega_{n}^{k+\tfrac{n}{2}}A_2(\omega_{n}^{2k+n})$

$A(\omega_{n}^{k}+\tfrac{n}{2})=A_1(\omega_{n}^{2k}\cdot\omega_{n}^{n})-\omega_{n}^{k}A_2(\omega_{n}^{2k}\cdot\omega_{n}^{n})$

$A(\omega_{n}^{k}+\tfrac{n}{2})=A_1(\omega_{n}^{2k})-\omega_{n}^{k}A_2(\omega_{n}^{2k})$

$A(\omega_{n}^k)=A_1(\omega_{\tfrac{n}{2}}^{k})-\omega_{n}^{k}A_2(\omega_{\tfrac{n}{2}}^{k})$

03快速傅里叶逆变换

$c_k=\sum_{i=0}^{n-1} [\sum_{j=0}^{n-1}a_j(\omega_n^i)^j] (\omega_n^{-k})^i$

$\quad \:=\sum_{i=0}^{n-1} \sum_{j=0}^{n-1}a_j(\omega_n^j)^i (\omega_n^{-k})^i$

$\quad \:=\sum_{j=0}^{n-1} a_j \sum_{i=0}^{n-1} (\omega_n^j \omega_n^{-k})^i$

$\quad \:=\sum_{j=0}^{n-1} a_j \sum_{i=0}^{n-1} (\omega_n^{j-k})^i$

$j \neq\ k$

$\sum_{i=0}^{n-1} (\omega_n^{j-k})^i=\frac{\omega_n^0 [1-(\omega_n^{j-k})^n]}{1-\omega_n^{j-k}}$

$\qquad \qquad \quad \: \: \:=\frac{1-(\omega_n^{j-k})^n}{1-\omega_n^{j-k}}$

$\qquad \qquad \quad \: \: \:=\frac{1-(\omega_n^n)^{j-k}}{1-\omega_n^{j-k}}$

$\qquad \qquad \quad \: \: \:=\frac{1-1^{j-k}}{1-\omega_n^{j-k}}$

$\qquad \qquad \quad \: \: \:=\frac{0}{1-\omega_n^{j-k}}$

$\qquad \qquad \quad \: \: \:=0$

$j = k$

$\sum_{i=0}^{n-1} (\omega_n^{j-k})^i = n$

$c_j=a_jn$

$a_j = \frac{c_j}{n}$

04实现

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#define MAXN 4000006
using namespace std;
class complex
{
public:
complex(){}
complex(double a,double b)
{
this->a=a;
this->b=b;
}
double a,b;
}a[MAXN],b[MAXN];
complex operator+ (complex x,complex y)
{
return complex(x.a+y.a,x.b+y.b);
}
complex operator- (complex x,complex y)
{
return complex(x.a-y.a,x.b-y.b);
}
complex operator* (complex x,complex y)
{
return complex(x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a);
}
const double pi=acos(-1.0);
void FFT(int l,complex *arr,int f)
{
if(l==1) return;
int dl=l>>1;
complex a1[dl],a2[dl];
for(int i=0;i<l;i+=2)
{
a1[i>>1]=arr[i];
a2[i>>1]=arr[i+1];
}
FFT(dl,a1,f);
FFT(dl,a2,f);
complex wn=complex(cos(2.0*pi/l),sin(2.0*pi/l)*f),w=complex(1.0,0.0);
for(int i=0;i<dl;i++,w=w*wn)
{
arr[i]=a1[i]+w*a2[i];
arr[i+dl]=a1[i]-w*a2[i];
}
}
int n,m,N;
int main()
{
scanf("%d%d",&n,&m);
for(int i=0;i<=n;i++)
scanf("%lf",&a[i].a);
for(int i=0;i<=m;i++)
scanf("%lf",&b[i].a);
N=1;
while(N<n+m+1) N<<=1;
FFT(N,a,1);
FFT(N,b,1);
for(int i=0;i<N;i++)
a[i]=a[i]*b[i];
FFT(N,a,-1);
for(int i=0;i<n+m+1;i++)
printf("%d ",(int)(a[i].a/N+0.5));
puts("");
return 0;
}

import numpy as np

pi = np.arccos(-1.0)

def get_numbers():
try:
get_numbers()
return 1
except:
return 0

if not get_numbers():
return 0
if not get_numbers():
return 0

class Complex:
# 复数类

def __init__(self, a=0.0, b=0.0):
self.a = a
self.b = b

return Complex(self.a + other.a, self.b + other.b)

def __sub__(self, other):
return Complex(self.a - other.a, self.b - other.b)

def __mul__(self, other):
return Complex(self.a * other.a - self.b * other.b, self.a * other.b + self.b * other.a)

def fft(num, f, args):
if num == 1:
return
div_num = num >> 1
a1 = []
a2 = []
for i in range(0, num, 2):
a1.append(args[i])
a2.append(args[i + 1])
fft(div_num, f, a1)
fft(div_num, f, a2)
wn = Complex(np.cos(2.0 * pi / num), np.sin(2.0 * pi / num) * f)
w = Complex(1.0, 0.0)

for i in range(0, div_num):
args[i] = a1[i] + w * a2[i]
args[i + div_num] = a1[i] - w * a2[i]
w = w * wn

aa = []
bb = []
for j in range(0, n + 1):
for j in range(0, m + 1):

nn = 1
while nn < n + m + 1:
nn <<= 1

for j in range(n + 1, nn):
aa.append(Complex(0.0, 0.0))
for j in range(m + 1, nn):
bb.append(Complex(0.0, 0.0))

fft(nn, 1, aa)
fft(nn, 1, bb)

for j in range(0, nn):
aa[j] = aa[j] * bb[j]
fft(nn, -1, aa)

for j in range(0, n + m + 1):
print(int(aa[j].a / nn + 0.5), end=' ')

05结语

posted @ 2021-12-02 23:28  Roor  阅读(198)  评论(0编辑  收藏  举报