Biomedical Signal processing matlab 信号处理函数 Zhongguo Liu Biomedical Engineering School of Control Science and Engineering, Shandong University 2019/1/1 1 Zhongguo Liu_Biomedical Engineering_Shandong Univ.

关于MATLAB MATLAB是美国MathWorks公司开发的一种功能极其强大的高技术计算语言和内容极其丰富的软件库。它以矩阵和向量的运算以及运算结果的可视化为基础，把广泛应用于各个学科领域的数值分析、矩阵计算、函数生成、信号、图形及图象处理、建模与仿真等诸多强大功能集成在一个便于用户使用的交互式环境之中，为使用者提供了一个高效的编程工具及丰富的算法资源。

MATLAB与信号处理直接有关的工具箱 （ Toolbox） Signal Processing (信号处理工具箱） Wavelet (小波工具箱) Image Processing（图象处理工具箱） Higher-Order Spectral Analysis （高阶谱分析工具箱）

与信号处理间接有关的工具箱： Control System (控制系统) Communication （通信） System Identification （系统辨识） Statistics （统计） Neural Network (神经网络)

例： z=peaks; surf(z)；

如何改变 的方差为P 与第二章内容有关的MATLAB文件 ? 1. rand.m 用来产生均值为0.5、幅度在 0~1之间均匀分布的伪白噪声: u=rand(N,1) (rand(N)生成N阶矩阵) 方差： 方差函数var(u) 标准差函数std(u) 如何改变 的方差为P ? 即如何确定a使au的方差为P？ 将au代替u带入上面方差公式可得

与第二章内容有关的MATLAB文件 1. rand.m 用来产生均值为0.5、幅度在 0~1之间均匀分布的伪白噪声: u=rand(N,1) (rand(N)生成N阶矩阵) randn.m 用来产生均值为零、方差为1 服从高斯（正态）分布的白噪声信号 u=randn(1, N) x=randn(1000,1) y=randn(1000,1) 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 x=[3 4 5]; y=xcorr(x,h) y=24 53 86 65 38 10 -0 两序列的相关运算 MATLAB实现：y=xcorr(x1,x2)。

5．xcorr: 其互相关和自相关。格式是：(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)都除以（N－abs(m)）；若’flag’缺省，则rx不定标。M和‘flag’同样适用于求互相关。

第三章 Z变换 F=ztrans( f ) 对f(n)进行Z变换，其结果为F(z) . 在MATLAB语言中有专门对信号进行正反Z变换的函数ztrans( ) 和itrans( )。其调用格式分别如下： F=ztrans( f )      对f(n)进行Z变换，其结果为F(z) F=ztrans(f,w)   对f(n)进行Z变换，其结果为F(w) F=ztrans(f,k,w)    对f(k)进行Z变换,其结果为F(w) f=itrans ( F )    对F(z)进行Z反变换,其结果为f(n) f=itrans(F,k)   对F(z)进行Z反变换，其结果为f(k) f=itrans(F,w,k)  对F(w)进行Z反变换,其结果为f(k) 注意： 在调用函数ztran( )及iztran( )之前，要用syms命令对所有需要用到的变量（如t,u,v,w）等进行说明，即要将这些变量说明成符号变量

Z变换 例①.求数列 fn=e-n的Z变换及其逆变换。命令如下： syms n z fn=exp(-n); Fz=ztrans(fn,n,z) %求fn的Z变换 f=iztrans(Fz,z,n) %求Fz的逆Z变换 Fz= z/(z - 1/exp(1)) f = (1/exp(1))^n 例② ．用MATLAB求出离散序列f=0.5k的Z变换MATLAB程序如下： syms k z f=0.5^k;          %定义离散信号 Fz=ztrans(f)       %对离散信号进行Z变换 运行结果如下： Fz = 2*z/(2*z-1)

Z变换 例③ ．已知一离散信号的Z变换式为 , 求出它所对应的离散信号f(k). MATLAB程序如下： syms k z Fz=2*z/(2*z-1);       %定义Z变换表达式 fk=iztrans(Fz,k)        %求反Z变换 运行结果如下：fk = (1/2)^k 阶跃序列符号 例④: 求序列的Z变换. syms k hk=sym('kroneckerDelta(k, 1) + kroneckerDelta(k, 2)+ kroneckerDelta(k, 3)') Hz=ztrans(hk) Hz=simplify(Hz) 运行结果如下：Fz= (z^2 + z + 1)/z^3

符号变换 Fourier变换及其反变换 Fw=fourier(ft,t,w) 求“时域”函数ft的Fourier变换 ft=ifourier(Fw,w,t) 求“频域”函数Fw的Fourier 反变换 Laplace变换及其反变换 Fs=laplace(ft,t,s) 求“时域”函数ft的Laplace变换 ft=ilaplace(Fs,s,t) 求“频域”函数Fs的Laplace

【例 】求 的Fourier变换。 （1）求Fourier变换 syms t w ut=heaviside(t); UT=fourier(ut) UT = pi*dirac(w)-i/w （2）求Fourier反变换 Ut=ifourier(UT,w,t) Ut =heaviside(t)

【例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),t^2*exp(-t)]; MS=laplace(Mt,t,s) MS = [ exp(-s*a), exp(-s*b)/s] [ 1/b/((s+a)^2/b^2+1), 2/(s+1)^3]

与系统响应、逆Z变换 相关的matlab 函数 1．filter.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 函数 x=[1,2,3,4]; y=[3,4,6] 1．filter.m x=[1,2,3,4]; y=[3,4,6] z_conv= conv(x,y) % x , y 为输入和单位脉冲响应时输出 z_conv_= conv(y, x) % x , y为输入和单位脉冲响应时输出 z_filter=filter(y,1,x) % x为输入, y为FIR单位脉冲响应时输出 z_filter_=filter(x,1,y) % y为输入, x为FIR单位脉冲响应时输出 z_conv = 3 10 23 36 34 24  z_conv_=3 10 23 36 34 24 z_filter = 3 10 23 36 z_filter_= 3 10 23 可见，conv(x,y)总是等于conv(y,x)。而filter(x,1,y)却不一定等于filter(y,1,x)，但是它们总是conv(x,y)截短的结果，截短的长度等于length(filter的第三个参数)

与逆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) 的有理分式分解成简单有理分式的和，因此可用来求逆Z变换。调用格式： [r,p,k]= residuez(b,a) 假如知道了向量r, p和k，利用residuez.m还可反过来求出多项式A(z)、B(z)。格式是 [b,a]= residuez(r,p,k)。

[H,w]=freqz(b,a,N,'whole',Fs) 4.频率响应函数：freqz.m 已知A(z)、B(z), 求系统的频率响应。基本的调用格式是： [H,w]=freqz(b,a,N,'whole',Fs) N是频率轴的分点数，建议N为2的整次幂；w是返回频率轴座标向量，绘图用；Fs是抽样频率，若Fs＝1，频率轴给出归一化频率；’whole’指定计算的频率范围是从0～FS,缺省时是从0～FS/2. 幅频响应：Hr=abs(H); Hphase=angle(H); Hphase=unwrap(Hphase);

5．filtfilt.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 作零相位滤波； 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');

5．filtfilt.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 作零相位滤波； 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'); 23

6．grpdelay.m 求系统的群延迟。调用格式 [gd w]=grpdelay(B, A, N) , 或 [gd F]=grpdelay(B, A, N, FS) 式中B和A仍是 的分子、分母多项式，gd是群延迟，w、F是频率分点，二者的维数均为N；FS为抽样频率，单位为Hz。

7．deconv.m ：实现系统的反卷积，其调用格式： [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');ylabel(' q(n)'); subplot(326);stem(q1,'.r');ylabel(' q1(n)'); clear all; % 实现多项式除法q=(4*x^3+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.5]

（1) tf2zpk.m, 调用[z,p,k]=tf2zpk(b,a) 适用于z-1的升幂排列 与极－零点有关的MATLAB函数 下面几个文件用于转移函数与极－零点之 间的相互转换及极－零点的排序： （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 (z,p,k) (3)roots.m， 调用 Z1=roots(b) (4) poly.m, 调用b =poly (Z1)

显示离散系统的极零图函数：zplane.m 本文件可用来显示离散系统的极－零图。其调用格式是： zplane(z,p), 或 zplane(b,a), 前者是在已知系统零点的列向量z和极点的列向量p的情况下画出极－零图，后者是在仅已知A(z)、B(z) 的情况下画出极－零图。

例1： 系统 解： 求极零点并绘极零分布图，部分画出幅频及相频响应： 转移函数： z=2; 2 p=0; 0 b=[1 -4 4]; K=1 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) r =[]; 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： Hphase1=unwrap(Hphase); 解： 相位的卷绕 (wrapping) 解卷绕

例： 给定系统 ? 解： 求：极－零图 单位抽样响应 频率响应 部分分式展开 stem(t,h,'.');grid on; b=[.1836 .7344 1.1016 .7374 .1836]/100 a =[1 -3.0544 3.8291 -2.2925 .55075] zplane(b,a); ? [h,t]=impz(b,a,40); stem(t,h,'.');grid on; Hphase=unwrap(Hphase); [r,p,k]= residuez(b,a) [b,a]= residuez(r,p,k)

zplane(b,a); 极－零图

[h,t]=impz(b,a,40); stem(t,h,'.');grid on; 单位抽样响应

频率响应 Hr=abs(H); subplot(311), plot(Hr) subplot(312), plot(Hphase) plot(Hphaseu) Hphaseu=unwrap(Hphase); 频率响应

与信号流图有关的MATLAB函数 下面几个文件用于转移函数、极－零点与二阶子系统sos(Second-Order Section)级联之间的相互转换： （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) sos是一L×6的矩阵，每行元素如下排列：

与第七章内容有关的MATLAB文件 1．buttord.m 确定 LP DF、或 LP AF的阶次； （1） [N, Wn] = buttord(Wp, Ws, Rp, Rs)； （2）[N, Wn] = buttord(Wp, Ws, Rp, Rs,‘s’)： (1)对应数字滤波器。其中 Wp, Ws分别是通带和阻带的截止频率，其值在 0～1 之间，1对应抽样频率的一半(归一化频率)。对低通和高通，Wp, Ws都是标量，对带通和带阻，Wp, Ws是1×2的向量。Rp, Rs 分别是通带和阻带的衰减(dB)。N是求出的相应低通滤波器的阶次，Wn是求出的3dB频率，它和Wp稍有不同。 (2)对应模拟滤波器，各变量含意和（1）相同，但Wp, Ws及Wn的单位为弧度/秒，它们实际上是频率。

2．buttap.m 设计模拟低通(Butt)原型滤波器。 [z, p, k]=buttap(N): N是欲设计的低通原型滤波器的阶次，z, p, k是设计出的极点、零点及增益。 buttap.m sets Ωc to 1 for a normalized result.

3．lp2lp.m、lp2hp.m、lp2bp.m, lp2bs.m 将模拟低通原型转换为实际的低通、高通、带通及带阻滤波器。 [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) b, a 是AF LP 的分子、分母的系数向量，B, A是转换后的的分子、分母的系数向量；(1)中，Wo是低通或高通滤波器的截止频率； (2)中Wo是带通或带阻滤波器中心频率，Bw是其带宽。

4．bilinear.m ：双线性变换，由模拟滤波器 得到数字滤波器。 [Bz, Az]=bilinear(B, A, Fs) 式中B, A分别是G（s）的分子、分母多项式 的系数向量，Bz, Az分别是H(z)的分子、分 母多项式的系数向量，Fs是抽样频率。

5．butter. m 用来直接设计Butterworth数字滤波器，实际上它把 buttord. m, buttap. m, lp2lp 5．butter.m 用来直接设计Butterworth数字滤波器，实际上它把 buttord.m, buttap.m, lp2lp.m, bilinear.m等文件都包含了进去，从而使设计过程更简捷。 [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’) 格式(1)～(3) 用来设计数字滤波器，B，A分别是H(z)的分子、分母多项式的系数向量，Wn是通带截止频率，范围在0～1之间。若Wn是标量，(1)用来设计低通数字滤波器，若Wn是1×2的向量，则(1) 用来设计数字带通滤波器； (2)用来设计数字高通滤波器； (3) 用来设计数字带阻滤波器，显然，这时的Wn是1×2的向量；格式(4) 用来设计模拟滤波器。

6．cheb1ord.m 求Cheb-Ⅰ型滤波器的阶； 7．cheb1ap.m 设计原型低Cheb-I型模拟滤波器； 8．cheby1.m 直接设计数字Cheb-Ⅰ滤波器。 以上三个文件的调用格式和对应的 Butterworth滤波器的文件类似。

9．cheb2ord.m; 10. cheb2ap.m; 11. cheby2.m; 12. ellipord.m; 13. ellipap.m; 14. ellip.m; 对应 Cheby-II、椭圆 IIR 滤波器 15．impinvar.m 用冲激响应不变法实现频率转换；

数字高通, 带通及带阻滤波器的设计 数字高通滤波器设计步骤 给出 得到 得到 模拟低通 的技术要求 数字高通 的技术要求 模拟高通 设计出 最后得到 数字高通转移 函数 得到 模拟高通转移 函数 数字高通滤波器设计步骤

对 带通（BP）、带阻（BS）数字滤波器的设计，只需改变图中 Step2 和 Step4：

例 ： 设计一 IIR BP DF,要求： 要求： 通带频率范围 ： 300Hz~ 400Hz ; 阻带频率范围 ：200Hz、500Hz 按上述转换办法，可以求出：

例7.1 设计 IIR LP DF， clear all; fp=100;fs=300;Fs=1000;rp=3;rs=20; wp=2*pi*fp/Fs;ws=2*pi*fs/Fs; Fs=Fs/Fs; % let Fs=1 % frequency prewarping ; wap=2*Fs*tan(wp/2);was=2*Fs*tan(ws/2); % [n,wn]=buttord(wap,was,rp,rs,'s') % Note: 's'! [z,p,k]=buttap(n); % [bp,ap]=zp2tf(z,p,k) % [bs,as]=lp2lp(bp,ap, wn) % [bz,az]=bilinear(bs,as,Fs) %% s=(2/Ts)(z-1)/(z+1); [h,w]=freqz(bz,az,256,Fs*1000); plot(w, 20*log10(abs(h)));grid on;

clear all; wp=.2*pi;ws=.6*pi;Fs=1000; rp=3;rs=20; 例7.1 设计 IIR LP DF， clear all; wp=.2*pi;ws=.6*pi;Fs=1000; rp=3;rs=20; % frequency prewarping; wap=2*Fs*tan(wp/2);was=2*Fs*tan(ws/2); [n,wn]=buttord(wap,was,rp,rs,'s');% Note: 's'! [z,p,k]=buttap(n); [bp,ap]=zp2tf(z,p,k); [bs,as]=lp2lp(bp,ap,wap) w1=[0:499]*2*pi; h1=freqs(bs,as,w1); [bz,az]=bilinear(bs,as,Fs) % z=(2/ts)(z-1)/(z+1); [h2,w2]=freqz(bz,az,500,Fs); plot(w1/2/pi,abs(h1),w2,abs(h2),'k');grid on; 模拟滤波器和数字滤波器幅频特性比较

例7.1 设计 IIR LP DF， clear all; wp=.2*pi;ws=.6*pi;Fs=1000; rp=3;rs=20; [n,wn]=buttord(wp/pi,ws/pi,rp,rs); [bz,az]=butter(n,wp/pi) [bz1,az1]=butter(n,wn) [h,w]=freqz(bz,az,128,Fs); [h1,w1]=freqz(bz1,az1,128,Fs); plot(w, 20*log10(abs(h) ),w1, 20*log10(abs(h1) ),'g.'); grid on; 保证通带上限指标满足 保证阻带下限指标满足

产生窗函数的文件有八个： bartlett（三角窗）; 2. blackman（布莱克曼窗） ; 两端为零 3. boxcar（矩形窗）; 4. hamming（哈明窗）; 5. hanning（汉宁窗）; 6. triang（三角窗）; 7. chebwin（切比雪夫窗）; 8 .kaiser（凯赛窗）; 两端为零 调用方式都非常简单请见help文件 两端不为零 稍为复杂

9．fir1.m 用“窗函数法”设计FIR DF。调用格式： （1）b = fir1(N,Wn); (2) b = fir1(N,Wn,‘high’); (3) b = fir1(N,Wn, ‘stop’); N：阶次，滤波器长度为N＋1；Wn：通带截止频率，其值在0～1之间，1对应 Fs/2； b: 滤波器系数。 格式（2）用来设计高通滤波器， 格式（3）用来设计带阻滤波器。 格式（1），若Wn为标量，则设计低通滤波器，若 Wn是1×2的向量，则用来设计带通滤波器，若Wn是 1×L的向量，则可用来设计L带滤波器。

对格式（1），若Wn为标量，则设计低通滤波器，若 Wn是1×2的向量，则用来设计带通滤波器，若Wn是 1×L的向量，则可用来设计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.2 用matlab设计一FIR滤波器，要求频率指标如下： Solution: clc;clear all; close all wp=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.8π','π'}) axis([0,pi,0.97,1.03]); Hanning窗

10．fir2.m 本文件采用“窗函数法”设计具有任意幅 频相应的FIR 数字滤波器。其调用格式是： b = fir2(N, F, M); F是频率向量，其值在0～1之间，M是和F相对应 的所希望的幅频相应。如同fir1, 缺省时自动选用 Hamming窗。高通DF，N为偶数，不能为奇数 例7.3 ：设计一多带滤波器，要求频率在0.2~0.3, 0.6~0.8 之间为1，其余处为零。设计结果如下： f = [0 0.2 0.2 0.3 0.3 0.6 0.6 0.8 0.8 1]; m = [0 0 1 1 0 0 1 1 0 0]; N1=30 N2=90 b1 = 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))

N=30 N=90 N=30,90时幅频响应响应及理想幅频响应；

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是设计的滤波器的系数，其长度为N＋1；F是频率向量，A是对应F的各频段上的理想幅频响应，W是各频段上的加权向量。

F、A及W的指定方式和例7.4.1和7.4.2所讨论过 的一样，唯一的差别是F的范围为0～1，而非0～0.5, 1对应抽样频率的一半。需要指出的是，若b的长度为偶数，设计高通和带阻滤波器时有可能出现错误，因此，最好保证b的长度为奇数，也即N应为偶数。

例7.4.1： 设计低通 FIR DF: 调整通带、阻带的加权及滤波器的长度。 b=remez(N, F, A, W) 调整N或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')

例7.1.1.设计低通 DF FIR, 令截止频率 ＝0. 25π, 取 M＝10, 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 b2(n)=0; b2(n)=win(n)*(-1)^(n-1-N/2)/(n-1-N/2); 例7.1.2： 理想差分器及其设计 M=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-')

12．remezord.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的长度为奇数。

与第八、九章有关的 MATLAB 文件 与本章内容有关的MATLAB文件主要是fft, ifft和 czt.m。顾名思义，fft实现快速傅立叶变换，ifft实现快速傅立叶反变换，czt.m 用来实现线性调频Z变换。 fft的调用格式是： X=fft(x), 或 X=fft(x, N)。 以高频pi为频谱中心 fftshift（X）：移位使以零频为频谱中心

信号x X=fft(x) fftshift(X)

fftfilt.m 用叠接相加法实现长序列卷积。格式是: y=fftfilt(h,x) 或 y=fftfilt(h, x，N) 记 的长度为 ， 的长度为 。 若采用第一个调用方式，程序自动地确定对 分段的长度 及做FFT的长度 ， 显然， 是最接近 的2的整次幂。分的段数为 。 采用第二个调用方式，使用者可自己指定做FFT的长度。建议使用第一个调用方式。

例1. 设x(n)是两个正弦信号及白噪声的叠加，进行频谱分析。 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: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 tight

信号x X=fft(x) fftshift(X) x’=ifft(X)

例2 令x(n)为一正弦加白噪声信号，长度为500，h(n)是用fir1. m文件设计的一个低通FIR滤波器，长度为11 例2 令x(n)为一正弦加白噪声信号，长度为500，h(n)是用fir1.m文件设计的一个低通FIR滤波器，长度为11.试用fftfilt实现长序列的卷积 clear; % 用叠接相加法，计算滤波器系数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 noise s=sin(2*pi*f*[0:N-1]); % s:sine signal x=u(1:N)+s; % x: a long sequence; y=fftfilt(h,x); % Overlap-add method of FIR filtering using FFT % y=filter(h,1,x); % y=x*h subplot(211); plot(x); subplot(212); plot(y); figure; subplot(211); plot(abs(fft(x))); subplot(212); plot(abs(fft(y)));

信号x 滤波输出信号y 信号x频谱 信号y频谱

clear; % 按照overlap-add算法编程实现长信号滤波输出 % 生成滤波器系数h和混有高斯白噪的正弦信号x h=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)]; end subplot(211);plot(x) subplot(212);plot(y(1:N))

czt.m 调用格式是： X＝czt(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);

例3 设x(n)由三个实正弦组成，频率分别是8HZ， 8.22HZ 和9HZ， 抽样频率是40HZ ，时域取128点。 用CZT计算的DFT 用FFT计算的DFT 7~(7+M×0.05) HZ 的CZT变换

程序 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;

% DFT Y2=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;

hilbert.m 文件用来计算信号Hilbert变换。调用 的格式是: y=hilbert(x),y的实部就 所以，y 实际上是 x 的解析信号。