Transcript 实验指导-信号处理MATLAB函数介绍

Biomedical Signal processing
matlab 信号处理函数
Zhongguo Liu
Biomedical Engineering
School of Control Science and Engineering, Shandong
University
2020/4/30
1
Zhongguo Liu_Biomedical Engineering_Shandong Univ.
关于MATLAB
MATLAB是美国MathWorks公司开发的
一种功能极其强大的高技术计算语言和
内容极其丰富的软件库。它以矩阵和向
量的运算以及运算结果的可视化为基础，
把广泛应用于各个学科领域的数值分析、
矩阵计算、函数生成、信号、图形及图
象处理、建模与仿真等诸多强大功能集
成在一个便于用户使用的交互式环境之
中，为使用者提供了一个高效的编程工
具及丰富的算法资源。
2
MATLAB与信号处理直接有关的工具箱
（ Toolbox）
Signal Processing (信号处理工具箱）
Wavelet
(小波工具箱)
Image Processing（图象处理工具箱）
Higher-Order Spectral Analysis
（高阶谱分析工具箱）
3
与信号处理间接有关的工具箱：
Control System
Communication
(控制系统)
（通信）
System Identification （系统辨识）
Statistics
Neural Network
4
（统计）
(神经网络)
例：
5
z=peaks; surf(z)；
与第二章内容有关的MATLAB文件
1. rand.m 用来产生均值为0.5、幅度在
0~1之间均匀分布的伪白噪声: u=rand(N,1)
(rand(N)生成N阶矩阵)
1
方差：  
N
2
u
2
N 1
 u ( n )   
n 0
方差函数var(u)
u
1
,  
12
2
u
标准差函数std(u)
如何改变 u(n) 的方差
u(n )
6
与第二章内容有关的MATLAB文件
1. rand.m 用来产生均值为0.5、幅度在
0~1之间均匀分布的伪白噪声: u=rand(N,1)
(rand(N)生成N阶矩阵)
2. randn.m 用来产生均值为零、方差为1
服从高斯（正态）分布的白噪声信号
u=randn(1, N)
x=randn(1000,1)
y=randn(1000,1)
v=var(x)
7
h=std(y)
u(n )
3.sinc ：用来产生 “sinc” 函数：
sinc函数定义为：
t 0
1

sin c(t )  0
t  k
sin( t )  t t为其它

t=-4:0.1:4;
x4=sinc(t);
plot(t,x4)
8
%
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];
x=[3 4 5];
h=[2 6 7 8]; h=[2 6 7 8];
y=conv(x,h)
y=xcorr(x,h)
运行结果:
y= 6 26 55 82 67 40
y =24 53 86 65 38 10 -0
y (m)   x1 (n) x2 (n  m)
两序列的相关运算
n
MATLAB实现：y=xcorr(x1,x2)。
9
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’同样适用于求互相关。
10
第三章 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=itrans ( F )
对f(u)进行Z变换,其结果为F(v)
对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）等进
行说明，即要将这些变量说明成符号变量
11
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变换
例② ．用MATLAB求出离散序列 的Z变换MATLAB
程序如下：
syms k z
f=0.5^k;
%定义离散信号
Fz=ztrans(f)
%对离散信号进行Z变换
运行结果如下：
Fz
=
2*z/(2*z-1)
12
Z变换
2z
F ( z) 
2z 1
例③ ．已知一离散信号的Z变换式为
,
求出它所对应的离散信号f(k). MATLAB程序如下：
syms k z
Fz=2*z/(2*z-1);
%定义Z变换表达式
fk=iztrans(Fz,k)
%求反Z变换
运行结果如下：fk = (1/2)^k
例④：求序列的Z变换. f (k )   (k  1)   (t  4)
syms n
hn=sym('kroneckerDelta(n, 1) +
kroneckerDelta(n, 2)+ kroneckerDelta(n, 3)')
Hz=ztrans(hn)
Hz=simplify(Hz)
运行结果如下：Fz= (z^2 + z + 1)/z^3
13
与逆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都是向量。
14
与逆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开始。
15
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)。
16
频率响应函数：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);
17
解卷绕：
B( z )
H ( z) 
A( z )
与极－零点有关的MATLAB函数
下面几个文件用于转移函数与极－零点之
间的相互转换及极－零点的排序：
（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)
18
显示离散系统的极零图函数：zplane.m
本文件可用来显示离散系统的极－零图。
其调用格式是：
zplane(z,p), 或
zplane(b,a),
前者是在已知系统零点的列向量z和极点的
列向量p的情况下画出极－零图，后者是在
仅已知A(z)、B(z) 的情况下画出极－零图。
19
例1： 系统 y (n)  x(n)  4 x(n  1)  4 x(n  2)
用极零分析大致画出幅频响应及相频响应：
解： 转移函数：
1
H ( z)  1  4 z  4 z
z=2; 2  ( z 2  4 z1  4) / z 2
b=[1 -4 4]; p=0; 0
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)
2
 ( z  2)2 / z 2
1
Imaginary Part
0.5
2
0
2
-0.5
r =[]; p=[];
k=[1 -4 4];
-1
-1
20
-0.5
0
0.5
Real Part
1
1.5
2
例1： 系统 y (n)  x(n)  4 x(n  1)  4 x(n  2)
用极零分析大致画出幅频响应及相频响应：
Amplitude Freq. Res.
[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)
21
2
 ( z 2  4 z1  4) / z 2  ( z  2)2 / z 2
wrap Phase Res.
b=[1 -4 4];
a=[1];
H ( z)  1  4 z  4 z
unwrap Phase Res.
解： 转移函数：
1
10
5
0
0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
0.6
0.8
1
5
0
-5
0
-5
-10
-15
例2：H ( z )  z  1
z 1
j
| H (e ) | 1

 ( )   ,   0  2
b=[1];
a=[0,1];
[H,w]=freqz(b,a,256,'whole',1);
Hr=abs(H);
1
1
1
1
Hphase=angle(H);
1
1
相位的卷绕
(wrapping)
0
50
100
150
200
250
300
4
2
Hphase1=unwrap(Hphase);
0
-2
-4
0
50
100
150
200
250
300
0
50
100
150
200
250
300
8
解卷绕
6
4
2
0
22
例： 给定系统
-1
-2
-3
-4
1 .1836+.7344z +1.1016z +.7374z +.1836z
H ( z) 
-1
-2
-3
-4
100 1-3.0544z +3.8291z -2.2925z +.55075z
求：极－零图 解：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;
[H,w]=freqz(b,a,256,'whole',1);
Hr=abs(H);
23
1
0.8
0.6
Imaginary Part
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
-1.5
24
-1
-0.5
0
Real Part
zplane(b,a); 极－零图
0.5
1
0.25
0.2
0.15
0.1
0.05
0
-0.05
-0.1
0
5
[h,t]=impz(b,a,40);
10
15
20
25
30
单位抽样响应
stem(t,h,'.');grid on;
25
35
40
Amplitude Freq. Res.
[H,w]=freqz(b,a,256,'whole',1);
Hr=abs(H);
1.5
1
0.5
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
unwrap Phase Res.
wrap Phase Res.
subplot(311), 0
0
plot(Hr)
5
subplot(312),
plot(Hphase) 0
subplot(313),
-5
plot(Hphaseu) 0
0
-10
-20
0
Hphase=angle(H);
26
频率响应
与信号流图有关的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的矩阵，每行元素如下排列：
[  k ,0  k ,1  k ,2 1 ak ,1 ak ,2 ], k  1,
 k ,0   k ,1 z   k ,2 z
1
H k ( z) 
27
1
1  ak ,1 z  ak ,2 z
2
,L
2
,
k  1,
,L
与本章内容有关的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的单位为弧度/秒，它们实际上是频率。
28
2．buttap.m 设计模拟低通(Butt)原型滤波器。
[z, p, k]=buttap(N): N是欲设计的低通原型滤波
器的阶次，z, p, k是设计出的极点、零点及增益。
3．lp2lp.m、lp2hp.m、lp2bp.m, lp2bs.m
将模拟低通原型转换为实际的低通、高通、带通及带阻
滤波器。
[B, A]=lp2bp(b, a, Wo , Bw)
[B, A]=lp2lp(b, a, Wo)
(2)
(1)
[B, A]=lp2bs(b, a, Wo , Bw)
[B, A]=lp2hp(b, a, Wo)
b, a 是AF LP 的分子、分母的系数向量，B, A是转换后的
的分子、分母的系数向量；(1)中，Wo是低通或高通滤波
器的截止频率；
(2)中Wo是带通或带阻滤波器中心频率，Bw是其带宽。
29
4．bilinear.m ：双线性变换，由模拟滤波器
得到数字滤波器。
[Bz, Az]=bilinear(B, A, Fs)
式中B, A分别是G（s）的分子、分母多项式
的系数向量，Bz, Az分别是H(z)的分子、分
母多项式的系数向量，Fs是抽样频率。
30
5．butter.m 用来直接设计Butterworth数字滤波器，实际
上它把 buttord.m, buttap.m, lp2lp.m, bilinear.m等文件都包
含了进去，从而使设计过程更简捷。
(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’)
格式(1)～(3) 用来设计数字滤波器，B，A分别是H(z)的
分子、分母多项式的系数向量，Wn是通带截止频率，
范围在0～1之间。若Wn是标量，(1)用来设计低通数字
滤波器，若Wn是1×2的向量，则(1) 用来设计数字带通
滤波器； (2)用来设计数字高通滤波器； (3) 用来设计数
字带阻滤波器，显然，这时的Wn是1×2的向量；格式(4)
用来设计模拟滤波器。
31
 s  20dB
DF : f p  100 Hz, f s  300 Hz, Fs  1000 Hz
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
% Firstly to finish frequency prewarping ;
wap=tan(wp/2);was=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) %
% Note: s=(2/Ts)(z-1)/(z+1);Ts=1,that is
2Fs=1,Fs=0.5;
32
例6.7.1(例6.5.1)
设计 IIR LP DF， p  3dB,
例6.7.1(例6.5.1)
设计 IIR LP DF， p  3dB,
 s  20dB
f s  300 Hz, Fs  1000 Hz
DF : f p  100 Hz,
clear all;
wp=.2*pi;ws=.6*pi;Fs=1000; rp=3;rs=20;%
% Firstly to finish 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)
% Note: z=(2/ts)(z33
例6.7.1(例6.5.1)
DF :
设计 IIR LP DF， p  3dB,
f p  100 Hz,
 s  20dB
f s  300 Hz, Fs  1000 Hz
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,abs(h),w1,abs(h1),'g.');grid on;
34
非线性关系
2
  tan( / 2)
Ts
又称为频率的预变形（Freq. Warping) 。例如 ：
DF :
f p  100 Hz,
f s  300 Hz, Fs  1000 Hz
 p  0.2 ,
 s  0.6 ,
抽样频率
2
AF :  p  tan( p / 2)  685.8  2 109(Hz)
Ts
 2 . f p
2
 s  tan( s / 2)  2452.76  2  438(Hz)
Ts
 2 . f s
设计的 AF 并不是按给定的技术指标，但当再
由 s 变回 z 后，保证了 DF的技术要求。
35
7.6 数字高通, 带通及带阻滤波器的设

计 给出


得到
得到

数字高通 H dhp (z )   tan(

的技术要求
step1
 p ,  s , p , s
2
p
) 模拟高通 H ahp (s )

的技术要求
数字高通转移
函数
H dhp (z )
s
z 1
得到
z  1 模拟高通转移
step5
函数 H ahp (s )

p
 p , s , p ,  s
q  s / p
p
s
step4
设计出
G( p)
数字高通滤波器设计步骤
36
模拟低通 G ( p )
的技术要求
step2
 p , s ,  p ,  s
q  1/ p
最后得到
1
step3
对 带通（BP）、带阻（BS）数字滤波器的设
计，只需改变图中 Step2 和 Step4：
 BW  3  1
 BW  3  1
    BW
    BW
  2


2
s  13
p
s (3  1 )
2
37
2
带
通

 2
2
  2
s (3  1 )
p 2
s  13
带
阻
例6.6.2： 设计一 IIR BP DF,要
求：通带频率范围 f1 ~ f 3 ： 300Hz~ 400Hz ;
阻带频率范围 f sl ~ f sh ：200Hz、500Hz
要求：  p  3dB,  s  18dB, Fs  2000 Hz
N 2
按上述转换办法，可以求出：
1
G( p)  2
p  2 p 1
( z  1)  2 ( z  1)
p
 BW ( z 2  1)
2
2
0.0201(1  2 z 2  z 4 )
H dbp ( z ) 
1  1.637 z 1  2.237 z 2  1.307 z 3  0.641z 4
38
2
与本章内容有关的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的单位为弧度/秒，它们实际上是频率。
39
2．buttap.m 设计模拟低通(Butt)原型滤波器。
[z, p, k]=buttap(N): N是欲设计的低通原型滤波
器的阶次，z, p, k是设计出的极点、零点及增益。
3．lp2lp.m、lp2hp.m、lp2bp.m, lp2bs.m
将模拟低通原型转换为实际的低通、高通、带通及带阻
滤波器。
[B, A]=lp2bp(b, a, Wo , Bw)
[B, A]=lp2lp(b, a, Wo)
(2)
(1)
[B, A]=lp2bs(b, a, Wo , Bw)
[B, A]=lp2hp(b, a, Wo)
b, a 是AF LP 的分子、分母的系数向量，B, A是转换后的
的分子、分母的系数向量；(1)中，Wo是低通或高通滤波
器的截止频率；
(2)中Wo是带通或带阻滤波器中心频率，Bw是其带宽。
40
4．bilinear.m ：双线性变换，由模拟滤波器
得到数字滤波器。
[Bz, Az]=bilinear(B, A, Fs)
式中B, A分别是G（s）的分子、分母多项式
的系数向量，Bz, Az分别是H(z)的分子、分
母多项式的系数向量，Fs是抽样频率。
41
5．butter.m 用来直接设计Butterworth数字滤波器，实际
上它把 buttord.m, buttap.m, lp2lp.m, bilinear.m等文件都包
含了进去，从而使设计过程更简捷。
(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’)
格式(1)～(3) 用来设计数字滤波器，B，A分别是H(z)的
分子、分母多项式的系数向量，Wn是通带截止频率，
范围在0～1之间。若Wn是标量，(1)用来设计低通数字
滤波器，若Wn是1×2的向量，则(1) 用来设计数字带通
滤波器； (2)用来设计数字高通滤波器； (3) 用来设计数
字带阻滤波器，显然，这时的Wn是1×2的向量；格式(4)
用来设计模拟滤波器。
42
6．cheb1ord.m
求Cheb-Ⅰ型滤波器的阶；
7．cheb1ap.m
设计原型低Cheb-I型模拟滤波器；
8．cheby1.m
直接设计数字Cheb-Ⅰ滤波器。
以上三个文件的调用格式和对应的
Butterworth滤波器的文件类似。
43
9．cheb2ord.m;
10. ellipord.m;
11.cheb2ap.m;
12. ellipap.m;
13.besselap.m;
对应 Cheby-II、
椭圆 IIR 滤波器
14. cheby2.m;
15. ellip.m;
16.besself.m
17．impinvar.m 用冲激响应不变法实现频率转换；
44
产生窗函数的文件有八个：
两端为零
1. bartlett（三角窗）;
2. blackman（布莱克曼窗） ;
调用方
式都非
常简单
请见help
文件
3. boxcar（矩形窗）;
4. hamming（哈明窗）;
5. hanning（汉宁窗）;
6. triang（三角窗）;
两端不为零
7. chebwin（切比雪夫窗）;
8 .kaiser（凯赛窗）;
45
稍为复杂
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带滤波器。
46
对格式（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)); 指定矩形窗
47
10．fir2.m
幅
本文件采用“窗函数法”设计具有任意
频相应的FIR 数字滤波器。其调用格式是：
b = fir2(N, F, M);
F是频率向量，其值在0～1之间，M是和F相对应
的所希望的幅频相应。如同fir1, 缺省时自动选用
Hamming窗。
例 ：设计一多带滤波器，要求频率在0.2~0.3,
0.6~0.8 之间为1，其余处为零。
设计结果如下：
48
0.5
N=30
0
-0.5
0
5
10
15
20
25
30
35
0.5
N=90
0
-0.5
h( n)
0
20
40
60
80
100
1
j
0.5
0
H (e )
0
0.1
0.2
0.3
0.4
0.5
N=30,90时幅频响应响应及理想幅频响应；
49
13 firls.m 用最小平方法设计线性相位FIR滤波器，
可设计任意给定的理想幅频响应；
14 fircls.m用带约束的最小平方法设计线性相位
FIR滤波器，可设计任意给定的理想幅频响应；
15 fircls1.m 用带约束的最小平方方法设计线性
相位FIR低通和高通滤波器。
16 sgolay.m 用来设计 Savitzky-Golay FIR 平滑滤
波器，其原理见9.1.1节
17 firrcos.m 用来设计低通线性相位FIR滤波器，
其过渡带为余弦函数形状。
50
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带滤波器。
51
对格式（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)); 指定矩形窗
52
 c ＝0. 25π,
例7.1.1.设计低通 DF FIR, 令截止频率
取 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);
2
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])
53
H d (e j )
1


c
c

2
clear all;N=40;n=0:N;
例7.1.2： 理想差
b1=fir1(N,0.25,boxcar(N+1));
分器及其设计
b2=fir1(N,0.25,hamming(N+1));
win=hamming(N+1);
h(n)  hd (n  M / 2) w(n)
for n=1:N+1
if (n-1-N/2)==0;
(1) n  M /2

w(n)
b1(n)=0;
nM /2
else
M=128;
b1(n)=(-1)^(n-1-N/2)/(n-1-N/2);
h1=freqz(b1,1,M);
end
h2=freqz(b2,1,M);
end
% h=freqz(b,1,M);
for n=1:N+1
f=0:0.5/M:0.5-0.5/M;
if (n-1-N/2)==0;
hd=2*pi*f;
b2(n)=0;
plot(f,abs(h1),f,abs(h2),f,hd,'k-')
else
b2(n)=win(n)*(-1)^(n-1-N/2)/(n-1-N/2);
end
end
54
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是各频段上的
加权向量。
55
F、A及W的指定方式和例7.4.1和7.4.2所讨论过
的一样，唯一的差别是F的范围为0～1，而非
0～0.5, 1对应抽样频率的一半。需要指出的是，
若b的长度为偶数，设计高通和带阻滤波器时
有可能出现错误，因此，最好保证b的长度为
奇数，也即N应为偶数。
56
例7.4.1： 设计低通 FIR DF:  p  0.6 ,  s  0.7
调整通带、阻带的加权及滤波器的长度。
b=remez(N, F, A, W)
clear all;
0.7
f=[0 .6 .7 1];% 给定频率轴分点；
0.6
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);
调整N或W的结果
figure(1);stem(b,'.');grid;
figure(2);plot(w,h);grid;
57
1
例7.4.2： 设计多阻带滤波器，抽样频率500Hz,
在50Hz、 100Hz 及150Hz处陷波。
通带加
权为8，
阻带为1
-17dB
通带、
阻带加
权都是1
-25dB
58
例7.4.2： 设计多阻带滤波器，抽样频率500Hz,
在50Hz、 100Hz 及150Hz处陷波。
0.14 0.26 0.34 0.46 0.54 0.66
0.18 0.22 0.38 0.42 0.58 0.62
250Hz
1
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')
59
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的长度为奇数。
60
与本章有关的 MATLAB 文件
fftfilt.m 用叠接相加法实现长序列卷积。格式是:
y=fftfilt(h,x) 或
y=fftfilt(h, x，N)
记
的长度为
，
的长度为
。 若采
用第一个调用方式，程序自动地确定对
分段
的长度 及做FFT的长度 ， 显然， 是最接近
的2的整次幂。分的段数为
。
采用第二个调用方式，使用者可自己指定做FFT
的长度。建议使用第一个调用方式。
61
例3.9.1令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
N=500;p=0.05;f=1/16;
%
of a 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);
% y=x*h
subplot(211)
plot(x);
subplot(212)
plot(y);
62
2
1
0
-1
-2
0
50
100
150
200
250
300
350
400
450
500
0
50
100
150
200
250
300
350
400
450
500
2
1
0
-1
-2
63
clear;
% 生成滤波器系数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)
64
subplot(212);plot(y(1:N))
2
1
0
-1
-2
0
50
100
150
200
250
300
350
400
450
500
0
50
100
150
200
250
300
350
400
450
500
2
1
0
-1
-2
65
hilbert.m
文件用来计算信号Hilbert变换。调用
的格式是: y=hilbert(x),y的实部就
是
，虚部是的Hilbert变换
所以，y 实际上是 x 的解析信号。
66
。
•
czt.m 调用格式是： X＝czt(x, M, W, A) 。x
是待变换的时域信号，其长度设为N，M是
变换的长度，W确定变换的步长，A确定变
换的起点。若M=N, A=1, 则CZT变成DFT。
j0
Ae , W e
 j0
A=exp(j*2*pi*f0/fs);
W=exp(-j*2*pi*DELf/fs);
67
设x(n)是两个正弦信号及白噪声的叠加，进行
频谱分析。
•例4.7.2
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(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)));
68
10
0
-10
0
10
20
30
40
50
60
70
600
400
200
0
-0.5
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5
10
0
-10
69
0
10
20
30
40
50
60
70
例4.7.2 设x(n)由三个实正弦组成，频率分别
是8HZ， 8.22HZ 和9HZ， 抽样频率是40HZ ，时域
取128点。
用CZT计算
的DFT
100
50
0
用FFT计算
的DFT
4
6
8
10
12
14
16
18
20
0
2
4
6
8
10
12
14
16
18
20
50
100
50
0
70
2
100
0
7~(7+M×0.05) HZ
的CZT变换
0
7
7.5
8
8.5
9
9.5
10
10.5
•
程序
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;
71
% DTFT
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;
72
与本章内容有关的MATLAB文件
1．filtfilt.m 本文件实现零相位滤波。其调用格
式是：y=filtfilt(B, A, x) 。式中B是 H ( z ) 的分子
多项式，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');
73
2．grpdelay.m 求系统的群延迟。调用格式
[gd w]=grpdelay(B, A, N) , 或
[gd F]=grpdelay(B, A, N, FS)
式中B和A仍是 H ( z ) 的分子、分母多项式，gd
是群延迟，w、F是频率分点，二者的维数均
为N；FS为抽样频率，单位为Hz。
3
Amplitude Freq. Res.
2.5
2
1.5
1
0.5
0
74
0
0.5
1
1.5
2
2.5
3
3.5
3．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=(x^3+1)/(x+1)
y=[1 0 0 1]; x=[1 1];
[q,r]=deconv(y,x);
q
75
3．tf2latc.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。
76
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是
 (n) 则输出y是 H ( z ) 的系数；
 ( N 1)
1
z
H
(
z
) 的系数。
g是
77
5. tf2ss.m 和 ss2tf.m 实现转移函数和相应状
态变量之间的转换。二者的调用格式分别是：
[A, B, C, D]=tf2ss(b, a), [b, a]=ss2tf(A, B, C, D)。
H ( z)
式中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的说明。
78