第4章 MATLAB在信号处理中的应用 4.1 信号及其表示4.2 信号的基本运算 4.3 信号的能量和功率 4.4 线性
时不变系统4.5 线性时不变系统的响应4.6 线性时不变系统的频率响应 4.7 傅里叶(Fourier)变换4.
8 IIR数字滤波器的设计方法4.9 FIR数字滤波器设计4.1 信号及其表示4.1.1连续时间信号的表示 连续
时间信号：时间变化连续。如y=x(t) 离散时间信号(序列)：时间离散，如x(nT)=x(t)|t=nT.4.1.2工具箱中的
信号产生函数4.1.3离散时间信号的表示 在MATLAB中，离散时间信号x(n)的表示：需用一个向量x表示序列幅值，
用另一个等长的定位时间变量n，才能完整地表示一个序列。 [例4-10] 绘制离散时间信号的棒状图。其中x(-1)=-1, x(0)
=1, x(1)=2, x(2)=1, x(3)=0, x(4)=-1。MATLAB源程序为：n=-3:5; %定位时间
变量x=[0,0,-1,1,2,1,-1,0,0];stem(n,x); grid; % 绘制棒状图line([-3,5],[0
,0]); %画x轴线xlabel(''n''); ylabel(''x[n]'')运行结果如图4.11所示。图 4.11 离散时间
信号图形4.1.4几种常用离散时间信号的表示1．单位脉冲序列直接实现：x=zeros(1,N); x(1,n0)=1;2．单位阶跃
序列 直接实现：n=[ns:nf]; x=[(n-n0)>=0];3．实指数序列直接实现：n=[ns:nf]; x=a.^n;4．
复指数序列直接实现：n=[ns:nf]; x=exp((sigema+jw)n);5．正(余)弦序列直接实现：n=[ns:nf]
; x=cos(wn+sita);4.2 信号的基本运算4.2.1信号的相加与相乘 y(n)=x
1(n)+x2(n) y(n)=x1(n)×x2(n) MATLAB实现：y=x1+x2; y=x1.x24
.2.2序列移位与周期延拓运算序列移位：y(n)=x(n-m)。MATLAB实现：y=x; ny=nx-m序列周期延拓：y(n)=
x((n))M，MATLAB实现：ny=nxs:nxf；y=x(mod(ny,M)+1)4.2.3 序列翻褶与序列累加运算序列翻褶
：y(n)=x(-n)。MATLAB可实现: y=fliplr(x)序列累加的数学描述为： MATLAB实现：y=cumsum(x
)4.2.4 两序列的卷积运算两序列卷积运算： MATLAB实现：y=conv(x1,x2)。序列x1(n)和x2(n)必须长度有
限。 4.2.5 两序列的相关运算两序列相关运算： 。MATLAB实现：y=xcorr(x1,x2)。4.3 信号的能量和功率1.
信号能量数字定义：MATLAB实现: E=sum(x.conj(x)); 或 E=sum(abs(x).^2);数字定义：2.
信号功率MATLAB实现: P=sum(x.conj(x))/N; 或 E=sum(abs(x).^2)/N;4.4 线性时不
变系统4.4.1 系统的描述1．常系数线性微分/差分方程2．系统传递函数3．零－极点增益模型连续系统： 连续系统： 离散系统： 离
散系统： 4．极点留数模型离散系统： 连续系统： 5．二次分式模型连续系统： 离散系统： 6．状态空间模型连续系统： 离散系统：
4.4.2 系统模型的转换函数 在MATLAB中，用sos、ss、tf、zp分别表示二次分式模型、状态空间模型、传递
函数模型和零－极点增益模型。其中sos表示二次分式，g为比例系数，sos为L×6的矩阵，即 （4－15） 1．ss2tf
函数 格式：[num, den]=ss2tf(A,B,C,D,iu)功能：将指定输入量iu的线性系统(A,B,C,D)转换为传递函
数模型[num,den]。2．zp2tf函数格式：[num,den]=zp2tf(z,p,k)功能：将给定系统的零－极点增益模型转
换为传递函数模型，z、p、k分别为零点列向量、极点列向量和增益系数。 线性系统模型的变换函数[例4-18] 求离散时间系统的零、
极点向量和增益系数。在命令窗口输入：>> num=[2,3]; den=[1,0.4,1];>> [num,den]=eqtfle
ngth(num,den); %使长度相等>> [z,p,k]=tf2zp(num,den)屏幕显示为z = 0 -1.50
00p = -0.2000 + 0.9798i -0.2000 - 0.9798ik = 24.4.3 系统互联与系统结
构MATLAB实现函数series( ) 格式：[A,B,C,D]=series(A1,B1,C1,D1,A2,B2,C2,D2)
或 [num,den]=series(num1,den1,num2,den2)将系统1、系统2级联，可得到级联连
接的传递函数形式为：1. 系统的级联MATLAB实现函数parallel( )格式：[A,B,C,D]=parallel(A1,B
1,C1,D1,A2,B2,C2,D2) 或 [num,den]=parallel(num1,den1,num2,
den2)2. 系统的并联将系统1、系统2并联，可得到并联连接的传递函数形式为：3. 两个系统的反馈连接函数feedback格式：
[A,B,C,D]=feedback(A1,B1,C1,D1,A2,B2,C2,D2,sign) 或 [nu
m,den]=feedback(num1,den1,num2,den2,sign)将系统1和系统2进行反馈连接，sign表示反馈方
式(默认值为-1)； 当sig=+1时表示正反馈；当sig=-1时表示负反馈。[例4-19] 求两个单输入单
输出子系统的级联、并联和反馈后系统的传递函数。MATLAB源程序为： num1=1; den1=[1,1]; %系统1 num2=
2; den2=[1,2]; %系统2 [nums,dens]=series(num1,den1,num2,den2) %实现
两个系统级联 [nump,denp]=parallel(num1,den1,num2,den2) %实现两个系统并联
[numf,denf]=feedback(num1,den1,num2,den2) %实现两个系统反
馈程序运行结果为：nums = 0 0 2 ； dens = 1 3 2nump =
0 3 4 ； denp = 1 3 2numf = 0 1 2
； denf = 1 3 4因此，各系统的传递函数分别为：4.5 线性时不变系统的响应4.5.1 线性
时不变系统的时域响应1．连续LTI系统的响应2．离散LTI系统的响应用MATLAB中的卷积函数conv( )来实现。 用MATLA
B中的卷积函数conv( )来实现。 格式：[y,x]=lsim(a,b,c,d,u,t)功能：返回连续LTI系统 (2)对任意输
入的离散LTI系统响应函数dlsim( )格式：[y,x]=dlsim(a,b,c,d,u)功能：返回离散LTI系统 对任意输入时
系统的输出响应y和状态记录x，其中u给出每个输入的时序列，一般情况下u为一个矩阵；t用于指定仿真的时间轴，它应为等间隔。对输入序列
u的响应y和状态记录x。3．时域响应函数(1)对任意输入的连续LTI系统响应函数lsim( ) 4.5.2 LTI系统的单位冲激响
应1. 求连续LTI系统的单位冲激响应函数impulse( )格式：[Y,T] = impulse(sys) 或impulse(s
ys)功能：返回系统的响应Y和时间向量T，自动选择仿真的时间范围。其中sys可为系统传递函数、零极增益模型或状态空间模型。2. 求
离散系统的单位冲激响应函数dimpulse( )格式：[y,x]=dimpulse(num,den)功能：返回项式传递函数的单位冲
激响应y向量和时间状态历史记录x向量。4.5.3 时域响应的其它函数1. 求连续LTI系统的零输入响应函数initial( )格式
：[y,t,x]=initial(a,b,c,d,x0)功能：计算出连续时间LTI系统由于初始状态x0所引起的零输入响应y。其中x
为状态记录，t为仿真所用的采样时间向量。2. 求离散系统的零输入响应函数dinitial( )格式：[y,x,n]=dinitia
l(a,b,c,d,x0)功能：计算离散时间LTI系统由初始状态x0所引起的零输入响应y和状态响应响应x，取样点数由函数自动选取。
n为仿真所用的点数。3. 求连续系统的单位阶跃响应函数step( )格式：[Y,T] = step(sys) 功能：返回系统的单位
阶跃响应Y和仿真所用的时间向量T，自动选择仿真的时间范围。其中sys可为系统传递函数(TF)、零极增益模型(ZPK)或状态空间模型
(SS)。4. 求离散系统的单位阶跃响应函数dstep( )格式：[y,x]= dstep (num,den)功能：返回多项式传递
函数G(z)=num(z)/den(z)表示的系统单位阶跃响应。4.6线性时不变系统的频率响应1．求模拟滤波器Ha(s)的频率响应
函数freqs( )格式：H＝freqs(B,A,W) 功能：计算由向量W(rad/s)指定的频率点上模拟滤器系统函数Ha(s)的
频率响应Ha(jΩ)，结果存于H向量中。[例4-31] 已知某模拟滤波器的系统函数求该模拟滤波器的频率响应。MATLAB源程序如下
。B=1;A=[1 2.6131 3.4142 2.6131 1];W=0:0.1:2pi5;freqs(B,A,W)图4.3
0 模拟滤波器的频率响应 [例4-32] 已知某滤波器的系统函数为求该滤波器的频率响应。MATLAB源程序为：B=[1
0 0 0 0 0 0 0 –1];A=1;freqz(B,A)该程序运行所绘出的幅频与相频性曲线如图4.31所示。图4.31滤波
器幅度和相位曲线 2．求数字滤波器H(z)的频率响应函数freqz( )格式：H=freqz(B,A,W)功能：计算由向量W(ra
d)指定的数字频率点上(通常指在H(z)的频率响应H(ejw )。 3．滤波函数filter格式：y=filter(B,A,x)功
能：对向量x中的数据进行滤波处理，即差分方程求解，产生输出序列向量y。B和A分别为数字滤波器系统函数H(z)的分子和分母多项式系数
向量。[例4-33] 设系统差分方程为MATLAB源程序为：B=1; A=[1,-0.8];N=0:31; x=0.8.^n;y=
filter(B,A,x);subplot(2,1,1);stem(x)subplot(2,1,2);stem(y)该程序运行所得
结果如图4.32所示。，求该系统对信号的响应。图4.32系统对信号的响应 4.7傅里叶(Fourier)变换4.7.1连续时间、连
续频率－傅里叶变换4.7.2 连续时间、离散频率－傅里叶级数正变换： 逆变换： 正变换： 逆变换： 4.7.3 时间离散、连续频率
－序列傅里叶变换4.7.4 离散时间、离散频率－离散傅里叶级数4.7.5离散时间、离散频率－离散傅里叶变换(DFT)正变换： 逆变
换： 正变换： 逆变换： 正变换： 逆变换： 1．一维快速正傅里叶变换函数fft格式：X=fft(x, N)功能：采用FFT算法计
算序列向量x的N点DFT变换， 当N缺省时，fft函数自动按x的长度计算DFT。当N为2整数次幂时，fft按基-2算法计算，否则用
混合算法。2．一维快速逆傅里叶变换函数ifft格式：x=ifft(X, N)功能：采用FFT算法计算序列向量X的N点IDFT变换。
[例4-36] 用快速傅里叶变换FFT计算下面两个序列的卷积。, 并测试直接卷积和快速卷积的时间。图4.35 快速卷积框图MA
TLAB程序(部分)：%线性卷积xn= sin(0.4[1:15]); %对序列x(n)赋值, M=15hn= 0.9.^(1:
20); %对序列h(n)赋值, N=20yn=conv(xn,hn); % 直接调用函数conv计算卷积%园周卷
积L=pow2(nextpow2(M+N-1)); Xk=fft(xn,L); Hk=fft(hn,L); Yk=Xk.
Hk; yn=ifft(Yk,L); 图4.36 x(n),h(n)及其线性卷
积波形4.8 IIR数字滤波器的设计方法1. 数字滤波器的频率响应函数幅度响应：相位响应：图4.37 理想低通、高通、带通、带阻
数字滤波器幅度特性 2. 滤波器的技术指标 幅度响应指标、相位响应指标 图4.38 数字低通滤波器的幅度特性
通带要求： 阻带要求： 通带最大衰减：阻带最小衰减：4.8.1冲激响应不变法2.MATLAB信号处理工箱中的专用函数impinv
ar( )：格式：[BZ,AZ] =impinvar(B,A,Fs) 功能：把具有[B,A]模拟滤波器传递函数模型转换成采样频率为
Fs(Hz)的数字滤波器的传递函数模型[BZ,AZ]。采样频率Fs的默认值为Fs=1。1. 冲激响应不变法设计IIR数字滤波器的基
本原理： [例4-37] MATLAB源程序如下：num=[1]; %模拟滤波器系统函数的分子den=[1,sqrt(5),
2,sqrt(2),1]; %模拟滤波器系统函数的分母[num1,den1]=impinvar(num,den) %求数字低
通滤波器的系统函数程序的执行结果如下：num1 = -0.0000 0.0942 0.2158 0.0311
den1 = 1.0000 -2.0032 1.9982 -0.7612 0.1069MATLAB信号处
理工具箱中的专用双线性变换函数bilinear( )格式：[numd,dend]＝bilinear(num,den,Fs)功能：把
模拟滤波器的传递函数模型转换成数字滤波器的传递函数模型。4.8.2双线性变换法双线性变换利用频率变换关系： [例4-38] MAT
LAB源程序如下： num=1; %模拟滤波器系统函数的分子 den=[1,sqrt(3),sqrt(2),
1]; %模拟滤波器系统函数的分母 [num1,den1]=bilinear(num,den,1) %求数字滤波器的
传递函数运算的结果如下： num1 = 0.0533 0.1599 0.1599 0.0533den1 =
1.0000 -1.3382 0.9193 -0.15464.8.3 IIR数字滤波器的频率变换设计法1.
IIR数字滤波器的频率变换设计法的基本原理 根据滤波器设计要求，设计模拟原型低通滤波器，然后进行频率变换，将其转换为相应的模
拟滤波器(高通、带通等)，最后利用冲激响应不变法或双线性变换法，将模拟滤波器数字化成相应的数字滤波器。 图4.39 IIR数字滤
波器MATLAB设计步骤流程图 1．MATLAB的典型设计利用在MATLAB设计IIR数字滤波器可分以下几步来实现 (1)按一定规
则将数字滤波器的技术指标转换为模拟低通滤波器的技术指标；(2)根据转换后的技术指标使用滤波器阶数函数，确定滤波器的最小阶数N和截止
频率Wc；(3)利用最小阶数N产生模拟低通滤波原型；(4)利用截止频率Wc把模拟低通滤波器原型转换成模拟低通、高通、带通或带阻滤波
器；(5)利用冲激响应不变法或双线性不变法把模拟滤波器转换成数字滤波器。[例4-39] 设计一个数字信号处理系统，它的采样率为Fs
＝100Hz，希望在该系统中设计一个Butterworth型高通数字滤波器，使其通带中允许的最小衰减为0.5dB，阻带内的最小衰减
为40dB，通带上限临界频率为30Hz，阻带下限临界频率为40Hz。 MATLAB源程序设计如下： %把数字滤波器的
频率特征转换成模拟滤波器的频率特征 wp=302pi;ws=402pi;rp=0.5;rs=40;Fs=10
0; [N,Wc]=buttord(wp,ws,rp,rs,''s''); %选择滤波器的最小阶数
[Z,P,K]=buttap(N); %创建Butterworth低通滤波器原型 [A
,B,C,D]=zp2ss(Z,P,K); %零－极点增益模型转换为状态空间模型 [AT
,BT,CT,DT]=lp2hp(A,B,C,D,Wc); %实现低通向高通的转变 [num1,den1]=s
s2tf(AT,BT,CT,DT); %状态空间模型转换为传递函数模型 %运用双线性变换法把模拟滤波器转换
成数字滤波器 [num2,den2]=bilinear(num1,den1,100);
[H,W]=freqz(num2,den2); %求频率响应 plot(WFs/(2
pi),abs(H));grid; %绘出频率响应曲线 xlabel(''频率/Hz''); yla
bel(''幅值'')程序运行结果如图4.40所示。 2．MATLAB的直接设计图4.39 IIR数字滤波器MATLAB设计步骤流程
图 [例4-41] 试设计一个带阻IIR数字滤波器，其具体的要求是：通带的截止频率：wp1＝650Hz、wp2＝850Hz；阻带的
截止频率：ws1＝700Hz、ws2＝800Hz；通带内的最大衰减为rp＝0.1dB；阻带内的最小衰减为rs＝50dB；采样频率为
Fs＝2000Hz。MATLAB源程序设计如下： wp1=650;wp2=850;ws1=700;ws2=800;rp=0.1;r
s=50;Fs=2000; wp=[wp1,wp2]/(Fs/2);ws=[ws1,ws2]/(Fs/2); %利用Nyquis
t频率频率归一化 [N,wc]=ellipord(wp,ws,rp,rs,''z'');
%求滤波器阶数 [num,den]=ellip(N,rp,rs,wc,''stop''); %求滤波器传递函数 [H,W]=
freqz(num,den); %绘出频率响应曲线 plot(WFs/(2pi),abs(H));grid;
xlabel(''频率/Hz'');ylabel(''幅值'')该程序运行后的幅频响应曲线如图4.42所示。 4.9 FIR数字滤波器设
计格式：w = boxcar(M) 功能：返回M点矩形窗序列。MATLAB信号处理工具箱中的窗函数法设计FIR数字滤波器的专用命令
fir1( )。格式：B＝fir1(N,wc)功能：设计一个具有线性相位的N阶(N点)的低通FIR数字滤波器，返回的向量B为滤波器
的系数(单位冲激响应序列)，其长度为N+1。4.9.1窗函数设计法窗函数设计的基本原理：
h(n)=w(n)hd(n) w(n)为窗函数， hd(n)理想数字滤波器的单位冲激响应。
在MATLAB信号处理工具箱中为用户提供了Boxcar (矩形)、Bartlet(巴特利特)、Hanning(汉宁)等窗函数，这
些窗函数的调用格式相同。 FIR数字滤波器的单位冲激响应h(n)满足偶(奇)对称 h(n)=
h(N-n-1) 或 h(n)=-h(N-n-1)FIR数字滤波器具有线性相位：或[例4-43] 用矩形窗设计线性相位FI
R低通滤波器。该滤波器的通带截止频率wc=pi/4，单位脉冲响h(n)的长度M=21。并绘出h(n)及其幅度响应特性曲线。MATL
AB源程序为：M=21; wc=pi/4; % 理想低通滤波器参数n=0:M-1; r=(M-1)/2;nr
=n-r+eps((n-r)==0);hdn=sin(wcnr)/pi./nr; % 计算理想低通单位脉冲响应hd
(n)if rem(M,2)~=0, hdn(r+1)=wc/pi; end % M为奇数时，处理n=r点的0/0型wn1=bo
xcar(M); % 矩形窗hn1=hdn.wn1''; % 加窗subplot(2,1,1);stem(
n,hn1,''.''); line([0,20],[0,0]);xlabel(''n''),ylabel(''h(n)''),title(''
矩形窗设计的h(n)'');hw1=fft(hn1,512); w1=2[0:511]/512; %求频谱subplot(2,
1,2), plot(w1,20log10(abs(hw1)))xlabel(''w/pi''), ylabel(''幅度(dB)'')
; title(''幅度特性(dB)'');程序运行结果如图4.44所示。4.9.2频率抽样法1. 频率抽样法的基本原理
对所期望的滤波器的频率响应Hd(ejw)进行等间隔采样获得H(k)，利用h(n)=IDFT[H(k)]求得FIR的单位冲激响应。
2. MATLAB信号处理工具箱中的频率抽样法专用函数命令fir2( )格式：B＝fir2(N,F,A)功能：设计一个N阶的FIR数字滤波器，其频率响应由向量F和A指定，滤波器的系数（单位冲激响应）返回在向量B中，长度为N+1。向量F和A分别指定滤波器的采样点的频率及其幅值，F中的频率必须在0.0到1.0之间，1.0对应于采样频率的一半。它们必须按递增的顺序从0.0开始到1.0为结束。 [例4-47] 试用频率抽样法设计一个FIR低通滤波器，该滤波器的截止频率为0.5pi，频率抽样点数为33。MATLAB源程序为：N=32;F=[0:1/32:1]; %设置抽样点的频率，抽样频率必须含0和1。A=[ones(1,16),zeros(1,N-15)]; %设置抽样点相应的幅值B=fir2(N,F,A);freqz(B); %绘制滤波器的幅相频曲线figure(2);stem(B,''.''); %绘制单位冲激响应的实部line([0,35],[0,0]);xlabel(''n'');ylabel(''h(n)'');图4.49滤波器的频率响应和单位冲激响应序列 4.9.3 MATLAB的其它相关函数1．最小二乘逼近法设计线性相位FIR滤波器函数fircls( )2．有限制条件的最小二乘逼近法设计低通和高通FIR数字滤波器函数fircls1( )3．最小二乘逼近法设计线性相位FIR数字滤波器函数firls( )4．升余弦FIR滤波器设计函数firrcos( )5．Parks-McClellan优化等纹波FIR滤波器设计函数remez( ) |
|
