数字信号处理实验全部程序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'）
（注：
文档可能无法思考全面，请浏览后下载，供参考。
）
升级会员 