信號(hào)與信號(hào)處理實(shí)驗(yàn)參考答案_新0_20070611.doc

實(shí)驗(yàn)一 熟悉MATLAB環(huán)境2、（2）粗略描繪下列各函數(shù)的波形說(shuō)明：MATLAB中有函數(shù) f(t)=sint/tt=-3*pi:0.01*pi:3*pi; t1=t/pi;y=sinc(t1); plot(t,y);%輸出一個(gè)波形，以t為縱坐標(biāo)，y函數(shù)為橫坐標(biāo)。 hold on%不用的話 前面的圖會(huì)消失的，要在第一個(gè)波形后面用。plot(t,0)%輸出一個(gè)參考坐標(biāo)軸。在自己的工作目錄work下創(chuàng)建Heaviside函數(shù)的M文件，該文件如下：function f=Heaviside(t)f=(t0) %t0時(shí)f為1，否則f為0在命令窗口輸入如下語(yǔ)句，就能繪出u(t)的波形。t=-1:0.01:3;f=Heaviside(t);plot(t,f)%繪圖的加不加分號(hào)都行。axis(-1 3 0.2 1.2)%只顯示前面兩個(gè)是橫坐標(biāo)范圍，后面兩個(gè)縱坐標(biāo)范圍內(nèi)的圖形。t=-1:0.01:2;g=Heaviside(t)-Heaviside(t-1); plot(t,g);axis(-1 2 -0.2 1.2)hold onplot(t,0)4、 分別用for和while循環(huán)結(jié)構(gòu)編寫程序，求出s=1+2+22+23+262+263并考慮一種避免循環(huán)的簡(jiǎn)潔方法來(lái)進(jìn)行求和。 程序如下： s=1;for k=1:63%么有逗號(hào)。 s=s+2k;ends運(yùn)行結(jié)果是：s =1. 8447e+019（2）s=1;k=1;while k=63%沒(méi)有分號(hào)。 s=s+2k; k=k+1;ends運(yùn)行結(jié)果：s =1. 8447e+019(3) k=0:63;s=sum(2.k)實(shí)驗(yàn)二 信號(hào)的卷積與系統(tǒng)的響應(yīng)和階躍響應(yīng)1 n=0:20;hn=0.9.n;xn=0,0 ones(1,8),0,0;%離散的所以是用矩陣賦值，連續(xù)的用HEAVISIDE。yn=conv(hn,xn);%求卷積，系統(tǒng)在前，輸入在后。stem(yn)%離散輸出波形。3. 利用MATLAB繪制下列信號(hào)的卷積積分f1(t)*f2(t)的時(shí)域波形。（1） f1(t)=2u(t+1)-u(t-1), f2(t)=u(t+2)-u(t-2)（2） f1(t)=tu(t), f2(t)=u(t)（3） f1(t)=u(t)-u(t-4), f2(t)=sin(t)u(t);（4） f1(t)=e-2tu(t), f2(t)=e-tu(t)(1) 先編寫實(shí)現(xiàn)連續(xù)信號(hào)卷積的通用函數(shù)sconv(),程序如下：functionf,k=sconv(f1,f2,k1,k2,p)%計(jì)算連續(xù)信號(hào)卷積積分f(t)=f1(t)*f2(t)%f:卷積積分f(t)對(duì)應(yīng)的非零樣值向量%k:f(t)的對(duì)應(yīng)時(shí)間向量%f1:f1(t)非零樣值向量%f2:f2(t)非零樣值向量%k1:f1(t)的對(duì)應(yīng)時(shí)間向量%k2:序列f2(t)的對(duì)應(yīng)時(shí)間向量%p:取樣時(shí)間間隔f=conv(f1,f2); %計(jì)算序列f1與f2的卷積ff=f*p;k0=k1(1)+k2(1); %計(jì)算序列f的非零樣值的起點(diǎn)位置k3=length(f1)+length(f2)-2; %計(jì)算卷積和f的非零樣值的寬度k=k0:p:(k3-(0-k0)/p)*p); %確定卷積和f非零樣值的時(shí)間向量subplot(2,2,1)plot(k1,f1)title(f1(t)xlabel(t)ylabel(f1(t)subplot(2,2,2)plot(k2,f2)title(f2(t)xlabel(t)ylabel(f2(t)subplot(2,2,3)plot(k,f);h=get(gca,position);h(3)=2.5*h(3);set(gca,position,h)title(f(t)=f1(t)*f2(t)xlabel(t)ylabel(f(t)p=0.01;k1=-1:p:1;f1=2*ones(1,length(k1);k2=-2:p:2;f2=ones(1,length(k2);f,k=sconv(f1,f2,k1,k2,p)（2）p=0.01;k1=0:p:10;f1=k1k2=0:p:10;f2=ones(1,length(k2);f,k=sconv(f1,f2,k1,k2,p)第（2）題圖上，第（3）題圖下，5已知描述某連續(xù)系統(tǒng)的微分方程為：y(t)+5y(t)+8y(t)=3f(t)+2f(t)繪出系統(tǒng)的沖激響應(yīng)波形，求出t=0.5s, 1s, 1.5s, 2s系統(tǒng)沖激響應(yīng)的數(shù)值解。6、step()函數(shù)的調(diào)用格式與impulse()函數(shù)類似。已知描述某連續(xù)系統(tǒng)的微分方程為：2y(t)+y(t)+8y(t)=f(t)試用MATLAB繪出該系統(tǒng)的沖激響應(yīng)與階躍響應(yīng)的波形。5. 程序如下：a=1 5 8;b=3 0 2;impulse(b,a)y=impulse(b,a,0:0.5:2)y = -15.0000 -1.3292 0.6302 0.3941 0.12656. a=2 1 8;b=1; impulse(b,a) hold on step(b,a)legend(沖激響應(yīng),階躍響應(yīng))7. b=1;a=1 1 0.9; n=-10:50;impz(b,a,n)實(shí)驗(yàn)三 利用DFT分析模擬信號(hào)頻譜1T=0.1;n=0:20;xn=exp(-2*n.*T);k=0:31Xk=fft(xn,32);XM=T.*XkMagXk=abs(XM)plot(k,XM)xlabel(k);ylabel(XM)矩形脈沖的頻譜T=0.1;n=-10:10;xn=Heaviside(n*T+1)-Heaviside(n*T-1)Xk=fft(xn,256);Xk=fftshift(T*abs(Xk);plot(Xk)2fs=400;n=0:160;xn=cos(2*pi*100*n/fs)+0.75*cos(2*pi*110*n/fs);Xk=fft(xn,256);Xk=(abs(Xk)/fs);k=0:255;plot(fs/256*k,Xk)另一個(gè)：fs=400;n=0:160;xn=cos(2*pi*100*n/fs)+0.75*cos(2*pi*110*n/fs);Xk=fft(xn,256);Xk=(abs(Xk)/fs);Xk=fftshift(Xk);k=0:255;plot(-fs/2+fs/256*k,Xk)3t=-0.2:0.005:0.2;x=square(2*pi*20*t,50); %產(chǎn)生頻率為20Hz的方波subplot(1,2,1);plot(t,x);axis(-0.2 0.2 -1.5 1.5)T=0.0025;N=20;n=0:N-1;xn=square(2*pi*20.*n*T,50)yn=0.5*(xn+1)Yn=fft(yn);Yn=fftshift(abs(Yn);Ym=1/N*Ynsubplot(1,2,2);stem(1/T)/N*n,Ym)實(shí)驗(yàn)五 離散信號(hào)與離散系統(tǒng)分析基礎(chǔ)參考答案functionz,p=ljdt(A,B)%繪離散系統(tǒng)零極點(diǎn)圖的函數(shù)p=roots(A); %求系統(tǒng)極點(diǎn)z=roots(B); %求系統(tǒng)零點(diǎn)p=p; %將極點(diǎn)列向量轉(zhuǎn)置為行向量z=z; %將零點(diǎn)列向量轉(zhuǎn)置為行向量x=max(abs(p z 1); %確定縱坐標(biāo)范圍x=x+0.1;y=x; %確定橫坐標(biāo)范圍clchold onaxis(-x x -y y) %確定坐標(biāo)軸顯示范圍w=0:pi/300:2*pi;t=exp(i*w); %畫單位圓plot(t)axis(square)plot(-x x,0 0) %畫橫坐標(biāo)軸plot(0 0,-y y) %畫縱坐標(biāo)軸text(0.1,x,jImz)text(y,0.1,Rez)plot(real(p),imag(p),x) %畫極點(diǎn)plot(real(z),imag(z),o) %畫零點(diǎn) title(離散系統(tǒng)零極點(diǎn)圖) %標(biāo)注標(biāo)題hold off2. 離散時(shí)間系統(tǒng)的系統(tǒng)函數(shù)如下所示，繪制系統(tǒng)的零極點(diǎn)圖，并判斷系統(tǒng)的穩(wěn)定性。(1) (2)A=1 3 2 1;B=3 2 2 5; ljdt(A,B)（2）A= 1 2 -4 1; B=1 0 2;ljdt(A,B)3B=1 2 1; A=1 -0.5 -0.005 0.3; ljdt(A,B) H,w=freqz(B,A); figure;subplot(1,2,1);plot(w/pi,abs(H) subplot(1,2,2);plot(w/pi,angle(H)4.A=1 0 0.61;B=1 -2 1;H,w=freqz(B,A);subplot(1,2,1)plot(w/pi,abs(H)xlabel(頻率/pi)ylabel(頻率響應(yīng)幅度)grid onsubplot(1,2,2)plot(w/pi,angle(H)xlabel(頻率/pi)ylabel(頻率響應(yīng)相位)grid on5.B=1 1 1;A=1 1 -2 0;r p k=residue(B,A)r = 0.5000 1.0000 -0.5000p = -2 1 0k =實(shí)驗(yàn)六 離散傅里葉變換參考程序1 dft1.m: 用for循環(huán)計(jì)算DFTfunction Am,pha=dft1(x)N=length(x);w=exp(-j*2*pi/N);for k=1:N sum=0; for n=1:N sum=sum+x(n)*w(k-1)*(n-1); end Am(k)=abs(sum); pha(k)=angle(sum);end調(diào)用dft1.m函數(shù)文件的M程序文件：x=ones(1,8),zeros(1,248);t=cputime;Am1,pha1=dft1(x);t1=cputime-t;n=0:(length(x)-1);w=(2*pi/length(x)*n;figure(1)subplot(2,1,1),plot(w,Am1,b);title(Magnitude fart);xlabel(frequency in radians);ylabel(|X(exp(jw)|)subplot(2,1,2),plot(w,pha1,r);grid;gtext(Phase part);xlabel(frequency in radians);ylabel(arg(Xexp(jw)/radians);2dft2.m: 用MATLAB矩陣運(yùn)算計(jì)算DFTfunction Am,pha=dft2(x)N=length(x);n=0:N-1;k=0:N-1;w=exp(-j*2*pi/N);nk=n*k;wnk=w.nk;Xk=x*wnk;Am=abs(Xk);pha=angle(Xk);3. dft3.m：調(diào)用FFT庫(kù)函數(shù)計(jì)算DFTfunction Am,pha=dft3(x) Xk=fft(x); Am=abs(Xk); pha=angle(Xk)t1 =0. 32802采樣16點(diǎn)，做16點(diǎn)DFT,程序如下：fs=32000;n=1:16;x=cos(2*pi*6.5*1000*n/fs)+cos(2*pi*7*1000*n/fs)+cos(2*pi*9*1000*n/fs);subplot(2,2,1);stem(n,x);title(signal x(n)Xk=fft(x,16);magXk=abs(Xk);k1=0:15;w1=2*pi/16*k1;subplot(2,1,2);plot(w1,magXk);title(16點(diǎn)DTFT的幅值);xlabel(w(單位：) 采樣16點(diǎn)，做256點(diǎn)DFT,程序如下：fs=32000;n=1:16;x=cos(2*pi*6.5*1000*n/fs)+cos(2*pi*7*1000*n/fs)+cos(2*pi*9*1000*n/fs);subplot(2,2,1);stem(n,x);title(signal x(n)Xk=fft(x,256);magXk=abs(Xk);k1=0:255;w1=2*pi/256*k1;subplot(2,1,2);plot(w1,magXk);title(16點(diǎn)DTFT的幅值);xlabel(w(單位：pi)結(jié)果如圖二（3）采集數(shù)據(jù)長(zhǎng)度N=256點(diǎn)，做256點(diǎn)DFT,結(jié)果如圖三fs=32000;n=0:255;x=cos(2*pi*6.5*1000*n/fs)+cos(2*pi*7*1000*n/fs)+cos(2*pi*9*1000*n/fs);subplot(2,2,1);stem(n,x);title(signal x(n)Xk=fft(x,256);magXk=abs(Xk);k1=0:255;w1=2*pi/256*k1;subplot(2,1,2);plot(w1,magXk);title(256點(diǎn)DTFT的幅值);xlabel(w(單位：pi)3參考程序和運(yùn)行結(jié)果n=0:1:15;xn=0:1:15;subplot(2,2,1)stem(n,xn);xlabel(n);ylabel(x(n);m=0:1:7;hn=ones(1,8);subplot(2,2,2)stem(m,hn)xlabel(n);ylabel(h(n);N=length(n)+length(m)-1;Xk=fft(xn,N);Hk=fft(hn,N);Yk=Xk.*Hk;yn=ifft(Yk,N);if all(imag(xn)=0)&(all(imag(hn)=0) %實(shí)序列的循環(huán)卷積仍為實(shí)序列 yn=real(yn)endsubplot(2,2,3)x=0:N-1;stem(x,yn)xlabel(n);ylabel(y(n)=x(n)*h(n)實(shí)驗(yàn)八 IIR數(shù)字濾波器的設(shè)計(jì)1wp=0.2*pi; %digital Passband freq in radws=0.3*pi; %digital Stopband freq in radFs=1000; wp=0.2*pi*Fs;ws=0.3*pi*Fs; Rp=1;Rs=15; N,Wn=buttord(wp,ws,Rp,Rs,s); 返回模擬濾波器的最小階數(shù)和截止頻率Z,P,K=buttap(N); 模擬低通濾波器原型Bap,Aap=zp2tf(Z,P,K);b,a=lp2lp(Bap,Aap,Wn);bz,az=bilinear(b,a,Fs);H,W=freqz(bz,az);plot(W,abs(H);gridxlabel(頻率/弧度) ylabel(頻率響應(yīng)幅度)H,W=freqz(bz,az,50)figurestem(W,abs(H)23. 參考程序fs=1000;wp=400*2/fs;ws=200*2/fs;rp=3;rs=15;Nn=128;N,wn=buttord(wp,ws,rp,rs)b,a=butter(N,wn,high)freqz(b,a,Nn,fs)運(yùn)行結(jié)果N = 2wn = 0.6630b = 0.1578 -0.3155 0.1578a = 1.0000 0.6062 0.23734參考程序fs=1000;wp=100 250*2/fs;ws=50 300*2/fs;rp=3;rs=30;Nn=128;N,wn=cheb2ord(wp,ws,rp,rs)b,a=cheby2(N,rp,wn,bandpass)freqz(b,a,Nn,fs)運(yùn)行結(jié)果實(shí)驗(yàn)九 FIR數(shù)字濾波器的設(shè)計(jì)1參考程序clear all;Wp=0.2*pi;Ws=0.4*pi;tr_width=Ws-Wp; %過(guò)渡帶寬度N=ceil(6.6*pi/tr_width)+1 %濾波器長(zhǎng)度（選擇海明窗）n=0:1:N-1;Wc=(Ws+Wp)/2; %理想低通濾波器的截止頻率alpha=(N-1)/2;m=n-alpha+eps;hd=sin(Wc*m)./(pi*m); %理想低通濾波器的單位沖激響應(yīng)w_ham=(hamming(N); %海明窗h=hd.*w_ham; %截取得到實(shí)際的單位脈沖響應(yīng)H,w=freqz(h,1,1000,whole);H=(H(1:1:501);w=(w(1:1:501);mag=abs(H);db=20*log10(mag+eps)/max(mag); %相對(duì)幅值響應(yīng)pha=angle(H);delta_w=2*pi/1000;Ap=-(min(db(1:1:Wp/delta_w+1) %實(shí)際通帶文波As=-round(max(db(Ws/delta_w+1:1:501) %實(shí)際阻帶文波subplot(221)stem(n,hd)title(理想單位脈沖響應(yīng)hd(n)subplot(222)stem(n,w_ham)title(海明窗w(n)subplot(223)stem(n,h)title(實(shí)際單位脈沖響應(yīng)h(n)subplot(224)plot(w/pi,db)title(幅度響應(yīng)(dB)axis(0,1,-100,10)運(yùn)行結(jié)果N = 34Ap = 0.0477As = 5291另一程序clear all;Wp=0.2*pi;Ws=0.4*pi;tr_width=Ws-Wp; %過(guò)渡帶寬度N=ceil(6.6*pi/tr_width)+1 %濾波器長(zhǎng)度（選擇海明窗）n=0:1:N-1;Wc=(Ws+Wp)/2; %理想低通濾波器的截止頻率h=fir1(N,Wc/pi,hamming(N+1)subplot(221)stem(h)H,w=freqz(h,1,512);mag=abs(H);db=20*log10(mag+eps)/max(mag); %相對(duì)幅值響應(yīng)subplot(222)plot(w/pi,db)2window=hanning(16); b=fir1(15,0.3 0.5,window); freqz(b,1,512)window=hanning(46); b=fir1(45,0.3 0.5,window); freqz(b,1,512)4. N=16;f=0:8/8;a=1,1,1,1,0.389,0,0,0,0;h=fir2(N,f,a);freqz(h,1)實(shí)驗(yàn)十 信號(hào)濾波：參考程序1 IIR濾波器實(shí)現(xiàn)clcclearFs=1000;T=1/10;t=0:1/Fs:1;x=square(2*pi*10*t);%產(chǎn)生幅度為-1,1的方波plot(t,x)title(矩形波)axis(0,0.2,-1.5,1.5)N=T*Fs; %方波的頻譜t1=0:1/Fs:(0.1-1/Fs);xn=square(2*pi*10*t1);fx=fftshift(fft(x,N)/N;figuref=(-N/2:(N/2-1)*Fs/N;stem(f,abs(fx)title(方波信號(hào)的頻譜)n,Wn=buttord(40/500,50/500,3,15); %確定butterworth濾波器階次b,a=butter(n,Wn); %確定濾波器傳遞函數(shù)figurefreqz(b,a); %濾波器頻率特性title(低通濾波器頻率特性)y=filter(b,a,x); %濾波figureplot(t,y)title(矩形波通過(guò)低通濾波器后的波形)axis(0,0.2,-1.5,1.5)N=512; %濾波后的頻譜fx=fftshift(fft(y,N)*(1/Fs);figuref=(-N/2:(N/2-1);f=f*Fs/N;plot(f,abs(fx)title(濾波后信號(hào)的頻譜)s=randn(1,1001); %信號(hào)疊加噪聲xs=x+s;figureplot(t,xs)title(信號(hào)疊加噪聲的波形)%axis(0,0.5,-2,2)n,Wn=buttord(80/500,100/500,1,15); %確定butterworth濾波器階次b,a=butter(n,Wn); %確定濾波器傳遞函數(shù)figurefreqz(b,a); %濾波器頻率特性title(低通濾波器頻率特性)y1=filter(b,a,xs); %采用filter函數(shù)濾波figureplot(t,y1) %濾波后的波形title(濾波后的波形)%axis(0,0.5,-2,2)參考程序2 FIR濾波器實(shí)現(xiàn)%采用FIR濾波器clcclearFs=1000;T=1/10;t=0:1/Fs:1;x=square(2*pi*10*t);%產(chǎn)生幅度為-1,1的方波figureplot(t,x)title(矩形波)axis(0,0.5,-1.5,1.5)N=T*Fs; %方波的頻譜t1=0:1/Fs:(0.1-1/Fs);xn=square(2*pi*10*t1);fx=fftshift(fft(x,N)/N;figuref=(-N/2:(N/2-1)*Fs/N;stem(f,abs(fx)title(方波信號(hào)的頻譜)n=48win=chebwin(n+1,40);%窗函數(shù)法設(shè)計(jì)FIR濾波器(切比雪夫窗），也可以采用其他窗bt=fir1(n,45/500,win); %截止頻率45hzfigurefreqz(bt,1); %濾波器頻率特性title(低通濾波器頻率特性)y=filter(bt,1,x); %濾波figureplot(t,y)title(矩形波通過(guò)低通濾波器后的波形)axis(0,0.5,-1.5,1.5)N=512; %濾波后的頻譜fx=fftshift(fft(y,N)*(1/Fs);figuref=(-N/2:(N/2-1);f=f*Fs/N;plot(f,abs(fx)title(濾波后信號(hào)的頻譜)s=randn(1,1001); %信號(hào)疊加噪聲xs=x+s;figureplot(t,xs)title(信號(hào)疊加噪聲的波形)axis(0,0.5,min(xs),max(xs)n=48win=chebwin(n+1,40);%窗函數(shù)法設(shè)計(jì)FIR濾波器(切比雪夫窗），也可以采用其他窗bt=fir1(n,80/500,win); %截止頻率80hz（可自己調(diào)整）figurefreqz(bt,1); %濾波器頻率特性title(低通濾波器頻率特性)y1=filter(bt,1,xs); %采用filter函數(shù)濾波figureplot(t,y1) %濾波后的波形title(濾波后的波形)axis(0,0.5,min(y1),max(y1)fxs=fftshift(fft(y1,N)*(1/Fs); %濾波后波形的頻譜figuref=(-N/2:(N/2-1);f=f*Fs/N;plot(f,abs(fxs) title(濾波后波形的頻譜)二參考程序clcclearmusic=wavread(d:matlab6p5toolboxdspblksdspdemosdspafsx_mono.wav);%可根據(jù)自己電腦情況調(diào)整music=music;a=size(music)%利用media player觀察音樂(lè)持續(xù)時(shí)間為3s，信號(hào)的采樣點(diǎn)數(shù)為60000，可知信號(hào)的采樣頻率為20KHzmusic=decimate(music,2);%音樂(lè)的采樣頻率降低為10KHzt=0:0.0001:(3-0.0001);x=sin(2*pi*1000*t)+sin(2*pi*3000*t);mx=music+x;xa=0.1;%顯示x軸上限，可調(diào)整figuresubplot(3,1,1),