数字信号处理之快速傅里叶变换 2021-10-24

数字信号处理

---

前言

本章主要讲述快速傅里叶变换。

一、概述

DFT的面临的最大问题就是运算量过大，导致无法在实时性要求高的通信系统中发挥作用。
求出N点$X(k)需要N^2$次复数乘法及N(N-1)次复数加法。实现一次复数乘法需要四次实数乘两次实数加，实现一次复数加法需要2次实数加法。当N很大时，计算量非常大。对于2D的DFT计算量更是惊人。因此想要将DFT的应用到实际通信系统中，需要更简单快速的算法。FFT便应运而生了。

DFT的运算过程中实际包括了大量的重复运算。如何利用这一特性称为实现FFT的关键。1965年，J.W.Cooley和J.W.Tukey首先提出了一种高效的算法，将DFT的计算量由$N^2次降为\frac{N}{2}lo{g}_{2}N次。N为1024时，计算量降为普通DFT的4.88$
自从Cooley-Tukey的算法出现后，新算法不断涌现。总的来说由两个主要方向：一是N等于2的整数次幂的算法，如基2、基4算法等；二是N不等于2的整数次幂的算法，它是以Winograd为代表的一类算法。（素因子算法，Winograd算法–WFTA）
split-radix算法是当前第一类算法中最理想的一种。
WFTA算法所需的乘法次数较第一类算法少，但理论复杂，编程困难，未来也许会有发展前景。

二、时间抽取基2的FFT算法和频域抽取基2的FFT算法

N点序列x(n)，N=2^M，M为正整数。将x(n)按奇偶分成两组，n=2r及n=2r+1，r = 0,1,… , n/2-1，于是：
$$X(k)=\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}$$

$$X(k)=A(k)+W_N^{k}B(k)$$
式中：
$$A(k)=\sum _{r=0}^{N/2-1}x(2r)W_{N/2}^{rk}，B(k)=\sum _{r=0}^{N/2-1}x(2r+1)W_{N/2}^{rk}，W_{N/2}=e^{-j4\pi /N}$$
所以用A(k)和B(k)可以完整的表示X(k)。A(k)和B(k)仍是高复合数(N/2)的DFT，可以按照上述方法继续分解，直到变为2点的DFT。所以N点DFT需要分解为M级，进行计算。这种将时间按照奇偶分开的方法称为时间抽取法(Decimation in time,DIT)

类似的有频域抽取算法（DIF）
无论是DIF还是DIT都可以看作是将N×N的W矩阵做分解来实现。

三、分裂基FFT算法

分裂基算法（split-radix）的实际是2基和4基的混合算法。N≥64时，分裂基算法需要的实数乘的次数约是2基算法的2/3。

四、线性调频Z变换（CZT）

$$X(Z_r)=CZT[X(n)]=\sum _{n=0}^{\infty }x(n)z_r^{-n}=\sum _{n=0}^{\infty }x(n)A^{-n}{W}^{nr}$$
当r=0时，$z_0=A_0{e}^{j{\theta }_0}，该点在Z平面上的幅度为A_0，相位为{\theta }_0$，是CZT的起点。
当r=1时，$z_1=A_0{W}_0^{-1}{e}^{j{\theta }_0+{\phi }_0}，该点在Z平面上的幅度变为为A_0{W}_0^{-1}，相位增加为{\phi }_0，z_0,z_1,z_2,...构成了CZT变换的路径。当A_0大于1时，在单位圆外，反之在单位圆内；当W_0大于1时，螺旋线内旋，反之螺旋线外旋。当A_0=1，W_0=1时，CZT的变换路径为单位圆上的一段弧；当A_0=1，W_0=1，{\theta }_0=0，M=N时，CZT变成普通DFT。$

五、离散时间系统的相位、结构与状态变量描述

1、离散时间系统的相频响应

定义系统的相位延迟(Phase Delay)：它表示的是输出相对于输入的时间延迟。
$${\tau }_p(\omega )=-\frac{\phi (\omega )}{\omega }$$
定义系统的群延迟(Group Delayl)为：
$${\tau }_g(\omega )=-\frac{d\phi (\omega )}{d\omega }$$
它表示的是信号中各个频率延迟的关系。如果系统具有线性相位，则它的群延迟是一个常数。如果群延迟和频率有关则系统不具有线性相位。也就是宽带信号经过系统后各个频率延迟不同，从而导致信号畸变。

2、FIR系统的线性相位特性

（1）如果FIR系统的单位冲击响应满足：$h(n)=\pm h(N-1-n)$，也就是系统满足对称性，则该系统具有线性相位。式中存在奇偶两种对称可能，N也有奇数、偶数的可能。对于奇对称情况，相当于在滤波器前先经过了一个90°移相器，幅频特性类似于差分器，希尔伯特变换器和差分器的单位抽样响应一般也是奇对称的。因此，一般用途的滤波器通常选用偶对称形式，长度N也往往取奇数。
$\bullet$类型I滤波器：h(n)为偶对称、N为奇数
$\bullet$类型II滤波器：h(n)为偶对称、N为偶数
$\bullet$类型III滤波器：h(n)为奇对称、N为奇数
$\bullet$类型IV滤波器：h(n)为奇对称、N为偶数
（2）由(1)的条件，可得到其对应的Z变换：
$$H(Z)=\pm z^{-(N-1)}H(z^{-1})$$
由上式可以看出，H(Z)的零点也是$H(z^{-1})$的零点，反之亦然。设$H(z)$的一个零点$Z_k=r_ke^{j\varphi_k}:
$\bullet 当z_k在单位圆内时，H(Z)是四阶系统，有四个对称的零点：z_k,z_k^{-1},(z_k^{\ast }{)}^{-1},z_k^{\ast }$
$\bullet 当z_k在实轴上时，H(Z)是二阶系统，有2个镜像对称的零点：z_k=r_k,z_k^{-1}=1/r_k$
$\bullet 当z_k在单位圆上时，H(Z)是二阶系统，有2个共轭对称的零点：z_k=e^{j\phi },z_k^{\ast })=e^{-j\phi }$
$\bullet 当z_k在单位圆外时，H(Z)是一阶系统，没有对称的零点。$

3、全通系统和最小相位系统

（1）全通系统：因果系统的幅频响应对所有频率都等于1或者常数，则此系统成为全通系统。
最简单的一阶全通系统：$H_{ap}(z)=z^{-k}$，另一个简单的全通系统$H_{ap}(z)=\frac{1-{\lambda }^{-1}{z}^{-1}}{1-\lambda z^{-1}}，|\lambda |<1$
一个N阶全通系统可表示为：
$$H_{ap}(z)=\pm \prod _{k=1}^N[\frac{z^{-1}-{\lambda }_k^{\ast }}{1-{\lambda }_k{z}^{-1}}]，|\lambda |<1$$
上式经过变换可表示为：
$$h_{ap}(z)=\pm \frac{z^{-N}A(z^{-1})}{A(z)}$$
由上式可以看出，全通系统的分子多项式和分母多项式是互为镜像的多项式。
$\bullet$全通系统一般是IIR系统；
$\bullet$全通系统的零、极点数相同
$\bullet$极点在单位圆内，零点在单位圆外
$\bullet$极点与零点以单位圆镜像对称
$\bullet$当$\omega 由0变到\pi 时，相频响应\phi (\omega )单调递减$
$\bullet$全通系统的群延迟始终为正值。
（2）最小相位系统：如果一个离散系统的极点和零点都在单位圆内，则称系统为最小相位系统；如果零点都在单位圆外，则称系统为最大相位系统；
性质1：在一组具有相同幅频响应的因果的且稳定的滤波器集合中，最小相位滤波器具有对于$\omega$轴（即0相位）具有最小的相位偏移；
性质2：$h_{min}(n)$是最小相位系统的单位抽样响应，则：
$$\sum _{n=0}^M{h}_{min}(n)⩾\sum _{n=0}^{M}h(n)=E(M)，0⩽M<\infty$$
式中，，E(M)为M的函数，表示单位抽样响应在M个点的累计能量。上式说明，最小相位系统的单位抽样响应能力集中在n为较小值的范围内，即最小相位离散系统的抽样响应具有最小的延迟。因此$h_{min(n)}$也叫做最小延迟序列。
性质3：最小相位系统H(z)可表示为$H(z)=|H(z)|e^{jarg[H(z)]}，令\check{H}=ln(H(z)=ln|H(z)|+jarg[H(z)]$
则${\check{H}}_R(e^{j\omega })和{\check{H}}_I(e^{j\omega })$构成一对希尔伯特变换，即：
$${\check{H}}_R(e^{j\omega })=\check{H}(0)+\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{\check{H}}_I(e^{j\omega })cot(\frac{\omega -\theta }{2})d\theta$$
$${\check{H}}_I(e^{j\omega })=-\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{\check{H}}_R(e^{j\omega })cot(\frac{\omega -\theta }{2})d\theta$$
性质4：给定一个稳定的、因果的系统H(z)=D(z)/N(z)，定义其逆滤波器$H_{IV}(z)=1/H(z)=D(z)/N(z)$
当且仅当H(z)是最小相位系统时，$H_{IV}(z))$才是稳定的、因果的，也就是物理可实现的。
性质5：任何一个非最小相位系统都可以由一个最小相位系统和一个全通系统级联而成。
$$H(z)=H_{min}(z)H_{ap}(z)=H_{min}(z)\frac{z^{-1}-z_0}{1-z_0^{\ast }{z}^{-1}}$$
上式中后面部分的$H_{min}(z)$是与H(z)具有相同的幅频响应，这给我们提供了一种能够由非最小相位系统构成最小相位系统的方法。

4.谱分解

重要性质：一个线性相位系统可以有两种级联形式：
$$P(z)=H(z)H(z^{-1})或P(z)=H_0(z)H_1(z)$$
对于前者，零点的位置和数量决定是否H(z)和H(1/Z)具有相同的幅频响应。
对于后者，如果把单位圆内的零点都赋给$H_0(z)$，则其实最小相位的，而 $H_1(z)$就是最大相位的。而如果不这样分配则它们都是混合相位的。

5.FIR结构

FIR系统可由差分方程或转移函数来描述：
$$y(n)=\sum _{r=0}^{M}b(r)x(n-r)orH(z)=\sum _{r=0}^{M}b(r)z^{-r}$$
FIR系统可以直接实现，也可以级联实现。


线性相位的FIR系统的结构：

FIR系统的频率抽样实现：
$$H(z)=\sum _{n=0}^{N-1}h(n)Z^{-n},H(k)=\sum _{n=0}^{N-1}h(n)W_N^{-nk}$$
H(k)实际上是$H(z){|}_{z=e^{j\omega }}=H(e^{j\omega })$在单位圆上的N个值，即H(k)实际上是$H(e^{j\omega })$在频域上的抽样，因此可以用H(k)来表示H(z)。
$$H(z)=\frac{1-z^{-N}}{N}\sum _{n=0}^{N-1}\frac{H(k)}{1-exp(j\frac{2\pi }{N}k)z^{-1}}$$

这种结构中，整个系统是N个$H_{2,k}$并联后再和$H_1(z)$级联的结果。$H_1(z)$有N个零点，分别抵消掉$H_{2,k}$的一个极点，这样总的转移函数是N个(N-1)阶$z^{-1}$多项式的和，仍是一个FIR系统。
$$H_{2,k}(z)=\frac{H(k)}{1-exp(j\frac{2\pi }{N}k)z^{-1}}=\frac{H(k)}{1-W_N^{-k}{z}^{-1}}$$
$$H_1(z)=\frac{1-z^{-N}}{N}$$
用硬件实现上图的结构时由于误差存在，$H_{2,k}$的极点和$H_1(z)$的零点不一定完全抵消，这样在单位圆上存在极点将使系统变得不稳定。为了客服这个问题，通过在半径为r（r<1,但接近于1）的圆上对H(z)抽样来实现。即用$r{z}^{-1}$来代替$z^{-1}$可以消除不稳定问题。

---

总结

本节主要学习了FFT，CZT，离散时间系统的相位、结构与状态变量描述。
FFT主要是了解一下基本算法，最后离散时间系统的相位结构状态描述是比较关键的，对于后续滤波器的设计有重要知道意义。