Matlab在数字信号处理中的应用.ppt

第 7章 在数字信号处理中的应用,7.1 离散信号的产生及时域处理,时域离散信号用x(n)表示,时间变量n（表示采样位置）只能取整数。因此,x(n)是一个离散序列,以后简称序列。用一个向量x不足以表示序列值x(n)。必须再用另一个等长的定位时间变量n。x和n同时使用才能完整地表示一个序列,由于n序列是按整数递增的,可简单地用其初值ns决定,因为它的终值nf取决于ns 和x的长度length(x),故可写成： n = ns:nf 或 n = ns: nslength(x)1,例7.1 序列的相加和相乘,给出两个序列x1(n)和x2(n)。 x1 = 0,1,2,3,4,3,2,1,0; n1 = -2:6,; x2 = 2,2,0,0,0,-2,-2; n2 = 2:8; 要求它们的和ya及乘积yp。 解：编程的思路是把序列长度延拓到覆盖n1和n2的范围，这样才能把两序列的时间变量对应起来，然后进行对应元素的运算。,例7.2 序列的合成和截取,用例6.13的结果编写产生矩形序列RN(n)的程序。序列起点为n0，矩形序列起点为n1，长度为N（n0，n1，N由键盘输入）。并用它截取一个复正弦序列exp(jn/8) . 解：建模:矩形序列可看成两个阶跃序列之差。 用MATLAB逻辑关系产生矩形序列x2(n)。而用它截取任何序列相当于元素群相乘x2.*x，也称为加窗运算。序列的合成和截取就是相加和相乘。,例7.3 序列的移位和周期延拓,已知 ，利用MATLAB生成并图示 表示x(n)以8为周期的延拓）和 解：方法1，利用矩阵乘法和冒号运算 x=1 2 3 4;y=x*ones(1，3); 方法2，采用求余函数mod， y = x(mod(n, M)+1)可实现对x(n) 以M为周期的周期延拓。加1是因为MATLAB向量下标只能从1开始，,例7.4 离散系统对信号的响应,本题给定6阶低通数字滤波器的系统函数，求它在下列输入序列x(n)下的输出序列y(n)。 解：本题的计算原理见例6.14，在这里用工具箱函数filter来解。如果已知系统函数H(z)=B(z)/A(z)，则filter函数可求出系统对输入信号x(n)的响应y(n)。 y = filter(B, A, x) 由差分方程可得到H(z)的分子和分母多项式系数向量A和B,再给出输入向量x即可。,例7.5 系统线性性质验证,设系统差分方程为 y(n) = x(n) + 0.8y (n-1) 要求用程序验证系统的线性性质。 解：产生两种输入序列，分别乘以常数后 1。分别激励系统，再求输出之和； 2。先相加，再激励系统求输出； 对两个结果进行比较，方法是求它们之差，按误差的绝对值是否极小进行判断。,例7.6 离散序列的卷积计算,给出两个序列 和 ,计算其卷积y(n)，并图示各输入输出序列。 解：在例6.4中，已经给出了直接调用MATLAB的卷积函数conv的方法，也给出了自编卷积计算程序的方法，要注意的是本例时间变量的设定和移位方法。在本例中，设定n为从零开始，向量x和h的长度分别为Nx=20和Nh=10；结果向量y的长度为length(y)=Nx+Nh-1。,求z的逆变换的方法,对于z变换分式 可以用部分分式法或长除法求其反变换。 用函数residuez可以求出它的极点留数分解 其中 r, p, k = residuez (B, A) 其反变换为：,例7.7 有限序列的z和逆z变换,两序列x1 = 1,2,3， n1 = -1:1 及x2 = 2,4,3,5，n2 = -2:1，求出x1与x2及其卷积x的z变换。 解：其z变换可写成 两个多项式乘积 可用conv函数来求得。n数组要自己判别。n的起点ns = ns1 + ns2 = 3，终点nf = nf1 + nf2 = 2。 n=ns:nf。由x和n即可得出X(z)。,例7.8 求z多项式分式的逆变换,设系统函数为 ， 输入例7.7中的x2信号，用z变换计算输出y(n) 解：由例7.7可知 ，故 Y(z)=X(z)W(z)= 其中nsy = 分母分子中z的最高幂次之差。 调用 r, p, k = residuez(B, A)，可由B，A求出r，p，k，进而求逆z变换，得,例7.8 z多项式分式逆变换(续),由程序算出nsy =-1，留数、极点分别为 r = -57.7581 和 204.7581 p = 0.7791 和 0.3209 k = -150 -30 代入 得,例7.9 离散时间傅里叶变换,取周期的正弦信号，作8点采样，求它的连续频谱。然后对该信号进行N个周期延拓，再求它的连续频谱。把N无限增大，比较分析其结果。 解： 先求离散傅立叶变换的MATLAB子程序 最后得到X = x*exp(-j*w*n)。 有了子程序，本例就没有什么难度了。,例7.9 离散时间傅里叶变换2, 程序运行结果 执行程序q709并按提示键入N = 4，所得图形如图7.10所示。N取得愈大，其峰值愈大，宽度愈窄。当N取得很大时，会出现内存不足的问题，这是用矩阵乘法做傅里叶变换的缺点。另外，因为那时峰值点处的宽度很窄，也会出现所选频点对不上峰值点的问题。所以对于N无限增大的情况，必须用fft函数来求。这时用连续频谱也没有意义了。这里用同样的横坐标把几种频谱进行对比，使读者更好地理解其关系。,例7.10 时域采样频率与频谱混叠,分别以采样频率fs=1000Hz，400Hz和200Hz对xa(t)进行等间隔采样，计算并图示三种采样频率下的采样信号及其幅频特性 解：程序分别设定4种采样频率fs=10kHz，1kHz，400Hz和200Hz，对xa(t)进行采样，得到采样序列xa(t)，xa1(n)，xa2(n)，xa3(n)，画出其幅度频谱。采样时间区间均为0.1秒。为了便于比较，画出了幅度归一化的幅频曲线，如图7.11所示。,例7.10 采样频率与频谱混叠(续),由于 由以上关系式可见，采样信号的频谱函数是原模拟信号频谱函数的周期延拓，延拓周期为2/T。如果以频率f为自变量（ = 2f），则以采样频率fs = 1/T为延拓周期。对频带限于fc的模拟信号xa(t)，只有当fs2fc时，采样后 才不会发生频谱混叠失真。这就是著名的采样定理,例7.11 由离散序列恢复模拟信号,用时域内插公式 其中 模拟用理想低通滤波器恢复的过程，观察恢复波形，计算出最大恢复误差。 解：这个公式与卷积公式相像，可以用向量和矩阵乘法来解决。,例7.11 由离散序列恢复模拟信号,xa = x * g (TNM) = x * G 其中 G = sinc(Fs * TNM) M表示在两个采样点之间增加的间隔数，使输出更密，更接近模拟信号。,例7.12 梳状滤波器零极点和幅特性,梳状滤波器系统函数有如下两种类型。 FIR型： IIR型： freqz 数字滤波器频率特性计算和绘制函数 zplane H(z)的零-极点图绘制。 解：调用函数freqz和zplane 很容易写出程序q712.m。,例7.13 低通滤波及时域卷积定理,输入信号 x(n) = cos(0.04n) + cos(0.08n) + cos(0.4n) + 0.3(n)，0n63 通过低通滤波器，计算滤波器对x(n)的响应输出y(n)，并图示x(n)和y(n)，观察滤波效果。 解：如前所述，只要求出H(z)=B(z)/A(z)的分子和分母多项式系数向量B和A，则可调用滤波器直接型实现函数filter对输入信号x(n)进行滤波。 y = filter(B, A, x),例7.14 用符号运算工具箱解z变换,解：无限长度时间序列的z变换和逆z变换都属于符号运算的范围。MATLAB的symbolic（符号运算）工具箱已提供了这种函数。如果读者已在计算机上安装了这个工具箱，可以键入以下程序。 MATLAB程序q714.m 其特点是程序的开始要指定符号自变量 syms z n a N w0 % 规定z,n,a为符号变量,7.3 离散傅里叶变换（DFT）,定义DFT： 用类似于例7.9中的方法，可把（7.3）式写成矩阵乘法运算。 其中，xn为序列行向量，Wnk是一NN阶方阵， 而 称为旋转因子。,7.3 离散傅里叶变换（DFT）,用矩阵乘法计算N点DFT的程序如下。 MATLAB程序q73a.m %用矩阵乘法计算N点DFT clear;close all xn=input(请输入序列x= ); N = length(xn); % n=0:N-1; k=n; nk=n*k; % 生成NN方阵 WN=exp(-j*2*pi/N); Wnk=WN.nk; % 产生旋转因子矩阵 Xk=xn*Wnk; % 计算N点DFT,例7.15 序列的离散傅立叶变换,求复正弦序列 余弦序列 正弦序列 的离散傅立叶变换，分别按N =16和N =8进行计算。绘出幅频特性曲线， 进行比较讨论。,例7.15 序列的离散傅立叶变换,在截取16点时，得到的是完整的余弦波形；而截取8点时，得到的是半截的余弦波形，当然有大量的谐波成分。,例7.16 验证N点DFT的物理意义,（1） ， 绘出幅频曲线和相频曲线。 （2）计算并图示x(n)的8点DFT。 （3）计算并图示x(n)的16点DFT。 解: 序列x(n)的点DFT的物理意义是 在0，2上进行点等间隔采样。 程序先密集采样，绘制出幅频曲线图。然后再分别做8点和16点DFT来验证这个采样关系。,例7.17 频域与时域采样对偶性,（1）产生三角波序列 （2）对M = 40，计算x(n)的64点DFT，并图示x(n)和X(k) = DFTx(n)，k = 0, 1, , 63。 （3）对（2）中所得X(k)在 0，2 上进行32点抽样得 （4）求 的32点IDFT，即 （5）绘出 的波形图，评述它与x(n)的关系。,例7.17 频域与时域采样对偶性,由于频域在0，2上的采样点数N（N = 32）小于x(n)的长度M（M = 40），所以，产生时域混叠现象，不能由X1(k)恢复原序列x(n)。只有满足NM时，可由频域采样X1(k)得到原序列x(n)。 这就是频域采样定理。对NM的情况，请读者自己编程验证。,例7.18 快速卷积,快速卷积就是根据DFT的循环卷积性质，将时域卷积转换为频域相乘，最后再进行IDFT得到时域卷积序列y(n)。其中时域和频域之间的变换均用FFT实现，所以使卷积速度大大提高。框图如下:,例7.19 用DFT求连续信号频谱,在计算机上用DFT对模拟信号进行谱分析时，只能以有限大的采样频率fs对模拟信号采样有限点样本序列（等价于截取模拟信号一段进行采样）作DFT变换，得到模拟信号的近似频谱。其误差主要来自以下因素： 截断效应（频谱泄露和谱间干扰） 频谱混叠失真 因素使谱分辨率（能分辨开的两根谱线间的最小间距）降低，并产生谱间干扰； 因素使折叠频率（fs / 2）附近的频谱产生较大失真。,例7.19 用DFT求连续信号频谱,加大截取长度Tp可提高频率分辨率；选择合适的窗函数可降低谱间干扰；而频谱混叠失真要通过提高采样频率fs和（或）预滤波（在采样之前滤除折叠频率以外的频率成分）来改善。 编写程序q719.m验证截断效应及加窗的改善作用，先选取以下参数： 采样频率fs = 400Hz，T = 1/fs 采样信号序列 对x(n)作4096点DFT作为的近似频谱Xa(jf)。,例7.19 用DFT求连续信号频谱,如图7.19所示。图中X1(jf)，X4(jf)和X8(jf)分别表示Tp = 0.04s，0.04*4s和0.04*8s时的谱分析结果。由图可见，由于截断使原频谱中的单频谱线展宽（又称之为泄漏），Tp越大泄漏越小，频率分辨率越高。Tp = 0.04s时，25Hz与50Hz两根谱线已分辨不清了。所以实际谱分析的截取时间Tp是由频率分辨率决定的。 另外，在本应为零的频段上出现了一些参差不齐的小谱包（称为谱间干扰）。谱间干扰的大小取决于加窗的类型。用矩形窗比用Hamming窗的频率分辨率高（泄漏小），但谱间干扰刚好相反。,例7.20 IIR滤波器直接型的转换,程序调用了信号处理工具箱函数tf2sos和扩展函数dir2par， sos, g = tf2sos(B, A) 实现从直接型到级联型（二阶分割形式）的转换。g为式中的增益，sos为L6阶矩阵，表示式中的系数。,例7.20 IIR滤波器直接型的转换, Cp, Bp, Ap = dir2par(B, A) 实现从直接型到并联型的转换。B为直接型H(z)的分子多项式系数向量，A为直接型H(z)的分母多项式系数向量；Cp，Bp，Ap的含义与扩展函数dir2par中的C，B，A相同。 dir2par中又调用了复共轭对比较函数cplxcomp。由于dir2par和cplxcomp是文献7中开发的，不属于MATLAB工具箱函数，所以将其M文件清单附在程序q720.m之后。,例7.20 IIR滤波器直接型的转换,根据计算结果, 级联型H(z)表达式： 级联型结构图。,例7.20 IIR滤波器直接型的转换,并联型结构H(z)表达式 并联型结构图,例7.21 直接型结构到格型梯形结构, tf2latc函数实现直接型到格型转换 K, C = tf2latc (B, A) 求出零-极点IIR系统格型梯形结构的格型参数向量K和梯形参数向量C（用A(1)归一化）。注意，当系统函数在单位圆上有极点时发生错误。 K = tf2latc (1, A) 求出全极点IIR系统的格型结构参数向量K。如果使用格式K, C = tf2latc (1, A)，返回的系数C为标量。 K = tf2latc (B) 求出FIR系统的格型梯形结构参数（反射系数）向量K（用H(z)的常数项B(1)归一化）。,例7.21 直接型结构到格型梯形结构,直接型系统函数 转换为格型梯形结构,例7.22 FIR滤波器直接型到其他型,系统函数为 调用信号处理工具箱函数tf2sos和tf2latc,给变元A赋值1,B=2,13/12,5/4,2/3即可. 级联型结构系数 sos = 1.0000 0.5360 0 1.0000 0 0 1.0000 0.0057 0.6219 1.0000 0 0 g = 2 格型结构系数（反射系数）： K = 0.2500 0.5000 0.3333,例7.22 FIR滤波器直接型到其他型,得出级联型为 格型结构为,例7.23 FIR格型到直接型转换,给定K = 2,1/4,1/2,1/3 ,用函数 latc2tf即可 由B=latc2tf(K)得到B,写出直接型结构,例7.24 系统函数的计算机推导,数字滤波器的网络结构图实际上也是一种信号流图。因此【例6.20】中介绍的方法和公式同样可以用来求离散域的数字滤波器的系统函数。不同的地方仅仅在于节点方程中出现了作为系数的符号变量z1，它将出现在系数矩阵中。MATLAB是不能处理上标变量的，因此在程序中设q=z 1，在计算完成后再人工地把结果中的q恢复为z1。,例7.24的结构图与方程,例7.24的方程的矩阵形式,由此可以求出系统函数,例7.24解出的系统函数,程序运行的结果 如果加入一个激励x(n)A(n)，则得出,7.5 FIR数字滤波器设计,滤波器的特性指标 用绝对值1,2表示; 用分贝Rp,Rs表示,（1）窗函数法设计FIR滤波器,先根据c和N求出相应的理想滤波器单位脉冲响应hd(n)。 第二步要选择合适的窗函数w(n)来截取hd(n)的适当长度（即阶数），以保证实现要求的阻带衰减；最后得到FIR滤波器单位脉冲响应h(n)=hd(n).*w(n),即其系数。,（2）等波纹最佳一致逼近法,（2）等波纹最佳一致逼近法: 信号处理工具箱采用remez算法实现线性相位FIR数字滤波器的等波纹最佳一致逼近设计。其优点是，设计指标相同时，使滤波器阶数最低；或阶数相同时，使通带最平坦，阻带最小衰减最大；通带和阻带均为等波纹形式，最适合设计片段常数特性的滤波器。其调用格式如下： b = remez(N, f, m, w, ftype) 其中N由remezord函数求出: N, fo, mo, w = remezord(f, m, dev, Fs) 输入变元dev为各逼近频段允许的波纹振幅。 remez函数可直接调用remezord返回的参数如下： b=remez(N, fo, mo, w),例7.25 窗函数法设计数字滤波器,分别用矩形窗和Hamming窗设计线性相位FIR低通滤波器。要求通带截止频率c = /4，单位脉冲响应h(n)的长度N = 21。绘出h(n)及其幅频响应特性曲线。 先求出相应的理想滤波器（本例应为理想低通）单位脉冲响应hd(n)，再根据阻带最小衰减选择合适的窗函数w(n)，最后得到FIR滤波器单位脉冲响应h(n)=hd(n).*w(n)。,例7.25 窗函数法设计数字滤波器,本题中，c = /4，N = 21，所以线性相位理想低通滤波器的单位脉冲响应为: 为了满足线性相位FIR滤波器条件h(n) = h(N-1-n)，要求 = (N-1)/2 = 10。 信号处理工具箱中有窗生成函数boxcar，hamming，hanning和blackman等。,例7.25 窗函数法设计数字滤波器,对两种窗函数的设计结果分别如右图7.25-1和图7.25-2所示。,工具箱设计函数fir1和fir2,MATLAB提供了基于窗函数法的FIR滤波器设计函数fir1和fir2,其功能及用法如下。 fir1功能：标准频率响应形状。 格式：b=fir1(N, wc, ftype, window)。 当wc=wc1，wc2时，是的带通滤波器。 当ftype=high时，设计高通FIR滤波器； 当ftype=stop时，设计带阻FIR滤波器。 fir2功能：任意频率响应形状。 格式：b = fir2(N, f, m, window),例7.26 窗函数法设计带通滤波器,使用fir1函数b = fir1(N, wc, window) 编程 参数c为行向量c = lp/，hp/ 根据阻带最小衰减Rs = 60dB选择窗函数类型和阶次。可以查上面列出的“窗函数设计滤波器时的阶数选择表”。选blackman窗，其滤波器阻带最小衰减可达到74dB，其窗口长度M由过渡带宽度B = 0.15 决定，Blackman窗设计的滤波器过渡带宽度为12/M，故M取80。因M = N+1，所以滤波器阶数N = 79。,例7.27 用remez函数低通滤波器,解：先由题意计算设计参数 f = 1/4,5/16，m = 1,0; dev的计算稍复杂一些，由于 所以 有了这几个参数就可以调用remezord和remez函数了.,例7.27 用remez函数低通滤波器,横线为-3dB，两条竖线分别位于频率/4和5/16。显然，通带指标稍有富裕，过渡带宽度和阻带最小衰减刚好满足指标要求。,程序输出的幅频特性,例7.28 remez函数设计高通滤波器,观察等波纹逼近法中加权系数w()及滤波器阶数N的作用和影响。期望逼近的滤波器通带为3/4，阻带为0，23/32。 解:在滤波器设计中，技术指标越高，实现滤波器的阶数也就越高。另外，对固定的阶数，通带与阻带指标可以互换，过渡带宽度与通带波纹和阻带衰减指标可以互换。 取f = 0, 3/4, 23/32, 1，m = 0, 0, 1, 1。其余参数分三种情况进行设计，N = 30，w = 1, 1；N = 30，w = 1, 5; N = 60，w = 1, 1。,例7.28 remez函数设计高通滤波器,程序运行结果如图 由图可见，w较大的频段逼近精度较高；w较小的频段逼近精度较低。N较大时逼近精度较高,N较小时逼近精度较低 。,7.6 IIR数字滤波器设计,IIR数字滤波器设计的主要方法是先设计低通模拟滤波器，进行频率变换，将其转换为相应的（高通、带通等）模拟滤波器，再转换为高通、带通或带阻数字滤波器。对设计的全过程的各个步骤，MATLAB都提供了相应的工具箱函数，使IIR数字滤波器设计变得非常简单。本节主要结合例题介绍这些IIR滤波器设计的工具箱函数。 IIR数字滤波器的设计步骤由以下的流程图来表示。下面以巴特沃斯滤波器设计函数为典型，介绍此流程图中函数的功能和用法。,IIR数字滤波器设计流程图,模拟低通滤波器原型设计 Buttap,cheb1ap,cheb2ap besselap，ellipap函数,频率变换(变为高通，带通，带阻等) lp2lp，lp2hp，lp2bp，lp2bs,模拟数字变换 bilinear impinvar,合为一步的设计函数 butter，cheb1，cheb2，ellip， besself,求最小阶数N Buttord, cheb1ord Cheb2ord,ellipord,巴特沃斯滤波器设计流程,（1）求最小阶数N的函数buttord N, wc = buttord (wp, ws, Rp, Rs, s) 根据滤波器指标wp，ws，Rp，Rs，求出巴特沃斯模拟滤波器的阶数N及频率参数wc，此处wp，ws及wc均以弧度/秒为单位。 （2）得到N后,调用设计函数buttap z，p，k = buttap(N) 得到z, p, k后，很容易求出滤波器系数B，A。 （3）调用模拟频率变换函数lp2lp Bt, At = lp2lp(B, A, wo) （4）调用模拟数字变换函数 Bd, Ad = bilinear (B, A, Fs),集成的数字滤波器设计函数,把（2）、（3）、（4）合为一步的数字滤波器设计函数butter(N, wc, ftype) B, A = butter (N, wc) 设计低通或带通数字滤波器系数B，A（当为带通滤波器时，第(1)类函数由wp = wp1, wp2会自动生成wc = w1, w2）。 B, A = butter (N, wc, high) 设计高通数字滤波器系数B，A。 B, A = butter (N, wc, stop) 设计带阻数字滤波器系数B，A。 butter(N, wc, ftype)还有零极增益和状态空间形式，读者可用help命令查阅。,例7.29 巴特沃斯模拟滤波器设计,设计一个低通巴特沃斯模拟滤波器，指标如下。 通带频率：fp = 3400Hz，最大衰减：Rp = 3dB 阻带频率：fs = 4000Hz，最小衰减：As = 40dB 解:它的系统函数完全由阶数N和3dB截止频率c决定。而N和c是由滤波器设计指标决定的。 取c = c1，通带指标刚好，阻带指标富裕； 取c = c2，则阻带指标刚好，通带指标富裕。 MATLAB工具箱函数buttord，butter就是根据以上公式编写的。因此就无需再记忆这些公式了。,模拟转换为数字：脉冲响应不变法,模拟滤波器离散化的基本方法有脉冲响应不变法和双线性变换法。 脉冲响应不变法及impinvar函数 单极点的N阶模拟滤波器Ha(s)，用部分分式展开为 脉冲响应不变法的数字化结果为 工具箱函数impinvar可实现以上计算，格式为 Bz, Az = impinvar(B, A, Fs),模拟转换为数字：双线性变换法, 双线性变换法函数bilinear 双线性变换法的原理是用 代换Ha(s)中的s值，得到H(z)。bilinear函数用来实现这个转换。其使用格式为 Bz, Az = bilinear(B, A, Fs) 脉冲响应不变法的缺点是存在频率混叠失真。双线性变换法可完全消除频率混叠失真，缺点是存在非线性频率失真。,例7.30 模拟低通转换为数字低通,已知一模拟滤波器的系统函数为 分别用脉冲响应不变法和双线性变换法将Ha(s)转换成数字滤波器系统函数H(z)，并图示Ha(s)和H(z)的幅频响应曲线。 程序中的核心语句是以下两条： d,c=impinvar(b,a,Fs) % 用impinvar 函数离散化 f,e=bilinear(b,a,Fs) % 用bilinear 函数离散化,例7.30 模拟低通转换为数字低通,图形结果如图7.30