第4章 快速傅里叶变换(FFT)-补充

第4章快速傅里叶变换补充

N为复合数的FFT算法4.6N为复合数的FFT算法上面讨论的是序列的点数N为2的幂次(即N=2M)情况下，按时间抽取和按频率抽取的基-2FFT算法的基本原理。这种基-2FFT算法在实际中使用得最多，因为它的程序简单，效率高，使用方便。但实际上无法保证总是处理长度为2的整数幂次的序列。若不满足N=2M，处理方法有以下几种:

（1）可将x(n)增补一些零值点，以使N增长到最邻近的一个2M数值。有限长序列补零之后，并不影响其频谱X(ejω)，只不过其频谱的抽样点数增加了，所造成的结果是增加了计算量而已。但是，有时计算量增加太多，浪费较大。例如,x(n)的点数N=300，则须补到N=29=512，要补212个零值点，因而人们才研究N≠2M时的FFT算法。有限长度序列补零后并不影响其频谱x(ejw)，只是频谱的采样点数增加了,当然，采样点的位置也有相应的变化。上例中由30点增加到32点，采样点则由2π/30的整数倍变为2π/32的整数倍，对应的绝对频率也由fs/30的整数倍变为fs/32的整数倍。在许多场合这种处理是可接受的，不适应对某些特殊频率点有特别要求的场合。

（2）若N是一个复合数，即它可以分解成一些因子的乘积，则可以用FFT的一般算法，即混合基FFT算法，如库利-图基（CooleyTukey）算法，而基-2算法只是这种一般算法的特例。总之，不管采用什么方法，计算DFT的高效算法是把计算长度为N的序列的DFT逐次分解成计算长度较短序列的DFT。这是很多高效算法的标准方法和基本原理。（3）若N为一个素数，要求准确的求出N点DFT可用线性调频Z变换（Chirp-Z变换）算法一、整数的多基多进制表示形式1.二进制（）数n<二进制数倒位序2.r进制（）数n<r进制数倒位序3.多基多进制(混合基）数n<二进制数倒位序Bb11-14二、的快速算法原理k=0,1,…,N-1N为一个复合数2、利用3、4、利用5、利用1、计算量4.7基－4FFT算法 在的特殊情况下,即当时,混合基算法变成基-4FFT算法步骤当时,采用分裂基算法可进一步降低运算量,4.8分裂基FFT算法n为偶数部分直接采用基---2FFT算法而n为奇数部分则用基-4FFT算法.接着N/2点DFT可以按以上同样办法分解成三部分(一个N/4点,两个N/8点DFT)第一次的两个N/4点DFT也分别可同样分成三个部分(每个N/4点可分成一个N/8点,两个N/16点DFT)如此进行到不能再分解为止.说明推导分裂基FFT时,是按n的奇偶划分的,所以其输入信号排序是基--2算法的倒位序.分裂基的基本蝶形运算:(倒L形)运算量：一个分裂基蝶形有两个复乘分裂基蝶形个数Xl11-134.9线性调频Z变换（Chirp-Z变换）算法利用FFT算法，可以很快地计算出有限长序列的DFT值，也即是Z变换在单位圆上的全部等间隔采样值。然而，在许多场合，并不一定需要计算全部频谱值，而仅需对某一频带内的信号频谱作较密集的分析。另外，采样也不一定局限于单位圆上，而需要计算出某一螺旋线上的等角度间隔的采样值。例如，在语音信号分析时K采用靠近语音信号序列Z变换的极点的螺旋线上进行采样，可以使语音信号的共振峰变得更尖锐，便于精确确定共振峰频率。Chirp-Z变换算法，在螺旋线上进行采样，可适用于更一般情况下，由x(n)求的快速算法，这种变换称为线性调频Z变换(简称CZT).图4-17单位圆与非单位圆采样（a）沿单位圆采样;(b)沿AB弧采样4.6.1算法基本原理

已知x(n)（0≤n≤N-1）是有限长序列，其Z变换为（4-28）为适应z可以沿Z平面更一般的路径取值，故沿Z平面上的一段螺线作等分角的采样，z的这些采样点zk为zk=AW-k

k=0,1,…,M-1（4-29）M为所要分析的复频率的点数，即采样点的总数，不一定等于N；A和W都是任意复数，可表示为:将式（4-30）与式（4-31）代入式（4-29），可得（4-30）（4-31）（4-32）因此有：图4-18螺线采样采样点在Z平面上所沿的周线如图4-18所示。由以上讨论和图4-18可以看出：（1）A0表示起始采样点z0的矢量半径长度，通常A0≤1;否则z0将处于单位圆|z|=1的外部。（2）θ0表示起始采样点z0的相角，它可以是正值或负值。（3）φ0表示两相邻采样点之间的角度差。φ0为正时，表示zk的路径是逆时针旋转的；φ0为负时，表示zk的路径是顺时针旋转的。（4）W0的大小表示螺线的伸展率。W0>1时，随着k的增加螺线内缩;W0<1时，则随k的增加螺线外伸;W0=1时，表示是半径为A0的一段圆弧。若又有A0=1,则这段圆弧是单位圆的一部分。当M=N，A=A0ejθ0=1，W=W0·e-jφ0=（W0=1，φ0=2π/N）这一特殊情况时，各zk就均匀等间隔地分布在单位圆上，这就是求序列的DFT。将式（4-29）的zk代入变换表达式（4-28），可得

0≤k≤M-1(3-33)直接计算这一公式，与直接计算DFT相似，总共算出M个采样点，需要NM次复数乘法与(N-1)M次复数加法。当N，M很大时，这个量很大，这就限制了运算速度。但是，下面我们将看到，通过一定的变换，以上运算可以转换为卷积形式，从而可以采用FFT算法，这样就可以大大提高运算速度。nk可以用以下表达式来替换将式（4-34）代入式（4-33），可得如果定义:n=0,1,…,N-1则它们的卷积为式中，k=0,1,…,M-1。式（4-38）正好是式（4-35）的一部分，因此式（4-35）又可以用卷积的形式表示为k=0,1,…,M-1由式（4-39）可看出，如果我们对信号按式（4-36）先进行一次加权处理，加权系数为;然后，通过一个单位脉冲响应为h(n)的线性系统即求g(n)与h(n)的线性卷积;最后，对该系统的前M点输出再做一次加权，这样就得到了全部M点螺线采样值X(zn)(n=0,1,…,M-1)。这个过程可以用图4-19表示。从图中我们看到，运算的主要部分是由线性系统来完成的。由于系统的单位脉冲响应 可以想象为频率随时间(n)呈线性增长的复指数序列。在雷达系统中,这种信号称为线性调频信号（ChirpSignal），因此，这里的变换称为线性调频Z变换。图4-1９Chirp-Z变换的线性系统表示k=0,1,…,M-1n=0,1,…,N-14.6.2Chirp-Z变换（CZT）的实现步骤由式（4-37）可看出，线性系统h(n)是非因果的，当n的取值为0到N-1，k的取值为0,1,…,M-1时，h(n)是在n=-(N-1)到n=M-1取值。也就是说，h(n)是一个有限长序列，点数为N+M-1，见图4-20（a）。输入信号g(n)也是有限长序列，点数为N。g(n)*h(n)的点数为2N+M-2，因而用圆周卷积代替线性卷积且不产生混叠失真条件是圆周卷积的点数应大于或等于2N+M-2。但是，由于我们只需要前M个值X（zk）(k=0,1,…,M-1)，对以后的其他值是否有混叠失真并不感兴趣，这样可将圆周卷积的点数缩减到最小为N+M-1。当然，为了进行基-2FFT运算，圆周卷积的点数应取为L≥N+M-1，同时又满足L=2m的最小L。这样可将h(n)先补零值点，补到点数等于L，也就是从n=M开始补L-(N+M-1)个零值点，补到n=L-N处，或补L-(N+M-1)个任意序列值，然后将此序列以L为周期而进行周期延拓，再取主值序列，从而得到进行圆周卷积的一个序列，如图4-20（b）所示。进行圆周卷积的另一个序列只需要将g(n)补上零值点，使之成为L点序列即可，如图4-20（c）所示。图4-20Chirp-Z变换的圆周卷积图（M≤n≤L-1时，h(n)和g(n)的圆周卷积不代表线性卷积）这样，我们可以列出CZT运算的实现步骤：（1）选择一个最小的整数L，使其满足L≥N+M-1，同时满足L=2m，以便采用基-2FFT算法。（2）将 （见图4-20（c））补上零值点，变为L点序列，因而有：

0≤n≤N-1N-1≤n≤L-1（4-40）并利用FFT法求此序列的L点DFT

0≤r≤L-1（4-41）（3）形成L点序列h(n)，如上所述，在n=0到M-1一段 ，n=M到L-N段取h(n)为任意值（一般为零），在L=N+M-1到L-1段取h(n)为 的周期延拓序列 ，即有:

0≤n≤M-1

L-N+1≤n≤L-1M≤n≤L-N（4-42）此h(n)可见图4-20（b）。实际上它就是图4-20（a）的序列 以L为周期的周期延拓序列的主值序列。对式（4-42）定义的h(n)序列，用FFT法求其L点DFT0≤r≤L-1（4-43）（4）将H(r)和G(r)相乘，得Q(r)=H(r)G(r)，Q(r)为L点频域离散序列。（5）用FFT法求Q(r)的L点IDFT，得h(n)和g(n)的圆周卷积L（4-44）式中，前M个值等于h(n)和g(n)的线性卷积结果［h(n)*g(n)］;n≥M的值没有意义，不必去求。g(n)*h(n)即g(n)与h(n)圆周卷积的前M个值可见图4-20（d）。（6）最后求X(zk)：（4-45）4.6.3运算量的估计

CZT的算法求X(zk)比直接求X(zk)的算法有效得多。CZT所需的乘法如下：（1）形成L点序列 ，但只有其中N点序列值，需要N次复乘，而系数 可事先准备好，不必在实时分析时计算。（2）形成L点序列h(n)。由于它是由 在-(N-1)≤n≤M-1的序列值构成的，而 是偶对称序列。如果设N>M，则只需要求得0≤n≤M-1一段N点序列值即可;h(n)也可事先准备好，不必实时分析时计算，因此，可不用考虑其计算量。同时，h(n)的L点FFT即H(r)也可预先计算好。（3）计算G(r)，q(n),需要二次L点FFT（或IFFT），共需要L

lbL次复乘。（4）计算Q(r)=G(r)H(r),需要L次复乘。（5）计算X(zk)=（0≤k≤M-1），需要M次复乘。综上所述，CZT总的复数乘法次数为mc=L

lbL+N+M+L

（4-46）

前面说过，直接计算式（4-33）的X(zk)需要NM次复数乘法。可以看出，当N,M都较大时（例如N,M都大于50时），CZT的FFT算法比直接算法的运算量要小得多。4.8FFT的其他应用4.8.1线性卷积的FFT算法——快速卷积

1.FFT的快速卷积法以FIR滤波器为例，因为它的输出等于有限长单位脉冲响应h(n)与有限长输入信号x(n)的离散线性卷积。

设x(n)为L点，h(n)为M点，输出y(n)为

y(n)也是有限长序列，其点数为L+M-1点。下面首先讨论直接计算线性卷积的运算量。由于每一个x(n)的输入值都必须和全部的h(n)值相乘一次，因而总共需要LM次乘法，这就是直接计算的乘法次数，以md表示为md=LM

用FFT算法也就是用圆周卷积来代替这一线性卷积时，为了不产生混叠，其必要条件是使x(n)，h(n)都补零值点，补到至少N=M+L-1，即:（4-65）0≤n≤L-1L≤n≤N-10≤n≤M-1M≤n≤N-1然后计算圆周卷积N这时，y(n)就能代表线性卷积的结果。用FFT计算y(n)的步骤如下：①求H(k)=DFT［h(n)］，N点；②求X(k)=DFT［x(n)］,N点；③计算Y(k)=X(k)H(k);④求y(n)=IDFT［Y(k)］，N点。步骤①、②、④都可用FFT来完成。此时的工作量如下：三次FFT运算共需要 次相乘，还有步骤③的N次相乘，因此共需要相乘次数为（4-66）比较直接计算线性卷积(简称直接法)和FFT法计算线性卷积(简称FFT法)这两种方法的乘法次数。设式（4-65）与式（4-66）的比值为Km，则（4-67）分两种情况讨论如下：（1）x(n)与h(n)点数差不多。例如，M=L，则N=2M-1≈2M，则这样可得下表：M=L8163264128256512102420484096Km0.5720.9411.62.785.928.821629.2454.999.9当M=8时，FFT法的运算量大于直接法；当M=32时，二者相当；当M=512时，FFT法运算速度可快16倍；当M=4096时，FFT法约快100倍。可以看出，当M=L且M超过32以后，M越长，FFT法的好处越明显。因而将圆周卷积称为快速卷积。（2）当x(n)的点数很多时，即当L>>M。通常不允许等x(n)全部采集齐后再进行卷积;否则，使输出相对于输入有较长的延时。此外，若N=L+M-1太大，h(n)必须补很多个零值点，很不经济，且FFT的计算时间也要很长。这时FFT法的优点就表现不出来了，因此需要采用分段卷积或称分段过滤的办法。即将x(n)分成点数和h(n)相仿的段，分别求出每段的卷积结果，然后用一定方式把它们合在一起，便得到总的输出，其中每一段的卷积均采用FFT方法处理。有两种分段卷积的办法:重叠相加法和重叠保留法。在这一公式中，已经把每一段的时间原点放在该段的起始点，而不是x(n)的原点。这种分段方法示于图4-28中，每段xi(n)和h(n)的圆周卷积结果以yi(n)表示，如图4-28（b）所示，图中已标出每一段输出段开始的(M-1)个点，0≤n≤M-2部分舍掉不用。把相邻各输出段留下的序列衔接起来，就构成了最后的正确输出，即（4-73）式中:M-1≤n≤N-1其他n

（4-74）这时，每段输出的时间原点放在yi

′(n)的起始点，而不是y(n)的原点。线性相关的FFT算法利用FFT来计算相关函数,也就是利用圆周相关来代替线性相关,常称之为快速相关.这与利用FFT的快速卷积相类似(它是利用圆周卷积代替线性卷积),也要利用补零值点的办法来避免混叠失真.方法设x(n)的长度为L点,y(n)的长度为M点,要求线性相关

利用FFT法来求线性相关,是用圆周相关来代替线性相关,选择,令4.8.2信号消噪

假设信号在传输过程中，受到噪声的干扰，则在接收端得到的信号由于受到噪声的干扰，信号将难以辨识。用FFT方法消噪就是对含噪信号的频谱进行处理，将噪声所在频段的X(k)值全部置零后，再对处理后的X(k)进行离散傅里叶反变换（IFFT）可得原信号的近似结果。这种方法要求噪声与信号的频谱不在同一频段，否则，将很难处理。

例4-5

将上述消噪原理用于语音消噪。图4-30是一个实际例子，语音信号受到了强烈的啸叫噪声干扰，无法听清语音意思，如图4-30(a)，信号淹没在噪声中（信噪比只有-10dB）。试用FFT方法消噪。先作FFT分析，得到其功率谱如图4-30(b)，可见在频率2.5kHz附近有一极强分量。这就是啸叫噪声干扰。图4-30(b)中频率在30～800Hz范围是语音信号。对频谱进行修正，去除噪声频段，即将大于2.5kHz部分的X(k)值全部置为零，图4-30(c)是去噪后的功率谱。再由反变换（IFFT）重构信号得到原语音信号如图4-30(d)。这时信噪比为14dB，提高了24dB。这就是早期的数字式录音音乐中所采用的消噪方法。图4-30语音信号消噪过程信号淹没在啸叫噪声中；(b)信号与噪声的功率谱；(c)去噪后的功率谱；(d)重构原语音信号4.8.3FFT在双音多频（DTMF）信号中的应用双音多频信号DTMF是按键式电话信令中的一般名称，它等效于在贝尔系统内部正在使用的按钮式拨号系统。近几年DT