Matlab中的信号处理函数

1、12021-11-21Zhongguo Liu_Biomedical Engineering_Shandong Univ.Biomedical Signal processingmatlab 信号处理函数信号处理函数Zhongguo LiuBiomedical EngineeringSchool of Control Science and Engineering, Shandong University2MATLAB是美国MathWorks公司开发的一种功能极其强大的高技术计算语言和内容极其丰富的软件库。它以矩阵和向量的运算以及运算结果的可视化为基础，把广泛应用于各个学科领域的数值分析、矩阵

2、计算、函数生成、信号、图形及图象处理、建模与仿真等诸多强大功能集成在一个便于用户使用的交互式环境之中，为使用者提供了一个高效的编程工具及丰富的算法资源。 关于MATLAB3 MATLAB与信号处理直接有关的工具箱与信号处理直接有关的工具箱 （ Toolbox） Signal Processing (信号处理工具箱） Wavelet (小波工具箱) Image Processing（图象处理工具箱） Higher-Order Spectral Analysis （高阶谱分析工具箱）4 与信号处理间接有关的工具箱：与信号处理间接有关的工具箱：Control System (控制系统)Communi

3、cation （通信）System Identification （系统辨识）Statistics （统计）Neural Network (神经网络)5例例： z=peaks; surf(z)；61. rand.m 用来产生均值为0.5、幅度在 01之间均匀分布的伪白噪声: u=rand(N,1) (rand(N)生成N阶矩阵) 2122011( ),12Nuuunu nN方差：( )u n与第二章内容有关的MATLAB文件如何改变 的方差为P( )u n方差函数var(u)标准差函数std(u)即如何确定即如何确定a使使au的方差为P？将将au代替u带入上面方差公式可得22uaP71. ra

4、nd.m 用来产生均值为0.5、幅度在 01之间均匀分布的伪白噪声: u=rand(N,1) (rand(N)生成N阶矩阵) 2. randn.m 用来产生均值为零、方差为1 服从高斯（正态）分布的白噪声信号 u=randn(1, N) ( )u n与第二章内容有关的MATLAB文件x=randn(1000,1)y=randn(1000,1)v=var(x)h=std(y)83.sinc ：用来产生 “sinc” 函数：sinc函数定义为： 10sin ( )0sin()tc ttkttt为其它t=-4:0.1:4;x4=sinc(t); %产生抽样函数plot(t,x4) 94. conv.

5、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两序列的相关运算两序列的相关运算12()( )()ny mxn xnmMATLAB实现：实现：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:

6、其互相关和自相关。格式是：(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( )。其调用格式分别如

7、下：。其调用格式分别如下：uF=ztrans( f ) 对对f(n)进行进行Z变换，其结果为变换，其结果为F(z)uF=ztrans(f,w) 对对f(n)进行进行Z变换，其结果为变换，其结果为F(w)uF=ztrans(f,k,w) 对对f(k)进行进行Z变换变换,其结果为其结果为F(w)uf=itrans ( F ) 对对F(z)进行进行Z反变换反变换,其结果为其结果为f(n)uf=itrans(F,k) 对对F(z)进行进行Z反变换，其结果为反变换，其结果为f(k)uf=itrans(F,w,k) 对对F(w)进行进行Z反变换反变换,其结果为其结果为f(k)u注意：注意： 在调用函数在调

8、用函数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求出离散序列求出离散序列f=0.5k的的Z变换变换MATLAB程序如下：程序如下：sy

9、ms k zf=0.5k; %定义离散信号定义离散信号Fz=ztrans(f) %对离散信号进行对离散信号进行Z变换变换u运行结果如下：运行结果如下： Fz = 2*z/(2*z-1)Fz= z/(z - 1/exp(1)f = (1/exp(1)nsyms k hk=sym(kroneckerDelta(k, 1) + kroneckerDelta(k, 2)+ kroneckerDelta(k, 3)Hz=ztrans(hk)Hz=simplify(Hz)u运行结果如下：运行结果如下：Fz= (z2 + z + 1)/z313Z变换u例例 已知一离散信号的已知一离散信号的Z变换式为变换式为

10、 , 求出它所对应的离散信号求出它所对应的离散信号f(k). MATLAB程序如下：程序如下：syms k zFz=2*z/(2*z-1); %定义定义Z变换表达式变换表达式fk=iztrans(Fz,k) %求反求反Z变换变换u运行结果如下：运行结果如下：fk = (1/2)ku例例: 求序列的求序列的Z变换变换.( )(1)(4)f kkk 2( )21zF zz(1)(2)(3)kkk(1)(4)u ku k 阶跃序列符号阶跃序列符号14符号变换符号变换 Fourier变换及其反变换变换及其反变换uFw=fourier(ft,t,w) 求求“时域时域”函数函数ft的的Fourier变换变

11、换uft=ifourier(Fw,w,t) 求求“频域频域”函数函数Fw的的Fourier 反变换反变换 Laplace变换及其反变换变换及其反变换uFs=laplace(ft,t,s) 求求“时域时域”函数函数ft的的Laplace变换变换uft=ilaplace(Fs,s,t) 求求“频域频域”函数函数Fs的的Laplace 反变换反变换15【例例 】求求 的的Fourier变换。变换。0001)(tttf（1）求）求Fourier变换变换syms t wut=heaviside(t);UT=fourier(ut) UT = pi*dirac(w)-i/w （2）求）求Fourier反变换

12、反变换Ut=ifourier(UT,w,t)Ut =heaviside(t)16【例例2.5-4】求求 的的Laplace变换。变换。 syms t s; syms a b positive; %不限定则无结果 Dt=dirac(t-a); Ut=heaviside(t-b); Mt=Dt,Ut;exp(-a*t)*sin(b*t),t2*exp(-t); MS=laplace(Mt,t,s) MS = exp(-s*a), exp(-s*b)/s 1/b/(s+a)2/b2+1), 2/(s+1)3 tatetbtebtuat2sin)()(17 1filter.m本文件用来求离散系统的输出

13、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 函数函数18 1filter.mx=1,2,3,4; y=3,4,6z_conv= conv(x,y) % x , y 为输入和单位脉冲响应时输出为输入和单位脉冲响应时输出z_conv_= conv(y, x) % x , y为输入和单位

14、脉冲响应时输出为输入和单位脉冲响应时输出z_filter=filter(y,1,x) % x为输入为输入, y为为FIR单位脉冲响应时输出单位脉冲响应时输出z_filter_=filter(x,1,y) % y为输入为输入, x为为FIR单位脉冲响应时输出单位脉冲响应时输出与逆与逆Z变换变换 相关的相关的matlab 函数函数可见，conv(x,y)总是等于conv(y,x)。而filter(x,1,y)却不一定等于filter(y,1,x)，但是它们总是conv(x,y)截短的结果，截短的长度等于length(filter的第三个参数) z_conv = 3 10 23 36 34 24 z

15、_conv_=3 10 23 36 34 24z_filter = 3 10 23 36z_filter_= 3 10 2319与逆与逆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开始。20 3. residuez.m将X(z) 的有理分式分解成简单有理分式的和，因此可用来求逆Z变换。调用格式： r,p,k= residuez(b,a)假如知道了向量r, p和k，利用resi

16、duez.m还可反过来求出多项式A(z)、B(z)。格式是 b,a= residuez(r,p,k)。 214.频率响应函数：频率响应函数：freqz.m已知A(z)、B(z), 求系统的频率响应。基本的调用格式是： H,w=freqz(b,a,N,whole,Fs)N是频率轴的分点数，建议N为2的整次幂；w是返回频率轴座标向量，绘图用；Fs是抽样频率，若Fs1，频率轴给出归一化频率；whole指定计算的频率范围是从0FS,缺省时是从0FS/2.( )( )( )B zH zA z幅频响应：幅频响应：Hr=abs(H);相频响应：相频响应： Hphase=angle(H);解卷绕：解卷绕：Hp

17、hase=unwrap(Hphase);225filtfilt.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,.);grid on; xlabel(n);b=0.06745 0.1349 0.06745; a=1 -1.143 0.4128;y=filtfilt(b,a,x); % 用给定系统(b,a)对信号 x 作零相位滤波；

18、y1=filter(b,a,x); % 用给定系统(b,a)对信号 x 作低通滤波；subplot(312);stem(n,y,.);grid on; xlabel(n);subplot(313);stem(n,y1,.);grid on; xlabel(n);( )H z235filtfilt.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);st

19、em(n,x,.);grid on; xlabel(n);b=0.06745 0.1349 0.06745; a=1 -1.143 0.4128;y=filtfilt(b,a,x); % 用给定系统(b,a)对信号 x 作零相位滤波；y1=filter(b,a,x); % 用给定系统(b,a)对信号 x 作低通滤波；subplot(312);stem(n,y,.);grid on; xlabel(n);subplot(313);stem(n,y1,.);grid on; xlabel(n);( )H z-20-1001020304050-202nOriginal signal-20-10010

20、20304050-202nZero-phase Filtering-20-1001020304050-202nConventional Filtering24 6grpdelay.m 求系统的群延迟。调用格式 gd w=grpdelay(B, A, N) , 或 gd F=grpdelay(B, A, N, FS)式中B和A仍是 的分子、分母多项式，gd是群延迟，w、F是频率分点，二者的维数均为N；FS为抽样频率，单位为Hz。( )H z00.511.522.533.500.511.522.53 Amplitude Freq. Res.257deconv.m ：实现系统的反卷积，其调用格式：

21、q,r=deconv(y,x); 也用来实现多项式除法。clear all;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);

22、ylabel( q(n);subplot(326);stem(q1,.r);ylabel( q1(n);clear all;% 实现多项式除法q=(4*x3+1)/(2*x+1)y=4 0 0 1; x=2 1; q,r=deconv(y,x); q=2 -1 0.5, r=0 0 0 0.526下面几个文件用于转移函数与极零点之 间的相互转换及极零点的排序：（1) tf2zpk.m, 调用z,p,k=tf2zpk(b,a) 适用于z-1的升幂排列 tf2zp.m, 调用z,p,k=tf2zp(b,a) 适用于适用于z的降幂排列的降幂排列 (2) zp2tf.m, 调用 b,a=zp2tf (

23、z,p,k) (3)roots.m， 调用 Z1=roots(b) (4) poly.m, 调用b =poly (Z1)与极零点有关的与极零点有关的MATLAB函数函数27 显示离散系统的极零图显示离散系统的极零图函数函数：zplane.m本文件可用来显示离散系统的极零图。其调用格式是： zplane(z,p), 或 zplane(b,a),前者是在已知系统零点的列向量列向量z和极点的列向量列向量p的情况下画出极零图，后者是在仅已知A(z)、B(z) 的情况下画出极零图。28-1-0.500.511.52-1-0.500.5122Real PartImaginary Part1221222(

24、)1 44(44)/(2) /H zzzzzzzz 求极零点并绘极零分布图，部分画出幅频及相频响应：( )( )4 (1)4 (2)y nx nx nx n例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;291221222( )1 44(44)/(2) /H zzzzzzzz 画出幅频响应及相频响应：( )( )4 (1)4 (2)y nx nx nx n例1： 系统解：转

25、移函数：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)00.810510 Amplitude Freq. Res.00.81-505 wrap Phase Res.00.81-15-10-50unwrap Phase Res.30|()| 1jH e11( )H zzz例

26、例2：( ),02 相位的卷绕相位的卷绕 (wrapping) 解卷绕解卷绕 b=1;a=0,1;H,w=freqz(b,a,256,whole,1);Hr=abs(H);Hphase=angle(H);Hphase1=unwrap(Hphase); 050100150200250300111111050100150200250300-4-202405010015020025030002468画出幅频响应及相频响应：解：31-1-2-3-4-1-2-3-4.1836+.7344z +1.1016z +.7374z +.183116z( )1-3.0544z +3.8291z -2.2925z

27、+.55075z00H z 例：例： 给定系统给定系统求：极零图 单位抽样响应 频率响应 部分分式展开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解：r,p,k= residuez(b,a)b,a= residuez(r,p,k)32-1.5-1-0

28、.500.51-1-0.8-0.6-0.4-0.60.81Real PartImaginary Part极零图极零图 zplane(b,a); 330510152025303540-0.1-0.0500.00.25单位抽样响应h,t=impz(b,a,40);stem(t,h,.);grid on;34频率响应0100.511.5 Amplitude Freq. Res.01-505 wrap Phase Res.00.40.5

29、0.91-20-100unwrap Phase Res.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)35下面几个文件用于转移函数、极零点与二阶子系统二阶子系统sos(Second-Order Section)级联之间的相互转换：（1) tf2sos.m, 调用sos,G=tf2sos(b,a) (2) sos2tf.m,

30、 调用 b,a=sos2tf (sos,G) (3) sos2zp.m, 调用z,p,k= sos2zp (sos,G) (4) zp2sos.m, 调用 sos,G=zp2sos (z,p,k)与信号流图有关的与信号流图有关的MATLAB函数函数12,0,1,212,1,2( ),1,1kkkkkkzzHzkLazaz,0,1,2,1,21,1,kkkkkaakLsos是一是一L6的矩阵，每行元素如下排列：的矩阵，每行元素如下排列：361buttord.m 确定 LP DF、或 LP AF的阶次；（1） N, Wn = buttord(Wp, Ws, Rp, Rs)；（2）N, Wn = b

31、uttord(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的单位为弧度/秒，它们实际上是频率。372buttap.m 设计模拟低通(Butt)原型滤波器。 z, p, k=buttap(N

32、): N是欲设计的低通原型滤波器的阶次，z, p, k是设计出的极点、零点及增益。2211NccHjjjbuttap.m sets c to 1 for a normalized result. 383lp2lp.m、lp2hp.m、lp2bp.m, lp2bs.m将模拟低通原型转换为实际的低通、高通、带通及带阻滤波器。b, a 是AF LP 的分子、分母的系数向量，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,

33、 A=lp2bs(b, a, Wo , Bw)(2)(2)中Wo是带通或带阻滤波器中心频率，Bw是其带宽。394bilinear.m ：双线性变换，由模拟滤波器 得到数字滤波器。 Bz, Az=bilinear(B, A, Fs) 式中B, A分别是G（s）的分子、分母多项式的系数向量，Bz, Az分别是H(z)的分子、分母多项式的系数向量，Fs是抽样频率。405butter.m 用来直接设计Butterworth数字滤波器，实际上它把 buttord.m, buttap.m, lp2lp.m, bilinear.m等文件都包含了进去，从而使设计过程更简捷。格式(1)(3) 用来设计数字滤波器

34、，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); (3) B,A=butter(N,Wn,stop); (4) B,A=butter(N,Wn,s)416cheb1ord.m 求Cheb-型滤波器的阶；7cheb1ap.m 设计原型低Ch

35、eb-I型模拟滤波器；8cheby1.m 直接设计数字Cheb-滤波器。以上三个文件的调用格式和对应的Butterworth滤波器的文件类似。429cheb2ord.m; 10. cheb2ap.m;11. cheby2.m; 12. ellipord.m; 13. ellipap.m; 14. ellip.m;15impinvar.m 用冲激响应不变法实现频率转换； 对应 Cheby-II、椭圆 IIR 滤波器43 给出给出数字高通数字高通的技术要求的技术要求)(zHdhpspsp, 得到得到模拟高通模拟高通的技术要求的技术要求)(sHahpspsp,)2tan(1step得到得到模拟低通模

36、拟低通的技术要求的技术要求)(pGspsp,p12step设计出设计出)(pG3stepspp4step 得到得到模拟高通转移模拟高通转移函数函数)(sHahp11zzs5step最后得到最后得到数字高通转移数字高通转移函数函数)(zHdhp数字高通滤波器设计步骤数字高通, 带通及带阻滤波器的设计psqpq/144对 带通（BP）、带阻（BS）数字滤波器的设计，只需改变图中 Step2 和 Step4：22231BWBW 带阻31213()sps 22231BWBW 带通21331()sps 452N 121)(2pppG) 1() 1() 1(22222zzzpBW432142641. 03

37、07. 1237. 2637. 11)21 (0201. 0)(zzzzzzzHdbps3dB,18dB,F2000Hzps要求：按上述转换办法，可以求出：例 ： 设计一 IIR BP DF,要求：通带频率范围 ： 300Hz 400Hz ;阻带频率范围 ：200Hz、500Hz31 ffslshff46例例7.1uclear all; fp=100;fs=300;Fs=1000;rp=3;rs=20;uwp=2*pi*fp/Fs;ws=2*pi*fs/Fs;uFs=Fs/Fs; % let Fs=1 u% frequency prewarping ;uwap=2*Fs*tan(wp/2);w

38、as=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, wn) %ubz,az=bilinear(bs,as,Fs) % s=(2/Ts)(z-1)/(z+1);uh,w=freqz(bz,az,256,Fs*1000);uplot(w, 20*log10(abs(h);grid on;设计 IIR LP DF，3dB,20dBpssDF:100Hz,300Hz,1000HzpsffF47uclear all;

39、 wp=.2*pi;ws=.6*pi;Fs=1000; rp=3;rs=20;u% frequency prewarping;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) % z=(2/ts)(z-1)/(z+1);uh2,w2=freqz(b

40、z,az,500,Fs);uplot(w1/2/pi,abs(h1),w2,abs(h2),k);grid on;设计 IIR LP DF，3dB,20dBpssDF:100Hz,300Hz,1000HzpsffF例例7.1模拟滤波器和数字滤波器幅频特性比较48uclear all; 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,

41、az1,128,Fs);uplot(w, 20*log10(abs(h) ),w1, 20*log10(abs(h1) ),g.);ugrid on;设计 IIR LP DF，3dB,20dBpssDF:100Hz,300Hz,1000HzpsffF例例7.13dBp保证通带上限指标满足保证阻带下限指标满足20dBs49产生窗函数的文件有八个：1. bartlett（三角窗）; 2. blackman（布莱克曼窗） ; 3. boxcar（矩形窗）; 4. hamming（哈明窗）; 5. hanning（汉宁窗）; 6. triang（三角窗）;7. chebwin（切比雪夫窗）; 8 .k

42、aiser（凯赛窗）; 两端为零两端不为零调用方式都非常简单请见help文件稍为复杂509fir1.m 用“窗函数法”设计FIR DF。调用格式： （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带滤波器。51对格式（1

43、），若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); 指定矩形窗52例例7.2 用matlab设计一FIR滤波器，要求频率指标如下：0.1 0.25jwH ew1 0.05

44、 00.15jwH ew 0.1spwww10lg20log 0.0526pdBlg20sdB Solution:clc;clear all; close allwp=0.15*pi;ws=0.25*pi;wdelta=ws-wp;N=ceil(8*pi/wdelta);Wn=(0.15+0.25)*pi/2;b=fir1(N,Wn/pi,hanning(N+1);figure% freqz(b,1,512) H,w=freqz(b,1,512) ;plot(w,abs(H)set(gca,XTick,0:0.2*pi:pi)set(gca,XTickLabel,0,0.2,0.4,0.6,0

45、.8,)axis(0,pi,0.97,1.03);Hanning窗8wM8Mwlgminimum(,)26psdB 5310fir2.m 本文件采用“窗函数法”设计具有任意幅 频相应的FIR 数字滤波器。其调用格式是： b = fir2(N, F, M); F是频率向量，其值在01之间，M是和F相对应 的所希望的幅频相应。如同fir1, 缺省时自动选用 Hamming窗。高通DF，N为偶数，不能为奇数例例7.3 ：设计一多带滤波器，要求频率在0.20.3, 0.60.8 之间为1，其余处为零。设计结果如下：f = 0 0.2 0.2 0.3 0.3 0.6 0.6 0.8 0.8 1;m =

46、0 0 1 1 0 0 1 1 0 0;N1=30N2=90b1 = fir2(N1,f,m);b2 = fir2(N2,f,m);H1,w=freqz(b1,1,512); H2,w=freqz(b2,1,512); subplot(311)stem(b1)subplot(312)stem(b2)subplot(313)plot(f,m,w/pi,abs(H1),w/pi,abs(H2)5405101520253035-0.60102030405060708090100-0.4-0100.511.5

47、Comparison of Frequency Response Magnitudes IdealN=30 fir filterN=90 fir filterN=30,90时幅频响应响应及理想幅频响应；N=30N=90( )h n()jH e5511. 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是给定的

48、滤波器的阶次，b是设计的滤波器的系数，其长度为N1；F是频率向量，A是对应F的各频段上的理想幅频响应，W是各频段上的加权向量。56F、A及W的指定方式和例7.4.1和7.4.2所讨论过的一样，唯一的差别是F的范围为01，而非00.5, 1对应抽样频率的一半。需要指出的是，若b的长度为偶数，设计高通和带阻滤波器时有可能出现错误，因此，最好保证b的长度为奇数，也即N应为偶数。57例7.4.1： 设计低通 FIR DF:0.6p,0.7sb=remez(N, F, A, W)0.60.71clear all;f=0 .6 .7 1;% 给定频率轴分点；A=1 1 0 0;% 频率分点上理想幅频响应；

49、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的结果58例7.4.2： 设计多阻带滤波器，抽样频率500Hz, 在50Hz、 100Hz 及150Hz处陷波。 通带加权为8，阻带为1-17dB通带、阻带加权都是1-25dB59例7.4.2： 设计多阻带滤波器，抽样频率500Hz, 在

50、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,FontSi

51、ze,14,Color,r)60.220.340.380.420.460.580.620.540.66250Hz60例7.1.1.设计低通 DF FIR, 令截止频率 0. 25, 取 M10, 20，40，观察其幅频响应的特点.cclear 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/

52、M:0.5-0.5/M;plot(f,abs(h1),f,abs(h3),f,abs(h4);grid; axis(0 0.5 0 1.2)()jdHe22cc161clear 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; elseb

53、2(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-)/2( )(/ 2) ( )( 1)( )/ 2dn Mh nh nMw nw nnM例7.1.2： 理想差分器及其设计6212remezord.m 本文件用来确定在用Chebyshev最佳一致逼近设计FIR滤波器时所需要的滤波器阶次。其调用格式是： N, Fo, Ao,

54、 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的长度为奇数。 63u与本章内容有关的MAT

55、LAB文件主要是fft, ifft和 czt.m。顾名思义，fft实现快速傅立叶变换，ifft实现快速傅立叶反变换，czt.m 用来实现线性调频Z变换。ufft的调用格式是： X=fft(x), 或 X=fft(x, N)。u以高频pi为频谱中心 ufftshift（X）：移位使以零频为频谱中心与第八、九章有关的 MATLAB 文件64X=fft(x)fftshift(X)信号x65fftfilt.m 用叠接相加法实现长序列卷积。格式是: y=fftfilt(h,x) 或 y=fftfilt(h, x，N)记 的长度为 ， 的长度为 。 若采用第一个调用方式，程序自动地确定对 分段的长度 及做

56、FFT的长度 ， 显然， 是最接近 的2的整次幂。分的段数为 。( )x nxN( )h nM( )x n()LMNL/xNLN采用第二个调用方式，使用者可自己指定做FFT的长度。建议使用第一个调用方式。66例例。clear all; % 产生两个正弦加白噪声； 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(4,1,1); plot(x(1:N/4); axis tight f1=0:2*pi/N:

57、2*pi-pi/N; f=-pi:2*pi/N:pi-pi/N; X=fft(x); y=ifft(X); subplot(4,1,2); plot(f1, abs(X);axis tight subplot(4,1,3); plot(f, fftshift(abs(X);axis tight subplot(4,1,4); plot(real(y(1:N/4); axis tight67102030405060-5050123456200400-3-2-10123200400102030405060-505X=fft(x)fftshift(X)信号信号xx=ifft(X)68clear;% 用叠接相加法，计算滤波器系数用叠接相加法，计算滤波器系数h和输入信号和输入信号x的卷积的卷积% 其中其中h为为10阶阶hanning窗，窗，x是带有高斯白噪的正弦信号是带有高斯白噪的正弦信号h=fir1(10,0.3,hanning(11); % h: is the impulse response of a N=500;p=0.05;f=1/16; % 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