第3章有限长序列离散变换和快速算法

本章讨论有限长离散序列的正交变换，主要关注离散傅里叶变换(DFT)，详细讨论DFT的定义和性质以及DFT的快速算法FFT。
DFT及其快速算法FFT是数字信号处理中最重要的工具之一，应用非常广泛。本章集中在DFT的定义和性质、DFT的快速算法FFT。有关DFT的各项应用将在后续章节陆续介绍，第4章集中讨论以DFT为核心的数字频谱分析技术，第5章包括了DFT在滤波器实现中的应用，第6章介绍DFT在FIR滤波器设计中的应用，第8章给出了DFT在滤波器组的应用等。
本章也简要介绍了其他的正交变换，例如离散余弦变换(DCT)等。
3.1离散正交变换
假设x［n］只在0≤n<N区间内取非零值，该区间之外为零，这样的序列称为N点有限长序列。对于N点有限长序列，最多只有N个非零值，故只有N个自由度，可以在N维矢量空间里完整地表示该信号。在N维矢量空间里，定义N个相互正交的基序列(相当于基矢量)，对于0≤k≤N-1，每个基序列为

ak［n］,n=0,1,…,N-1(3.1)

每个基序列也只定义在区间0≤n<N内。为了由基序列对有限长序列定义正交变换，要求基序列满足正交性和完备性。
正交性

∑N-1n=0ak［n］a*k′［n］=Cδ［k-k′］(3.2)

完备性

∑N-1k=0ak［n］a*k［n′］=Cδ［n-n′］(3.3)

C是一个常数，当C=1时正交基是归一化的。
由于基序列是正交的，对每个基序列定义一个变换系数T［k］为

T［k］=∑N-1n=0x［n］ak［n］，0≤k≤N-1(3.4)

由完备性，可以证明由全部变换系数T［k］可以重构有限长信号为

x［n］=1C∑N-1k=0T［k］a*k［n］,0≤n≤N-1(3.5)

为得到式(3.5)，对式(3.4)两侧同乘a*k［m］，并对k从0~N-1求和，得到

∑N-1k=0T［k］a*k［m］=∑N-1k=0a*k［m］∑N-1n=0ak［n］x［n］=∑N-1n=0x［n］∑N-1k=0ak［n］a*k［m］

=C∑N-1n=0x［n］δ［n-m］=Cx［m］

上式中将m换为n，得到式(3.5)。




如前所述，由式(3.4)给出了变换的定义式，然后由完备性导出信号重构公式(3.5)(反变换公式)。容易证明若从式(3.5)出发，两边同乘ak′(n)对n求和，利用正交性式(3.2)，可得到式(3.4)，因此，式(3.4)和式(3.5)构成正、反变换对。
满足条件式(3.2)和式(3.3)的基序列有很多。给出每一种基序列，用式(3.4)和式(3.5)可定义一种离散变换，但不是每一种基序列和每一种离散变换都能得到广泛应用。其中有几种离散变换已得到广泛应用，应用最广泛的当属离散傅里叶变换(DFT)。对于按如下定义的基序列

ak［n］=e-j2πNkn，k=0,1,…,N-1； n=0,1,…,N-1

其正交性表示为

∑N-1n=0e-j2πNkne-j2πNrn*=Nδ［k-r］

完备性表示为

∑N-1k=0ej2πNknej2πNkl*=Nδ［n-l］

用X［k］表示变换系数的DFT变换对定义为


X［k］=DFTx［n］=∑N-1n=0x［n］e-j2πNkn，k=0,1,…,N-1(3.6)

x［n］=IDFTX［k］=1N∑N-1k=0X［k］ej2πNkn，n=0,1,…,N-1(3.7)

这是利用有限长序列正交展开引出的DFT变换的定义。由于DFT应用十分广泛，为了理解DFT定义的多方面性质，3.2节我们从其他不同方面重新引入DFT的定义。当然这些来自于不同方面对DFT的定义式是相同的，但却帮助读者从不同侧面理解DFT。
3.2离散傅里叶变换
在2.4.3节讨论周期序列的傅里叶变换时我们已经看到，周期序列的傅里叶变换是线谱形式，每个冲激的幅度由对单周期信号的DTFT在频点ωk=2πNk的采样值确定。由此可以得到结论，只需单周期信号DTFT在频点ωk=2πNk的采样值就可以准确表示一个周期信号。那么很自然，它也能够准确表示这个有限长的单周期信号。因此，对于有限长信号，我们并不需要DTFT的全部，而是只需要对DTFT在一些离散频点的采样值。
假设x［n］是有限长信号，只有在0≤n<N区间内取非零值，这样，只在0≤n<N区间内使用式(2.59)，并结合式(2.61)得到如下的变换对

X［k］=∑N-1n=0x［n］e-j2πNkn，k=0,1,…,N-1(3.8)

x［n］=1N∑N-1k=0X［k］ej2πNkn，n=0,1,…,N-1(3.9)

我们称这一变换对为离散傅里叶变换对(DFT)，其中式(3.8)是离散傅里叶变换对的正变换(DFT)，式(3.9)是离散傅里叶变换对的反变换(IDFT)。
当x［n］是有限长信号，只有在0≤n<N区间内取非零值，此区间之外取零值时，其DFT和DTFT关系为

X［k］=X(ejω)ω=2πN k=∑N-1n=0x［n］e-jωnω=2πN k(3.10)

由式(3.10)，我们可以看到，DFT变换系数X［k］有清晰的物理意义。
3.2.1DFT作为对DTFT的频域离散采样
为了进一步理清一个离散信号的DFT和DTFT之间的关系、明确DFT的物理意义和限制条件，以下我们从对一般离散信号的DTFT采样入手，重新审视DFT定义的含义和限制。
从一般信号x［n］出发，讨论对x［n］的DTFT X(ejω)采样问题。这里先不限制x［n］为有限长序列，在频点ωk=2πNk上对X(ejω)进行采样，在2π的一个周期内可以采得N个点，其他的采样值都是周期重复的。采样值写为

X［k］=X(ejω)ωk=2πN k=∑+∞n=-∞x［n］e-j2πNkn(3.11)

为了使用DTFT的反变换公式且仅用DTFT的采样值求取原信号，我们定义采样离散时间傅里叶变换为

X^(ejω)=∑+∞k=-∞2πNX［k］δω-2πNk(3.12)

将X^(ejω)代入DTFT反变换公式，仅利用DTFT的采样值得到的时域离散信号用x^［n］表示，有

x^［n］=12π∫2π0X^ejωejωndω

=12π∫2π0∑+∞k=-∞2πNX［k］δω-2πNkejωndω

=1N∑N-1k=0X［k］∫2π0δω-2πNkejωndω

=1N∑N-1k=0X［k］ej2πNkn(3.13)

根据冲激函数的抽样性质，式(3.12)的X^(ejω)也可以写成

X^(ejω)=X(ejω)2πN∑+∞k=-∞δω-2πNk(3.14)

由卷积定理和周期冲激串的傅里叶变换得

x^［n］=x［n］*∑+∞r=-∞δn-rN=∑+∞r=-∞xn-rN(3.15)

从式(3.15)可以看出，对于无穷长序列x［n］，仅由DTFT的采样值X［k］是无法得到原信号的，得到的是原信号的重叠相加后的一个周期化序列，这种现象称为时域混叠。只有当x［n］是一个有限长序列，其不为零的值集中在0≤n<M且M≤N时，式(3.15)中各xn-rN项对不同的r互不重叠，此时没有时域混叠发生。只要取出x^(n)中0≤n<N对应的一个周期的值，就可得到x［n］，即

x［n］= x^［n］RN［n］

总结如上的讨论： 当x［n］是长度不超过N的有限长序列，由对DTFT的采样所获得的X［k］可完全恢复信号x［n］。此时变换对式(3.8)和式(3.9)成立，称为DFT变换对。其中X［k］是对X(ejω)的采样值，由式(3.10)确定。当x［n］是长于N的序列时，由采样值X［k］所恢复的信号是如式(3.15)所示的混叠后的周期化序列，周期是N。
DFT变换有几个明显的特点。其一，它是可直接计算的。不像DTFT，由于ω取值的连续性，无法编写一个程序或通过一个数字系统直接计算出所有DTFT的值。当x［n］是N点有限长序列，直接通过式(3.8)计算出所有变换系数X［k］（k=0,1,…,N-1）需要N2次复数乘法和N(N-1)次复数加法。进一步的研究表明，DFT存在高效算法，这些高效算法统称快速傅里叶变换(FFT)，可以把计算量降低到Nlog2N的量级。
其二，DFT有明确的物理意义。它是对DTFT的采样，表示在频点ωk=2πNk处的频谱密度值。实际频谱分析的软件或仪器，大多采用DFT来计算频谱密度，再通过以X［k］为基点插值获得频谱密度X(ejω)。在这些软件或仪器中，是通过用X［k］插值来显示频谱图X(ejω)的。容易证明，对于N点序列x［n］，由X［k］通过插值可获得X(ejω)精确重构。插值公式为

Xejω=∑N-1k=0X［k］φω-2πNk(3.16)

其中

φω＝1Ne-jN-12ωsinN2ωsinω2

式(3.16)推导如下

X(ejω)=∑N-1n=0x［n］e-jωn=∑N-1n=01N∑N-1k=0X［k］ej2πNkne-jωn

=∑N-1k=0X［k］∑N-1n=01Nej2πNkne-jωn=∑N-1k=0X［k］1N
1-e-jNω-2πNk1-e-jω-2πNk


=∑N-1k=0X［k］1Ne-jN-12ω-2πNksinN2ω-2πNksin12ω-2πNk


由X［k］通过式(3.16)可精确重构X(ejω)，但在实际应用中，例如频谱分析的软件包或频谱分析仪中，显示X(ejω)的幅度谱和相位谱的显示器分辨率是有限的，并不需要完全精确的X(ejω)值，可以用一些更简单的插值方式由X［k］得到X(ejω)的近似值。例如，在N充分大时，可以用直线连接相邻X［k］的值，即用线性插值代替式(3.16)的精确插值(见图3.1)，或者将相邻几个X［k］的值用多项式插值来近似X(ejω)。



图3.1序列的DFT(垂直点线)和插值得到的DTFT(包络线)（图中只画出幅度谱）


其三，可以利用DFT来计算离散LTI系统的输出。在一定条件下，利用DFT的快速算法FFT，可获得更高效率的卷积计算方法。对该问题，在3.4节介绍DFT性质时再详细讨论，其在滤波器实现中的应用在第5章再做进一步分析。
3.2.2DFT的矩阵表示
DFT和IDFT可以用矩阵运算形式表示，我们定义信号矢量

x=x［0］,x［1］,…,x［N-1］T

变换系数矢量

X=X［0］,X［1］,…,X［N-1］T

变换矩阵

T=e-j2πNknN×N=11…1
1e-j2πN…e-j2πN(N-1)

1e-j2πN(N-1)…e-j2πN(N-1)(N-1)

反变换矩阵

U=1Nej2πNnkN×N=1N11…1
1ej2πN…ej2πN(N-1)

1ej2πN(N-1)…ej2πN(N-1)(N-1)

DFT的矩阵形式为

X=Tx(3.17)

IDFT的矩阵形式为

x=UX(3.18)

不难验证

UT=I(3.19)

和

U=1NTH(3.20)

式(3.19)表示了变换的完备性。通过式(3.17)的变换获得一组变换系数，通过式(3.18)的反变换重构原信号矢量，在这个过程中，能够通过变换和反变换准确重构原信号的条件就是式(3.19)。这可以通过将式(3.17)代入式(3.18)得到验证

x=UX=UTx(3.21)

由式(3.21)知，只有满足式(3.19)才能通过变换系数准确重构原信号矢量。
3.2.3DFT的实例
通过以下例子，进一步理解DFT的定义，尤其是通过例子说明DFT定义隐含的灵活性。
例3.2.1N点有限长信号

x［n］=RN［n］=1,0≤n<N

0,其他

做N点DFT，DFT系数记为X［k］，由DFT的定义得

X［0］=∑N-1n=0x［n］e-j2πNknk=0=∑N-1n=0x［n］=N

对于k≠0,1≤k<N，有

X［k］=∑N-1n=0x［n］e-j2πNkn=∑N-1n=0e-j2πNkn=0

其N点DFT系数除X［0］外均为零。由这组DFT系数，可以反变换得到原信号，对于0≤n<N，有

x［n］=1N∑N-1k=0X［k］ej2πNkn=1NX［0］ej2πN0n=1,0≤n<N

从另外一个角度理解DFT系数，2.4节得到x［n］的DTFT为

X(ejω)=e-j(N-1)ω/2sinNω/2sin(ω/2)

对DTFT在ω=2πNk处采样，首先考虑k=0的情况

X［0］=X(ejω)ω=0=N

对于k≠0,1≤k<N，有

X［k］=X(ejω)ω=2πN k=e-j(N-1)2π2N ksinN2π2Nksin2π2Nk=e-j(N-1)πN ksinπksinπNk=0

除了线性相位项e-j(N-1)ω/2外，X(ejω)重画于图3.2中，可以看到每个过零点位置对应Nω/2=kπ，即

ω=2πNk

处，这与DFT对DTFT的采样点重合，故除X［0］外，其余DFT系数均为0。


图3.2N点等值序列的DTFT示意图(忽略线性相位项)


在例3.2.1中，x［n］有N个非零值，至少要做N点DFT，才能由变换系数求得原序列的全部非零值。在实际中为了处理的方便，对于有N个非零值的序列，可以进行N1>N点的DFT。也就是说，在x［n］后面补上若干零，把x［n］看成是N1长度的序列进行DFT。这样，用N1点DFT系数除了可以计算出x［n］的全部非零值，还会把后面补的零值也计算出来，这在应用中没有什么危害。所以，对于只有N个非零值的序列，通过补零计算更长点数的DFT，是DFT应用中经常采用的一种办法。后续还会看到，DFT的这种灵活性，大大地扩宽了DFT的应用范围。
例3.2.2继续讨论例3.2.1中的x［n］序列，取N1=2N，对x［n］做N1点DFT。为了与上例对比，用x1［n］表示补零后的x［n］，用X1［k］表示对x［n］补零后的2N点DFT系数。因此

x1［n］=x［n］,0≤n<N

0,N≤n<2N

直接代入DFT定义式，得

X1［k］=∑2N-1n=0x1［n］e-j2π2Nkn=∑N-1n=0e-jπNkn

=N,k=0

1-e-jπNkN1-e-jπN k,1≤k<2N
=N,k=0

1-(-1)k1-e-jπN k,1≤k<2N

=N,k=0

0,k为偶数

e-jπ2-π2Nk1sinπ2Nk,k为奇数

请读者自行验证，用

X1［k］=X(ejω)ω=2π2N k

也得到相同的X1［k］表达式。对x［n］补上N个零，做2N点的DFT，相当于对X(ejω)在

ω=2π2Nk=πNk

处采样，采样密度比作N点DFT高了一倍，相当于对图3.2采样时，将两个过零点之间的峰值处也进行了采样，这正是X1［k］中k为奇数点的那些取值。

补零后做更长点数的DFT，一方面，可得到对DTFT更密集(更细致)的采样； 另一方面,通过反变换仍准确重构信号值，这给DFT带来灵活性。但也注意到，既然通过N点DFT就可准确重构信号和插值得到X(ejω)，因此，理论上补零做更长点数DFT不会带来更多的有用信息，但在实际应用中会带来很多方便性。至于DFT的这种灵活性带来的方便，本章后文还会多次遇到。
3.3DFT与周期序列傅里叶级数的关系
为讨论方便，重写离散周期序列的傅里叶级数(DFS)和DFT以便于比较。周期序列x～［n］=∑+∞r=-∞x［n-rN］，周期为N，其傅里叶级数系数为X～［k］。尽管只关心X～［k］的N个取值，但它也是周期为N的序列。用简化符号WN=e-j2πN，离散周期序列的傅里叶级数展开重写如下

X～［k］=∑N-1n=0x～［n］WknN(3.22)

x～［n］=1N∑N-1k=0X～［k］W-knN(3.23)

考虑有限长序列x［n］，长度为N，序号范围为0≤n≤N-1，其DFT重写为

X［k］=DFTx［n］=∑N-1n=0x［n］WknN，0≤k≤N-1(3.24)

x［n］=IDFTX［k］=1N∑N-1k=0X［k］W-knN，0≤n≤N-1(3.25)

DFS和DFT仅从数学形式上是一致的，但两者的目标不同。DFT的定义中，只关心有限长序列，它本质上是对DTFT的采样。由此，其可以看成是时域/频域双域采样得到的“纯离散变换”。尽管DFT定义时仅关心有限长序列，但是因为时域采样带来的频域周期拓展，使得DTFT自身就是周期的，作为DTFT采样的X［k］自然具有周期性(尽管不关心0≤k≤N-1之外的值)； X［k］作为对DTFT的采样，又引起由X［k］重构的信号构成周期拓展，尽管只关心由式(3.25)在0≤n≤N-1内重构x［n］，但式(3.25)也隐含地满足周期性，即在0≤n≤N-1之外按周期N重复x［n］的取值。所以，尽管DFT是仅就有限长序列定义的，但其存在自然的周期延拓性，即X［k］隐含周期性。
为后面讨论方便，这里给出将有限长序列周期化和从周期序列取出一个周期形成有限长序列的表达方式。由有限长序列x［n］，得到周期序列x～［n］可表示为

x～［n］=∑+∞r=-∞x［n+rN］=x((n))N(3.26)

这里((·))N表示对N取余数运算，例如x((r+kN))N=x［r］。由周期序列，取出一个周期表示为

x［n］= x～［n］RN［n］=x((n))NRN［n］(3.27)

我们在讨论DFT性质时，尤其是与反转、位移等相关联的性质时，有限长序列位移或反转后序号超出0≤n≤N-1，但由于DFT对0≤n≤N-1之外没有定义，不方便于这些性质的理解。为了处理这些性质，可先将有限长序列拓展成周期序列，进行反转或位移后再取其一个周期。将看到，在表述DFT的一些性质时，就是用的这种方式。
至此，已讨论了各种不同的傅里叶变换。各种傅里叶变换在时域和频域有一些不同的对偶关系，把理想采样信号的傅里叶变换也包括在内，几种傅里叶变换时域和频域的一些对偶关系总结在表3.1中，时域和频域的一般对偶性总结在表3.2中。



表3.1各种傅里叶变换时域和频域对偶关系



变 换 名 称时 域 性 质频 域 性 质

连续傅里叶级数周期连续离散，无限
连续傅里叶变换连续，无限连续，无限
采样信号离散，无限周期连续
DTFT离散，无限周期连续,频域归一化离散傅里叶级数离散，周期离散，周期DFT离散，有限长离散，有限长


表3.2时域和频域一般对偶性



时域 频域 
连续，无限连续，无限
离散，无限 周期连续
周期连续离散，无限 
离散，周期 离散，周期

3.4DFT的性质
DFT有许多性质，其中一些是简单和自明的，另一些性质却并不直观，需要仔细地解释。对于一些简单性质，只简单列出，证明留作练习。对于几个需要仔细解释的性质，给出详细的讨论。
在叙述DFT性质之前，为了在各性质叙述中节省重复声明的篇幅，假设X［k］=DFTx［n］，X1［k］=DFTx1［n］，X2［k］=DFTx2［n］。这里DFTx［n］表示对x［n］求离散傅里叶变换，若不加说明，表示N点DFT。
性质1线性性质
若x［n］=a1x1［n］+a2x2［n］，则X［k］=a1X1［k］+a2X2［k］。线性性质中，各序列均取N点DFT，若x1［n］ 和x2［n］长度不同，则N代表较长序列的长度，较短的序列通过补零也做N点DFT。
性质2反转性
对于定义在0≤n≤N-1的N点有限长序列，如何描述其反转性？信号的反转对应于DFT系数如何变化？对有限长序列的反转性，既要有“反转”之表现，又要保持在0≤n≤N-1区间。 为了研究反转性，首先将x［n］周期化，然后反转，再取其一个周期值，相当于

x((-n))NRN［n］=x［0］,n=0

x［N-n］,n≠0(3.28)

反转的表现为x［0］不变，其他值在取值范围内依次前后对换。为了使表示式更简洁，可将反转运算简记为

x((-n))NRN［n］=x［N-n］(3.29)

例3.4.1序列反转的示意图如图33所示，图中示出了由x［n］=5,4,3,2,1周期延拓得到x((n))N，然后反转得到x((-n))N，最后仅取第一个周期得到x((-n))NRN［n］，最后得到的反转序列如式(3.28)所示。


图3.3有限长序列反转的示意图


DFT变换的反转性叙述为

DFTx((-n))NRN［n］=X((-k))NRN［k］

或简记为

X［N-k］=DFTx［N-n］(3.30)

即信号反转对应其DFT系数反转。
证明

DFTx-nN=∑N-1n=0x-nNe-j2πNnk=x0+∑N-1n=1xN-ne-j2πNnk

=x0+∑N-1n′=1xn′e-j2πNn′((-k))N=∑N-1n′=0xn′e-j2πNn′((-k))N

=X-kN

上式的推导过程中，由于求和号已经限制了序列的取值范围，故省略了RN［n］和RN［k］。
对于N个非零值的有限长序列，也可通过补零把序列看作N1>N长的序列，若对该补零序列做N1点DFT，如上反转性质用N1替代N进行描述。


图3.4序列补零做反转

例3.4.2在例3.4.1中，若对序列补上5个零，将序列作为x［n］={5,4,3,2,1,0,0,0,0,0}的N1=2N长序列进行反转，其反转x((-n))10R10［n］如图3.4所示。

性质3对偶性(duality)
如果有限长序列xn的N点DFT为Xk，把Xn看成时域离散序列，则Xn的N点DFT为

DFTXn=∑N-1n=0Xne-j2πNnk=∑N-1n=0∑N-1m=0xme-j2πNmne-j2πNnk

=∑N-1m=0xm∑N-1n=0e-j2πNm+kn
=Nx-kN

可见，如果xn(时域)和Xk(频域)构成一对离散傅里叶变换对，则Xn(时域)和Nx-kNRN［k］之间构成一对离散傅里叶变换对。
利用对偶性，只要知道了某个有限长序列的DFT，则与其频域序列相同的时域序列的DFT也就得到了，无须再重新计算。

性质4循环位移性
循环位移性表述为下式的形式

DFTx((n-m))NRN［n］=WkmNX［k］(3.31)

这里x((n-m))NRN［n］称为循环位移。设m>0，x((n-m))N表示周期化后右移m点，RN［n］取出一个周期。由于x((n-m))N的周期性，不难看出，x((n-m))NRN［n］相当于x［n］依次从0≤n≤N-1的尾部移出一个数据，又从头部移进来，直到移动完m次。这相当于把0≤n≤N-1序列绕在一个圆柱上右移m次，所以称这种移位为循环移位。m<0时是相反的移动，数据从头部移出从尾部移入，是循环左移。不管移动方向如何，循环移位对应DFT系数乘以WkmN=e-j2πNkm的指数项。


图3.5循环位移

例3.4.3序列x［n］=5,4,3,2,1按N=5进行循环位移，位移量m=2的示意图如图3.5所示。
上述移位的过程之所以被称为“循环”移位，因为可以形象地用图3.6来表示。图3.6(a)中有一个圆桶，四周被均匀地划分成N个格，依次刻上有限长序列的值。把序列第一点x0对准0刻度压到一根直线轴上，随着圆桶转动，刻在桶上的序列被复制到直线轴上，就形成了周期延拓后的序列。
循环移位，就是先把圆桶旋转指定的步数，正延迟逆时针旋转，负延迟顺时针旋转，然后再把旋转后的圆桶对准0刻度压下去，随着圆桶转动，刻在桶上的循环移位后的序列被复制到直线轴上，就形成了延迟后的序列的周期延拓。图3.6(b)是延迟为2时的示意图。


图3.6循环移位的圆桶示意图


证明
取x1［n］=x((n-m))NRN［n］，直接代入DFT定义运算如下

DFTx1n=∑N-1n=0x1ne-j2πNnk=∑N-1n=0xn-mNe-j2πNnk

=∑((m))N-1n=0xn+N-((m))Ne-j2πNnk+∑N-1n=((m))Nxn-((m))Ne-j2πNnk

=∑N-1n′=N-((m))Nxn′e-j2πNn′ke-j2πN((m))Nk+∑N-((m))N-1n′=0xn′e-j2πNnke-j2πN((m))Nk

=∑N-1n′=0xn′e-j2πNn′ke-j2πNmk=Xke-j2πNmk

这里可能会问，DFT怎样反映自然移位性(即延迟特性)，即对N点有限长序列x［n］，x［n-m］对应的DFT系数怎样？若x［n］在0≤n≤N-1范围内的取值均非零，由于x［n-m］的部分数据已移出0≤n≤N-1范围，无法直接用N点DFT表示，故无法用N点DFT性质表述自然移位x［n-m］。但是，这并不是说DFT就无法表示自然移位了，可以通过DFT定义的灵活性解决这个问题。
设有一个有限长信号x［n］在0≤n≤N-1范围内的取值均非零。若应用中需要用到m<M的自然移位x［n-m］的DFT域表示，这只要把x［n］看成是长度为N1=N+M的序列，在N≤n<N+M范围x［n］=0，对x［n］做N1点DFT。这种做法在DFT应用中很常用，也就是对x［n］补M个零从而构成长度为N1=N+M的序列，再做N1点的DFT。这样，当延迟不超过M时，x［n-m］和循环移位x((n-m))N1RN1［n］相等，注意到，现在循环移位的循环周期是N1，这样，x［n-m］对应的DFT系数是WkmN1X［k］，这里的X［k］是x［n］补零后的N1点DFT系数。


图3.7补零后的循环位移

例3.4.4将例3.4.3中的序列补零，得到x［n］=5,4,3,2,1,0,0,0,0,0的序列，按N1=10进行循环位移，得到图3.7的位移序列，这种情况下，循环位移等于自然位移。
我们将会看到，DFT的基本性质大多建立在“循环变化”基础上，但通过补零，DFT可用于描述各种自然的信号处理算法。
性质5频域循环位移性
这是性质4的对偶性质，仅表示如下

WrnNx［n］X((k+r))NRN［k］(3.32)

性质6循环卷积性
假设x1［n］和x2［n］均为N点序列时，X1［k］和X2［k］分别是两个序列的N点DFT。若取X［k］=X1［k］X2［k］，其反变换得到的序列为

x［n］=IDFTX［k］=∑N-1m=0x1［m］x2((n-m))N(3.33)

式(3.33)的求和项称为循环卷积(和)。和式中的第二项是循环反转和位移运算，尽管只关心0≤n≤N-1范围的结果，但实际上该求和项的结果是周期性的。严格讲，式(3.33)应该写为

x［n］=IDFTX［k］=∑N-1m=0x1［m］x2((n-m))NRN(n)

由于我们清楚x［n］只取0≤n≤N-1范围的值，若为简单计可省略RN(n)项。


证明

x［n］=IDFTX1［k］X2［k］=1N∑N-1k=0X1［k］X2［k］ej2πNnk

=1N∑N-1k=0∑N-1m=0x1［m］e-j2πNmkX2［k］ej2πNnk

=∑N-1m=0x1［m］1N∑N-1k=0X2［k］e-j2πNmkej2πNnk

=∑N-1m=0x1［m］x2((n-m))N

注意到上式第3行相当于求X2［k］e-j2πNmk的反变换，故用性质4(循环移位性)，得到了第4行的结果。
两个序列DFT乘积的IDFT是两个序列的循环卷积，循环卷积过程可用图3.8形象地表示。图中，两个内外相嵌的圆桶分别被划分成N个格，一个序列逆时针排列在其中一个圆桶上，另一个序列顺时针排列在另一个圆桶上。把按照顺时针排列的圆桶按照逆时针方向每次转动一格(或者等价地把按照逆时针排列的圆桶按照顺时针方向每次转动一格)，对应格的序列值相乘然后求和，就得到循环卷积结果。


图3.8循环卷积的圆桶示意图


循环卷积实际上是两个序列以N为周期延拓之后在一个周期上的卷积和。式(3.33)的求和式称为循环卷积，用符号x1［n］○N x2［n］表示。如下讨论循环卷积与线性卷积的关系。为了区分两种卷积，用如下符号表示线性卷积

xL［n］=x1［n］*x2［n］=∑+∞m=-∞x1［m］x2［n-m］

=∑Nm=0x1［m］x2［n-m］(3.34)

假设x1［n］和x2［n］均为N点序列，这是为了使两个序列均可做N点DFT，实际上可以放宽到x1［n］真正非零长度为L≤N，x2［n］真正非零长度为P≤N，通过补零对x1［n］和x2［n］均做N点DFT。由第2章介绍的线性卷积性质知，两个有限长序列x1［n］和x2［n］的线性卷积xL［n］不为零的长度为L+P-1。可以证明： 循环卷积与线性卷积的关系为

x［n］=∑+∞r=-∞xL［n-rN］RN［n］(3.35)

这说明，N点循环卷积是线性卷积以周期为N的周期性延拓叠加结果。为了证明式(3.35)的正确性，对式(3.34)两侧做DTFT，得到

XL(ejω)=X1(ejω)X2(ejω)(3.36)

对式(3.36)两侧均在ωk=2πNk处采样，得

XL(ejω)ω=2πN k=X1(ejω)ω=2πN kX2(ejω)ω=2πN k

右侧的采样值分别对应X1［k］和X2［k］，故

X［k］=X1［k］X2［k］=XL(ejω)ω=2πN k

因此，X［k］相当于是对xL［n］的DTFT XL(ejω)在ωk=2πNk的采样值，由DFT定义时导出的结论式(3.15)得到式(3.35)。
先来看循环卷积和线性卷积最不同的情况。当x1［n］和x2［n］非零长度均为N，这时xL［n］是长度为2N-1的序列。既然X［k］相当于是对XL(ejω)的N点采样在ωk=2πNk的采样值，由DFT定义时的讨论得知，由X［k］无法重构2N-1长的序列xL［n］，而是得到混叠的结果。这个结果如式(3.35)所示，由X［k］做IDFT得到的x［n］中几乎所有值都是混叠了的(只有一点的值是对的)。
上述讨论说明，无法直接用DFT计算线性卷积。但是通过对序列补零，利用DFT的灵活性，是可以用DFT计算有限长序列的线性卷积的。这里讨论两种特殊情况。
第一种情况，设x1［n］和x2［n］长度不等，x1［n］长度为L，x2［n］长度为P，且L>P。取N=L，x1［n］直接做N点DFT，x2［n］补零后也做N点DFT，对X［k］=X1［k］X2［k］做N点IDFT得到x［n］。在这些条件下，通过式(3.35)xL［n］和x［n］的关系示于图3.9中。


图3.9线性卷积和与循环卷积和的关系


如图3.9所示，矩形框住的值就是IDFT的结果x［n］，而三角图形表示线性卷积xL［n］和它的各移位求和项。在该情况下，式(3.35)中只有xL［n］和xL［n+N］两项对矩形窗内的值有影响，两个三角形重叠的区域x［n］的取值由xL［n］和xL［n+N］混合而成。观察图3.9发现，x［n］在0≤n≤P-2范围内的P-1个值是混叠结果，而x［n］在P-1≤n≤N-1范围内的L-P+1个值是与线性卷积相等的。这种情况说明，当x1［n］和x2［n］长度不等时，若以长序列的点数做DFT，利用IDFT计算得到的循环卷积中保留了线性卷积的部分正确结果，即P-1≤n≤N-1范围内循环卷积保留了线性卷积的部分正确值。这个结论在FIR滤波器计算时会得到应用。

第二种情况，仍设x1［n］长度为L，x2［n］长度为P。取N≥L+P-1，对x1［n］和x2［n］均补零，做N点DFT，计算X［k］=X1［k］X2［k］的IDFT，得到循环卷积x［n］。在这种条件下，因满足

x［n］=∑+∞r=-∞xL［n-rN］RN［n］＝xL［n］

循环卷积等于线性卷积。由此，对于有限长序列，用循环卷积计算线性卷积的步骤为取N≥L+P-1，计算如下： 
(1) x1［n］在P~N-1补零，做N点DFT； 
(2) x2［n］在L~N-1补零，做N点DFT； 
(3) X［k］=X1［k］X2［k］，0≤k≤N-1； 
(4) xL［n］=x［n］=IDFTX［k］，0≤n≤N-1。
用循环卷积计算线性卷积有其实际意义。若把上述讨论中的x1［n］换作一个离散LTI系统的单位抽样响应，把x2［n］看作系统的输入信号，则线性卷积就是系统的输出。利用循环卷积也可计算系统的输出。由于DFT存在快速算法FFT，一定条件下用DFT计算系统输出可能更经济。
在实际应用中还会遇到LTI系统的单位抽样响应是有限长的，即FIR系统，但系统输入信号是无限长序列。这种情况下，可将输入序列分成有限长的段进行处理，具体算法将在第5章进一步讨论。

图3.10循环卷积

例3.4.5有两个信号，x1［n］=1,1,1,1,1，x2［n］=5,4,3,2,1，取N=5，进行循环卷积，首先利用式(3.33)的定义直接计算循环卷积，计算过程的示意图如图3.10所示，记循环卷积的结果为y［n］，显然y［n］=15,15,15,15,15。
图3.10(c)是x2((-k))5，注意到

y［0］=∑4k=0x1［k］x2((-k))5=15

对于1≤n<5，图3.10(c)进行循环位移，得到x1［k］x2((n-k))5并进行求和，但由于x1［k］=1，所以循环卷积结果总是x2［n］非零值之和，恒为15。
以上结果同样可以通过DFT求得。N=5时，例3.2.1已得到x1［k］的DFT系数为： 只有X1［0］=5非零，故也只需要求出X2［0］=15，因此

X1［k］X2［k］=75,0,0,0,0

求IDFT为

y［n］=15∑4k=0X1［k］X2［k］W-kn5=15×75W05=15

例3.4.6离散序列与例3.4.5相同，本例取N=10进行循环卷积，补零后两个信号为

x1［n］=1,1,1,1,1,0,0,0,0,0，x2［n］=5,4,3,2,1,0,0,0,0,0

直接计算循环卷积，示意图如图3.11所示。图3.11(b)和图3.11(c)给出了x2((n-k))10在n=0和n=2的示意图，图3.11(d)是循环卷积的结果。可见对本例，取N=10时，计算得到的循环卷积等于线性卷积。本例满足循环卷积等于线性卷积的最小N值是9。
同样的结果，可通过DFT求得，做N=10点DFT，得

X1［k］={5, 1 -j 3.0777, 0,  1 -j 0.7265, 0,  1,  0,  1 + j0.7265, 0, 1 +j 3.0777}

X2［k］={15, 7.7361 -j 7.6942, 2.5 -j 3.4410,  3.2639 -j 1.8164i, 

2.5 -j 0.8123i3,  2.5 +j 0.8123,  3.2639 + j1.8164,  2.5 +j 3.4410,  

7.7361 +j 7.6942}

相乘后做IDFT得

y［n］=IDFTX1［k］X2［k］={5,9,12,14,15,10,6,3,1,0}



图3.11补零后循环卷积等于线性卷积的例子


性质7频域卷积性
设序列x［n］=x1［n］x2［n］，其中x［n］、x1［n］、x2［n］均为N点有限长序列，并分别做N点DFT，则

X［k］=DFTx［n］=∑N-1l=0X1［l］X2((k-l))N

性质8共轭性质

DFTx*［n］=X*［N-k］

证明

DFTx*［n］=∑N-1n=0x*［n］e-j2πNkn=∑N-1n=0x［n］e-j2πN(-k)n*

=∑N-1n=0x［n］e-j2πN((-k))Nn*=∑N-1n=0x［n］e-j2πN(N-k)n*=X*［N-k］

性质9DFT的对称性
利用有限长序列的周期延拓，定义有限长序列的周期共轭对称性为

xep((n))NRN［n］=x*ep((-n))NRN［n］

记住，除零点自身互共轭(一定是实数)外，周期共轭对称关系可简写为

xep［n］=x*ep［N-n］

类似地，周期共轭反对称分量简写为

xop［n］=-x*op［N-n］

对于一个任意序列x［n］，可分解为周期共轭对称序列和周期共轭反对称序列之和。周期共轭对称序列和周期共轭反对称序列与x［n］的关系为

xep［n］=12x［n］+x*［N-n］

xop［n］=12x［n］-x*［N-n］

如果x［n］的DFT分别写成其实部和虚部之和，即X(k)=XR(k)+jXI(k)，则有

XR［k］=DFTxep［n］

jXI［k］=DFTxop［n］

对偶的关系表述为： 若x［n］=xr［n］+jxi［n］为复序列，xr［n］是实部，xi［n］为虚部，则

DFTxr［n］=Xep［k］=12X［k］+X*［N-k］

DFTjxi［n］=Xop［k］=12X［k］-X*［N-k］

例3.4.7利用对称性可以简化DFT的计算。例如，有两个实值信号x1［n］和x2［n］，用它们分别作为实部和虚部构成一个复信号x［n］=x1［n］+jx2［n］。直接计算得到复信号的DFT X［k］，利用对称性得到x1［n］和x2［n］的DFT为

X1［k］=DFTx1［n］=Xep［k］=12X［k］+X*［N-k］

X2［k］=DFTx2［n］=1jXop［k］=12jX［k］-X*［N-k］

性质10实序列DFT的对称性
这是性质9的特例，因为重要，故单独列为一条性质。若离散信号是实值信号，其信号自身就是其实部。这说明，实值信号的DFT自身满足周期共轭对称性，即

X［k］=X*［N-k］(3.37)

也就是说，实值有限长序列的DFT系数满足X［0］=X*［0］，故X［0］是实值； 其他DFT系数在1≤k≤N-1范围内，满足式(3.37)的周期共轭对称性质。若N是偶数，也可确定X［N/2］=X*［N/2］是实数。
对一个N点实序列，其N点DFT系数存在共轭对称性，仅由X［0］、X［N/2］和1≤k≤N/2-1的系数即可确定全部系数。这样，可以用N/2个复数存储单元即可存储N点实序列的DFT系数。由于存储一个复数需要两个实数存储单元，因此存储实序列和存储其DFT系数需要相同数目的存储单元。

性质11帕塞瓦尔定理
如果有限长序列xn的N点DFT为Xk，帕塞瓦尔定理(Parsevals Theorem)反映变换前后序列的能量关系

∑N-1n=0xn2=1N∑N-1k=0Xk2

证明

∑N-1n=0xn2=∑N-1n=0xnx*n=∑N-1n=01N∑N-1k=0Xkej2πNnkx*n

=1N∑N-1k=0Xk∑N-1n=0xne-j2πNnk*=1N∑N-1k=0XkX*k

=1N∑N-1k=0Xk2

DFT有许多性质，其中一些性质很独特。灵活使用DFT的各种性质，可以极大地方便DFT的各类应用。
3.5用DFT计算相关序列
如果一有限长离散信号x［n］，在0≤n≤N-1内非零，根据自相关的定义，自相关序列为

rxx［k］=∑+∞n=-∞x［n］x*［n-k］=∑N-1n=kx［n］x*［n-k］,0≤k≤N-1
∑N-1+kn=0x［n］x*［n-k］,-N+1≤k<0
0,k≥N

对于有限长序列来讲，信号总是能量有限的，故用如上的能量信号的自相关定义。自相关取值不为零的区间为k≤N-1，非零长度为2N-1。
第2章已经详细研究了自相关的DTFT，假设

Exx(ejω)=DTFTrxx［k］=X(ejω)2

由于自相关是2N-1长度的，只要取L≥2N-1，在ωk=2πLk对Exx(ejω)=X(ejω)2进行采样，由采样值可重构自相关序列。故得到用DFT计算自相关序列的算法如下。
算法1： 
(1) 取L≥2N-1，离散信号x［n］补零，得到

xL［n］=x［n］,0≤n≤N-1

0,N≤n≤L-1

(2) 对xL［n］做L点DFT，得到XL［k］，0≤k≤L-1；
(3) 计算Exx［k］=XL［k］2，0≤k≤L-1；
(4) 计算IDFT，得到rL［n］=IDFTXL［k］2，0≤n≤L-1；
(5) 得到自相关序列为

rxx［k］=rL［k］,0≤k≤N-1
rL［L+k］,-N+1≤k<0
0,k≥N

可通过类似的讨论，给出互相关的计算步骤如下。
算法2： 
(1) 取L≥2N-1，离散信号x［n］、y［n］补零，得到

xL［n］=x［n］,0≤n≤N-1

0,N≤n≤L-1

yL［n］=y［n］,0≤n≤N-1

0,N≤n≤L-1

(2) 对xL［n］做L点DFT，得到XL［k］，0≤k≤L-1，对yL［n］做L点DFT，得到YL［k］，0≤k≤L-1； 
(3) 计算Exy［k］=XL［k］Y*L［k］，0≤k≤L-1； 
(4) 计算IDFT，得到rxyL［n］=IDFTExy［k］，0≤n≤L-1； 
(5) 得到互相关序列为

rxy［k］=rxyL［k］,0≤k≤N-1

rxyL［L+k］,-N+1≤k<0

0,k≥N


在研究了DFT的快速算法FFT后，自相关和互相关计算中的DFT和IDFT都采用FFT算法高效实现。
3.6DFT的快速计算方法
由DFT的定义

X［k］=∑N-1n=0x［n］WnkN(3.38)

计算每一个DFT系数需要N次复数乘法和N-1次复数加法，本节不加说明时，乘法和加法次数均指复数运算。这样计算所有的DFT系数需要N2次复数乘法和N(N-1)次复数加法。一次复数乘法需要四次实数乘法和两次实数加法，一次复数加法需要2次实数加法。所以，直接计算DFT所需要的实数乘法次数为4N2，所需要的实数加法次数为N4N-2。这样的算法复杂度我们称为N2量级的，简记为ON2。在数据采集速率比较高，N取值较大时，实时计算DFT的运算量仍是很可观的。本节研究DFT的快速计算方法： 快速傅里叶变换(FFT)。
DFT定义中的WnkN项，有许多性质可用于简化DFT的计算。首先，一些特殊值可节省运算，如W0N=WkNN=1，WN/2N=-1，这些项都节省乘法运算； 另外一些项如WN/4N=-j，W3N/4N=j在复数乘法运算时，并不需要实际做乘法，只是交换实部和虚部，可统计为0次乘法，这些特殊项均可节省乘法次数。另外一些性质，主要是周期性、齐次性和对称性，与序列分解结合可显著节省运算量。
周期性为

WnkN=W(N+n)kN=Wn(k+N)N

齐次性为

WkN=Wk/mN/m

变换基的对称性表现为复共轭对称性(complex conjugate symmetry)

W(N-n)kN=W-nkN=WnkN*

利用变换基的周期性和对称性可以从一定程度上提高计算效率，但是并不是数量级上的改进。提高DFT计算效率的另一类方法是把长序列分解成短序列，先做短序列的DFT，然后再在此基础上进一步得到长序列的DFT。
应用变换基的这些性质可方便地将长序列分解为短序列进行计算。那些可直接节省乘法运算次数的项，对于长序列DFT的直接计算来讲，节省的乘法次数并不明显，但对于短序列可以明显节省乘法次数。例如N=2或N=4时，可不需要任何乘法运算。为了看清楚这一点，列出如下几个短序列直接做DFT的计算公式。
2点DFT的直接计算，不需要乘法，只需要两次加法，计算如下

X［k］=∑1n=0x［n］Wnk2=x［0］W02+x［1］Wk2，k=0,1

计算两个系数的矩阵形式为

X［0］
X［1］=W02W02
W02W12x［0］
x［1］=11
1-1x［0］
x［1］=x［0］+x［1］
x［0］-x［1］

这个计算过程可表示为图3.12的蝶形结构。
3点DFT的计算如下

X［k］=∑2n=0x［n］Wnk3=x［0］W03+x［1］Wk3+x［2］W2k3，k=0,1,2

矩阵形式为

X［0］
X［1］
X［2］=W03W03W03
W03W13W23
W03W23W43x［0］
x［1］
x［2］=111
1W13W23
1W23W13x［0］
x［1］
x［2］

3点DFT的蝶形结构示于图3.13。
4点DFT计算为

X［k］=∑3n=0x［n］Wnk4=x［0］W04+x［1］Wk4+x［2］W2k4+x［3］W3k4，k=0,1,2,3




图3.122点DFT的蝶形计算




图3.133点DFT的蝶形计算




矩阵形式为

X［0］
X［1］
X［2］
X［3］=W04W04W04W04
W04W14W24W34
W04W24W44W64
W04W34W64W94x［0］
x［1］
x［2］
x［3］=1111
1-j-1j
1-11-1
1j-1-jx［0］
x［1］
x［2］
x［3］

当DFT算法要求的变换点数N为组合数时，既然N为组合数，比如N=N1N2，那么n可以表示为

n=n2N1+n1,n1=0,1,…,N1-1;n2=0,1,…,N2-1

此时

WnkN=W(n2N1+n1)kN=Wn2N1kNWn1kN

而

Wn2N1kN=e-j2πNn2N1k=e-j2πN2n2k=Wn2kN2

是N2点DFT变换的基。变换基的以上性质对k同样适用。
利用这种组合数性质，可以把N点长序列的变换基转变成N1或N2点短序列的变换基。经过适当的分解与合成，将对长序列的DFT计算问题转换为计算短序列的DFT然后再组合，这样做可能大大降低运算量。尤其是短序列因子Ni为2或4时，运算效率更高。这类算法统称为快速傅里叶变换(fast fourier transform，FFT)。注意，FFT不是一类新变换，而是DFT快速计算算法的总称。
快速傅里叶变换可以把DFT的计算量减少到ONlgN量级。由于计算量小，FFT得到广泛应用。
但是标准FFT算法必须要把所有Xk都计算出来，对于某些只需要一部分谱的应用来说，特别是当点数N比较大时，FFT算法在效率方面并没有优势。此时采用卷积实现DFT的线性调频z变换(chirpz)算法可以减少无效计算。
在对序列做滑窗(slidingwindow)分析时，如果上次的FFT结果已知，则分析窗滑动一步，只有一个最老的值被滑出窗外，同时只有一个新值被滑进窗内，此时可以利用上次FFT结果递推得到新的FFT结果。这种算法可以把每一段变换的计算量减小到ON。
3.6.1按时间抽取基2FFT算法
既然目的是将长序列分解成短序列，先进行短序列的DFT，将短序列的DFT合成为长序列的DFT，很自然的一种分解方法是将序列分成两部分，偶数序号为一个新的短序列，奇数序号组成第二个短序列，即

e［n］=x［2n］,


f［n］=x［2n+1］,0≤n≤N2-1 (3.39)

这里e［n］和f［n］都是N/2点序列，对其分别DFT，记为E［k］和F［k］，均为N/2点DFT。考虑到利用短序列DFT合成长序列DFT时，序号k在0≤k≤N-1范围取值，利用DFT隐含的周期性，用E((k))N/2和F((k))N/2分别表示E［k］和F［k］的周期延拓。这样

X［k］=∑N-1n=0x［n］WnkN
=∑N2-1r=0x［2r］W2rkN+∑N2-1r=0x［2r+1］W(2r+1)kN

=∑N2-1r=0e［r］WrkN/2+WkN∑N2-1r=0f［r］WrkN/2

=E((k))N/2+WkNF((k))N/2(3.40)

上式第三行是N/2点DFT，考虑到k在0≤k≤N-1范围取值，第四行用了N/2点DFT E［k］和F［k］的周期延拓形式。
式(3.40)给出了用短序列变换合成长序列变换的公式，为了方便，取0≤k≤N2-1，重写式(3.40)最后一行为

X［k］=E［k］+WkNF［k］，

Xk+N2=E［k］-WkNF［k］,0≤k≤N2-1(3.41)

由式(3.41)更清楚地看到，由于WkNF(k)只需要计算一次，若已经计算出N/2点DFT E［k］和F［k］，由E［k］和F［k］得到X［k］,0≤k≤N-1，仅需要N次加法和N/2次乘法。图3.14中，以N=8为例，画出了由两个N/2点DFT E［k］和F［k］合成得到X［k］的示意图。


图3.148点序列第一次分解流图


图3.14中看到，只有奇序列的变换F［k］输出端需要乘因子WkN，该因子称为旋转因子。一次分解后总运算量为
乘法次数
2N22+N2=N2N+1≈N22
加法次数
2N2N2-1+N=N22
乘法次数里的近似符号对大的N成立，可以看到一次分解获得的效果是大约降低了一半的乘法次数和加法次数，因此，这种分解可以继续进行下去，将e［n］和f［n］，按照奇偶继续划分为

a［n］=e［2n］=x［4n］，
b［n］=e［2n+1］=x［4n+2］，
c［n］=f［2n］=x［4n+1］，
d［n］=f［2n+1］=x［4n+3］，0≤n≤N4-1

对以上序列分别做N/4点DFT，类似式(3.41)的推导，得到 

E［k］=A［k］+WkN/2B［k］=A［k］+W2kNB［k］，
Ek+N4=A［k］-WkN/2B［k］=A［k］-W2kNB［k］,0≤k≤N4-1(3.42)

F［k］=C［k］+WkN/2D［k］=C［k］+W2kND［k］,


Fk+N4=C［k］-WkN/2D［k］=C［k］-W2kND［k］,0≤k≤N4-1(3.43)

以N=8为例，图3.15画出第二次分解后的示意图。


图3.158点序列第二次分解流图


这个分解过程一直进行下去。若取N=2m，经过m-1次分解后，序列已经变成2点序列。2点序列的DFT只需要2次加法，直接实现，不需要再分解。对于N=8的情况，最后的2点DFT直接用图3.12的蝶形图实现，最后得到的完整图形如图3.16所示。图3.16中最左侧的蝶形中的旋转因子W0N是为保持图中各级运算单元的一致性加入的。


图3.168点序列按时间抽取FFT流图


图3.16中，每一级的基本运算单元都是一致的，是一种蝶形运算单元，重新画在图3.17中。每个蝶形运算单元需要1次乘法和2次加法运算。在N=2m的情况下，总共有 m级运算，每一级由N/2个蝶形运算单元构成。故总运算量为
乘法次数

mc=m×N2=N2log2(2m)=N2log2(N)(3.44)

加法次数

ac=m×N=Nlog2(2m)=Nlog2(N)(3.45)



图3.17按时间抽取FFT的蝶形计算结构


观察图3.16的计算流程图，发现图中最左侧信号的输入不再是自然的顺序。图中所示的这种顺序称为倒位序。为了理解倒位序的一般性顺序，需分析序号的二进制表示。设N=2m为有限长信号的长度，需要用m位二进制数表示样本的序号n，n的二进制表示为

n=bm-1bm-2…b1b02=bm-12m-1+bm-22m-2+…+b121+b020

上式中，bi只取0和1。序号n的倒位序 n-记为

n-=b0b1…bm-2bm-12=b02m-1+b12m-2+…+bm-221+bm-120

从二进制表示来讲，倒位序是由原序号的二进制表示按位进行次序反转得到。
不难验证，图3.16输入端的排列次序是3位二进制数的倒位序。例如4=1002，其倒位序为0012=1，这正是图中x［4］的位置。另一例子5=1012，其倒位序取值不变。
不难理解为什么按时间抽取FFT算法的输入顺序是倒位序排列。在序号的二进制表示中，最低位决定奇偶性，最高位决定前一半和后一半。在按时间抽取的分解时，用最低位(序号的奇偶性)决定了被分到上一半还是下一半，因此实际上由最低位和最高位进行了交换； 在下一次分解时，原最高位和最低位都不再起作用，次高位和次低位交换； 当分解过程一直进行下去，到最后一级时，二进制位完成了倒位序过程。
注意到，在分解过程中，FFT的输出顺序保持原来的顺序，称为正位序。
在用数字系统实现FFT时，有一些专用处理器内部带有进行倒位序运算的单元，方便FFT的编程实现。在用通用计算机编程实现FFT运算时，可以通过编写一个专用子程序巧妙地实现倒位序运算。
3.6.2按频率抽取基2FFT算法
与时间抽取FFT的推导不同，换一种思路，将长序列的前后各一半分开，观察会得到什么结果

X［k］=∑N2-1n=0x［n］WnkN+∑N-1n=N2x［n］WnkN

=∑N2-1n=0x［n］WnkN+∑N2-1n=0xn+N2Wn+N2kN(3.46)

=∑N2-1n=0x［n］+xn+N2(-1)kWnkN

在式(3.46)中，仅考虑k=2r的偶数序号DFT系数，得到

X［2r］=∑N2-1n=0x［n］+xn+N2WrnN2=∑N2-1n=0e［n］WrnN2=E［r］(3.47)

仅考虑k=2r+1的奇数序号DFT系数，得到

X［2r+1］=∑N2-1n=0x［n］-xn+N2WnNWrnN2=∑N2-1n=0f［n］WrnN2=F［r］(3.48)

在式(3.47)和式(3.48)中，对偶数序号和奇数序号DFT系数的计算，相当于先定义如式(3.49)和式(3.50)所表示的两个N/2点长的序列，然后做N/2点DFT，即

e［n］=x［n］+xn+N2，0≤n<N/2(3.49)

f［n］=x［n］-xn+N2WnN，0≤n<N/2(3.50)

以N=8为例，这个分解过程如图3.18所示。


图3.18按频率抽取FFT的第一次分解流图


与时间抽取情况类似，第一步分解可节省约一半运算量。可继续分解，第二级分解为

a［n］=e［n］+en+N4，0≤n<N/4

b［n］=e［n］-en+N4WnN/2=e［n］-en+N4W2nN，0≤n<N/4

c［n］=f［n］+fn+N4，0≤n<N/4

d［n］=f［n］-fn+N4WnN/2=f［n］-fn+N4W2nN，0≤n<N/4

显然

X［4r］=E［2r］=∑N4-1n=0a［n］WrnN4=A［k］

X［4r+2］=E［2r+1］=∑N4-1n=0b［n］WrnN4=B［k］

X［4r+1］=F［2r］=∑N4-1n=0c［n］WrnN4=C［k］

X［4r+3］=F［2r+1］=∑N4-1n=0d［n］WrnN4=D［k］

以N=8为例，第二级分解过程如图3.19所示。


图3.19按频率抽取FFT的第二次分解流图


这种分解过程一直持续下去，直到只剩下2点长序列，2点长序列直接实现即可。以N=8为例，第二级分解后只剩下2点序列，直接用2点DFT的蝶形实现，完整的计算过程如图3.20所示。


图3.208点序列按频率抽取FFT的流图



由于这种分解算法，在每一步都将DFT系数(频域序号)分成偶序号和奇数序号分别运算，因此称为按频率抽取的基2FFT算法。它的基本运算蝶形图如图3.21所示，显然其乘法和加法运算量与按时间抽取基2FFT算法一致。


图3.21按频率抽取FFT的蝶形计算结构


注意到，按频率抽取FFT流程中，输入是按正位序排列的，输出是倒位序。
注意到，不管是时间抽取或频率抽取，基2FFT的序号n和k互为倒序。还要注意到，不管是时间抽取还是频率抽取，蝶形运算结构都是从前一级的两个位置取出数据，经蝶形运算后，结果放置回相同的位置。这种运算结构称为同址运算，可有效节省运算过程中的存储器空间。
3.6.3基4和分裂基FFT
设N=4m，类似于基2分解的方式讨论基4FFT算法。设N=4N1，将序列分成4个子序列，这里只讨论按时间抽取的基4FFT，故按照序号对4取余，分成4个子序列。

n=4m0+n0(3.51)

这里0≤m0<N1,0≤n0<4，n0表示每个子序列，m0表示子序列中的序号。
输出序号k按下式分解

k=N1km-1+m-1(3.52)

这个式子是按时间抽取基2FFT输出合成的推广，4个子序列的N1点DFT系数，用序号m-1表示，且0≤m-1<N1，4个子序列DFT系数中，序号m-1取值相等的一组(4个)值构成一个4蝶形，4蝶形的输出序号用km-1，0≤km-1<4表示，4个输出值在最终DFT系数列中间隔为N1，故DFT输出序号表示为式(3.52)。
用稍紧凑一点的表示方法，描述第一次分解过程为

X［k］=∑3n0=0∑N1-1m0=0x［4m0+n0］W(4m0+n0)(N1km-1+m-1)N

=∑3n0=0∑N1-1m0=0x［4m0+n0］Wm0m-1N1Wn0m-1NWn0km-14(3.53)

=∑3n0=0Xm-1,n0Wn0m-1NWn0km-14

上式第二行中，小括号内是每个子序列的DFT，即

Xm-1,n0=∑N1-1m0=0x［4m0+n0］Wm0m-1N1，0≤m-1<N1

针对子序列标号不同，乘以不同的旋转因子Wn0m-1N，然后通过4蝶形运算合成N点DFT。
图3.22是以N=64为例的按时间抽取基4FFT的第一次分解的部分流程图。由于m=3，故图中m-1=2，4个16点DFT分别对应n0的4个取值，每个16点DFT输出序号用2表示。由于复杂性原因，在一张图上已经无法清晰地画出所有蝶形运算，只画出了一个完整4蝶形的示意图，图3.23单独画出了一个标准4蝶形运算。
与基2一样，式(3.53)中的每一个N1点DFT可继续分解，直到分解为4点为止。通过对每一级分解和每一个4蝶形运算量的统计，可以得到基4FFT的乘法运算量为

mc=3×N4×(m-1)≈34Nm=34Nlog4N

=34N12log2N=34N2log2N

可见，基4FFT比基2FFT进一步减少乘法运算次数。



图3.22N=64 的基4FFT部分流程图




图3.234蝶形运算

流程图



直接实现4蝶形需要12次加法，故加法运算量为

ac=12×N4×(m-1)≈3Nm=3Nlog4N=32Nlog2N

可见，基4FFT的加法运算量高于基2FFT。若对每一个基4蝶形重新安排一下运算(分解成2蝶形的级联，如同4点序列的基2FFT实现，见习题)则只需8次加法，这样基4FFT的加法运算与基2FFT的相等。
分裂基分裂基是基2和基4两种FFT的结合。尽管1965年就提出了基2FFT，基4FFT是基2的直接推广，但分裂基却晚了近20年。直到1984年由杜阿梅尔(P Dohamel)和霍尔曼(H Hollmann)提出。
在按时间抽取基2FFT的算法流程中，偶数序号构成的子序列产生的DFT系数没有加旋转因子，奇数序号子序列的DFT系数都加了旋转因子。由于基4FFT进一步减少乘法运算，以此为启发，将奇数序号进一步再分成两个子序列，这样就构成了3个子序列，即

x1［r］=x［2r］,0≤r≤N2-1
x2［l］=x［4l+1］,0≤l≤N4-1
x3［l］=x［4l+3］,0≤l≤N4-1

将原序列的DFT分解成三个子序列的DFT的组合，进行类似基2或基4的推导得到

X［k］=X1((k))N/2+WkNX2((k))N/4+W3kNX2((k))N/4,0≤k≤N-1

将上式分成4段，针对0≤k≤N4-1得到

X［k］=X1［k］+WkNX2［k］+W3kNX3［k］

Xk+N4=X1k+N4+Wk+N4NX2k+N4N4+W3k+N4NX3k+N4N4

=X1k+N4-j(WkNX2［k］-W3kNX3［k］)

Xk+N2=X1k+N2N2+Wk+N2NX2k+N2N4

+W3k+N2NX3k+N2N4

=X1［k］-WkNX2［k］-W3kNX3［k］

Xk+3N4=X1k+3N4N2+Wk+3N4NX2k+3N4N4

+W3k+3N4NX3k+3N4N4

=X1k+N4+j(WkNX2［k］-W3kNX3［k］)

稍加整理，重新写为


X［k］=X1［k］+WkNX2［k］+W3kNX3［k］,

Xk+N2=X1［k］-(WkNX2［k］+W3kNX3［k］),

Xk+N4=X1k+N4-j(WkNX2［k］-W3kNX3［k］),

Xk+3N4=X1k+N4+j(WkNX2［k］-W3kNX3［k］),0≤k≤N4-1 (3.54)



图3.24给出分裂基实现的流图。与基4情况一样，图中只完整画出了一个蝶形，分裂基的蝶形结构更一般地表示在图3.25中。
注意到分裂基的特点，它的蝶形是倒L型的，同时存在不同点数的DFT需要进一步分解。分裂基运算量的统计留作习题，乘法和加法运算量的近似结果分别为

mc=13Nlog2N
ac=Nlog2N

分裂基的乘法运算量比基2FFT节省33%，加法运算量相同。在所讨论的基2、基4和分裂基FFT中，分裂基运算量最少，实现结构比基2FFT稍复杂。分裂基要求N=2m而不是基4的N=4m，对序列长度的要求与基2相同，比基4算法灵活。综合来讲，分裂基在算法有效性和复杂性方面取得了好的平衡。



图3.24分裂基实现流图




图3.25分裂基的蝶形结构







3.6.4滑窗FFT算法
如果数据做DFT时是滑动处理的，即每次处理之后，观察窗移动一个采样点。起始做如下序列的FFT

x0,x1,x2,…,xN-1

滑动一步，做如下序列的FFT

x1,x2,…,xN-1,xN



滑动m步，做如下序列的FFT

xm,xm+1,…,xm+N-2,xm+N-1

现假设已计算出第m步的DFT结果Xmk，求下次滑动后的DFT结果。根据DFT的定义，第m步的DFT为

Xmk=∑N-1n=0xm+nWnkN

第m+1步的DFT为

Xm+1k=∑N-1n=0xm+1+nWnkN

=∑N-2n=-1xm+1+nWnkN+xm+NW-kN-xmW-kN

=W-kN∑N-1n′=0xm+n′Wn′kN+xm+NW-kN-xmW-kN

=W-kNXmk+xm+N-xm(3.55)

可见滑动一步之后的DFT可从上一步的DFT结果加上滑进窗内的新数据与滑出窗外的老数据之差，再乘以旋转因子W-kN得到。
滑动DFT算法的优点是计算量比直接做FFT还要小。如果起始完成一次FFT(假设采用基2算法)之后滑动M步，则总的复数乘法计算量为

12Nlog2N+MN

平均每一步滑动DFT的计算量为

N2Mlog2N+N

当M≈N时，平均每步DFT的计算量为ON。
滑动DFT的另一个优点是，滑动过程中可以不必计算所有频点的变换结果，而只计算所关心的频点的DFT结果，这对某些应用而言可以进一步减少计算量。
以上讨论了每次只滑动一步，若滑动M步则由M次一步滑动构成，也可以讨论一次直接滑动M步的情况，其详细推导留作练习。
由于滑动过程实际上是一个积累过程，因此滑动DFT存在误差稳定性的问题。一旦由于某种原因在滑动计算过程中引入某个误差，这个误差将保持在滑动计算结果中，直至下次重新计算一次起始FFT才能将这个误差清除。
*3.6.5组合数FFT算法简述
将FFT的分解思想应用到任意组合数的情况。假定N是一个组合数，可以表达成υ个因子的乘积

N=N1N2…Nυ

例如，N=3×2×5是一个组合数的例子。
为了简化表示，用Na~b表示下标从a到b的连乘，即

Na~b=Na…Nb

这样
N=N1~υ

则任意n=0,1,…,N-1可以表示为

n=n0+n1N1+n2N1~2+…+nυ-1N1~υ-1

其中

n0=0,1,…,N1-1

n1=0,1,…,N2-1



nυ-1=0,1,…,Nυ-1

任意k=0,1,…,N-1可以表示为

k=k0+k1Nυ+k2Nυ~υ-1+…+kυ-1Nυ~2

其中

k0=0,1,…,Nυ-1

k1=0,1,…,Nυ-1-1



kυ-1=0,1,…,N1-1

DFT的变换基可以表示为

WnkN=W(n0+n1N1+n2N1~2+…+nυ-1N1~υ-1)kN

=Wn0kNWn1N1kNWn2N1~2kN…Wnυ-1N1~υ-1·kN

=Wn0kNWn1kN/N1Wn2kN/N1~2…Wnυ-1kN/N1~υ-1

=Wn0kN1~υWn1kN2~υWn2kN3~υ…Wnυ-1kNυ

根据DFT的定义

Xk=∑N-1n=0xn=n0+n1N1+…+nυ-1N1~υ-1WnkN

=∑N1-1n0=0Wn0kN1~υ∑N2-1n1=0Wn1kN2~υ…∑Nυ-1nυ-1=0xnWnυ-1kNυ

可见，DFT可以分υ步逐步计算。第一步令

Gk0n′=∑Nυ-1nυ-1=0xn′+nυ-1N1~υ-1Wnυ-1k0Nυ,n′=0,1,…,N1~υ-1-1

这样的子序列一共有Nυ个，每个序列的长度为N1~υ-1。对每个这样的子序列再进行下一步计算，令

Gk1k0［n′］=∑Nυ-1-1nυ-2=0Gk0［n′+nυ-2N1~υ-2］Wnυ-2(k0+k1Nυ)Nυ~υ-1，n′=0,1,…,N1~υ-2-1

对每一个k0，这样的子序列共有Nυ~υ-1个，每个序列的长度为N1~υ-2。如此重复计算，第m步有

Gkm-1…k1k0n′=∑Nυ-m+1-1nυ-m=0Gkm-2…k1k0n′+nυ-mN1~υ-mWnυ-mkNυ~υ-m+1,n′=0,1,…,N1~υ-m-1

对每个给定的km-2…k1k0，这样的子序列共有Nυ~υ-m+1个，每个序列的长度为N1~υ-m。如此重复计算直至子序列长度为1止，一共要计算υ步。每一步蝶形的复数乘法次数为Nυ~υ-m+1N1~υ-mNυ-m+1=NNυ-m+1，复数加法次数为

N1~υ-mNυ~υ-m+1Nυ-m+1-1=NNυ-m+1-1

故总的计算量为： 复数乘法

∑υm=1NNυ-m+1=N∑υ-1m=0Nυ-r=N∑υm=1Nm

复数加法

∑υm=1NNm-1

计算量比直接计算DFT所需的N2要少很多，且k和n是位(digit)倒序关系(称为广义倒序关系)。
3.6.6快速傅里叶反变换
对于IDFT，导出快速算法的思路与DFT几乎一致。由于反变换和正变换的计算结构如此相似，不再单独导出快速傅里叶反变换(IFFT)算法，而是采用FFT实现IFFT。由IDFT的定义

x［n］=IDFTX［k］=1N∑N-1k=0X［k］W-knN=1N∑N-1k=0X*［k］WknN*

上式中，括号内的运算与DFT运算完全一致，只是需要以X*［k］为输入序列，故以FFT运算单元实现IFFT的过程如下： 
(1) 以X*［k］作为输入，执行FFT运算； 
(2) FFT输出取共轭并除以N，得到反变换x［n］。
3.7CZT算法
DFT和其标准FFT算法都有其局限性。当DFT的输入N点，输出N点，DFT系数均匀分布于z平面单位圆，两个系数的频域间隔为2πN。可见，DFT有如下局限性： 
(1) 当输入采样点数很少时，若希望DFT输出的频率点数较多，则需要补零做大点数FFT，增加了计算量;
(2) 对于窄带信号，希望通带内采样点密，带外疏，或根本不用计算，用DFT计算需计算全部密集采样点，再取所需的一部分； 
(3) 难以准确得到信号的“自然频率”位置，DFT系数代表的是信号频谱在2πNk，而不是任意感兴趣点的取值。
例3.7.1设一个带通信号，采样获得了N=150个点，希望只分析在π4~3π8之间的频谱，希望在这π8的间隔内，获得128点的采样密度，即两个DFT系数对应的频率间隔为π8×128，如果用DFT(通过FFT计算)，需要计算的总点数是L=128×2ππ/8=2048，由于采样点只有150点，因此，需要补上L-N=1898个零，然后做L=2048点FFT，FFT结束后取出相当于频率π4~3π8的128个系数，即X［256］~X［384］的点。
CZT变换相当于是一种广义的DFT，同时也可以进一步利用快速卷积方法给出CZT的快速算法。相比于DFT是对z变换在单位圆上的均匀采样，CZT变换相当于是在z平面一个更一般曲线上的采样值。如图3.26，左侧是DFT在z平面的部分采样点，右图是CZT在z平面的部分采样点，CZT的采样点是位于一条螺旋线上的。


图3.26DFT(左)和CZT(右)取样点比较



对于N点序列，CZT计算M个输出，N和M可以相等，也可以不等。CZT的一般定义为

XCZ(z)z=zk=∑N-1n=0x［n］z-nk,k=0,1,…,M-1(3.56)

其中

zk=AW-k,0≤k≤M-1

更细致地取

zk=(A0ejθ0)(W0e-j0)-k

实际上zk中的A=A0ejθ0代表第一个采样点在z平面的位置，W=W0e-j0表示每一个采样点在z平面的步进关系，0表示两个采样点对应的角度差，W0表示采样曲线的形式。W0＞1，螺线向中心弯曲； W0＜1，螺线远离中心弯曲； 若取W0=1，A0=1，则采样点位于单位圆上，但起始点由θ0确定。
例3.7.2若取θ0=0,0=2πN,W0=A0=1,N=M，显然CZT就退化为DFT； 若取θ0=π4,0=2π2048,W0=A0=1,N=150,M=128，用CZT即可只计算例3.7.1要求的π4~3π8之间的频谱。
如下导出利用卷积计算CZT的算法，由

XCZ(zk)=∑N-1n=0x［n］z-nk=∑N-1n=0x［n］A-nWnk，0≤k≤M-1

利用

nk=12［n2+k2-(k-n)2］

代入上式得到

XCZ(zk)=∑N-1n=0x［n］A-nW12n2W-12(k-n)2W12k2

=W12k2∑N-1n=0(x［n］A-nW12n2)W-12(k-n)2(3.57)

如果定义两个新序列

f［n］=x［n］A-nW12n2,0≤n<N(3.58)

h［n］=W-12n2,-(N-1)≤n≤M-1(3.59)

代入CZT表示式得

XCZ(zk)=W12k2∑N-1n=0f［n］h［k-n］=W12k2f［k］*h［k］(3.60)

此式将CZT表示成两个序列的卷积形式。以卷积方式实现CZT计算的方框图如图3.27所示。


图3.27CZT计算的卷积结构



由于两个有限长序列的卷积可用DFT来计算，故可用FFT实现CZT的卷积计算形式，若取FFT的长度L≥N+M-1，通过周期延拓得到h［n］在0≤n≤L-1的取值，如图3.28所示。


图3.28单位抽样响应的取值范围


用FFT计算CZT的方框图如图3.29所示，算法详细表述为如下步骤： 


图3.29CZT的FFT实现结构


(1) 选择可用FFT计算的点数，取L=2m； 
(2) 计算f［n］=x［n］A-nWn22,0≤n<N

0,N≤n<L； (N次乘法)
(3) 计算DFT{f［n］}=F［k］； 12Llog2L次乘法


(4) 计算

h［n］=W-n2/2,0≤n≤M-1
0,M≤n<L-N
W-(L-n)2/2,L-N+1≤n<L； 

(5) 计算DFT{h［n］}=H［k］，0≤k≤L-1； 
(6) 计算H［k］F［k］； (L次乘法)
(7) 计算y［k］=IDFT{H［k］F［k］}； 12Llog2L次乘法
(8) 取y［k］的L点中的前M点并乘权值： XCZ(zk)=W12k2y［k］,0≤k≤M-1。(M次乘法)
在该算法实现中，由于h［n］是预先已知的量，因此(4)和(5)两步可预先计算好，在实际操作时不需要运算量。其他各步的运算量放在算法表述的括号中。注意到，这里只统计了乘法次数，并假设用的是基2FFT，若采用分裂基等其他FFT算法可相应调整系数。集中起来总乘法运算量为

mcf=N+12Llog2L+L+12Llog2L+M

=L(log2L+1)+M+N(3.61)

不难统计，直接用卷积计算CZT，需要乘法运算量

mcc=NM(3.62)


例3.7.3对例3.7.1的计算任务，用CZT完成，如例3.7.2讨论的，N=150，M=128，用直接卷积实现CZT，需要运算量

mcc=150×128=19200

若使用FFT计算CZT，利用基2FFT算法，故取L=512，得到运算量

mcc=512×(9+1)+150+128=5398

若采用例3.7.1所示的2048点FFT，乘法运算量为

mc=2048×11/2=11264

在这个任务中，利用FFT的CZT算法，运算量最低。
3.8离散余弦变换及其快速算法
除了DFT外，数字信号处理领域还存在多种离散正交变换，其中有些变换也得到了重要应用。例如DCT变换特别适用于信号压缩问题，在当前的数字图像和数字视频压缩标准中被广泛采用。
3.8.1离散余弦变换
离散余弦变换(DCT)广泛应用于信息压缩中，一维DCT变换的基函数集是

ak［n］=1N，k=0,0≤n≤N-1
2Ncosπ(2n+1)k2N，1≤k≤N-1,0≤n≤N-1(3.63)

因此，矢量x=x［0］,x［1］,…,x［N-1］T的一维DCT定义为

XDCT［k］=［k］1N∑N-1n=0x［n］cosπ(2n+1)k2N,0≤k≤N-1(3.64)

这里

［0］=1

［k］=2,1≤k≤N-1

反变换(IDCT)为

x［n］=1N∑N-1k=0［k］XDCT［k］cosπ(2n+1)k2N(3.65)

很容易验证式(3.63)的基函数集满足正交性和完备性，因此式(3.64)和式(3.65)构成正交变换对。DCT的由来还可以从几个方面来理解。DCT可看作KL变换的一种逼近，KL变换可去除变换序列中的相关性，是信源编码中最理想的变换。但KL变换实现复杂，基函数集不定，限制了其实际应用。人们证明，DCT变换是KL变换的一种良好的逼近。正是如此，DCT在图像和视频编码标准中被广泛采用。由于DFT变换的物理概念明确，接下来讨论DCT与DFT的关系，以此帮助理解DCT的物理意义。
1. DCT与DFT的关系


图3.30延拓的对称序列

可以通过序列延拓建立DFT和DCT之间的关系，对于N点序列x［n］，通过对称性延拓得到序列x2［n］为

x2［n］=x［n］,0≤n<N

x［2N-n-1］,N≤n<2N-1(3.66)

序列x2［n］是长度2N的对称性序列。图3.30示出了N=5的序列对称延拓成N=10序列的例子。
对x2［n］做2N点DFT得到

X2［k］=∑2N-1n=0x2［n］Wkn2N

=∑N-1n=0x［n］Wkn2N+∑2N-1n=Nx［2N-n-1］Wkn2N

对第二个求和项做变量替换为n′=2N-n-1，得

X2［k］=∑N-1n=0x［n］Wkn2N+∑0n′=N-1x［n′］Wk(2N-n′-1)2N

=∑N-1n=0x［n］e-jπ2N·2kn+ejπ2N k(2n+2)

=∑N-1n=0x［n］ejπ2N ke-jπ2N k(2n+1)+ejπ2N k(2n+1)

=2ejπ2N k∑N-1n=0x［n］cosπk(2n+1)2N

对比式(3.64)有

X2［k］=2NXDCT［0］,k=0

2Nejπ2N kXDCT［k］,0<k<N(3.67)

或由2N点DFT系数X2［k］表示DCT系数为

XDCT［k］=12NX2［0］,k=0

12Ne-jπ2N kX2［k］,0<k<N

=［k］12Ne-jπ2N kX2［k］,0≤k<N(3.68)

除了比例因子外，N点序列的DCT系数可由对称延拓2N点序列DFT的前N个系数得到。也可以把DCT看成是2N点变换，利用式(3.64)可验证(留作习题)

XDCT［k］=-XDCT［2N-k］,0≤k<2N

并且把式(3.67)推广为

X2［k］=2NXDCT［0］,k=0

2Nejπ2N kXDCT［k］,0<k<N

0,k=N

-2Nejπ2N kXDCT［2N-k］,N≤k<2N(3.69)

因为DFT意味着时域的周期延拓，故2N点对称序列的DFT其实质对应了如图3.31所示的周期延拓序列。这种延拓相当于先做对称延拓，再以2N为周期进行周期延拓，与直接对序列x［n］以N为周期进行周期延拓相比，大大减少了周期延拓后序列在周期边界处的跳跃性。因此，X2［k］随k的衰减要快于直接对x［n］的DFT变换X［k］，由此使得XDCT［k］比X［k］随k的衰减更快。关于这方面的性质，在DCT的能量紧致特性中进一步讨论。


图3.31对称延拓序列的周期延拓



2. DCT的能量紧致性
对于有限长序列的正交变换，变换系数和信号序列是完全等价的，由一种表示可精确得到另外一种表示。但在一些应用中，需要用到变换的一种特殊性质，即能量紧致性。所谓能量紧致性指的是用变换系数的子集有效重构信号序列的能力。
信号紧致性的典型应用是信源的压缩编码，例如语音、图像或视频的压缩编码，用尽可能少的数据尽可能精确地恢复信号是信源编码的目标。目前最常用的编码技术是变换编码。首先对信号进行变换，在变换域只保留部分变换系数，对这部分变换系数进行存储或传输，在播放端或接收端用这些部分系数反变换重构信号。由于只保留部分变换系数，重构的信号必定存在误差。出于对信源质量的要求，对误差大小要有限制。也就是给定一个误差门限，在误差值不大于门限的情况下，保留最少的变换系数。
对于一种变换来讲，给出要求的误差门限，需要保留的系数数目越少，则称该变换的能量紧致性越好。人们发现，对于大多数实际信号，DCT变换的能量紧致性比DFT要好。也正是这个原因，图像和视频编码标准大都采用DCT变换。
设有N点有限长序列x［n］，分别进行N点DFT和DCT，系数分别记为X［k］和XDCT［k］，用Π表示保留的DFT系数的序号集合，用ΠDCT表示保留的DCT系数的序号集合，用Π表示集合中元素的个数，用X^［k］和X^DCT［k］表示只保留部分系数，其他系数置为零，即

X^［k］=X［k］,k∈Π

0,kΠ(3.70)

和

X^DCT［k］=XDCT［k］,k∈ΠDCT

0,kΠDCT(3.71)

利用部分DFT系数重构的信号为

x^［n］=IDFTX^［k］(3.72)

利用部分DCT系数重构的信号为

x^DCT［n］=IDCTX^［k］(3.73)

重构信号的平方误差和分别为

e=∑N-1n=0x［n］-x^［n］2(3.74)

和

eDCT=∑N-1n=0x［n］-x^DCT［n］2(3.75)

对于一般实际信号(实际语音、实际图像或视频数据)DCT比DFT的能量紧致性好是指如下两种情况之一： 若ΠDCT=Π，则eDCT<e; 或取eDCT=e，则有ΠDCT<Π。
与帕塞瓦尔定理类似，可以证明误差的时域和频域等价表示，即如下的两个等式，证明过程留作习题。

e［Π］=∑N-1n=0x［n］-x^［n］2=1N∑N-1k=0X［k］2-∑k∈ΠX［k］2

=1N∑kΠX［k］2(3.76)

和

eDCT［ΠDCT］=∑N-1n=0x［n］-x^DCT［n］2=∑N-1k=0XDCT［k］2-∑k∈ΠDCTXDCT［k］2

=∑kΠDCTXDCT［k］2(3.77)

注意以上两式的区别，DCT的基序列是归一化的，DFT不是归一化的，故需要除以N。
在变换编码中，怎样选择保留系数的集合Π和ΠDCT是有技巧的，进一步的细节讨论超出本书范围，有兴趣的读者可参考有关逼近论或信源编码的文献［30，58］。
例3.8.1对于一个N=32的实指数衰减序列，只保留部分系数，比较DFT和DCT的能量紧致性。
给出一个最简单的保留系数序号的选择，对于DFT，若只保留一个系数，就保留序号k=0的系数(直流分量)，若保留更多系数，因为实序列DFT的系数满足共轭对称性X［k］=X*［N-k］，保留的非零系数必须保持该性质，因此序号k,N-k必须同时保留。例如，若保留5个系数，则保留系数序号集为Π=0,1,2,N-2,N-1，对于DCT系数，问题简单得多，若保留M个系数，只需保留0≤k≤M-1序号的系数，若保留5个DCT系数，则ΠDCT=0,1,2,3,4。
若保留系数个数为Π=ΠDCT=N时，全部系数被保留，反变换重构的信号无失真，当Π=ΠDCT=1时，只保留一个直流系数，重构信号是常数，误差较大。改变保留系数数目，计算重构误差，结果示于图3.32中。注意图中横坐标是被置为零的系数个数，即N-Π，纵坐标是均方误差，即式(3.76)和式(3.77)误差项除以N得到的平均误差。可见，在此例中，DCT的紧致性比DFT好得多。


图3.32DCT和DFT能量紧致性的比较


注意到，能量紧致性不是从数学上严谨地证明的特性，而是一种统计结果。我们说DCT的紧致性比DFT好，意味着对大多数实际信号讲，这是成立的； 但也能找到一些信号，DFT的紧致性可能更好。例如，对于信号Acos2πNmn+Bcos2πNkn，这里m和k都是整数，对这种特殊信号，DFT的紧致性是最好的。
3.8.2离散余弦变换的快速算法
由于DCT获得广泛应用，DCT的快速算法得到深入的研究。这里简单介绍两类DCT的快速算法。第一类是基于FFT的间接算法，第二类是针对DCT直接导出的快速算法。
基于FFT的快速算法，利用DCT与一些特殊序列DFT的关系，直接利用FFT计算DFT系数后，得到DCT系数。显然，第一个算法是利用式(3.66)得到2N点对称序列，用FFT得到2N点DFT系数，再通过式(3.68)得到DCT系数。若做DCT反变换，由N点DCT系数通过式(3.69)得到2N点DFT系数通过IFFT得到2N点序列x2［n］，再取出前N个点。
一种稍简单的方法是定义2N点序列

x1［n］=x［n］,0≤n<N

0,N≤n<2N-1

对x1［n］做2N点DFT得到X1［k］，可以证明(留作习题)

XDCT［k］=［k］1NRee-jπ2N kX1［k］,0≤k<N

这里Re{·}指取实部。
DCT自身存在高效的快速算法，利用DCT的性质可构造专用的快速算法。这方面结果很多，进一步可以参考文献［59］。这里仅简要介绍Loeffler等提出的一个快速算法，它的结构规范，适用于硬件实现。对于8点DCT仅需要11次乘法和29次加法，乘法次数达到理论上的下限。这个算法的流程图如图3.33(a)所示，其中图3.33(a)中的主要运算单元如图3.33(b)所示。注意，在图(b)的运算单元中，用下列等式可以将4次乘法、2次加法运算变成3次乘法和3次加法

y0=ax0+bx1=(b-a)x1+a(x0+x1)


y1=-bx0+ax1=-(a+b)x0+a(x0+x1)(3.78)

由于a、b是常数，a+b、b-a可以预先算好存在寄存器中，不需要增加运算。对于快速反变换，只需要改换输入/输出方向。
在图像压缩标准JPEG或视频压缩标准MPEG中，对8×8的块图像数据进行DCT变换，可通过多次对8点序列的DCT完成，因此8点DCT是实际中最常用的。图3.33的算法流程图很有用。


图3.338点DCT的快速算法



图3.33(续)



*3.9一些其他离散变换简介
除DFT和DCT外，信号处理中还用到若干离散正交变换，这些变换在不同领域得到了应用，本节简述几种其他的离散变换。
3.9.1离散正弦变换
与离散余弦变换对应，存在一种离散正弦变换(DST)，一维DST变换的基函数集是

ak［n］=2N+1sinπ(n+1)(k+1)N,0≤k≤N-1,0≤n≤N-1(3.79)

因此，矢量x=x［0］,x［1］,…,x［N-1］T的DST定义为

XDST［k］=2N+1∑N-1n=0x［n］sinπ(k+1)(n+1)N+1,0≤k≤N-1(3.80)

反变换(IDST)为

x［n］=2N+1∑N-1k=0XDST［k］sinπ(k+1)(n+1)N+1,0≤n≤N-1(3.81)

DST可以通过对序列进行反对称延拓，建立起与DFT的联系。本质上，DST与DCT一样，也是KL变换的一种逼近。
有学者进一步定义了离散余弦变换和离散正弦变换族，包括四类余弦变换和四类正弦变换。其中3.8节讨论的也是最常用的这种离散余弦变换称为第二类余弦变换(DCTⅡ)。关于余弦和正弦变换族的进一步讨论，请参考文献［89］。
3.9.2Hadamard变换
除一个归一化因子外，Hadamard变换的基序列仅由1、-1组成，适宜于数字信号的合成与分解。

设变换序列N=2n，Hadamard变换矩阵Hn写成如下递推关系

Hn=12Hn-1Hn-1
Hn-1-Hn-1

由初始值

H1=1211
1-1

递推可以得到长度为8的变换矩阵为

H3=1811111111
1-11-11-11-1
11-1-111-1-1
1-1-111-1-11
1111-1-1-1-1
1-11-1-11-11
11-1-1-1-111
1-1-11-111-1

3.9.3Haar变换
Haar变换相当于用两个滤波器Hl=12,12和Hh=12,-12，分别对数据序列滤波，并进行2∶1降采样。得到的高频输出部分保持，而低频部分继续这样的分解过程，直到最后只有一个点。设N=2n，分解过程需进行n次，这个过程也等价为一个线性变换。N=8的Haar变换矩阵为

Hr3=1811111111
1111-1-1-1-1
22-2-20000
000022-2-2
2-2000000
002-20000
00002-200
0000002-2

第8章将会看到，Haar变换是一种最简单的小波变换。
3.9.4Slant变换
Slant变换的核函数也可以通过递推公式求出，此处从略，仅给出一个N=22的变换矩阵的例子

S2=121111
3515-15-35
1-1-11
15-3535-15

还有很多通过特别方法构造的正交变换，这些变换各有一些特殊性质，并获得一些应用。正交变换很多［58，59］，目前非正交的变换甚至突破基函数概念的框架变换、词典变换等的研究也很活跃，有兴趣的读者可参考有关文献［30］。
3.10与本章相关的MATLAB函数与实例
3.10.1相关的MATLAB函数简介
1. fft

功能介绍离散傅里叶变换(FFT)。
语法

y=fft(x)

y=fft(x,n)

y=fft(x,n,dim)
输入变量
y=fft(x)利用FFT算法计算矢量x的离散傅里叶变换。x表示原始信号，当x的长度小于n时，在x的尾部补零，构成n点数据； 当x的长度大于n时，x被截断成为n点数据。若x为矩阵，则fft函数作用于x的每一列。dim表示fft函数作用于x的维度。
输出内容
y代表输出信号。
2. ifft
功能介绍离散傅里叶反变换(IFFT)。
语法

y=ifft(x)

y=ifft(x,n)

y=ifft(x,n,dim)
输入变量
y=ifft(x)利用IFFT算法计算矢量x的离散傅里叶反变换。x表示原始变换系数，当x的长度小于n时，在x的尾部补零，构成n点数据； 当x的长度大于n时，x被截断成为n点数据。若x为矩阵，则ifft函数作用于x的每一列。dim表示ifft函数作用于x的维度。
输出内容
y代表输出信号。
3. fftshift
功能介绍对fft的输出重新排列，将零频分量移到频谱的中心。
语法

y=fftshift(x)
输入变量
y=fftshift(x)对x重新排列。x表示原始变换系数，当x为向量时，将x的左右部分交换； 当x为矩阵时，将x的第1,3象限以及2,4象限分别交换。
输出内容
y代表输出信号。
4. czt
功能介绍线性调频z变换。
语法

y=czt(x,m,w,a)
输入变量
y=czt(x,m,w,a)计算x的线性调频z变换。x表示原始信号，m为要分析的频谱点数，w为z平面取样轮廓线上各点之间的比率，a为轮廓线上的复数起点。
输出内容
y代表输出信号。
5. dct
功能介绍离散余弦变换。
语法

y=dct(x)

y=dct(x,n)
输入变量
y=dct(x)计算矢量x的离散余弦变换。x表示原始信号，当x的长度小于n时，在x的尾部补零，构成n点数据； 当x的长度大于n时，x被截断成为n点数据。若x为矩阵，则dct函数作用于x的每一列。
输出内容
y代表输出的变换系数。
6. idct
功能介绍离散余弦反变换。
语法

y=idct(x)

y=idct(x,n)
输入变量
y=idct(x)计算矢量x的离散余弦反变换。x表示原始变换系数，当x的长度小于n时，在x的尾部补零，构成n点数据； 当x的长度大于n时，x被截断成为n点数据。若x为矩阵，则idct函数作用于x的每一列。
输出内容
y代表输出信号。
3.10.2MATLAB例程
例3.10.1已知信号x［n］=［0,1,2,3,4,5,6,7］，计算其离散傅里叶变换X［k］，中心为零频的序列X^［k］，离散傅里叶反变换x^［n］，并画出相应的图形。

x=［0 1 2 3 4 5 6 7］;

y=fft(x);

y1=fftshift(y);

x1=ifft(y);

figure(1);

stem(x,'ko');

xlabel('samples');

ylabel('Amplitude');

figure(2);

stem(abs(y),'ko');

xlabel('samples');

ylabel('Amplitude');

figure(3);

stem(abs(y1),'ko');

xlabel('samples');

ylabel('Amplitude');

figure(4);

stem(x1,'ko');

xlabel('samples');

ylabel('Amplitude');

运行结果如图3.34所示。


图3.34例3.10.1的例程结果


例3.10.2设x［n］是由3个实正弦信号组成，频率分别为7Hz、8Hz、9Hz，采样频率为50Hz，采样点数为256，在6~10Hz频率段求其CZT。


t=0:1/50:255/50;

x=sin(2*pi*7*t)+sin(2*pi*8*t)+sin(2*pi*9*t);

fs=50;

f1=6; 

f2=10; 

m=50;

w=exp(-j*2*pi*(f2-f1)/(m*fs));

a=exp(j*2*pi*f1/fs);

y=czt(x,m,w,a);

fy=(0:(m-1))*(f2-f1)/m+f1;

plot(fy,abs(y),'ko');

title('CZT');

xlabel('Frequency');

ylabel('Amplitude');

运行结果如图3.35所示。



图3.35例3.10.2的例程结果



例3.10.3计算信号x［n］=2n+100cos(2πn/5),n=1,2,…,50的DCT和IDCT。

n=1:50;figure(2);

x=2*n+100*cos(2*pi*n/5);stem(x1,'ko')； 

y=dct(x);xlabel('samples')； 

x1=idct(y);ylabel('Amplitude')； 

figure(1);

stem(abs(y),'ko');

xlabel('samples');

ylabel('Amplitude');

运行结果如图3.36所示。


图3.36例3.10.3的例程结果

3.11本章小结
DFT和DCT等变换都是基于正交基的正交变换。正交变换的一个明显优势是正、反变换计算上的简洁性。在很多应用中，人们首先选择正交基，正交变换是信号表示中的一种最基本的变换方式。随着研究的深入，也出现很多非正交基，甚至基于“非基”序列集的变换。例如利用更广义的“框架”替代基序列集，构成的一般变换，在表示信号的稀疏性和时频局域性等方面显示出优势。本书主要讨论有关正交变换及其应用，属于“框架”变换的离散小波变换将在第8章给出简要介绍，有关信号变换的更广泛的讨论，可参考有关文献，例如Mallat的文献［30］。
习题
3.1计算如下信号的N点DFT（N为偶数）。
（1） x［n］=RN［n］e-j2πNn
（2） x［n］=δ［n-n1］-δ［n-n2］,0<n1,n2<N
（3） x［n］=1,n为偶数，0≤n<N

0,n为奇数，0≤n<N
（4） x［n］=1,0≤n<N/2

0,N/2≤0≤n<N-1

3.2已知信号x［n］=cosπn6，0≤n<12
0，其他。
（1） 求N=12点的DFTX［k］； 
（2） 若取x［n］的前24点，求N=24点的DFTX1［k］。

3.3如果x～［n］是周期为N的周期序列，那么它也是周期为2N的周期序列。如果把x～［n］看作周期为N的周期序列，令X～1［k］表示x～［n］的离散傅里叶级数的系数； 当把x～［n］看作周期为2N的周期序列，令对应的傅里叶级数的系数为X～2［k］。请利用X～1［k］来确定X～2［k］。

3.4若X［k］=DFT{x［n］}，试证明下列各式。
（1） DFT{x［n］}=NxN-k
（2） X［0］=∑N-1n=0x［n］
（3） x［0］=1N∑N-1k=0X［k］
（4） 若N为偶数，则XN2=∑N-1n=0(-1)nx［n］

3.5已知连续时间信号xa(t)=ej2πt+φ0为复单频信号，现以T为周期对它进行采样得到离散序列x［n］=ej2πnT+φ0。对该序列做N点DFT得到N个变换系数。如果要使得DFT的N个变换系数只有一个不为零而其他全部为零，T和N应满足什么条件？为什么？

3.6已知N点序列x［n］的N点DFT是X［k］。设y［n］定义如下： 

y［n］=x［n］，0≤n<N

0，N≤n<3N


对y［n］做3N点的DFT，请用X［k］表示Y［k］。

3.7已知N点序列x［n］的N点DFT是X［k］。对如下定义的N点序列求N点DFT，请用X［k］表示新的变换系数： 
（1） y1［n］=x［N-1-n］,0≤n<N
（2） y2［n］=(-1)nx［n］,0≤n<N

3.8已知N点序列x［n］的N点DFT是X［k］。对如下定义的2N点序列求2N点DFT，请用X［k］表示新的变换系数： 
(1) y1［n］=x［n］,0≤n<N

x［n-N］,N≤n<2N

0,其他
(2) y2［n］=x［n］,0≤n<N

0,N≤n<2N

0,其他
(3) y3［n］=x［n/2］,n为偶数,0≤n<2N

0,n为奇数,0≤n<2N

3.9一个长度为N的有限长序列，从0时刻起有连续N个非零值，该序列的z变换为X(z)。对该序列补2N个零，做3N点长度的DFT，用X(·)表示第k=5的DFT系数值。
3.10已知N点序列x［n］的N点DFT是X［k］，请证明： 
（1） 如果x［n］是实周期偶对称，即x［n］=xN-n，则X［k］也具有实偶对称性。
（2） 如果x［n］是实周期奇对称，即x［n］=-xN-n，则X［k］为纯虚数且奇对称。

3.11已知N点序列x［n］的N点DFT是X［k］。可通过如下三种方式将其扩展成rN点的序列。
（1） 在原序列后面补零，即

y1［n］=x［n］,0≤n<N

0,N≤n<rN


（2） 在原序列两点间插入r-1个零，即

y2［n］=x［n/r］,n=mr,0≤m<N

0,n≠mr,0≤m<N


（3） 将原序列进行r次的周期延拓，即

y3［n］=x［((n))N］RrN［n］


求上述三个rN点序列的rN点DFTY1［k］、Y2［k］和Y3［k］与X［k］之间的关系。

3.12设X［k］是一个8点长实序列x［n］的8点DFT，已知其中5个值
X［0］=X［1］=0,X［2］=4,X［3］=X［4］=0，求∑7n=0x2［n］。

3.13一个4点序列的DFT系数分别为： {2,1+j,0,1-j}，不加运算能否判断原序列是实序列还是复序列，为什么？序列的零序号值x［0］为多少？
3.14已知长度为9的序列x［n］={-2,3,1,4,-3,0,2,-1,6}，其DTFT为X(ejω)。
（1） 从ω=0开始在［0,2π)区间以间隔π/6对X(ejω)进行均匀采样，得到序列X1［k］。x1［n］为X1［k］对应的IDFT序列，不计算X(ejω)和X1［k］请直接写出序列x1［n］。
（2） 从ω=0开始在［0,2π)区间以间隔π/4对X(ejω)进行均匀采样，得到序列X2［k］。x2［n］为X2［k］对应的IDFT序列，不计算X(ejω)和X2［k］请直接写出序列x2［n］。

3.15已知无限长序列x［n］的z变换为X(z)，且其收敛域包含单位圆； 另一长度为N的序列y［n］的N点DFT为Y［k］。已知Y［k］和X(z)间满足如下关系

Y［k］=X(z)z=ej2πNk,0≤k<N


求y［n］和x［n］间的关系。

3.16已知5点序列x［n］=2δ［n］+δ［n-1］-δ［n-4］。
（1） 求对应的DFT X［k］； 
（2） 分别求x1［n］=IDFT{X2［k］}和x2［n］=IDFT{|X［k］|2}，并说明这两个结果的物理意义。

3.17已知长度为12的实序列x［n］的12点DFT为X［k］，X［k］的前7个样本取值为

X［k］={12,-18-21j,-10+4j,-6+7j,9+8j,19-16j,39},0≤k≤6。不计算X［k］的IDFT，其如下表达式的值。
（1） x［0］（2） x［6］（3） ∑11n=0x［n］（4） ∑11n=0|x［n］|2 （5） ∑11n=0ejπn/3x［n］

3.18已知12点序列x［n］={2,-3,4,-4,5,-3,3,0,-1,-1,-1,1}，其离散时间傅里叶变换为X(ejω)。要求通过M点DFT来计算Xej83π。
（1） 确定M可取的最小值； 
（2） 按照最小M取值给出通过DFT计算Xej83π的过程和结果。

3.19已知如题3.19图所示的5点长序列，画出下列各序列的图形


题3.19图


（1） x～1［n］=x［((-n))5］
（2） x～2［n］=x［((-n))3］
（3） x～3［n］=x［((-n))7］
（4） x4［n］=x［((n-2))5］R5［n］
（5） x5［n］=x［((2-n))3］R5［n］
（6） x6［n］=x［((n-2))7］R7［n］
（7） xep［n］和xop［n］

3.20已知两个序列分别为

x1［n］=δ［n］+δ［n-1］-δ［n-2］-δ［n-3］
x2［n］=δ［n］-δ［n-2］+δ［n-4］


（1） 求两个序列的线性卷积y1［n］=x1［n］*x2［n］； 
（2） 求两个序列的5点循环卷积y2［n］=x1［n］x2［n］,N=5； 
（3） 求两个序列的8点循环卷积y3［n］=x1［n］x2［n］,N=8。

3.21已知x1［n］是100点长的序列，其非零值范围为0≤n≤99； x2［n］是20点长的序列，其非零值范围为10≤n≤39。两序列做100点的循环卷积得到y［n］。请确定y［n］中哪些n值对应x1［n］和x2［n］线性卷积的结果。

3.22计算一个长度为5000点的序列与一个长度为100点序列的线性卷积，要求利用重叠相加法并通过256点的FFT和IFFT来实现。请回答至少需要多少次FFT和多少次IFFT，并请详细说明理由。

3.23已知x［n］是N（N为偶数）点的有限长序列，它的DFT是X［k］。
（1） 如果由X［k］中所有奇次谐波组合成新的N2点频谱X1［r］ 即X1［r］=X［2r+1］， 0≤r≤N2-1，试推导x1［n］=IDFT{X1［r］}0≤n≤N2-1与原序列x［n］的关系？
（2） 根据以上的推导，问对x［n］经过怎样的预处理，就可以用N2点的复数FFT模块得到X［k］中全部的奇次谐波频谱，以N=8为例画出处理的流程框图。

3.24在按时间抽取基2的128点FFT算法中，第一级中与x［47］组成同一个蝶形的是x［m］，请确定m的取值。

3.25对一个长度为N=2m的序列x［n］做N点DFT。如果只需要计算DFT的奇数序号值，即X［1］,X［3］,…,X［N-1］，能否用一个N2点的DFT处理器进行计算？如果可以，请写出N2点DFT处理器的输入序列的一般表达式。

3.26设有一个N=2m点的复数基2 FFT运算模块，要求利用它一次算得2N点的实序列的DFT。请说明计算方法，并画出相应的处理框图。

3.27序列x［n］的长度为N，其中N为偶数。x［n］的N点DFT可以按如下方式计算： 将x［n］拆成两个N2点序列g1［n］=x［2n］和g2［n］=x［2n+1］，分别计算这两个序列N2点DFT得到G1［k］和G2［k］，然后将其进行组合可得到X［k］。如果在构建g1［n］和g2［n］时出现了错误，使得g1［n］=x［2n+1］、g2［n］=x［2n］，但仍按原来方式进行计算，这样就会得到了一个错误的结果X⌒［k］。请利用X［k］给出X⌒［k］的表达式。


3.28设x［n］是长度为N的序列，其中N是偶数，且

x［n］=-xn+N2,n=0,1,…,N2-1


（1） 证明： 在x［n］的N点DFT X［k］中，当k为偶数时，X［k］=0； 
（2） 能否用如题3.28图所示的一个N2点的DFT处理器，通过级联系统A和系统B进行计算？图中y［n］和Y［k］长度均为N2点，Y［k］是y［n］的N2点DFT。如果可以，请推导出x［n］与y［n］的关系，以及X［k］与Y［k］的关系。


题3.28图

3.29一个序列长为N=16，按频率抽取进行基4分解计算其FFT。
（1） 请写出第一级分解的原理表达式； 
（2） 画出第一级自上而下的第2个蝶形和第二级分解的自上而下第3个蝶形的蝶形图（注意，蝶形图画标准的基4蝶形，不需要将基4蝶形分解成基2蝶形，标注蝶形的输入、输出及各支路系数）。

3.30对于长度N=64的序列，进行基4按时间抽取FFT，按从上到下，写出第二级分解对应的前8个旋转因子。

3.31一个长度为9的序列，按3×3的时间抽取进行FFT，推导它的按时间抽取方式分解的原理表达式，画出完整的蝶形图，分析其乘法运算量（其中，乘以±1,±j不计作一次乘法），比较它与直接补零用16点分裂基的运算量（使用分裂基算法时，不考虑一些输入信号是否是由补零产生）。
3.32在一些实际应用中，一边采集信号一边进行DFT运算，当信号采集到N个样本后，得到信号序列

x(1)={x［0］,x［1］,…,x［N-2］,x［N-1］}


用FFT进行计算，得到这第一组DFT系数X(1)［k］。若信号采集继续进行，信号存于固定长度的先入先出的存储器中，在继续采集了M个样本后，存储器里存的信号是

x(2)={x［M］,x［M+1］,…,x［N+M-2］,x［N+M-1］}


对这一组信号的DFT记为X(2)［k］，这样的过程，称为滑动DFT或离散短时傅里叶变换。常用两类方法计算滑动DFT，一种是对每一组信号都直接采用FFT，一种计算方法是通过已计算出的X(1)［k］计算X(2)［k］，
（1） 导出由X(1)［k］和一些采样值联合计算X(2)［k］的一般表达式； 
（2） 假设M=2m,M<N=2L，对如上导出的各组运算尽可能使用FFT，估计其运算量，并与直接进行FFT的算法进行比较。

3.33证明DCT系数满足对称关系： XDCT［k］=-XDCT［2N-k］,0≤k<2N。

3.34定义2N序列点序列

x1［n］=x［n］,0≤n<N

0,N≤n<2N-1


对x1［n］做2N点DFT得到X1［k］，证明DCT变换可由下式求得

XDCT［k］=［k］1NRee-jπ2NkX1［k］,0≤k<N


这里Re{·}指取实部。

3.35证明式（3.76）和式（3.77）。

MATLAB习题
3.1请用MATLAB画出习题3.2中序列x［n］的DTFT，然后利用MATLAB自带fft函数计算序列x［n］的12点和24点的DFT，并将计算结果与DTFT画在同一张图中。通过该练习理解“DFT作为对DTFT的频域离散采样”的含义。
3.2请利用MATLAB自带fft函数计算习题3.20中两个序列的线性卷积和循环卷积，并与按照定义直接计算的结果进行比对。
3.3请利用MATLAB编写任意2的整数次幂点数的基2 FFT的通用函数，并与MATLAB自带fft函数进行比对确认自己编写函数的正确性。
3.4请利用MATLAB按照3.7节中介绍的用FFT计算CZT的算法编写计算CZT的通用函数，并与MATLAB自带的CZT函数进行比对确认自己编写函数的正确性。
3.5请利用MATLAB验证3.4节所讨论的DFT各种性质。