数字信号处理matlab函数使用详细

11离散信号和系统分析的MATLAB实现1.9.1利用MATLAB产生离散信号用MATLAB表示一离散序列x[k]时，可用两个向量来表示。其中一个向量表示自变量k的取值范围，另一个向量表示序列x[k]的值。例如序列x[k]={2,1,1,-1,3,0,2}可用MATLAB表示为K=-2:4;x=[2,1,1,-1,3,0,2]可用stem(k,x)（k和X大小写无关系）画出序列波形。K个数必须和x个数相同图标：xlabel（‘要标的数’）；ylabel。。。。当序列是从k=0开始时，可以只用一个向量x来表示序列（stem中一个参数x就够了）。由于计算机内寸的限制，MATLAB无法表示一个无穷长的序列。例1-38利用MATLAB计算单位脉冲序44]2[kk在范围内各点的取值。解：%progran1_1产生单位脉冲序列Ks=-4;ke=4;n=2;K=[ks:ke];X=[(k-n)==0];Stem(k,x):（；也可）xlabel(‘k’);程序产生的序列波形如图1-49所示。例1-39利用MATLAB画出信号X[k]=10sin(0.02k)+n[k],1000k的波形。其中n[k]表示为均值为0方差为1的Gauss分布随机信号。解：2MALAB提供了两个产生（伪）随机序列的函数。Rand(1,N)产生1行N列的[0，1]均匀分布随机数。Randn(1,N)产生1行N列均值为0方差为1的Gauss分布随机数。%program1_2产生受噪声干扰的正弦信号N=100;k=0:N;X=10*sin(0.02*pi*k)+randn(1,N+1);Plot(k,x);Xlabel(‘k’);Ylabel(‘x[k]’);%大小写不影响程序程序产生序列如图1-50所示。1.9.2离散卷积的计算离散卷积是数字信号处理中的一个基本运算，MTLAB提供的计算两个离散序列卷积的函数是conv，其调用方式为y=conv(x,h)其中调用参数x,h为卷积运算所需的两个序列，返回值y是卷积结果。MATLAB函数conv的返回值y中只有卷积的结果，没有y的取值范围。由离散序列卷积的性质可知，当序列x和h的起始点都为k=0时，y的取值范围为k=0至length(x)+length(h)-2。当序列x或（和）h的起始点不是k=0时，由例1-3知，y的非零值范围可根据例1-4的结论进行计算。例1-40试用MATLAB函数conv计算例1-2中序列的卷积。解：program1_3计算离散卷积x=[-0.5,0,0.5,1];%序列x的值kx=-1:2;%序列x的取值范围3h=[1,1,1];kh=-2:0;y=conv(x,h);%计算卷积k=kx(1)+kh(1):kx(end)+kh(end);%计算y的取值范围若不加此句则当做kx和kh从0开始stem(k,y);xlabel(‘k’);ylabel(‘y’);程序的运行结果如图1-51所示。1.9.3离散LTI系统响应MATLAB求解许多离散LTI都可用如下的线性常系数的差分方程描述NnMnnkbnxnkany00][][(1-151)其中x[k]、y[k]分别系统的输入和输出。在已知差分方程的N个初始状态y[k],1kN和输入x[k],就可由下式迭代计算出系统的输出。y[k]=-NnMnnkxbbnnkyaan10][)0/(][)0/(利用MATLAB提供的filter函数，可方便地计算出上述差分方程的零状态响应。filter函数调用形式为y=filter(b,a,x)a对应y，b对应x其中a=[a0,a1,…,aN],b=[bo,b1,…,bM]分别表示差分方程系数。X表示输入序列，y表示输出序列。输出序列的长度和序列相同。例1-40试用M=9点滑动平均系统4y[k]=M110][Mnnkx处理例1-39中的受噪声干扰的正弦信号。解：由式（1-151）可知，M点滑动平均系统可看成是N=0的差分方程。在调用filter函数时，调用参数a=1,b为有M个元素的向量，b中每个元素的值均为1/M。%program1_4滑动平均去噪M=9;N=100;k=0:N;s=10*sin(0.02*pi*k);x=s+randn(1,N+1);b=ones(M,1)/M;a=1;y=filter(b,a,x);plot(k,y,s,’:’);%改成（k，y，k，s，‘：’用虚线表示）不能用stemxlabel(‘k’);ylabel（‘幅度’）legend(‘y[k]’,’s[k]’);%图标说明程序的运行结果如图1-52所示。图中的虚线表示未受噪声干扰的原信号s[k],实线为9点滑动的响应y[k]。比较图1-50的信号可以看出，系统滤出了受干扰信号中的大部分噪声，输出信号相对输入信号有4个样本的延迟。1.9.4利用MATLAB计算DTFT当序列的DTET可写成je的有理多项式时，可用MATLAB信号处理工具箱提供的freqz函数计算DTFT的抽样值。另外，可用MATLAB提供的abs、angle、real、imag等基本函数计算DTFT的幅度、相位、实部、虚部。若X（je）可表示为5X（je）=)()(jjeAeB=NjjMjjaNeeaabMeebbo...10...1则freqz的调用形式为X=freqz(b,a,w)(1-153)式（1-153）中的b和a分别是表示式（1-152）中分子多项式和分母多项式系数的向量，即b=[b0,b1,…,bM]a=[a0,a1,…,aN]w为抽样的频率点，在以式（1-153）形式调用freqz函数时，向量w的长度至少为2。返回值X就是DTFT在抽样点w上的值。注意一般情况下，函数freqz的返回值X是复数。例1-41已知离散系统的H（z）为H（z）=212121214167.08519.01)(4167.08519.01()5320.00640.15320.0)(5009.00019.15009.0(zzzzzzzz试画出该系统的幅度响应。解：%program1_5计算系统幅度响应b1=[0.5009-1.00190.5009];b2=[0.53201.06400.5320];a1=[1.0000-0.85190.4167];a2=[1.00000.85190.4167];b=conv(b1,b2);%卷积计算系数！！！a=conv(a1,a2);%计算H（z）的分母多项式w=linspace(0,pi,512);%默认是512点，常常在0到pi间6H=freqz(b,a,w);plot(w/pi,abs(H));ylabel(‘幅度’);xlabel(‘Normalizedfrequency’);系统幅度响应如图1-531.9.5部分分式法的MATLAB实现MATLAB的信号处理工具箱提供了一个计算X（z）的部分分式展开的函数，它的调用形式如下[r,p,k]=resifuez(b,a)其中调用参数b,a分别表示用z1表示X（z）的分子和分母多项式。如果X（z）的部分分式展开为X(z)=12111121)31(4313212111zkkzprzprzprzpr则residuez的返回参数r,p,k分别为r=[r1r2r3r4]p=[p1p2p3p3]k=[k1k2]同一极点p3在向量p中出现了两次，表示p3是一个二阶的重极点。Residuez也用于由r,p,k计算z1表示的X（z）的分子和分母多项式，其调用形式为[b,a]=residuez(r,p,k)例1-43试用MATLAB计算X（z）=6.0072.06.04.11144.02.1608.298.05.13214321zzzzzzzzz的部分分式展开。7解：计算部分分式展开的MATLAB程序如下%program1_6部分分式展开b=[1.5,0.98,-2.608,1.2,-0.144];a=[1,-1.4,0.6,-0.072];[r,p,k]=residuez(b,a);disp(‘留书’)；disp(r’)要直接显示的菜用‘’引起来disp(‘极点’)；disp(p’)disp(‘常数’)；disp(k’);程序运行结果为留数0.70000.50000.3000极点0.60000.60000.2000常数02部分分式展开的结果为X(z)=112112.013.0)6.01(5.06.017.0zzzzz反变换的结果为x[k]=(0.7x0.6k+0.5x(k+1)0.6k+0.3x0.2k)u[k]+2[k-1]1.9.6用MATLAB计算系统的零极点MATLAB信号处理工具箱提供的tf2zp和zp2tf函数可进行系统函数的不同表示形式的转换。z正幂有理多项式表示的系统函数为8)(...)2()1()(...)2()1()(11NazazaMbzbzbzHNNMM(1-154)用零点、极点和常数表示的一阶因子形式的系统函数为))())...(2())(1(())())...(2())(1(()(NpzpzpzMzzzzzzkzH（1-155）MATLAB函数[z,p,k]=tf2zp(b,a)将有理多项式表示的系统函数转换为一阶因子形式的系统函数。[b,a]=zp2tf(z,p,k)将一阶因子形式的系统函数转换为有理多项式表示的系统函数。例1-44试将下面的系统函数表示为一阶因子形式。3212128.016.08.0104.01)(zzzzzH（1-156）解：用z正幂有理多项式表示的系统函数为128.016.08.004.0)(233zzzzzzH%program1_7确定一阶因子形式的系统函数b=[100.040];%此处的0千万别少，少了阶数就少了，一定要考虑到z正幂的常数项即z负幂的最低次项以及所缺的项！！a=[1-0.80.16-0.128];[z,p,k]=tf2zp(b,a);disp(‘零点’);disp(z’);%此处有无分号都会显示disp(‘极点‘)；disp(p’);disp(‘常数’)；disp(k’);程序运行结果为零点00+0.2000i0-0.2000i9极点0.80000.0000+0.4000i0.0000-0.4000i常数1系统函数的一阶因子形式为)4.0)(4.0)(8.0()2.0)(2.0()(jzjzzjzjzzzH（1-157）为了在H（z）的表达式中不出现复数，也可将式（1-157）等价的写成二阶因子的形式)16.01)(8.01(04.01)(212zzzzH（1-158）当H（z）是表示为式（1-154）的形式时，可用[z,p,k]=tf2zp(b,a)求出系统的零极点，从而可根据极点的位置判断系统的稳定性，也可用函数zplane(b,a)（用于给出参数直接根据154式子画图！）画出z平面的零极点分布来判断系统的稳定性。例１-45已知一因果的ＬＴＩ系统的函数为3.0005.05.012)(232zzzzzzH试用函数zplane画出系统的零极点分布，并判断系统的稳定性。解：%program1_8系统零极点分布b=[121];a=[1-0.5-0.0050.3];zplane(b,a);图１－５４画出了系统零极点分布。图中符号表示零点。符号旁的数字表示零点的阶数，符号x表示极点，图中的虚线画的是单位圆。由图可知，该系统的极点全在单10位圆内，故系统是稳定的。１.９.７简单数字滤波器的设计例1-46已知信号x[k]中含有角频率分别为015.01和4.02的离散正弦信号。试设计一高通滤波器，滤出信号x[k]中低频分量，保留信号x[k]高频分量。解：为了简单起见，假设数字滤波器是一个具有如下形式的长度为３的ＦＩＲ系统h