[笔记]FFT算法

原 作者 殷健翔 | 发表于 | 分类于 2016级

前言

对于学通信的人来说，在学到数字信号处理时都会学到一个东东，叫做快速傅里叶变换(Fast Fourier Transform,简称FFT)。这东西真的挺有用的，但是只要有那么一点用的东西，就是特别难的。(现在也有很多不完整的地方，以后再补充~)

什么是FFT

FFT，即为快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步

FFT是一种用来计算DFT（离散傅里叶变换）和IDFT（离散傅里叶反变换）的一种快速算法。这种算法运用了一种高深的数学方式、把原来复杂度为 O ( n 2 )
的朴素多项式乘法转化为了 O ( n l o g n ) 的算法。

首先上DFT的公式，下面将由这个公式一步步推出FFT的算法。

X [ k ] = ∑ n = 0 N − 1 x [ n ] ∗ e − i 2 π k N n

它其实就是针对离散的采样点 x[n]的傅里叶变换，其意义还是将信号从时域转换到频域。那么这里的 X[k] 表示的就是频域信号了，下标k表示信号的频率， X[k] 的值就是该频率信号的幅值。这里可以看到，当k固定的时候，对应于频率为k的信号幅值是与 x[n]有关的一个级数的和，因此每一个频域的点实际上包含了所有时域点的信息。

有两点需要注意，一是因为 x[n]的长度为 N ，因此 DFT 后的频域信号 X[k] 的长度也是 N；二是为了进行FFT变换，N的取值一般是2的整数次幂。

将下标n分解为奇偶两部分，然后分别求和

X [ k ] = ∑ r = 0 N 2 − 1 x [ 2 r ] ∗ e − i 2 π k N 2 r + ∑ r = 0 N 2 − 1 x [ 2 r + 1 ] ∗ e − i 2 π k N ( 2 r + 1 )

这里进行一次代换，用 x 0 [ n ] 表示 x[n] 中的偶数项，用 x 1 [ n ] 表示 x[n] 中的奇数项

X [ k ] = ∑ n = 0 N 2 − 1 x 0 [ n ] ∗ e − i 2 π k N / 2 n + ∑ n = 0 N 2 − 1 x 1 [ n ] ∗ e − i 2 π k N / 2 n − i 2 π k N

第二项求和中， e − i 2 π k N 与 n 无关，因此可以单独提出

对 x 0 [ n ] 进行DFT
X [ k ] = ∑ n = 0 N 2 − 1 x 0 [ n ] e − i 2 π k N / 2 n

对 x 1 [ n ] 进行DFT
X [ k ] = e − i 2 π k N + ∑ n = 0 N 2 − 1 x 1 [ n ] e − i 2 π k N / 2 n

利用欧拉公式，将： e − i 2 π k N = cos ⁡ 2 π k N − i sin ⁡ 2 π k N 代入：

X [ k ] = D F T ( x 0 ) + ( cos ⁡ 2 π k N − i sin ⁡ 2 π k N ) D F T ( x 1 )

设： X 0 [ k ] = D F T ( x 0 ) ， X 1 [ k ] = D F T ( x 1 )

X [ k ] = X 0 [ k ] + ( cos ⁡ 2 π k N − i sin ⁡ 2 π k N ) X 1 [ k ]

因为这里 x 0 [ n ] 和 x 1 [ n ] 分别是 x[n]的偶数项和奇数项，因此它们的长度都是 N 2 ，所以 x 0 [ n ] 和 x 1 [ n ] 的长度也是 N 2 ，而 X [ k ] 是 x 0 [ n ] 和 x 1 [ n ] 的线性组合，因此其长度也是 N 2

前面已经提到，原式的变换中 X[k]的长度为 N ，经过一轮变换之后其长度变为了 N 2 ，那岂不是意味着丢失了一半的频率信息吗？下面就利用 DFT的对称性找回这一半丢失的信息。

在此之前先引入旋转因子的概念：

旋转因子表示为：
W N = e − i 2 π N = cos ⁡ 2 π N − i sin ⁡ 2 π N

旋转因子的对称性：

W N k + N 2 = e − i 2 π N ( k + N 2 ) = e − i 2 π N ∗ e − i π = − e − i 2 π N = − W N k

使用旋转因子表达将更简洁：

X [ k ] = X 0 [ k ] + W N k X 1 [ k ]

上面已经说到， X 0 [ k ] 和 X 1 [ k ] 的长度为 N 2 ，那么超出 N 2 的部分会出现什么情况呢？事实上超出的部分将会呈周期性变化，即： X 0 [ k + N 2 ] = X 0 [ k ] ，这一点将 k + N 2 代入 X 0 [ k ] 的表达式即可证明，这里就省去了。知道这一点后，我们用在上式中用 k + N 2 代替 k

X [ k + N 2 ] = X 0 [ k + N 2 ] + W N k + N 2 X 1 [ k + N 2 ] = X 0 [ k ] − W N k X 1 [ k ]

至此已经得出了整个频域上的表达式：

X [ k ] = X 0 [ k ] + W N k X 1 [ k ]

X [ k + N 2 ] = X 0 [ k ] − W N k X 1 [ k ]

这里写图片描述
下面举一个8点FFT的例子，假设使用x[0]~x[7]表示要进行FFT的8个点，用X[0]~X[7]表示对应频域上的点，那么根据以上思路，需先将0~7分为奇偶两组，即：{0,2,4,6}和{1,3,5,7}，然后继续分组：{0,4}和{2,6}、{1,5}和{3,7}，最后一次分组：{0}和{4}、{2}和{6}、{1}和{5}、{3}和{7}。经过三次分组后发现DFT的运算量已经是1了，这时再进行蝶形运算

这里写图片描述

FFT的基本步骤：
1. 根据序列长度 L ，计算出倒序的码位；
2. 算出加权序列 W u L ；
3. 计算 2 N − 1 个长度为2的离散傅里叶变换；
4. 基于前述的蝶形运算公式和分解流程，递归算出原序列共 2 N 个点的DFT。

FFT在数字信号处理作用

对于信号的频域分析，经常会用到FFT
当然，matlab这个东东都给我们封装好了库函数,不不，matlab叫做工具包~反正就是封装好了

Y = fft(y);
Y = fft(y,N);

参数
y是目标序列
N是所做FFT的点数
返回值
Y是序列的快速傅里叶变换
hint：
y可以是一向量或矩阵，若y为向量，则Y是y的FFT，并且与y具有相同的长度。若y为一矩阵，则Y是对矩阵的每一列向量进行FFT

随便举个栗子（水题）

一个由50Hz和120Hz正弦信号构成的信号，收到零均值的随机噪声干扰，数据采样频率为1000Hz，利用FFT分析其信号频率部分，画出时域波形以及信号的功率谱密度

直接贴代码了~太水了

t=0:0.001:0.8;
x=sin(2*pi*50*t)+cos(2*pi*120*t);%信号源
y = x + 1.5 * randn(1,length(t));%加入噪声
subplot(3,1,1);plot(t,x);%信号时域波形
subplot(3,1,2);plot(t,y);%掺杂后的时域波形

Y = fft(y,512);%进行FFT
P=Y.*conj(Y)/512;%利用共轭计算Y序列各点的模
f=1000*(0:255)/512;
subplot(3,1,3);
plot(f,P(1:256));%信号功率谱密度波形

这里写图片描述

FFT的其他地方的作用

理论

比如在信息竞赛中，FFT的应用也是很多的~比如FFT来求卷积，也就是多项式乘法
一般情况下，多项式乘法的时间复杂度为 o ( n 2 ) ，例如下面的方法

for(int i = 0;i < n;i++) {
   for(int j = 0;j < n;j++) {
       c[i+j] += a[i] * b[j];
   }
}

数据量变大之后，就会出现TLE

首先我们单独考虑一个 n 项（ n = 2 x )的多项式 A ( x ) ，其系数向量为 ( a 0 , a 1 , a 2 , . . . , a n − 1 ) 。我们将 n 次单位根的0 ~ n − 1 次幂分别带入 A ( x ) 得到其点值向量 ( A ( w n 0 ) , A ( w n 1 ) , A ( w n 2 ) , . . . , A ( w n n − 1 ) ) 。

所以我们必须要利用到单位根 ω 的特殊性质。

对于 A ( x ) = a 0 + a 1 ⋅ x 1 + a 2 ⋅ x 2 + a 3 ⋅ x 3 + . . . + a n − 1 ⋅ x n − 1

考虑将其按照奇偶分组

A ( x ) = ( a 0 + a 2 ⋅ x 2 + a 4 ⋅ x 4 . . . + a n − 2 ⋅ x n − 2 ) + ( a 1 ⋅ x 1 + a 3 ⋅ x 3 + a 5 ⋅ x 5 + . . . + a n − 1 ⋅ x n − 1 )

A ( x ) = ( a 0 + a 2 ⋅ x 2 + a 4 ⋅ x 4 . . . + a n − 2 ⋅ x n − 2 ) + x ⋅ ( a 1 + a 3 ⋅ x 2 + a 5 ⋅ x 4 + . . . + a n − 1 ⋅ x n − 2 )

令

A 1 ( x ) = ( a 0 + a 2 ⋅ x + a 4 ⋅ x 2 . . . + a n − 2 ⋅ x n − 2 2 )

A 2 ( x ) = ( a 1 + a 3 ⋅ x + a 5 ⋅ x 2 . . . + a n − 1 ⋅ x n − 2 2 )

则可得到

A ( x ) = A 1 ( x 2 ) + x ⋅ A 2 ( x 2 )

分类讨论

设 0 ≤ k ≤ n 2 − 1 ， k ∈ Z

A ( ω n k ) = A 1 ( ω n 2 k ) + ω n k ⋅ A 2 ( ω n 2 k )

折半引理

A ( ω n k ) = A 1 ( ω n 2 k ) + ω n k ⋅ A 2 ( ω n 2 k )

对于 n 2 ≤ k + n 2 ≤ n − 1

A ( ω n k + n 2 ) = A 1 ( ω n 2 k + n ) + ω n k + n 2 ⋅ A 2 ( ω n 2 k + n )

其中 ω n 2 k + n = ω n 2 k ⋅ ω n n = ω n 2 k = ω n 2 k

由消去引理 ω n k + n 2 = − ω n k

故 A ( ω n k + n 2 ) = A 1 ( ω n 2 k ) − ω n k ⋅ A 2 ( ω n 2 k )

注意， k 与 k + n 2 取遍了 [ 0 , n − 1 ] 中的 n 个整数，保证了可以由这 n 个点值反推解出系数.

如果已知了 A 1 ( x ) , A 2 ( x ) 分别在 ω n 2 0 , ω n 2 1 ,…, ω n 2 n 2 − 1 , 的取值，可以在 O ( n ) 的时间内求出 A ( x ) 的取值。

而 A 1 ( x ) , A 2 ( x ) 都是 A ( x ) 一半的规模，显然可以转化为子问题递归求解。

这里写图片描述

时间复杂度：

T ( n ) = 2 T ( n 2 ) + O ( n ) = O ( n l o g n )

实现

一个简单的模板 模板来源

struct Complex // 复数
{
    double r, i;
    Complex(double _r = 0, double _i = 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);
    }
};

递归实现

Complex* RecursiveFFT(Complex a[], int n)//n表示向量a的维数
{
    if(n == 1)
        return a;
    Complex wn = Complex(cos(2*PI/n), sin(2*PI/n));
    Complex w = Complex(1, 0);
    Complex* a0 = new Complex[n >> 1];
    Complex* a1 = new Complex[n >> 1];
    for(int i = 0; i < n; i++)
        if(i & 1) a1[(i - 1) >> 1] = a[i];
        else a0[i >> 1] = a[i];
    Complex *y0, *y1;
    y0 = RecursiveFFT(a0, n >> 1);
    y1 = RecursiveFFT(a1, n >> 1);
    Complex* y = new Complex[n];
    for(int k = 0; k < (n >> 1); k++)
    {
        y[k] = y0[k] + w*y1[k];
        y[k + (n >> 1)] = y0[k] - w*y1[k];
        w = w*wn;
    }
    delete a0;
    delete a1;
    delete y0;
    delete y1;
    return y;
}

非递归实现

void change(Complex y[], int len) // 二进制平摊反转置换 O(logn)  
{
    int i, j, k;
    for (i = 1, j = len / 2;i < len - 1;i++)
    {
        if (i < j)swap(y[i], y[j]);
        k = len / 2;
        while (j >= k)
        {
            j -= k;
            k /= 2;
        }
        if (j < k)j += k;
    }
}
void fft(Complex y[], int len, int on) //FFT:on=1; IFFT:on=-1
{
    change(y, len);
    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)
        for (int i = 0;i < len;i++)
            y[i].r /= len;
}

最后举个栗子

HDU4609 3-idiots
Time Limit: 10000/5000 MS (Java/Others) Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 7102 Accepted Submission(s): 2462

Problem Description
King OMeGa catched three men who had been streaking in the street. Looking as idiots though, the three men insisted that it was a kind of performance art, and begged the king to free them. Out of hatred to the real idiots, the king wanted to check if they were lying. The three men were sent to the king’s forest, and each of them was asked to pick a branch one after another. If the three branches they bring back can form a triangle, their math ability would save them. Otherwise, they would be sent into jail.
However, the three men were exactly idiots, and what they would do is only to pick the branches randomly. Certainly, they couldn’t pick the same branch - but the one with the same length as another is available. Given the lengths of all branches in the forest, determine the probability that they would be saved.

Input
An integer T(T≤100) will exist in the first line of input, indicating the number of test cases.
Each test case begins with the number of branches N(3≤N≤105).
The following line contains N integers a_i (1≤a_i≤105), which denotes the length of each branch, respectively.

Output
Output the probability that their branches can form a triangle, in accuracy of 7 decimal places.

Sample Input

2
4
1 3 3 4
4
2 3 3 4

Sample Output
0.5000000
1.0000000

Source
2013 Multi-University Training Contest 1

分析：我们考虑两边长度之和为n的方案数，设num[x]为长度为x的个数，那么 ∑ x = 1 n n u m [ n − x ] ∗ n u m [ x ] 即两边长度之和为n

的方案数。容易发现这这正是卷积！

然后我们就可以用FFT预处理出所有的两边长度之和为i的方案数。FFT求出来的第i项的系数就是方案数。由于FFT处理出来的是重复的，以及部分非法的（自己和自己构成两边之和），这些我们可以很容易的去除：设长度为x的方案数为ans，如果长度x为偶数，首先 a n s = a n s − n u m [ x 2 ] ，之后 a n s = a n s 2

接下来，我们枚举两边长度之和i，易知第三边长度 x ≥ i 的都是不合法的，这样我们减去就好了。这时候我很自然的想到了一种情况：如果i=5，其中包含情况x=1,y=4,x+y=i，此时如果第三条边长度为2，那不是不合法么？然而我们不需要担心。因为在x=1,y=2,x+y=3,z=4的时候我们就已经排除了，所以这样所有不合法的情况我们都恰好排除了。 

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <algorithm>
using namespace std ;

typedef long long LL ;

const int MAXN = 100005 ;
const double pi = acos ( -1.0 ) ;

struct Complex {
    double r , i ;
    Complex () {}
    Complex ( double r , double i ) : r ( r ) , i ( i ) {}
    Complex operator + ( const Complex& t ) const {
        return Complex ( r + t.r , i + t.i ) ;
    }
    Complex operator - ( const Complex& t ) const {
        return Complex ( r - t.r , i - t.i ) ;
    }
    Complex operator * ( const Complex& t ) const {
        return Complex ( r * t.r - i * t.i , r * t.i + i * t.r ) ;
    }
} ;

void FFT ( Complex y[] , int n , int rev ) {
    for ( int i = 1 , j , t , k ; i < n ; ++ i ) {
        for ( j = 0 , t = i , k = n >> 1 ; k ; k >>= 1 , t >>= 1 ) j = j << 1 | t & 1 ;
        if ( i < j ) swap ( y[i] , y[j] ) ;
    }
    for ( int s = 2 , ds = 1 ; s <= n ; ds = s , s <<= 1 ) {
        Complex wn ( cos ( rev * 2 * pi / s ) , sin ( rev * 2 * pi / s ) ) , w ( 1 , 0 ) , t ;
        for ( int k = 0 ; k < ds ; ++ k , w = w * wn ) {
            for ( int i = k ; i < n ; i += s ) {
                y[i + ds] = y[i] - ( t = w * y[i + ds] ) ;
                y[i] = y[i] + t ;
            }
        }
    }
    if ( rev == -1 ) for ( int i = 0 ; i < n ; ++ i ) y[i].r /= n ;
}

int num[MAXN] , sum[MAXN] ;
int n ;
Complex y[MAXN << 2] ;

void solve () {
    int x , n1 = 1 , maxv = 0 ;
    clr ( num , 0 ) ;
    scanf ( "%d" , &n ) ;
    for ( int i = 0 ; i < n ; ++ i ) {
        scanf ( "%d" , &x ) ;
        num[x] ++ ;
        maxv = max ( maxv , x ) ;
    }
    while ( n1 <= 2 * maxv ) n1 <<= 1 ;
    for ( int i = 0 ; i <= maxv ; ++ i ) y[i] = Complex ( num[i] , 0 ) ;
    for ( int i = maxv + 1 ; i < n1 ; ++ i ) y[i] = Complex ( 0 , 0 ) ;
    FFT ( y , n1 , 1 ) ;
    for ( int i = 0 ; i < n1 ; ++ i ) y[i] = y[i] * y[i] ;
    FFT ( y , n1 , -1 ) ;
    for ( int i = 1 ; i <= maxv ; ++ i ) sum[i] = sum[i - 1] + num[i] ;
    LL tot = ( LL ) n * ( n - 1 ) * ( n - 2 ) / 6 , ans = tot ;
    for ( int i = 2 ; i <= maxv ; ++ i ) {
        LL x = ( LL ) ( y[i].r + 0.5 ) ;
        if ( i % 2 == 0 ) x -= num[i / 2] ;
        x /= 2 ;
        ans -= x * ( n - sum[i - 1] ) ;
    }
    printf ( "%.7f\n" , ( double ) ans / tot ) ;
}

int main () {
    int T ;
    scanf ( "%d" , &T ) ;
    for ( int i = 1 ; i <= T ; ++ i ) solve () ;
    return 0 ;
}