《数字信号处理》课件-第4章

4.1引言

4.2提高DFT运算效率的基本途径4.3基2时分FFT算法

4.4基2频分FFT算法

4.5IDFT的快速算法

4.6实序列DFT的有效计算方法

4.7线性调频Z(Chirp-Z)变换算法习题与上机题

离散傅里叶变换(DFT)是信号分析与处理中的一种重要变换。因为直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大，所以在快速傅里叶变换(FastFourierTransform，FFT)出现以前，直接用DFT进行谱分析和信号的实时处理是不切实际的。4.1引言自从1965年图基(J.W.Tukey)和库利(J.W.Cooley)在《计算数学》(MathematicsComputation)杂志上发表了著名的《机器计算傅里叶级数的一种算法》论文后，桑德-图基等快速算法相继出现，很快形成了一套高效运算方法，这种算法使DFT的运算效率提高了几个数量级，为数字信号处理技术应用于各种信号的实时处理创造了良好的条件，大大推动了数字信号处理技术的发展。多年来，人们继续寻求更快、更灵活的好算法，目前除了基2、基4算法外，还有分裂基、混合基等算法。对于N=1024，直接计算需要复数乘法1048576次，而几种常用FFT算法使计算量大大降低，例如基2算法需要复数乘法5120次，基4算法需要复数乘法3840次，

分裂基算法需要复数乘法3413次。本章重点介绍基2快速傅里叶变换算法。设x(n)为N点有限长序列，其DFT为

(4.2.1)

其反傅里叶变换(IDFT)为

(4.2.2)4.2提高DFT运算效率的基本途径正变换与反变换的差别只在于的指数符号不同，以及差一个乘数因子，因此我们以下只讨论DFT的运算量，IDFT的运算量与DFT的基本相同。

一般来说，X(k)为复数(只有x(n)为实数且条件偶对称时，X(k)才为实数)，x(n)多数情况下也为复数(例如：信号经A/D采样，数字下变频芯片之后输出为I、Q正交两路，详见第9章)。因此，计算一点X(k)需要N次复数乘法和N-1次复数加法。由于有N点X(k)，因此完成N点序列的DFT共需要N2次复数乘法，即

CM=N2

(4.2.3)

以及N(N-1)次复数加法，即

CA=N(N-1)≈N2

(4.2.4)实际实现中复数乘法和加法都要依赖实数乘法和加法。由于一次复数乘法需要4次实数乘法和两次实数加法，一次复数加法需要2次实数加法，因此直接计算DFT的实数乘法次数为

RM=4N2

(4.2.5)

实数加法次数为

RA=2N(N-1)+2N2≈4N2

(4.2.6)所以，直接计算DFT时的运算量与N2成正比。设N=1024，则直接计算DFT需要100多万次(1048576次)复数乘法。在实时计算较宽带宽的频谱时，即使使用目前运算速度最快的DSP芯片也无法满足计算要求。

观察式(4.2.1)，要想提高DFT的运算效率，只有利用DFT定义式中旋转因子的特点。根据第3章所讲内容，我们可知有以下几个特点:

(1)共轭对称性：。

(2)周期性：。

(3)可约性：。利用共轭对称性，可以合并式(4.2.1)中的乘法项，从而减少了近一半的乘法次数。

利用旋转因子的周期性和可约性，可以将长序列的DFT分解成短序列的DFT。由于DFT的运算量与N2成正比，N越小，计算量越小。因此将长序列DFT变为短序列DFT是减少DFT运算量的根本途径。快速傅里叶变换算法正是基于此思路而发展起来的。

将长序列的DFT划分为短序列的DFT有两种划分方法：一种是将时域序列按奇偶划分，称之为时间抽取(DecimationinTime，DIT)算法或基2时分算法;一种是将频域DFT按奇偶划分，称之为频率抽取(DecimationinFrequency，DIF)算法或基2频分算法。4.3.1基2时分蝶式运算定理

设序列x(n)的长度为N，且满足N=2M，M为整数，其N点DFT为X(k)=DFT［x(n)］，0≤k，n≤N-1。按n的奇偶将x(n)分解为两个点的子序列

x1(r)=x(2r)，r=0，1，…，(4.3.1)

x2(r)=x(2r+1)，r=0，1，…，(4.3.2)4.3基2时分FFT算法若X1(k)=DFT［x1(r)］，X2(k)=DFT［x2(r)］，0≤k≤

-1，则有

(4.3.3a)

(4.3.3b)

证明x(n)的N点DFT为

利用旋转因子的可约性，即，

X(k)可以写成

由于X1(k)和X2(k)都隐含周期性，周期为，因此上式可以写为

当0≤k≤

-1时

当≤k≤N-1时

令，则有

由于，因此

最后

由证明过程可得，N点DFT可以分解为两个点的DFT，它们按式(4.3.3)又组合成一个N点的DFT。式(4.3.3)也称为蝶式运算公式。由蝶式运算定理可以看出，要完成一个蝶形运算(即计算一次蝶式运算定理)，需要一次复数乘法和两次复数加法，计算N点DFT共需计算两个点的DFT和个蝶形运算。计算点DFT需要次复数乘法和次复数加法。当N1时，计算N点DFT共需由此可见，仅经过一次分解，就使运算量减少近一半。如果N=2M，则我们可以继续分解，直至分解到1点DFT，运算量即可大大减少。4.3.2基2时分的蝶形流图与计算量分析

1.蝶形流图

基2时分FFT蝶式运算定理的运算公式可以用图4.3.1(a)所示的信号流图表示，由于这个图形呈蝶形，故称为蝶形流图，图(b)是其简化形式。一个蝶形流图是基2时分算法的一个基本单元。图4.3.1DIT-FFT蝶式运算信号流图(a)蝶形流图;(b)蝶形流图简化形式图中左面两支为输入，中间以一个点表示加、减运算，右上支为相加输出，右下支为相减输出。如果在某一支路上信号需要进行乘法运算，则在该支路上标以箭头，并将

相乘的系数标在箭头边。这样，蝶式运算定理的运算公式(4.3.3)所表示的运算，可用图4.3.1(b)所表示的“蝶形结”来表示。采用这种表示法，可将基2时分的分解过程用计算流图来表示。

【例4.3.1】将8点序列的DFT用基2时分FFT算法进行分解。

解

(1)第一次分解。

将8点序列x(n)分解为两个4点序列x1(n)和x2(n)，x1(n)为偶数序列，x2(n)为奇数序列，即

x1(n)={x(0)，x(2)，x(4)，x(6)}

x2(n)={x(1)，x(3)，x(5)，x(7)}

4点序列x1(n)和x2(n)分别做4点DFT得到X1(k)和X2(k)，由X1(k)和X2(k)通过蝶形构造获得8点DFTX(k)。第一次分解的旋转因子为(k=0，1，2，3)，运算流图如图4.3.2所示。

(2)第二次分解。将两个4点序列x1(n)和x2(n)分解为四个2点序列x3(n)、x4(n)、和。

图4.3.2

N点DFT基2时分第一次分解运算流图(N=8)

2点序列分别做2点DFT得到X3(k)、X4(k)、、

，由四个2点DFT通过蝶形构造获得两个4点DFT，再通过蝶形构造获得8点DFTX(k)。第二次分解的旋转因

子为，运算流图如图4.3.3所示。图4.3.3

N点DFT基2时分第二次分解运算流图(N=8)

(3)第三次分解。

将四个2点序列分解为8个1点序列，其中

由于1点序列的DFT值为其序列本身，因此在最后一次分解后，流图中已经没有直接计算DFT的环节。第三次分解旋转因子为，运算流图如图4.3.4所示。

由以上例子我们可以看出，由于每一次分解都是按输入序列在时域上的次序是偶数还是奇数来抽取的，最终分解成N个1点DFT，因此称为基2时分FFT。图4.3.4

N点DFT基2时分FFT运算流图(N=8)

2.计算量分析

对于任何一个2的整数幂N=2M，总是可以通过M次分解直到1点的DFT运算。这样的M次分解，就构成了从x(n)到X(k)的M级运算过程。从图4.3.4可以看出，每一级运算都由个蝶形运算构成。由于每一个蝶形运算需要一次复数乘法和两次复数加法，因此每一级分解所需的复数乘法和加法次数分别为

N点DFT(M级分解)所需总的复数乘法和加法次数分别为

(4.3.4)

(4.3.5)我们知道直接计算N点DFT需要的复数乘法为N2次，复数加法为N(N-1)次，当N1时，基2时分FFT算法将比直接计算DFT的运算量大大减少。基2时分算法与直接计算N点DFT的复数乘法次数比为

(4.3.6)复数加法次数比为

(4.3.7)若N=210，则直接计算需要的复数乘法和加法次数分别为CM=10242=1048576，CA=1024(1024-1)=1047552，基2时分需要的复数乘法和加法次数分别为CM=5120，CA=10240,复数乘法次数比为αM≈，即基2时分FFT的复数乘法运算效率提高了约200倍。图4.3.5显示了直接计算DFT与基2时分FFT算法所需乘法次数的比较曲线，显然，N越大，FFT算法的效率越高。图4.3.5直接计算DFT与基2时分FFT算法所需乘法次数的比较曲线

4.2.4基2时分FFT算法的运算规律及编程思想由基2时分FFT运算流图，可以看出这种算法的运算规律。

1.原位计算

由图4.3.4的基2时分FFT运算流图可以看出，同一级运算中，每个蝶形的两个输入数据只对计算本蝶形有用，而且每个蝶形的输入、输出数据结点又同在一条水平线上，这就意味着计算完一个蝶形后，所得输出数据可立即存入原输入数据所占用的存储单元。这样，经过M级运算后，原来存放输入序列数据的N个存储单元中便依次存放X(k)的N个值。

这种利用同一存储单元存储蝶形计算输入、输出数据的方法称为原位运算。

2.旋转因子的变化规律

在基2时分FFT运算流图中，旋转因子的变化和两个节点之间的距离都呈现一定的规律，p称为旋转因子的指数。设L表示从左到右的运算级数，即L=1，2，…，M。设N=23=8，L=1，2，3，第L级共有2L-1个不同的旋转因子，即

L=1时，;

L=2时，;

L=3时，。

对N=2M的一般情况，L=1，2，…，M，第L级的旋转因子为

(4.3.8)在实际编程实现中，旋转因子有三种计算方法。

(1)直接计算。在用高级编程语言实现FFT，且对运算时间没有要求的理论分析场合下，

一般采用直接计算求旋转因子。根据欧拉公式有

如果采用相同的旋转因子计算完后存储起来的办法，则共需计算个，即N个正弦函数。

(2)递推计算。根据旋转因子的特点，即

可以利用递推的办法求解旋转因子。根据上式，在每一级只需计算一个旋转因子，其余的旋转因子均可通过上式求出。因此N=2M点FFT共需计算2M个正弦函数。

(3)查表。在要求实时计算FFT的应用场合，用DSP汇编或FPGA编程实现FFT算法时，一般用查表方法。由于N点基2时分FFT算法中共有个不同的旋转因子，因此事先计算好个旋转因子表以供查询，可以省去计算正弦函数的时间。

3.输入序列的逆序

基2时分FFT算法输入序列的排序看起来似乎很乱，但仔细分析就会发现这种排序是

有规律的，我们称之为比特逆序。

对于N=2M，N个顺序数可以用M位二进制数(nM-1nM-2

…n1n0)2来表示，即(n)10=(nM-1nM-2…n0)2，定义

ρM(n)=(n0n1…nM-1)2

(4.3.9)

为n的M位比特逆序数。

【例4.3.2】计算n=2，12的二进制比特逆序数ρ4(n)。

解

造成基2时分FFT算法的输入序列顺序逆序的原因是将输入序列x(n)按下标n的奇偶不断分解造成的。表4.3.1列出了N=8的自然顺序数和比特逆序数。由表显而易见，只要

将自然顺序的二进制数按位倒置，就能得到比特逆序数。

在实际应用中，一般先按自然顺序数将输入序列存入存储单元A(m)，m=0，1，…，N-1。

为了得到比特逆序的排列，可以通过变址运算来完成，变址过程如图4.3.6所示。表4.3.1N=8的自然顺序数和比特逆序数图4.3.6比特逆序的变址过程(N=8)由表4.3.1可以看出，自然顺序的次序增加是在二进制数低位加1，二进制进位由低向高;比特逆序的次序增加在二进制数高位加1，进位由高向低。因此，实际求数

n(n=0，1，2，…，N-1)的比特逆序ρM(n)步骤如下:

(1)0的比特逆序为0，因此ρM(0)=0。

(2)求n+1的逆序ρM(n+1)。如ρM(n)的二进制数最高位为0，则将该位赋1，循环结束。

如ρM(n)的二进制数从高位开始数前几位均是1，则把为1的最低位的右边为0的位赋1，同时将为1的位赋0，循环结束。

综合以上求比特逆序步骤和图4.3.6，将自然顺序的输入序列通过比特逆序变址变为比特逆序输入次序的程序流图如图4.3.7所示。图4.3.7比特逆序程序流图图4.3.7所示的程序流图是求j的逆序，同时将存储单元内自然顺序排列的输入序列变为逆序排序，i存储的是j-1的逆序数。由于j=0，N-1的逆序既是数据本身，因此程序循

环从j=1到j=N-2。s存放的是二进制数每一位的权，每次j循环开始，s为二进制最高位的权。程序流图中虚线框部分是将为1的位赋0的过程。

4.编程思路与程序流图

观察基2时分FFT蝶形流图，可归纳出一些对编程有用的运算规律：第L级运算中，每个蝶形的两个输入数据相距2L-1个点，同一个旋转因子对应着间隔为2L点的2M-L个蝶形。总结上述运算规律，可以采用以下方法。首先从第一级(输入端)开始，逐级进行，共进行M级运算。在进行第L级运算时，依次求出2L-1个不同的旋转因子，每求出一个旋转因子，就计算完它对应的所有2M-L个运算蝶。因此基2时分FFT程序由逆序和三重循环组成。三重循环与流图的对应关系为：最外层循环(循环变量为L)对应蝶形流图的级数M，中间循环(循环变量为k)对应流图中每一级不同旋转因子对应的运算蝶数目2L-1，最内层循环(循环变量为J)对应每一级相同旋转因子运算蝶的顺序，程序流图如图4.3.8所示，C语言程序见附录。图4.3.8基2时分FFT算法程序流图

【例4.3.3】设x(n)为N点有限序列，N为偶数，其N点DFT为X(k)，试用X(k)表示如下序列的N点DFT。

解

(1)解法1。

(2)解法2。用基2时分蝶式运算定理将x2(n)分解为两个

点的序列x20(n)和x21(n)。

根据蝶式运算定理有

因此有

同理，将x(n)分解为两个点的序列x0(n)和x1(n)。

因为

所以

令，有

因为

所以

4.4.1基2频分蝶式运算定理

设序列x(n)的长度为N=2M，其N点DFT为X(k)=DFT［x(n)］，0≤k≤N-1，按k的奇偶将X(k)分解为两个点的子序列X1(k)和X2(k)。

4.4基2频分FFT算法若x1(n)=IDFT［X1(r)］，x2(n)=IDFT［X2(r)］，0≤n≤

-1，则

(4.4.1)

(4.4.2)证明由证明过程可得，N点序列x(n)按式(4.4.1)和式(4.4.2)分解为两个点的子序列，这两个点子序列的DFT正好是N点DFTX(k)的偶数组和奇数组。由蝶式运算定理可以看出，要完成一个蝶形运算，需要一次复数乘法和两次复数加法，计算N点DFT共需计算两个点DFT和次蝶形运算。计算点DFT需要次复数乘法和次复数加法。当N>>1时，计算N点DFT共需

次复数加法。由此可见，仅经过一次分解，就使运算量减少近一半。如果N=2M，则我们可以继续分解，直至分解

到1点DFT，运算量即可大大减少。4.4.2基2频分的蝶形流图与计算量分析

1.蝶形流图

基2频分FFT蝶式运算定理的运算公式可以用图4.4.1所示的信号流图表示，一个蝶形流图是基2频分算法的一个基本单元。图4.4.1DIF-FFT蝶式运算信号流图

【例4.4.1】将8点序列的DFT用基2频分FFT算法进行分解。

解

(1)第一次分解。

将8点序列x(n)通过蝶形构造分解为两个4点序列x1(n)和x2(n),分别做4点DFT得到X1(k)和X2(k)，即

X1(k)={X(0)，X(2)，X(4)，X(6)}

X2(k)={X(1)，X(3)，X(5)，X(7)}

由X1(k)和X2(k)获得8点DFTX(k)，旋转因子为(n=0，1，2，3)，运算流图如图4.4.2所示。图4.4.2

N点DFT基2频分第一次分解运算流图(N=8)

(2)第二次分解。

将两个4点序列分别通过蝶形构造分解为四个2点序列，分别做2点DFT得到

由2点DFT即可获得8点DFTX(k)，第二次分解的旋转因子为(n=0，1)，运算流图如图4.4.3所示。图4.4.3

N点DFT基2频分第二次分解运算流图(N=8)

(3)第三次分解。

将四个2点序列通过蝶形构造分解为8个1点序列，其中

由于1点序列的DFT值为其序列本身，因此，在最后一次分解后，流图中已经没有直接计算DFT的环节。第三次分解旋转因子为，运算流图如图4.4.4所示。图4.4.4

N点DFT基2频分FFT运算流图(N=8)由以上例子我们可以看出，由于每一次分解都是按输入序列在频域上的次序是偶数还是奇数来抽取的，最终分解成N个1点DFT，因此称为基2频分FFT。

2.计算量分析

对于任何一个长度为2的整数幂(N=2M)序列，其DFT总是可以通过M次分解到1点的DFT运算。这样的M次分解，就构成了从x(n)到X(k)的M级运算过程。从图4.4.4可以看出，每一级运算都由个蝶形运算构成。由于每一个蝶形运算需要一次复数乘法和两次复数加法，因此N点DFT(M级分解)所需总的复数乘法和加法次数分别为

(4.4.3)

CA=NM=Nlog2N

(4.4.4)由于基2时分FFT算法和基2频分FFT算法抽取的域不同，因此它们的蝶形构造发生的地方、基本运算和蝶形流图都不相同。基2时分是时域抽取，频域蝶形构造，输入逆序，输出顺序;基2频分是频域抽取，时域蝶形构造，输入顺序，输出逆序。但两者的运算量是相同的。具体的编程思想不再赘述。

在前面章节的学习中我们或者使用基2时分FFT算法，或者使用基2频分FFT算法将某一个序列的DFT计算出来。下面我们举个例子，在同一个序列中某些步骤用基2时分，某些步骤用基2频分。

【例4.4.2】在N=8的DFT中，分两步做快速算法，第一步用基2频分将8点DFT化为两个4点DFT，第二步用基2时分将这两个4点DFT化为四个2点DFT，并画出流图。

解根据题目要求，第一步使用基2频分，构造4点序列x1(n)和x2(n)，其DFT分别是

X1(k)={X(0)，X(2)，X(4)，X(6)}，

X2(k)={X(1)，X(3)，X(5)，X(7)}

运算流图如图4.4.5所示。图4.4.5例4.4.2第一步基2频分FFT第一次分解运算流图第二步使用基2时分，以x1(n)为例，将其按照n的奇偶分成2点序列x3(n)和x4(n),由2点DFTX3(k)和X4(k)通过蝶形运算即可构造出X1(k)，如图4.4.6所示。同理可得X2(k)。

最后，将前两步的结果合并起来，运算流图如图4.4.7所示。图4.4.6例4.4.2第二步基2时分FFT分解运算流图图4.4.7例4.4.2运算流图长度为N的有限长序列x(n)的离散傅里叶变换为

(4.5.1)

其离散傅里叶逆变换为

(4.5.2)4.5IDFT的快速算法比较以上两式，可见离散傅里叶变换对具有对称性。因此可以利用DFT的快速算法(FFT)来实现IDFT。具体来说，就是将FFT运算流图中的旋转因子取共轭(虚部取反)，最后输出乘以。此时，输入为X(k)，输出即为x(n)。前面讨论的序列x(n)的DFT快速算法都是复数运算，包括序列x(n)也认为是复数。在实际应用中，信号有可能是实序列，任何实数都可看成虚部为零的复数。例如，求某实信号x(n)的频谱，可认为是将实信号加上数值为零的虚部变成复信号(x(n)+j0)，再用FFT求其离散傅里叶变换。这种做法很不经济，因为把实序列变成复序列，存储器要增加一倍，且计算机运行时，即使虚部为零，也要进行涉及虚部的运算，浪费了运算量。

合理的解决方法是利用复数FFT对实数进行有效计算。4.6实序列DFT的有效计算方法

1.一个N点FFT同时计算两个N点实序列的DFT

【例4.6.1】设x1(n)和x2(n)是彼此独立的两个N点实序列，且X1(k)=DFT［x1(n)］，X2(k)=DFT［x2(n)］。

试通过一次FFT运算同时求得X1(k)和X2(k)。

解首先将x1(n)和x2(n)分别作为一个N点复序列x(n)的实部及虚部，即

x(n)=x1(n)+jx2(n)通过调用一个N点的FFT算法即可获得x(n)的DFT序列X(k)。

利用离散傅里叶变换的共轭对称性可获得X1(k)和X2(k)，即

下面分析计算量。通过调用N点FFT算法求X(k)需要的复数乘法次数和复数加法次数分别为

根据X(k)求X1(k)和X2(k)需要的复数乘法次数和复数加法次数分别为

CM2=N-1，CA2=2(N-1)

因此，该方法总共需要的复数乘法次数和复数加法次数分别为

如果调用两次N点FFT算法直接计算X1(k)和X2(k)，所需的运算量为

CM=Nlog2N，CA=2Nlog2N

设N=1024，则直接调用两次FFT计算需要的CM=10240，CA=20480。利用上述算法则有CM=6143，CA=12286，计算量减少近半。

2.一个N点FFT计算一个2N点实序列的DFT

【例4.6.2】设x(n)是2N点的实序列，试设计用一次N点的FFT算法完成计算X(k)的高效算法。

解将x(n)分解为偶数序列x1(n)和奇数序列x2(n)，即

x1(n)=x(2n)，n=0，1，…，N-1

x2(n)=x(2n+1)，n=0，1，…，N-1将x1(n)和x2(n)分别作为一个N点复序列y(n)的实部及虚部，即

y(n)=x1(n)+jx2(n)

通过N点FFT运算可得到y(n)的N点DFTY(k)，根据前面的讨论可得

根据基2时分蝶式运算定理得

下面分析计算量。通过调用N点FFT算法求Y(k)需要的复数乘法次数和复数加法次数分别为

根据Y(k)求X1(k)和X2(k)需要的复数乘法次数和复数加法次数分别为

CM2=N-1，

CA2=2(N-1)

利用基2时分蝶式运算定理求X(k)需要的复数乘法次数和复数加法次数分别为

CM3=N，

CA3=2N因此，该方法总共需要的复数乘法次数和复数加法次数分别为

直接计算2N点FFT计算量为

CM=Nlog22N，

CA=2Nlog22N

设N=1024，则直接计算CM=11264，CA=22528。

利用上述算法则有CM=7167，CA=14334。前面已经介绍过，采用FFT算法可以很快计算出全部N点DFT值，即Z变换X(z)在z平面单位圆上的全部等间隔取样值。而实际中也许不需要计算整个单位圆上Z变换的取样值。例如对于窄带信号，只需要对信号所在的一段频带进行分析，这时希望频谱的取样集中在这一频带内，以获得较高的分辨率，而频带以外的部分可不考虑。另外，实际中有时也对非

单位圆上的取样值感兴趣。4.7线性调频Z(Chirp-Z)变换算法例如需要知道Z变换的极点所在处的频率，如极点位置离单位圆较远，则其单位圆上的频谱就很平滑，这时很难从中识别出极点所在处的频率，此时，如果取样不是沿单位圆而是沿一条接近这些极点的弧线进行，则极点所在频率上的频谱将出现明显的尖峰，由此可较准确地测定极点频率。

Z变换采用螺旋线取样是一种适合于以上需要的变换，且可以采用FFT来快速计算，这种变换也称做线性调频Z变换(也称Chirp-Z变换或CZT)，它是适用于这种更为一般情况下由x(n)求X(zk)的快速算法。4.7.1Chirp-Z变换的基本原理

已知x(n)(0≤n≤N-1)为有限长序列，其Z变换为

(4.7.1)

为适应z可以沿z平面上的一段螺旋线做等分角的取样，令z的取样点为

zk=AW-k，k=0，1，…，M-1

(4.7.2)其中，M为所要分析的复频谱的点数(不一定等于N);A和W是任意复数，可表示为

(4.7.3)

(4.7.4)图4.7.1CZT在z平面的螺旋线取样将式(4.7.3)和式(4.7.4)带入式(4.7.2)中，可得

(4.7.5)

有

取样点在z平面上所沿的周线如图4.7.1所示。由以上讨论及图可见:

(1)A0表示起始取样点的矢量半径长度，通常A0≤1，否则z0将处于单位圆外。

(2)ω0表示起始取样点z0的相角，可以是正值也可以是负值。

(3)φ0表示两相邻取样点之间的等分角度，φ0为正，表示zk的路径是逆时针旋转的;φ0为负，表示zk的路径是顺时针旋转的。

(4)W0的大小表示螺旋线的伸展率，当W0<1时，随着k的增加螺旋线外伸;当W0>1时,随着k的增加螺旋线内缩；W0=1表示半径为A0的一段圆弧，若A0=1，这段圆弧则是单位圆的一部分。当

时，各个zk均匀等间隔地分布在单位圆上，此时CZT即为DFT。

将式(4.7.2)中的zk代入Z变换表达式(4.7.1)中，可得

(4.7.6)

显然，同直接计算DFT情况相仿，按照式(4.7.6)计算出全部M点取样值需要NM次复数乘法和(N-1)M次复数加法，当N及M较大时，计算量迅速增加。如果通过一定的变换，如采用布鲁斯坦(Bluestein)等式，将以上运算转换为卷积形式，则可采用FFT进行，这样即可大大提高计算速度。布鲁斯坦等式为

(4.7.7)将上式代入式(4.7.6)，可得

令

(4.7.8)

(4.7.9)

则有

(4.7.10)式(4.7.10)说明，如对信号x(n)先进行一次加权处理，加权系数为，然后通过一个单位脉冲响应为h(n)的线性时不变系统，最后对该系统的前M点输出再作一次的加权，就可得到全部M点螺旋线取样值。运算流程如图4.7.2所示。图4.7.2线性调频Z变换运算流程由于系统的单位脉冲响应与频率随时间成线性增加的线性调频信号(即ChirpSignal)相似，因此这种算法称为Chirp-Z变换。

4.7.2Chirp-Z变换的实现步骤

由式(4.7.10)可以看出，由于输入信号g(n)是长度为N的有限长序列，序列是无限长的，而计算0~M-1点卷积g(k)*h(k)所需要的h(n)是取值在n=-(N-1)~(M-1)那一部分的值，如图4.7.3(a)所示。因此，可认为h(n)是一个有限长序列，长度为L=N+M-1。所以，Chirp-Z变换为两个有限长序列的线性卷积

g(k)*h(k),可用循环卷积通过FFT来实现，即

g(k)*h(k)=IDFT［G(r)H(r)］根据第3章所学知识，g(k)*h(k)的长度为2N+M-2，因而用循环卷积求线性卷积时，循环卷积的长度最短为2N+M-2。由于我们只需要求0~M-1处的线性卷积结果，因此选循环卷积的长度为两个序列中最长的长度，即L=N+M-1点。另外，h(n)的主值序列(0≤n≤L-1)可由h(n)作周期延拓后取0≤n≤L-1部分值获得，如图4.7.3(b)所示。

将h(n)与g(n)作循环卷积后，其输出的前M个值就是CZT的M个值。这个循环卷积的过程可在频域上通过FFT实现。图4.7.3

h(n)的主值序列求解过程总之，CZT运算步骤如下：

(1)选择一个最小整数L，使其满足L≥N+M-1，同时满足L=2m，以便采用基2FFT算法。

(2)对x(n)加权并补零得到g(n)：

(4.7.11)并利用FFT算法求此序列的L点DFTG(r)。

(4.7.12)

(3)求h(n)的主值序列:

(4.7.13)并利用FFT算法求此序列的L点DFTHm(r):

(4.7.14)

(4)将Hm(r)和G(r)相乘，得到Q(r)=G(r)Hm(r)。

(5)用FFT算法求L点Q(r)的IDFT，得到hm(n)与g(n)的循环卷积。

(4.7.15)其中，q(n)的前M-1点等于h(n)与g(n)的线性卷积结果(与第3章所讲的线性卷积与循环卷积的关系不同，这里的hm(n)是h(n)以L为周期延拓后的主值序列)。

(6)求X(zk):

(4.7.16)4.7.3Chirp-Z变换运算量的估算

采用上小节的算法后，CZT所需的运算量如下:

(1)形成L点序列，由于补零关系，

g(n)只有N点有序列值，需要N次复数乘法。系数可以通过递推求得。令

则

初始条件。

因此通过递推求Cn需要2N次复数乘法。

(2)形成L点序列hm(n)，由于它是由在

-N+1≤n≤M-1段内的序列值构成的，是偶对称序列，如果设N≥M，则只需求0≤n≤N-1段内的N点序列值。求

也可以用递推方式求解，因此需要复数乘法2N次。

(3)计算Hm(r)、G(r)和q(n)共需要3次L点FFT算法，因此需要复数乘法次数为。

(4)计算Q(r)=G(r)Hm(r)需要L次复数乘法。

(5)计算X(zk)=，0≤k≤M-1需要M次复数乘法。综上所述，Chirp-Z变换总共需要的复数乘法次数为

(4.7.17)

前面讨论过直接计算X(zk)需要复数乘法为NM次，当N=925，M=100时，直接计算需要CM=92500次，采用快速算法后所需的复数乘法次数为CM=21109，计算量减少较多。总之，与DFT相比，CZT有以下特点：

(1)输入序列的长度N与输出序列的长度M不需要相等；(2)N及M不必是合成数，二者均可为素数；

(3)zk点的角度间隔φ0是任意的，因此频率分辨率也是任意的；

(4)周线不必是z平面上的圆；

(5)起始点z0可任意选定，因此可以从任意频率上开始对输入数据进行窄带高分辨率分析；

(6)若A=1，M=N，即使N为素数，也可用Chirp-Z变换计算DFT。4.7.4用Matlab计算Chirp-Z变换

Matlab信号处理工具箱提供了内部函数czt用于实现Chirp-Z变换，调用格式如下：

y=czt(x，M，W，A);

y=czt(x);

y=czt(x，M，W，A)用于计算指定参数M、W、A下的Chirp-Z变换，而y=czt(x)则使用默认参数进行Chirp-Z变换，默认参数为M=length(x)，W=exp(j*2*pi/m)，A=1，此时即计算DFT。

【例4.7.1】某一低通滤波器的零、极点分布如图4.7.4所示，利用Chirp-Z变换观察滤波器零点特性。

clear；fs=2000; ％采样频率为2000Hz

fp=150；

p1=0.9*exp(-j*2*pi*fp/fs);

p2=p1′;p=［p1，p2］′; ％极点在150Hz处

fz1=500;

fz2=300;

z1=1.2*exp(-j*2*pi*fz1/fs)；z2=z1′;

z3=1.2*exp(-j*2*pi*fz2/fs);z4=z3′;图4.7.4滤波器零、极点分布图

z=［z1，z2z3z4］′;

％零点在300Hz和500Hz处的单位圆外，半径为1.2

［b，a］=zp2tf(z，p，1);

ww=0：0.005*pi：2*pi;

h1=freqz(b，a，ww);

hn1=ifft(h1);hn=real(hn1)；

N=length(ww)； ％滤波器单位脉冲响应长度(401点)

f2=700;f1=50；M=201；

W=exp(-j*2*pi*(f2-f1)/(M*fs));

％在z平面ω=2πf1~2πf2之间取M点

A=1.19*exp(j*2*pi*f1/fs); ％取样圆弧半径为1.19

yzz=czt(hn，M，W，A); ％求CZT由于零点在单位圆外，半径为1.2处，因此CZT的圆弧线半径为1.19，起始点ω=2π×50,终止点为ω=2π×700，取样点数为201点，如图4.7.5所示。滤波器单位圆上的幅频响应如图4.7.6所示，由于单位圆距离零点较远，两个零点处的频率在幅频响应里不明显。由于CZT取样曲线接近零点，因此在300Hz和500Hz处的频率在CZT中非常明显，CZT如图4.7.7所示。由于CZT的取样圆弧离极点位置较远，因此极点处的频率在CZT曲线中已不明显。图4.7.5CZT取样螺旋线图图4.7.6滤波器幅频响应图图4.7.7滤波器CZT图用Matlab按快速算法的实现步骤编写的CZT程序如下(参数设置与上相同):

L=M+N-1;

gn=zeros(1，L);

gn(1：N)=A.^(-(0：N-1)).*W.^((0：N-1).^2./2).*hn;

gf=fft(gn，L);

hhn(1：M)=W.^(-(0：M-1).^2./2);

hhn(M+1：L)=W.^(-(L-(M：L-1)).^2./2);

hhnf=fft(hhn);

Qr=hhnf.*gf;

qn=ifft(Qr);

yzz=qn(1：M).*W.^((0：M-1).^2./2)

【例4.7.2】设信号为x(n)，它由4个频率分别为60Hz、64Hz、95Hz和100Hz的正弦序列组合而成。采样频率为600Hz，样点数为200点。试用CZT观察信号频谱。

解为比较起见，分别用200点的DFT和CZT观察信号频谱。

图4.7.8所示为x(n)的DFT，由于点数太少，60Hz与64Hz、95H