基于MATLAB的数字信号处理实验指导书Word格式.doc
《基于MATLAB的数字信号处理实验指导书Word格式.doc》由会员分享，可在线阅读，更多相关《基于MATLAB的数字信号处理实验指导书Word格式.doc（38页珍藏版）》请在冰点文库上搜索。
为克服上述困难，美国Mathwork公司于1967年推出“MatrixLaboratory”（缩写为MATLAB）软件包，并不断更新和扩充。
MATLAB软件包现已成为国际公认的最优秀的科技界应用软件，是一种面向科学和工程计算的高级语言，它强大的计算功能、计算结果的可视化以及极高的编程效率，是其它语言无与伦比的。
MATLAB包含的几十个工具箱，覆盖了通讯、自动控制、信号处理、图象处理等科学领域，汲取了当今世界这些领域的最新研究成果，已经成为从事科学研究和工程设计不可缺少的工具软件。
MATLAB作为一种计算机语言，要想运用自如，充分发挥它的威力，也需要先系统地学习它。
但由于使用MATLAB编程运算与人进行科学计算的思路和表达方完全一致，所以不象学习其他高级语言那样难于掌握。
实践证明，你可以在一个多小时的时间学会MATLAB的基础知识，在短短的几小时的使用中就能初步掌握它。
从而使你可以在短时间内高效地、出色地完成复杂的科学计算、能迅速地测试他们的构想、综合评测系统性能。
所以，在欧美高等院校中，MATLAB已成为大学生、研究生、教师的必备的基本技能。
近年来，国内高校也在大力推广应用MATLAB软件。
1.2MATLAB的基本操作
一、启动与退出
二、命令窗口与M文件编辑窗口
1．命令窗口：
输入一条MATLAB语句，就立即执行。
【例】
a=[1.51.11.3;
2.12.32.5;
3.33.53.1]
a=
1.50001.10001.3000
2.10002.30002.5000
3.30003.50003.1000
inv（a）%矩阵求逆
ans=
1.9565-1.37680.2899
-2.1014-0.43481.2319
0.28991.9565-1.3768
a'
%转置
1.50002.10003.3000
1.10002.30003.5000
1.30002.50003.1000
b=a*a%矩阵相乘
b=
8.85008.73008.7300
16.230016.350016.2300
22.530022.530022.6500
det（a）%a的行列式
-0.8280
eig（a）%a的特征值
6.9000
0.3464
-0.3464
rank（a）%a的秩
3
[LU]=lu（a）%a做KU分解
L=
0.45451.00000
0.6364-0.14811.0000
1.000000
U=
0-0.4909-0.1091
000.5111
但仅靠一条一条输入MATLAB语句，难以实现复杂功能。
为了设计诸如循环、条件分支等功能，MATLAB利用了M文件。
2．M文件的编辑、路径的设定
CommandWindow——File ——New——mfile,打开MATLABEditor/DebuggerWindow，编辑m文件。
a=[123;
456;
789];
b=a'
c=a*a
第二章信号的生成和信号的变换
2.1序列的表示及基本序列的生成
MATLAB是用向量表示序列的。
MATLAB向量的第一个元素位置是x
（1），不是x（0）。
为了清楚表示序列{x（n）}；
要用两个向量，其中一个向量n表示序列的位置，而另一个向量x表示序列值。
如：
n=[-2,-1,0,1,2,3];
x=[0,1,2,3,2,1];
stem（n,x）;
1．单位取样序列
%CREATEDELTASEQUENCE
n0=0;
n1=-10;
n2=10;
n=[n1:
n2];
nc=length（n）;
x=zeros（1,nc）;
fori=1:
nc
ifn（i）==n0
x（i）=1;
end
end
stem（n,x）
xlabel（'
n'
）;
ylabel（'
x（n）'
title（'
DeltaSequence'
grid
2．单位阶跃序列
%CREATESTEPSEQUENCE
x=[（n-n0）>
=0];
stepSequence'
3．正弦序列
%CREATESINESEQUENCE
n=[0:
100];
x=2*sin（0.05*pi*n+pi/4）;
SineSequence'
4．复指数序列
%CREATEAComplexPowerSEQUENCE
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=（180/pi）*angle（x）;
subplot（221）
stem（n,Real_x）;
RealPart'
subplot（222）
stem（n,Image_x）;
ImaginaryPart'
subplot（223）
stem（n,Mag_x）;
Magnitude'
subplot（224）
stem（n,Phase_x）;
Phase'
5．随机序列
rand（1,n）——产生[0，1]上均匀分布的随机序列，长度为n。
randn（1,n）——产生N（0,1）高斯随机序列，即白噪声，长度n。
6．序列的线性卷积
y（n）=x（n）*h（n）
h（n）
y（n）
x（n）
y=conv（x,h）
conv——一维卷积
con2——二维卷积
%MATLABPROGRAMofconvolution
b=0.15;
a=[1-0.8];
%Creatinputoffilter
x1=2*sin（0.05*pi*n）
x=x1+0.4*randn（1,101）;
%Calculateresponseaccordingtoconvolution
imp=[1;
zeros（100,1）];
h=filter（b,a,imp'
yc=conv（h,x）;
y=yc（1:
101）;
%Calculateresponseaccordingtomatlabfunction
y1=filter（b,a,x）;
subplot（311）,plot（n,x1）;
x'
subplot（312）,plot（n,x）;
x+noise'
subplot（313）,plot（n,y）;
y'
2.2序列的DFT、FFT
有限长序列x（n）
离散傅里叶变换X（k）也是一个长度为N的序列，所对应的离散频率在0~之间，且频率间隔相等，为。
可用函数fft和ifft实现离散傅里叶正变换和反变换。
调用格式为：
；
【例】已知序列x（n）=cos（0.4n）+cos（0.6n）,,试绘制x（n）波形及其傅里叶变换的波形。
【注：
】
（，，）
%CalculateDFTofx（n）
N=120;
n=0:
N-1;
xn1=cos（0.4*pi*n）+cos（0.5*pi*n）;
xn=xn1+rand（1,N）;
Xk=fft（xn,N）;
magXk=abs（Xk）;
subplot（211）,plot（n,xn）;
x（n）N=100'
f=n/N;
subplot（212）;
plot（f,magXk）;
frequency'
|X（k）|'
|X（k）|N=100'
2.3用FFT法求线性卷积
对于有限长序列，存在两种卷积方法：
线性卷积和循环卷积。
若x（n）是长度为nx的有限长序列，h（n）是长度为nh的有限长序列，要使N点循环卷积和线性卷积相等，且不发生混叠失真的条件是：
。
而循环卷积可以采用FFT求取。
【例】求序列x1与x2线性卷积
x1=[122];
x2=[1234];
disp（'
LinearConvolution'
）
yl=conv（x1,x2）
CircularConvolutionusingFFTN=4'
N=4;
x1=[x1,zeros（1,N-length（x1））];
x2=[x2,zeros（1,N-length（x2））];
Xk1=fft（x1,N）;
Xk2=fft（x2,N）;
Yk=Xk1.*Xk2;
y=ifft（Yk,N）
CircularConvolutionusingFFTN=6'
N=6;
在CommandWindow中显示的结果如下：
LinearConvolution
yl=14914148
CircularConvolutionusingFFTN=4
y=15.000012.0000-0.0000i9.000014.0000+0.0000i
CircularConvolutionusingFFTN=6
y=1.0000-0.0000i4.0000-0.0000i9.0000-0.0000i14.0000+0.0000i
14.0000+0.0000i8.0000+0.0000i
第三章数字滤波器的结构
3.1直接型——传递函数形式
用分子和分母多项式的系数两个向量来表示。
Num=[b
（1）b
（2）……b（nb+1）]
Den=[a
（1）a
（2）……a（na+1）]
num=[1-311-2718];
den=[16122-4-1];
3.2零极点增益形式
用零点向量z、极点向量p以及增益K表示
[z,p,k]=tf2zp（num,den）
z=0.0000+3.0000i
0.0000-3.0000i
2.0000
1.0000
p=-0.5000+0.5000i
-0.5000-0.5000i
0.5000
-0.2500
k=0.0625
3.3级联型——二阶因子级联形式
用的数组SOS表示
sos=zp2sos（z,p,k）
sos=0.0625-0.18750.12501.0000-0.2500-0.1250
1.0000-0.00009.00001.00001.00000.5000
3.4并联型——部分分式展开式形式
用函数residuez（）把传递函数形式转换为部分分式形式。
[r,p,c]=residuez（num,den）
r=
-5.0250-1.0750i
-5.0250+1.0750i
0.9250
27.1875
p=
-0.5000+0.5000i
-0.5000-0.5000i
0.5000
-0.2500
c=-18
然后将共轭极点两两合并成如下的一阶节和二阶节的并联。
MATLABSignalProcessingToolbox中的系统模型之间的相互转换函数有：
[z,p,k]=tf2zp（b,a）
[b,a]=zp2tf（z,p,k）
[z,p,k]=sos2zp（sos）
[b,a]=sos2tf（sos）
[r,p,k]=residuez（b,a）
[b,a]=residuz（r,p,k）
第四章IIR数字滤波器设计
4.1MATLAB中模拟滤波器设计函数介绍
设计流程图：
给定模拟滤波器的技术指标
采用真实角频率
（弧度/秒）
计算模拟滤波器的阶次N和截止频率Wc。
（利用buttord,cheb1ord,cheb2ord,ellipord等）
设计模拟低通滤波器原型G（p）。
（利用buttap,cheb1ap,cheb2ap,ellipap等）
模拟滤波器完全设计函数H（s）
（butter,cheby1,cheby2,ellip）
由模拟低通原型经频率变换获得所需要的低通、高通、带通和带阻滤波器H（s）。
（利用lp2lp,lp2hp,lp2bp,lp2bs）
1.阶次计算函数；
（1）Butterworth模拟滤波器
[N,Wn]=buttord（Wp,Ws,Rp,Rs,’s’）；
（2）ChebyshevI型模拟滤波器
[N,Wn]=cheb1ord（Wp,Ws,Rp,Rs,’s’）
（3）ChebyshevII型模拟滤波器
[N,Wn]=cheb2ord（Wp,Ws,Rp,Rs,’s’）
（4）椭圆滤波器
[N,Wn]=ellipord（Wp,Ws,Rp,Rs,’s’）
其中：
Wp,Ws分别是通带和阻带的截止频率，单位为弧度/秒。
对低通和高通，Wp,Ws都是标量，对带通和带阻，Wp,Ws是1×
2的向量。
Rp,Rs分别是通带和阻带的衰减（dB）。
N是求出的相应低通滤波器的阶次，Wn是求出的3dB频率，它和Wp稍有不同。
2．模拟低通原型滤波器设计函数（注意：
得到的是归一化的传递函数）
[z,p,k]=buttap（N）；
滤波器传递函数具有如下形式：
[z,p,k]=cheb1ap（N,Rp）
[z,p,k]=cheb2ap（N,Rs）
[z,p,k]=ellipap（N,Rp,Rs）
其中，N是欲设计的低通原型滤波器的阶次，z,p,k是设计出的极点、零点及增益。
3．频率变换函数
将模拟低通原型转换为实际的低通、高通、带通及带阻滤波器。
（1）[B,A]=lp2lp（b,a,wo）,
（2）[B,A]=lp2hp（b,a,wo）
（3）[B,A]=lp2bp（b,a,wo,Bw）,
（4）[B,A]=lp2bs（b,a,wo,Bw）
b,a是AFLP的分子、分母的系数向量，B,A是转换后的的分子、分母的系数向量；
在
（1）、
（2）中，wo是低通或高通滤波器的截止频率；
在（3）、（4）中，wo是带通或带阻滤波器的中心频率，Bw是其带宽。
，
4．模拟滤波器完全设计函数（=2+3）
[b,a]=butter（N,Wn,’s’）%低通或带通滤波器
[b,a]=butter（N,Wn,’high’,’s’）%高通滤波器
[b,a]=butter（N,Wn,’stop’,’s’）%带阻滤波器
（2）ChebyshevI型模拟滤波器
[b,a]=cheby1（N,Rp,Wn,’s’）%低通或带通滤波器
[b,a]=cheby1（N,Rp,Wn,’high’,’s’）%高通滤波器
[b,a]=cheby1（N,Rp,Wn,’stop’,’s’）%带阻滤波器
[b,a]=cheby2（N,Rs,Wn,’s’）%低通或带通滤波
[b,a]=cheby2（N,Rs,Wn,’high’,’s’）%高通滤波器
[b,a]=cheby2（N,Rs,Wn,’stop’,’s’）%带阻滤波器
（4）椭圆滤波器
[b,a]=ellip（N,Rp,Rs,Wn,’s’）%低通或带通滤波器
[b,a]=ellip（N,Rp,Rs,Wn,’high’,’s’）%高通滤波器
[b,a]=ellip（N,Rp,Rs,Wn,’stop’,’s’）%带阻滤波器
【例1】设计一个Butterworth模拟低通滤波器，通带截至频率fp=5KHz，阻带截至频率fs=10KHz；
通带最大衰减Ap=3，阻带最小衰减As=30。
画出该滤波器的频率响应。
%DesignaButterworthAnaloglowpassfilter
%Specifications
Wp=5000*2*pi;
Ws=10000*2*pi;
Ap=3;
As=30;
%ComputeorderandCutofffrequency
[N,Wc]=buttord（Wp,Ws,Ap,As,'
s'
%designanaloglowpassfilterprototype
[z,p,k]=buttap（N）%归一化的传递函数G（p）
[b,a]=zp2tf（z,p,k）;
%frequencytransform,得到反归一化传递函数G（s）
[bt,at]=lp2lp（b,a,Wp）;
w=0:
20*2*pi:
20000*2*pi;
H=freqs（bt,at,w）;
subplot（211）,plot（w/（2*pi）,abs（H））;
Frequency（Hz）'
ButterworthAnaloglowpassfilter'
%Usingfilterwholedesignfuction
[bb,aa]=butter（N,Wc,'
HH=freqs（bb,aa,w）;
subplot（212）,plot（w/（2*pi）,abs（HH））;
【例2】设计一个Chebyshev模拟带通滤波器，通带截至频率分别为1000Hz和2000Hz，阻带截至频率分别为500Hz和2500Hz；
通带最大衰减Ap=1，阻带最小衰减As=100。
%DesignaChebyshevAnalogbandpassfilter
Wp=[10002000]*2*pi;
Ws=[50
升级会员 