神经网络第二章-山东大学课件

BiomedicalSignalprocessingmatlab信号处理函数ZhongguoLiuBiomedicalEngineeringSchoolofControlScienceandEngineering,ShandongUniversity2024/7/131ZhongguoLiu_BiomedicalEngineering_ShandongUniv.MATLAB是美国MathWorks公司开发的一种功能极其强大的高技术计算语言和内容极其丰富的软件库。它以矩阵和向量的运算以及运算结果的可视化为基础，把广泛应用于各个学科领域的数值分析、矩阵计算、函数生成、信号、图形及图象处理、建模与仿真等诸多强大功能集成在一个便于用户使用的交互式环境之中，为使用者提供了一个高效的编程工具及丰富的算法资源。关于MATLABMATLAB与信号处理直接有关的工具箱（Toolbox）SignalProcessing(信号处理工具箱）Wavelet(小波工具箱)ImageProcessing（图象处理工具箱）Higher-OrderSpectralAnalysis（高阶谱分析工具箱）与信号处理间接有关的工具箱：ControlSystem(控制系统)Communication（通信）SystemIdentification（系统辨识）Statistics（统计）NeuralNetwork(神经网络)例：例：z=peaks;surf(z)；1.rand.m用来产生均值为用来产生均值为0.5、幅度在、幅度在01之之 间间 均均 匀匀 分分 布布 的的 伪伪 白白 噪噪 声声:u=rand(N,1)(rand(N)生成生成N阶矩阵阶矩阵)方差：如何改变 的方差与第二章内容有关的MATLAB文件方差函数var(u)标准差函数std(u)1.rand.m用来产生均值为用来产生均值为0.5、幅度在、幅度在01之之 间间 均均 匀匀 分分 布布 的的 伪伪 白白 噪噪 声声:u=rand(N,1)(rand(N)生成生成N阶矩阵阶矩阵)2.randn.m用来产生均值为零、方差为用来产生均值为零、方差为13.服从高斯（正态）分布的白噪声信号服从高斯（正态）分布的白噪声信号u=randn(1,N)与第二章内容有关的MATLAB文件x=randn(1000,1)y=randn(1000,1)v=var(x)h=std(y)3.sinc：用来：用来产产生生“sinc”函数：函数：sinc函数定函数定义为义为：t=-4:0.1:4;x4=sinc(t);%产生抽样函数plot(t,x4)4.conv.m用来实现两个离散序列的线性卷用来实现两个离散序列的线性卷积。其调用格式是：积。其调用格式是：y=conv(x,h).若若x(n)和和y(n)的长度分别为的长度分别为M和和N，则返回值则返回值是长度为是长度为M+N-1的序列。的序列。例例 x(n)=3 4 5;h(n)=2 6 7 8，求求其其线线性性卷卷积。积。MATLAB语句如下：语句如下：x=3 4 5;h=2 6 7 8;y=conv(x,h)运行结果运行结果:y=6 26 55 82 67 40 u两序列的相关运算两序列的相关运算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 -05xcorr:其其互互相相关关和和自自相相关关。格格式式是是：(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同样适用于求互相关。同样适用于求互相关。第三章 Z变换.在在MATLAB语语言中有言中有专门对专门对信号信号进进行正反行正反Z变换变换的函数的函数ztrans()和和itrans()。其。其调调用格式分用格式分别别如下：如下：F=ztrans(f)对对f(n)进进行行Z变换变换，其，其结结果果为为F(z)F=ztrans(f,v)对对f(n)进进行行Z变换变换，其，其结结果果为为F(v)F=ztrans(f,u,v)对对f(u)进进行行Z变换变换,其其结结果果为为F(v)f=itrans(F)对对F(z)进进行行Z反反变换变换,其其结结果果为为f(n)f=itrans(F,u)对对F(z)进进行行Z反反变换变换，其，其结结果果为为f(u)f=itrans(F,v,u)对对F(v)进进行行Z反反变换变换,其其结结果果为为f(u)注意：注意：在在调调用函数用函数ztran()及及iztran()之前，要用之前，要用syms命令命令对对所有所有需要用到的需要用到的变变量（如量（如t,u,v,w）等）等进进行行说说明，即要将明，即要将这这些些变变量量说说明成符号明成符号变变量量Z变换例例.求数列求数列fn=e-n的的Z变换及其逆及其逆变换。命令如下：。命令如下：symsnzfn=exp(-n);Fz=ztrans(fn,n,z)%求求fn的的Z变换f=iztrans(Fz,z,n)%求求Fz的逆的逆Z变换u例例用用MATLAB求出离散序列求出离散序列的的Z变换MATLAB程序如下：程序如下：usymskzuf=0.5k;%定定义离散信号离散信号uFz=ztrans(f)%对离散信号离散信号进行行Z变换u运行运行结果如下：果如下：uFz=2*z/(2*z-1)Z变换u例例已知一离散信号的已知一离散信号的Z变换式式为,求出它所求出它所对应的离散信号的离散信号f(k).MATLAB程序如下：程序如下：usymskzuFz=2*z/(2*z-1);%定定义Z变换表达式表达式ufk=iztrans(Fz,k)%求反求反Z变换u运行运行结果如下：果如下：fk=(1/2)ku例例：求序列的：求序列的Z变换.usymsnuhn=sym(kroneckerDelta(n,1)+kroneckerDelta(n,2)+kroneckerDelta(n,3)uHz=ztrans(hn)uHz=simplify(Hz)u运行运行结果如下：果如下：Fz=(z2+z+1)/z31filter.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函数函数与逆与逆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开始。开始。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)。频率响应函数：频率响应函数：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);下面几个文件用于转移函数与极零点之间的相互转换及极零点的排序：（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函数函数显示离散系统的极零图函数：zplane.m本文件可用来显示离散系统的极零图。其调用格式是：zplane(z,p),或zplane(b,a),前者是在已知系统零点的列向量z和极点的列向量p的情况下画出极零图，后者是在仅已知A(z)、B(z)的情况下画出极零图。用极零分析大致画出幅频响应及相频响应：例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;用极零分析大致画出幅频响应及相频响应：例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)例例2：相位的卷绕相位的卷绕(wrapping)解卷绕解卷绕 b=1;a=0,1;H,w=freqz(b,a,256,whole,1);Hr=abs(H);Hphase=angle(H);Hphase1=unwrap(Hphase);例：例：给定系统给定系统求：极零图单位抽样响应频率响应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解：极零图极零图 zplane(b,a);单位抽样响应h,t=impz(b,a,40);stem(t,h,.);grid on;频率响应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)下面几个文件用于转移函数、极零点与二阶子系统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的矩阵，每行元素如下排列：的矩阵，每行元素如下排列：1buttord.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都是标量，对带通和带阻，Wp,Ws是12的向量。Rp,Rs分别是通带和阻带的衰减(dB)。N是求出的相应低通滤波器的阶次，Wn是求出的3dB频率，它和Wp稍有不同。(2)对应模拟滤波器，各变量含意和（1）相同，但Wp,Ws及Wn的单位为弧度/秒，它们实际上是频率。2buttap.m设计模拟低通(Butt)原型滤波器。z,p,k=buttap(N):N是欲设计的低通原型滤波器的阶次，z,p,k是设计出的极点、零点及增益。3lp2lp.m、lp2hp.m、lp2bp.m,lp2bs.m将模拟低通原型转换为实际的低通、高通、带通及带阻滤波器。b,a是AFLP的分子、分母的系数向量，B,A是转换后的的分子、分母的系数向量；(1)中，Wo是低通或高通滤波器的截止频率；B,A=lp2lp(b,a,Wo)B,A=lp2hp(b,a,Wo)(1)B,A=lp2bp(b,a,Wo,Bw)B,A=lp2bs(b,a,Wo,Bw)(2)(2)中Wo是带通或带阻滤波器中心频率，Bw是其带宽。4bilinear.m：双线性变换，由模拟滤波器得到数字滤波器。Bz,Az=bilinear(B,A,Fs)式中B,A分别是G（s）的分子、分母多项式的系数向量，Bz,Az分别是H(z)的分子、分母多项式的系数向量，Fs是抽样频率。5butter.m用来直接设计Butterworth数字滤波器，实际上它把buttord.m,buttap.m,lp2lp.m,bilinear.m等文件都包含了进去，从而使设计过程更简捷。格式(1)(3)用来设计数字滤波器，B，A分别是H(z)的分子、分母多项式的系数向量，Wn是通带截止频率，范围在01之间。若Wn是标量，(1)用来设计低通数字滤波器，若Wn是12的向量，则(1)用来设计数字带通滤波器；(2)用来设计数字高通滤波器；(3)用来设计数字带阻滤波器，显然，这时的Wn是12的向量；格式(4)用来设计模拟滤波器。(1)B,A=butter(N,Wn);(2)B,A=butter(N,Wn,high);(2)(3)B,A=butter(N,Wn,stop);(4)B,A=butter(N,Wn,s)例6.7.1(例6.5.1)uclearall;fp=100;fs=300;Fs=1000;rp=3;rs=20;uwp=2*pi*fp/Fs;ws=2*pi*fs/Fs;uFs=Fs/Fs;%letFs=1u%Firstlytofinishfrequencyprewarping;uwap=tan(wp/2);was=tan(ws/2);%un,wn=buttord(wap,was,rp,rs,s)%Note:s!uz,p,k=buttap(n);%ubp,ap=zp2tf(z,p,k)%ubs,as=lp2lp(bp,ap,wap)%u%Note:s=(2/Ts)(z-1)/(z+1);Ts=1,thatis2Fs=1,Fs=0.5;ubz,az=bilinear(bs,as,Fs/2)%uh,w=freqz(bz,az,256,Fs*1000);uplot(w,abs(h);gridon;设计IIRLPDF，例6.7.1(例6.5.1)uclearall;uwp=.2*pi;ws=.6*pi;Fs=1000;rp=3;rs=20;%u%Firstlytofinishfrequencyprewarping;uwap=2*Fs*tan(wp/2);was=2*Fs*tan(ws/2);un,wn=buttord(wap,was,rp,rs,s);%Note:s!uz,p,k=buttap(n);ubp,ap=zp2tf(z,p,k);ubs,as=lp2lp(bp,ap,wap)uw1=0:499*2*pi;uh1=freqs(bs,as,w1);ubz,az=bilinear(bs,as,Fs)%Note:z=(2/ts)(z-1)/(z+1);uh2,w2=freqz(bz,az,500,Fs);uplot(w1/2/pi,abs(h1),w2,abs(h2),k);gridon;设计IIRLPDF，例6.7.1(例6.5.1)uclearall;uwp=.2*pi;ws=.6*pi;Fs=1000;urp=3;rs=20;un,wn=buttord(wp/pi,ws/pi,rp,rs);ubz,az=butter(n,wp/pi)ubz1,az1=butter(n,wn)uh,w=freqz(bz,az,128,Fs);uh1,w1=freqz(bz1,az1,128,Fs);uplot(w,abs(h),w1,abs(h1),g.);gridon;设计IIRLPDF，非线性关系设计的AF并不是按给定的技术指标，但当再由变回后，保证了DF的技术要求。又称为频率的预变形（Freq.Warping)。例如：抽样频率给出给出数字高通数字高通的技术要求的技术要求得到得到模拟高通模拟高通的技术要求的技术要求得到得到模拟低通模拟低通的技术要求的技术要求设计出设计出得到得到模拟高通转移模拟高通转移函数函数最后得到最后得到数字高通转移数字高通转移函数函数数字高通滤波器设计步骤7.6 数字高通,带通及带阻滤波器的设计对带通（BP）、带阻（BS）数字滤波器的设计，只需改变图中Step2和Step4：带阻带通要求：按上述转换办法，可以求出：例6.6.2：设计一IIRBPDF,要求：通带频率范围：300Hz400Hz;阻带频率范围：200Hz、500Hz1buttord.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都是标量，对带通和带阻，Wp,Ws是12的向量。Rp,Rs分别是通带和阻带的衰减(dB)。N是求出的相应低通滤波器的阶次，Wn是求出的3dB频率，它和Wp稍有不同。(2)对应模拟滤波器，各变量含意和（1）相同，但Wp,Ws及Wn的单位为弧度/秒，它们实际上是频率。2buttap.m设计模拟低通(Butt)原型滤波器。z,p,k=buttap(N):N是欲设计的低通原型滤波器的阶次，z,p,k是设计出的极点、零点及增益。3lp2lp.m、lp2hp.m、lp2bp.m,lp2bs.m将模拟低通原型转换为实际的低通、高通、带通及带阻滤波器。b,a是AFLP的分子、分母的系数向量，B,A是转换后的的分子、分母的系数向量；(1)中，Wo是低通或高通滤波器的截止频率；B,A=lp2lp(b,a,Wo)B,A=lp2hp(b,a,Wo)(1)B,A=lp2bp(b,a,Wo,Bw)B,A=lp2bs(b,a,Wo,Bw)(2)(2)中Wo是带通或带阻滤波器中心频率，Bw是其带宽。4bilinear.m：双线性变换，由模拟滤波器得到数字滤波器。Bz,Az=bilinear(B,A,Fs)式中B,A分别是G（s）的分子、分母多项式的系数向量，Bz,Az分别是H(z)的分子、分母多项式的系数向量，Fs是抽样频率。5butter.m用来直接设计Butterworth数字滤波器，实际上它把buttord.m,buttap.m,lp2lp.m,bilinear.m等文件都包含了进去，从而使设计过程更简捷。格式(1)(3)用来设计数字滤波器，B，A分别是H(z)的分子、分母多项式的系数向量，Wn是通带截止频率，范围在01之间。若Wn是标量，(1)用来设计低通数字滤波器，若Wn是12的向量，则(1)用来设计数字带通滤波器；(2)用来设计数字高通滤波器；(3)用来设计数字带阻滤波器，显然，这时的Wn是12的向量；格式(4)用来设计模拟滤波器。(1)B,A=butter(N,Wn);(2)B,A=butter(N,Wn,high);(2)(3)B,A=butter(N,Wn,stop);(4)B,A=butter(N,Wn,s)6cheb1ord.m求Cheb-型滤波器的阶；7cheb1ap.m设计原型低Cheb-I型模拟滤波器；8cheby1.m直接设计数字Cheb-滤波器。以上三个文件的调用格式和对应的Butterworth滤波器的文件类似。9cheb2ord.m;10.ellipord.m;11.cheb2ap.m;12.ellipap.m;13.besselap.m;14.cheby2.m;15.ellip.m;16.besself.m17impinvar.m用冲激响应不变法实现频率转换；对应Cheby-II、椭圆IIR滤波器产生窗函数的文件有八个：bartlett（三角窗）;2.blackman（布莱克曼窗）;3.boxcar（矩形窗）;4.hamming（哈明窗）;5.hanning（汉宁窗）;6.triang（三角窗）;7.chebwin（切比雪夫窗）;8.kaiser（凯赛窗）;两端为零两端不为零调用方式都非常简单请见help文件稍为复杂9fir1.m用“窗函数法”设计FIRDF。调用格式：（1）b=fir1(N,Wn);(2)b=fir1(N,Wn,high);(3)b=fir1(N,Wn,stop);N：阶次，滤波器长度为N1；Wn：通带截止频率，其值在01之间，1对应Fs/2；b:滤波器系数。格式（2）用来设计高通滤波器，格式（3）用来设计带阻滤波器。格式（1），若Wn为标量，则设计低通滤波器，若Wn是12的向量，则用来设计带通滤波器，若Wn是1L的向量，则可用来设计L带滤波器。对格式（1），若Wn为标量，则设计低通滤波器，若Wn是12的向量，则用来设计带通滤波器，若Wn是1L的向量，则可用来设计L带滤波器。这时，格式（1）要改为:b=fir1(N,Wn,DC-1),或b=fir1(N,Wn,DC-0)前者保证第一个带为通带，后者保证第一个带为阻带。在上述所有格式中，若不指定窗函数的类型，fir1自动选择Hamming窗。指定窗函数格式：（4）b=fir1(N,Wn,wind);例b=fir1(N,Wn,boxcar(N+1);指定矩形窗10fir2.m本文件采用“窗函数法”设计具有任意幅频相应的FIR数字滤波器。其调用格式是：b=fir2(N,F,M);F是频率向量，其值在01之间，M是和F相对应的所希望的幅频相应。如同fir1,缺省时自动选用Hamming窗。例：设计一多带滤波器，要求频率在0.20.3,0.60.8之间为1，其余处为零。设计结果如下：N=30,90时幅频响应响应及理想幅频响应；N=30N=9013firls.m用最小平方法设计线性相位FIR滤波器，可设计任意给定的理想幅频响应；14fircls.m用带约束的最小平方法设计线性相位FIR滤波器，可设计任意给定的理想幅频响应；15fircls1.m用带约束的最小平方方法设计线性相位FIR低通和高通滤波器。16sgolay.m用来设计Savitzky-GolayFIR平滑滤波器，其原理见9.1.1节17firrcos.m用来设计低通线性相位FIR滤波器，其过渡带为余弦函数形状。9fir1.m用“窗函数法”设计FIRDF。调用格式：（1）b=fir1(N,Wn);(2)b=fir1(N,Wn,high);(3)b=fir1(N,Wn,stop);N：阶次，滤波器长度为N1；Wn：通带截止频率，其值在01之间，1对应Fs/2；b:滤波器系数。格式（2）用来设计高通滤波器，格式（3）用来设计带阻滤波器。格式（1），若Wn为标量，则设计低通滤波器，若Wn是12的向量，则用来设计带通滤波器，若Wn是1L的向量，则可用来设计L带滤波器。对格式（1），若Wn为标量，则设计低通滤波器，若Wn是12的向量，则用来设计带通滤波器，若Wn是1L的向量，则可用来设计L带滤波器。这时，格式（1）要改为:b=fir1(N,Wn,DC-1),或b=fir1(N,Wn,DC-0)前者保证第一个带为通带，后者保证第一个带为阻带。在上述所有格式中，若不指定窗函数的类型，fir1自动选择Hamming窗。指定窗函数格式：（4）b=fir1(N,Wn,wind);例b=fir1(N,Wn,boxcar(N+1);指定矩形窗例7.1.1.设计低通DFFIR,令截止频率0.25,取M10,20，40，观察其幅频响应的特点.clear all;N=10;b1=fir1(N,0.25,boxcar(N+1);b3=fir1(2*N,0.25,boxcar(2*N+1);b4=fir1(4*N,0.25,boxcar(4*N+1);M=128;h1=freqz(b1,1,M);h3=freqz(b3,1,M);h4=freqz(b4,1,M);f=0:0.5/M:0.5-0.5/M;plot(f,abs(h1),f,abs(h3),f,abs(h4);grid;axis(0 0.5 0 1.2)例7.1.2：理想差分器及其设计clear all;N=40;n=0:N;b1=fir1(N,0.25,boxcar(N+1);b2=fir1(N,0.25,hamming(N+1);win=hamming(N+1);for n=1:N+1 if(n-1-N/2)=0;b1(n)=0;else b1(n)=(-1)(n-1-N/2)/(n-1-N/2);end endfor n=1:N+1 if(n-1-N/2)=0;b2(n)=0;elseb2(n)=win(n)*(-1)(n-1-N/2)/(n-1-N/2);end endM=128;h1=freqz(b1,1,M);h2=freqz(b2,1,M);%h=freqz(b,1,M);f=0:0.5/M:0.5-0.5/M;hd=2*pi*f;plot(f,abs(h1),f,abs(h2),f,hd,k-)11.remez.m设计Chebyshev最佳一致逼近FIR滤波器、Hilbert变换器和差分器。调用格式是：(1)b=remez(N,F,A);(2)b=remez(N,F,A,W);（3）b=remez(N,F,A,W,Hilbert);(4)b=remez(N,F,A,W,differentiator)N是给定的滤波器的阶次，b是设计的滤波器的系数，其长度为N1；F是频率向量，A是对应F的各频段上的理想幅频响应，W是各频段上的加权向量。F、A及W的指定方式和例7.4.1和7.4.2所讨论过的一样，唯一的差别是F的范围为01，而非00.5,1对应抽样频率的一半。需要指出的是，若b的长度为偶数，设计高通和带阻滤波器时有可能出现错误，因此，最好保证b的长度为奇数，也即N应为偶数。例7.4.1：设计低通FIRDF:b=remez(N,F,A,W)clear all;f=0.6.7 1;%给定频率轴分点；A=1 1 0 0;%频率分点上理想幅频响应；weigh=1 10;%频率分点上的加权；b=remez(32,f,A,weigh);%设计出切比雪夫最佳一致逼近滤波器;h,w=freqz(b,1,256,1);h=abs(h);h=20*log10(h);figure(1);stem(b,.);grid;figure(2);plot(w,h);grid;调整通带、阻带的加权及滤波器的长度。调整N或W的结果例7.4.2：设计多阻带滤波器，抽样频率500Hz,在50Hz、100Hz及150Hz处陷波。通带加权为8，阻带为1-17dB通带、阻带加权都是1-25dB例7.4.2：设计多阻带滤波器，抽样频率500Hz,在50Hz、100Hz及150Hz处陷波。clear all;f=0.14.18.22.26.34.38.42.46.54.58.62.66 1;A=1 1 0 0 1 1 0 0 1 1 0 0 1 1;weigh1=8 1 8 1 8 1 8;b1=remez(64,f,A,weigh1);h1,w1=freqz(b1,1,256,1);h1=abs(h1);h1=20*log10(h1);subplot(211);plot(w1,h1);grid;axis(0 0.5-60 10)title(N=32,weight=8 1 8 1 8 1 8,FontSize,14,Color,r)250Hz12remezord.m本文件用来确定在用Chebyshev最佳一致逼近设计FIR滤波器时所需要的滤波器阶次。其调用格式是：N,Fo,Ao,W=remezord(F,A,DEV,Fs)。F、A的含意同文件remez，DEV是通带和阻带上的偏差；输出的是适合要求的滤波器阶次N、频率向量Fo、幅度向量Ao和加权向量W。若设计者事先不能确定要设计的滤波器的阶次，那么，调用remezord后，就可利用这一族参数调用remez,即b=remez(N,Fo,Ao,W)，从而设计出所需要滤波器。因此，remez和remezord常结合起来使用。需要说明的是，remezord给出的阶次N有可能偏低，这时适当增加N即可；另外，最好判断一下，若N为奇数，就令其加一，使其变为偶数，这样b的长度为奇数。fftfilt.m用叠接相加法实现长序列卷积。格式是:y=fftfilt(h,x)或y=fftfilt(h,x，N)与本章有关的 MATLAB 文件记的长度为，的长度为。若采用第一个调用方式，程序自动地确定对分段的长度及做FFT的长度，显然，是最接近的2的整次幂。分的段数为。采用第二个调用方式，使用者可自己指定做FFT的长度。建议使用第一个调用方式。clear;%用叠接相加法，计算滤波器系数用叠接相加法，计算滤波器系数h和输入信号和输入信号x的卷积的卷积%其中其中h为为10阶阶hanning窗，窗，x是带有高斯白噪的正弦信号是带有高斯白噪的正弦信号h=fir1(10,0.3,hanning(11);%h:is the impulse responseN=500;p=0.05;f=1/16;%of a low-pass filter.u=randn(1,N)*sqrt(p);%u:white noises=sin(2*pi*f*0:N-1);%s:sine signalx=u(1:N)+s;%x:a long sequence;y=fftfilt(h,x);%y=x*hsubplot(211)plot(x);subplot(212)plot(y);例3.9.1令x(n)为一正弦加白噪声信号，长度为500，h(n)是用fir1.m文件设计出的一个低通FIR滤波器，长度为11.试用fftfilt实现长序列的卷积clear;%生成滤波器系数h和混有高斯白噪的正弦信号xh=fir1(10,0.3,hanning(11);N=500;p=0.05;f=1/16;u=randn(1,N)*sqrt(p);%s=sin(2*pi*f*0:N-1);x=u(1:N)+s;%将x分为长度为L的小段L=20;M=length(h);y=zeros(1,N+M-1);tempy=zeros(1,M+L-1);tempX=zeros(1,L);for k=0:N/L-1 tempx(1:L)=x(k*L+1:(k+1)*L);tempy=conv(tempx,h);y=y+zeros(1,k*L),tempy,zeros(1,N-(k+1)*L);endsubplot(211);plot(x)subplot(212);plot(y(1:N)hilbert.m 文件用来文件用来计计算信号算信号Hilbert变换变换。调调用用的格式是的格式是:y=hilbert(x),y的的实实部就部就是是 ，虚部是的，虚部是的Hilbert变换变换 。所以，所以，y 实际实际上是上是 x 的解析信号。的解析信号。czt.m调用格式是：Xczt(x,M,W,A)。x是待变换的时域信号，其长度设为N，M是变换的长度，W确定变换的步长，A确定变换的起点。若M=N,A=1,则CZT变成DFT。A=exp(j*2*pi*f0/fs);W=exp(-j*2*pi*DELf/fs);例例4.7.2设设x(n)是是两两个个正正弦弦信信号号及及白白噪噪声声的的叠叠加加，进进行行频频谱分析。谱分析。clearall;%产生两个正弦加白噪声；产生两个正弦加白噪声；N=256;f1=.1;f2=.2;fs=1;a1=5;a2=3;w=2*pi/fs;x=a1*sin(w*f1*(0:N-1)+a2*sin(w*f2*(0:N-1)+randn(1,N);%应用应用FFT求频谱；求频谱；subplot(3,1,1);plot(x(1:N/4);f=-0.5:1/N:0.5-1/N;X=fft(x);y=ifft(X);subplot(3,1,2);plot(f,fftshift(abs(X);subplot(3,1,3);plot(real(y(1:N/4);例例4.7.2设x(n)由三个由三个实正弦正弦组成，成，频率分率分别是是8HZ，8.22HZ和和9HZ，抽抽样频率是率是40HZ，时域取域取128点。点。用用用用CZTCZT计算计算计算计算的的的的DFTDFT用用用用FFTFFT计算计算计算计算的的的的DFTDFT7(7+M0.05)7(7+M0.05)HZHZ的的的的CZTCZT变换变换变换变换 程序clear all；%构造三个不同频率的正弦信号的叠加作为试验信号N=128;f1=8;f2=8.22;f3=9;fs=40;stepf=fs/N;n=0:N-1;t=2*pi*n/fs;n1=0:stepf:fs/2-stepf;x=sin(f1*t)+sin(f2*t)+sin(f3*t);M=N;W=exp(-j*2*pi/M);%A=1时的czt变换A=1;Y1=czt(x,M,W,A);subplot(311)plot(n1,abs(Y1(1:N/2);grid on;%DTFTY2=abs(fft(x);subplot(312)plot(n1,abs(Y2(1:N/2);grid on;%详细构造A后的czt M=60;f0=7.2;DELf=0.05;A=exp(j*2*pi*f0/fs);W=exp(-j*2*pi*DELf/fs);Y3=czt(x,M,W,A);n2=f0:DELf:f0+(M-1)*DELf;subplot(313);plot(n2,abs(Y3);grid on;1filtfilt.m本文件实现零相位滤波。其调用格式是：y=filtfilt(B,A,x)。式中B是的分子多项式，A是分母多项式，x是待滤波信号，y是滤波后的信号。clear;N=32;n=-N/2:N+N/2;w=0.1*pi;x=cos(w*n)+cos(2*w*n);subplot(311);stem(n,x,.);gridon;xlabel(n);b=0.067450.13490.06745;a=1-1.1430.4128;y=filtfilt(b,a,x);%用给定系统(b,a)对信号x作零相位滤波；y1=filter(b,a,x);%用给定系统(b,a)对信号x作低通滤波；subplot(312);stem(n,y,.);gridon;xlabel(n);subplot(313);stem(n,y1,.);gridon;xlabel(n);与本章内容有关的MATLAB文件2grpdelay.m求系统的群延迟。调用格式gdw=grpdelay(B,A,N),或gdF=grpdelay(B,A,N,FS)式中B和A仍是的分子、分母多项式，gd是群延迟，w、F是频率分点，二者的维数均为N；FS为抽样频率，单位为Hz。3deconv.m：实现系统的反卷积，其调用格式：q,r=deconv(y,x);也用来实现多项式除法。clearall;k=0:1:7;x=k+1;h=ones(1,4);y=conv(h,x);%y=x*h;q,r=deconv(y,x);%由y,x作反卷积，求出h；q1,r1=deconv(y,h);%由y,h作反卷积，求出x；subplot(321);stem(h,.b);ylabel(h(n);subplot(322);stem(x,.b);ylabel(x(n);subplot(323);stem(y,.b);ylabel(y(n);subplot(325);stem(q,.r);ylabel(q(n);subplot(326);stem(q1,.r);ylabel(q1(n);clearall;%实现多项式除法q=(x3+1)/(x+1)y=1001;x=11;q,r=deconv(y,x);q3tf2latc.m和latc2tf.m：实现转移函数和Lattice系数之间的相互转换。tf2latc的调用格式是：(1)k=tf2latc(b),(2)k=tf2latc(1,a),(3)k,c=tf2latc(b,a),其中（1）对应全零系统，（2）对应全极系统，（3）对应极零系统。latc2tf的调用格式和tf2latc正好相反。需要说明的是，tf2latc求出的Lattice系数k和本书求出的k差一个负号，这是由于我们在图中用的是k。4.latcfilt.m用来实现Lattice结构下的信号滤波。调用格式是：(1)y,g=latcfilt(k,x)：对应全零系统(2)y,g=latcfilt(k,1,x)：对应全极系统(3)y,g=latcfilt(k,c,x)：对应极零系统x是待滤波的信号，y是用Lattice结构作正向滤波的输出，g是作反向滤波的输出。若输入x是则输出y是的系数；g是的系数。5.tf2ss.m和ss2tf.m实现转移函数和相应状态变量之间的转换。二者的调用格式分别是：A,B,C,D=tf2ss(b,a),b,a=ss2tf(A,B,C,D)。式中b,a分别是分子、分母多项式的系数向量，A,B,C及D的定义见书。5.sos2ss.m实现由转移函数的二阶级联形式转换为状态变量表示。调用格式：A,B,C,D=sos2ss(sos,g),A,B,C,D的定义 见 书。有 关 sos和 g的 说 明 见 2.8节 关 于tf2sos.m的说明。