本文是“离散时间信号处理”系列的第 09 章，主题为“离散傅里叶变换的计算”。

上一篇：离散傅里叶变换 · 下一篇：利用离散傅里叶变换的信号傅里叶分析

离散傅里叶变换的计算

离散傅里叶变换的高效计算

直接计算DFT的运算次数
\[ \begin{aligned} X[k]=&\sum_{n=0}^{N-1}x[n]W_N^{kn},0\leq k\leq N-1\\ =&\sum_{n=0}^{N-1}[(\mathcal{R}e\{x[n]\}\mathcal{R}e\{W_N^{kn}\}-\mathcal{I}m\{x[n]\}\mathcal{I}m\{W_N^{kn}\})+\\ &j(\mathcal{R}e\{x[n]\}\mathcal{I}m\{W_N^{kn}\}+\mathcal{I}m\{x[n]\}\mathcal{R}e\{W_N^{kn}\})],k=0,1,\cdots,N-1 \end{aligned} \]

(1) 计算DFT的每一个值需要\(N\)次复数乘法和\((N-1)\)次复数加法.

(2) 计算DFT的\(N\)个值需要\(N^2\)次复数乘法和\(N(N-1)\)次复数加法.

(3) 每个复数乘法需要4次实数乘法和2次实数加法, 每个复数加法需要2次实数加法.

(4) 序列\(x[n]\)的傅里叶变换的直接计算法需要\(4N^2\)次实数乘法和\(N(4N-2)\)次实数加法.

按时间抽取的FFT算法

基2-DIT-FFT算法

将\(x[n]\)分解为偶数点和奇数点, 可得
\[ \begin{aligned} X[k]&=\sum_{n\text{为偶数}}x[n]W_N^{nk}+\sum_{n\text{为奇数}}x[n]W_N^{nk}\\ &=\sum_{r=0}^{(N/2)-1}x[2r]W_N^{2rk}+\sum_{r=0}^{(N/2)-1}x[2r+1]W_N^{(2r+1)k}\\ &=\sum_{r=0}^{(N/2)-1}x[2r]W_{N/2}^{rk}+W_N^k\sum_{r=0}^{(N/2)-1}x[2r+1]W_{N/2}^{rk}\\ &=G[k]+W_N^kH[k],k=0,1,\cdots,N-1\\ \end{aligned} \]

其中\(G[k]\)是原序列偶数点的\(\dfrac{N}{2}\)点DFT, \(H[k]\)是原序列奇数点的\(\dfrac{N}{2}\)点DFT.

将上式化简后得到\(N\)点DFT的表示为
\[ \begin{cases} X[k]=G[k]+W_N^kH[k]\\[2pt] X\left[k+\dfrac{N}{2}\right]=G[k]-W_N^kH[k] \end{cases},k=0,1,\cdots,\frac{N}{2}-1 \]

(1) 对于没有利用对称性的直接计算, 需要\(N^2\)次复数乘法和加法(\(N-1\)近似为\(N\)).

(2) 分解成两个\(\dfrac{N}{2}\)点DFT计算, 需要\(\left(N+\dfrac{N^2}{2}\right)\)次复数乘法和加法.

把一个 \(N\) 点 DFT 计算分解成两个 \(\dfrac{N}{2}\) 点 DFT 计算的按时间抽取流图

进一步分解\(\dfrac{N}{2}\)点DFT直至将计算简化为2点DFT的计算

按时间抽取将8点DFT计算完全分解的流图

2点DFT的流图

(1) 对于\(N\)点的DFT, 可以分解成\(\upsilon=\log_2N\)级的运算.

(2) 每一级的蝶形计算次数为\(\dfrac{N}{2}\), 因此总的蝶形计算次数为\(\dfrac{N\log_2N}{2}\).

(3) 每次蝶形计算需要2次复数乘法和加法, 因此总共需要\(N\log_2N\)次复数乘法和加法.

利用系数\(W_N^r\)的对称性和周期性进一步减少计算量
\[ W_N^{r+N/2}=W_N^{N/2}W_N^r=-W_N^r \]

(a) 基本蝶形计算流图(b) 只需要1次复数乘法的简化蝶形计算流图

化简后每次蝶形计算需要1次复数乘法和2次复数加法, 因此总共需要\(\dfrac{N\log_2N}{2}\)次复数乘法和\(N\log_2N\)次复数加法.

基2-DIT-FFT算法实现8点DFT的最终流图

同址运算

同址运算的流图

(1) 若将\(X_m[p]\)和\(X_m[q]\)分别存放在原存放\(X_{m-1}[p]\)和\(X_{m-1}[q]\)的同一存储寄存器中, 每一行节点共用一个寄存器, 则实现全部计算只需要一列存储\(N\)个复数的寄存器, 这种计算称为同址计算.

(2) 同址运算必须满足蝶形计算结构, 1个\(N\)点的基2-DIT-FFT结构需要\(N\)个寄存器.

倒位序

为了实现同址运算, 输入数据的存储和读取需要按照倒位序存储.对于8点流图, 用二进制形式来标注整个数据, 则可以得到

\(X_0[000]=x[000]\)	\(X_0[100]=x[001]\)
\(X_0[001]=x[100]\)	\(X_0[101]=x[101]\)
\(X_0[010]=x[010]\)	\(X_0[110]=x[011]\)
\(X_0[011]=x[110]\)	\(X_0[111]=x[111]\)

若\((n_2,n_1,n_0)\)为序列\(x[n]\)中标号的二进制表示, 则序列值\(x[n_2,n_1,n_0]\)存放在数列\(X_0[n_0,n_1,n_2]\)的位置上.

按频率抽取的FFT算法

基2-DIF-FFT算法

将\(X[k]\)分解为偶序号频率样本和奇序号频率样本, 可得
\[ \begin{aligned} X[2r]&=\sum_{n=0}^{N-1}x[n]W_N^{2rn}\\ &=\sum_{n=0}^{(N/2)-1}x[n]W_N^{2rn}+\sum_{n=N/2}^{N-1}x[n]W_N^{2rn}\\ &=\sum_{n=0}^{(N/2)-1}x[n]W_N^{2rn}+\sum_{n=0}^{(N/2)-1}x\left[n+\dfrac{N}{2}\right]W_{N}^{2r(n+N/2)}\\ &=\sum_{n=0}^{(N/2)-1}\left(x[n]+x\left[n+\dfrac{N}{2}\right]\right)W_{N/2}^{rn},r=0,1,\cdots,\dfrac{N}{2}-1 \end{aligned} \]
\[ \begin{aligned} X[2r+1]&=\sum_{n=0}^{N-1}x[n]W_N^{n(2r+1)}\\ &=\sum_{n=0}^{(N/2)-1}x[n]W_N^{n(2r+1)}+\sum_{n=N/2}^{N-1}x[n]W_N^{n(2r+1)}\\ &=\sum_{n=0}^{(N/2)-1}x[n]W_N^{n(2r+1)}+\sum_{n=0}^{(N/2)-1}x\left[n+\dfrac{N}{2}\right]W_{N}^{[n+(N/2)](2r+1)}\\ &=\sum_{n=0}^{(N/2)-1}\left(x[n]-x\left[n+\dfrac{N}{2}\right]\right)W_N^nW_{N/2}^{nr},r=0,1,\cdots,\dfrac{N}{2}-1 \end{aligned} \]

按频率抽取分解的2点DFT的流图

8点DFT计算按频率抽取完全分解的流图

按频率抽取的流图和按时间抽取的流图互为转置关系.

补充知识点

用1次\(N\)点FFT算法来计算2个\(N\)点实序列DFT的方法

已知\(N\)点实序列\(x_1[n]\)和\(x_2[n]\), 为了只用1次FFT算法得到这两个序列的DFT, 令\(x[n]=x_1[n]+jx_2[n]\), 经过FFT计算后得到的序列为\(X[k]\), 则有
\[ \begin{aligned} x_1[n]=\frac{x[n]+x^*[n]}{2}\Leftrightarrow X_1[k]&=\frac{X[k]+X^*[((-k))_N]}{2}\\ &=X_{ep}[k]\\ &=\frac{1}{2}(\mathcal{R}e\{X[k]\}+\mathcal{R}e\{X[N-k]\})+\frac{j}{2}(\mathcal{I}m\{X[k]\}-\mathcal{I}m\{X[N-k]\}) \end{aligned} \]
\[ \begin{aligned} x_2[n]=\frac{x[n]-x^*[n]}{2j}\Leftrightarrow X_2[k]&=\frac{X[k]-X^*[((-k))_N]}{2j}\\& =-jX_{op}[k]\\ &=\frac{1}{2}(\mathcal{I}m\{X[k]\}+\mathcal{I}m\{X[N-k]\})-\frac{j}{2}(\mathcal{R}e\{X[k]\}-\mathcal{R}e\{X[N-k]\}) \end{aligned} \]

(扩展一) 可以将1个\(2N\)点实序列分解成2个\(N\)点实序列(偶数点和奇数点), 则得到用1次\(N\)点FFT算法来计算1个\(2N\)点实序列DFT的方法, 其中
\[ x_1[n]=x[2n]\quad\quad x_2[n]=x[2n+1] \]

构造\(y[n]=x_1[n]+jx_2[n]\), 其对应的DFT为\(Y[k]\), 则\(X_1[k]\)和\(X_2[k]\)分别可表示为
\[ X_1[k]=\dfrac{1}{2}\{Y[k]+Y^*[N-k]\} \]
\[ X_2[k]=\dfrac{1}{2j}\{Y[k]-Y^*[N-k]\} \]

因此\(X[k]\)的\(2N\)点DFT可表示为
\[ \begin{cases} X[k]=X_1[k]+W_{2N}^kX_2[k]\\ X[k+N]=X_1[k]-W_{2N}^kX_2[k] \end{cases},k=0,1,\cdots,N-1 \]

(扩展二) 可以将实序列\(x_1[n]\)和\(x_2[n]\)进一步分解成对称和反对称序列, 则得到用1次\(N\)点FFT算法来计算4个\(N\)点实对称或反对称序列DFT的方法, 其中
\[ x_1[n]=y_1[n]+y_3[n]\quad\quad x_2[n]=y_2[n]+y_4[n] \]

则4个DFT分别可表示为
\[ y_1[n]=y_1[N-n]\Leftrightarrow Y_1[k]=\mathcal{R}e\{X_1[k]\}=\frac{1}{2}(\mathcal{R}e\{X[k]\}+\mathcal{R}e\{X[N-k]\}) \]
\[ y_2[n]=y_2[N-n]\Leftrightarrow Y_2[k]=\mathcal{R}e\{X_2[k]\}=\frac{1}{2}(\mathcal{I}m\{X[k]\}+\mathcal{I}m\{X[N-k]\}) \]
\[ y_3[n]=-y_3[N-n]\Leftrightarrow Y_3[k]=j\mathcal{I}m\{X_1[k]\}=\frac{j}{2}(\mathcal{I}m\{X[k]\}-\mathcal{I}m\{X[N-k]\}) \]
\[ y_4[n]=-y_4[N-n]\Leftrightarrow Y_4[k]=j\mathcal{I}m\{X_2[k]\}=\frac{-j}{2}(\mathcal{R}e\{X[k]\}-\mathcal{R}e\{X[N-k]\}) \]

已知一个实对称序列\(x[n]\), 考虑序列\(u[n]\)
\[ u[n]=x[((n+1))_N]-x[((n-1))_N] \]

\(u[n]\)满足以下关系
\[ \begin{aligned} u[N-n]&=x[((N-n+1))_N]-x[((N-n-1))_N]\\ &=x[((n-1))_N]-x[((n+1))_N]\\ &=-u[n] \end{aligned} \]

因此\(u[n]\)是一个反对称序列.

用FFT算法来计算IDFT的方法

(方法一) 已知IDFT的计算公式为
\[ x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]W_N^{-kn},0\leq n\leq N-1 \]

等式两边取共轭
\[ x^*[n]=\frac{1}{N}\sum_{k=0}^{N-1}X^*[k]W_N^{kn}\Rightarrow Nx^*[n]=\sum_{k=0}^{N-1}X^*[k]W_N^{kn} \]

因此通过FFT算法来计算IDFT的流程为
\[ X[k]\xrightarrow{\text{共轭}}X^*[k]\xrightarrow{\text{FFT}}Nx^*[n]\xrightarrow{\frac{1}{N}}x^*[n]\xrightarrow{\text{共轭}}x[n] \]

(方法二) 直接将\(X[k]\)输入FFT得到的输出为
\[ g[n]=\sum_{k=0}^{N-1}X[k]W_N^{kn} \]

\(g[n]\)进行翻转后的DFT为
\[ \begin{aligned} g[((N-n))_N]&=\sum_{k=0}^{N-1}X[k]W_N^{k(N-n)}\\ &=\sum_{k=0}^{N-1}X[k]W_N^{-kn} \end{aligned} \]

因此通过FFT算法来计算IDFT的流程为
\[ X[k]\xrightarrow{\text{FFT}}g[n]\xrightarrow{\text{翻转}}g[((N-n))_N]\xrightarrow{\frac{1}{N}}x[n] \]

旋转因子\(W_N^{kn}\)的性质

(1) 对称性
\[ (W_N^{kn})^*=W_N^{-kn}=W_N^{k(N-n)} \]

(2) 周期性
\[ W_N^{kn}=W_N^{k(n+N)}=W_N^{(k+N)n}\quad\text{(}n\text{和}k\text{均以}N\text{为周期)} \]

(3) 可约性
\[ W_N^{kn}=W_{N/2}^{kn/2}=W_{2N}^{2kn} \]

(4) 特殊值
\[ W_N^0=1\quad\quad W_N^{N/2}=W_2^1=-1 \]
\[ W_N^{N/4}=W_4^1=-j\quad\quad W_N^{k+N/2}=-W_N^k \]

---

上一篇：离散傅里叶变换 · 下一篇：利用离散傅里叶变换的信号傅里叶分析