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

第四章

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

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

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

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

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

模擬信號時域采樣周期延拓主值區(qū)間時域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)系

:采樣插補

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

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

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

functiony=myditfft(x)

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

%若x的長度不是2的冪，補零到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的初始值

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

對每次作DFT運算

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

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

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

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

t=y(kp)*u;%蝶形運算的乘積項

y(kp)=y(k)-t;%蝶形運算的加法項

y(k)=y(k)+t;%蝶形運算的減法項

end

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

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

兩端取共軛：

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

244.3用FFT計算序列的頻譜

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

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

25用FFT算序列頻譜例4.3.1考慮長度為5的有限序列:x=[1,3,5,3,1],要求用FFT來計算其頻譜。解：因為要求的是頻率域中的連續(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));%求對稱位置向量X=fftshift(fft(x,N));%求對稱的FFT序列值subplot(1,2,1),plot(k*D,abs(X),'o:')%畫幅特性subplot(1,2,2),plot(k*D,angle(X),'o:')%相特性圖26用FFT計算序列的頻譜曲線如下。注意plot語句具有線性內(nèi)插作用，可近似得出連續(xù)頻譜DTFT.27用FFT計算序列的頻譜用FFT計算了頻譜時，其頻率樣本數(shù)等于序列樣本數(shù)，上例中于N=5，所以相鄰兩個頻率樣本點的間距為這個頻率間距非常大，分辨率很差。改善分辨率的最好方法是增加數(shù)據(jù)的長度N。如果沒有更多數(shù)據(jù)，也可以給輸入序列補零。補了1019個零以后，得出連續(xù)光滑的頻譜。不過改善的只是頻譜的視在分辨率。

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

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

在編寫程序時，要準(zhǔn)備給出不同的N進(jìn)行比較。所以程序hc432的寫法應(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的對稱位置FFTsubplot(1,2,r),plot(k*D,X)%畫實頻特性圖end30用FFT計算序列的頻譜由圖，N=20給出了頻譜的大體形狀，但分辨率差，N=1024的結(jié)果則接近于精確的結(jié)果。

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

32用FFT計算無限序列的頻譜當(dāng)把一個無限序列截斷成為有限序列，會導(dǎo)致泄漏和波動，算出的頻譜與無限序列的頻譜有誤差。因此要確定截斷的長度至少應(yīng)有多長，才能保證其頻譜足夠準(zhǔn)確。

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

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

36比較N1=N和N2=2N的計算頻點。用N1點算出的頻譜X1位于下列頻點上

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

(4.3.5)令 ，解出k2=2k1。必須在的這些頻點上比較X1和X2的幅度。按照這個原理來編寫程序hc424。

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

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

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

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

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

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

誤差b%24.015.76480.33230.00111.22e-8394.4

時域采樣中的頻譜變換

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

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

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

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

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

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

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

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

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

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

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

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

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

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

理想方法：

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

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

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

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

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

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

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

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

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

78連續(xù)信號的頻譜計算為了用一個程序而，用不同的參數(shù)畫出圖4.5.1的多幅子圖。可以用for循環(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)計算六次

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

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

程序運行結(jié)果說明，N取11時分辨率太低，取T=0.5,用1024點FFT,(相當(dāng)于后面補零到1024）可得到較好結(jié)果。84連續(xù)信號的頻譜計算85連續(xù)信號的頻譜計算周期信號的頻譜計算：從理論上說，如果一個連續(xù)時間信號是周期性的，它必定有無限的長度，計算出的頻譜將不會收斂。然而，用計算機計算時，序列總是有限長的。因為我們遇到的都是經(jīng)過截斷了的周期信號。從算法上看，完全可以用上節(jié)中的方法和程序，只是在計算結(jié)果上，頻譜會有向有限個頻點上集中，并成為脈沖的趨勢。

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

解：在計算機計算中，正余弦函數(shù)必須截斷為有限長度L。信號＝cos5t的帶寬限制于5。純理論地看，只要采樣周期小于π/5=0.63秒，就不會發(fā)生頻率混疊。然而如果把cos5t截斷為長L的信號，則它的頻譜就不再是有限帶寬了，所以必須采用更小的采樣周期，任選T＝0.1，并選N＝50，得到L＝TN＝5。按此來截斷信號。

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

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

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

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

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

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

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

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

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

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

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

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

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

X(jω)

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

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

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

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

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

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

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

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

程序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);%對X求IFFT,得出循環(huán)對稱的xstem(T*[0:N-1],x)%繪圖107從頻譜計算離散時間序列108從頻譜計算離散時間序列當(dāng)N=3時，因為N小于時間序列長度，產(chǎn)生了時間混疊，其標(biāo)志是序列中沒有取零的點或點列，結(jié)果不正確；當(dāng)N=10>4時，結(jié)果是正確的，其標(biāo)志則是在算出的IDFT序列x(n)中，包含一段如(4.6.8)式所指的全零的部分，

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

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

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從頻譜計算離散時間序列112從頻譜計算離散時間序列N=8時得到左圖，其中的8個數(shù)據(jù)中沒有一個為零，可見存在著頻率采樣所造成的時間混疊，再取N=20進(jìn)行計算。結(jié)果畫在右圖中，算得的序列中都有1