數(shù)字信號(hào)處理

1、第一次討論課內(nèi)容（1）時(shí)域信號(hào)（a）如何由模擬信號(hào)產(chǎn)生時(shí)域離散信號(hào)；模擬信號(hào)：信號(hào)的自變量和函數(shù)值都為連續(xù)值。一段連續(xù)的時(shí)間間隔內(nèi)，其代表信息的特征量可以在任意瞬間呈現(xiàn)為任意數(shù)值的信號(hào)。時(shí)域離散信號(hào)：自變量取離散值，函數(shù)值取連續(xù)值。離散信號(hào)是一個(gè)序列，即其自變量是“離散”的。這個(gè)序列的每一個(gè)值都可以被看作是連續(xù)信號(hào)的一個(gè)采樣。模擬信號(hào)，以采樣間隔T對(duì)它進(jìn)行等間隔采樣，得到時(shí)域離散信號(hào)。即： (取整數(shù)， )采樣是將時(shí)間上、幅值上都連續(xù)的模擬信號(hào)，在采樣脈沖的作用，轉(zhuǎn)換成時(shí)間上離散（時(shí)間上有固定間隔）、但幅值上仍連續(xù)的離散模擬信號(hào)。所以采樣又稱為波形的離散化過程。 對(duì)模擬信號(hào)進(jìn)行采樣可以看作一個(gè)模

2、擬信號(hào)通過一個(gè)電子開關(guān)S。設(shè)電子開關(guān)的作用等效成一寬度為，周期為T的矩形脈沖串，采樣信號(hào)就是與相乘的結(jié)果。理想采樣：0，脈沖串變?yōu)閱挝粵_激串。中每個(gè)單位沖激處在采樣點(diǎn)上，強(qiáng)度為1，理想采樣是與相乘的結(jié)果。即： 式中是單位沖激信號(hào)，上式只有在時(shí)，才可能是非零值，所以可寫成:由模擬信號(hào)經(jīng)采樣產(chǎn)生時(shí)域離散信號(hào)的MATLAB程序： 一個(gè)連續(xù)的周期性三角波信號(hào)頻率為50HZ，信號(hào)幅度在0+2V之間，在窗口顯示2個(gè)周期信號(hào)波形，對(duì)信號(hào)的一個(gè)周期進(jìn)行16點(diǎn)采樣來獲取離散信號(hào)。代碼：f=50;Um=1;ts=2; %輸入信號(hào)的頻率、振幅、顯示周期N=16; %信號(hào)一個(gè)采樣周期的采樣點(diǎn)數(shù)為16T=1/f; %信

3、號(hào)周期Tdt=T/N; %采樣時(shí)間間隔n=0:ts*N-1; %建立離散時(shí)間的時(shí)間序列列tn=n*dt; %確定時(shí)間序列樣點(diǎn)在時(shí)間軸上的位置y=Um*sawtooth(2*f*pi*tn,1/2)+1; %三角波subplot(2,1,1),stem(tn,y); %顯示采樣后的信號(hào)title('離散信號(hào)');subplot(2,1,2),plot(tn,y); %顯示原連續(xù)信號(hào)title('連續(xù)信號(hào)');運(yùn)行結(jié)果： 一個(gè)連續(xù)的周期性正弦信號(hào)頻率為200Hz，信號(hào)幅度在-1+1V之間，在窗口上顯示2個(gè)周期信號(hào)波形，用Fs=4KHZ的頻率對(duì)連續(xù)信號(hào)進(jìn)行采樣，試顯示

4、連續(xù)信號(hào)和采樣獲得的離散信號(hào)波形。代碼：f=200;Um=1;nt=2; %輸入信號(hào)的頻率、振幅、顯示周期Fs=4000; %采樣頻率N=Fs/f; %采樣點(diǎn)數(shù)T=1/f; %信號(hào)周期 dt=T/N; %采樣時(shí)間間隔n=0:nt*N-1; %建立離散時(shí)間的時(shí)間序列tn=n*dt; %確定時(shí)間序列樣點(diǎn)在時(shí)間軸上的位置y=Um*sin(2*f*pi*tn); %正弦波subplot(2,1,1),stem(tn,y); %顯示采樣后的信號(hào)title('離散信號(hào)');subplot(2,1,2),plot(tn,y); %顯示原連續(xù)信號(hào)title('連續(xù)信號(hào)');運(yùn)行

5、結(jié)果：抽樣定理是連續(xù)時(shí)間信號(hào)和離散時(shí)間信號(hào)之間的橋梁，在時(shí)域該系統(tǒng)實(shí)現(xiàn)了輸入信號(hào)與抽樣序列的相乘，完成了時(shí)間軸的離散，在頻域?qū)崿F(xiàn)了原信號(hào)頻譜的周期延拓。在奈奎斯特抽樣定理的條件下(抽樣頻率不小于被抽樣帶限信號(hào)最高頻率的兩倍)，一個(gè)連續(xù)時(shí)間信號(hào)完全可以用該信號(hào)在等時(shí)間間隔點(diǎn)上的樣本來表示，在頻率軸上實(shí)現(xiàn)了原信號(hào)頻譜無混疊的周期化。采樣定理：設(shè)連續(xù)信號(hào)的的最高頻率為，如果采樣頻率，那么采樣信號(hào)可以唯一的恢復(fù)出原連續(xù)信號(hào)，否則會(huì)造成采樣信號(hào)中的頻譜混疊現(xiàn)象，不可能無失真地恢復(fù)原連續(xù)信號(hào)。驗(yàn)證采樣定理的MATLAB程序：畫出連續(xù)時(shí)間信號(hào)的時(shí)域波形及其幅頻特性曲線，信號(hào)為 對(duì)信號(hào)進(jìn)行采樣，得到采樣序列，

6、繪制采樣頻率分別為80Hz，120 Hz，150 Hz時(shí)的采樣序列波形；、對(duì)不同采樣頻率下的采樣序列進(jìn)行頻譜分析，繪制其幅頻曲線，對(duì)比各頻率下采樣序列和的幅頻曲線有無差別。由采樣序列恢復(fù)出連續(xù)時(shí)間信號(hào)，畫出其時(shí)域波形，對(duì)比與原連續(xù)時(shí)間信號(hào)的時(shí)域波形。代碼：%實(shí)現(xiàn)采樣頻譜分析繪圖函數(shù)function fz=caiyang(fy,fs) %第一個(gè)輸入變量是原信號(hào)函數(shù)，信號(hào)函數(shù)fy以字符串的格式輸入%第二個(gè)輸入變量是采樣頻率fs0=10000; tp=0.1; t=-tp:1/fs0:tp; k1=0:999; k2=-999:-1; m1=length(k1); m2=length(k2); f=

7、fs0*k2/m2,fs0*k1/m1; %設(shè)置原信號(hào)的頻率數(shù)組w=-2*pi*k2/m2,2*pi*k1/m1; fx1=eval(fy); FX1=fx1*exp(-j*1:length(fx1)'*w); %求原信號(hào)的離散時(shí)間傅里葉變換figure % 畫原信號(hào)波形subplot(2,1,1),plot(t,fx1,'r') title('原信號(hào)'), xlabel('時(shí)間t (s)') axis(min(t),max(t),min(fx1),max(fx1)% 畫原信號(hào)幅度頻譜subplot(2,1,2),plot(f,abs(F

8、X1),'r') title('原信號(hào)幅度頻譜') , xlabel('頻率f (Hz)') axis(-100,100,0,max(abs(FX1)+5) % 對(duì)信號(hào)進(jìn)行采樣Ts=1/fs; %采樣周期t1=-tp:Ts:tp; %采樣時(shí)間序列f1=fs*k2/m2,fs*k1/m1; %設(shè)置采樣信號(hào)的頻率數(shù)組t=t1; %變量替換fz=eval(fy);%獲取采樣序列FZ=fz*exp(-j*1:length(fz)'*w); %采樣信號(hào)的離散時(shí)間傅里葉變換figure % 畫采樣序列波形subplot(2,1,1),stem(t,f

9、z,'.'), title('取樣信號(hào)') , xlabel('時(shí)間t (s)') line(min(t),max(t),0,0) % 畫采樣信號(hào)幅度頻譜subplot(2,1,2),plot(f1,abs(FZ),'m') title('取樣信號(hào)幅度頻譜') , xlabel('頻率f (Hz)')function fh=huifu(fz,fs)%第一個(gè)輸入變量是采樣序列%第二個(gè)輸入變量是得到采樣序列所用的采樣頻率T=1/fs; dt=T/10; tp=0.1; t=-tp:dt:tp; n=-t

10、p/T:tp/T; TMN=ones(length(n),1)*t-n'*T*ones(1,length(t); fh=fz*sinc(fs*TMN); % 由采樣信號(hào)恢復(fù)原信號(hào)k1=0:999; k2=-999:-1; m1=length(k1); m2=length(k2); w=-2*pi*k2/m2,2*pi*k1/m1; FH=fh*exp(-j*1:length(fh)'*w); % 恢復(fù)后的信號(hào)的離散時(shí)間傅里葉變換figure % 畫恢復(fù)后的信號(hào)的波形subplot(2,1,1),plot(t,fh,'g'), st1=sprintf('由

11、取樣頻率fs=%d',fs); st2='恢復(fù)后的信號(hào)' st=st1,st2; title(st) , xlabel('時(shí)間t (s)') axis(min(t),max(t),min(fh),max(fh) line(min(t),max(t),0,0) % 畫重構(gòu)信號(hào)的幅度頻譜f=10*fs*k2/m2,10*fs*k1/m1; %設(shè)置頻率數(shù)組subplot(2,1,2),plot(f,abs(FH),'g') title('恢復(fù)后信號(hào)的頻譜') , xlabel('頻率f (Hz)') axis(-

12、100,100,0,max(abs(FH)+2);%主函數(shù)f1='sin(2*pi*60*t)+cos(2*pi*25*t)+cos(2*pi*30*t)'%輸入一個(gè)信號(hào)fs0=caiyang(f1,80);%欠采樣fr0=huifu(fs0,80); fs1=caiyang(f1,120);%臨界采樣fr1=huifu(fs1,120); fs2=caiyang(f1,150);%過采樣fr2=huifu(fs2,150); 運(yùn)行結(jié)果：(1) 采樣頻率時(shí)，為原信號(hào)的欠采樣信號(hào)和恢復(fù)，采樣頻率不滿足時(shí)域采樣定理，那么頻移后的各相臨頻譜會(huì)發(fā)生相互重疊，這樣就無法將他們分開，因而也

13、不能再恢復(fù)原信號(hào)。頻譜重疊的現(xiàn)象被稱為混疊現(xiàn)象。欠采樣信號(hào)的離散波形及頻譜、恢復(fù)后信號(hào)見下圖。(2) 采樣頻率 時(shí)，為原信號(hào)的臨界采樣信號(hào)和恢復(fù)，下圖為其采樣的離散波形和頻譜，從下圖恢復(fù)后信號(hào)和原信號(hào)先對(duì)比可知，只恢復(fù)了低頻信號(hào)，高頻信號(hào)未能恢復(fù)。（3）采樣頻率時(shí)，為原信號(hào)的過采樣信號(hào)和恢復(fù)，由圖采樣信號(hào)離散波形和頻譜，可以看出采樣信號(hào)的頻譜是原信號(hào)頻譜進(jìn)行周期延拓形成的，從圖采樣恢復(fù)后的波形和頻譜，可看出與原信號(hào)誤差很小了，說明恢復(fù)信號(hào)的精度已經(jīng)很高。(b) 線性時(shí)不變系統(tǒng)輸入和輸出之間的關(guān)系設(shè)系統(tǒng)的輸入用表示，表示成單位脈沖序列移位加權(quán)和為： 則系統(tǒng)輸出為： 根據(jù)線性系統(tǒng)疊加性質(zhì)：又根據(jù)時(shí)

14、不變性質(zhì) 即：線性時(shí)不變系統(tǒng)的輸出等于輸入序列和該系統(tǒng)的單位脈沖響應(yīng)的卷積。驗(yàn)證例子：設(shè)模擬正弦信號(hào)：，現(xiàn)以T為周期，對(duì)其進(jìn)行采樣，得到正弦序列為： 其中為數(shù)字頻率，為數(shù)字角頻率。產(chǎn)生信號(hào)： N取大于500設(shè)離散時(shí)間線性時(shí)不變系統(tǒng)為：令通過此系統(tǒng)MATLAB程序如下：n=0:500;xa=sin(0.0364*pi*n);figure(1);stem(n,xa,'.');xlabel('n');ylabel('xa(n)');hn(1:50)=1;y1=conv(xa,hn);figure(2);stem(y1,'.');xlab

15、el('n');ylabel('y1(n)');y2=filter(hn,1,xa);figure(3);stem(y2,'.');axis(0,600,-20,20)xlabel('n');ylabel('y2(n)');運(yùn)行結(jié)果：的時(shí)域波性卷積conv方法：通過系統(tǒng)的輸出波形濾波Filter函數(shù)方法：通過系統(tǒng)的輸出波形兩種函數(shù)運(yùn)行結(jié)果相同，證明了線性時(shí)不變系統(tǒng)的輸出等于輸入序列和該系統(tǒng)的單位脈沖響應(yīng)的卷積。（2）離散傅立葉變換(a) 離散傅立葉變換的性質(zhì)；設(shè)是一個(gè)長(zhǎng)度為M的有限長(zhǎng)序列，則定義的N點(diǎn)離散傅里葉變換

16、為： 1、 線性性質(zhì)如果兩個(gè)有限長(zhǎng)序列分別為和x2(n)，長(zhǎng)度分別為N1和N2，且y(n)ax1(n)bx2(n) (a、b均為常數(shù))則該y(n)的N點(diǎn)DFT為Y(k)DFTy(n)aX1(k)bX2(k) 0kN1其中：NmaxN1，N2，X1(k)和X2(k)分別為x1(n)和x2(n)的N點(diǎn)DFT。驗(yàn)證線性性質(zhì)例子：已知x1(n)0，1，2，4，x2(n)1，0，1，0，1，(1)y(n)2x1(n)3x2(n)，再由y(n)的N點(diǎn)DFT獲得Y1(k)；(2)由x1(n)、x2(n)求X1(k)、X2(k)，再求Y2(k)2X1(k)3X2(k)。MATLAB程序：xn1=0,1,2,4

17、; %建立xn1序列xn2=1,0,1,0,1;%建立xn2序列N1=length(xn1);N2=length(xn2);N=max(N1,N2); %確定Nif N1>N2 xn2=xn2,zeros(1,N1-N2);%對(duì)長(zhǎng)度短的序列補(bǔ)0elseif N2>N1 xn1=xn1,zeros(1,N2-N1);endyn=2*xn1+3*xn2;%計(jì)算ynn=0:N-1;k=0:N-1;Yk1=yn*(exp(-j*2*pi/N).(n'*k) %求yn的N點(diǎn)DFTXk1=xn1*(exp(-j*2*pi/N).(n'*k); %求xn1的N點(diǎn)DFTXk2=xn

18、2*(exp(-j*2*pi/N).(n'*k); %求xn2的N點(diǎn)DFTYk2=2*Xk1+3*Xk2%由Xk1、Xk2求Yksubplot(3,2,1),stem(n,xn1);title('x1(n)');subplot(3,2,2),stem(n,Xk1);title('X1(k)');subplot(3,2,3),stem(n,xn2);title('x2(n)');subplot(3,2,4),stem(n,Xk2);title('X1(k)');%subplot(4,2,5),stem(n,yn);%titl

19、e('yn');subplot(3,2,6),stem(n,Yk2);axis(-1,5,-20,30)title('2*Xk1+3*Xk2');subplot(3,2,5),stem(n,Yk1);axis(-1,5,-20,30)title('DFTy(n)');運(yùn)行結(jié)果：Yk1 = 23.0000 + 0.0000i -7.5902 + 1.5388i 3.5902 - 0.3633i 3.5902 + 0.3633i -7.5902 - 1.5388iYk2 = 23.0000 + 0.0000i -7.5902 + 1.5388i 3.

20、5902 - 0.3633i 3.5902 + 0.3633i -7.5902 - 1.5388i由運(yùn)算結(jié)果和波形圖都可知，Y1(K)=Y2(K)，即y(n)2x1(n)3x2(n)，DFTy(n)= 2X1(k)3X2(k),可驗(yàn)證線性性質(zhì)。2、 循環(huán)移位性質(zhì)如果有限長(zhǎng)序列為x(n)，長(zhǎng)度為N，將x(n)的循環(huán)移位定義為：x(n)左移m位的過程可由以下步驟獲得：(1)將x(n)以N為周期進(jìn)行周期延拓，得到 ；(2)將左移m位，得到 ；(3)取的主值序列，得到x(n)循環(huán)移位序列y(n)。驗(yàn)證例子：已知有限長(zhǎng)序列x(n)1，2，3，4，5，6，求x(n)左移2位成為新的向量y(n)，并畫出循環(huán)

21、移位的中間過程。MATLAB程序：xn=1,2,3,4,5,6; %建立xn序列Nx=length(xn);nx=0:Nx-1;nx1=-Nx:2*Nx-1;%設(shè)立周期延拓的范圍x1=xn(mod(nx1,Nx)+1);%建立周期延拓序列ny1=nx1-2;y1=x1;%將x1左移2位,得到y(tǒng)1RN=(nx1>=0)&(nx1<Nx);%在x1的位置向量nx1上設(shè)置主值窗RN1=(ny1>=0)&(ny1<Nx);%在y1的位置向量ny1上設(shè)置主值窗subplot(4,1,1),stem(nx1,RN.*x1);%畫出x1的主值部分subplot(4,1

22、,2),stem(nx1,x1);%畫出x1subplot(4,1,3),stem(ny1,y1);%畫出y1subplot(4,1,4),stem(ny1,RN1.*y1); %畫出y1的主值部分運(yùn)行結(jié)果：可見，循環(huán)移位的實(shí)質(zhì)是將x(n)左移m位時(shí)，移出主值區(qū)的序列值又依次從右端進(jìn)入主值區(qū)。有限長(zhǎng)序列的移位也稱為循環(huán)移位。3、時(shí)域和頻域循環(huán)卷積特性離散傅里葉變換的循環(huán)卷積特性也稱為圓周卷積，分為時(shí)域卷積和頻域卷積兩類。1)時(shí)域循環(huán)卷積假定有限長(zhǎng)序列長(zhǎng)度分別為N1、N2，。 2)頻域循環(huán)卷積利用時(shí)域和頻域的對(duì)稱性，可以得到頻域卷積特性。若 ，則 下面重點(diǎn)討論時(shí)域循環(huán)卷積。時(shí)域循環(huán)卷積的方法有多

23、種：方法1：直接使用時(shí)域循環(huán)卷積。由于有限長(zhǎng)序列可以看成是周期序列的主值，因此，時(shí)域圓周卷積的結(jié)果可以由對(duì)應(yīng)的周期序列卷積和取主值部分獲得。方法2：用頻域DFT相乘再求逆變換。即先分別求x1(n)、x2(n)的DFTX1(k)、X2(k)，再求Y(k)的IDFT獲得y(n)。方法3：用FFT和IFFT進(jìn)行循環(huán)卷積?；舅悸吠椒?，但直接使用了MATLAB提供的fft和ifft子函數(shù)來實(shí)現(xiàn)。舉例：已知的兩個(gè)時(shí)域周期序列分別取主值，得到x11，1，1，0，0，0，x20，1，2，3，0，0，求時(shí)域循環(huán)卷積y(n)并用圖形表示。xn1=0,1,2,3,0,0; %建立x1(n)序列xn2=1,1,

24、1,0,0,0; %建立x2(n)序列N=length(xn1);n=0:N-1;k=0:N-1;Xk1=xn1*(exp(-j*2*pi/N).(n'*k); %由x1(n)的DFT求X1(k)Xk2=xn2*(exp(-j*2*pi/N).(n'*k); %由x2(n)的DFT求X2(k)Yk=Xk1.*Xk2; %Y(k)=X1(k)X2(k)yn=Yk*(exp(j*2*pi/N).(n'*k)/N; %由Y(k)的IDFT求y(n)yn=abs(yn) %取模值,消除DFT帶來的微小復(fù)數(shù)影響subplot(2,3,1),stem(n,xn1);title(&#

25、39;x1(n)')subplot(2,3,2),stem(n,xn2);title('x2(n)')subplot(2,3,3),stem(n,yn);title('y(n)')subplot(2,3,4),stem(n,Xk1);title('X1(k)')subplot(2,3,5),stem(n,Xk2);title('X2(k)')subplot(2,3,6),stem(n,Yk);title('Y(k)')運(yùn)行結(jié)果：得到：yn = 0.0000 1.0000 3.0000 6.0000 5.00

26、00 3.00004、共軛對(duì)稱性由于序列x(n)及其離散傅里葉變換X(k)的定義在主值為0N1的區(qū)間，因此DFT的循環(huán)對(duì)稱性對(duì)時(shí)間序列是指關(guān)于n0和nN/2的對(duì)稱性，對(duì)頻譜序列是關(guān)于數(shù)字頻率為0和p的對(duì)稱性。本實(shí)驗(yàn)重點(diǎn)分析實(shí)序列的循環(huán)對(duì)稱性。實(shí)序列x(n)可以分解為循環(huán)偶序列xe(n)和循環(huán)奇序列xo(n)：x(n)xe(n)xo(n) 0nN1其中：設(shè)DFTx(n)X(k)ReX(k)j*ImX(k)，則有即實(shí)序列中的偶序列xe(n)對(duì)應(yīng)于x(n)的離散傅里葉變換X(k)的實(shí)部，而實(shí)序列中的奇序列xo(n)對(duì)應(yīng)于x(n)的離散傅里葉變換X(k)的虛部。舉例：已知一個(gè)定義在主值區(qū)間的實(shí)序列xo

27、nes(1，4)，zeros(1，4)，試將其分解成為偶對(duì)稱序列和奇對(duì)稱序列，并求它們的DFT，驗(yàn)證離散傅里葉變換的循環(huán)對(duì)稱性。MATLAB代碼：x=ones(1,5),zeros(1,5) %建立x(n)序列N=length(x);n=0:N-1;k=0:N-1;xr=x(mod(-n,N)+1); %求x(-n)xe=0.5*(x+xr) %求x(n)的偶序列xo=0.5*(x-xr) %求x(n)的奇序列X=x*(exp(-j*2*pi/N).(n'*k); %由x(n)的DFT求X(k)Xe=xe*(exp(-j*2*pi/N).(n'*k); %由xe(n)的DFT求

28、Xe(k)Xo=xo*(exp(-j*2*pi/N).(n'*k); %由xo(n)的DFT求Xo(k)error1=(max(abs(real(X)-Xe) %計(jì)算X(k)的實(shí)部與Xe(k)的差值error2=(max(abs(j*imag(X)-Xo) %計(jì)算X(k)的虛部與Xo(k)的差值subplot(2,4,1),stem(n,x);title('x(n)')subplot(2,4,2),stem(n,xr);title('x(-n)')subplot(2,4,3),stem(n,xe);title('xe(n)')subplo

29、t(2,4,4),stem(n,xo);title('xo(n)')subplot(2,4,5),stem(n,real(X);title('X(k)的實(shí)部')subplot(2,4,6),stem(n,imag(X);title('X(k)的虛部') subplot(2,4,7),stem(n,Xe);title('Xe(k)=DFT(xe(n)')subplot(2,4,8),stem(n,Xo);title('Xo(k)=DFT(xo(n)')x = 1 1 1 1 1 0 0 0 0 0xe = 1.000

30、0 0.5000 0.5000 0.5000 0.5000 0 0.5000 0.5000 0.5000 0.5000xo = 0 0.5000 0.5000 0.5000 0.5000 0 -0.5000 -0.5000 -0.5000 -0.5000error1 = 3.8778e-15error2 = 3.9475e-15(b) 如何利用離散傅立葉變換對(duì)信號(hào)進(jìn)行譜分析；譜分析的誤差及改進(jìn)方法；DFT的高分辨率頻譜與高密度頻譜之間的區(qū)別。利用離散傅立葉變換對(duì)信號(hào)進(jìn)行譜分析：時(shí)域連續(xù)信號(hào)離散傅里葉分析的基本步驟如下圖所示：其中低通濾波器LPF（預(yù)濾波器）的引入，是為了消除或減少時(shí)域連續(xù)信號(hào)轉(zhuǎn)

31、換成序列時(shí)，可能出現(xiàn)的頻譜混疊的影響。實(shí)際工作中，時(shí)域離散信號(hào)x(n)的時(shí)寬是很長(zhǎng)的甚至是無限長(zhǎng)。由于DFT的需要，必須把x(n)限制在一定的時(shí)間間隔之內(nèi)，即進(jìn)行數(shù)據(jù)截?cái)?。?shù)據(jù)的截?cái)嘞喈?dāng)于加窗處理。因此，在計(jì)算x(n)的DFT之前，用一個(gè)時(shí)域有限的窗函數(shù)w(n)加到x(n)上。通過A/D變換器轉(zhuǎn)換成采樣序列x(n)時(shí)域離散序列。其頻譜用表示，它是頻率的周期函數(shù)，即 由于進(jìn)行DFT的需要，對(duì)序列x(n)進(jìn)行加窗處理，即v(n)=x(n).w(n)。加窗對(duì)頻域的影響，用周期卷積表示。 最后進(jìn)行DFT運(yùn)算。加窗后的DFT是 例子：先用MATLAB產(chǎn)生數(shù)字信號(hào)：然后逐個(gè)用DFT進(jìn)行譜分析，分別取DFT

32、的長(zhǎng)度N16，32，畫出信號(hào)的幅譜圖，分析實(shí)驗(yàn)結(jié)果。MATLAB程序：例2：先用MATLAB產(chǎn)生出如下模擬信號(hào)：設(shè)采樣頻率Hz，然后用DFT進(jìn)行譜分析，分別取DFT的長(zhǎng)度N16，32，64，畫出信號(hào)的幅譜圖，橫軸打印模擬頻率。討論DFT長(zhǎng)度與頻率分辨率的關(guān)系，要在頻譜圖上分離出上述3個(gè)頻率，DFT長(zhǎng)度至少取多大？代碼：N1=16;N2=32;N3=64;n=0:N1-1;Fs=64;t=n/Fs;x4=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t);figure,subplot(231),stem(n,x4);gridX4=abs(fft(x4,N1);subplo

33、t(232),stem(n,X4);gridn=0:N2-1;Fs=64;t=n/Fs;x4=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t);subplot(233),stem(n,x4);gridX4=abs(fft(x4,N2);subplot(234),stem(n,X4);gridn=0:N3-1;Fs=64;t=n/Fs;x4=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t);subplot(235),stem(n,x4);gridaxis(0 31 0 4);X4=abs(fft(x4,N3);subplot(236),stem

34、(n,X4);gridaxis(0 31 0 40);運(yùn)行結(jié)果：可見，只有當(dāng)N=64時(shí)，3個(gè)頻率在頻譜圖上被分離出來。分析： 譜分析的誤差及改進(jìn)方法：（2）柵欄效應(yīng)（3）截?cái)嘈?yīng)DFT的高分辨率頻譜與高密度頻譜之間的區(qū)別:高分辨率譜通過增加信號(hào)的記錄長(zhǎng)度得到，頻率分辨率F與信號(hào)的實(shí)際長(zhǎng)度成反比。信號(hào)越長(zhǎng)，F(xiàn)值越小，即分辨率（頻域取樣中兩相鄰點(diǎn)間的頻率間隔）越高。高密度譜是指當(dāng)信號(hào)的時(shí)間域長(zhǎng)度不變時(shí)，在頻域內(nèi)對(duì)它的頻譜進(jìn)行提高采樣頻率而得到的高密度譜，它只可以更細(xì)化當(dāng)前分辨率下的頻譜，克服柵欄效應(yīng)，但不能改變DFT的分辨率。另外，采用尾部補(bǔ)零的方法不能提高DFT的高分辨率。即增大M,在原序列后面

35、加0，采樣頻率不變。區(qū)別1、高密度譜呈許多譜線型,補(bǔ)充0越多,譜線也越密集，但不改變頻譜分辨率。2、通過增加采樣點(diǎn)數(shù)N，可提高物理分辨率，得到高分辨率譜。實(shí)例分析：用兩個(gè)不同頻率的正弦信號(hào)開研究分辨率的大小，能分辨的的兩個(gè)正弦信號(hào)頻率越靠近，表明分辨率越高。設(shè)序列，序列有兩個(gè)主要頻率0.48和0.52。畫出其10點(diǎn)DFT、64點(diǎn)DFT及10點(diǎn)后序列補(bǔ)零至64點(diǎn)的DFT。10點(diǎn)DFT程序如下：clearn=0:9;xn=cos(0.48*pi*n)+cos(0.52*pi*n);xn1=xn(1:1:10);N=10;n=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);n

36、k=n'*k;WNnk=WN.nk;Xk=xn1*WNnk;subplot(121)stem(n,xn1);xlabel('n');ylabel('x(n)');axis(0 9 -2.5 2.5);title('正弦序列');n1=0:1:9;subplot(122)plot(n1,abs(Xk(1:1:10)title('幅度譜');xlabel('k');ylabel('|X(k)|'); 10點(diǎn)DFT正弦序列及頻譜圖64點(diǎn)DFT程序：clearn=0:1:63;xn=cos(0.48

37、*pi*n)+cos(0.52*pi*n);xn1=xn(1:1:64);N=64;n=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.nk;Xk=xn1*WNnk;subplot(121)stem(n,xn1);xlabel('n');ylabel('x(n)');axis(0 63 -2.5 2.5);title('正弦序列');n1=0:1:63;subplot(122)plot(n1,abs(Xk(1:1:64)title('幅度譜');xlabel('

38、;k');ylabel('|X(k)|');64點(diǎn)DFT正弦序列及頻譜圖10點(diǎn)序列補(bǔ)零值64點(diǎn)DFT的MATLAB程序：clearn=0:1:9;xn=cos(0.48*pi*n)+cos(0.52*pi*n);xn1=xn(1:1:10),zeros(1,54);N=64;n=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.nk;Xk=xn1*WNnk;subplot(121)stem(n,xn1);xlabel('n');ylabel('x(n)');axis(0 63

39、-2.5 2.5);title('正弦序列');n1=0:1:63;subplot(122)plot(n1,abs(Xk(1:1:64)title('幅度譜');xlabel('k');ylabel('|X(k)|');10點(diǎn)序列補(bǔ)零值64點(diǎn)DFT正弦序列及頻譜圖 :DFT的高分辨率頻譜與高密度頻譜之間的區(qū)別:高分辨率譜通過增加信號(hào)的記錄長(zhǎng)度得到，頻率分辨率F與信號(hào)的實(shí)際長(zhǎng)度成反比。信號(hào)越長(zhǎng)，F(xiàn)值越小，即分辨率（頻域取樣中兩相鄰點(diǎn)間的頻率間隔）越高。高密度譜是指當(dāng)信號(hào)的時(shí)間域長(zhǎng)度不變時(shí)，在頻域內(nèi)對(duì)它的頻譜進(jìn)行提高采樣頻率而得到的高

40、密度譜，它只可以更細(xì)化當(dāng)前分辨率下的頻譜，克服柵欄效應(yīng)，但不能改變DFT的分辨率。另外，采用尾部補(bǔ)零的方法不能提高DFT的高分辨率。即增大M,在原序列后面加0，采樣頻率不變。區(qū)別1、高密度譜呈許多譜線型,補(bǔ)充0越多,譜線也越密集，但不改變頻譜分辨率。2、通過增加采樣點(diǎn)數(shù)N，可提高物理分辨率，得到高分辨率譜。已知序列： x(n)=cos(0.82n)+2sin(0.43n) ，試確定一合適樣本數(shù)N，利用MATLAB的fft函數(shù)分析計(jì)算信號(hào)()xn的頻譜。 解：序列x(n)是一個(gè)周期序列。為了說明高密度譜和高分辨率譜之間的區(qū)別，分以下幾種情況進(jìn)行討論： 先取x(n)的前10個(gè)樣本，10點(diǎn)DFT的M

41、ATLAB程序如下：n=0:1:9;x = cos(0.82*pi*n)+2*sin(0.43*pi*n); subplot(2,1,1); stem(n,x,'.'); title('x(n), 0 <= n <= 9'); xlabel('n'); ylabel('x(n)'); axis(0,10,-2.5,2.5); Xk = fft(x); magXk = abs(Xk(1:1:6); k1 = 0:1:5; w1 = 2*pi/10*k1; subplot(2,1,2); stem(w1/pi,magXk,

42、'.'); title('Samples of DTFT Magnitude'); xlabel('frequency in pi units'); ylabel('|X(k)|'); axis(0,1,0,10); 由于樣本數(shù)不足，難以獲得足夠的信息而得到正確的結(jié)論。即從頻譜圖無法觀測(cè)到原復(fù)合余弦信號(hào)x(n)的w=0.43和w=0.82兩個(gè)頻率分量. 在先前()xn的前10個(gè)樣本后補(bǔ)90個(gè)零，以期得到一個(gè)更高密度的頻譜。補(bǔ)零后程序：n = 0:1:9; x = cos(0.82*pi*n)+2*sin(0.43*pi*n); n

43、1 = 0:1:99; x1 = x zeros(1,90); subplot(2,1,1); stem(n1,x1,'.'); title('x(n), 0 <= n <= 9 + 90 zeros'); xlabel('n'); ylabel('x(n)'); axis(0,100,-2.5,2.5); Xk = fft(x1); magXk = abs(Xk(1:1:51); k1 = 0:1:50; w1 = 2*pi/100*k1; 結(jié)果表明，序列在w=0.4處存在一個(gè)主頻，而原序列并未提供這一信息（原序列僅

44、含w=0.43和w=0.82兩個(gè)頻率分量）。補(bǔ)零運(yùn)算只是獲得一個(gè)更高密度的頻譜，改善了柵欄現(xiàn)象而已為了獲得更多的頻譜信息，現(xiàn)采用()xn的前120個(gè)樣本，以得到一個(gè)高分辨率頻譜。120點(diǎn)DFT的MATLAB程序如下subplot(2,1,2); stem(w1/pi,magXk,'.'); title('Samples of DTFT Magnitude'); xlabel('frequency in pi units'); ylabel('|X(k)|'); axis(0,1,0,10);n = 0:1:119; x = cos

45、(0.82*pi*n)+2*sin(0.43*pi*n); subplot(2,1,1); stem(n,x,'.'); title('x(n), 0 <= n <= 119'); xlabel('n'); ylabel('x(n)'); axis(0,120,-2.5,2.5); Xk = fft(x); magXk = abs(Xk(1:1:61); k1 = 0:1:60; w1 = 2*pi/120*k1; subplot(2,1,2); stem(w1/pi,magXk,'.'); title('Samples of DTFT Magnitude'); xlabel('frequency in pi units'); ylabel('|X(k)|'); axis(0,1,0,60); 現(xiàn)在可從頻譜圖上清楚地觀測(cè)到w=0.43和w=0.82兩個(gè)頻率分量，這才是 x(n)的高分辨率頻譜。(c) 線性卷積和圓周卷積的關(guān)系；如何利用圓周卷積求系統(tǒng)響應(yīng)。L點(diǎn)圓周卷積等于線性卷積以L為周期的周期延拓序列