数字信号处理实验全部程序MATLAB.docx
- 文档编号:17435457
- 上传时间:2023-04-24
- 格式:DOCX
- 页数:21
- 大小:168.94KB
数字信号处理实验全部程序MATLAB.docx
《数字信号处理实验全部程序MATLAB.docx》由会员分享,可在线阅读,更多相关《数字信号处理实验全部程序MATLAB.docx(21页珍藏版)》请在冰豆网上搜索。
数字信号处理实验全部程序MATLAB
实验一熟悉MATLAB环境
一、实验目的
(1)熟悉MATLAB的主要操作命令。
(2)学会简单的矩阵输入和数据读写。
(3)掌握简单的绘图命令。
(4)用MATLAB编程并学会创建函数。
(5)观察离散系统的频率响应。
二、实验内容
认真阅读本章附录,在MATLAB环境下重新做一遍附录中的例子,体会各条命令的含义。
在熟悉了MATLAB基本命令的基础上,完成以下实验。
上机实验内容:
(1)数组的加、减、乘、除和乘方运算。
输入A=[1234],B=[3456],求C=A+B,D=A-B,E=A.*B,F=A./B,G=A.^B并用stem语句画出A、B、C、D、E、F、G。
实验程序:
A=[1234];
B=[3456];
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');
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,'fill');xlabel('时间序列n');ylabel('A.^B');
运行结果:
(2)用MATLAB实现以下序列。
a)x(n)=0.8n0≤n≤15
实验程序:
n=0:
15;x=0.8.^n;
stem(n,x,'fill');xlabel('时间序列n');ylabel('x(n)=0.8^n');
b)x(n)=e(0.2+3j)n0≤n≤15
实验程序:
n=0:
15;x=exp((0.2+3*j)*n);
stem(n,x,'fill');xlabel('时间序列n');ylabel('x(n)=exp((0.2+3*j)*n)');
运行结果:
a)的时间序列b)的时间序列
c)x(n)=3cos(0.125πn+0.2π)+2sin(0.25πn+0.1π)0≤n≤15
实验程序:
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)');
运行结果:
d)将c)中的x(n)扩展为以16为周期的函数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)扩展为以10为周期的函数x10(n)=x(n+10),绘出四个周期
实验程序:
n=0:
1:
39;
x=3*cos(0.125*pi*rem(n,10)+0.2*pi)+2*sin(0.25*pi*rem(n,10)+0.1*pi);
stem(n,x,'fill');xlabel('时间序列n');ylabel('x10(n)');
运行结果:
d)的时间序列e)的时间序列
(3)x(n)=[1,-1,3,5],产生并绘出下列序列的样本。
a)x1(n)=2x(n+2)-x(n-1)-2x(n)
实验程序:
n=0:
3;
x=[1-135];
x1=circshift(x,[0-2]);x2=circshift(x,[01]);x3=2*x1-x2-2*x;
stem(x3,'fill');xlabel('时间序列n');ylabel('x1(n)=2x(n+2)-x(n-1)-2x(n)');
b)
实验程序:
n=0:
3;
x=[1-135];
x1=circshift(x,[01]);x2=circshift(x,[02]);x3=circshift(x,[03]);
x4=circshift(x,[04]);x5=circshift(x,[05]);
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-4)+5x(n-5)');
运行结果:
a)的时间序列b)的时间序列
(4)绘出时间函数的图形,对x轴、y轴图形上方均须加上适当的标注。
a)x(t)=sin(2πt)0≤t≤10sb)x(t)=cos(100πt)sin(πt)0≤t≤4s
实验程序:
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)');
subplot(2,1,2);
plot(t2,xb);xlabel('t');ylabel('x(t)');title('x(t)=cos(100*pi*t2).*sin(pi*t2)');
运行结果:
(5)编写函数stepshift(n0,n1,n2)实现u(n-n0),n1 实验程序: clc; n1=input('请输入起点: '); n2=input('请输入终点'); n0=input('请输入阶跃位置'); n=n1: n2; x=[n-n0>=0]; stem(n,x,'fill');xlable('时间序列n');ylable('u(n-n0)'); 请输入起点: 2 请输入终点: 8 请输入阶跃位置: 6 运行结果: (5)运行结果(6)运行结果 (6)给一定因果系统 求出并绘制H(z)的幅频响应与相频响应。 实验程序: a=[1-0.670.9]; b=[1sqrt (2)1]; [hw]=freqz(b,a); fp=20*log(abs(h)); subplot(2,1,1); plot(w,fp);xlabel('时间序列t');ylabel('幅频特性'); xp=angle(h); subplot(2,1,2); plot(w,xp);xlabel('时间序列t');ylabel('相频特性'); 运行结果: (右上图) 常用典型序列 单位采样序列 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); fori=1: nc ifn(i)==n0 x(i)=1 end end stem(n,x); title('电信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') 实数指数 n=[0: 10]; x=0.9.^n; stem(n,x); title('电信1201') 复数指数序列 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=angle(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') 例: 求出下列波形 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,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') 已知某一系统方程为: y[n]-y[n-1]+0.9y[n-2]=x[n]计算并画出脉冲响应h(n),n=(-20,100) n=(-20: 100); num= (1);den=[1-10.9]; x=impseq(-20,100,0); h=filter(num,den,x); stem(n,h) xlabel('时间序号N'); ylabel('脉冲响应h'); title('脉冲响应电信1201'); 实现两个序列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') 自己给的数 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) %信号处理的改进卷积程序 nyb=nx (1)+nh (1); 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);stem(ny,y);title('电信1201') 实验二 functionxk=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)' IDFS functionxn=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)是复指数,它满足周期性,我们将在两个周期中的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=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,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*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的正弦信号组成,数据采样频率 fs=100; N=128; n=0: N-1; t=n/fs; x=0.5*sin(2*pi*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; subplot(2,1,1); plot(y); title('电信1201') subplot(2,1,2); plot(f,P(1: 256)); title('电信1201') 对信号进行DFT,对其结果进行IDFT,并将IDFT的结果和原信号进行比较。 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('originalsignal电信1201') y=fft(x,N); mag=abs(y); f=(0: length(y)-1)'*fs/length(y); subplot(2,2,2) plot(f,mag) title('FFTtooriginalsignal电信1201') xifft=ifft(y); magx=real(xifft); ti=[0: length(xifft)-1]/fs; subplot(2,2,3) plot(ti,magx); title('signalfromIFFT电信1201') yif=fft(xifft,N); mag=abs(yif); subplot(2,2,4) plot(f,mag) title('FFTtosignalfromIFFT电信201') 用函数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) 实验四 例: 用冲激响应不变法设计Butterworth低通数字滤波器,通带波纹小于1dB,阻带在 wp=0.2*pi;ws=0.3*pi; rp=1;rs=15;ts=0.01;Nn=128; Wp=wp/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') 设计一个Butterworth高通数字滤波器,通带边界频率为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,Nn,fs) title('电信1201') 设计一个24阶FIR带通滤波器,通带频率为 wp=[0.35,0.65]; N=24; b=fir1(2*N,wp); freqz(b,1,512) title('电信1201') f=0: 0.002: 1; m(1: 201)=1;m(202: 301)=0; m(302: 351)=0.5;m(352: 401)=0;m(402: 501)=1; hold plot(f,m,'r: ') b=fir2(64,f,m); [h,f1]=freqz(b); f1=f1./pi; plot(f1,abs(h)) title('N=64电信1201') f=0: 0.002: 1; m(1: 201)=1;m(202: 301)=0; m(302: 351)=0.5;m(352: 401)=0;m(402: 501)=1; hold plot(f,m,'r: ') b=fir2(256,f,m); [h,f1]=freqz(b); f1=f1./pi; plot(f1,abs(h)) title('N=256电信1201') 针对一个含有5Hz、15Hz和30Hz的混和正弦波信号,设计一个FIR带通滤波器。 fc1=10;fc2=20;fs=100; [n,Wn,beta,ftype]=kaiserord([7131723],[010],[0.010.010.01],100); w1=2*fc1/fs;w2=2*fc2/fs; window=kaiser(n+1,beta);%使用kaiser窗函数 b=fir1(n,[w1w2],window);%使用标准频率响应的 freqz(b,1,512);%数字滤波器频率响应 t=(0: 100)/fs; s=sin(2*pi*t*5)+sin(2*pi*t*15)+sin(2*pi*t*30); sf=filter(b,1,s);%对信号s进行滤波 figure subplot(2,1,1);plot(t,s);title('电信1201') subplot(2,1,2);plot(t,sf) title('电信1201') (注: 文档可能无法思考全面,请浏览后下载,供参考。 )
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数字信号 处理 实验 全部 程序 MATLAB
![提示](https://static.bdocx.com/images/bang_tan.gif)