《數(shù)字信號(hào)處理教程-MATLAB釋義與實(shí)現(xiàn)》第四章

第四章

信號(hào)頻譜的高效計(jì)算1信號(hào)頻譜的高效計(jì)算4.1各種傅立葉變換及其關(guān)系4.2快速付立葉變換(FFT)4.3用FFT計(jì)算序列的頻譜4.4時(shí)域采樣中的頻譜變換4.5連續(xù)信號(hào)的頻譜計(jì)算4.6用反變換從頻譜計(jì)算信號(hào)4.7用FFT計(jì)算能量4.8小結(jié)2（1）時(shí)域的周期性對(duì)應(yīng)于頻域的離散化。（2）

時(shí)域的離散化對(duì)應(yīng)于頻域的周期性。其主周期就是乃奎斯特頻率范圍[-π/T,π/T]。（3）

周期性的離散序列將對(duì)應(yīng)于離散并周期性的頻譜。即離散傅立葉級(jí)數(shù)(DFS)。離散有利于數(shù)值計(jì)算，但信號(hào)無(wú)窮延伸又不利于計(jì)算。（4）把時(shí)域和頻域數(shù)據(jù)長(zhǎng)度都限于主周期，并且使之相等，形成離散傅立葉變換(DFT)。它既離散，又長(zhǎng)度有限，適合于計(jì)算機(jī)數(shù)值計(jì)算。

4.1各種傅立葉變換及其關(guān)系3各種傅立葉變換的特點(diǎn)變換名稱時(shí)域信號(hào)（傅立葉反變換）頻譜曲線（傅立葉變換）(連續(xù))傅立葉變換(CFT)

連續(xù)信號(hào)連續(xù)頻譜(連續(xù))傅立葉級(jí)數(shù)(CFS)周期性，連續(xù)信號(hào)離散頻譜離散時(shí)間傅立葉變換(DTFT)離散信號(hào)周期性，連續(xù)頻譜離散傅立葉級(jí)數(shù)(DFS)周期性，離散信號(hào)周期性，離散頻譜離散傅立葉變換(DFT)有限長(zhǎng)離散信號(hào)(隱含周期)有限長(zhǎng)離散頻譜(隱含周期)4對(duì)離散傅立葉變換(DFT)，人們開(kāi)發(fā)了可以高效地進(jìn)行計(jì)算的方法，稱為快速傅立葉變換(FFT)。人們想盡量利用FFT來(lái)解決其他類型信號(hào)的頻譜計(jì)算問(wèn)題。所以要充分弄清各種傅立葉變換之間的關(guān)系。本章就討論這個(gè)主題。先把主要結(jié)果列出，其中有些結(jié)論已經(jīng)討論過(guò)，與采樣定理有關(guān)的結(jié)論將在后面討論，讀者可先接受下來(lái)。目的是走通下圖的路線，完成時(shí)頻域傅立葉變換的數(shù)值計(jì)算。各種傅立葉變換及其相互關(guān)系5各種傅立葉變換及其相互關(guān)系

模擬信號(hào)時(shí)域采樣周期延拓主值區(qū)間時(shí)域xa(t)——x(n)————x(n)·RN(n)|

|(FFT)頻域Xa(jΩ)——X(jΩ)————X(k)·RN(k) CTFTDTFTDFS DFT ICTFT IDTFT 頻域采樣

IDFT6各種傅立葉變換及其相互關(guān)系(1).DFT與DFS的關(guān)系： 或 (4.1.7)(2).DFT與DTFT的關(guān)系

:采樣插補(bǔ)

7各種傅立葉變換及其相互關(guān)系(3).DTFT與CTFT的關(guān)系:用采樣定理或乃奎斯特定理建立關(guān)系。由CTFT求DTFT由DTFT求CTFT:頻率泄漏可以忽略不計(jì)時(shí)近似解。

8各種傅立葉變換及其相互關(guān)系DFT與CTFT的相互關(guān)系:由CTFT求DFT只是取出其一部分，沒(méi)有誤差問(wèn)題。而由DFT求CTFT要先經(jīng)過(guò)由DFT求DTFT，再經(jīng)過(guò)由DTFT求CTFT的兩步。由離散求連續(xù)不一定有解，只有在信號(hào)的頻譜足夠窄，頻率泄漏可以忽略不計(jì)的條件才可能近似由下式求得。

94.2快速付立葉變換(FFT)計(jì)算DFT的運(yùn)算次數(shù)按N2快速增長(zhǎng)。設(shè)N可以被2整除，把x(n)分成兩個(gè)子序列x1(n)和x2(n)，則原序列的DFT可寫(xiě)成(設(shè)N1=N/2)：10快速付立葉變換(FFT)設(shè)它們的傅立葉變換分別為X1(m)和X2(m)，其周期為N1=N/2：則x(n)的傅立葉變換X(m)可表為：如果N=8,則N1=4,故X(m)周期是8,而X1(m)和X2(m)周期是4,即X1(4)=X1(0),X1(5)=X1(1),…。依次取m=0,1,…7，上式對(duì)應(yīng)于右方的運(yùn)算圖。11快速付立葉變換(FFT)12快速付立葉變換(FFT)X1(m)和X2(m)是N1點(diǎn)的DFT，它們的計(jì)算又可以用類似方法化為兩個(gè)更短的N2=N1/2點(diǎn)的DFT,一直分解下去，直到2為止，這就構(gòu)成了上述FFT的全部運(yùn)算圖。粗算其中每一根線條代表一次乘法，線條的合成點(diǎn)代表一次加法。則每一級(jí)要N次乘法和加法。N=8時(shí)，需log28=3級(jí)，故共要24次乘法和加法。原來(lái)要N2=64次。若N=1024，需要10242次乘法，而用FFT,需分解為log21024=9級(jí)，只需1024×9次乘法，加快了100多倍。13快速付立葉變換(FFT)把上述運(yùn)算次數(shù)的估計(jì)精確化：(1).每次分解的乘法次數(shù)為N，共log2N次， 乘法次數(shù)＝N·log2N(2).考慮到蝶形運(yùn)算，又把乘法次數(shù)減半。在FFT運(yùn)算圖中，基本單元為如右圖所示的蝶形結(jié)構(gòu)，實(shí)際上不需要四次乘法，而只需要兩次乘法即可完成。14蝶形運(yùn)算節(jié)省一半乘法考慮到在算時(shí)，把m和N/2+m兩點(diǎn)成組來(lái)進(jìn)行，即構(gòu)成上圖的蝶形運(yùn)算，就節(jié)省一半乘法15快速付立葉變換(FFT)用這些措施后總的乘法次數(shù)約為（當(dāng)N很大時(shí)）：當(dāng)N=1024時(shí)，結(jié)果為5120，與10242＝1048576相比，減小了近200倍。FFT除了提高速度，還要減少計(jì)算時(shí)所用的內(nèi)存。理想的方法是實(shí)現(xiàn)原位計(jì)算，即每次運(yùn)算的結(jié)果就放在輸入數(shù)據(jù)的位置上，最后輸出結(jié)果就放在原輸入數(shù)據(jù)的位置上。這樣所需的內(nèi)存數(shù)目就是N個(gè)。為此，在具體實(shí)現(xiàn)FFT時(shí)，要遵循一些規(guī)則和技巧：16FFT如何節(jié)省內(nèi)存(排序）(1).輸入數(shù)據(jù)和輸出數(shù)據(jù)的下標(biāo)按二進(jìn)制排序：時(shí)域抽取法(Decimate-in-time—DIT)—輸入數(shù)據(jù)遵循倒序排列，則輸出按順序排列。如N=8,二進(jìn)制順序?yàn)閇000,001,010,011,100,101,110,111]：即：0,1,2,3,4,5,6,7二進(jìn)制倒序?yàn)閇000,100,010,110,001,101,011,111]： 即：0,4,2,6,1,5,3,7故DIT輸入數(shù)據(jù)按x(0),x(4),x(2),x(6),x(5),x(3),x(3),x(7)排列，輸出下標(biāo)即為順序排列。17快速付立葉變換(FFT)(2).運(yùn)算結(jié)構(gòu)圖按蝶形運(yùn)算成對(duì)安排；每一段蝶形運(yùn)算完成后，左邊的一對(duì)數(shù)據(jù)就無(wú)用了，它們的位置就用來(lái)保存右邊的一對(duì)數(shù)據(jù)。(3).運(yùn)算從左到右，按段進(jìn)行。每一段運(yùn)算完成后，左一列數(shù)據(jù)就用不到了，全用右一列數(shù)據(jù)來(lái)取代。(4).各段全部運(yùn)算完畢后，N個(gè)輸入數(shù)據(jù)就換成了N個(gè)輸出數(shù)據(jù)。其下標(biāo)序號(hào)則由倒序變成了順序。18快速付立葉變換(FFT)快速付立葉變換的其他變型和討論：以上討論的方法稱為基2-DIT-FFT方法(1).仔細(xì)分析可以發(fā)現(xiàn)，N不必分解到2(基2)，基4的計(jì)算速度是最快的。(2).DIF-FFT方法與DIT-FFT形成對(duì)偶；(3).N不是2的冪次時(shí)可按素?cái)?shù)分解,構(gòu)成其他FFT算法.(4).FFT算法通常都用匯編語(yǔ)言編寫(xiě)，以進(jìn)一步提高運(yùn)算速度，所以各種軟件系統(tǒng)都有FFT模塊可供調(diào)用。一般不需自己編程。19快速付立葉變換(FFT)用MATLAB模仿FFT算法

functiony=myditfft(x)

m=nextpow2(x);N=2^m;%長(zhǎng)度x對(duì)應(yīng)的2的冪次m

%若x的長(zhǎng)度不是2的冪，補(bǔ)零到2的整數(shù)冪

iflength(x)<N

x=[x,zeros(1,N-length(x))];end

%求1:N數(shù)列的倒序

nxb=dec2bin([1:N]-1,m);%求二進(jìn)制順序值

nxd=bin2dec(fliplr(nxd))+1;求十進(jìn)制倒序值

y=x(nxd);%將x倒序排列作為y的初始值

（接下頁(yè)）20快速付立葉變換(FFT)formm=1:m%將DFT作m次基2分解,

對(duì)每次作DFT運(yùn)算

Nmr=2^mm;u=1;%旋轉(zhuǎn)因子u初始化

WN=exp(-i*2*pi/Nmr);%基本DFT因子

forj=1:Nmr/2%本次間隔內(nèi)的各次蝶形運(yùn)算

fork=j:Nmr:N%本次蝶形運(yùn)算跨越間隔 Nmr=2^mm;kp=k+Nmr/2;%蝶形運(yùn)算的下標(biāo)

t=y(kp)*u;%蝶形運(yùn)算的乘積項(xiàng)

y(kp)=y(k)-t;%蝶形運(yùn)算的加法項(xiàng)

y(k)=y(k)+t;%蝶形運(yùn)算的減法項(xiàng)

end

u=u*WN;%修改旋轉(zhuǎn)因子,

end,end21快速付立葉反變換(IFFT)快速付立葉反變換(IFFT):它和正變換是互成對(duì)偶的并有相似的形式:

兩端取共軛：

可見(jiàn)x*(n)是X*(k)/N的傅立葉變換，故要求X(k)的反變換x(n),可先求X*(k)/N的FFT,再把求得的結(jié)果取共軛。22用MATLAB計(jì)算FFT和IFFTMATLAB中FFT函數(shù)的調(diào)用：正變換：X=fft(x);或 X=fft(x,N)N指定了FFT的長(zhǎng)度。在默認(rèn)條件下，它取x的長(zhǎng)度。因?yàn)楫?dāng)N取2的冪次時(shí)，計(jì)算的速度最快，所以當(dāng)x的長(zhǎng)度不是2的冪次時(shí)，應(yīng)盡量指定N。反變換用的是ifft函數(shù)，其調(diào)用方法與fft函數(shù)類似x=ifft(X);或 x=ifft(x,N)23數(shù)字信號(hào)處理器(DSP)中的FFTFFT進(jìn)行的是重復(fù)性的乘積累加計(jì)算（MAC—Multiply-Accumulate-Calculation），普通PC機(jī)用的通用處理器需要用許多個(gè)指令周期才能完成一個(gè)MAC。所以開(kāi)發(fā)了專門(mén)的DSP處理器，其主要特點(diǎn)是能高效地實(shí)現(xiàn)MAC運(yùn)算。例如最新的DSP處理器可以在一個(gè)指令周期內(nèi)完成一條MAC，并且在硬件中實(shí)現(xiàn)位倒序。在所有的DSP芯片開(kāi)發(fā)系統(tǒng)中，F(xiàn)FT也都是標(biāo)準(zhǔn)的固件。

244.3用FFT計(jì)算序列的頻譜

有限長(zhǎng)序列的頻譜計(jì)算：設(shè)已知位于n1~n2區(qū)間的x(n)，要求其頻譜。在使用FFT時(shí)，必須用主值區(qū)間n=0,…,N-1的數(shù)據(jù)，F(xiàn)FT函數(shù)會(huì)產(chǎn)生N個(gè)數(shù)據(jù)，它們應(yīng)該定位在下列頻點(diǎn)上。

如果要求頻譜關(guān)于零頻率對(duì)稱，可以利用X1=fftshift(X)函數(shù)。此時(shí)X1應(yīng)分布在下列頻點(diǎn)上。

25用FFT算序列頻譜例4.3.1考慮長(zhǎng)度為5的有限序列:x=[1,3,5,3,1],要求用FFT來(lái)計(jì)算其頻譜。解：因?yàn)橐蟮氖穷l率域中的連續(xù)頻譜，故求出x的DFT和DTFT,進(jìn)行比較。程序hc431rx=[1,3,5,3,1];nx=0:4;%給定原始數(shù)據(jù)N=length(x);D=2*pi/N;%求出及頻率分辨率k=floor((-(N-1)/2):((N-1)/2));%求對(duì)稱位置向量X=fftshift(fft(x,N));%求對(duì)稱的FFT序列值subplot(1,2,1),plot(k*D,abs(X),'o:')%畫(huà)幅特性subplot(1,2,2),plot(k*D,angle(X),'o:')%相特性圖26用FFT計(jì)算序列的頻譜曲線如下。注意plot語(yǔ)句具有線性內(nèi)插作用，可近似得出連續(xù)頻譜DTFT.27用FFT計(jì)算序列的頻譜用FFT計(jì)算了頻譜時(shí)，其頻率樣本數(shù)等于序列樣本數(shù)，上例中于N=5，所以相鄰兩個(gè)頻率樣本點(diǎn)的間距為這個(gè)頻率間距非常大，分辨率很差。改善分辨率的最好方法是增加數(shù)據(jù)的長(zhǎng)度N。如果沒(méi)有更多數(shù)據(jù)，也可以給輸入序列補(bǔ)零。補(bǔ)了1019個(gè)零以后，得出連續(xù)光滑的頻譜。不過(guò)改善的只是頻譜的視在分辨率。

28用FFT計(jì)算序列的頻譜例4.3.2考慮長(zhǎng)度為11的矩形窗函數(shù)序列，計(jì)算其頻譜。解：假如選N=20作為重復(fù)周期，則要在序列后面補(bǔ)9個(gè)零。使用FFT時(shí)，我們必須按N=20的周期延拓序列中取從n=0到19的主值部分，因此FFT的輸入為

x=[ones(1,6),zeros(1,N-11),ones(1,5)]

在編寫(xiě)程序時(shí)，要準(zhǔn)備給出不同的N進(jìn)行比較。所以程序hc432的寫(xiě)法應(yīng)能適用于不同的N。29序列頻譜的程序hc432rC=[20,1024];forr=[1,2]N=C(r);D=2*pi/N;x=[ones(1,6),zeros(1,N-11),ones(1,5)];%主值區(qū)間的xk=floor(-N/2+0.5:N/2-0.5);%頻率位置向量X=fftshift(fft(x,N));%求出x的對(duì)稱位置FFTsubplot(1,2,r),plot(k*D,X)%畫(huà)實(shí)頻特性圖end30用FFT計(jì)算序列的頻譜由圖，N=20給出了頻譜的大體形狀，但分辨率差，N=1024的結(jié)果則接近于精確的結(jié)果。

31用FFT計(jì)算序列的頻譜長(zhǎng)N的有限序列的頻譜可以用N點(diǎn)FFT方便地求出。存在的問(wèn)題是：頻譜是定義在奈氏頻率范圍內(nèi)所有的ω上的，而FFT只計(jì)算樣本點(diǎn)上的值。因此，要由插值才能求出全部點(diǎn)的頻譜。較簡(jiǎn)單的插值方法是在輸入數(shù)據(jù)的尾部補(bǔ)零，以改善頻率分辨率。如果一個(gè)有限序列位置向量不在主值區(qū)間，比如在某些 ，則這些負(fù)時(shí)間的部分必須移到補(bǔ)零后的尾部，使位置向量位于[0:N-1],然后再用FFT。否則得到的相頻特性將是錯(cuò)誤的。

32用FFT計(jì)算無(wú)限序列的頻譜當(dāng)把一個(gè)無(wú)限序列截?cái)喑蔀橛邢扌蛄校瑫?huì)導(dǎo)致泄漏和波動(dòng)，算出的頻譜與無(wú)限序列的頻譜有誤差。因此要確定截?cái)嗟拈L(zhǎng)度至少應(yīng)有多長(zhǎng)，才能保證其頻譜足夠準(zhǔn)確。

但原無(wú)限序列的頻譜又是未知的，怎么知道誤差？所以就采用另一種思路：不斷地倍增截?cái)嗟拈L(zhǎng)度，把相鄰兩次計(jì)算的結(jié)果進(jìn)行比較。如果其間的誤差已遠(yuǎn)小于算出的頻譜峰值，則認(rèn)為該結(jié)果已經(jīng)接近于（收斂于）精確值。33用FFT計(jì)算無(wú)限序列的頻譜例4.3.3計(jì)算下列序列x(n)的頻譜,精度要求為頻譜峰值的1%。解：這是理想題，它的頻譜已經(jīng)解析地算出為：現(xiàn)在要問(wèn)，應(yīng)該截取多長(zhǎng)的序列，才能使計(jì)算出的頻譜達(dá)到給定的精度。以后再考慮不知道頻譜解析解的情況。34無(wú)限序列頻譜的算例4.3.3rN=32;dw=2*pi/N;n=0:N-1;x=0.7.^n;k=n;X=fft(x); %求出32點(diǎn)數(shù)據(jù)的DFT數(shù)值解%解析法得出的32點(diǎn)頻譜XtXt=1.0./(1-0.7.*exp(-j*k*dw));%求數(shù)值解和解析解的最大誤差e=max(abs(abs(X)-abs(Xt)));Xm=max(abs(X));%求序列的最大幅度pe=e/Xm*100 %求出最大相對(duì)誤差（百分?jǐn)?shù)）

35無(wú)限序列的頻譜計(jì)算程序算出的百分?jǐn)?shù)誤差pe是0.0011,相當(dāng)于十萬(wàn)分之一?？梢?jiàn)取32個(gè)樣本足夠了?，F(xiàn)在設(shè)不知道準(zhǔn)確頻譜。要選定最小的N值，使相對(duì)誤差不超過(guò)β%。這時(shí)，要逐次倍增N值，比較相鄰兩次計(jì)算結(jié)果的誤差，以這個(gè)誤差作為判斷的標(biāo)準(zhǔn)。使得用N=2a個(gè)數(shù)據(jù)點(diǎn)計(jì)算出的頻譜，與用N/2個(gè)數(shù)據(jù)點(diǎn)的計(jì)算結(jié)果的誤差小于峰值的β%，比較頻譜的誤差必須在相同的頻點(diǎn)上進(jìn)行。因?yàn)辄c(diǎn)數(shù)N變了，同樣的頻點(diǎn)在相鄰兩次計(jì)算中的下標(biāo)不同。

36比較N1=N和N2=2N的計(jì)算頻點(diǎn)。用N1點(diǎn)算出的頻譜X1位于下列頻點(diǎn)上

(4.3.4)而用N2點(diǎn)算出的頻譜X2，則位于下列頻點(diǎn)上

(4.3.5)令 ，解出k2=2k1。必須在的這些頻點(diǎn)上比較X1和X2的幅度。按照這個(gè)原理來(lái)編寫(xiě)程序hc424。

無(wú)限序列的頻譜計(jì)算37計(jì)算程序例hc434ara=1;b=100;beta=1;%給定初始數(shù)據(jù)whileb>beta%判斷是否應(yīng)結(jié)束循環(huán)運(yùn)算

N1=2^a;n1=0:N1-1;%確定數(shù)據(jù)長(zhǎng)度N1x1=0.7.^n1;X1=fft(x1);%長(zhǎng)度N1的x1及其FFT1N2=2*N1;n2=0:N2-1;%數(shù)據(jù)長(zhǎng)度加倍

x2=0.7.^n2;X2=fft(x2);%長(zhǎng)度N2的x2及其FFTk1p=0:N1/2-1;k2p=2*k1p;%兩序列對(duì)應(yīng)k2=2k1d=max(abs(X1(k1p+1)-X2(k2p+1)));%對(duì)應(yīng)點(diǎn)的誤差

mm=max(abs(X1(k1p+1)));%X1的最大幅值

b=d/mm*100;a=a+1;%求相對(duì)誤差,長(zhǎng)度加倍end38在此程序中，如果設(shè)β=1，程序?qū)⒌玫絅2=32，b=0.3323，滿足了b<β的要求。這時(shí)的b基本上反映的是N=16時(shí)的誤差，也就是說(shuō)取N=16就夠了。不過(guò)既然已經(jīng)按N=32算出了較精確的結(jié)果，也沒(méi)必要再退回去取N=16重算了。比較相鄰兩個(gè)N值的計(jì)算誤差列表如下：

程序hc434的計(jì)算結(jié)果N1/N24/88/1616/3232/6464/128

誤差b%24.015.76480.33230.00111.22e-8394.4

時(shí)域采樣中的頻譜變換

時(shí)域采樣定理(Nyquist定理)需要弄清經(jīng)采樣后的序列與原始的模擬信號(hào)在時(shí)域和頻譜方面的數(shù)學(xué)關(guān)系。先把采樣的過(guò)程看作模擬信號(hào)通過(guò)一個(gè)電子開(kāi)關(guān)S；再考慮增益而改為乘法器。見(jiàn)圖4.1.1(a)和(b)。把開(kāi)關(guān)信號(hào)等價(jià)成一寬度為τ，周期為T(mén)、面積為1（因而幅度為1/τ）的矩形脈沖串 ，采樣信號(hào)就是xa(t)與 相乘的結(jié)果。如果讓電子開(kāi)關(guān)合上時(shí)間τ→0，則形成理想采樣，此時(shí)上面的脈沖串變成單位沖激串40連續(xù)信號(hào)采樣中的頻譜變換采樣過(guò)程用電子開(kāi)關(guān)模型(a)和乘法器模型(b)來(lái)表示。后者將采樣過(guò)程看作把xa(t)與寬度為τ,周期為T(mén),面積為1，因而幅度為1/τ的矩形脈沖串PT(t)相乘，當(dāng)寬度

τ趨向于零時(shí)，這些脈沖趨向于δ函數(shù)。41連續(xù)信號(hào)采樣中的頻譜變換因此連續(xù)信號(hào)與其采樣序列間的時(shí)域關(guān)系表為：對(duì)等式兩端求傅立葉變換，由于 而即單個(gè)脈沖具有無(wú)限寬的均勻頻譜。當(dāng)單個(gè)脈沖構(gòu)成等間隔脈沖序列時(shí)，42連續(xù)信號(hào)采樣中的頻譜變換它的頻譜是一個(gè)等幅度、周期性的脈沖頻譜。因?yàn)闀r(shí)域相乘對(duì)應(yīng)于頻域卷積：δ函數(shù)的一個(gè)特性是：它與任何函數(shù)的卷積就等于該函數(shù)。一串間隔為Ωs的δ函數(shù)與X(jΩ)卷積的結(jié)果是43連續(xù)信號(hào)采樣中的頻譜變換這是采樣序列頻譜和原模擬信號(hào)頻譜的關(guān)系，注意兩者也差了一個(gè)時(shí)間量綱T，與時(shí)域相同。上式表明，采樣信號(hào)的頻譜是一串原模擬信號(hào)的頻譜無(wú)限次平移的疊加。平移的間隔為采樣角頻率Ωs=2πFs。可以看出，知道Xa(jΩ)求X(jΩ)是有唯一解的，而知道X(jΩ)求Xa(jΩ)則不是唯一的。44連續(xù)信號(hào)采樣中的頻譜變換設(shè)模擬信號(hào)xa(t)是帶限信號(hào)，最高截止頻率為Ωc，其頻譜Xa(jΩ)如

(a)圖。Pδ(t)的頻譜Pδ(jΩ)如圖(b)。是一連串間隔為Ωs的脈沖頻譜串。x(t)的頻譜X(jΩ)是(a)和(b)的卷積。Ωc<Ωs/2時(shí)如圖

(c)，原模擬頻譜互相分開(kāi)。Ωc>Ωs/2時(shí)如圖

(d)，原模擬頻譜互相交疊。45連續(xù)信號(hào)采樣中的頻譜變換設(shè)模擬信號(hào)xa(t)是帶限信號(hào)，最高截止頻率為Ωc，其頻譜如圖4.4.2(a)所示。Pδ(t)的頻譜如圖4.4.2(b)所示。那么按照（4.4.7）式，x(t)的頻譜如圖4.4.2(c)所示。如果基帶頻譜較窄，滿足Ωc<Ωs/2，基帶譜與它的周期延拓的譜沒(méi)有重疊，此時(shí)的序列頻譜就如圖4.4.2(c)所示的周期頻譜。如果讓這個(gè)采樣信號(hào)經(jīng)過(guò)理想低通濾波器G(jΩ)，就能得到輸出y(t)，不失真地提取原模擬信號(hào)，如圖4.4.3所示。此時(shí)，(4.4.7)式成為：46連續(xù)信號(hào)采樣中的頻譜變換對(duì)離散序列進(jìn)行低通濾波[如圖

(a)]復(fù)原連續(xù)信號(hào)的可行性：如果原模擬信號(hào)的頻譜比采樣頻率之半低，即Ωc<Ωs/2，如圖

(b);取低通濾波器的帶寬等于乃奎斯特頻率，如圖

(c);則濾波輸出的頻譜就是原模擬信號(hào)的頻譜如圖

(d);47連續(xù)信號(hào)采樣中的頻譜變換(4.4.8)式說(shuō)明，在沒(méi)有混疊的情況下，序列的頻譜等于原始模擬信號(hào)的頻譜除以采樣周期。如果選擇采樣頻率Ωs低，或者說(shuō)信號(hào)最高截止頻率Ωc高，以至Ωc>Ωs/2，或者fc>Fs/2，則基帶譜與把它作Ωs的周期延拓形成的譜就會(huì)發(fā)生重疊，稱為頻譜混疊現(xiàn)象，如圖4.4.2(d)所示。這種情況下，再用圖4.4.3所示的理想低通濾波器對(duì)它進(jìn)行濾波，得到的將是失真了的模擬信號(hào)。

48連續(xù)信號(hào)采樣中的頻譜變換采樣定理(Nyquist定理)：

采樣信號(hào)頻譜是原連續(xù)信號(hào)的頻譜以采樣頻率Ωs為周期進(jìn)行周期延拓并疊加形成的，用公式（4.4.7）表示。

Ωs/2稱為折疊頻率或乃奎斯特頻率。如信號(hào)帶寬大于折疊頻率，Ωc>Ωs/2，則采樣信號(hào)x(t)的頻譜會(huì)是原模擬信號(hào)頻譜折疊相加，稱為混疊或泄漏，造成頻譜的失真。實(shí)際中需根據(jù)模擬信號(hào)的最高頻率Ωc，按Ωs>2Ωc的要求選擇采樣頻率Ωs(Fs=Ωs/2π

）。49模擬信號(hào)的重構(gòu)重構(gòu)模擬信號(hào)是采樣的逆過(guò)程——從離散到連續(xù)。設(shè)序列x(n)的頻譜為X(jω)，根據(jù)公式(4.4.7)，它與模擬頻譜Xa(jΩ)關(guān)系如下：如果信號(hào)帶寬小于折疊頻率，則可以用帶寬為[-π/T,π/T]的理想濾波器完全無(wú)失真地恢復(fù)原來(lái)的低頻模擬頻譜。在時(shí)域也可完全無(wú)失真地恢復(fù)原來(lái)的低頻模擬信號(hào)波形。50模擬信號(hào)的重構(gòu)由離散序列重構(gòu)模擬信號(hào)的過(guò)程可以分解成兩步。(1).首先經(jīng)過(guò)一個(gè)脈沖變換器，將離散序列乘以T，使它變成一個(gè)等價(jià)的模擬信號(hào)；(2).然后再讓它通過(guò)一個(gè)低通濾波器。濾波器的輸出就是模擬信號(hào)。下面將分理想方法和實(shí)際方法兩種情況討論。

51模擬信號(hào)的重構(gòu)

理想方法：

(1).脈沖變換器把每個(gè)數(shù)字樣本x(n)變成相應(yīng)權(quán)值的δ函數(shù)，即產(chǎn)生如下的時(shí)間函數(shù)xs(t)(2).然后xs(t)讓通過(guò)一個(gè)幅頻特性如圖4.4.3(c)的理想重構(gòu)濾波器。其數(shù)學(xué)形式為：

52模擬信號(hào)的重構(gòu)對(duì)上述理想頻率特性作傅立葉反變換（例3.2.8),得到它的脈沖響應(yīng)：右圖是它的波形，下圖是使它右移一拍。可以看出，它是非因果的，所以是不可實(shí)現(xiàn)的理想濾波器。53模擬信號(hào)的重構(gòu)濾波后的輸出ya(t)應(yīng)為xs(t)和h(t)的卷積：它是理想濾波器的脈沖響應(yīng)平移加權(quán)的疊加，如右圖所示。54模擬信號(hào)的重構(gòu)由于理想脈沖響應(yīng)h(t-nT)在t=nT點(diǎn)都等于零，疊加的結(jié)果對(duì)除中心樣本點(diǎn)外其它樣本點(diǎn)的值沒(méi)有影響。所以輸出的模擬信號(hào)ya(t)將經(jīng)過(guò)所有的樣本點(diǎn)x(n)，它就是樣本序列的包絡(luò)。理想方法的兩個(gè)步驟都是實(shí)際上無(wú)法實(shí)現(xiàn)的。第一步中的脈沖函數(shù)δ(t)在工程中無(wú)法做到；第二步中的理想濾波器的脈沖響應(yīng)又是非因果過(guò)程。所以理想方法只有理論意義，它給出了一種理論計(jì)算的方法和努力追求的目標(biāo)。55模擬信號(hào)的重構(gòu)實(shí)際方法：(1).用零階保持器（zero-order-holder—ZOH）代替脈沖變換器。它把序列保持一個(gè)采樣周期T，直到下一個(gè)樣本的到來(lái)。（圖1.1.1)這樣，就相當(dāng)于把信號(hào)乘了T，并使之成為模擬信號(hào)；

(2).零階保持器自身也有低通濾波作用，因?yàn)樗拿}沖響應(yīng)為：56模擬信號(hào)的重構(gòu)用傅立葉變換求得并畫(huà)出零階保持器的頻率特性：57模擬信號(hào)的重構(gòu)零階保持器的幅頻特性的通帶約為π，即乃奎斯特頻率。這是很好的特性。從時(shí)間波形而言，它大體相當(dāng)于圖2.1.1(b)，但保持波形右移了半拍。所以它基本反映了樣本序列的包絡(luò)，不過(guò)在時(shí)間上延遲了半個(gè)采樣周期。這也可以從它的相頻特性上看出。作為濾波器，零階保持器的過(guò)渡帶太緩，高頻區(qū)又有較大的波動(dòng)。反映在波形中有高頻的尖峰，通常還不能滿足工程的指標(biāo)，往往還要加上后濾波器。58預(yù)濾波作用的定量分析按兩種情況作對(duì)比分析：(1)在圖4.4.7(a)中，連續(xù)信號(hào)xa(t)經(jīng)過(guò)預(yù)濾波后為xp(t)，然后把它采樣后變成x1(nT)。再用一個(gè)帶理想低通濾波器的D/A變換器把它轉(zhuǎn)換回連續(xù)信號(hào)xa1(t)。(2)在圖4.4.6(b)中，連續(xù)信號(hào)xa(t)不經(jīng)過(guò)預(yù)濾波而被采樣，得到的x2(nT)也立即再用同樣的D/A變換器把它轉(zhuǎn)換成連續(xù)信號(hào)xa2(t)。比較誤差xa1(t)-xa(t)和xa2(t)-xa(t)

哪個(gè)大？59預(yù)濾波作用的定量分析60預(yù)濾波作用的定量分析用誤差平方積分E1和E2和帕塞瓦爾定理：由于x1(n)沒(méi)有混疊，在奈奎斯特主頻區(qū)間有經(jīng)過(guò)理想重構(gòu)后，只保留主頻區(qū)間頻譜不變，別處為零61預(yù)濾波作用的定量分析平方積分E1中[-π,π]的一段積分為零：可寫(xiě)成對(duì)于E2,，由于存在著頻率混疊，在|Ω|≤π/T的范圍內(nèi)，Xa2(Ω)將不等于Xa(Ω)，因此平方積分E2可分成三段來(lái)寫(xiě)：62預(yù)濾波作用的定量分析與E1相比，其首尾兩段相同，所以可寫(xiě)成這證實(shí)了信號(hào)經(jīng)過(guò)預(yù)濾波后再作數(shù)字處理的比不經(jīng)過(guò)預(yù)濾波的誤差要小。因?yàn)檫@個(gè)誤差主要取決于采樣后造成的頻率混疊。預(yù)濾波器的帶寬應(yīng)當(dāng)取得比乃奎斯特頻率π/T小。634.5連續(xù)信號(hào)的頻譜計(jì)算對(duì)于實(shí)際中遇到的大多數(shù)連續(xù)信號(hào)，都不能采用解析方法求頻譜。而必須采用計(jì)算機(jī)輔助的方法。第一步是采樣，把連續(xù)信號(hào)離散化；第二步是截?cái)?，取有限的?shù)據(jù)。不然就用不成計(jì)算機(jī)。所以，不僅是實(shí)時(shí)處理信號(hào)要用數(shù)字信號(hào)處理的理論，非實(shí)時(shí)地分析研究大部分連續(xù)信號(hào)照樣要用數(shù)字信號(hào)處理的理論和方法。64設(shè)xa(t)為絕對(duì)可積的連續(xù)時(shí)間信號(hào)，x(nT)是它的采樣序列，如果T足夠小，有數(shù)字頻譜X(jΩT)=X(jω)可以用FFT求得，所以就可由上式求得連續(xù)時(shí)間信號(hào)的頻譜。計(jì)算中的關(guān)鍵問(wèn)題是如何選擇采樣周期T、數(shù)據(jù)長(zhǎng)度N和/或信號(hào)持續(xù)時(shí)間L。三個(gè)變量之間以L=TN相關(guān)聯(lián)，只有兩個(gè)獨(dú)立可選。連續(xù)信號(hào)的頻譜計(jì)算65連續(xù)信號(hào)的頻譜計(jì)算在參數(shù)選擇時(shí)，要考慮如下三個(gè)問(wèn)題：1．頻率混疊：如果不是有限帶寬的，則采樣周期T的選擇的主要依據(jù)是，取足夠小的T，使得由時(shí)間采樣所造成的頻率混疊小到可以忽略。2．頻率分辨率：連續(xù)信號(hào)頻譜Xa(Ω)是定義在所有Ω上的，當(dāng)使用FFT來(lái)計(jì)算時(shí)，只計(jì)算出了它在N個(gè)數(shù)字頻率上的樣本值X(jω)，ω=2kπ/N,k=0,1,…,N-1。數(shù)字頻率分辨率為dω=2π/N，選擇大的N可以改善數(shù)字頻率分辨率，獲得頻譜細(xì)節(jié)的更多信息。663．截?cái)嘈?yīng)：如果信號(hào)是無(wú)限長(zhǎng)的，就必須把它截?cái)嗟介L(zhǎng)度L=NT，截?cái)嘞喈?dāng)于加窗，它會(huì)引起吉布斯效應(yīng)（波動(dòng)），也會(huì)把窗函數(shù)的頻譜引入信號(hào)頻譜，造成混疊，需要考慮其誤差的問(wèn)題。此外，模擬頻率分辨率為D=dω/T=2π/L，它與L成反比，選擇大的信號(hào)持續(xù)時(shí)間L=NT可以改善模擬頻率分辨率，得知頻譜的細(xì)節(jié)。

連續(xù)信號(hào)的頻譜計(jì)算67連續(xù)信號(hào)的頻譜計(jì)算由于上述三個(gè)問(wèn)題，當(dāng)用FFT來(lái)計(jì)算的頻譜時(shí)必須小心。很明顯，采樣周期T應(yīng)該選得小些以減少混疊，而N要選得足夠大來(lái)提高分辨率；如果N是確定的，為了減小頻率混疊就要減小T；為了減小截?cái)嘈?yīng)，需要增加L=NT，所以要增大T；這是兩個(gè)互相矛盾的要求。此外，取很小的T和很大的N往往又是不必要的。因此建議在計(jì)算與頻譜時(shí)，采用以下的步驟來(lái)選擇T和N。68連續(xù)信號(hào)的頻譜計(jì)算1．先初步選擇時(shí)間紀(jì)錄長(zhǎng)度L，使得0到L之間包括了大部分非零的xa(t)，然后用逐次減小T和加大N的步驟來(lái)選擇周期T，使得時(shí)間采樣造成的頻率混疊可以忽略不計(jì)。2．在選定T以后，進(jìn)一步選擇L。即把L(和N)增加一倍。并與用原來(lái)的結(jié)果進(jìn)行比較。因?yàn)閮蓚€(gè)頻譜所用的T相同，因此，兩個(gè)頻譜之間的差別是由于截?cái)嘈?yīng)不同（即L不同）引起的。如果兩者的誤差不大，這個(gè)L就可以接受；反之，再把L加倍，直到用兩個(gè)相繼的L算出的頻譜非常接近為止。69連續(xù)信號(hào)的頻譜計(jì)算例4.5.1考慮連續(xù)時(shí)間信號(hào):用FFT計(jì)算其頻譜。解：先選tmax=L=10, ，所以L=10并沒(méi)有完全包括信號(hào)的主要部分。但時(shí)間區(qū)間[0,10]還是能大體反映t>0的信號(hào)的主要頻率分量。再隨意地選擇N=5，得到T=L/N=2。然后編寫(xiě)計(jì)算頻譜的程序hc451,核心語(yǔ)句為：70連續(xù)信號(hào)頻譜計(jì)算例4.5.1T=？;N=？;D=2*pi/(N*T);n=0:N-1;x=exp(-0.1*n*T); ％求出離散序列X=fftshift(fft(x));％用求對(duì)稱離散頻譜，Xa=T*X； ％乘以T求連續(xù)頻譜k=floor(-(N-1)/2:((N-1)/2);％位置向量也取對(duì)稱Xa(1) %求乃奎斯特頻率處的頻譜值plot(k*D,abs(Xa)) %繪圖

程序中有幾個(gè)地方需加說(shuō)明：71連續(xù)信號(hào)頻譜計(jì)算例4.5.1(1).從模擬信號(hào)求連續(xù)頻譜，經(jīng)過(guò)的步驟是：連續(xù)信號(hào)xa(t)→(采樣)離散序列x(n)=xa(nT)→(FFT)序列頻譜X(k)=fft(x)→(乘T)連續(xù)信號(hào)頻譜Xa(kD)=TX(k)，這反映在2～4行程序中.(2).檢驗(yàn)混疊，用的方法是在折疊頻率上的頻譜幅度足夠小。以判斷參數(shù)T選擇是否合理。故求對(duì)稱頻譜向量的第一項(xiàng)X(1)的值。(3).選擇N,用不同T和N算得到的幅頻特性畫(huà)在圖中。注意六張子圖上的頻率坐標(biāo)是相同的，但它們的乃奎斯特頻率邊界[-π/T,π/T]不同。72連續(xù)信號(hào)的頻譜計(jì)算73連續(xù)信號(hào)的頻譜計(jì)算T的選擇：可以用乃奎斯特頻率處的幅特性來(lái)評(píng)價(jià)混疊的嚴(yán)重程度。圖4.5.1(a)中T=2,乃奎斯特頻率邊界是

[-1.57,1.57];，圖(b)的T=1,邊界是[-3.14,3.14];（超出了圖中橫軸的邊界）?？梢钥闯?，在圖4.5.1(b)的邊界(±3)處，有一個(gè)相當(dāng)大的非零值Xa(1)=0.3319。說(shuō)明T=1時(shí)，頻率混疊仍不能忽略。繼續(xù)多次把T減半和把N加倍（保持L不變）進(jìn)行計(jì)算，直到T=0.1,N=100，這時(shí)的乃奎斯特頻率已經(jīng)是±31.4，其上的Xa(1)值已減小到0.0318，峰值是6，相對(duì)誤差小于0.032/6=0.4%。這說(shuō)明，選T=0.1，頻率混疊可以忽略不計(jì)。

744.5連續(xù)信號(hào)的頻譜計(jì)算N的選擇：圖4.5.1(e)用實(shí)線畫(huà)出了T=0.1,N=200的計(jì)算頻譜，與圖4.5.1(d)中T=0.1,N=100的結(jié)果進(jìn)行比較。這里只畫(huà)出了[-3,3]的頻率范圍，可見(jiàn)它們之間的區(qū)別很大，說(shuō)明截?cái)嘈?yīng)相當(dāng)明顯，有必要采用更大的L。圖4.5.1(f)給出了用T=0.1,N=400的計(jì)算頻譜。與圖4.5.1(e)相比仍有相當(dāng)差別。如果再畫(huà)出T=0.01,N=800的結(jié)果，它們的誤差在[-3，3]頻段上就很難分辨了。因此取L=0.1×800=80秒已經(jīng)大得足以避免截?cái)嘈?yīng)的影響。

75截?cái)嚅L(zhǎng)度和補(bǔ)零的影響Xa是t=0到L=40，采樣周期為0.1秒所得的全部400個(gè)數(shù)據(jù)，即xa=exp(-0.1*(1:N)*T);xa1是只取xa的前N2=100個(gè)數(shù)據(jù),而后面N-N2個(gè)數(shù)據(jù)補(bǔ)以零的數(shù)組xa1=[exp(-0.1*[1:N2]*T),zeros(1,N-N2)];xa和xa1的波形及其頻譜見(jiàn)圖4.5.2。由于截?cái)啵l譜Xa1有波動(dòng)。即窗函數(shù)的吉布斯效應(yīng)。解的結(jié)果為：可以從中看出波動(dòng)分量。76截?cái)嚅L(zhǎng)度和補(bǔ)零的影響77截?cái)嚅L(zhǎng)度和補(bǔ)零的影響可見(jiàn)，序列補(bǔ)零可以得到高的分辨率，經(jīng)FFT得到很密的頻譜，但頻譜的值不見(jiàn)得精確，因?yàn)樗鼪](méi)有增加信息。在實(shí)用中，如果可以得到的更多數(shù)據(jù)時(shí)，應(yīng)該補(bǔ)充數(shù)據(jù)而不是補(bǔ)零。補(bǔ)零可能反而給出了頻譜中不正確的細(xì)節(jié)。在本例中，時(shí)域截?cái)嗍怯?jì)算中人為地引入的，其頻域波動(dòng)也就是人為引入的。所以首先要解決的是原始時(shí)域數(shù)據(jù)改善的問(wèn)題，而不是去補(bǔ)零來(lái)提高分辨率。

78連續(xù)信號(hào)的頻譜計(jì)算為了用一個(gè)程序而，用不同的參數(shù)畫(huà)出圖4.5.1的多幅子圖?？梢杂胒or循環(huán)如下：T0=[2,1,0.5,0.1,0.1,0.1];%各次的T,編成向量L0=[10,10,10,10,20,40];%各次的N,編成向量forr=1:6%循環(huán)計(jì)算六次

T=T0(r);N=L0(r)/T0(r);%順序調(diào)用T及N（下接程序hc451核心語(yǔ)句）subplot(3,2,r),%在第r個(gè)子圖位置plot(k*D,abs(Xa))%繪圖end79連續(xù)信號(hào)的頻譜計(jì)算例4.5.2用FFT計(jì)算下列連續(xù)信號(hào)的頻譜。解：本例的計(jì)算原理同上，因?yàn)樾盘?hào)由兩個(gè)頻率很接近的正弦信號(hào)組成，所以重點(diǎn)是計(jì)算時(shí)保證頻譜的分辨率D,它取決于截?cái)嘈盘?hào)的長(zhǎng)度L,D=2π/L。計(jì)算的程序?yàn)閔c452,80連續(xù)信號(hào)的頻譜計(jì)算T0=[0.6,0.15,0.15,0.15];N0=[256,256,256,2048];forr=1:4(類似于hc451的計(jì)算連續(xù)信號(hào)頻譜核心語(yǔ)句）subplot(2,2,r),plot(k*D,abs(Xa))%畫(huà)幅頻特性end上述四組T0和N0對(duì)應(yīng)于L取[153.6,38.4,307.2]三種數(shù)據(jù).其頻率分辨率為D=[0.0409,0.1636,0.0205]。由于兩個(gè)正弦頻率之差為0.1，所以只有T0=0.15和N=2048能同時(shí)保證較小的混疊和較高的分辨率。8182連續(xù)信號(hào)的頻譜算例例4.5.3計(jì)算下列信號(hào)的頻譜設(shè)a=5，這個(gè)信號(hào)是有限時(shí)間信號(hào)，它是實(shí)的和偶的，因此它的頻譜也一定是實(shí)的和偶的。先任選T=1，則在|t|≤5的h(nT)中總共有11個(gè)樣本點(diǎn)。因此，把時(shí)間樣本數(shù)選成奇數(shù)，并且要按循環(huán)對(duì)稱的概念排列在0到N-1的主值下標(biāo)區(qū)間內(nèi)。編出程序hc45383連續(xù)信號(hào)的頻譜計(jì)算T=1;M=5/T;N=2*M+1;n=1:M;%原始數(shù)據(jù)D=2*pi/(N*T);%頻率分辨率hp=sin(2*n*T)./(pi*n*T);%生成正序列h=[2/pihpfliplr(hp)];%構(gòu)成循環(huán)對(duì)稱序列H=T*fftshift(fft(h));H(1)%求xa的對(duì)稱FFT，k=floor(-(N-1)/2:(N-1)/2);%頻率下標(biāo)向量

subplot(1,2,1),plot(k*D,H)%繪圖

程序運(yùn)行結(jié)果說(shuō)明，N取11時(shí)分辨率太低，取T=0.5,用1024點(diǎn)FFT,(相當(dāng)于后面補(bǔ)零到1024）可得到較好結(jié)果。84連續(xù)信號(hào)的頻譜計(jì)算85連續(xù)信號(hào)的頻譜計(jì)算周期信號(hào)的頻譜計(jì)算：從理論上說(shuō)，如果一個(gè)連續(xù)時(shí)間信號(hào)是周期性的，它必定有無(wú)限的長(zhǎng)度，計(jì)算出的頻譜將不會(huì)收斂。然而，用計(jì)算機(jī)計(jì)算時(shí)，序列總是有限長(zhǎng)的。因?yàn)槲覀冇龅降亩际墙?jīng)過(guò)截?cái)嗔说闹芷谛盘?hào)。從算法上看，完全可以用上節(jié)中的方法和程序，只是在計(jì)算結(jié)果上，頻譜會(huì)有向有限個(gè)頻點(diǎn)上集中，并成為脈沖的趨勢(shì)。

86連續(xù)信號(hào)的頻譜計(jì)算例4.5.4計(jì)算定義在全部t上的xa(t)=cos5t的頻譜，它的理論頻譜是 ,它包含了權(quán)重為π的位于Ω=±5上的兩個(gè)脈沖函數(shù)。

解：在計(jì)算機(jī)計(jì)算中，正余弦函數(shù)必須截?cái)酁橛邢揲L(zhǎng)度L。信號(hào)＝cos5t的帶寬限制于5。純理論地看，只要采樣周期小于π/5=0.63秒，就不會(huì)發(fā)生頻率混疊。然而如果把cos5t截?cái)酁殚L(zhǎng)L的信號(hào)，則它的頻譜就不再是有限帶寬了，所以必須采用更小的采樣周期，任選T＝0.1，并選N＝50，得到L＝TN＝5。按此來(lái)截?cái)嘈盘?hào)。

87連續(xù)信號(hào)的頻譜計(jì)算程序hc454。N=input('N=');T=0.1;n=1:N;%原始數(shù)據(jù)D=2*pi/(N*T);%頻率分辨率xa=sin(5*n*T);%生成有限長(zhǎng)的正弦序列Xa=T*fftshift(fft(xa));Xa(1)%求xa的對(duì)稱FFT，k=floor(-(N-1)/2:(N-1)/2);%頻率下標(biāo)向量

plot(k*D,abs(Xa)) %繪圖

程序運(yùn)行的結(jié)果見(jiàn)下圖。88連續(xù)信號(hào)的頻譜計(jì)算89連續(xù)信號(hào)的頻譜計(jì)算把N作為用戶自選的參數(shù)。選N=50,100,500及628所得的計(jì)算結(jié)果在圖中畫(huà)出，在±5處都出現(xiàn)了兩個(gè)尖峰。結(jié)論是：如果用FFT計(jì)算出的幅頻譜，包含了窄的尖峰，而N每加一倍，它也大體上加倍，這就說(shuō)明這個(gè)連續(xù)時(shí)間信號(hào)包含一個(gè)周期分量。因?yàn)榧夥逭艘院?，很難準(zhǔn)確采樣在峰值處，所以圖上顯示的峰值，并非真正的最大值。從理論上說(shuō)，只有長(zhǎng)度L恰好是周期的整數(shù)倍時(shí)，可以得到準(zhǔn)確的峰值。

90連續(xù)信號(hào)的頻譜計(jì)算例4.5.5求下列周期信號(hào)的頻譜xa(t)=-1+2sin0.2πt-3cosπt解：它的最高頻率是Ωm=π弧度／秒，因此采樣周期T必須小于π／Ωm＝1。任選T＝0.2，在圖（a）中畫(huà)出N＝50時(shí)的幅頻特性，而在圖(b)中用N＝1000。

從圖(b)中也可大致看出，在ω＝0，±0.628和±3.14的三個(gè)高度之比為2：2：3，由于信號(hào)的直流分量是它的零階傅立葉級(jí)數(shù)的一半，可見(jiàn)原信號(hào)的三個(gè)頻率分量幅度之比為1：2：3，91連續(xù)信號(hào)的頻譜計(jì)算92連續(xù)信號(hào)的頻譜計(jì)算用FFT從連續(xù)信號(hào)計(jì)算其頻譜的步驟概括如下：連續(xù)信號(hào)xa(t)―→(采樣)

離散序列x=xa(nT)―→(FFT)

序列頻譜X(k)=fft(x)―→(乘T)信號(hào)頻譜Xa=T?X1.采樣前的xa本來(lái)是解析式，但在MATLAB中和x同樣用數(shù)組表示，所以在程序中用同一變量表示；2.序列樣本經(jīng)FFT后變?yōu)樾蛄蓄l譜，其頻率混疊誤差由折疊頻率處的幅特性決定。該處的幅特性愈小愈好。3.控制混疊誤差和頻率分辨率的參數(shù)是T,N和L。4.序列頻譜必須乘T才等于信號(hào)頻譜。這也表明信號(hào)和它的采樣序列物理上不等價(jià)。93循環(huán)計(jì)算中對(duì)應(yīng)頻點(diǎn)的確定

當(dāng)要靠計(jì)算機(jī)自動(dòng)選擇參數(shù)N,T或L時(shí)，需要比較相鄰兩次計(jì)算中對(duì)應(yīng)頻點(diǎn)上幅特性的誤差。本節(jié)討論在N,T或L發(fā)生變換的條件下，如何找到對(duì)應(yīng)的頻點(diǎn)Ωk？

因?yàn)楫?dāng)相鄰兩次計(jì)算中，N、T或L發(fā)生變化，必須要在同樣的Ωk下進(jìn)行比較。94循環(huán)計(jì)算中對(duì)應(yīng)頻點(diǎn)的確定令得到：只有滿足(4.5.7)式的整數(shù)k1和k2，才是可以進(jìn)行頻譜比較的對(duì)應(yīng)點(diǎn)。這是從數(shù)學(xué)上說(shuō)，如果要得到更形象的結(jié)果，可以用由程序hc456生成的圖4.5.7來(lái)觀察。

95循環(huán)計(jì)算中對(duì)應(yīng)頻點(diǎn)的確定96循環(huán)計(jì)算中對(duì)應(yīng)頻點(diǎn)的確定以子圖一作為標(biāo)準(zhǔn)，T=0.1秒，N=16，信號(hào)長(zhǎng)度L=1.6秒，乃奎斯特頻率范圍為 ；內(nèi)分16個(gè)點(diǎn)。把T減小一半，N不變，L也減小一半，得到子圖二。其乃奎斯特頻率范圍擴(kuò)大了一倍，成為 內(nèi)分16個(gè)點(diǎn)，間隔大了一倍。故只有8個(gè)點(diǎn)能與標(biāo)準(zhǔn)情況對(duì)應(yīng)。把T減小一半，同時(shí)把N加大一倍，因此L不變，得到子圖三。其乃奎斯特頻率范圍也擴(kuò)大了一倍，分成32個(gè)點(diǎn)，故中間16個(gè)能與標(biāo)準(zhǔn)情況對(duì)應(yīng)。T不變，N加大一倍，L也加大一倍，得到子圖四。其乃奎斯特頻率范圍不變，分成32個(gè)點(diǎn)，間隔D縮小一半，每隔一個(gè)點(diǎn)有一個(gè)點(diǎn)與標(biāo)準(zhǔn)情況對(duì)應(yīng)。974.6用反變換從頻譜計(jì)算信號(hào)

從連續(xù)頻譜計(jì)算信號(hào)的步驟應(yīng)為求頻譜過(guò)程的逆：信號(hào)頻譜Xa(Ω)→(頻率域采樣除T)

離散序列頻譜X(k)=Xa/T→(IFFT)時(shí)間序列x(n)=fft(X)→(取包絡(luò))連續(xù)信號(hào)xa信號(hào)頻譜必須采樣并除T才等于序列頻譜；序列頻譜經(jīng)IFFT后變?yōu)闀r(shí)間序列，其時(shí)間混疊誤差由時(shí)間序列中幅度最小的區(qū)域決定。該區(qū)的幅度愈小愈好。3.控制混疊誤差和時(shí)間分辨率的參數(shù)是T,N和L。4.

xa是x的包絡(luò)，可以用plot(nT,x)畫(huà)出。但作為變量，xa在MATLAB中和x用同一變量表示。98頻率域采樣定理

X(jω)

的IDTFT是有限長(zhǎng)序列x(n)；其頻率域采樣序列X(k)(以后記作Y(k))的IDFT是無(wú)限長(zhǎng)周期序列的主值y(n)。計(jì)算機(jī)不能做IDTFT，只能做IDFT，因而只能得到y(tǒng)(n)，得不到x(n)。從概念上說(shuō)，X(jω)中包含的信息是完全的，Y(k)只是其部分樣本，是不完全的。所以，知道X可以確定Y；知道Y未必能完全確定X。對(duì)應(yīng)于時(shí)域信號(hào)。知道x(n)可以確定y(n)；知道y(n)未必能完全確定x(n)。頻率域采樣定理從數(shù)學(xué)上求出兩者的數(shù)學(xué)關(guān)系。99證明的基本思路是：寫(xiě)出x(n)=IDTFT[X(jω)]；寫(xiě)出y(n)=IDFT[Y(k)]；通過(guò)Y(k)=X(jωk)=X(j2πk/N)把兩者關(guān)聯(lián)起來(lái)；頻率域采樣定理的證明100注意 ,即周期為N,把x(n)分成 長(zhǎng)N的段。段號(hào)為r，段內(nèi)下標(biāo)為n1?[0,N-1]，則x(n)=x(n1+rN)，r=0,±1,±2,…

（4.6.4）中n1相同的點(diǎn)的乘子 相同，可以合并得此式與(4.6.2)式恰成傅立葉變換對(duì)，故有可見(jiàn)，y(n)是無(wú)限個(gè)x(n)平移相加的結(jié)果。101由（4.6.7)，如果已知x(n)，可以將它無(wú)限次平移N相加而得到y(tǒng)(n)。如果已知y(n)，則未必可能求出x(n)。只有當(dāng)x(n)在[0,N-1)主值區(qū)間外取值為零，即此時(shí)x(n)平移后不互相重疊,即無(wú)時(shí)間泄漏。

x(n)=y(n)(n=0,1,…,N-1)這就要把N取得足夠大，即頻率域采樣足夠密。頻率域采樣定理102從頻譜計(jì)算離散時(shí)間序列為了利用IFFT函數(shù)，輸入的頻率樣本點(diǎn)必須合乎IFFT的要求。(1).IFFT輸入的N個(gè)樣本應(yīng)是數(shù)字頻率ω＝ΩT,它的最大值不得超過(guò)2π。假如頻譜的最大頻率為Ωm，采樣周期T必須滿足：式中，F(xiàn)s為采樣頻率(赫茲)。T或Fs確定后，才能把給定的頻譜X（Ω）映射到數(shù)字頻率坐標(biāo)上來(lái)，成為X（ω）。頻率分辨率為

Δω=2π/N或D=Δω/T=2π/(NT)103從頻譜計(jì)算離散時(shí)間序列(2)通常序列的DTFT是用對(duì)稱于原點(diǎn)的正負(fù)頻率表示，頻率下標(biāo)取為

k＝floor(-(N-1)/2:(N-1)/2)

而IFFT的輸入必須從k=0到N-1，因此，頻率樣本點(diǎn)中k<0的部分（用kn表示）必須平移到k≥0部分（用kp表示）的右面。

(3)若待求的時(shí)間序列是實(shí)的，則輸入頻譜X(Ω)應(yīng)具有共軛對(duì)稱性。X(k)應(yīng)具有循環(huán)共軛對(duì)稱性。所以把k作上述移動(dòng)不影響其對(duì)稱性。

104從頻譜計(jì)算離散時(shí)間序列由上節(jié)，頻率采樣會(huì)引入時(shí)間混疊。如果時(shí)間混疊不能忽略，那樣算出的結(jié)果是無(wú)用的。檢查時(shí)間混疊的方法如下：如果算出x1(n)在所有的n處都明顯地不等于零，說(shuō)明時(shí)間混疊相當(dāng)大，必須另選大些的N。當(dāng)算出的x1(n)在n的某個(gè)范圍內(nèi)實(shí)際上等于零了，即IDFT存在著(4.6.8)式中含零的部分，就認(rèn)可這個(gè)N。105從頻譜計(jì)算離散時(shí)間序列例4.6.1設(shè)某時(shí)間序列的頻譜為要求用IDFT數(shù)值計(jì)算求出它的時(shí)間序列。解:因?yàn)轭}中用的是數(shù)字頻率而沒(méi)有涉及T，這就是默認(rèn)T=1.先取N=3，輸入頻率數(shù)組必須是把2π等分為3份的數(shù)字頻率。然后再取N=10，用同一程序，對(duì)計(jì)算結(jié)果進(jìn)行比較。106從頻譜計(jì)算離散時(shí)間序列

程序hc461T=1;N=3;D=2*pi/N; %N及頻率分辨率kn=floor(-(N-1)/2:-1/2);%負(fù)頻率下標(biāo)向量kp=floor(0:(N-1)/2); %正頻率下標(biāo)向量w=[kp,kn]*D;%形成主值區(qū)頻率排序%按新的頻率排序輸入數(shù)字頻譜X=2-exp(-j*w)+exp(-j*2*w)+exp(-j*3*w);x=ifft(X);%對(duì)X求IFFT,得出循環(huán)對(duì)稱的xstem(T*[0:N-1],x)%繪圖107從頻譜計(jì)算離散時(shí)間序列108從頻譜計(jì)算離散時(shí)間序列當(dāng)N=3時(shí)，因?yàn)镹小于時(shí)間序列長(zhǎng)度，產(chǎn)生了時(shí)間混疊，其標(biāo)志是序列中沒(méi)有取零的點(diǎn)或點(diǎn)列，結(jié)果不正確；當(dāng)N=10>4時(shí)，結(jié)果是正確的，其標(biāo)志則是在算出的IDFT序列x(n)中，包含一段如(4.6.8)式所指的全零的部分，

109從頻譜計(jì)算離散時(shí)間序列例4.6.2給出在頻率范圍[-62.8,62.8]弧度/秒間的頻譜。求出它的時(shí)間序列。解：本題給出的是模擬頻率Ω，Ωm=62.8取T≤π/Ωm=0.05秒。即采樣頻率為20赫茲或40π。N取8和20兩種情況。

110從頻譜計(jì)算離散時(shí)間序列wm=62.83;T=pi/wm; %選擇采樣周期TN0=[8,20]; %設(shè)定兩種N值forr=1:2%作兩次循環(huán)計(jì)算

N=N0(r);D=2*pi/(N*T);%取N,求頻率分辨率

k=[0:N-1]+eps;%求頻率下標(biāo)，偏移微量

X=sin(0.275*k*D)./sin(0.025*k*D);算出頻譜

x=ifft(X);%求X的IDFT，得出x(n)subplot(1,2,r),stem([0:N-1]*T,x)%繪圖end111從頻譜計(jì)算離散時(shí)間序列112從頻譜計(jì)算離散時(shí)間序列N=8時(shí)得到左圖，其中的8個(gè)數(shù)據(jù)中沒(méi)有一個(gè)為零，可見(jiàn)存在著頻率采樣所造成的時(shí)間混疊，再取N=20進(jìn)行計(jì)算。結(jié)果畫(huà)在右圖中，算得的序列中都有1