山东大学MATLAB编程指导课件第7章 matlab 信号处理函数

2022/5/191Zhongguo Liu_Biomedical Engineering_Shandong Univ.Biomedical Signal processingmatlab信号处理函数信号处理函数ZhongguoLiuBiomedicalEngineeringSchoolofControlScienceandEngineering,ShandongUniversity2MATLAB是美国MathWorks公司开发的一种功能极其强大的高技术计算语言和内容极其丰富的软件库。它以矩阵和向量的运算以及运算结果的可视化为基础，把广泛应用于各个学科领域的数值分析、矩阵计算、函数生成、信号、图形及图象处理、建模与仿真等诸多强大功能集成在一个便于用户使用的交互式环境之中，为使用者提供了一个高效的编程工具及丰富的算法资源。关于MATLAB3MATLAB与信号处理直接有关的工具箱与信号处理直接有关的工具箱（Toolbox）SignalProcessing(信号处理工具箱）Wavelet(小波工具箱)ImageProcessing（图象处理工具箱）Higher-OrderSpectralAnalysis（高阶谱分析工具箱）4与信号处理间接有关的工具箱：与信号处理间接有关的工具箱：ControlSystem(控制系统)Communication（通信）SystemIdentification（系统辨识）Statistics（统计）NeuralNetwork(神经网络)5例例：z=peaks; surf(z)；61. rand.m用来产生均值为0.5、幅度在01之 间 均 匀 分 布 的 伪 白 噪 声 :u=rand(N,1)(rand(N)生成N阶矩阵)方差：如何改变 的方差与第二章内容有关的MATLAB文件方差函数var(u)标准差函数std(u)71. rand.m用来产生均值为0.5、幅度在01之 间 均 匀 分 布 的 伪 白 噪 声 :u=rand(N,1)(rand(N)生成N阶矩阵)2.randn.m 用来产生均值为零、方差为1服从高斯（正态）分布的白噪声信号u=randn(1,N)与第二章内容有关的MATLAB文件x=randn(1000,1)y=randn(1000,1)v=var(x)h=std(y)83.sinc ：用来产生“sinc”函数：sinc函数定义为： t=-4:0.1:4;x4=sinc(t); %产生抽样函数plot(t,x4) 94. conv.m用来实现两个离散序列的线性卷积。其调用格式是：y=conv(x,h).若x(n)和y(n)的长度分别为M和N， 则返回值是长度为M+N-1的序列。例例 x(n)=345;h(n)=2678，求其线性卷积。MATLAB语句如下： x=3 4 5; h=2 6 7 8; y=conv(x,h)运行结果:y=62655826740u两序列的相关运算两序列的相关运算MATLAB实现：实现：y=xcorr(x1,x2)。x=3 4 5;h=2 6 7 8;y=xcorr(x,h)y =24 53 86 65 38 10 -010 5xcorr: 其互相关和自相关。格式是：(1)rxy=xcorr(x,y)： 求 x,y的 互 相 关 ； (2)rx=xcorr(x,M,flag）：求x的自相关，M：rx的单边长度，总长度为2M+1；flag是定标标志，若flag=biased，则表示是“有偏”估 计 ， 需 将 rx(m)都 除 以 N， 若flag=unbiased,则表示是“无偏”估计，需将rx(m)都除以（Nabs(m)）；若flag缺省，则rx不定标。M和flag同样适用于求互相关。11第三章 Z变换. 在在MATLAB语言中有专门对信号进行正反语言中有专门对信号进行正反Z变换的函数变换的函数ztrans()和和itrans()。其调用格式分别如下：。其调用格式分别如下：uF=ztrans(f)对对f(n)进行进行Z变换，其结果为变换，其结果为F(z)uF=ztrans(f,v)对对f(n)进行进行Z变换，其结果为变换，其结果为F(v)uF=ztrans(f,u,v)对对f(u)进行进行Z变换变换,其结果为其结果为F(v)uf=itrans(F)对对F(z)进行进行Z反变换反变换,其结果为其结果为f(n)uf=itrans(F,u)对对F(z)进行进行Z反变换，其结果为反变换，其结果为f(u)uf=itrans(F,v,u)对对F(v)进行进行Z反变换反变换,其结果为其结果为f(u)u注意：注意：在调用函数在调用函数ztran()及及iztran()之前，要用之前，要用syms命令对所有需要用到的变量（如命令对所有需要用到的变量（如t,u,v,w）等进）等进行说明，即要将这些变量说明成符号变量行说明，即要将这些变量说明成符号变量 12Z变换例例.求数列求数列 fn=e-n的的Z变换及其逆变换。命令如下：变换及其逆变换。命令如下：syms n zfn=exp(-n);Fz=ztrans(fn,n,z) %求求fn的的Z变换变换f=iztrans(Fz,z,n) %求求Fz的逆的逆Z变换变换u例例用用MATLAB求出离散序列求出离散序列的的Z变换变换MATLAB程序如下：程序如下：symskzf=0.5k;%定义离散信号定义离散信号Fz=ztrans(f)%对离散信号进行对离散信号进行Z变换变换u运行结果如下：运行结果如下：Fz=2*z/(2*z-1)13Z变换u例例已知一离散信号的已知一离散信号的Z变换式为变换式为,求出它所对应的离散信号求出它所对应的离散信号f(k).MATLAB程序如下：程序如下：symskzFz=2*z/(2*z-1);%定义定义Z变换表达式变换表达式fk=iztrans(Fz,k)%求反求反Z变换变换u运行结果如下：运行结果如下：fk=(1/2)ku例例：求序列的：求序列的Z变换变换.symsnhn=sym(kroneckerDelta(n,1)+kroneckerDelta(n,2)+kroneckerDelta(n,3)Hz=ztrans(hn)Hz=simplify(Hz)u运行结果如下：运行结果如下：Fz=(z2+z+1)/z314 1filter.m本文件用来求离散系统的输出y(n)。若系统的h(n)已知，由y(n)=x(n)*h(n)，用conv.m文件可求出y(n)。y=conv(x,h)filter文件是在A(z)、B(z)已知，但不知道h(n)的情况下求y(n)的。调用格式是:y=filter(b,a,x)x,y,a和b都是向量。与逆与逆Z变换变换相关的相关的matlab函数函数15与逆与逆Z变换变换相关的相关的matlab函数函数2.impz.m在A（z）、B（z）已知情况下,求系统的单位抽样响应h(n)。调用格式是：h = impz(b, a, N) 或h,t=impz(b,a,N) N是所需的的长度。前者绘图时n从1开始，而后者从0开始。16 3. residuez.m将X(z)的有理分式分解成简单有理分式的和，因此可用来求逆变换。调用格式：r,p,k= residuez(b,a)假如知道了向量r,p和k，利用residuez.m还可反过来求出多项式A(z)、B(z)。格式是b,a= residuez(r,p,k)。17频率响应函数：频率响应函数：freqz.m已知A(z)、B(z),求系统的频率响应。基本的调用格式是：H,w=freqz(b,a,N,whole,Fs)N是频率轴的分点数，建议N为2的整次幂；w是返回频率轴座标向量，绘图用；Fs是抽样频率，若Fs1，频率轴给出归一化频率；whole指定计算的频率范围是从0FS,缺省时是从0FS/2.幅频率响应：幅频率响应：Hr=abs(H);相频响应：相频响应： Hphase=angle(H);解卷绕：解卷绕：Hphase=unwrap(Hphase);18下面几个文件用于转移函数与极零点之间的相互转换及极零点的排序：（1)tf2zpk.m,调用z,p,k=tf2zpk(b,a)(2)zp2tf.m,调用b,a=zp2tf(z,p,k)(3)roots.m，调用Z1=roots(b)(4)poly.m,调用b=poly(Z1)与极零点有关的与极零点有关的MATLAB函数函数19显示离散系统的极零图显示离散系统的极零图函数函数：zplane.m本文件可用来显示离散系统的极零图。其调用格式是：zplane(z,p),或zplane(b,a),前者是在已知系统零点的列列向向量量z和极点的列列向向量量p的情况下画出极零图，后者是在仅已知A(z)、B(z)的情况下画出极零图。20用极零分析大致画出幅频响应及相频响应：例1： 系统解：转移函数：b=1 -4 4;a=1;z,p,k=tf2zpk(b,a)zplane(b,a)zplane(z,p)r,p,k=residuez(b,a)b,a=residuez(r,p,k)z=2; 2p=0; 0K=1r =; p=;k=1 -4 4;21用极零分析大致画出幅频响应及相频响应：例1： 系统解：转移函数：b=1 -4 4;a=1;H,w=freqz(b,a,128,whole,1)Hr=abs(H);Hphase=angle(H);Hphaseu=unwrap(Hphase);subplot(311),plot(Hr)subplot(312),plot(Hphase)subplot(313),plot(Hphaseu)22例例2：相位的卷绕相位的卷绕 (wrapping) 解卷绕解卷绕 b=1;a=0,1;H,w=freqz(b,a,256,whole,1);Hr=abs(H);Hphase=angle(H);Hphase1=unwrap(Hphase); 23例：例： 给定系统给定系统求：极零图单位抽样响应频率响应H,w=freqz(b,a,256,whole,1); Hr=abs(H); Hphase=angle(H); Hphase=unwrap(Hphase); h,t=impz(b,a,40);stem(t,h,.);grid on;zplane(b,a);b=.1836 .7344 1.1016 .7374 .1836/100a =1 -3.0544 3.8291 -2.2925 .55075解：24极零图极零图 zplane(b,a); 25单位抽样响应h,t=impz(b,a,40);stem(t,h,.);grid on;26频率响应Hphase=angle(H); Hphaseu=unwrap(Hphase); H,w=freqz(b,a,256,whole,1); Hr=abs(H); subplot(311),plot(Hr)subplot(312),plot(Hphase)subplot(313),plot(Hphaseu)27下面几个文件用于转移函数、极零点与二二阶阶子子系系统统sos(Second-OrderSection)级联之间的相互转换：（1)tf2sos.m,调用sos,G=tf2sos(b,a)(2)sos2tf.m,调用b,a=sos2tf(sos,G)(3)sos2zp.m,调用z,p,k=sos2zp(sos,G)(4)zp2sos.m,调用sos,G=zp2sos(z,p,k)与信号流图有关的与信号流图有关的MATLAB函数函数sos是一是一L6的矩阵，每行元素如下排列：的矩阵，每行元素如下排列：281buttord.m确定LPDF、或LPAF的阶次；（1）N,Wn=buttord(Wp,Ws,Rp,Rs)；（2）N, Wn = buttord(Wp, Ws, Rp, Rs,s)：与本章内容有关的MATLAB文件(1)对应数字滤波器。其中Wp,Ws分别是通带和阻带的截止频率，其值在01之间，1对应抽样频率的一半(归一化频率)。对低通和高通，Wp,Ws都是标量，对