第4章 MATLAB在信号处理中的应用 4.1 信号及其表示 4.2 信号的基本运算 4.3 信号的能量和功率 4.4 线性时不变系统 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) 4.1.2工具箱中的信号产生函数 4.1 信号及其表示 4.1.1连续时间信号的表示 连续时间信号：时间变化连续。如y=x(t) 离散时间信号(序列)：时间离散，如x(nT)=x(t)|t=nT. 4.1.2工具箱中的信号产生函数 函数名 功能 sawtooth 产生锯齿波或三角波信号 pulstran 产生冲激串 square 产生方波信号 rectpule 产生非周期的方波信号 sinc 产生sinc函数波形 tripuls 产生非周期的三角波信号 chirp 产生调频余弦信号 diric 产生Dirichlet或周期sinc函数 gauspuls 产生高斯正弦脉冲信号 gmonopuls 产生高斯单脉冲信号 vco 电压控制振荡器

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(w*n+sita);

4.2 信号的基本运算 4.2.1信号的相加与相乘 4.2.2序列移位与周期延拓运算 4.2.3 序列翻褶与序列累加运算 y(n)=x1(n)+x2(n) y(n)=x1(n)×x2(n) MATLAB实现：y=x1+x2; y=x1.*x2 4.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 两序列的卷积运算 4.2.5 两序列的相关运算 两序列卷积运算： MATLAB实现：y=conv(x1,x2)。序列x1(n)和x2(n)必须长度有限。 4.2.5 两序列的相关运算 两序列相关运算： 。MATLAB实现：y=xcorr(x1,x2)。

4.3 信号的能量和功率 1.信号能量 2. 信号功率 数字定义： 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.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分别为零点列向量、极点列向量和增益系数。

线性系统模型的变换函数 函数名 功能说明 ss2tf 状态空间模型转换为传递函数模型 zp2tf 零－极点增益模型转换为传递函数模型 ss2zp 状态空间模型转换为零－极点增益模型 zp2ss 零－极点增益模型转换为状态空间模型 ss2sos 状态空间模型转换为二次分式模型 zp2sos 零－极点增益模型转换为二次分式模型 tf2ss 传递函数模型转换为状态空间模型 sos2tf 二次分式模型转换为传递函数模型 tf2zp 传递函数模型转换为零－极点增益模型 sos2zp 二次分式模型转换为零－极点增益模型 tf2sos 传递函数模型转换为二次分式模型 sos2ss 二次分式模型转换为状态空间模型

[例4-18] 求离散时间系统 的零、极点向量和增益系数。 在命令窗口输入： >> num=[2,3]; den=[1,0.4,1]; >> [num,den]=eqtflength(num,den);　%使长度相等 >> [z,p,k]=tf2zp(num,den) 屏幕显示为 z =　 0 -1.5000 p = -0.2000 + 0.9798i 　 -0.2000 - 0.9798i k = 2

[num,den]=series(num1,den1,num2,den2) 4.4.3 系统互联与系统结构 1. 系统的级联 MATLAB实现函数series( ) 格式：[A,B,C,D]=series(A1,B1,C1,D1,A2,B2,C2,D2) 或 [num,den]=series(num1,den1,num2,den2) 将系统1、系统2级联，可得到级联连接的传递函数形式为：

[num,den]=parallel(num1,den1,num2,den2) 2. 系统的并联 MATLAB实现函数parallel( ) 格式：[A,B,C,D]=parallel(A1,B1,C1,D1,A2,B2,C2,D2) 或 [num,den]=parallel(num1,den1,num2,den2) 将系统1、系统2并联，可得到并联连接的传递函数形式为： 3. 两个系统的反馈连接 函数feedback 格式：[A,B,C,D]=feedback(A1,B1,C1,D1,A2,B2,C2,D2,sign) 或 [num,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 2 nump = 0 3 4 ； denp = 1 3 2 numf = 0 1 2 ； 　denf = 1 3 4 因此，各系统的传递函数分别为：

4.5 线性时不变系统的响应 4.5.1 线性时不变系统的时域响应 1．连续LTI系统的响应 4.5.1　线性时不变系统的时域响应 1．连续LTI系统的响应 用MATLAB中的卷积函数conv( )来实现。 2．离散LTI系统的响应 用MATLAB中的卷积函数conv( )来实现。

3．时域响应函数 (1)对任意输入的连续LTI系统响应函数lsim( ) 格式：[y,x]=lsim(a,b,c,d,u,t) 功能：返回连续LTI系统 对任意输入时系统的输出响应y和状态记录x，其中u给出每个输入的时序列，一般 情况下u为一个矩阵；t用于指定仿真的时间轴，它应为等间隔。 (2)对任意输入的离散LTI系统响应函数dlsim( ) 格式：[y,x]=dlsim(a,b,c,d,u) 功能：返回离散LTI系统 对输入序列u的响应y和状态记录x。

4.5.2 LTI系统的单位冲激响应 1. 求连续LTI系统的单位冲激响应函数impulse( ) 格式：[Y,T] = impulse(sys)　或impulse(sys) 功能：返回系统的响应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]=dinitial(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:2*pi*5; freqs(B,A,W) 图4.30 模拟滤波器的频率响应

2．求数字滤波器H(z)的频率响应函数freqz( ) 格式：H=freqz(B,A,W) 功能：计算由向量W(rad)指定的数字频率点上(通常指在H(z)的频率响应H(ejw )。 [例4-32] 已知某滤波器的系统函数为 求该滤波器的频率响应。MATLAB源程序为： B=[1 0 0 0 0 0 0 0 –1]; A=1; freqz(B,A) 该程序运行所绘出的幅频与 相频性曲线如图4.31所示。 图4.31滤波器幅度和相位曲线

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.5离散时间、离散频率－离散傅里叶变换(DFT) 4.7.3 时间离散、连续频率－序列傅里叶变换 正变换： 逆变换： 4.7.4 离散时间、离散频率－离散傅里叶级数 正变换： 逆变换： 4.7.5离散时间、离散频率－离散傅里叶变换(DFT) 正变换： 逆变换：

功能：采用FFT算法计算序列向量X的N点IDFT变换。 格式：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 快速卷积框图

MATLAB程序(部分)： %线性卷积 xn= sin(0.4*[1:15]);　%对序列x(n)赋值, M=15 hn= 0.9.^(1:20);　　 %对序列h(n)赋值, N=20 yn=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冲激响应不变法 1. 冲激响应不变法设计IIR数字滤波器的基本原理： 2.MATLAB信号处理工箱中的专用函数impinvar( )： 格式：[BZ,AZ] =impinvar(B,A,Fs) 功能：把具有[B,A]模拟滤波器传递函数模型转换成采样频率为Fs(Hz)的数字滤波器的 传递函数模型[BZ,AZ]。采样频率Fs的默认值为Fs=1。 [例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.1069

4.8.2双线性变换法 双线性变换利用频率变换关系： MATLAB信号处理工具箱中的专用双线性变换函数bilinear( ) 格式：[numd,dend]＝bilinear(num,den,Fs) 功能：把模拟滤波器的传递函数模型转换成数字滤波器的传递函数模型。 [例4-38] MATLAB源程序如下： num=1; %模拟滤波器系统函数的分子 den=[1,sqrt(3),sqrt(2),1]; %模拟滤波器系统函数的分母 [num1,den1]=bilinear(num,den,1) %求数字滤波器的传递函数 运算的结果如下：　 num1 = 0.0533 0.1599 0.1599 0.0533 den1 = 1.0000 -1.3382 0.9193 -0.1546

4.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=30*2*pi;ws=40*2*pi;rp=0.5;rs=40;Fs=100; [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]=ss2tf(AT,BT,CT,DT); %状态空间模型转换为传递函数模型 %运用双线性变换法把模拟滤波器转换成数字滤波器 [num2,den2]=bilinear(num1,den1,100);　 [H,W]=freqz(num2,den2);　 %求频率响应 plot(W*Fs/(2*pi),abs(H));grid; %绘出频率响应曲线 xlabel('频率/Hz'); ylabel('幅值') 程序运行结果如图4.40所示。

wp1=650;wp2=850;ws1=700;ws2=800;rp=0.1;rs=50;Fs=2000; 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;rs=50;Fs=2000; wp=[wp1,wp2]/(Fs/2);ws=[ws1,ws2]/(Fs/2);　 %利用Nyquist频率频率归一化 [N,wc]=ellipord(wp,ws,rp,rs,'z'); %求滤波器阶数 [num,den]=ellip(N,rp,rs,wc,'stop');　　　　　%求滤波器传递函数 [H,W]=freqz(num,den); 　　　　　　　　　%绘出频率响应曲线 plot(W*Fs/(2*pi),abs(H));grid; xlabel('频率/Hz');ylabel('幅值') 该程序运行后的幅频响应曲线如图4.42所示。

4.9 FIR数字滤波器设计 4.9.1窗函数设计法 FIR数字滤波器的单位冲激响应h(n)满足偶(奇)对称 h(n)=h(N-n-1) 或 h(n)=-h(N-n-1) FIR数字滤波器具有线性相位： 或 4.9.1窗函数设计法 窗函数设计的基本原理： h(n)=w(n)hd(n) w(n)为窗函数， hd(n)理想数字滤波器的单位冲激响应。 在MATLAB信号处理工具箱中为用户提供了Boxcar (矩形)、Bartlet(巴特利特)、Hanning(汉宁)等窗函数，这些窗函数的调用格式相同。 格式：w = boxcar(M) 功能：返回M点矩形窗序列。 MATLAB信号处理工具箱中的窗函数法设计FIR数字滤波器的专用命令fir1( )。 格式：B＝fir1(N,wc) 功能：设计一个具有线性相位的N阶(N点)的低通FIR数字滤波器，返回的向量B 为滤波器的系数(单位冲激响应序列)，其长度为N+1。

[例4-43] 用矩形窗设计线性相位FIR低通滤波器。该滤波器的通带截止频率wc=pi/4，单位脉冲响h(n)的长度M=21。并绘出h(n)及其幅度响应特性曲线。MATLAB源程序为： M=21; wc=pi/4; 　　　　　　　 % 理想低通滤波器参数 n=0:M-1; r=(M-1)/2; nr=n-r+eps*((n-r)==0); hdn=sin(wc*nr)/pi./nr; 　 % 计算理想低通单位脉冲响应hd(n) if rem(M,2)~=0, hdn(r+1)=wc/pi; end　 % M为奇数时，处理n=r点的0/0型 wn1=boxcar(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,20*log10(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( ) 3．最小二乘逼近法设计线性相位FIR数字滤波器函数firls( ) 4．升余弦FIR滤波器设计函数firrcos( ) 5．Parks-McClellan优化等纹波FIR滤波器设计函数remez( )