MATLAB在信号处理中的应用.ppt

第 4 章 MATLAB在信号处理中的应用 第4章 MATLAB在信号处理中的应用 4.1 信号及其表示 4.2 信号的基本运算 4.3 信号的能量和功率 4.4 线性时不变系统 4.5 线性时不变系统的响应 4.6 线性时不变系统的频率响应 4.7 傅里叶(Fourier)变换 4.8 IIR数字滤波器的设计方法 4.9 FIR数字滤波器设计 第 4 章 MATLAB在信号处理中的应用 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 章 MATLAB在信号处理中的应用 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(xn) 运行结果如图4.11所示。 图 4.11 离散时间信号图形 第 4 章 MATLAB在信号处理中的应用 4.1.4几种常用离散时间信号的表示 1单位脉冲序列 直接实现：x=zeros(1,N); x(1,n0)=1; 2单位阶跃序列 直接实现：n=ns:nf; x=(n-n0)=0; 第 4 章 MATLAB在信号处理中的应用 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 章 MATLAB在信号处理中的应用 4.2 信号的基本运算 4.2.1信号的相加与相乘 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 章 MATLAB在信号处理中的应用 4.2.4 两序列的卷积运算 两序列卷积运算： MATLAB实现：y=conv(x1,x2)。序列x1(n)和x2(n)必须长度有限。 4.2.5 两序列的相关运算 两序列相关运算： 。MATLAB实现：y=xcorr(x1,x2)。 第 4 章 MATLAB在信号处理中的应用 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 章 MATLAB在信号处理中的应用 4.4 线性时不变系统 4.4.1 系统的描述 1常系数线性微分/差分方程 2系统传递函数 3零极点增益模型 连续系统： 连续系统： 离散系统： 离散系统： 第 4 章 MATLAB在信号处理中的应用 4极点留数模型 离散系统： 连续系统： 5二次分式模型 连续系统： 离散系统： 6状态空间模型 连续系统： 离散系统： 第 4 章 MATLAB在信号处理中的应用 4.4.2 系统模型的转换函数 在MATLAB中，用sos、ss、tf、zp分别表示二次分式模型、状态空间 模型、传递函数模型和零极点增益模型。其中sos表示二次分式，g为比 例系数，sos为L6的矩阵，即 （415） 1ss2tf函数 格式：num, den=ss2tf(A,B,C,D,iu) 功能：将指定输入量iu的线性系统(A,B,C,D)转换为传递函数模型num,den。 2zp2tf函数 格式：num,den=zp2tf(z,p,k) 功能：将给定系统的零极点增益模型转换为传递函数模型，z、p、k分别为零 点列向量、极点列向量和增益系数。 第 4 章 MATLAB在信号处理中的应用 线性系统模型的变换函数 函数名功能说明函数名功能说明 ss2tf 状态态空间间模型转换为转换为 传递传递 函数模型 zp2tf零极点增益模型转换为传递转换为传递 函数模型 ss2zp状态态空间间模型转换为转换为 零 极点增益模型 zp2ss零极点增益模型转换为转换为 状态态 空间间模型 ss2sos状态态空间间模型转换为转换为 二 次分式模型 zp2sos零极点增益模型转换为转换为 二次 分式模型 tf2ss传递传递 函数模型转换为转换为 状 态态空间间模型 sos2tf二次分式模型转换为传递转换为传递 函数 模型 tf2zp传递传递 函数模型转换为转换为 零 极点增益模型 sos2zp二次分式模型转换为转换为 零极点 增益模型 tf2sos传递传递 函数模型转换为转换为 二 次分式模型 sos2ss二次分式模型转换为转换为 状态态空间间 模型 第 4 章 MATLAB在信号处理中的应用 例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 第 4 章 MATLAB在信号处理中的应用 4.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. 系统的级联 第 4 章 MATLAB在信号处理中的应用 MATLAB实现函数parallel( ) 格式：A,B,C,D=parallel(A1,B1,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) 或 num,den=feedback(num1,den1,num2,den2,sign) 将系统1和系统2进行反馈连接，sign表示反馈方式(默认值为-1)； 当sig=+1时表示正反馈；当sig=-1时表示负反馈。 第 4 章 MATLAB在信号处理中的应用 例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 章 MATLAB在信号处理中的应用 4.5 线性时不变系统的响应 4.5.1 线性时不变系统的时域响应 1连续LTI系统的响应 2离散LTI系统的响应 用MATLAB中的卷积函数conv( )来实现。 用MATLAB中的卷积函数conv( )来实现。 第 4 章 MATLAB在信号处理中的应用 格式：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 章 MATLAB在信号处理中的应用 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 章 MATLAB在信号处理中的应用 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 章 MATLAB在信号处理中的应用 4.6线性时不变系统的频率响应 1求模拟滤波器Ha(s)的频率响应函数freqs( ) 格式：Hfreqs(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 模拟滤波器的频率响应 第 4 章 MATLAB在信号处理中的应用 例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(rad)指定的数字频率点上(通常指在H(z)的频率响应H(ejw )。 第 4 章 MATLAB在信号处理中的应用 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 章 MATLAB在信号处理中的应用 4.7傅里叶(Fourier)变换 4.7.1连续时间、连续频率傅里叶变换 4.7.2 连续时间、离散频率傅里叶级数 正变换： 逆变换： 正变换： 逆变换： 第 4 章 MATLAB在信号处理中的应用 4.7.3 时间离散、连续频率序列傅里叶变换 4.7.4 离散时间、离散频率离散傅里叶级数 4.7.5离散时间、离散频率离散傅里叶变换(DFT) 正变换： 逆变换： 正变换： 逆变换： 正变换： 逆变换： 第 4 章 MATLAB在信号处理中的应用 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 快速卷积框图 第 4 章 MATLAB在信号处理中的应用 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 章 MATLAB在信号处理中的应用 4.8 IIR数字滤波器的设计方法 1. 数字滤波器的频率响应函数 幅度响应： 相位响应： 图4.37 理想低通、高通、带通、带阻数字滤波器幅度特性 第 4 章 MATLAB在信号处理中的应用 2. 滤波器的技术指标 幅度响应指标、相位响应指标 图4.38 数字低通滤波器的幅度特性 通带要求： 阻带要求： 通带最大衰减： 阻带最小衰减： 第 4 章 MATLAB在信号处理中的应用 4.8.1冲激响应不变法 2.MATLAB信号处理工箱中的专用函数impinvar( )： 格式：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.1069 第 4 章 MATLAB在信号处理中的应用 MATLAB信号处理工具箱中的专用双线性变换函数bilinear( ) 格式：numd,dendbilinear(num,den,Fs) 功能：把模拟滤波器的传递函数模型转换成数字滤波器的传递函数模型。 4.8.2双线性变换法 双线性变换利用频率变换关系： 例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 章 MATLAB在信号处理中的应用 4.8.3 IIR数字滤波器的频率变换设计法 1. IIR数字滤波器的频率变换设计法的基本原理 根据滤波器设计要求，设计模拟原型低通滤波器，然后进行频率变换，将其转 换为相应的模拟滤波器(高通、带通等)，最后利用冲激响应不变法或双线性变换法， 将模拟滤波器数字化成相应的数字滤波器。 图4.39 IIR数字滤波器MATLAB设计步骤流程图 第 4 章 MATLAB在信号处理中的应用 1MATLAB的典型设计 利用在MATLAB设计IIR数字滤波器可分以下几步来实现 (1)按一定规则将数字滤波器的技术指标转换为模拟低通滤波器的技术指标； (2)根据转换后的技术指标使用滤波器阶数函数，确定滤波器的最小阶数N和截止 频率Wc； (3)利用最小阶数N产生模拟低通滤波原型； (4)利用截止频率Wc把模拟低通滤波器原型转换成模拟低通、高通、带通或带阻 滤波器； (5)利用冲激响应不变法或双线性不变法把模拟滤波器转换成数字滤波器。 例4-39 设计一个数字信号处理系统，它的采样率为Fs100Hz，希望在该系统中 设计一个Butterworth型高通数字滤波器，使其通带中允许的最小衰减为0.5dB，阻 带内的最小衰减为40dB，通带上限临界频率为30Hz，阻带下限临界频率为40Hz。 第 4 章 MATLAB在信号处理中的应用 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所示。 第 4 章 MATLAB在信号处理中的应用 2MATLAB的直接设计 图4.39 IIR数字滤波器MATLAB设计步骤流程图 例4-41 试设计一个带阻IIR数字滤波器，其具体的要求是：通带的截止频率：wp1 650Hz、wp2850Hz；阻带的截止频率：ws1700Hz、ws2800Hz；通带内的 最大衰减为rp0.1dB；阻带内的最小衰减为rs50dB；采样频率为Fs2000Hz。 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 章 MATLAB在信号处理中的应用 4.9 FIR数字滤波器设计 格式：w = boxcar(M) 功能：返回M点矩形窗序列。 MATLAB信号处理工具箱中的窗函数法设计FIR数字滤波器的专用命令fir1( )。 格式：Bfir1(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 章 MATLAB在信号处理中的应用 例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.