您的位置：首页 > 其它

初探 FFT/DFT

2015-02-07 08:29 197 查看

有用的学习链接&书籍

傅立叶变化-维基百科

离散傅立叶变化-维基百科·长整数与多项式乘法

维基百科看英文的更多内容&有趣的图

快速傅立叶变化-百度百科，注意其中的图！

组合数学（第4版） Page 287~291（讲得挺详细）

FFT/DFT是个什么东西？

说实话，我也不知道，不过根据维基百科上面的图，就可以略窥一二了：




傅里叶变换将函数的时域（红色）与频域（蓝色）相关联。频谱中的不同成分频率在频域中以峰值形式表示。

——“Fourier transform time and frequency domains (small)”作者Lucas V. Barbosa - 自己的作品。来自维基共享资源 - http://commons.wikimedia.org/wiki/File:Fourier_transform_time_and_frequency_domains_(small).gif#mediaviewer/File:Fourier_transform_time_and_frequency_domains_(small).gif根据公有领域授权


（注，下面的话都是我乱编的，不具有很强的科学性）

它的意思是，红色的周期函数可以通过若干的波（图中蓝色部分）来表示。

FFT的作用就是通过红色的那一部分求出蓝色那一部分。在oi中几乎都是离散函数（你可以认为数组的下标就是定义域），所以oi中用的几乎都是DFT。

它的最基本函数是(摘自维基百科)：



它们的过程十分相似，我们只要会了DFT，就能套入IDFT了。

注：\(e^{i\theta} = cos \theta + i sin \theta\)，没错，这是复数运算，c++中有
complex<double>
。

如果我们按照上面的公式计算，时间复杂度为\(O(n^2)\)。

加快计算

令\(\omega = e^{2 \pi i}\)，\(\omega_n^p= e^{ \frac {2 \pi i \cdot p} {n}}\)。

组合数学一书中通过举例子，发现了计算规律（通过化简式子，比如\(\omega^6_4=\omega^2_4\)），发现规律的过程就不展开了。

下面的图片摘自百度百科：



最左边的\(8\)个黑点是原数列，最右边的\(8\)个是输出，没错，输出的下标是有顺序的，但是原数列是无序的（但有规律，看第\(4\)行，左边是\(6=1100_2\)，右边是\(3=0011_2\)，二进制的字符顺序恰好是反的）。

计算过程可以从代码中看出：

C表示，\(1\)为DFT，\(-1\)为IDFT，
MyType
即
complex<double>
，
cnt
是数列大小，
unit
是\(\omega_{cnt}\)。

void FFT(MyType A[], int cnt, int C) {
if (cnt == 1) return;

MyType* L = newMemory(cnt >> 1);
MyType* R = newMemory(cnt >> 1);

for (int i = 0; i < cnt; i += 2) {
L[i >> 1] = A[i];
R[i >> 1] = A[i + 1];
}

FFT(L, cnt >> 1, C);
FFT(R, cnt >> 1, C);

MyType w(1, 0);
MyType unit(cos(C * PI * 2 / cnt), sin(PI * 2 / cnt * C));
//also unit(cos(    PI * 2 / cnt), sin(PI * 2 / cnt    ) * C)

for (int i = 0; i < cnt >> 1; i ++) {
A[i] = L[i] + w * R[i];
A[i + (cnt >> 1)] = L[i] - w * R[i];
w *= unit;
}
}

时间复杂度\(O(n * logn)\)。

oi里用DFT来求两个多项式的乘积

算法（此处仅有算法，我不知道它的正确性如何证明）：

对于两个多项式\(A\)和\(B\)，我们要算出结果\(C\)：

对\(A\)进行DFT得到\(\hat{A}\)

对\(B\)进行DFT得到\(\hat{B}\)

令\(\hat C_i = \hat A_i * \hat B_i\)

对\(\hat C\)进行逆DFT（IDFT）得到\(C\)

后

这篇随笔写得不清不楚的，其实我是想我自己以后看了自己明白。

听说可以写模意义下的DFT（没有误差），找个时间看看？？

lyl所说的“总结”图：



附上求两个多项式乘积的代码，注complex的两个值调用分别是real和imag，看这里：

#include <cstdio>
#include <cmath>
#include <algorithm>
#include <iostream>
#include <complex>
#include <assert.h>
using namespace std;

const int MAXLOGN = 18;
const int MAXN = 1 << MAXLOGN;
const int MAXMEMORY = MAXLOGN * MAXN;
const double PI = acos(-1.);

typedef complex<double> MyType;

int n;

MyType mem[MAXMEMORY];
MyType* cur_mem;

MyType* newMemory(int size) {
cur_mem += size;
return cur_mem - size;
}

void FFT(MyType A[], int cnt, int C) {
if (cnt == 1) return;

MyType* L = newMemory(cnt >> 1);
MyType* R = newMemory(cnt >> 1);

for (int i = 0; i < cnt; i += 2) {
L[i >> 1] = A[i];
R[i >> 1] = A[i + 1];
}

FFT(L, cnt >> 1, C);
FFT(R, cnt >> 1, C);

MyType w(1, 0);
MyType unit(cos(C * PI * 2 / cnt), sin(PI * 2 / cnt * C));

for (int i = 0; i < cnt >> 1; i ++) {
A[i] = L[i] + w * R[i];
A[i + (cnt >> 1)] = L[i] - w * R[i];
w *= unit;
}
}

int main() {
static MyType A[MAXN], B[MAXN];
int cntA, cntB;
scanf("%d%d", &cntA, &cntB);
cntA ++, cntB ++;
for (int i = 0; i < cntA; i ++)
scanf("%lf", &A[i].real());
for (int i = 0; i < cntB; i ++)
scanf("%lf", &B[i].real());

n = 1;
while (n < cntA + cntB - 1) n <<= 1;

cur_mem = mem; FFT(A, n, 1);
cur_mem = mem; FFT(B, n, 1);
for (int i = 0; i < n; i ++) A[i] *= B[i];
cur_mem = mem; FFT(A, n, -1);

for (int i = 0; i < cntA + cntB - 1; i ++)
printf("%d ", (int) floor(A[i].real() / n + 0.1));
printf("\n");
return 0;
}
内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 点击这里给我发消息
标签： 
相关文章推荐
新的分享
章节导航