第 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 序列的合成和截取 用例的结果编写产生矩形序列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)。 解：本题的计算原理见例，在这里用工具箱函数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(ny (n-1)要求用程序验证系统的线性性质。解：产生两种输入序列，分别乘以常数后1。分别激励系统，再求输出之和；2。先相加，再激励系统求输出；对两个结果进行比较，方法是求它们之差，按误差的绝对值是否极小进行判断。例7.6 离散序列的卷积计算给出两个序列 和 ,计算其卷积y(n)，并图示各输入输出序列。解：在例中，已经给出了直接调用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多项式分式的逆变换 设系统函数为，输入例中的x2信号，用z变换计算输出y(n)解：由例可知，故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 和 p = 0.7791 和 k = -150 -30代入得例7.9 离散时间傅里叶变换取周期的正弦信号，作8点采样，求它的连续频谱。然后对该信号进行N个周期延拓，再求它的连续频谱。把N无限增大，比较分析其结果。解： 先求离散傅立叶变换的MATLAB子程序最后得到X = x*exp(-j*w*n)。有了子程序，本例就没有什么难度了。例7.9 离散时间傅里叶变换2 程序运行结果执行程序q709并按提示键入N = 4，所得图形如图所示。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)，画出其幅度频谱。采样时间区间均为秒。为了便于比较，画出了幅度归一化的幅频曲线，如图所示。 例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 很容易写出程序。 例7.13 低通滤波及时域卷积定理 输入信号 x(nnnn(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：用类似于例中的方法，可把（）式写成矩阵乘法运算。其中，xn为序列行向量，Wnk是一NN阶方阵， 而 称为旋转因子。7.3 离散傅里叶变换（DFT）用矩阵乘法计算N点DFT的程序如下。 MATLAB程序%用矩阵乘法计算N点DFTclear;close allxn=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和（或）预滤波（在采样之前滤除折叠频率以外的频率成分）来改善。编写程序验证截断效应及加窗的改善作用，先选取以下参数： 采样频率fs = 400Hz，T = 1/fs 采样信号序列 对x(n)作4096点DFT作为的近似频谱Xa(jf)。例7.19 用DFT求连续信号频谱如图所示。图中X1(jf)，X4(jf)和X8(jf)分别表示Tp，0.04*4s和0.04*8s时的谱分析结果。由图可见，由于截断使原频谱中的单频谱线展宽（又称之为泄漏），Tp越大泄漏越小，频率分辨率越高。Tp时，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文件清单附在程序之后。 例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 0g = 2格型结构系数（反射系数）：例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。例的结构图与方程例的方程的矩阵形式由此可以求出系统函数例解出的系统函数程序运行的结果如果加入一个激励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 窗函数法设计数字滤波器对两种窗函数的设计结果分别如右图和图所示。 工具箱设计函数fir1和fir2MATLAB提供了基于窗函数法的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 决定，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,cheb2apbesselap，ellipap函数频率变换(变为高通，带通，带阻等)lp2lp， lp2hp， lp2bp，lp2bs模拟数字变换bilinearimpinvar合为一步的设计函数 butter， cheb1， cheb2，ellip， besself求最小阶数NButtord, cheb1ordCheb2ord,ellipord巴特沃斯滤波器设计流程（1）求最小阶数N的函数buttordN, wc = buttord (wp, ws, Rp, Rs, s) 根据滤波器指标wp，ws，Rp，Rs，求出巴特沃斯模拟滤波器的阶数N及频率参数wc，此处wp，ws及wc均以弧度/秒为单位。（2）得到N后,调用设计函数buttapz，p，k = buttap(N) 得到z, p, k后，很容易求出滤波器系数B，A。 （3）调用模拟频率变换函数lp2lpBt, 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 模拟低通转换为数字低通图形结果如图所示。由图（b）可见，对脉冲响应不变法，采样频率Fs越高（T越小），混叠越小；由图（c）可见，对双线性变换法，无频率混叠，但存在非线性失真。例7.31 切比雪夫数字滤波器设计 解：切比雪夫型滤波器通带内为等波纹，阻带内单调下降；切比雪夫型滤波器通带内为单调下降，阻带内等波纹。调用cheb2ord函数和cheby2函数使切比雪夫型设计变得非常简单。先用N, wc = Cheb2ord(wp, ws, Rp, Rs)求出N和 wc，提供函数cheby2的输入变元，再由B, A = cheby2(N, Rp, wc)设计切比雪夫型数字滤波器。B和A分别为H(z)的分子和分母多项式系数。 对切比雪夫型滤波器，同样有相应的工具箱函数cheb1ord和cheby1。 例7.32 椭圆带通数字滤波器设计设计一个椭圆带通数字滤波器，要求如下 Rp = 1 dB， Rs = 60 dB解：调用ellipord函数和ellip函数，代入参数wp，0.45；ws，0.55；Rp = 1；Rs = 60;即可得到本题的程序 。例7.33 高通和带通数字滤波器设计用双线性变换法从Ha(s)设计四阶带通巴特沃斯数字滤波器，并图示 （设采样周期T=1s）。指标如下：解：本题主要涉及三个步骤：（1）由数字滤波器指标求相应的模拟滤波器指标；（2）模拟滤波器频率变换（因为已给定阶数和模拟滤波器归一化低通原型）；（3）由相应的模拟滤波器到数字滤波器转换（双线性变换法）。例7.33 高通和带通数字滤波器设计首先要设计出与该指标相应的四阶Butterworth模拟滤波器。然后，调用bilinear函数将其转换成数字滤波器。对双线性变换法，由数字边界频率求相应的模拟边界频率时，一定要考虑预畸变矫正。这样，最终设计结果才能满足所给指标。步骤（1）： 设计高通滤波器时，模拟高通3dB截止频率为 设计带通滤波器时，模拟带通3dB截止频率为 例7.33 高通和带通数字滤波器设计步骤（2）：可调用MATLAB频率变换函数lp2lp，lp2hp，lp2bp，分别实现从模拟低通到模拟低通、高通、带通、带阻的频率变换。本题用到的lp2hp和lp2bp的格式简要说明如下： Bt，At = lp2hp(B, A, wc) 将系数向量为B和A的模拟滤波器归一化低通原型（3dB截止频率为1rad/s）变换成3dB截止频率为wc的高通模拟滤波器，返回高通模拟滤波器系数向量Bt和At。步骤（3）：调用bilinear函数将滤波器系数向量Bt和At转换到数字滤波器Bz和Az。