数字信号处理入门:离散傅里叶变换 DFT 与快速傅里叶变换 FFT

从 DTFT 的连续频谱出发,系统整理 DFT 的定义、周期延拓、圆周移位、循环卷积与 FFT 的蝶形分解,帮助建立离散频域分析的整体框架。

数字信号处理入门:离散傅里叶变换 DFT 与快速傅里叶变换 FFT

这篇文章继续沿着数字信号处理的主线,讲清楚离散傅里叶变换 DFT 与快速傅里叶变换 FFT。前面学习过 DTFT、Z 变换和系统频响,但真正到了计算机实现时,会遇到一个现实问题:计算机不能直接处理连续频率曲线,也不能承受过大的计算量。因此,DFT 要解决两个核心问题:频谱如何离散化,以及离散化以后怎样快速计算。

本文按讲义顺序展开,但去掉了原 PPT 页码和模板式标题,改写为适合博客发布的章节。读者只需要知道序列、单位冲激、卷积、DTFT 的基本定义,就可以从 0 开始理解 DFT 和 FFT。

傅里叶工具关系图

为什么需要 DFT

DTFT 的形式是

X(ejω)=n=x(n)ejωn,X(e^{j\omega})=\sum_{n=-\infty}^{\infty}x(n)e^{-j\omega n},

反变换是

x(n)=12πππX(ejω)ejωndω.x(n)=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(e^{j\omega})e^{j\omega n}d\omega.

它的特点是:时域是离散的,但频域仍然是连续变量 ω\omega。连续频率曲线不适合直接交给计算机逐点处理。计算机更喜欢有限个数值,例如 X(0),X(1),,X(N1)X(0),X(1),\cdots,X(N-1)。所以我们希望把频域也离散化。

这里有一个重要规律:时域周期化会导致频域离散化;时域离散化会导致频域周期化。DTFT 已经做到“时域离散”,所以频域是周期的;若进一步想让频域也离散,就要在时域引入周期延拓。DFT 正是建立在这个思想上。

DTFT到DFT的频域采样

从 DTFT 到 DFT:有限长序列的周期延拓

x(n)x(n) 是长度为 MM 的有限长序列。为了得到 NN 点 DFT,通常要求 NMN\ge M。若 N>MN>M,就要在序列末尾补零,使它成为 NN 点序列。DFT 的思想不是孤立地看这 NN 个点,而是把它看成一个周期序列的主值区间。

这句话很关键:DFT 只对有限长序列作周期延拓或周期序列成立。在 DFT 中,有限长序列总是隐含周期性。

周期延拓与主值区间

所谓主值区间,就是从 n=0n=0n=N1n=N-1 的第一个周期。周期延拓后的序列可以写作

x~(n)=m=x(n+mN)=x((n))N,\tilde{x}(n)=\sum_{m=-\infty}^{\infty}x(n+mN)=x((n))_N,

其中 ((n))N((n))_N 表示对 NN 取模后的余数。比如

((25))9=7,((4))9=5.((25))_9=7,\qquad ((-4))_9=5.

余数运算的目的,是把任意下标重新折回到 0N10\sim N-1 的主值区间内。

DFT 的定义

x(n)x(n) 是长度为 MM 的有限长序列,取 NMN\ge M,则 NN 点 DFT 定义为

X(k)=DFT[x(n)]N=n=0N1x(n)ej2πNkn=n=0N1x(n)WNkn,X(k)=DFT[x(n)]_N=\sum_{n=0}^{N-1}x(n)e^{-j\frac{2\pi}{N}kn} =\sum_{n=0}^{N-1}x(n)W_N^{kn},

其中

WN=ej2πNW_N=e^{-j\frac{2\pi}{N}}

称为旋转因子。

对应的反变换为

x(n)=IDFT[X(k)]N=1Nk=0N1X(k)ej2πNkn=1Nk=0N1X(k)WNkn.x(n)=IDFT[X(k)]_N=\frac{1}{N}\sum_{k=0}^{N-1}X(k)e^{j\frac{2\pi}{N}kn} =\frac{1}{N}\sum_{k=0}^{N-1}X(k)W_N^{-kn}.

从定义可以看出,DFT 把 NN 个时域样本 x(0),x(1),,x(N1)x(0),x(1),\cdots,x(N-1) 变成 NN 个频域样本 X(0),X(1),,X(N1)X(0),X(1),\cdots,X(N-1)。IDFT 则从这 NN 个频域样本恢复主值区间内的时域序列。

旋转因子为什么重要

旋转因子 WN=ej2π/NW_N=e^{-j2\pi/N} 是 DFT 与 FFT 的核心。它本质上是复平面单位圆上的一个等角度旋转。

它有四个常用性质。

周期性:

WNn=WNn+iN,iZ.W_N^n=W_N^{n+iN},\qquad i\in\mathbb{Z}.

共轭对称性:

WNn=(WNn).W_N^n=(W_N^{-n})^*.

正交性:

1Nk=0N1WN(nm)k={1,nm=iN0,nmiN.\frac{1}{N}\sum_{k=0}^{N-1}W_N^{(n-m)k}=\begin{cases}1,& n-m=iN\\0,& n-m\ne iN\end{cases}.

可约性:

WiNin=WNn.W_{iN}^{in}=W_N^n.

正交性保证了 DFT 和 IDFT 能互相恢复;可约性和周期性则是 FFT 能减少计算量的基础。

矩形序列的 8 点与 16 点 DFT

x(n)=R4(n),x(n)=R_4(n),

也就是 x(0)=x(1)=x(2)=x(3)=1x(0)=x(1)=x(2)=x(3)=1,其余点为 0。

当取 8 点 DFT 时:

X(k)=n=03ej2π8kn=ej3π8ksin(πk/2)sin(πk/8),0k7.X(k)=\sum_{n=0}^{3}e^{-j\frac{2\pi}{8}kn} =e^{-j\frac{3\pi}{8}k}\frac{\sin(\pi k/2)}{\sin(\pi k/8)},\qquad 0\le k\le 7.

当取 16 点 DFT 时:

X(k)=n=03ej2π16kn=ej3π16ksin(πk/4)sin(πk/16),0k15.X(k)=\sum_{n=0}^{3}e^{-j\frac{2\pi}{16}kn} =e^{-j\frac{3\pi}{16}k}\frac{\sin(\pi k/4)}{\sin(\pi k/16)},\qquad 0\le k\le 15.

同一个时域序列,取 8 点 DFT 和 16 点 DFT 的结果不同,这是因为频率采样点数不同。16 点 DFT 不是“换了信号”,而是在频域上采得更密。

DFT 与 Z 变换、DTFT 的关系

有限长序列的 Z 变换为

X(z)=n=0N1x(n)zn.X(z)=\sum_{n=0}^{N-1}x(n)z^{-n}.

DFT 为

X(k)=n=0N1x(n)WNkn.X(k)=\sum_{n=0}^{N-1}x(n)W_N^{kn}.

比较可知,当

z=WNk=ej2πNkz=W_N^{-k}=e^{j\frac{2\pi}{N}k}

时,Z 变换在单位圆上的取值就是 DFT:

X(k)=X(z)z=ej2πNk.X(k)=X(z)\bigg|_{z=e^{j\frac{2\pi}{N}k}}.

所以,DFT 可以理解为对 Z 变换在单位圆上的 NN 点等间隔采样。

对 DTFT 而言,如果 x(n)x(n) 是长度为 NN 的有限长序列,则

X(ejω)=n=0N1x(n)ejωn.X(e^{j\omega})=\sum_{n=0}^{N-1}x(n)e^{-j\omega n}.

ω=2πNk,\omega=\frac{2\pi}{N}k,

就得到

X(k)=X(ejω)ω=2πNk.X(k)=X(e^{j\omega})\bigg|_{\omega=\frac{2\pi}{N}k}.

因此,DFT 也是对 DTFT 连续频谱在 [0,2π)[0,2\pi) 上的 NN 点等间隔采样。

DFT是DTFT采样

解答题例 1:由 DTFT 采样得到 4 点 DFT

例如,有限长序列

x(n)=[1,0,3,1]x(n)=[1,0,3,-1]

的 DTFT 为

X(ejω)=1+3e2jωe3jω.X(e^{j\omega})=1+3e^{-2j\omega}-e^{-3j\omega}.

它的 4 点 DFT 是

X(k)={3,2j,5,2+j}.X(k)=\{3,-2-j,5,2+j\}.

这 4 个值就是连续频谱在 4 个等间隔频率点上的样本。

解答

这里的关键不是重新推一遍 DTFT,而是直接使用

X(k)=X(ejω)ω=2πNkX(k)=X(e^{j\omega})\bigg|_{\omega=\frac{2\pi}{N}k}

这条采样关系。对长度为 4 的序列做 4 点 DFT,本质上就是把连续频谱在

ω=0,π2,π,3π2\omega=0,\frac{\pi}{2},\pi,\frac{3\pi}{2}

四个频率点代入求值,因此得到

X(k)={3,2j,5,2+j}.X(k)=\{3,-2-j,5,2+j\}.

这个例题说明:DFT 不是另一套脱离 DTFT 的定义,而是 DTFT 在有限频率点上的离散采样。

DFT 的线性性质

DFT[x1(n)]=X1(k),DFT[x2(n)]=X2(k),DFT[x_1(n)]=X_1(k),\qquad DFT[x_2(n)]=X_2(k),

DFT[ax1(n)+bx2(n)]=aX1(k)+bX2(k).DFT[ax_1(n)+bx_2(n)]=aX_1(k)+bX_2(k).

这就是 DFT 的线性性质。若两个序列长度不同,应取

N=max(N1,N2),N=\max(N_1,N_2),

再把短序列补零到 NN 点。

圆周移位:DFT 中的“绕圈移动”

由于 DFT 隐含周期性,所以普通移位要变成圆周移位。对有限长序列 x(n)x(n),其圆周移位定义为

y(n)=x((n+m))NRN(n).y(n)=x((n+m))_N R_N(n).

操作过程可以分成三步:先把 x(n)x(n) 作周期延拓;再对周期延拓序列移位;最后截取主值区间。

圆周移位示意

圆周移位的直观含义是:如果某个样本从主值区间一端移出去,它会从另一端重新移进来。把 x(0),x(1),,x(N1)x(0),x(1),\cdots,x(N-1) 排在一个 NN 等分圆周上,圆周移位就相当于沿圆周旋转。

时域圆周移位的 DFT 性质为

xm(n)=x((n+m))NRN(n),x_m(n)=x((n+m))_N R_N(n),

Xm(k)=WNkmX(k).X_m(k)=W_N^{-km}X(k).

这说明:有限长序列的圆周移位只引入相位因子,不改变频谱幅度。

频域圆周移位也有对应的时域意义。若频域发生圆周移位,则时域相当于乘上复指数,也就是调制。常见形式为

IDFT[X((k+l))NRN(k)]=WNnlx(n).IDFT[X((k+l))_N R_N(k)]=W_N^{nl}x(n).

如果时域乘余弦或正弦,频域会出现左右两个圆周移位分量。这就是调制会搬移频谱的数学原因。

循环卷积:DFT 世界里的卷积

h(n)h(n)x(n)x(n) 的长度分别为 NNMM。它们的 LL 点循环卷积定义为

yc(n)=h(n)x(n)=[m=0L1h(m)x((nm))L]RL(n),y_c(n)=h(n)\otimes x(n)=\left[\sum_{m=0}^{L-1}h(m)x((n-m))_L\right]R_L(n),

其中 Lmax(N,M)L\ge \max(N,M)

循环卷积和普通线性卷积的区别在于:循环卷积中的下标按 LL 取模,超出范围的部分会绕回主值区间。

循环卷积可以用矩阵法计算。把一个序列写成循环倒相序列,再逐行向右循环移位,就可以形成循环卷积矩阵。若 x(n)x(n)h(n)h(n) 的长度小于 LL,要先在序列末尾补零到 LL 点。循环卷积满足交换律:

x(n)h(n)=h(n)x(n).x(n)\otimes h(n)=h(n)\otimes x(n).

线性卷积与循环卷积

解答题例 2:4 点与 8 点循环卷积的差别

例如

h(n)={1,2,3,4},x(n)={1,1,1,1}.h(n)=\{1,2,3,4\},\qquad x(n)=\{1,1,1,1\}.

4 点循环卷积为

yc(n)={10,10,10,10}.y_c(n)=\{10,10,10,10\}.

若补零后作 8 点循环卷积,则结果为

yc(n)={1,3,6,10,9,7,4,0}.y_c(n)=\{1,3,6,10,9,7,4,0\}.

这里 8 点循环卷积刚好等于线性卷积结果,因为线性卷积长度为

4+41=7,4+4-1=7,

L=8L=8 已经足够避免混叠。

解答

4 点循环卷积把结果长度固定在 4 个样本内,线性卷积尾部会按模 4 折回前面,所以得到

yc(n)={10,10,10,10}.y_c(n)=\{10,10,10,10\}.

当两个序列都补零到 8 点以后,再做 8 点循环卷积,此时

L=84+41=7,L=8 \ge 4+4-1=7,

已经满足“循环卷积等于线性卷积”的条件,因此结果恢复为

yc(n)={1,3,6,10,9,7,4,0}.y_c(n)=\{1,3,6,10,9,7,4,0\}.

这个例题真正要记住的是判据:若想用 DFT 计算线性卷积,循环长度至少要满足

LN+M1.L \ge N+M-1.

循环卷积定理

DFT[x1(n)]N=X1(k),DFT[x2(n)]N=X2(k),DFT[x_1(n)]_N=X_1(k),\qquad DFT[x_2(n)]_N=X_2(k),

DFT[x1(n)x2(n)]N=X1(k)X2(k).DFT[x_1(n)\otimes x_2(n)]_N=X_1(k)X_2(k).

这就是时域循环卷积定理。它告诉我们:时域循环卷积可以转化为频域相乘。

反过来,若时域逐点相乘,则频域是循环卷积:

DFT[x1(n)x2(n)]N=1NX1(k)X2(k).DFT[x_1(n)x_2(n)]_N=\frac{1}{N}X_1(k)\otimes X_2(k).

用DFT计算循环卷积

解答题综合例:频域相乘对应时域循环卷积

x(n)=δ(n)+2δ(n5),x(n)=\delta(n)+2\delta(n-5),

求它的 10 点 DFT。

由定义得

X(k)=1+2W105k=1+2ejπk=1+2(1)k.X(k)=1+2W_{10}^{5k}=1+2e^{-j\pi k}=1+2(-1)^k.

Y(k)=ej2k2π10X(k)=W102kX(k),Y(k)=e^{j2k\frac{2\pi}{10}}X(k)=W_{10}^{-2k}X(k),

则这相当于 x(n)x(n) 向左圆周移位 2 点,因此

y(n)=x((n+2))10R10(n)=2δ(n3)+δ(n8).y(n)=x((n+2))_{10}R_{10}(n)=2\delta(n-3)+\delta(n-8).

Y(k)=X(k)W(k),Y(k)=X(k)W(k),

其中 W(k)W(k) 是序列

w(n)={1,0n60,其他w(n)=\begin{cases}1,&0\le n\le 6\\0,&\text{其他}\end{cases}

的 10 点 DFT,则频域相乘对应时域循环卷积。讲义给出的结果为

y(n)={3,3,1,1,1,3,3,2,2,2}.y(n)=\{3,3,1,1,1,3,3,2,2,2\}.
解答

先求原序列的 10 点 DFT:

X(k)=1+2W105k=1+2(1)k.X(k)=1+2W_{10}^{5k}=1+2(-1)^k.

若频域乘上

W102k,W_{10}^{-2k},

就对应时域左移 2 点的圆周移位,所以立即得到

y(n)=x((n+2))10R10(n)=2δ(n3)+δ(n8).y(n)=x((n+2))_{10}R_{10}(n)=2\delta(n-3)+\delta(n-8).

若频域改成

Y(k)=X(k)W(k),Y(k)=X(k)W(k),

则要改用循环卷积定理理解,因为频域逐点相乘对应时域循环卷积。于是输出由 x(n)x(n)w(n)w(n) 的 10 点循环卷积给出,结果正是

y(n)={3,3,1,1,1,3,3,2,2,2}.y(n)=\{3,3,1,1,1,3,3,2,2,2\}.

这个综合例把两条性质串起来了:频域乘相位因子对应时域圆周移位,频域逐点相乘对应时域循环卷积。

复共轭序列的 DFT

DFT[x(n)]N=X(k),DFT[x(n)]_N=X(k),

则有

DFT[x(n)]N=X(Nk),DFT[x^*(n)]_N=X^*(N-k),

以及

DFT[x(Nn)]N=X(k).DFT[x^*(N-n)]_N=X^*(k).

这些公式说明,取共轭、反转与 DFT 结果之间存在对应关系。它们后面会用于推导实序列 DFT 的对称性。

DFT 的共轭对称性

DTFT 中的对称性通常是关于原点讨论的,而 DFT 的序列定义在 0N10\sim N-1,所以有限长序列的对称性通常是关于 N/2N/2 点讨论的。

有限长共轭对称序列满足

xep(n)=xep(Nn),0nN1.x_{ep}(n)=x^*_{ep}(N-n),\qquad 0\le n\le N-1.

有限长共轭反对称序列满足

xop(n)=xop(Nn),0nN1.x_{op}(n)=-x^*_{op}(N-n),\qquad 0\le n\le N-1.

任意有限长序列都可以拆成共轭对称分量和共轭反对称分量:

x(n)=xep(n)+xop(n),x(n)=x_{ep}(n)+x_{op}(n),

其中

xep(n)=12[x(n)+x(Nn)],x_{ep}(n)=\frac{1}{2}[x(n)+x^*(N-n)], xop(n)=12[x(n)x(Nn)].x_{op}(n)=\frac{1}{2}[x(n)-x^*(N-n)].

同理,频域函数也可以分解成共轭对称分量和共轭反对称分量。

如果把序列分为实部和虚部:

x(n)=Re[x(n)]+jIm[x(n)],x(n)=Re[x(n)]+jIm[x(n)],

则有限长序列实部的 DFT 等于 X(k)X(k) 的共轭对称分量;有限长序列虚部乘 jj 后的 DFT 等于 X(k)X(k) 的共轭反对称分量。

反过来,如果把序列分成共轭对称分量和共轭反对称分量,则

DFT[xep(n)]=Re[X(k)],DFT[x_{ep}(n)]=Re[X(k)], DFT[xop(n)]=jIm[X(k)].DFT[x_{op}(n)]=jIm[X(k)].

实际中经常对实序列进行 DFT,利用这些对称性可以减少运算量。若 x(n)x(n) 是实序列,则

X(k)=X(Nk).X(k)=X^*(N-k).

x(n)x(n) 是纯虚序列,则

X(k)=X(Nk).X(k)=-X^*(N-k).

x(n)x(n) 是实偶序列,则 X(k)X(k) 是实偶对称;若 x(n)x(n) 是实奇序列,则 X(k)X(k) 是纯虚奇对称。

这一性质还能设计高效算法:构造

x(n)=x1(n)+jx2(n),x(n)=x_1(n)+jx_2(n),

只做一次 NN 点 DFT 得到 X(k)X(k),就可以分离出两个实序列的 DFT:

X1(k)=12[X(k)+X(Nk)],X_1(k)=\frac{1}{2}[X(k)+X^*(N-k)], X2(k)=j2[X(k)X(Nk)].X_2(k)=-\frac{j}{2}[X(k)-X^*(N-k)].

频域采样定理

若序列 x(n)x(n) 的长度为 MM,频域采样点数为 NN。只有当

NMN\ge M

时,才能由频域采样 X(k)X(k) 通过 IDFT 恢复原序列:

xN(n)=IDFT[X(k)]=x(n).x_N(n)=IDFT[X(k)]=x(n).

如果 N<MN<M,就会产生时域混叠。直观理解是:频域采样点太少,相当于时域周期太短,原序列不同部分会折叠到同一个主值区间内。

用 DFT 计算循环卷积

设有限长序列 x1(n)x_1(n)x2(n)x_2(n) 的长度分别为 NNMM,要求 LL 点循环卷积。步骤如下:

先计算

X1(k)=DFT[x1(n)]L,X_1(k)=DFT[x_1(n)]_L,

再计算

X2(k)=DFT[x2(n)]L,X_2(k)=DFT[x_2(n)]_L,

然后频域相乘:

Y(k)=X1(k)X2(k).Y(k)=X_1(k)X_2(k).

最后作 IDFT:

y(n)=IDFT[Y(k)]L.y(n)=IDFT[Y(k)]_L.

这样就把时域循环卷积转换成了频域乘法。

线性卷积与循环卷积的关系

x1(n)x_1(n) 长度为 NNx2(n)x_2(n) 长度为 MM。它们的线性卷积非零长度为

N+M1.N+M-1.

循环卷积可以看作线性卷积的周期延拓后再取主值序列。因此,如果循环长度 LL 太小,线性卷积的尾部会折叠回前面,产生混叠。

若希望 LL 点循环卷积等于线性卷积,必须满足

LN+M1.L\ge N+M-1.

课堂练习中的结论也正是这个条件。若 y(n)=x(n)h(n)y(n)=x(n)*h(n),而 w(n)w(n)LL 点循环卷积,则 w(n)=y(n)w(n)=y(n) 的条件是

LM+N1.L\ge M+N-1.

例如,对 6 点序列

{5,1,3,0,5,2}\{5,1,3,0,5,2\}

向左 2 点圆周移位,得到

{3,0,5,2,5,1}.\{3,0,5,2,5,1\}.

这类题本质上都是主值区间内“首尾相接”的移动。

DFT、循环卷积与线性系统输出

若线性时不变系统的单位脉冲响应和输入分别为

h(n)=R4(n),x(n)=R8(n),h(n)=R_4(n),\qquad x(n)=R_8(n),

则线性卷积输出为

y(n)={1,2,3,4,4,4,4,4,3,2,1}.y(n)=\{1,2,3,4,4,4,4,4,3,2,1\}.

若对 x(n)x(n)h(n)h(n) 作 12 点 DFT,得到 X(k)X(k)H(k)H(k),令

Y1(k)=H(k)X(k),Y_1(k)=H(k)X(k),

再作 IDFT,则得到

y1(n)={1,2,3,4,4,4,4,4,3,2,1,0}.y_1(n)=\{1,2,3,4,4,4,4,4,3,2,1,0\}.

这里取 12 点是安全的,因为线性卷积长度为

4+81=11,4+8-1=11,

12 点不会混叠,只是在末尾多了一个 0。

为什么需要 FFT

直接计算 DFT 时,每一个 X(k)X(k) 都需要对 n=0N1n=0\sim N-1 求和。计算全部 NNX(k)X(k),需要约

N2N^2

次复数乘法和

N(N1)N(N-1)

次复数加法。当 N=1024N=1024 时,复数乘法数量约为 1048576 次,实时处理会非常困难。

FFT 的基本思想是:不断把长序列的 DFT 分解成短序列的 DFT,并利用旋转因子的周期性、对称性和可约性减少重复计算。FFT 不是另一种变换,它是 DFT 的快速算法。

DFT与FFT运算量比较

FFT 主要分两类:按时间抽取 DIT-FFT 和按频率抽取 DIF-FFT。

基2 DIT-FFT:按时间抽取

基2 FFT 要求

N=2M.N=2^M.

若长度不足,可以补零到 2 的整数次幂。

DIT 的第一步是按时域下标奇偶分组:

x1(r)=x(2r),r=0,1,,N21,x_1(r)=x(2r),\qquad r=0,1,\cdots,\frac{N}{2}-1, x2(r)=x(2r+1),r=0,1,,N21.x_2(r)=x(2r+1),\qquad r=0,1,\cdots,\frac{N}{2}-1.

将 DFT 按偶数项和奇数项拆开,可得

X(k)=X1(k)+WNkX2(k),X(k)=X_1(k)+W_N^kX_2(k),

其中 X1(k)X_1(k)X2(k)X_2(k) 是两个 N/2N/2 点 DFT。

后一半频谱由同一组 X1(k),X2(k)X_1(k),X_2(k) 得到:

X(k+N2)=X1(k)WNkX2(k).X\left(k+\frac{N}{2}\right)=X_1(k)-W_N^kX_2(k).

这两个式子构成蝶形运算:

{X(k)=X1(k)+WNkX2(k),X(k+N/2)=X1(k)WNkX2(k).\begin{cases} X(k)=X_1(k)+W_N^kX_2(k),\\ X(k+N/2)=X_1(k)-W_N^kX_2(k). \end{cases}

FFT蝶形运算

N=8N=8 为例,先把序列分成偶数下标和奇数下标:

x1={x(0),x(2),x(4),x(6)},x_1=\{x(0),x(2),x(4),x(6)\}, x2={x(1),x(3),x(5),x(7)}.x_2=\{x(1),x(3),x(5),x(7)\}.

分别作 4 点 DFT 后,再用蝶形运算合成 X(0)X(7)X(0)\sim X(7)。进一步分解时,4 点 DFT 又可以拆成 2 点 DFT,直到每个子问题只剩下 2 点 DFT。

DIT-FFT分解树

FFT 的运算量

N=2MN=2^M,则基2 FFT 共有 M=log2NM=\log_2N 级蝶形运算。每一级有 N/2N/2 个蝶形。每个蝶形需要 1 次复乘和 2 次复加。

因此,FFT 的复乘次数为

N2log2N,\frac{N}{2}\log_2N,

复加次数为

Nlog2N.N\log_2N.

直接 DFT 的复乘次数是 N2N^2,复加次数是 N(N1)N(N-1)。当 NN 很大时,FFT 的优势非常明显。

课堂练习中,一个蝶形运算需要 1 次复数乘法和 2 次复数加法;对于 N=2MN=2^M 的基2 DIT-FFT,总共需要

MN2\frac{MN}{2}

次复数乘法和

MNMN

次复数加法。

原位运算与旋转因子规律

FFT 的程序实现中常用原位运算。所谓原位运算,就是输入数据、中间结果和最后输出都使用同一组存储单元。原因是每个蝶形的两个输入在计算完两个输出后不再需要保留,因此可以直接覆盖。

若第 LL 级运算中,蝶形两个输入相距

B=2L1,B=2^{L-1},

则可以用形式化表达:

AL(J)AL1(J)+AL1(J+B)WNp,A_L(J)\leftarrow A_{L-1}(J)+A_{L-1}(J+B)W_N^p, AL(J+B)AL1(J)AL1(J+B)WNp.A_L(J+B)\leftarrow A_{L-1}(J)-A_{L-1}(J+B)W_N^p.

其中第 LL 级的旋转因子指数满足

p=J2ML,J=0,1,,2L11.p=J\cdot 2^{M-L},\qquad J=0,1,\cdots,2^{L-1}-1.

程序设计时,通常用级数 LL 作最外层循环,再安排每一级中的蝶形组和旋转因子。

倒位序输入

DIT-FFT 的一个特点是:输入通常按倒位序排列,输出按自然顺序排列。

N=8N=8 为例,自然顺序与倒位序如下。

倒位序表

自然顺序

0,1,2,3,4,5,6,70,1,2,3,4,5,6,7

对应二进制

000,001,010,011,100,101,110,111.000,001,010,011,100,101,110,111.

把二进制位反过来,就得到倒位序

0,4,2,6,1,5,3,7.0,4,2,6,1,5,3,7.

因此,DIT-FFT 的输入次序常写成

x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7).x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7).

形成这种顺序的原因,是算法不断按下标奇偶抽取。

DIF-FFT:按频率抽取

DIF-FFT 与 DIT-FFT 思想相似,但分解对象不同。DIF 先把 DFT 的输出 X(k)X(k)kk 的奇偶分组。

从 DFT 定义出发,可以把序列前半段和后半段配对:

x(n),x(n+N2).x(n),\qquad x\left(n+\frac{N}{2}\right).

利用

WNN/2=1,W_N^{N/2}=-1,

可构造两个 N/2N/2 点序列:

x1(n)=x(n)+x(n+N2),x_1(n)=x(n)+x\left(n+\frac{N}{2}\right), x2(n)=[x(n)x(n+N2)]WNn.x_2(n)=\left[x(n)-x\left(n+\frac{N}{2}\right)\right]W_N^n.

其中 x1(n)x_1(n)N/2N/2 点 DFT 给出偶数频率点,x2(n)x_2(n)N/2N/2 点 DFT 给出奇数频率点。

DIF 的流程是:先蝶形运算,后分解成小 DFT;DIT 的流程是:先分解成小 DFT,后蝶形合成。

DIT与DIF比较

二者相同点是都能原位运算,运算量也相同:复乘约 (N/2)log2N(N/2)\log_2N,复加约 Nlog2NN\log_2N。不同点是:DIT 输入为倒位序、输出为自然顺序;DIF 输入为自然顺序、输出为倒位序。

IDFT 如何用 FFT 实现

IDFT 定义为

x(n)=1Nk=0N1X(k)WNnk.x(n)=\frac{1}{N}\sum_{k=0}^{N-1}X(k)W_N^{-nk}.

与 DFT 相比,区别是旋转因子指数符号相反,并且多了系数 1/N1/N。因此可以直接修改 FFT 中的旋转因子符号,再乘 1/N1/N

另一种更常见的方法是利用共轭:

先将 X(k)X(k) 取共轭;再对 X(k)X^*(k) 作 FFT;最后把结果再取共轭并乘以 1/N1/N。公式可写成

x(n)=1N[DFT(X(k))].x(n)=\frac{1}{N}\left[DFT(X^*(k))\right]^*.

这个方法常用于工程编程,因为可以复用已有 FFT 程序。

FFT 的应用:实数序列的高效计算

前面讨论的 FFT 默认是复数序列。但现实信号大多数是实信号。如果直接把实序列当成虚部为 0 的复序列,会浪费一半存储和一部分计算量。

一种高效方法是:用一次 NN 点复数 FFT 同时计算两个 NN 点实序列的 DFT。

x(n)x(n)y(n)y(n) 是两个实序列,构造复序列

g(n)=x(n)+jy(n).g(n)=x(n)+jy(n).

g(n)g(n) 作一次 FFT,得到 G(k)G(k)。再利用 DFT 的共轭对称性,把 X(k)X(k)Y(k)Y(k) 分离出来。这样一次复数 FFT 就可以完成两个实序列的 DFT,效率提高接近一倍。

实序列FFT技巧

这也是为什么理解共轭对称性不仅是为了做题,更是为了写高效算法。

结语

这一讲的核心可以归纳为四条线索。

第一,DFT 是为了让计算机能够处理频谱。它把 DTFT 的连续频谱变成有限个等间隔频率样本。

第二,DFT 隐含周期性。有限长序列在 DFT 中被看成周期序列的一个主值区间,因此会出现圆周移位、循环卷积、取模下标等概念。

第三,循环卷积与线性卷积不同。循环卷积是线性卷积周期延拓后的主值序列;只有当 LN+M1L\ge N+M-1 时,循环卷积才等于线性卷积。

第四,FFT 是 DFT 的快速算法。它不是新的变换,而是利用旋转因子的周期性、对称性和可约性,把 N2N^2 级别的计算量降到 Nlog2NN\log_2N 级别。

理解 DFT 和 FFT 的关键不是死记公式,而是抓住一句话:DFT 把有限长序列看作周期序列的一个周期,在频域上取有限个样本;FFT 则利用这种结构进行快速计算。

Back to archive

Discussion

Comments

Post

Share questions, corrections, or extra notes about this post.