數(shù)字信號(hào)處理實(shí)驗(yàn)指導(dǎo)書(shū)

1、實(shí)驗(yàn)一 離散系統(tǒng)的時(shí)域分析一. 實(shí)驗(yàn)?zāi)康模?）掌握求系統(tǒng)響應(yīng)的方法。（2）掌握時(shí)域離散系統(tǒng)的時(shí)域特性。二.實(shí)驗(yàn)原理與方法在時(shí)域中，描寫(xiě)系統(tǒng)特性的方法是差分方程和單位脈沖響應(yīng)，在頻域可以用系統(tǒng)函數(shù)描述系統(tǒng)特性。已知輸入信號(hào)可以由差分方程、單位脈沖響應(yīng)或系統(tǒng)函數(shù)求出系統(tǒng)對(duì)于該輸入信號(hào)的響應(yīng)，本實(shí)驗(yàn)僅在時(shí)域求解。在計(jì)算機(jī)上適合用遞推法求差分方程的解，最簡(jiǎn)單的方法是采用MATLAB語(yǔ)言的工具箱函數(shù)filter函數(shù)。也可以用MATLAB語(yǔ)言的工具箱函數(shù)conv函數(shù)計(jì)算輸入信號(hào)和系統(tǒng)的單位脈沖響應(yīng)的線性卷積，求出系統(tǒng)的響應(yīng)。三實(shí)驗(yàn)內(nèi)容及結(jié)果（1）編制程序，包括產(chǎn)生輸入信號(hào)、單位脈沖響應(yīng)序列的子程序，用fi

2、lter函數(shù)或conv函數(shù)求解系統(tǒng)輸出響應(yīng)的主程序。程序中要有繪制信號(hào)波形的功能。給定一個(gè)低通濾波器的差分方程為： 輸入信號(hào)： a) 分別求出系統(tǒng)對(duì)兩個(gè)輸入信號(hào)的響應(yīng)序列，并畫(huà)出其波形。 b) 求出系統(tǒng)的單位沖響應(yīng)，畫(huà)出其波形。% A=1,-0.9;% B=0.05,0.05;% x1n=1 1 1 1 1 1 1 1 zeros(1,25);% x2n=ones(1,32);% y1n=impz(B,A,33);% n=0:length(y1n)-1;% figure(1);% stem(n,y1n);% title('1102303005 y1(n)');% box on;

3、% xlabel('n');% ylabel('y1(n)');% x2n=ones(1,33);% y2n=filter(B,A,x2n);% n=0:length(y2n)-1;% figure(2);% stem(n,y2n);% title('1102303005 y2(n)');% box on;% xlabel('n');% ylabel('y2(n)');% hn=filter(B,A,x1n);% n=0:length(hn)-1;% figure(3);% stem(n,hn);% title(&

4、#39;1102303005 h(n)');% box on% xlabel('n');%ylabel('h(n)');（2）給定系統(tǒng)的單位脈沖響應(yīng)為：用線性卷積法分別求兩個(gè)系統(tǒng)對(duì)輸入信號(hào)的輸出響應(yīng)，并畫(huà)出波形。% x1n=1 1 1 1 1 1 1 1 zeros(1,24);% h1n=1 1 1 1 1 1 1 1 1 1 zeros(1,22);% h2n=1 2.5 2.5 1 zeros(1,28);% y21n=conv(h1n,x1n);% y22n=conv(h2n,x1n);% n=0:length(y21n)-1;% figure(

5、4);% stem(n,y22n);% title('1102303005 x1(n)');% box on;% xlabel('n');% ylabel('x1(n)');% n=0:length(y21n)-1;% figure(5);% stem(n,y21n);% title('1102303005 y21(n)');% box on;% xlabel('n');% ylabel('y21(n)');% n=0:length(y22n)-1;% figure(6);% stem(n,y22n

6、);% title('1102303005 y22(n)');% box on;% xlabel('n');% ylabel('y22(n)');（3）測(cè)量人耳辨別回聲的最小時(shí)間設(shè)計(jì)一個(gè)混合聲音數(shù)字系統(tǒng)，利用MATLAB語(yǔ)言實(shí)現(xiàn)聲音的混合。A=1,-0.9;B=0.05,0.05;x1n=1 1 1 1 1 1 1 1 zeros(1,25);x2n=ones(1,32);xn=1 zeros(1,32);y1n=filter(B,A,33);h=filter(B,A,xn);load mtlb.matyn=conv(h,mtlb)n=0:len

7、gth(yn)-1;figure(1);stem(n,yn);title('1102303005 yn)');box on;xlabel('n');ylabel('yn');四思考題（1）簡(jiǎn)述在時(shí)域求系統(tǒng)響應(yīng)的方法。有兩種，一種是通過(guò)解差分方程求得系統(tǒng)的輸出，注意要合理的選擇初始條件，一種是已知系統(tǒng)的單位脈沖響應(yīng)，通過(guò)球輸出信號(hào)和系統(tǒng)的單位脈沖的線性卷積求得性輸出。（2）簡(jiǎn)述判斷系統(tǒng)穩(wěn)定性的方法。如果信號(hào)經(jīng)過(guò)低通濾波器，把信號(hào)的高頻部分濾去，時(shí)域信號(hào)的的劇烈變化將被平滑。實(shí)驗(yàn)二 離散系統(tǒng)的頻域分析一、實(shí)驗(yàn)?zāi)康?（1）熟悉對(duì)離散系統(tǒng)的頻率響應(yīng)分析方

8、法。 （2）加深對(duì)零、極點(diǎn)分布的概念理解。二、實(shí)驗(yàn)原理 離散系統(tǒng)的時(shí)域方程為： 其變換域分析方法如下： 頻域: 系統(tǒng)的頻率響應(yīng)為: Z域: 系統(tǒng)函數(shù)為: 分解因式: ， 其中和稱(chēng)為零、極點(diǎn)。三、實(shí)驗(yàn)內(nèi)容及結(jié)果 （1）求差分方程 所對(duì)應(yīng)的系統(tǒng)的頻率響應(yīng)。程序：k=256;%采樣頻率為256Hznum=0.8 -0.44 0.36 0.02;den=1 0.7 -0.45 -0.6;w=0:pi/k:pi;h=freqz(num,den,w);subplot(2,2,1);plot(w/pi,real(h);grid;title(1102303005 實(shí)部)；xlabel('omega/p

9、i');ylabel('幅度');subplot(2,2,2);plot(w/pi,imag(h);grid;title('1102303005 虛部)；xlabel('omega/pi');ylabel('Amplitude');subplot(2,2,3);plot(w/pi,abs(h);grid;title(1102303005 幅度譜)；xlabel('omega/pi');ylabel('幅值');subplot(2,2,4);plot(w/pi,angle(h);grid;title(

10、1102303005 相位譜)；xlabel('omega/pi');ylabel('弧度'); 結(jié)果：（2）求系統(tǒng)的零、極點(diǎn)和幅度頻率響應(yīng)和相位響應(yīng)。程序：clc;clear ;num=0.0528 0.0797 0.1295 0.1295 0.0797 0.0528;den=1 -1.8107 2.4947 -1.8801 0.9537 -0.2336;figure(1)zplane(num,den);grid;z,p,k=tf2zp(num,den);disp('零點(diǎn)=');disp(z);disp('極點(diǎn)=');disp(

11、p);disp('增益系數(shù)=');disp(k);k=256; %采樣頻率為K=256Hzw=0:pi/k:pi;h=freqz(num,den,w);figure(2)subplot(2,2,1);plot(w/pi,real(h);grid;title('實(shí)部');xlabel('omega/pi');ylabel('幅度');subplot(2,2,2);plot(w/pi,imag(h);grid;title('虛部');xlabel('omega/pi');ylabel('Ampl

12、itude');subplot(2,2,3);plot(w/pi,abs(h);grid;title('幅度譜');xlabel('omega/pi');ylabel('幅值');subplot(2,2,4);plot(w/pi,angle(h);grid;title('相位譜');xlabel('omega/pi');ylabel('弧度'); 結(jié)果：的零、極點(diǎn)和幅度頻率響應(yīng)和相位響應(yīng)。實(shí)驗(yàn)三 FFT算法的應(yīng)用一、實(shí)驗(yàn)?zāi)康?（1）加深對(duì)離散信號(hào)的DFT的理解。 （2）熟悉FFT算法的應(yīng)用。

13、二、實(shí)驗(yàn)原理 N點(diǎn)序列的DFT和IDFT變換定義式如下： , 利用旋轉(zhuǎn)因子具有周期性，可以得到快速算法（FFT）。 在MATLAB中，可以用函數(shù)和 計(jì)算N點(diǎn)序列的DFT正、反變換。三、實(shí)驗(yàn)內(nèi)容及結(jié)果 （1）觀察一句話的頻譜特點(diǎn)：將自己錄制的一段話均勻分段，并用DFT計(jì)算每段信號(hào)的頻譜。程序：Fs=8000;T=4;x=wavrecord(T*Fs,Fs);L=length(x);t=(1:L)/Fs;subplot(411);plot(t,x);xlabel('t/s');ylabel('x(t)');title('1102303005)；axis(0

14、4 -0.7 0.7);subplot(413);specgram(x,256,Fs,256,250);xlabel('t/s');ylabel('f/Hz');title('1102303005)；結(jié)果： （2）2N點(diǎn)實(shí)數(shù)序列 N=64。一次算出，并繪出 的圖形。程序：N=64;n=0:2*N-1;x1=cos(2*pi*7*n/N);x2=0.5*cos(2*pi*19*n/N);x=x1+x2;subplot(2,1,1);stem(n,x);xlabel('t/T');ylabel('x(n)');title(&#

15、39;1102303005)；xk=fft(x,N);disp(xk);xk=abs(xk);k=0:N-1;subplot(2,1,2);stem(k,xk);xlabel('k');ylabel('X(k)');title('1102303005)；結(jié)果： （3）已知某序列在單位圓上的N=64等分樣點(diǎn)的Z變換為:。用N點(diǎn)IFFT程序計(jì)算出。程序：N=64;k=0:N-1;xk=1./(1-0.8*exp(2*pi*k/N)*-1*j);subplot(2,1,1)stem(k,xk)xlabel('k');ylabel('X(

16、k)');title('1102303005)；axis(0 63 0 6)x=ifft(xk,N);subplot(2,1,2)stem(k,x)xlabel('t/T');ylabel('x(n)');axis(0 63 0 1.5)title('1102303005)；結(jié)果：四、思考題（1）如果周期信號(hào)的周期預(yù)先不知道，如何用FFT進(jìn)行譜分析？ 答：周期信號(hào)的周期預(yù)先不知道時(shí),可先截取M點(diǎn)進(jìn)行DFT,再將截取長(zhǎng)度擴(kuò)大1倍截取,比較結(jié)果,如果二者的差別滿足分析誤差要求,則可以近似表示該信號(hào)的頻譜,如果不滿足誤差要求就繼續(xù)將截取長(zhǎng)度加倍

17、,重復(fù)比較,直到結(jié)果滿足要求。（2）如何選擇FFT的變換區(qū)間？（包括非周期信號(hào)和周期信號(hào)） 答：一、對(duì)于非周期信號(hào)：有頻譜分辨率F，而頻譜分辨率直接和FFT的變換區(qū)間有關(guān)，因?yàn)镕FT能夠?qū)崿F(xiàn)的頻率分辨率是2/N.因此有最小的N>2/F。就可以根據(jù)此式選擇FFT的變換區(qū)間。 二、對(duì)于周期信號(hào)，周期信號(hào)的頻譜是離散譜，只有用整數(shù)倍周期的長(zhǎng)度作FFT，得到的離散譜才能代表周期信號(hào)的頻譜。實(shí)驗(yàn)四 IIR數(shù)字濾波器的設(shè)計(jì)1、實(shí)驗(yàn)?zāi)康模?）掌握雙線性變換法及脈沖相應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器的具體設(shè)計(jì)方法;（2）熟悉用雙線性變換法及脈沖響應(yīng)不變法設(shè)計(jì)低通、高通和帶通IIR數(shù)字濾波器的計(jì)算機(jī)編程。2、

18、實(shí)驗(yàn)原理在MATLAB中，可以用下列函數(shù)輔助設(shè)計(jì)IIR數(shù)字濾波器：（1）利用buttord和cheb1ord確定低通原型巴特沃斯和切比雪夫?yàn)V波器的階數(shù)和截止頻率；如：求階數(shù)N, Wn = cheb1ord (Wp, Ws, Rp, Rs,s)選擇項(xiàng)說(shuō)明：high類(lèi)別。缺省為low； s模/數(shù)，缺省為數(shù)Rp即 ap, Rs即as ； WnChebyshev自然頻率（3dB頻率），數(shù)字設(shè)計(jì)： Wpwp/p Wsws/p 。（2）num,den=butter（N,Wn）(巴特沃斯) num,den=cheby1（N,Wn）,num,den=cheby2（N,Wn）(切比雪夫1型和2型)（3）lp2h

19、p，lp2bp，lp2bs可以完成低通濾波器到高通、帶通、帶阻濾波器的轉(zhuǎn)換；（4）使用bilinear可以對(duì)模擬濾波器進(jìn)行雙線性變換，求得數(shù)字濾波器的傳輸函數(shù)系數(shù)；（5）利用impinvar可以完成脈沖響應(yīng)不變法的模擬濾波器到數(shù)字濾波器的轉(zhuǎn)換。3、實(shí)驗(yàn)內(nèi)容利用MATLAB編程，用脈沖響應(yīng)不變法和雙線性變換法設(shè)計(jì)一個(gè)數(shù)字低通濾波器，指標(biāo)要求如下：通帶邊緣頻率：，通帶峰值起伏：； 阻帶邊緣頻率：，最小阻帶衰減：。畫(huà)出兩種設(shè)計(jì)方法下的損耗函數(shù)和幅頻特性。%脈沖響應(yīng)不變法Fs=400;T=1/Fs;wp=0.3*pi;ws=0.8*pi;Wp=wp*Fs;Ws=ws*Fs;Ap=1;As=40;N,W

20、c=buttord(Wp,Ws,Ap,As,'s');B,A=butter(N,Wc,'s');D,C=impinvar(B,A,Fs);Hz=freqz(D,C,W);f=(W*Fs)/(2*pi);figure(1)plot(f,abs(Hz);grid on;title('1102303005脈沖響應(yīng)不變法、幅頻特性');xlabel('頻率/Hz');ylabel('歸一化幅值');figure(2)a=20*log10(abs(Hz);plot(f,a);grid on;title('110230

21、3005脈沖響應(yīng)不變法、衰減函數(shù)')xlabel('頻率/Hz');ylabel('幅度/dB');%雙線性變化Fs=400;T=1/Fs;wp=0.3*pi;ws=0.8*pi;Wp=(2/T)*tan(wp/2);Ws=(2/T)*tan(ws/2);Ap=1;As=40;N,Wc=buttord(Wp,Ws,Ap,As,'s');B,A=butter(N,Wc,'s');D,C=bilinear(B,A,Fs);Hz,W=freqz(D,C);f=(W*Fs)/(2*pi);figure(1)plot(f,abs(H

22、z);grid on;title('1102303005雙線性變換法、幅頻特性');xlabel('頻率/Hz');ylabel('歸一化幅值');figure(2)a=20*log10(abs(Hz);plot(W/pi,a);grid on;title('1102303005雙線性變換法、幅頻特性');xlabel('頻率/Hz');ylabel('幅度/dB');實(shí)驗(yàn)五 FIR數(shù)字濾波器的設(shè)計(jì)1、實(shí)驗(yàn)?zāi)康模?（1）加深對(duì)數(shù)字濾波器的常用指標(biāo)理解。 （2）學(xué)習(xí)數(shù)字濾波器的設(shè)計(jì)方法。 2、實(shí)驗(yàn)原理

23、：（1）低通濾波器的常用指標(biāo)：（a）通帶邊緣頻率；（b）阻帶邊緣頻率；（c）通帶起伏；（d）通帶峰值起伏，（e）阻帶起伏，最小阻帶衰減。（2）在MATLAB中，熟悉函數(shù)fir1、kaiserord 、remezord、remez的使用； B = fir1(n,Wn,'high','noscale')設(shè)計(jì)濾波器；n,Wn,beta,ftype = kaiserord(f,a,dev)估計(jì)濾波器階數(shù)；（默認(rèn)情況下，過(guò)濾器的第一通帶中心為精確的一個(gè)窗口。使用一個(gè)尾隨'noscale'參數(shù)來(lái)防止這種縮放remezord、remez一起使用設(shè)計(jì)等波紋線性相

24、位FIR濾波器）3、實(shí)驗(yàn)內(nèi)容：利用MATLAB編程,用窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器，指標(biāo)要求如下：通帶邊緣頻率：，通帶峰值起伏：。阻帶邊緣頻率：，最小阻帶衰減：。clear allclcwp1=0.45*pi;wp2=0.65*pi;ws1=0.3*pi;ws2=0.75*pi;width1 =wp1-ws1;width2 =ws2-wp2;width=min(width1,width2)N=ceil(6.62*pi/width)n =0: (N-1);a =(N-1)/2+0.00001;wc1=(ws1+wp1)/2;wc2=(ws2+wp2)/2;hd=(sin(wc2*(n-a)-si

25、n(wc1*(n-a)./(pi*(n-a);win =(hanning(N)'h=hd.*win;freqz(h);subplot(2,1,1)title('1102303005 幅頻特性')xlabel('頻率')ylabel('db')subplot(2,1,2)title('1102303005 相頻特性')xlabel('頻率')ylabel('db')實(shí)驗(yàn)六 數(shù)字信號(hào)處理在雙音多頻撥號(hào)系統(tǒng)中的應(yīng)用1、 DTMF的組成在電話中，數(shù)字09的中每一個(gè)都用兩個(gè)不同的單音頻傳輸，所用的8個(gè)

26、頻率分成高頻帶和低頻帶兩組，低頻帶有四個(gè)頻率：679Hz,770Hz,852Hz和941Hz；高頻帶也有四個(gè)頻率：1209Hz,1336Hz,1477Hz和1633Hz.。每一個(gè)數(shù)字均由高、低頻帶中各一個(gè)頻率構(gòu)成，例如1用697Hz和1209Hz兩個(gè)頻率，信號(hào)用表示，其中。這樣8個(gè)頻率形成16種不同的雙頻信號(hào)。2、 電話中的雙音多頻（DTMF）信號(hào)的產(chǎn)生與檢測(cè)（1）雙音多頻信號(hào)的產(chǎn)生假設(shè)時(shí)間連續(xù)的DTMF信號(hào)用表示，式中f1和f2是按照上表選擇的兩個(gè)頻率，f1代表低頻帶中的一個(gè)頻率，f2代表高頻帶中的一個(gè)頻率。采用數(shù)字方法產(chǎn)生DTMF信號(hào)，規(guī)定用8KHz對(duì)DTMF信號(hào)進(jìn)行采樣，得到時(shí)域離散信號(hào)

27、：因?yàn)椴蓸宇l率是8000Hz，因此要求每125ms輸出一個(gè)樣本，得到的序列再送到D/A變換器和平滑濾波器，輸出便是連續(xù)時(shí)間的DTMF信號(hào)。DTMF信號(hào)通過(guò)電話線路送到交換機(jī)。（2）雙音多頻信號(hào)的檢測(cè)在接收端，要對(duì)收到的雙音多頻信號(hào)進(jìn)行檢測(cè)，檢測(cè)兩個(gè)正弦波的頻率是多少，以判斷所對(duì)應(yīng)的十進(jìn)制數(shù)字或者符號(hào)。顯然這里仍然要用數(shù)字方法進(jìn)行檢測(cè)，因此要將收到的時(shí)間連續(xù) DTMF信號(hào)經(jīng)過(guò)A/D變換，變成數(shù)字信號(hào)進(jìn)行檢測(cè)。檢測(cè)的方法有兩種，一種是用一組濾波器提取所關(guān)心的頻率，根據(jù)有輸出信號(hào)的2個(gè)濾波器判斷相應(yīng)的數(shù)字或符號(hào)。另一種是用DFT（FFT）對(duì)雙音多頻信號(hào)進(jìn)行頻譜分析，由信號(hào)的幅度譜，判斷信號(hào)

28、的兩個(gè)頻率，最后確定相應(yīng)的數(shù)字或符號(hào)。3、 根據(jù)實(shí)驗(yàn)原理利用數(shù)字信號(hào)處理方法進(jìn)行雙音多頻信號(hào)的產(chǎn)生與解調(diào)。clear all;clc; tm=1,2,3,65;4,5,6,66;7,8,9,67;42,0,35,68;N=205;K=18,20,22,24,31,34,38,42; f1=697,770,852,941;f2=1209,1336,1477,1633;TN=input('鍵入10學(xué)號(hào)=');TNr=0;for l=1:10;d=fix(TN/10(10-l);TN=TN-d*10(10-l); for p=1:4; for q=1:4; if tm(p,q)=ab

29、s(d); break,end end if tm(p,q)=abs(d); break,end end n=0:1023; x = sin(2*pi*n*f1(p)/8000) + sin(2*pi*n*f2(q)/8000); sound(x,8000); pause(0.1) X=goertzel(x(1:205),K+1); val = abs(X); subplot(5,2,l); stem(K,val,'.');grid;title('1102303005');xlabel('k');ylabel('|H(k)|');axis(10 50 0 120) limit = 80;