快速傅里叶变换

2014.12.06 15:26 Sat| 15 visits oi_2015| 2015_算法笔记| Text

知识铺垫

时域和频域

时域

时域(Time domain)描述数学函数或物理信号对时间的关系。例如一个信号的时域波形可以表达信号随着时间的变化。

若考虑离散时间,时域中的函数或信号,在各个离散时间点的数值均为已知。若考虑连续时间,则函数或信号在任意时间的数值均为已知。

频域

频域(Frequency domain)是指在对函数或信号进行分析时,分析其和频率有关部份,而不是和时间有关的部份,和时域一词相对。

以信号为例,信号在时域下的图形可以显示信号如何随着时间变化,而信号在频域下的图形(一般称为频谱)可以显示信号分布在哪些频率及其比例。频域的表示法除了有各个频率下的大小外,也会有各个频率的相位,利用大小及相位的资讯可以将各频率的弦波给予不同的大小及相位,相加以后可以还原成原始的信号。

函数或信号可以透过一对数学的运算子在时域及频域之间转换。例如傅里叶变换可以将一个时域信号转换成在不同频率下对应的振幅及相位,其频谱就是时域信号在频域下的表现,而反傅里叶变换可以将频谱再转换回时域的信号。

复数

复数的基本概念

复数为实数的延伸,它使任一多项式方程式都有根。复数当中有个“虚数单位” $i$ ,它是-1的一个平方根,即 $i ^2 = -1$ 。任一复数都可表达为 $x + yi$ ,其中 $x$ 及 $y$ 皆为实数,分别称为复数之“实部”和“虚部”。

如果复数 $z = a + ib$ ,则实部 $a$ 被指示为$\Re(z)$ ,而虚部 $b$ 被指示为 $\Im(z)$ 。

所有复数的集合记为 $\Bbb C$ 。通过把实数的所有成员看作复数 $a = a + 0i$ ,实数 $\Bbb R$ 可以被当作 $\Bbb C$ 的子集。

复数中的虚数是无法比较大小的,即两个虚数只有相等和不等两种等量关系。

复数的基本运算

通过形式上应用代数的结合律、交换律和分配律,再加上等式 $i^2 = -1$ ,定义复数的加法、减法、乘法和除法:

  • 加法: $(a + bi) + (c + di) = (a + c) + (b + d)i$
  • 减法: $(a + bi) - (c + di) = (a - c) + (b - d)i$
  • 乘法: $(a + bi) (c + di) = ac + bci + adi + bd i^2 = (ac - bd) + (bc + ad)i$
  • 除法: $\frac{(a + bi)}{(c + di)} = \frac{(a+bi)(c-di)}{(c+di)(c-di)} =\frac{ac+bci-adi-bd i^2}{c^2 -(di)^2}=\frac{(ac+bd)+(bc-ad)i}{c^2+d^2}=\left({ac + bd \over c^2 + d^2}\right) + \left( {bc - ad \over c^2 + d^2} \right)i$

欧拉公式、复数的模和幅角

  • 欧拉公式: $\cos \theta + i \sin \theta = e ^{i \theta}$
  • 复数 $z = re^{i\phi}$ 的模值是 $z$ 。如果 $z=a+bi$ ,则 $|z| = \sqrt{a^2+b^2}$ 。
  • 复数 $z = re^{i\phi}$ 的幅角为 $\phi$ 。此值对模 $2\pi$而言是唯一的。

复平面、复数的两种表示形式

复平面(complex plane)是用水平的实轴与垂直的虚轴建立起来的复数的几何表示。

  • 笛卡尔形式
    建立二维坐标系。将复数 $z=a十bi$ 写做 $(a,b)$ ,则可以用笛卡尔坐标指定。复数的笛卡尔坐标表示叫做复数的“笛卡尔形式”、“直角形式”或“代数形式”。

  • 极坐标形式
    建立极坐标系。将复数 $re ^{i \theta}$ 写做 $(r,\theta)$ ,则可以用极坐标指定。复数的极坐标表示叫做复数的“极坐标形式”。

  • 从极坐标形式到笛卡儿坐标形式的转换

周期信号的表示

信号分析和其他领域使用复数可以方便的表示周期信号。模值 $|z|$ 表示信号的幅度,辐角 $arg(z)$ 表示给定频率的正弦波的相位。

利用傅里叶变换可将实信号表示成一系列周期函数的和。这些周期函数通常用形式如下的复函数的实部表示:

其中 $ω$ 对应角频率,复数 $z$ 包含了幅度和相位的信息。

单位根

方程 $z^n = 1 \qquad (n = 1, 2, 3, \cdots )$ 的复数根 $z$ 为 $n$ 次单位根。

$n$ 次单位根有 $n$ 个: $e^{\frac{2 \pi k {\mathrm{i}} }{n} } \qquad (k = 0, 1, 2, \cdots, n - 1)$ 。它们位于复平面的单位圆上,构成正 $n$ 边形的顶点,其中一个顶点是1。

单位根的性质

下面,我们用 $\omega_n$ 来表示 n 次主单位根 $e^{-i\frac{2\pi}{n}}$。

$\omega_n$ 的性质:

  • 周期性: $\omega_n$ 具有周期 $n$ ,即 $\omega_n^{k+n}=\omega_n^{k}$ 。
  • 对称性: $\omega_n^{k+\frac{n}{2}}=-\omega_n^{k}$ 。
  • 若 m 是 n 的约数, $\omega_n^{mkn}=W_{\frac{n}{m}}^{kn}$ 。

卷积

在泛函分析中,卷积是通过两个函数 $f$ 和 $g$ 生成第三个函数的一种数学算子,表征函数 $f$ 与经过翻转和平移的 $g$ 的重叠部分的面积。

定义

函数 $f$ 与 $g$ 的卷积记作 $f * g$ ,它是其中一个函数翻转并平移后与另一个函数的乘积的积分,是一个对平移量的函数。

积分区间取决于 $f$ 与 $g$ 的定义域。

对于定义在离散域的函数,卷积定义为

利用定义计算卷积

以上过程总时间复杂度为 $\mathcal O(MN)$ 。

傅里叶变换

傅里叶变换(Fourier transform)是一种线性的积分变换,常在将信号在时域(或空域)和频域之间变换时使用,在物理学和工程学中有许多应用。因其基本思想首先由法国学者约瑟夫·傅里叶系统地提出,所以以其名字来命名以示纪念。

约瑟夫·傅里叶男爵(法语:Joseph Fourier,1768年3月21日-1830年5月16日),法国数学家、物理学家,提出傅里叶级数,并将其应用于热传导理论与振动理论,傅里叶变换也以他命名。他被归功为温室效应的发现者。

经过傅里叶变换而生成的函数 $\hat f$ 称作原函数 $f$ 的傅里叶变换。傅里叶变换是可逆的,即可通过$\hat f$确定其原函数 $f$ 。通常情况下,$f$是实数函数,而 $\hat f$ 则是复数函数,用一个复数来表示振幅和相位。

傅里叶变换家族

傅里叶变换家族较为庞大。在数学领域的谐波分析中,连续傅里叶变换(continuous Fourier transform, CFT)与傅里叶级数 (Fourier series, FS)有非常微妙的关系。而且连续傅里叶变换也与离散时间傅里叶变换(discrete time Fourier transform, DTFT)和离散傅里叶变换(discrete Fourier transform, DFT)有很近的关系。傅里叶变换家族通常就是指这四种变换。观察这四种变换的特性可以发现:函数在时(频)域的离散对应于其像函数在频(时)域的周期性.反之连续则意味着在对应域的信号的非周期性。

从DFT到FFT

离散傅里叶变换

离散傅里叶变换(Discrete Fourier Transform),是傅里叶变换在时域和频域上都呈离散的形式,将信号的时域采样变换为其离散时间傅里叶变换的频域采样。在形式上,变换两端(时域和频域上)的序列是有限长的,而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号作DFT,也应当将其看作其周期延拓的变换。在实际应用中通常采用快速傅里叶变换计算DFT。

对于 $n$ 点序列 $a$ ,它的离散傅里叶变换(DFT)为

其中 $e$ 是自然对数的底数, $i$ 是虚数单位。通常以符号 $\mathcal{F}$ 表示这一变换,即

离散傅里叶变换的逆变换(IDFT)为:

可以记为:

快速傅里叶变换

库利-图基算法是最常见的FFT算法。这一方法以分治法为策略递归地将长度为 $N=N_1 N_2$ 的离散傅里叶变换分解为长度为 $N_1$ 的 $N_2$ 个较短串行的离散傅里叶变换,以及与 $\mathcal O(N)$ 个旋转因子的复数乘法。

这种方法以及FFT的基本思路在1965年J. W. Cooley和J. W. Tukey合作发表An algorithm for the machine calculation of complex Fourier series之后开始为人所知。但后来发现,实际上这两位作者只是重新发明了高斯在1805年就已经提出的算法(此算法在历史上数次以各种形式被再次提出)。

设计思想

为了简单起见,我们下面设待变换串行长度 $n=2^r$ 。 根据上面单位根的对称性,求级数 $y_k=\sum_{n=0}^{n-1} \omega_n^{kn}x_n$ 时,可以将求和区间分为两部分:

由此式可以计算出 $y_k$ 的前 $n/2$ 个点。对于后 $n/2$ 个点,注意 $\mathcal F_{odd}(k)$ 和 $\mathcal F_{even}(k)$ 都是周期为 $n/2$ 的函数,由单位根的对称性,于是有以下变换公式:

这样,一个 $n$ 点变换就分解成了两个 $n/2$ 点变换,可递归求解,此时算法的时间复杂度为 $\mathcal O(n\log n)$ 。

实际应用

频谱分析

DFT是连续傅里叶变换的近似,因此可以对连续信号 x(t) 均匀采样并截断以得到有限长的离散序列 ${x_n}$ ,对这一序列作离散傅里叶变换,可以分析连续信号 x(t) 频谱的性质。

数据压缩

由于人类感官的分辨能力存在极限,因此很多有损压缩算法利用这一点将语音、音频、图像、视频等信号的高频部分除去。高频信号对应于信号的细节,滤除高频信号可以在人类感官可以接受的范围内获得很高的压缩比。这一去除高频分量的处理就是通过离散傅里叶变换完成的。将时域或空域的信号转换到频域,仅储存或传输较低频率上的系数,在解压缩端采用逆变换即可重建信号。

解偏微分方程

离散傅里叶变换及其多维形式在偏微分方程的求解中也有应用。此时DFT被看作傅里叶级数的近似。傅里叶级数将函数在复指数 einx 上展开,这正是微分算子的特征方程: d/dx einx = in einx 。因此,通过傅里叶级数的形式,线性常微分方程被转换为代数方程,而后者是很容易求解的。此时得到的结果是偏微分方程解的级数表示,只要通过DFT逆变换即可得到其一般表示。这种方法被称作谱方法或级数解法。

正交频分复用

正交频分复用(Orthogonal frequency-division multiplexing)在宽带无线通信中有重要的应用。这种技术将带宽分为N个等间隔的子载波,可以证明这些子载波相互正交。尤其重要的是,OFDM调制可以由IDFT实现,而解调可以由DFT实现。OFDM还利用DFT的移位性质,在每个帧头部加上循环前缀(Cyclic Prefix),使得只要信道延时小于循环前缀的长度,就能消除信道延时对传输的影响。

FFT在信息学竞赛中的应用

卷积定理

卷积定理指出,函数卷积的傅里叶变换是函数傅里叶变换的乘积。即一个域中的卷积对应于另一个域中的乘积,例如时域中的卷积对应于频域中的乘积。

其中 $\mathcal{F}(f)$ 表示 $f$ 的傅里叶变换。下面这种形式也成立:

借由傅里叶逆变换 $\mathcal{F}^{-1}$ ,也可以写成

长整数与多项式乘法

目前长整数或多项式乘法最快速的算法是基于离散傅里叶变换的。由于整数(或多项式)乘法是逐位(或逐项)乘累加的形式,因此整数(或多项式)乘积的数字(或系数)可以用乘数数字(或乘式系数)的卷积表示。利用卷积定理,只要将数字(或系数)序列通过离散傅里叶变换变到频域,就可以将逐个乘累加的卷积变为对位的乘法,从而减少计算量,再以一次逆变换便可以得到乘法结果。需要注意整数乘法还有进位的问题。下面以多项式乘法为例说明这一应用:

假设需要计算多项式乘法 $c(x) = a(x)·b(x)$ ,这一乘积结果的各项系数的表达式实际上是线性卷积的形式。由于离散傅里叶变换对应于圆周卷积,因此需要将这两个乘式的系数按照次数升序排列,序列后“补零”,使它们的长度d大于两式最高项次数的和,然后作圆周卷积:

其中 $c$ 就是 $c(x)$ 系数的向量。由于圆周卷积的DFT是乘积,有:

利用快速傅里叶变换,这一算法的运算复杂度为$\mathcal O(n \log n)$。

对于以上过程的理解

每一个多项式都可以表示成点-值形式,即其在点集 ${x}$ 中每个点处的取值。可以证明,有且仅有一个次数小于 $n$ 的多项式满足“点-值”序列 ${(x_0,y_0),(x_1,y_1),\ldots,(x_{n-1},y_{n-1})}$ 。因此若两个多项式的点-值表示中点集相同,则将它们对应的值相乘,所得的结果即为 $A(x)$ 与 $B(x)$ 的乘积的点-值表达式。此时我们只要快速计算出两个次数均不超过 $n$ 的单变量多项式 $A(x)$ 和 $B(x)$ 在 $2n$ 个点处的取值,便可在$\mathcal O(n)$的时间内求解它们之间的乘积在这些点处的取值。将这一取值转化为系数表示,便是两多项式间的乘积。这里FFT可以帮助我们快速进行一个多项式的系数表示与它在单位根处的点-值表示之间的转换。

具体实现步骤

  1. 补零
    在两个多项式的最前面补零,得到两个 $2n$ 次多项式,设系数向量分别为 $v_1$ 和 $v_2$ 。

  2. 求值
    用FFT计算 $\hat{f_1}=\mathcal{F}v_1$ 和 $\hat{f_2}=\mathcal{F}v_2$ 。这里得到的 $\hat{f_1}$ 和 $\hat{f_2}$ 分别是两个多项式在 $2n$ 次单位根处的各个取值(即点-值表示)。

  3. 乘法
    把两个向量 $\hat{f_1}$ 和 $\hat{f_2}$ 的每一位对应相乘,得到向量 $\hat{f}$ 。它对应输入多项式乘积的点-值表示。

  4. 插值
    用FFT计算 $v=\mathcal{F}^{-1}\hat{f}$ ,其中 $v$ 就是乘积的系数向量。

代码实现

递归实现

void FFT(complex<double> x[],int n,int type)
{
    if(n==1) return;
    complex<double> l[n>>1],r[n>>1];
    for(int i=0;i<n;i+=2)
        l[i>>1]=x[i], r[i>>1]=x[i+1];
    FFT(l,n>>1,type), FFT(r,n>>1,type);
    complex<double> wn(cos(type*2*pi/n),sin(type*2*pi/n)),w(1,0);
    for(int i=0;i<n>>1;w*=wn,i++)
        x[i]=l[i]+w*r[i], x[i+(n>>1)]=l[i]-w*r[i];
}

非递归实现

利用FFT非递归实现长整数乘法完整代码如下:

#include <bits/stdc++.h>
using namespace std;
#define N 131075

int n,c[N];
complex<double> a[N],b[N],p[N];
const double pi=3.14159265358979323;

void FFT(complex<double> x[],int n,int type)
{
    for(int i=0,t=0;i<n;i++)
    {
        if(i>t) swap(x[i],x[t]);
        for(int j=n>>1;(t^=j)<j;j>>=1);
    }
    for(int s=2;s<=n;s<<=1)
    {
        complex<double> wn(cos(type*2*pi/s),sin(type*2*pi/s));
        for(int i=0;i<n;i+=s)
        {
            complex<double> w(1,0),t;
            for(int j=0;j<s>>1;w*=wn,j++)
                t=w*x[i+j+(s>>1)],
                x[i+j+(s>>1)]=x[i+j]-t, x[i+j]+=t;
        }
    }
}

int main()
{
    cin>>n;
    getchar(); for(int i=n-1;i>=0;i--)
        a[i]=getchar()-'0';
    getchar(); for(int i=n-1;i>=0;i--)
        b[i]=getchar()-'0';
    for(int t=n,i=1;i>>2<t;i<<=1) n=i;
    FFT(a,n,1), FFT(b,n,1);
    for(int i=0;i<n;i++) p[i]=a[i]*b[i];
    FFT(p,n,-1);
    int len=0;
    for(int i=0;i<n;i++)
        c[i]=p[i].real()/n+0.1;
    for(int i=0;i<n;i++)
        if(c[i]) len=i, c[i+1]+=c[i]/10, c[i]%=10;
    for(int i=len;i>=0;i--)
        printf("%d",c[i]);
    return 0;
}

注:以上代码中 n 为两因数的位数。