數(shù)字信號(hào)處理實(shí)驗(yàn)報(bào)告

《數(shù)字信號(hào)處理實(shí)驗(yàn)報(bào)告》由會(huì)員分享，可在線閱讀，更多相關(guān)《數(shù)字信號(hào)處理實(shí)驗(yàn)報(bào)告（39頁(yè)珍藏版）》請(qǐng)?jiān)谘b配圖網(wǎng)上搜索。

1、數(shù)字信號(hào)處理實(shí)驗(yàn)報(bào)告（一） 實(shí)驗(yàn)一、 <1>A=[1.2,3,5.1,-1.5,2.4],n=[-2,-1,0,1,2];B=0.5m,m=0~10; 請(qǐng)編程完成下面操作并畫圖表示結(jié)果 AB，A-B，A.*B，A（n-2） 解答： 1、AB 這是矩陣的叉乘，所以要保證矩陣乘法的規(guī)則。即A的列數(shù)等于B的行數(shù)。故A要補(bǔ)零A2=[1.2,3,5.1,-1.5,2.4];---------A3=[1.2,3,5.1,-1.5,2.4,0,0,0,0,0,0]; 同時(shí)要轉(zhuǎn)置成列矩陣。A1=A3; C=A1*B; 2、 A-B D=A3-B; A和B要保持長(zhǎng)度一致，所

2、以用補(bǔ)零的A3來(lái)和B作差。 3、 A.*B E=A3.*B; 點(diǎn)乘要求長(zhǎng)度一致，計(jì)算時(shí)對(duì)應(yīng)的元素相乘。 4、 A（n-2） A右移兩位 F=[0,0,A2(1:end)]; 程序： %1、A=[1.2,3,5.1,-1.5,2.4],n=[-2,-1,0,1,2];B=0.5m,m=0~10; %請(qǐng)編程完成下面操作并畫圖表示結(jié)果 %AB，A-B，A.*B，A（n-2） m=0:10; % A1=[1.2;3;5.1;-1.5;2.4;0;0;0;0;0;0]; A2=[1.2,3,5.1,-1.5,2.4]; A3=[1.2,3,5

3、.1,-1.5,2.4,0,0,0,0,0,0];A1=A3; n1=[-2,-1,0,1,2]; B=0.5.^m; C=A1*B; D=A3-B; subplot(311); stem(m,D); title(A-B); E=A3.*B; subplot(312); stem(m,E); F=[0,0,A2(1:end)]; subplot(313); u=[-2,-1,0,1,2,3,4]; stem(u,F); <2>p138 2.28，2.34. 解答： 1、求頻率響應(yīng)的幅度相應(yīng)和相位響應(yīng) 可以直接用freqz來(lái)求，也可以是用公式計(jì)算

4、出來(lái) clear all; n=0:1:1000; h=0.9.^n; w=-pi:2*pi/N:pi; H=freqz(h,1,w);%繪制定制區(qū)域,則直接給定范圍,該范圍為給定自變量向量,例如 freqz(b,a,[-pi:2*pi/n:pi]) plot(w,abs(H)); 注意：用freqz做dtft時(shí)需要求特定范圍的dtft，這時(shí)可以在freqz中直接寫上給定的范圍，比如我們這次要求-pi到+pi范圍內(nèi)的dtft，所以在freqz函數(shù)中直接寫上范圍就可以了。 這是一種方法。 x=exp(

5、i*w)./(exp(i*w)-0.9*ones(1,1000)); y=exp(i*w)./(exp(i*w)-0.3*ones(1,1000))+exp(i*w)./(exp(i*w)-0.5*ones(1,1000)); 也可以像上面那樣直接算出來(lái) 。 2、 當(dāng)輸入為x（n）時(shí)，可以用conv和filter兩種方法。 Conv在計(jì)算卷積時(shí)，得到的序列的長(zhǎng)度是length MAX([LENGTH(A)+LENGTH(B)-1,LENGTH(A),LENGTH(B)]). 而filter卷積時(shí)得到的序列長(zhǎng)度是MAX(LENGTH(A),LENGTH(B))-1）； 可以

6、說(shuō)filter做卷積是conv卷積的截?cái)唷? 程序： w=linspace(-pi,pi,1000); % LINSPACE(X1, X2, N) generates N points between X1 and X2. x=exp(i*w)./(exp(i*w)-0.9*ones(1,1000)); %這是公式 ONES(M,N) or ONES([M,N]) is an M-by-N matrix of ones. % A./B denotes element-by-

7、element division. A and B subplot(4,2,1); plot(w/pi,abs(x)); xlabel(頻率f); title(幅度特性); ylabel(幅值) grid; %GRID ON adds major grid lines to the current axes. %GRID OFF removes major and minor grid lines from the current axes.

8、 subplot(4,2,2); plot(w/pi,angle(x)); grid; xlabel(頻率f); title(相位特性); ylabel(相位值) y=exp(i*w)./(exp(i*w)-0.3*ones(1,1000))+exp(i*w)./(exp(i*w)-0.5*ones(1,1000)); subplot(4,2,3); plot(w/pi,abs(y)); xlabel(頻率f); title(幅度特性);

9、 ylabel(幅值) grid; %GRID ON adds major grid lines to the current axes. %GRID OFF removes major and minor grid lines from the current axes. subplot(4,2,4); plot(w/pi,angle(y)); grid; xlabel(頻率f);

10、 title(相位特性); ylabel(相位值) n=0:1:50; h1=0.9.^n; h2=0.3.^n+0.5.^n; x1=3*cos(0.5*pi*n+0.25*pi)+2*sin(0.3*pi*n); y1=conv(x1,h1); subplot(4,2,5); stem(y1); grid; title(conv卷積h1); y2=conv(x1,h2); subplot(4,2,6); stem(y2); grid; title(conv卷積h2); y3=filter(h1,1,x1); subpl

11、ot(4,2,7); stem(y3); grid; title(filter卷積h1); y4=filter(h2,1,x1); subplot(4,2,8); stem(y4); grid; title(filter卷積h2); 解答：已知系統(tǒng)函數(shù)求輸出y（n），用x（n）卷積h（n）。 y1=filter(B1,A1,x); 在畫零極點(diǎn)圖是用zplane（B，A）； 程序： clear all; n=0:0.5:39; x=2*cos(pi*n/5); subplot(331);

12、 stem(n,x); B1=[1,1,0];A1=[1,-0.6]; y1=filter(B1,A1,x); subplot(332); stem(n,y1); subplot(333); zplane(B1,A1); subplot(334); stem(n,x); B2=[1,1,1];A2=[1,0.5,-0.25]; y2=filter(B2,A2,x); subplot(335); stem(n,y2); subplot(336); zplane(B2,A2); subplot(337); stem(n,x); B3=[1,0,1];A3=[1,

13、-6,9]; y3=filter(B3,A3,x); subplot(338); stem(n,y3); subplot(339); zplane(B3,A3); 結(jié)論： Filter卷積之后的長(zhǎng)度沒(méi)變，因?yàn)槌绦蛑衳和y在畫圖時(shí)用的是同一個(gè)n。 實(shí)驗(yàn)二 <1> 將DTFT和DFT畫在同一張圖上，這樣可以比較兩者之間的關(guān)系。 解答： 1、 用freqz做DTFT，用fft做DFT 畫在同一張圖上，注意兩個(gè)圖的x坐標(biāo)要一致，不一致的話沒(méi)有比較的意義。 N=100; x=[2,1,4,2,3]; [Y1,w]=freqz(x,1

14、,N,whole); %[H,W] = FREQZ(B,A,N) 因?yàn)閿?shù)字分析，所以我們只能取N盡量大，看起來(lái)H（e^j^w）是連續(xù)的。此時(shí)得到的w是（0-2π）； N1=5; n1=0:(2*pi/5):2*pi*(1-1/5); Y2=fft(x,N1); 這里也可以用freqz來(lái)求dft Y2=fft(x,N1); % Y2=freqz(x,1,N1,whole); 注意fft做的是0-2pi范圍內(nèi)的，所以freqz必須加‘whole’； Fft是有限長(zhǎng)的序列的傅里葉變換，所以N取序列的長(zhǎng)度5. 但是為了讓兩者對(duì)應(yīng)，所以讓兩者x坐標(biāo)保持一致，所以要加

15、上n1=0:(2*pi/5):2*pi*(1-1/5); stem(n1,Y2); 使得在畫圖時(shí)橫坐標(biāo)保持一致。 程序： clear all; N=100; x=[2,1,4,2,3]; [Y1,w]=freqz(x,1,N,whole); %[H,W] = FREQZ(B,A,N) subplot(221); plot(w,Y1); title(DTFT1); xlabel(n); ylabel(DTFT); hold on; N1=5; n1=0:(2*pi/5):2*pi*(1-1/5); Y2=fft(x,N1); % s

16、ubplot(222); stem(n1,Y2); title(DFT1); xlabel(n); ylabel(DFT); x1=[2,1,4,2,3,0,0,0]; [Y3,w]=freqz(x1,1,N,whole); subplot(212); plot(w,Y3); title(DTFT2); xlabel(n); ylabel(DTFT); hold on; N2=8; n2=0:(2*pi/8):2*pi*(1-1/8); Y4=fft(x1,N2); stem(n2,Y4); title(DFT2); xlabel(n); ylabel(D

17、FT); 2、 結(jié)論 由上圖可知，在x（n）尾部補(bǔ)零時(shí)，對(duì)DTFT沒(méi)有影響，因?yàn)樵谟?jì)算公式中 DTFT和n無(wú)關(guān)，但是DFT和k有關(guān)，而K的取值范圍是（0-N-1）；所以補(bǔ)零后改變了DFT的點(diǎn)數(shù)，從而改變了抽樣的點(diǎn)數(shù)。 同時(shí)還可以看出DFT是對(duì)DTFT的抽樣。 <2> 解答： 1、選擇合適的變換區(qū)間N，第一個(gè)由n為0~7，所以變換區(qū)間N=8.然后利用Y1=fft(x1,N1); magY1=abs(Y1); angY1=angle(Y1); stem(n,magY1); 畫出圖來(lái)即可。 2、 由n為-7~7，所以N選擇15.

18、 先將15個(gè)值算出來(lái)，然后再畫圖時(shí)stem(n,magY2); 程序： clear all; N1=8;n=0:1:N1-1; x1=3*cos(0.25*pi*n); Y1=fft(x1,N1); magY1=abs(Y1); angY1=angle(Y1); subplot(221); stem(n,magY1); title(x1幅頻特性); subplot(222); stem(n,angY1); title(x1相頻特性); clear all; N2=15;n2=-7:1:7; x2=0.8.^abs(

19、n2); Y2=fft(x2,N2); magY2=abs(Y2); angY2=angle(Y2); subplot(223); stem(n2,magY2); title(x2幅頻特性); subplot(224); stem(n2,angY2); % grid; title(x2相頻特性); <3> 解答： 1、n取0~10時(shí)，N為11。用fft算出DFT即可。 N=11;n=0:1:10; x=cos(0.48*pi*n)+cos(0.52*pi*n); % DFT Y1=fft(x,N); subplot(221);

20、 plot(n,Y1); 2、 補(bǔ)零用x1=[x,zeros(1,90)]; 注意此時(shí)的N變?yōu)?01 3、 不補(bǔ)零，直接取N=101 N2=101;n2=0:1:100; x2=cos(0.48*pi*n2)+cos(0.52*pi*n2); Y3=fft(x2,N2); subplot(223); plot(n2,abs(Y3)); 程序： clear all; N=11;n=0:1:10; x=cos(0.48*pi*n)+cos(0.52*pi*n); % DFT Y1=fft(x,N); subplot(221); plot(n,Y1);

21、%補(bǔ)零 90 % clear all; N1=101;n1=0:1:100; x1=[x,zeros(1,90)]; Y2=fft(x1,N1); subplot(222); plot(n1,Y2); % 高信號(hào)分辨率 但是并沒(méi)有區(qū)別 ？？？？？ % clear all; N2=101;n2=0:1:100; x2=cos(0.48*pi*n2)+cos(0.52*pi*n2); Y3=fft(x2,N2); subplot(223); plot(n2,abs(Y3)); 結(jié)論：在時(shí)域抽樣點(diǎn)數(shù)n不變的情況下，在有效數(shù)據(jù)后增加零值點(diǎn)可以減小柵欄效應(yīng)，

22、提高DFT的計(jì)算分辨率，但是并沒(méi)有提高頻率分辨率。只有增加有效數(shù)據(jù)的位數(shù)才可以提高頻率分辨率。 <1> 解答： 1、 由于題目本身給的就是級(jí)聯(lián)型的形式，所以我們只要將每一個(gè)因式歸一化后即可直接劃出級(jí)聯(lián)型的流圖。 2、直接型的流圖可以用sos2tf函數(shù)將級(jí)聯(lián)型轉(zhuǎn)化成直接型。 3、直接II型根據(jù)直接I型就可以畫出。 4、使用[C,B,A]=tf2par(6*b,a); 將直接型轉(zhuǎn)化成并聯(lián)型。 程序： clear all; % SOS = [ b01 b11 b21 1 a11 a21 %

23、 b02 b12 b22 1 a12 a22 % ... % b0L b1L b2L 1 a1L a2L ] sos=[1,0,1,1,-0.6,0.36;1,-1/3,0,1,-0.65,0;1,2,1,1,0,0.49]; [b,a] = sos2tf(sos); %級(jí)聯(lián)型到直接型的轉(zhuǎn)化 不需要考慮系統(tǒng)增益g %b為直接型多項(xiàng)式分子系數(shù)向量 a為直接型多項(xiàng)式分母系數(shù)向量 [C,B,A]=residuez(6*b,a);%直接型到并聯(lián)型的轉(zhuǎn)化

24、 % C為直接型多項(xiàng)式系數(shù)向量 B 為包含各個(gè)bi的K乘二維實(shí)系數(shù)矩陣 A為包含ai的k乘三維實(shí)系數(shù)矩陣 % b為直接型多項(xiàng)式分子系數(shù)向量 a為直接型多項(xiàng)式分母系數(shù)向量 <2> 解答： 思路是 將題目給的模擬頻率直接轉(zhuǎn)化成數(shù)字頻率。%ΩT=w Ω=2πfp1 等等 Fs=25000;wp1=2*5000/Fs;wp2=2*7000/Fs;ws1=2*3500/Fs;ws2=2*8500/Fs; wp=[wp1,wp2];ws=[ws1,ws2];Rp=0.5;As=

25、45; 然后用[N,wc]=ellipord(wp,ws,Rp,As); %用ellipord求階次N和wc 求出階次N和wc 再用[b,a]=ellip(N,Rp,As,wc); %求系統(tǒng)函數(shù)分子和分母向量 求出系統(tǒng)函數(shù)的分子和分母向量 最后用freqz來(lái)求頻率響應(yīng)，并畫出圖。 程序1： %根據(jù)題目給的模擬頻率直接轉(zhuǎn)化成數(shù)字頻率 clear all;clc; %ΩT=w Ω=2πfp1 等等 Fs=25000;wp1=2*5000/Fs;wp2=2*7000/Fs;ws1=2*3500/Fs;ws2=2*8500/Fs; wp=[

26、wp1,wp2];ws=[ws1,ws2];Rp=0.5;As=45; [N,wc]=ellipord(wp,ws,Rp,As); %用ellipord求階次N和wc [b,a]=ellip(N,Rp,As,wc); %求系統(tǒng)函數(shù)分子和分母向量 hn=dimpulse(b,a); %求離散時(shí)間線性系統(tǒng)的脈沖響應(yīng) w0=[ws1*pi,wp1*pi,wp2*pi,ws2*pi]; Hx=freqz(b,a,w0); %求在w0時(shí)的H（exp（iw））的值 [H,w]=freqz(b,a); dbHx=-2

27、0*log10(abs(H)/max(abs(H))); [db,mag,pha,grd,w]=freqz_m(b,a); subplot(211);plot(12.5*w/pi,db,linewidth,2);title(橢圓函數(shù)數(shù)字帶通幅度響應(yīng)（db）); xlabel(f(KHz));ylabel(dB);axis([0,12.5,-50,3]);%axis([xmin xmax ymin ymax])設(shè)置坐標(biāo)的范圍 set(gca,xtickmode,manual,xtick,[1.25,3.5,5,7,8.5,10,12.5]); set(gca,ytickmode,man

28、ual,ytick,[-50,-As,-10,0]);grid; subplot(2,1,2);plot(w/pi,pha,linewidth,2);title(橢圓函數(shù)數(shù)字帶通相位響應(yīng)); xlabel(\omega/\pi);%w/π ylabel(rad);grid;axis([0,1,-4,4]); figure(2);stem(hn,.,linewidth,2);title(系統(tǒng)沖激響應(yīng));xlabel(n);axis([0,121,-0.15,0.15]); N=4; 用FDAtool工具設(shè)計(jì) 將參數(shù)填寫進(jìn)去即可。 2、 同理可得 程序2

29、%根據(jù)題目給的模擬頻率直接轉(zhuǎn)化成數(shù)字頻率 clear all;clc; %ΩT=w Ω=2πfp1 等等 Fs=25000;wp1=2*5000/Fs;wp2=2*7000/Fs;ws1=2*3500/Fs;ws2=2*8500/Fs; wp=[wp1,wp2];ws=[ws1,ws2];Rp=0.5;As=45; [N,wc]=cheb1ord(wp,ws,Rp,As); [b,a]=cheby1(N,Rp,wc); hn=dimpulse(b,a); %求h（n） w0=[ws1*pi,wp1*pi,wp2*pi,ws2*pi]; Hx=freqz(b,a,w0)

30、; %求在w0時(shí)的H（exp（iw））的值 [H,w]=freqz(b,a); dbHx=-20*log10(abs(H)/max(abs(H))); [db,mag,pha,grd,w]=freqz_m(b,a); subplot(211);plot(12.5*w/pi,db,linewidth,2);title(切比雪夫I型數(shù)字帶通幅度響應(yīng)（db）); xlabel(f(KHz));ylabel(dB);axis([0,12.5,-50,3]);%axis([xmin xmax ymin ymax])設(shè)置坐標(biāo)的范圍 set(gca,xtickmode,

31、manual,xtick,[1.25,3.5,5,7,8.5,10,12.5]); set(gca,ytickmode,manual,ytick,[-50,-As,-10,0]);grid; subplot(2,1,2);plot(w/pi,pha,linewidth,2);title(切比雪夫I型數(shù)字帶通相位響應(yīng)); xlabel(\omega/\pi);%w/π ylabel(rad);grid;axis([0,1,-4,4]); figure(2);stem(hn,.,linewidth,2);title(系統(tǒng)沖激響應(yīng));xlabel(n);axis([0,121,-0.15,

32、0.15]); N=5; 3、 程序3 %根據(jù)題目給的模擬頻率直接轉(zhuǎn)化成數(shù)字頻率 clear all;clc; %ΩT=w Ω=2πfp1 等等 Fs=25000;wp1=2*5000/Fs;wp2=2*7000/Fs;ws1=2*3500/Fs;ws2=2*8500/Fs; wp=[wp1,wp2];ws=[ws1,ws2];Rp=0.5;As=45; [N,wc]=cheb2ord(wp,ws,Rp,As); [b,a]=cheby2(N,As,wc); hn=dimpulse(b,a); %求h（n） w0=[ws1*pi,wp1*pi,wp2*

33、pi,ws2*pi]; Hx=freqz(b,a,w0); %求在w0時(shí)的H（exp（iw））的值 [H,w]=freqz(b,a); dbHx=-20*log10(abs(H)/max(abs(H))); [db,mag,pha,grd,w]=freqz_m(b,a); subplot(211);plot(12.5*w/pi,db,linewidth,2);title(切比雪夫II型數(shù)字帶通幅度響應(yīng)（db）); xlabel(f(KHz));ylabel(dB);axis([0,12.5,-50,3]);%axis([xmin xmax ymin yma

34、x])設(shè)置坐標(biāo)的范圍 set(gca,xtickmode,manual,xtick,[1.25,3.5,5,7,8.5,10,12.5]); set(gca,ytickmode,manual,ytick,[-50,-As,-10,0]);grid; subplot(2,1,2);plot(w/pi,pha,linewidth,2);title(切比雪夫II型數(shù)字帶通相位響應(yīng)); xlabel(\omega/\pi);%w/π ylabel(rad);grid;axis([0,1,-4,4]); figure(2);stem(hn,.,linewidth,2);title(系統(tǒng)沖激響

35、應(yīng));xlabel(n);axis([0,121,-0.15,0.15]); N=5; 4、 程序4 %根據(jù)題目給的模擬頻率直接轉(zhuǎn)化成數(shù)字頻率 clear all;clc; %ΩT=w Ω=2πfp1 等等 Fs=25000;wp1=2*5000/Fs;wp2=2*7000/Fs;ws1=2*3500/Fs;ws2=2*8500/Fs; wp=[wp1,wp2];ws=[ws1,ws2];Rp=0.5;As=45; [N,wc]=buttord(wp,ws,Rp,As); [b,a]=butter(N,wc); hn=dimpulse(b,a); %求

36、h（n） w0=[ws1*pi,wp1*pi,wp2*pi,ws2*pi]; Hx=freqz(b,a,w0); %求在w0時(shí)的H（exp（iw））的值 [H,w]=freqz(b,a); dbHx=-20*log10(abs(H)/max(abs(H))); [db,mag,pha,grd,w]=freqz_m(b,a); subplot(211);plot(12.5*w/pi,db,linewidth,2);title(巴特沃斯型數(shù)字帶通幅度響應(yīng)（db）); xlabel(f(KHz));ylabel(dB);axis([0,12.5,-50,3])

37、;%axis([xmin xmax ymin ymax])設(shè)置坐標(biāo)的范圍 set(gca,xtickmode,manual,xtick,[1.25,3.5,5,7,8.5,10,12.5]); set(gca,ytickmode,manual,ytick,[-50,-As,-10,0]);grid; subplot(2,1,2);plot(w/pi,pha,linewidth,2);title(巴特沃斯型數(shù)字帶通相位響應(yīng)); xlabel(\omega/\pi);%w/π ylabel(rad);grid;axis([0,1,-4,4]); figure(2);stem(hn,.,

38、linewidth,2);title(系統(tǒng)沖激響應(yīng));xlabel(n);axis([0,121,-0.15,0.15]); N=7; <1> 解答： 1、基本的編程思路是 ①確定濾波器的設(shè)計(jì)指標(biāo) ②設(shè)計(jì)文博參數(shù)δ1和δ2 ③調(diào)用firpmord求解M ④調(diào)用firpm求出系統(tǒng)沖激響應(yīng) ⑤調(diào)用freqz求解濾波器頻率響應(yīng) 程序： clear all; clc; Rp=1;As=75; Fs=100000;f=[25000,30000];a=[1,0]; dev=[(10^(Rp/2

39、0)-1)/(10^(Rp/20)+1),10^(-As/20)]; [M, f0, a0, weights] = firpmord(f,a,dev,Fs); %用firpmord函數(shù)求M h=firpm(M,f0,a0,weights); %用firpm求解系統(tǒng)沖激響應(yīng) [H,f]=freqz(h,1,1024,Fs); %用freqz求解頻率響應(yīng) subplot(211);plot(f/Fs,20*log10(abs(H)),linewidth,2); title(幅度響應(yīng)（dB）); xlabel(f/Fs);ylabel(20log|H(e^j^\om

40、ega)|(dB)); axis([0,0.5,-140,20]); set(gca,xtickmode,manual,xtick,[0,0.1,0.2,0.25,0.3,0.4,0.5]); set(gca,ytickmode,manual,ytick,[-140,-110,-70,-30,0,20]);grid; subplot(212); plot(f/Fs,angle(H),linewidth,2);grid; title(相位響應(yīng));xlabel(f/Fs);ylabel(arg[H(e^j^\omega)]); 8.18 解答： clear

41、all; clc; Rp=1;As=50; f=[0.1,0.2,0.4,0.5,0.6,0.7,0.8,0.9];a=[0,1,0,1,0]; delta1=(10^(Rp/20)-1)/(10^(Rp/20)+1); delta2=(1+delta1)*(10^(-As/20)); dev=[delta2,delta1,delta2,delta1,delta2]; [M, f0, a0, weights] = firpmord(f,a,dev); %用firpmord函數(shù)求M h=firpm(M+4,f0,a0,weights); %用firpm求解系統(tǒng)沖激響應(yīng) 經(jīng)

42、驗(yàn)算取M+4才能滿足設(shè)計(jì)要求 [H,f]=freqz(h,1,1024); %用freqz求解頻率響應(yīng) subplot(211);plot(f/pi,20*log10(abs(H)),linewidth,2);grid; title(幅度響應(yīng)（dB）);xlabel(\omega/\pi);ylabel(20log|H(e^j^\omega)|(dB)); axis([0,1,-100,10]); set(gca,xtickmode,manual,xtick,[0:0.1:1]);grid on; set(gca,ytickmode,manual,ytick,[-100,

43、-50,-1,10]);grid on; subplot(212); plot(f/pi,angle(H),linewidth,2);grid; title(相位響應(yīng));xlabel(\omega/\pi);ylabel(arg[H(e^j^\omega)]); figure(2);stem(h,.,linewidth,2);title(系統(tǒng)沖激響應(yīng));xlabel(n);axis([0,40,-0.4,0.6]); <2>頻域抽樣定理的驗(yàn)證 給定信號(hào) 。 1．利用DTFT計(jì)算信號(hào)的頻譜，一個(gè)周期內(nèi)角頻率離散為M=1000點(diǎn)，畫出頻譜圖，標(biāo)明坐標(biāo)軸。 解答：先畫圖表

44、示出給定的序列，利用for循環(huán) n=0:100; %先將給定的信號(hào)表示出來(lái)，即 for n1=0:13 x(n1+1)=n1+1; %不需要從x（0）開(kāi)始 只要畫圖時(shí)n從0開(kāi)始畫就可以 end for n2=14:26; x(n2+1)=27-n2; end for n3=27:100 x(n3+1)=0; end figure(1); stem(n,x); xlabel(n); ylabel(x(n)); title(給定的信號(hào)序列); 然后利用freqz做DTFT運(yùn)算， M=1000; w=0:2*pi/

45、M:(2-1/M)*pi; %為了下面的方便 要畫出0~2pi X=freqz(x,1,M,whole); %[H,f]=freqz(h,1,1024); %用freqz求解頻率響應(yīng) %其中h是單位沖擊響應(yīng) X1=abs(X);figure(2);plot(w,X1);xlabel(w);ylabel(X); title(x(n)的DTFT幅頻響應(yīng)); %2．分別對(duì)信號(hào)的頻譜在(0~2pi)區(qū)間上等間隔抽樣16點(diǎn)和32點(diǎn)， %得到x32(k)和x16(k)。離散傅里葉反變換

46、后得到時(shí)域信號(hào)和。 2．分別對(duì)信號(hào)的頻譜在區(qū)間上等間隔抽樣16點(diǎn)和32點(diǎn)，得到和。離散傅里葉反變換后得到時(shí)域信號(hào)和。 解答： N1=16; w1=0:2*pi/16:(1-1/16)*2*pi; Xk16=freqz(x,1,N1,whole); 實(shí)現(xiàn)等間隔抽樣 % Xk16=fft(x,N1); figure(3); stem(w1,abs(Xk16)); xlabel(w);ylabel(X16（k）);title(fft后求出的X16（k）); N2=32; w2=0:2*pi/32:(1-1/32)*2*pi;

47、 Xk32=freqz(x,1,N2,whole); %對(duì)X做16點(diǎn)抽樣 figure(4); stem(w2,abs(Xk32)); xlabel(w);ylabel(X32（k）);title(fft求出的X32（k）); 3．畫出信號(hào)和的圖形，計(jì)算與和的均方誤差。從時(shí)域角度上進(jìn)行對(duì)比和分析，驗(yàn)證頻域抽樣定理。 解答：利用ifft求fft的反變換，其中注意求出來(lái)的時(shí)域信號(hào)，長(zhǎng)度變化了。 %ifft求反 N3=16; w3=0:1:15; xk16=ifft(Xk16,N3); %對(duì)X做16點(diǎn) figure(5); stem(w3,x

48、k16); xlabel(w);ylabel(x16（n）);title(ifft后求出的x16（n）); N4=32; w4=0:1:31; xk32=ifft(Xk32,N4); %對(duì)X做16點(diǎn)抽樣 figure(6); stem(w4,xk32); xlabel(w);ylabel(x32（n）);title(ifft后求出的x32（n）); 從圖形上來(lái)看，16點(diǎn)的抽樣做ifft之后還原回來(lái)的信號(hào)時(shí)域上已經(jīng)發(fā)生了嚴(yán)重的混疊失真；而32點(diǎn)的抽樣還原回來(lái)的時(shí)域信號(hào)和原來(lái)的信號(hào)具有較好的一致性，從而證明了頻域抽樣定理，即對(duì)于M點(diǎn)的有限長(zhǎng)序列，頻域抽樣不失

49、真的條件是頻域抽樣點(diǎn)數(shù)N要大于或者等于M。 計(jì)算均方誤差，直接根據(jù)公式寫出來(lái)就可以。 %計(jì)算均方誤差 for c=17:27 xk16(c)=12; %根據(jù)周期性將剩余的點(diǎn)補(bǔ)齊 值均為12 end sum=0; for n=1:27 sum=sum+(abs(x(n))-abs(xk16(n)))^2; end d=sum/27; %計(jì)算32點(diǎn)的均方誤差 sum1=0; for n=1:27 sum1=sum1+(abs(x(n))-abs(xk32(n)))^2; end d1=sum1/27; 結(jié)

50、果： >> d d = 37.4815 >> d1 d1 = 3.6521e-31 計(jì)算得到的16點(diǎn)抽樣還原信號(hào)與原來(lái)信號(hào)的方差是37.4815 32點(diǎn)抽樣還原信號(hào)與原來(lái)信號(hào)的方差是3.6521e-31 可以近似為0。 4．利用內(nèi)插公式，由和分別得到的估計(jì)值，計(jì)算均方誤差，從頻域角度驗(yàn)證頻率抽樣定理。 解答： 利用內(nèi)插公式，即 在程序中直接寫出來(lái)就可以。 %內(nèi)插公式 M=16; L=1000; for b=0:L-1 wn=2*pi*b/L; sum=0; for k1=0:M-1 if si

51、n((wn-(2*pi*k1/M))/2)~=0 sum=sum+(1/M)*Xk16(k1+1)*(sin((wn-(2*pi*k1/M))*M/2)/sin((wn-(2*pi*k1/M))/2))*exp(-j*(M-1)/2*(wn-(2*pi*k1/M))); else sum=sum+Xk16(k1+1); end end Xejw16(b+1)=sum; end k2=0:L-1; w2=2*pi*k2/L; figure(7); plot(w2,abs(Xejw16)); xlabel(w/

52、rad); ylabel(|Xejw|); title(16點(diǎn)DTFT抽樣還原Xejw); %內(nèi)插 32點(diǎn) M=32; L=1000; for b=0:L-1 wn=2*pi*b/L; sum=0; for k1=0:M-1 if sin((wn-(2*pi*k1/M))/2)~=0 sum=sum+(1/M)*Xk32(k1+1)*(sin((wn-(2*pi*k1/M))*M/2)/sin((wn-(2*pi*k1/M))/2))*exp(-j*(M-1)/2*(wn-(2*pi*k1/M))); else

53、 sum=sum+Xk32(k1+1); end end Xejw32(b+1)=sum; end k2=0:L-1; w2=2*pi*k2/L; figure(8); plot(w2,abs(Xejw32)); xlabel(w/rad); ylabel(|Xejw|); title(32點(diǎn)DTFT抽樣還原Xejw); 從圖像上看，與原來(lái)信號(hào)做對(duì)比 原信號(hào)： 顯然16點(diǎn)的抽樣還原后有明顯的混疊失真，而32點(diǎn)的則吻合的很好。從而從頻域證明了頻域抽樣定理，即對(duì)于M點(diǎn)的有限長(zhǎng)序列，頻域抽樣不失真的條件是頻

54、域抽樣點(diǎn)數(shù)N要大于或者等于M。 求均方誤差 %均方誤差 sum=0; for c=1:L sum=sum+(abs(X(c))-abs(Xejw16(c)))^2; end d3=sum/L; %d3 第三個(gè)均方誤差 sum=0; for c=1:L sum=sum+(abs(X(c))-abs(Xejw32(c)))^2; end d4=sum/L; 結(jié)果： >> d3 d3 = 230.2132 >> d4 d4 = 1.0918e-26 可以看出32點(diǎn)抽

55、樣的方差幾乎為零。 <3>音頻信號(hào)處理 1、語(yǔ)音信號(hào)的采集 　 利用Windows 附件中的錄音機(jī)，錄制一段自己的話音，時(shí)間在1 s內(nèi)。在Matlab軟件平臺(tái)下，利用函數(shù)wavread對(duì)語(yǔ)音信號(hào)進(jìn)行采樣，記住采樣頻率和采樣點(diǎn)數(shù)。通過(guò)wavread函數(shù)的使用，理解采樣頻率、采樣點(diǎn)數(shù)等概念。 解答： clc;clear;close all [z1,fs,bits]=wavread(d:\123.wav);%調(diào)用函數(shù)畫出時(shí)域信號(hào) 利用wavread函數(shù)得到的音頻文件存在z1中，同時(shí)還得到了視頻文件的抽樣頻率fs和抽樣點(diǎn)數(shù)bits。 >> fs fs = 44100

56、 >> bits bits =16 2、語(yǔ)音信號(hào)的頻譜分析 　 畫出語(yǔ)音信號(hào)的時(shí)域波形，然后對(duì)語(yǔ)音號(hào)進(jìn)行快速傅里葉變換，得到原始模擬信號(hào)的頻譜特性，畫出頻譜圖，標(biāo)注坐標(biāo)軸。 解答： 利用fft快速算法直接對(duì)得到的時(shí)域信號(hào)做運(yùn)算。 figure(1) plot(z1); title(原始語(yǔ)音信號(hào)); y1=z1(1:length(z1)); Y1=fft(y1);%直接調(diào)用得到原始頻域信號(hào) n=1:length(Y1); figure(2) plot(n,abs(Y1)); title(原始信號(hào)的頻譜特性);n1=0:31; 注意要對(duì)z1

57、抽取。 3、用濾波器對(duì)信號(hào)進(jìn)行濾波 　在離散時(shí)間域，使信號(hào)通過(guò)沖激響應(yīng)為 的低通濾波器，得到系統(tǒng)的輸出。 解答： h=0.5.*(1-cos(2*pi.*n1./31));%通過(guò)系統(tǒng)及時(shí)域卷積 z2=conv(y1,h); 利用conv做時(shí)域的卷積運(yùn)算，得到低通濾波后的時(shí)域信號(hào)。 4、比較濾波前后語(yǔ)音信號(hào)的波形及頻譜 畫出濾波前后的時(shí)域波形及頻譜，進(jìn)行比較。 解答： n2=0:length(z1)+30; figure(3) plot(n2,z2); title(濾波后時(shí)域波形); yyy=fft(z2); figure(4) plot(n2,abs(yyy));title(濾波后頻域波形); 注意卷紙之后時(shí)域信號(hào)的長(zhǎng)度變化，L=N1+N2-1；所以n2要加上30； 5、回放語(yǔ)音信號(hào) 在Matlab中，函數(shù)sound可以對(duì)聲音進(jìn)行回放，調(diào)用格式：sound(x，fs，bits)，可以感覺(jué)濾波前后的聲音有變化，音量較原來(lái)變大了，并且伴有雜音。 sound(z1,44100,16);%原來(lái)語(yǔ)音 sound(z2,44100,16);%濾波后語(yǔ)音