一步一步的无障碍理解快速傅立叶变换

时间：2015-08-05 20:29:58 阅读：2936 评论：0 收藏：0 [点我收藏+]

///////////////////////////////////////////////////////////////////////////////////////////////////////
作者：tt2767
声明：本文遵循以下协议自由转载-非商用-非衍生-保持署名|Creative Commons BY-NC-ND 3.0
查看本文更新与讨论请点击：http://blog.csdn.net/tt2767
链接被删请百度： CSDN tt2767
///////////////////////////////////////////////////////////////////////////////////////////////////////

一. 傅立叶变换的原理

任何一种变换方法都是为了我们能够更简单的解决问题，傅立叶变换亦是如此。我们平时是如何去分解一个复杂的问题呢？一个经典的方法就是把这个复杂的问题分解成为多个简单的可操作的子问题。

让我详细说来：

1.傅立叶级数

对于任何一个函数的曲线，傅立叶认为任何一个复杂的曲线都可以分解为正弦函数与余弦函数的叠加！

$s(n) = \dfrac{a_0}2+Σ_{n=1}^\infty[a_n\cos(nx)+b_n\sin(nx)]$
$\qquad——wiki（实值函数的傅立叶级数）$

函数的叠加：

而对于傅立叶级数的思想我们可以把一些复杂的曲线分解成这样：

技术分享

2.时域与频域

刚才就说道傅立叶变换是一种变换，那到底由什么变成了什么呢？其实，傅里叶变换就是把时域中的函数转换成了频域中的函数，以方便我们做数据的处理，处理完之后，用同样的原理再做一次傅立叶逆变换就得到了我们所需的结果！

域的意思是范围也是一种视角，就如高中学过的定义域一样。时域的意思就是随着时间变化的视角观察到的变化，可以理解为直角坐标系中x轴的意义为时间；频域就是频率域，意义也同样如此。

对于上图矩形波分解出来的函数们，他们在频域中是这个样子的：

这两者之间的转换就是傅立叶变换！

3.再说傅立叶级数

还记得高中老师讲过的三角函数么？

动图请看：File:Fourier series square wave

我们通常认为在时域中的基本单位是 1 ，而由上图可以看到正弦波其实是圆在直线上的投影，所以我们可以认为在频域中的基本单位是一个圆扫过的周期。

他们两个的联系就是这样的：

从左面叠加看过去就是合成的时域图像，而从右面看过去投射出的就是频域图像了。

这样的好处就是用一个个的线段把每个波表示了出来，要是需要从复杂波中去掉某一个波的时候就十分方便了。

而在更深的层次，我们可以理解为把连续的函数抽象为了单个离散的点与线的集合！

4.连续傅立叶变换

在傅立叶级数中，我们把一个在时域连续的周期函数转化成了在频域离散的非周期函数，而连续傅立叶变换是把一个在时域连续的非周期函数转化为一个频域连续的非周期函数。（说着绕嘴，一个字一个字的读）

还记得傅立叶级数中分解的那幅图么：

技术分享

被分解的曲线密度越来越密，让我想想一下，当曲线密集到一定程度的时候，一个曲线的左面无穷小处与右面无穷小处都存在另一个曲线的话，那么我们认为他们之间是相互连续的，他们的叠加就构成了连续傅立叶变换，不过这时由于无穷小处也存在需叠加的曲线，那么计算时的+号就应该变为积分符号了：

$f(x)=Σ_{n=-\infty}^{\infty}c_n e^{2\pi i(n/T)x}$
$\qquad——wiki（连续傅立叶变换公式）$

二. 离散傅立叶变换（DFT）

傅立叶变换存在4种变体，它们都是线性积分变换：傅里叶级数，连续傅里叶变换，离散傅里叶变换，离散时域傅里叶变换。

我们来看一下4种变换的图样与特性：

大家都知道，计算机系统是离散的，所以离散傅立叶变换在计算机领域得到了很大的应用与改进。
DFT的图像上不是一条连续的曲线了，而是一个一个点组成的图像，所以说它的公式是这样的：

$y_k=Σ_{j=0}^{n-1}x_j(e^{-2\pi i\dfrac{k}{n}})^j$

这里也给出逆离散傅立叶变换公式(IDFT)：

$x_j = \frac1nΣ_{k=0}^{n-1}y_k(e^{-2\pi i\dfrac{j}{n}})^k$

它最大的作用之一就是求卷积！

1.卷积

百度百科上讲的卷积很抽象，这里有一个“通俗“的例子：假如我打了你一拳，这一拳的痛感会在1个小时之内消失，但你这一个小时之内的感觉也不一样。我们记使用最轻的力打你一拳时，你在这一个小时内的疼痛感觉图像为2h（t），那么用两倍的力打你就是2h（t），3倍力3h（t）……f（n）倍力就为f（n）h（t）。如果我每秒打你一次，连续打你半分钟，那么你在一个小时零半分钟之内都会感到疼（每一次有新痛感来临时上一次痛感都会变化），那么怎么计算你在这一时间段内某一秒时的痛感呢？就是把之前的所有痛感叠加起来啊！Y(t) = 0 + …… + f（n）h（t-n），然而当你每次被打一拳的前无穷小时刻与后无穷小时刻都被打了一拳，那么我就认为打你是一件连续的行为，这时候就要把加和改写为积分了：

$Y(t) = \int_{-\infty}^{+\infty}f(n)h(t-n)dt$

而卷积定理指出：函数卷积的傅里叶变换是函数傅里叶变换后的乘积！
也就是说要求两个函数的卷积，我们就可以先把两个函数做傅立叶变换，之后算出它们的乘积，再做一次傅立叶逆变换就得出结果了！大学计算瞬间变成小学算术啊有木有！

2.点值表示法

可以粗略的理解为卷积是两个相关函数的乘积（这么说并不准确）

求多项式乘法的本质就是求卷积，像这样：
$（8x^3+6x^2+7）\times (4x^8+x^4+2x)$
可是如果单单就是这样相乘就没有什么可优化的了，可是让我们想一下多项式其实也是一个函数，每一个x值都对应一个y值，那么就如两个点确定一条直线一样，我们也能用几个点确定一个多项式，这就是点值表示法。
如果多项式的最高次项的次数为n，那么在常数已知的情况下，我们只需要知道n个不同点就能联立求得多项式了，那么我就可以用这n个点的集合表示这个多项式了。只要把点值对应相乘就可以了，复杂度为O（n）。
那么问题来了，如果单纯的代值进去，运算量会非常大，复杂度为O（n2），这时候就需要用到快速傅立叶变换了。

三. 快速傅立叶变换（FFT）

快速傅立叶变换能在O（nlogn）时间内求出卷积的结果。

1.单位复根

FFT正是利用了单位复根作为旋转因子把DFT二分才获得了如此高的效率。

定义：若 $W^n = 1$ ，则称W为1的复根，也称作单位复根。单位复根必有n过，它们是 $\sum_{k=0}^{n-1}e^{2\pi i\dfrac{k}{n}}$ ,为了方便我们定义： $W_n^k = e^{2\pi i\dfrac{k}{n}}$

它的性质为以下这些：

$W_n^n = 1$
根据欧拉公式 $e^{i\pi}+1=0$ ,
当 $W_n^k$ 中k=n时， $W_n^n = e^{2\pi i} = (e^{\pi i})^2 = 1$
相消引理
$W_{dn}^{dk} = e^{2\pi i\dfrac{k}{n}} = W_n^k$
折半引理
$（W_n^k）^2 = e^{2\pi i\dfrac{2k}{n}} = e^{2\pi i \dfrac{k}{\frac{n}{2}}} = W_\frac{n}{2}^k$

$W_n^{k +\dfrac{n}{2}} = e^{2\pi i \dfrac{2k+n}{2n}}= e^{(2\pi i\dfrac{k}{n})}*e^{\pi i} = - W_n^k$

$(W_n^{k +\dfrac{n}{2}})^2 = (-W_n^k)^2 = W_\frac{n}{2}^k$
欧拉展开
$W_n^k = e^{2\pi i \dfrac{k}{n}} = \cos(2\pi \frac{k}{n})+i * \sin(2\pi \frac{k}{n})$

Ok，基础知识补充完成，看不太懂的话无所谓，接着往下看

2.FFT原理

回顾一下DFT的公式：
$y_k=Σ_{j=0}^{n-1}x_j(e^{-2\pi i\dfrac{k}{n}})^j$

当x[]数组存放的是复数，我们可以用如下代码表示：

complex x[10],y[10]; //复数数组
for(int j = 0 ; j < n ; j++)
        y[k] += x[j] * pow（e，-2πijk/n）; //实际运算的时候用欧拉展开后的形式

注：根据公式，调整以下 $W_n^k =e^{-2\pi i\dfrac{k}{n}}$ (原理其实是一样)

首先让我们把整个公式拆掉：所有数 = 奇数 + 偶数（j = 2r + 2r+1）

$y_k=\sum_{r=0}^{\frac{n}{2}-1}x_{2r}(e^{-2\pi i\dfrac{2r}{n}})^k+Σ_{r=0}^{\frac{n}{2}-1}x_{2r+1}(e^{-2\pi i\dfrac{2r+1}{n}})^k$

再把把奇数项拆开：记x0[r] = x[2r] , x1[r] = x[2r+1]

$y_k=\sum_{j=0}^{\frac{n}{2}-1}x0_{r}((e^{-2\pi i\dfrac{r}{n}})^2)^k + e^{(-2\pi i\dfrac{r}{n})}* \sum_{j=0}^{\frac{n}{2}-1}x1_{r}((e^{-2\pi i\dfrac{r}{n}})^2)^k$

计算量就从（n-1）变成了（n/2-1）瞬间减少了一半啊！

然后我们把它抽象出来就是：

$FFT(x)[k] = DFT(x0)[k] + W_n^r * DFT(x1)[k]$

仔细观察上述公式，有没有发现x0[r] 与x1[r]里面的值为 $(W_n^r)^{2k}$ !

所以：

$FFT(x)[k] = DFT( (W_n^r)^{2k} ) + W_n^r * DFT((W_n^r)^{2k})$
$FFT(x)[k] = DFT( (W_{\frac{n}{2}}^r)^k) + W_n^r * DFT((W_{\frac{n}{2}}^r)^k)$
这里r = （0，…… ， $\frac{n}{2}-1$ ）

的确变快了但是这也就是说我们只能求出一半的值来，那后一半要怎么计算？

回顾一下单位复根的折半引理：

$W_n^{k +\dfrac{n}{2}} = e^{2\pi i \dfrac{2k+n}{2n}}= e^{(2\pi i\dfrac{k}{n})}*e^{\pi i} = - W_n^k$

$(W_n^{k +\dfrac{n}{2}})^2 = (-W_n^k)^2 = W_\frac{n}{2}^k$

我们就会发现当r=( $\frac{n}{2}$ , …… ，n-1 )时：

$FFT(x)[k] = DFT( (W_n^r)^{2k} ) - W_n^r * DFT((W_n^r)^{2k})$
$FFT(x)[k] = DFT( (W_{\frac{n}{2}}^r)^k) - W_n^r * DFT((W_{\frac{n}{2}}^r)^k)$

这两半居然是对称的！！

由于k对应着变换后的（0，……，n-1），我们就可以把FFT公式写成：

$FFT(X)[k] = DFT(x0)[r] + W_n^r * DFT(x1)[r] (r\leq\frac{n}{2}-1)$
$FFT(X)[k] = DFT(x0)[r] - W_n^r * DFT(x1)[r] (r>\frac{n}{2}-1)$

这个就是快速傅立叶变换的原理了，这样二分下去就可以在O（nlogn）复杂度下算出结果。

逆快速傅里叶变换同理，但要记得根据逆DFT，结果要除以n。

3.FFT的具体计算方法

FFT的计算方法被称为蝶形计算：

经过箭头代表乘上箭头的值，汇聚代表相加，这样正好对应了我们的公式。

以8点FFT（8个数的）为例：

它的运行顺序有些奇怪是 0,4,2,6,1,5,3,7 ，看看第一个蝶形运算，我们把它们变形：
0,2,4,6
1,3,5,7
两个两个看！是不是找到规律了，其实就是奇偶分开并且让前一半的开头与后一半的开头对齐运算。
我们把这种方法叫做倒序运算，这里的倒叙指的是每个数按位倒叙：
0 → 000 → 000 → 0
1 → 001 → 100 → 4
2 → 010 → 010 → 2
3 → 011 → 110 → 6
4 → 100 → 001 → 1
5 → 101 → 101 → 5
6 → 110 → 011 → 3
7 → 111 → 111 → 7

但是要是把每个数按位逆序的话写成程序计算量就会很大了，我们接着找规律：
仔细看，是不是只有1,4与3,6互换了位置而其他都不变？
0与n-1这两个数逆序与正序是相同的，这也限制了FFT中的n必须为 $2^n$ 。

对于中间的数（1~n-2），一组一组的看它们二进制的每一位，从低到高看，可以看出每一组是按高位值为0或者1 来决定分组的。我们设两个指针i=1，i是可能被置换的数，一直把他历遍到结尾，它的二进制是000……1；j = $\frac{n}{2}$ ，j是可能要置换的值，我们通过一系列的操作去维护它，它的二进制恰好是1……000（不信试试看）。

接着往下找，以后的j 要比i大才有意义，刚才看到了一个数加上 $\frac{n}{2}$ 也就是可以把最高位由0→1 ，同理加 $\frac{n}{4}$ 是次高位由0→1，加 $\frac{n}{8}$ 是次次高位由0→1……，so，我们每一次从右往左寻找，遇1变0，遇0跳出并且把该位变成1，这样当 $i<j$ 的时候置换他们，否则就接着找，这种方法叫做叫做二进制原地逆序置换。
让我们试一下16点的离散傅立叶变换：
00 → 0000 → 0000 → 0
01 → 0001 → 1000 → 8
02 → 0010 → 0100 → 4
03 → 0011 → 1100 → 12
04 → 0100 → 0010 → 2
05 → 0101 → 1010 → 10
06 → 0110 → 0110 → 6
07 → 0111 → 1110 → 14
08 → 1000 → 0001 → 1
09 → 1001 → 1001 → 9
10 → 1010 → 0101 → 5
11 → 1011 → 1101 → 13
12 → 1100 → 0011 → 3
13 → 1101 → 1011 → 11
14 → 1110 → 0111 → 7
15 → 1111 → 1111 → 15

其实是一样的！

这段程序可以这样写：
先给出复数结构体


struct complex
{
    double r,i;
    complex(double _r = 0.0,double _i = 0.0)
    {
        r = _r; i = _i;
    }
    complex operator +(const complex &b)
    {
        return complex(r+b.r,i+b.i);
    }
    complex operator -(const complex &b)
    {
        return complex(r-b.r,i-b.i);
    }
    complex operator *(const complex &b)
    {
        return complex(r*b.r-i*b.i,r*b.i+i*b.r);
    }
};

进行倒序运算：

void change(complex y[],int len)
{
    int i,j,k;
    for(i = 1, j = len/2;i < len-1; i++)
    {
        if(i < j)
            swap(y[i],y[j]);
        //交换互为小标反转的元素，i<j保证交换一次
        //i做正常的+1，j左反转类型的+1,始终保持i和j是反转的

        k = len/2;                  //求j的下一个倒位序
        while( j >= k)           //如果 j>=k ,表示j的最高位为1
        {
            j -= k;                    //把最高位变成0
            k /= 2;                   //k/2，比较次高位，依次类推，
                                       //逐个比较，直到某个位为0
        }
        j += k;                  //把0改为1
     }

倒序之后我们开始按照蝶形计算的方法计算：

void fft(complex y[],int len,int on) //on == 1 为FFT on==-1为 IFFT
{   
    change(y,len);                  //这3个循环非常精妙
    for(int h = 2; h <= len; h <<= 1)   //大家自己把每个数带进去手算一遍
    {                               //就明白，它是怎么进行蝶形计算的了
        complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
        for(int j = 0;j < len;j+=h)
        {
            complex w(1,0);
            for(int k = j;k < j+h/2;k++)
            {
                complex u = y[k];
                complex t = w*y[k+h/2];
                y[k] = u+t;
                y[k+h/2] = u-t;
                w = w*wn;
            }
        }
    }
    if(on == -1)                    //计算IFFT时别忘了除len
        for(int i = 0;i < len;i++)
            y[i].r /= len;           // y[i].r代表实部
}

这样我们就能快速的对函数进行离散傅立叶变换了！

4.快速傅立叶变换到底有多快？

粗略认为最高次数相同时：
DFT计算时多项式对应相乘，复杂度为：变换O（n2）+ O（n2）,频域相乘O(n)，逆运算O（n2），总复杂度为O（n2）；
FFT计算时，变换O（nlogn）+ O（nlogn）,频域相乘O(n)，逆运算O（nlogn），总复杂度为O（nlogn）；
作出图来是这样：

当1024点时FFT的计算量大约是DFT的1%！

四. 复习一下

傅立叶变换把函数从时域变换到频域
时域卷积对应着频域乘积
点值表示法把函数离散化了
离散傅立叶变换可以把点值变换到频域
点值在频域的乘积结果，逆离散傅立叶变换后得到时域函数的卷积结果
快速傅立叶变换巧妙利用了计算顺序使计算结果重复利用减少了计算次数

五. 在ACM中的应用

1.超大整数相乘

超大整数模拟手算的话复杂度是很高的，我们仔细想一想，其实超大整数相乘就是多项式相乘的一种特殊情况，所以我们可以用FFT来处理：
Eg：[HDU-1402  A * B Problem Plus](http://acm.hdu.edu.cn/showproblem.php?pid=1402)与

大数运算讲解题集中的题目

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>

using namespace std;

typedef  long long int LL;
const int INF = 99999999 ;
const long double PI = acos(0.0) * 2.0;
const int N = 200010;
complex x1[N],x2[N];
char a[N>>1],b[N>>1],res[N]={‘\0‘};
char * BigMul(char * a,char * b);

int main()
{
    while(~scanf("%s%s",a,b))
        puts(BigMul(a,b));
    return 0;
}
char * BigMul(char* a ,char*b)
{
        //保存结果的正负
        bool sign,sign_a,sign_b,flag;
        sign=sign_a=sign_b=flag=false;
        if(a[0] == ‘-‘) {a++ ; sign_a = true; }
        if(b[0] == ‘-‘) {b++ ; sign_b = true ; }
        if(sign_a != sign_b)   flag = true;

        int la = strlen(a);
        int lb = strlen(b);
        int len = 1;
        //计算点数，保证为2^n
        while(len < la*2 || len < lb*2)     len<<=1;
        //用前导0补足点数长度
        for(int i = 0;i < la;i++)
            x1[i] = complex(a[la-1-i]-‘0‘,0);
        for(int i = la;i < len;i++)
            x1[i] = complex(0,0);
        for(int i = 0;i < lb;i++)
            x2[i] = complex(b[lb-1-i]-‘0‘,0);
        for(int i = lb;i < len;i++)
            x2[i] = complex(0,0);

        FFT(x1,len,1); //FFT变换
        FFT(x2,len,1);
        for(int i = 0;i < len;i++) //频域对应相乘
            x1[i] = x1[i]*x2[i];
        FFT(x1,len,-1); //IFFT变换

        long long int num[N] = {0};
        memset(num,0,sizeof(num));
        for(int i = 0;i < len;i++) //四舍五入补足精度
            num[i] = (long long int)(x1[i].r+0.5);
        for(int i = 0;i < len;i++)//按位化为十进制数
        {
            num[i+1]+=num[i]/10;
            num[i]%=10;
        }
        len = la+lb-1; //找到最大长度
        while(num[len] <= 0 && len > 0)//删除前导0
            len--;
        int k = 0;  // 逆序转化为字符串，这步也可直接输出
        if(flag) res[k++] = ‘-‘;
        for(int i = len;i >= 0;i--)
            res[k++] = (int)num[i]+‘0‘;
        res[k] = ‘\0‘;
        if(res[0] == ‘-‘ && res[1] == ‘0‘ && res[2] ==‘\0‘)
            res[0] = ‘0‘ , res[1] = ‘\0‘;   //把-0改成0

       return res;
}

2.组合数学计数

直接看例子讲：hdu4609
题意：给出n个边的长度，求任取3条能组成三角形的概率

题中的例子a[4] = {1,3,3,4} 把每个长度的边的个数记下来 num[5]={0,1,0,2,1}，
长度为3的边有两条，为2的有0条。然后把它们两两组合起来{0,1,0,2,1}×{0,1,0,2,1}，
看出来了么，是不是就是求卷积！得出的结果num[4*2+1]={0,0,1,0,4,2,4,4,1},它的意思也是计数个数组合后长度为0 的有0 个，长度为1的有0个，……，长度为7的有4个，为8的有1个。

为了计算方便，我们把输入的数组a排序，把每一个边都当作最长边去算，这样的话虽然简便了但还是会有重复的值，把他们删去就得到最终的结果了。
具体的解释边看代码边说（已AC）：

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>

using namespace std;

#define INF 99999999
#define LL long long int
const long double PI = acos(0.0) * 2.0;
const int N = 100005;
int a[N];
LL num[N<<2],sum[N<<2];
complex x[N<<2];

LL Comb(LL x , LL y); //计算组合数
void change(complex y[],int len);
void FFT(complex y[],int len,int on);

int main()
{

    int Case,n;
    scanf("%d",&Case);
    while(Case--)
    {
        memset(num,0,sizeof(num));
        scanf("%d",&n);
        for(int i = 0 ; i < n ; i++)  //num[x]代表长度为x的边有几条
        {
            scanf("%d",&a[i]);
            num[ a[i] ]++;
        }
        sort(a,a+n);
        int len =1,l = a[n-1]+1;  // l为num的长度
        while(len < l<<1) len<<=1;  //获取len，注意为2^n;
        int i = 0;
        for(; i < l ;i++)    x[i] = complex(num[i],0);
        for(;i < len ; i++)     x[i] = complex(0,0);

        FFT(x,len,1);
        for(i = 0 ; i < len ; i++)  x[i] = x[i] * x[i];
        FFT(x,len,-1);
        for(i = 0 ; i <len ;i++) //补充精度
            num[i] = (LL) (x[i].r+0.5);

        len = a[n-1]<<1;
        //对于前n个，自己和自己的组合（a,a）应该减去
        for(i = 0 ; i <n ;i++)     num[ a[i]<<1 ]--;
        //对于所有长度组合，（a,b）与（b,a）重复了
        for(i = 1 ; i <=len ;i++)    num[i] /= 2;

        sum[0] = 0 ;
        for(i = 1 ; i <= len ; i++)   //计算前缀和
            sum[i] = sum[i-1] + num[i];

        LL res = 0;
        for(i = 0 ; i < n ; i++)
        {
            res += sum[len] - sum[ a[i] ];  //即为num[ (a[i]+1) ] +……+num[len]
            res -= Comb(i,1)*Comb(n-1-i,1); //减掉一个取小，一个取大的组合
            res -= n-1;//减掉一个取本身，另外一个取其它的组合，因为前面减过一次，所以为n-1
            res -=Comb(n-i-1,2);//减掉大于它的取两个的组合
        }
        printf("%.7f\n",(double)res/Comb(n,3));
    }
    return 0;
}
LL Comb(LL x , LL y)
{
    LL u=1,v=1,len;
    if(y==0||x==y) return 1;
    if(y==1)  return x;
    if(x < y ) return 0;

    if(y > x/2) y=len = x-y;
    else len = y;

    while(len--)
    {
        u*= x--;
        v*= y--;
    }
    return  (LL)(u/v);
}

六. 参考文献

一步一步的无障碍理解快速傅立叶变换

标签：acm fft 快速傅里叶变换

原文地址：http://blog.csdn.net/tt2767/article/details/47301849

踩

(0)

(2)

评论一句话评论（0）

分享档案

更多>

2021年07月29日 (22)
2021年07月28日 (40)
2021年07月27日 (32)
2021年07月26日 (79)
2021年07月23日 (29)
2021年07月22日 (30)
2021年07月21日 (42)
2021年07月20日 (16)
2021年07月19日 (90)
2021年07月16日 (35)

周排行