數(shù)字信號處理實驗全部程序MATLAB

1、實驗一 熟悉MATLAB環(huán)境一、實驗?zāi)康?1)熟悉MATLAB的主要操作命令。(2)學(xué)會簡單的矩陣輸入和數(shù)據(jù)讀寫。(3)掌握簡單的繪圖命令。(4)用MATLAB編程并學(xué)會創(chuàng)建函數(shù)。(5)觀察離散系統(tǒng)的頻率響應(yīng)。二、實驗內(nèi)容認(rèn)真閱讀本章附錄，在MATLAB環(huán)境下重新做一遍附錄中的例子，體會各條命令的含義。在熟悉了MATLAB基本命令的基礎(chǔ)上，完成以下實驗。上機實驗內(nèi)容：(1)數(shù)組的加、減、乘、除和乘方運算。輸入A=1 2 3 4，B=3 4 5 6，求C=A+B，D=A-B，E=A.*B，F(xiàn)=A./B，G=A.B并用stem語句畫出A、B、C、D、E、F、G。實驗程序：A=1 2 3 4;B=3

2、 4 5 6;n=1:4;C=A+B;D=A-B;E=A.*B;F=A./B;G=A.B;subplot(4,2,1);stem(n,A,'fill');xlabel ('時間序列n');ylabel('A');subplot(4,2,2);stem(n,B,'fill');xlabel ('時間序列n ');ylabel('B');subplot(4,2,3);stem(n,C,'fill');xlabel ('時間序列n ');ylabel('A+B

3、9;);subplot(4,2,4);stem(n,D,'fill');xlabel ('時間序列n ');ylabel('A-B');subplot(4,2,5);stem(n,E,'fill');xlabel ('時間序列n ');ylabel('A.*B');subplot(4,2,6);stem(n,F,'fill');xlabel ('時間序列n ');ylabel('A./B');subplot(4,2,7);stem(n,G,'f

4、ill');xlabel ('時間序列n ');ylabel('A.B');運行結(jié)果：(2)用MATLAB實現(xiàn)以下序列。a）x(n)=0.8n 0n15實驗程序：n=0:15;x=0.8.n;stem(n,x,'fill'); xlabel ('時間序列n ');ylabel('x(n)=0.8n');b）x(n)=e(0.2+3j)n 0n15實驗程序：n=0:15;x=exp(0.2+3*j)*n);stem(n,x,'fill'); xlabel ('時間序列n ');

5、ylabel('x(n)=exp(0.2+3*j)*n)');運行結(jié)果：a）的時間序列 b）的時間序列c）x(n)=3cos(0.125n+0.2)+2sin(0.25n+0.1) 0n15實驗程序：n=0:1:15;x=3*cos(0.125*pi*n+0.2*pi)+2*sin(0.25*pi*n+0.1*pi);stem(n,x,'fill'); xlabel('時間序列n ');ylabel('x(n)=3*cos(0.125*pi*n+0.2*pi)+2*sin(0.25*pi*n+0.1*pi)');運行結(jié)果：d）將c

6、)中的x(n)擴(kuò)展為以16為周期的函數(shù)x16(n)=x(n+16),繪出四個周期實驗程序：n=0:1:63;x=3*cos(0.125*pi*rem(n,16)+0.2*pi)+2*sin(0.25*pi*rem(n,16)+0.1*pi);stem(n,x,'fill'); xlabel ('時間序列n ');ylabel('x16(n)');e）將c)中的x(n)擴(kuò)展為以10為周期的函數(shù)x10(n)=x(n+10),繪出四個周期實驗程序：n=0:1:39;x=3*cos(0.125*pi*rem(n,10)+0.2*pi)+2*sin(0.2

7、5*pi*rem(n,10)+0.1*pi);stem(n,x,'fill'); xlabel ('時間序列n ');ylabel('x10(n)');運行結(jié)果：d）的時間序列 e）的時間序列(3)x(n)=1,-1,3,5，產(chǎn)生并繪出下列序列的樣本。a）x1(n)=2x(n+2)-x(n-1)-2x(n)實驗程序：n=0:3;x=1 -1 3 5;x1=circshift(x,0 -2);x2=circshift(x,0 1);x3=2*x1-x2-2*x;stem(x3,'fill'); xlabel ('時間序列n

8、');ylabel('x1(n)=2x(n+2)-x(n-1)-2x(n)');b）實驗程序：n=0:3;x=1 -1 3 5;x1=circshift(x,0 1);x2=circshift(x,0 2);x3=circshift(x,0 3);x4=circshift(x,0 4);x5=circshift(x,0 5);xn=1*x1+2*x2+3*x3+4*x4+5*x5;stem(xn,'fill'); xlabel('時間序列n ');ylabel('x2(n)=x(n-1)+2x(n-2)+3x(n-3)+4x(n-

9、4)+5x(n-5)');運行結(jié)果： a）的時間序列b）的時間序列(4)繪出時間函數(shù)的圖形，對x軸、y軸圖形上方均須加上適當(dāng)?shù)臉?biāo)注。 a) x(t)=sin(2t) 0t10s b) x(t)=cos(100t)sin(t) 0t4s實驗程序：clc;t1=0:0.001:10;t2=0:0.01:4;xa=sin(2*pi*t1);xb=cos(100*pi*t2).*sin(pi*t2);subplot(2,1,1);plot(t1,xa);xlabel ('t');ylabel('x(t)');title('x(t)=sin(2*pi*t)

10、');subplot(2,1,2);plot(t2,xb);xlabel ('t');ylabel('x(t)');title('x(t)=cos(100*pi*t2).*sin(pi*t2)');運行結(jié)果：(5)編寫函數(shù)stepshift(n0,n1,n2)實現(xiàn)u(n-n0),n1<n0<n2,繪出該函數(shù)的圖形，起點為n1，終點為n2。實驗程序：clc;n1=input('請輸入起點：');n2=input('請輸入終點');n0=input('請輸入階躍位置');n=n1:n

11、2;x=n-n0>=0;stem(n,x,'fill');xlable('時間序列n');ylable('u(n-n0)');請輸入起點：2請輸入終點：8請輸入階躍位置：6運行結(jié)果： （5）運行結(jié)果 （6）運行結(jié)果(6)給一定因果系統(tǒng)求出并繪制H(z)的幅頻響應(yīng)與相頻響應(yīng)。實驗程序：a=1 -0.67 0.9;b=1 sqrt(2) 1;h w=freqz(b,a);fp=20*log(abs(h);subplot(2,1,1);plot(w,fp);xlabel('時間序列t');ylabel('幅頻特性'

12、);xp=angle(h);subplot(2,1,2);plot(w,xp);xlabel('時間序列t');ylabel('相頻特性');運行結(jié)果：(右上圖)常用典型序列單位采樣序列function x,n=impseq(n1,n2,n0) n=n1:n2; x=(n-n0)=0;x,n=impseq(-2,8,2);stem(n,x);title('電信1201')n0=-2;n=-10:10;nc=length(n);x=zeros(1,nc);for i=1:nc if n(i)=n0 x(i)=1 endendstem(n,x) ;t

13、itle('電信1201 采樣序列第二種方法')單位階躍序列function x,n=stepseq(n1,n2,n0) n=n1:n2; x=(n-n0)>=0;x,n=stepseq(-2,8,2);stem(n,x);title('電信1201')實數(shù)指數(shù)n=0:10;x=0.9.n;stem(n,x);title('電信1201')復(fù)數(shù)指數(shù)序列n=-10:10;alpha=-0.1+0.3*j;x=exp(alpha*n);real_x=real(x); image_x=imag(x);mag_x=abs(x); phase_x=a

14、ngle(x);subplot(2,2,1); stem(n,real_x) ;title('電信1201')subplot(2,2,2); stem(n,image_x);title('電信1201')subplot(2,2,3); stem(n,mag_x);title('電信1201')subplot(2,2,4); stem(n,phase_x);title('電信1201')正余弦n=0:10;x=3*cos(0.1*pi*n+pi/3);stem(n,x);title('電信1201')例： 求出下列波

15、形x1(n)=2x(n-5)-3x(n+4)function y,n=sigadd(x1,n1,x2,n2) n=min(min(n1),min(n2) : max(max(n1),max(n2); y1=zeros(1,length(n); y2=y1; y1(find(n>=min(n1)&(n<=max(n1)=1)=x1; y2(find(n>=min(n2)&(n<=max(n2)=1)=x2; y=y1+y2;function y,n=sigshift(x,m,n0) n=m+n0; y=x; n=-2:10;x=1:7,6:-1:1;x11

16、,n11=sigshift(x,n,5);x12,n12=sigshift(x,n,-4);x1,n1=sigadd(2*x11,n11,-3*x12,n12);stem(n1,x1);title('電信1201')已知某一系統(tǒng)方程為：yn-yn-1+0.9yn-2=xn計算并畫出脈沖響應(yīng)h(n),n=(-20,100)n=(-20:100);num=(1); den=1 -1 0.9;x=impseq(-20,100,0);h=filter(num,den,x);stem(n,h)xlabel('時間序號N');ylabel('脈沖響應(yīng)h');

17、title('脈沖響應(yīng) 電信1201');實現(xiàn)兩個序列a,b的卷積x=3,11,7,0,-1,4,2; h=2,3,0,-5,2,1; c=conv(x,h); stem(c); title('電信1201')自己給的數(shù)x=0.5,1,1.5,0,0; h=1,1,1;c=conv(x,h);stem(c);title('電信1201')試求卷積C(t)=f1(t)*f2(t),并繪制出f1、f2、及卷積以后的波形。function y,ny= conv_m(x,nx,h,nh,p) % 信號處理的改進(jìn)卷積程序 nyb=nx(1)+nh(1);

18、nyc=nx(length(x)+nh(length(h); ny=(nyb :p: nyc); y=conv(x , h);p=0.1;t1= 0:p:1;f1=t1.*(t1>0);t2= -1:0.1:2;f2=t2.*exp(t2).*(t2>=0)+exp(t2).*(t2<0);y,ny=conv_m(f1,t1,f2,t2,p);subplot(3,1,1);stem(t1,f1);title('電信1201')subplot(3,1,2);stem(t2,f2);title('電信1201')subplot(3,1,3); st

19、em(ny,y);title('電信1201')實驗二function xk=dfs(xn,N) n=(0:1:N-1); k=n; WN=exp(-j*2*pi/N); nk=n'*k; WNnk=WN.nk; xk=xn* WNnk;xn=0,1,2,3;N=4;xk=dfs(xn,N)' IDFSfunction xn=idfs(xk,N) n=(0:1:N-1); k=n; WN=exp(-j*2*pi/N); nk=n'*k; WNnk=WN.(-nk); xn=xk*WNnk/N;分析：因為x(n)是復(fù)指數(shù)，它滿足周期性，我們將在兩個周期中的

20、401個頻點上作計算來觀n=0:10;x=(0.9*exp(j*pi/3).n;k=-200:200;w=(pi/100)*k;X=x*(exp(-j*pi/100).(n'*k);magX=abs(X); angX=angle(X);subplot(2,1,1); plot(w/pi,magX); subplot(2,1,2); plot(w/pi,angX/pi); title('電信1201')檢驗頻移特性n=0:100; x=cos(pi*n/2);k=-100:100; w=(pi/100)*k;X=x*(exp(-j*pi/100).(n'*k);y

21、=exp(j*pi*n/4).*x;Y=y*(exp(-j*pi/100).(n'*k);subplot(2,2,1);plot(w/pi,abs(X); axis(-1,1,0,60);title('幅值 電信1201');xlabel('以pi為單位的頻率');ylabel('絕對值X');subplot(2,2,2);plot(w/pi,angle(X)/pi);title('相位 電信1201');xlabel('以pi為單位的頻率');ylabel('絕對值X');axis(-1,

22、1,-1,1);subplot(2,2,3);plot(w/pi,abs(Y);title('幅值 電信201');xlabel('以pi為單位的頻率');ylabel('絕對值Y');axis(-1,1,0,60);subplot(2,2,4);plot(w/pi,angle(Y)/pi); title('相位 電信201');xlabel('以pi為單位的頻率');ylabel('絕對值Y');axis(-1,1,-1,1);b=1;a=1,-0.8;n=0:100;x=cos(0.05*pi*

23、n);y=filter(b,a,x);subplot(2,1,1); stem(n,x)xlabel('n');ylabel('x(n)');title('輸入序列 電信1201');subplot(2,1,2); stem(n,y)xlabel('n');ylabel('y(n)');title('輸出序列 電信1201');實驗三例：已知信號由15Hz幅值0.5的正弦信號和40Hz幅值2的正弦信號組成，數(shù)據(jù)采樣頻率fs=100;N=128;n=0:N-1;t=n/fs;x=0.5*sin(2*p

24、i*15*t)+2*sin(2*pi*40*t);y=fft(x,N);f=(0:length(y)-1)'*fs/length(y);mag=abs(y);stem(f,mag);title('N=128點 電信1201')已知帶有測量噪聲信號f1=50Hz,f2=1為均值為零、方差為1的隨機信號，采樣頻率為1000Hz，t=0:0.001:0.6;x=sin(2*pi*50*t)+sin(2*pi*120*t);y=x+2*randn(1,length(t);Y=fft(y,512);P=Y.*conj(Y)/512; %求功率f=1000*(0:255)/512;

25、subplot(2,1,1);plot(y);title('電信1201')subplot(2,1,2);plot(f,P(1:256);title('電信1201')對信號 進(jìn)行DFT，對其結(jié)果進(jìn)行IDFT，并將IDFT的結(jié)果和原信號進(jìn)行比較。 fs=100; N=128; n=0:N-1; t=n/fs;x=sin(2*pi*40*t)+sin(2*pi*15*t);subplot(2,2,1)plot(t,x)title('original signal 電信1201')y=fft(x,N);mag=abs(y);f=(0:length(y

26、)-1)'*fs/length(y);subplot(2,2,2)plot(f,mag)title('FFT to original signal 電信1201')xifft=ifft(y);magx=real(xifft);ti=0:length(xifft)-1/fs;subplot(2,2,3)plot(ti,magx);title('signal from IFFT 電信1201')yif=fft(xifft,N);mag=abs(yif);subplot(2,2,4)plot(f,mag)title('FFT to signal fro

27、m IFFT 電信201')用函數(shù)conv和FFT計算同一序列的卷積，比較其計算時間。L=5000; N=L*2-1; n=1:L;x1=0.5*n; x2=2*n;t0=clock; yc=conv(x1,x2);conv_time=etime(clock,t0)t0=clock;yf=ifft(fft(x1,N).*fft(x2,N);fft_time=etime(clock,t0)實驗四例：用沖激響應(yīng)不變法設(shè)計Butterworth低通數(shù)字濾波器，通帶波紋小于1dB，阻帶在 wp=0.2*pi; ws=0.3*pi;rp=1;rs=15;ts=0.01;Nn=128;Wp=wp/

28、ts; Ws=ws/ts;N,Wn=buttord(Wp,Ws,rp,rs,'s');z,p,k=buttap(N);Bap,Aap=zp2tf(z,p,k);b,a=lp2lp(Bap,Aap,Wn);bz,az=impinvar(b,a,1/ts);freqz(bz,az,Nn,1/ts)title('電信1201')設(shè)計一個Butterworth高通數(shù)字濾波器，通帶邊界頻率為300Hz，阻帶邊界頻率為200Hz，通帶波紋小于1dB，阻帶衰減大于20dB，采樣頻率為1000Hz。fs=1000;wp=300/(fs/2); ws=200/(fs/2);rp=1; rs=15; Nn=128;N,Wn=buttord(wp,ws,rp,rs);b,a=butter(N,Wn,'high');freqz(b,a,