数字信号处理实验全部程序MATLAB.docx
《数字信号处理实验全部程序MATLAB.docx》由会员分享,可在线阅读,更多相关《数字信号处理实验全部程序MATLAB.docx(19页珍藏版)》请在冰豆网上搜索。
数字信号处理实验全部程序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')