数字信号处理实验指导书(正1)

数字信号处理 实验指导书

长沙理工大学电气与信息工程学院电子信息工程教研室

二零零柒年柒月

目 录

实验一：离散时间信号的时域分析........................................................................ 3 实验二：离散时间系统的时域分析........................................................................ 6 实验三：离散时间信号的频域分析........................................................................ 9 实验四：线性时不变离散时间系统的频域分析.................................................. 13 实验五： IIR数字滤波器的设计 .......................................................................... 17 实验六： FIR数字滤波器的设计 ......................................................................... 24 附录A MATLAB系统的常用概念 ....................................................................... 28 附录B 信号处理工具箱函数 .......................................................................... 33

数字信号处理研究数字序列信号的表示方法，并对信号进行运算，以提取包含在其中的特殊信息。近几十年来，由于在研究及应用两方面均取得了进展，数字信号处理领域已日趋成熟。由于计算机的大量使用，从而很容易向学生提供实际环境，以验证所学的概念和算法。

本指导书编程语言是MATLAB，它广泛应用于高性能数值计算和可视化。本书假定读者已具备MATLAB基础知识。前面的一些实验帮助学生理解信号处理的重要概念，后面以设计性实验项目为主，有利于加强对原理的理解并且加强对技术的应用。附录中给出了本书中用到的MATALB函数及简要解释。

实验一：离散时间信号的时域分析

一．实验目的

1．熟悉MATLAB中产生信号和绘制信号的基本命令。

2．熟悉序列的简单运算，如：加法、标量乘法、时间反转、延时、乘法等。

二．实验相关知识准备

1．用到的MATLAB命令 运算符：

： . + - * / ; %

基本矩阵：

i ones pi rand randn zeros 基本函数：

cos exp imag real 数据分析： sum 二维图形：

axis grid legend plot subplot stairs stem title xlable ylable clf

工具箱：sawtooth square

三．实验内容

1．序列的产生

(1) 程序1－1：单位抽样序列的产生和绘制

% Program P1_1

% Generation of a Unit Sample Sequence clf;

% Generate a vector from -10 to 20 n = -10:20;

% Generate the unit sample sequence u = [zeros(1,10) 1 zeros(1,20)]; % Plot the unit sample sequence stem(n,u);

xlabel('Time index n');ylabel('Amplitude'); title('Unit Sample Sequence'); axis([-10 20 0 1.2]);

(2) 程序1－2：正弦序列的产生和绘制

% Program P1_2

% Generation of a sinusoidal sequence n = 0:40; f = 0.1; A = 1.5;

phase = 0;

arg = 2*pi*f*n - phase; x = A*cos(arg); clf;

% Clear old graph

% Plot the generated sequence

stem(n,x); grid;

title('Sinusoidal Sequence'); xlabel('Time index n'); ylabel('Amplitude'); axis;

axis([0 40 -2 2]);

2．序列的运算（信号的平滑）

数字信号处理应用的一个常见例子是从被加性噪声污染的信号中移除噪声。假定信号s[n]被噪声d[n]所污染，得到一个含有噪声的信号x[n]=s[n]+d[n]。我们需要对x[n]进行运算，产生一个合理的逼近s[n]，对时刻n的样本求平均，产生输出信号是一种简单有效的方法。如：三点滑动平均的信号。

程序1－3 实现三点滑动平均的信号运算：

y[n]=(x(n-1)+x(n)+x(n+1))/3 程序如下:

% Signal Smoothing by Averaging clf; R = 51;

d = 0.8*(rand(R,1) - 0.5); % 产生随机噪声 m = 0:R-1;

s = 2*m.*(0.9.^m); % 产生未被污染的信号 x = s + d'; %产生被噪声污染的信号 subplot(2,1,1);

plot(m,d','r-',m,s,'g--',m,x,'b-.'); xlabel('Time index n');ylabel('Amplitude'); legend('d[n] ','s[n] ','x[n] '); x1 = [0 0 x];x2 = [0 x 0];x3 = [x 0 0]; y = (x1 + x2 + x3)/3; subplot(2,1,2);

plot(m,y(2:R+1),'r-',m,s,'g--'); legend( 'y[n] ','s[n] ');

xlabel('Time index n');ylabel('Amplitude');

3．实验问题回答

(1) 命令clf,axis,title,xlable和ylable的作用是什么？

(2) 程序1－2中正弦序列的频率是多少？怎样可以改变它？哪个参数控制该

序列的振幅？该序列周期是多少？

(3) 程序1－3中加性噪声d[n]是什么样的形式？语句x=s+d代表什么？信号

x1,x2,x3与信号x之间的关系是什么？

四．实验报告要求

1．按照实验内容要求完成相关实验程序，并得出相关的实验结果（包括图形结果）。

2．回答实验中提出的问题。

3．总结本次实验结果，按照实验报告格式要求，书写实验报告。

五．实验设备

PC机，MATLAB软件

实验二：离散时间系统的时域分析

一．实验目的

通过MATLAB仿真一些简单的离散时间系统，并研究它们的时域特性。

二．实验相关知识准备

1．用到的MATLAB命令 语言构造与调试：

break end for if input 基本函数： abs num2str 多项式和内插函数： conv 工具箱：

filter impz

三．实验内容

1．线性和非线性系统 例2－1 设系统为

y[n]-0.4y[n-1]+0.75y[n-2]=2.2403x[n]+2.4908x[n-1]+2.2403x[n-2]

要求用MATLAB程序仿真系统，输入三个不同的输入序列x1(n),x2(n)和 x(n)=a.x1(n)+b.x2(n)，计算并求出相应的输出响应y1[n],y2[n]和y[n]。

% Generate the input sequences clf; n = 0:40; a = 2;b = -3;

x1 = cos(2*pi*0.1*n); x2 = cos(2*pi*0.4*n); x = a*x1 + b*x2;

num = [2.2403 2.4908 2.2403]; den = [1 -0.4 0.75];

ic = [0 0]; % Set zero initial conditions

y1 = filter(num,den,x1,ic); % Compute the output y1[n] y2 = filter(num,den,x2,ic); % Compute the output y2[n] y = filter(num,den,x,ic); % Compute the output y[n]

yt = a*y1 + b*y2;

d = y - yt; % Compute the difference output d[n] % Plot the outputs and the difference signal subplot(3,1,1) stem(n,y);

ylabel('Amplitude');

title('Output Due to Weighted Input: a \\cdot x_{1}[n] + b \\cdot x_{2}[n]'); subplot(3,1,2) stem(n,yt);

ylabel('Amplitude');

title('Weighted Output: a \\cdot y_{1}[n] + b \\cdot y_{2}[n]'); subplot(3,1,3) stem(n,d);

xlabel('Time index n');ylabel('Amplitude'); title('Difference Signal');

2．时不变系统和时变系统

实验内容2－2：因果系统的时不变特性研究

编写MATLAB程序仿真例2.1给出的系统，以产生两个不同的输入序列x(n)（在例2－1中x(n)=a.x1(n)+b.x2(n)）和x(n+D)，计算并画出相应的输出序列y1[n],y2[n+D]和y1[n]-y2[n+D]

参数：n = 0:40; D = 10;a = 3.0;b = -2;

实验内容2－3：线性时不变系统的冲激响应计算

设系统为y[n]-0.4y[n-1]+0.75y[n-2]=2.2403x[n]+2.4908x[n-1]+2.2403x[n-2] 用y=impz(num,den,N)计算系统的冲激响应的前N个样本。

3. 实验问题回答

(1) 修改程序2－1的相关语句，对由加权输入得到的y[n]与相同权系数下输

出y1[n]和y2[n]相加得到的yt[n]进行比较，这两个序列是否相等？该系统是线性系统吗？

(2) 比较实验内容2－2中输出序列y[n]和yd[n-10]，该系统是时不变系统吗？ (3) 计算并画出实验内容2－3的前40个样本。

四．实验报告要求

1．按照实验内容要求完成相关实验程序，并得出相关的实验结果（包括图

形结果）。

2．回答实验中提出的问题。

3．总结本次实验结果，按照实验报告格式要求，书写实验报告。

五．实验设备

PC机，MATLAB软件

实验三：离散时间信号的频域分析

一．实验目的

1．在学习了离散时间信号的时域分析的基础上，对这些信号在频域上进行分析，从而进一步研究它们的性质。

2．熟悉离散时间序列的3种表示方法：离散时间傅立叶变换（DTFT）,离散傅立叶变换（DFT）和Z变换。

二．实验相关知识准备

1．用到的MATLAB命令 运算符和特殊字符：

< > .* ^ .^ 语言构造与调试：

error function pause 基本函数：

angle conj rem 数据分析和傅立叶变换函数： fft ifft max min 工具箱：

freqz impz residuez zplane

三．实验内容

1． 离散傅立叶变换

在MATLAB中，使用fft可以很容易地计算有限长序列x[n]的离散傅立叶变换。此函数有两种形式： y=fft(x)

y=fft(x,n) 求出时域信号x的离散傅立叶变换

n为规定的点数，n的默认值为所给x的长度。当n取2的整数幂时变换的速度最快。通常取大于又最靠近x的幂次。（即一般在使用fft函数前用n=2^nextpow2(length(x))得到最合适的n)。

当x的长度小于n时，fft函数在x的尾部补0，以构成长为n点数据。

当x的长度大于n时，fft函数将序列x截断，取前n点。

一般情况下，fft求出的函数多为复数，可用abs及angle分别求其幅度和相位。 注意：栅栏效应，截断效应（频谱泄露和谱间干扰），混叠失真

例3－1： fft函数最通常的应用是计算信号的频谱。考虑一个由100hz和200hz正弦信号构成的信号，受零均值随机信号的干扰，数据采样频率为1000hz。通过fft函数来分析其信号频率成分。

t=0:0.001:1;%采样周期为0.001s，即采样频率为1000hz

x=sin(2*pi*100*t)+sin(2*pi*200*t)+1.5*rand(1,length(t));%产生受噪声污染的正弦波信号 subplot(2,1,1);

plot(x(1:50));%画出时域内的信号 y=fft(x,512);%对x进行512点的fft

f=1000*(0:256)/512;%设置频率轴（横轴）坐标，1000为采样频率 subplot(2,1,2);

plot(f,y(1:257));%画出频域内的信号

实验内容3－2：混叠现象与频率分辨率

假设现有含有三种频率成分的信号x(n)=sin(2nπf1/fs)+sin(2nπf2/fs)+sin(2nπf3/fs),其中f1=20hz,f2=20.5hz,f3=40hz,采样频率fs=100hz.

（1）求x(n)在0～128之间的DFT的x(k)；

（2）求把（1）中的x(n)以补0方式使其加长到0～512之间后的DFT的x(k)； （3）求x(n)在0～512之间的DFT的x(k)。最后分析三种形式的DFT结果图给出结论。

提示：根据采样定理，在信号采样的过程中，如果采样频率小于信号最高频率的两倍，那么信号的各次调制频谱就会相互重叠起来，有些频率部分的幅值就与原始情况不同，因而不能分开和恢复这些部分，取样造成了信号的损失，这种频谱重叠的出现，就成为“混叠现象”。

令fs为采样频率，fh为信号的最高频率，为了避免混叠现象，则必须使

fs2fh，而

FF又为

fs 。所以，在fs和fh两个参数中，保持其中一个不变而增加另一个的唯一N方法，就是增加在一记录长度内的点数N。如果fs与fh都已给定，则N必须满足N2fh F2． 一维逆快速傅立叶变换 y=ifft(x) y=ifft(x,n)

实验内容3-3：频域采样定理的验证 （1） 产生一个三角波序列x(n)

0nM/2n x(n)MnM/2nM(2)对M=40，计算x(n)的64点DFT，并图示x(n)和X(k)=DFT[x(n)],k=0,1,…,63。 (3)对（2）中所得X(k)在[0,2pi]上进行32点抽样得

X1(k)=X(2k), k=0,1,…,31。并图示。

（4）求X1（k）的32点IDFT,并图示。即x1(n)=IDFT[X1(k)], k=0,1,…,31 (5)采用周期延拓的方法绘出x1((n))32的波形，评述x1((n))32与x(n)的关系，并根据频域采样理论加以解释。

3． 线性调频Z变换

离散傅立叶变换（DFT）可以看作信号在Z域上沿单位圆的均匀采样。但在实际应用中，并非整个单位圆上的频谱都有意义。一些情况下，如对于窄带信号，只希望分析信号所在的一段频带等，采样点的轨迹是一条弧线或圆周。这种需求，就导致了线性调频Z变换（Chirp z变换）的出现。

Chirp z变换与DFT计算整个频谱的算法不同，它是一种更为灵活的计算频谱的算法，可以用来计算单位圆上任一段曲线的Z变换，作频谱分析时输入的点数和输出的点数可以不相等，从而达到频域“细化”的目的。 y=czt(x,m,w,a) y=czt(x)

例3-4：利用Chirp z变换计算滤波器h（且h=fir1(30,125/500,boxcar(31))）在100Hz~200Hz的频率特性，并用图形文件比较CZT函数和FFT函数。 h=fir1(30,125/500,boxcar(31)); fs=1000; f1=100;f2=200; m=1024;

y=fft(h,1024); fy=fs*(0:1023)/1024; subplot(2,1,1); plot(fy,abs(y)); axis([0,500,0,1.5]);

w=exp(-j*2*pi*(f2-f1)/(m*fs)); a=exp(j*2*pi*f1/fs);

z=czt(h,m,w,a);%产生1024个z值 fz=(f2-f1)*(0:1023)/1024+f1;subplot(2,1,2); plot(fz,abs(z))

4． 实验问题回答

（1）完成例3－1的程序并观察信号序列x(t)在时域和频域上分析的特点。 （2）在实验内容3-2和3-3中按要求完成程序，并回答3-2和3-3中提出的问题。

（2） 完成例3－4的程序并观察Chirp z变换与fft变换的特点。

四．实验报告要求

1．按照实验内容要求完成相关实验程序，并得出相关的实验结果（包括图形结果）。

2．回答实验中提出的问题。

3．总结本次实验结果，按照实验报告格式要求，书写实验报告。

五．实验设备

PC机，MATLAB软件

实验四：线性时不变离散时间系统的频域分析

一．实验目的

1．在学习了离散时间系统的时域分析的基础上，对这些系统在频域上进行分析，掌握LTI系统的频域分析法。

2．熟悉传输函数和频率响应。

二．实验相关知识准备

1．用到的MATLAB命令

MATLAB的工具箱，为同一离散系统的不同系统模型间的多种表达方式的变换，提供了很多功能丰富快捷的MATLAB函数，如下表所示：

系统原模型 传递函数 传递函数 零极点增益 零极点增益 零极点增益 状态空间 状态空间 状态空间 二次分式 二次分式 多项式 Lattice结构 基本函数：

abs angle imag log10 real zplane可以查看零极点 poly2rc可以直接检测稳定性

系统转换模型 状态空间 零极点增益 状态空间 传递函数 二次分式 传递函数 零极点增益 二次分式 零极点增益 状态空间 Lattice结构 多项式 函数 Tf2ss Tf2zp Zp2ss Zp2tf Zp2sos Ss2tf Ss2zp Ss2sos Sos2tf Sos2zp Sos2ss Poly2rc 三．实验内容

1. 例4－1 求截短的理想低通滤波器的冲激响应

理想低通滤波器的冲激响应是无限的，因此实际上不可能实现。在实际应用中，将冲激响应截短到一个有限值，可以近似地实现该滤波器。然而，将冲激响应截短后滤波器表现为一个非因果的滤波器。将经过截短后的滤波器的冲激响应向右移动N/2个样本，可得因果上的近似，如下所示：

h(n)滤波器的长度为N+1。

sincnNnN22 , 0nN

下面这个使用函数sinc的程序，可用来计算上述近似滤波器的冲激响应。

clf; fc = 0.25; n = [-6.5:1:6.5];

y = 2*fc*sinc(2*fc*n);k = n+6.5;

stem(k,y);title('N = 13');axis([0 13 -0.2 0.6]); xlabel('Time index n');ylabel('Amplitude');grid;

实验内容4－2

①运行例4－1的程序计算并画出近似理想低通滤波器的冲激响应。 ②在例4－1中的低通有限冲激响应滤波器的长度是多少？该程序中哪个语句确定滤波器的长度？哪个参数控制截止频率？

③修改例4－1中的程序，计算长度为20，截止角频率为0.45的有限冲激响应低通滤波器的冲激响应。（h[n]不变）

2．例4－3 线性相位有限冲激响应滤波器的零点的位置研究。

下列程序先产生四类滤波器的冲激响应，然后产生零极点图，最终显示零点的位置。

clf;

b = [1 -8.5 30.5 -63]; num1 = [b 81 fliplr(b)]; num2 = [b 81 81 fliplr(b)]; num3 = [b 0 -fliplr(b)]; num4 = [b 81 -81 -fliplr(b)]; n1 = 0:length(num1)-1;

n2 = 0:length(num2)-1;

subplot(2,2,1); stem(n1,num1);

xlabel('Time index n');ylabel('Amplitude'); grid; title('Type 1 FIR Filter'); subplot(2,2,2); stem(n2,num2);

xlabel('Time index n');ylabel('Amplitude'); grid; title('Type 2 FIR Filter'); subplot(2,2,3); stem(n1,num3);

xlabel('Time index n');ylabel('Amplitude'); grid; title('Type 3 FIR Filter'); subplot(2,2,4); stem(n2,num4);

xlabel('Time index n');ylabel('Amplitude'); grid; title('Type 4 FIR Filter'); pause

subplot(2,2,1); zplane(num1,1); title('Type 1 FIR Filter'); subplot(2,2,2); zplane(num2,1); title('Type 2 FIR Filter'); subplot(2,2,3); zplane(num3,1); title('Type 3 FIR Filter'); subplot(2,2,4); zplane(num4,1); title('Type 4 FIR Filter');

disp('Zeros of Type 1 FIR Filter are'); disp(roots(num1));

disp('Zeros of Type 2 FIR Filter are'); disp(roots(num2));

disp('Zeros of Type 3 FIR Filter are'); disp(roots(num3));

disp('Zeros of Type 4 FIR Filter are'); disp(roots(num4));

实验内容4－4

①运行例4－3的程序，生成每一类线性相位有限冲激响应滤波器的冲激响应。并回答每一个有限冲激响应滤波器的长度是多少？

②在例4－3中验证冲激响应序列的对称性，每一个滤波器的零点位置及它们的线性相位特性。并回答这些滤波器的群延迟是多少？

四．实验报告要求

1．按照实验内容要求完成相关实验程序，并得出相关的实验结果（包括

图形结果）。

2．回答实验中提出的问题。

3．总结本次实验结果，按照实验报告格式要求，书写实验报告。

五．实验设备

PC机，MATLAB软件

实验五： IIR数字滤波器的设计

一．实验目的

1．加深对IIR滤波器的理解

2．学会用直接法和间接法设计满足响应要求的IIR数字滤波器

二．实验相关知识准备 1．模拟滤波器的设计函数

为了利用模拟滤波器设计数字滤波器，先必须设计对应的模拟滤波器，常用的模拟滤波器有：bessel滤波器，butterworth滤波器，chebyshevI型滤波器，chebyshevII型滤波器及椭圆函数滤波器。 1》．模拟滤波器的设计函数

1) 设计bessel模拟低通滤波器 [z,p,k]=besselap(n)

2）设计butterworth模拟低通滤波器 [z,p,k]=buttap(n)

3）设计chebyshevI型模拟低通滤波器 [z,p,k]=cheb1ap(n,Rp)

Rp：通带内的波纹系数，单位分贝 4）设计chebyshevII型模拟低通滤波器 [z,p,k]=cheb2ap(n,Rs)

Rs：阻带内的波纹系数低于通带Rs分贝 5）设计椭圆模拟滤波器 [z,p,k]=ellipap(n,Rp.Rs) 2》滤波器阶数的选择

下列函数除了能选择模拟滤波器的阶数外，同时也能选择数字滤波器的阶数。

1） 选择butterworth滤波器阶数 数字域：[n,Wn]=buttord(Wp,Ws,Rp,Rs) 模拟域：[n,Wn]=buttord(Wp,Ws,Rp,Rs,’s’) 2）选择chebyshevI型滤波器阶数 数字域：[n,Wn]=cheb1ord(Wp,Ws,Rp,Rs)

模拟域：[n,Wn]=cheb1ord(Wp,Ws,Rp,Rs,’s’) 3）选择chebyshevII型滤波器阶数 数字域：[n,Wn]=cheb2ord(Wp,Ws,Rp,Rs) 模拟域：[n,Wn]=cheb2ord(Wp,Ws,Rp,Rs,’s’) 4) 选择椭圆滤波器阶数

数字域：[n,Wn]=ellipord(Wp,Ws,Rp,Rs) 模拟域：[n,Wn]=ellipord(Wp,Ws,Rp,Rs,’s’)

注意：

n：返回符合要求性能指标的数字滤波器或模拟滤波器的最小阶数 Wn：滤波器的截至频率（即3db频率）

Wp：通带的截至频率，Ws：阻带的截至频率，单位rad/s。且均为规一化频率，即0Wp(Ws)1。1对应π弧度。

频率规一化：信号处理工具箱中使用的频率为莱奎斯特频率，根据香农定理，它为采样频率的一半，在滤波器设计中的截止频率均使用莱奎斯特频率进行规一化。规一化频率转换为角频率，则将规一化频率乘以π。如果将规一化频率转换为Hz，则将规一化频率乘以采样频率的一半。 3》．模拟频率变换

1）低通到低通模拟滤波器变换 [bt,at]=lp2lp(b,a,Wo)

lp2lp函数将截至频率为1rad/s（规一化截至频率）的模拟低通滤波器原型变换成截至频率为Wo的低通滤波器。 2）低通到带通模拟滤波器变换 [bt,at]=lp2bp(b,a,Wo,Bw)

lp2bp函数可将截止频率为1rad/s的模拟低通滤波器原型变换成具有指定带宽Bw和中心频率Wo的带通滤波器。其中Wo=sqrt(W1*W2)，Bw=W2-W1 ，W2高端截止频率，W1低端截止频率。 3）低通到高通模拟滤波器变换 [bt,at]=lp2hp(b,a,Wo)

4）低通到带阻模拟滤波器变换 [bt,at]=lp2bs(b,a,Wo,Bw)

5）besself模拟滤波器设计 [b,a]=besself(n,Wn) [b,a]=besself(n,Wn,’ftype’) [z,p,k]=besself(n,Wn) [z,p,k]=besself(n,Wn,’ftype’)

2．由模拟滤波器变换成等效的数字滤波器

方法一：双线性变换法实现模拟滤波器到数字滤波器的变换

双线性变换将S域映射成Z域，从而将模拟滤波器变换成等效的数字滤波器 [zd,pd,kd]=bilinear(z,p,k,Fs)为零极点增益表示的bilinear函数

其中z,p,k为S域传递函数的零点、极点、和增益，Fs为取样频率，zd,pd,kd为双线性变换后Z域传递函数的零点、极点和增益。

[numd,dend]=bilinear(num,den,Fs,Fp)为传输函数表示的bilinear函数 方法二：冲激响应不变法实现模拟到数字滤波器的变换 [bz,az]=impinvar(b,a,Fs)

[bz,az]=impinvar(b,a)采用默认值为1hz的Fs。 小结：

设计IIR数字滤波器的一般步骤：

1）把给出的数字滤波器的性能指标转换为模拟滤波器的性能指标

2）根据转换后的性能指标，通过滤波器阶数选择函数，来确定滤波器的最小阶数N和固有频率Wn

3）由最小阶数N得到低通滤波器原型

4）由固有频率Wn把模拟低通滤波器原型转换为低通、高通、带通、带阻滤波器

5）运用脉冲响应不变法或双线性变换法把模拟滤波器转换成数字滤波器

MATLAB信号处理工具箱提供了几个用于直接设计IIR数字滤波器的函数，这些函数把上面这些复杂的步骤融合为一个整体，为设计滤波器带来了极大的方便。

3．直接设计IIR数字滤波器

1》Butterworth模拟和数字滤波器设计

数字域：[b,a]=butter(n,Wn)可设计出截止频率为Wn的n阶butterworth滤波器

[b,a]=butter(n,Wn,’ftype’)当ftype=high时，可设计出截止频率为Wn的高通滤波器；当ftype＝stop时，可设计出带阻滤波器 [z,p,k]=butter(n,Wn) [zp,k]=buter(n,Wn,’ftype’) [A,B,C,D]=butter(n,Wn) [A,B,C,D]=butter(n,Wn,’ftype’)

模拟域：[b,a]=butter(n,Wn,’s’)可设计出截止频率为Wn的n阶模拟butterworth滤波器，

其余形式类似于数字域的。 2》chebyshevI型滤波器（通带等波纹）设计

数字域：[b,a]=cheby1(n,Rp,Wn)可设计出n阶chebyshevI滤波器,其截止频率由Wn确定，通带内的波纹由Rp确定

[b,a]=cheby1(n,Rp,Wn,’ftype’)当ftype=high时，可设计出截止频率为Wn的高通滤波器；当ftype＝stop时，可设计出带阻滤波器 [z,p,k]=cheby1(n,Rp,Wn) [zp,k]= cheby1 (n,Rp,Wn,’ftype’) [A,B,C,D]= cheby1 (n,Rp,Wn) [A,B,C,D]= cheby1 (n,Rp,Wn,’ftype’)

模拟域：[b,a]= cheby1 (n,Rp,Wn,’s’)可设计出截止频率为Wn的n阶chebyshevI型模拟滤波器，其余形式类似于数字域的。 3》chebyshevII型滤波器（阻带等波纹）设计

数字域：[b,a]=cheby2(n,Rs,Wn)可设计出n阶chebyshevI滤波器,其截止频率由Wn确定，阻带内的波纹由Rs确定

[b,a]=cheby2(n,Rs,Wn,’ftype’)当ftype=high时，可设计出截止频率为Wn的高通滤波器；当ftype＝stop时，可设计出带阻滤波器 [z,p,k]=cheby2(n,Rs,Wn) [zp,k]= cheby2(n,Rs,Wn,’ftype’) [A,B,C,D]= cheby2 (n,Rs,Wn) [A,B,C,D]= cheby2(n,Rs,Wn,’ftype’)

模拟域：[b,a]= cheby2(n,Rs,Wn,’s’)可设计出截止频率为Wn的n阶chebyshevII型模拟滤波器，其余形式类似于数字域的。

4》椭圆滤波器设计

数字域：[b,a]=ellip(n,Rp,Rs,Wn)可设计出n阶椭圆滤波器,其截止频率由Wn确定，通带内的波纹由Rp确定,阻带内的波纹由Rs确定。

[b,a]=ellip(n,Rp,Rs,Wn,’ftype’)当ftype=high时，可设计出截止频率为Wn的高通滤波器；当ftype＝stop时，可设计出带阻滤波器 [z,p,k]=ellip(n,Rp,Rs,Wn) [zp,k]=ellip(n,Rp,Rs,Wn,’ftype’) [A,B,C,D]=ellip(n,Rp,Rs,Wn) [A,B,C,D]=ellip(n,Rp,Rs,Wn,’ftype’)

模拟域：[b,a]=ellip(n,Rp,Rs,Wn,’s’)可设计出截止频率为Wn的n阶椭圆型模拟滤波器，其余形式类似于数字域的。 5》基于幅度特性设计IIR数字滤波器

[b,a]=yulewalk(n,f,m)可得到以b,a表示的n阶IIR滤波器，其幅频特性逼近由频率矢量f和幅度矢量m所给定的特性。矢量f和n的长度必须相同。 6》基于频率特性设计IIR数字滤波器

[b,a]=invfreqz(h,w,nb,na)可得到以b,a表示的IIR滤波器，其幅频特性逼近由频率矢量w和复频率响应矢量h所给定的特性。矢量w和h的长度必须相同，nb,na分别为指定的分子和分母多项式的阶数。

三．实验内容 1． 例5－1

1》设计一butterworth低通滤波器，通带范围0～100hz，通带波纹小于3db，阻带为-30db，并利用最小的阶数来实现。采样频率为1000hz。 wp=100/500;ws=200/500;%频率归一化 [n,wn]=buttord(wp,ws,3,30); [b,a]=butter(n,wn); freqz(b,a,512,1000)

2》要求与1》相同，但采用椭圆函数滤波器实现。 wp=100/500;ws=200/500;%频率归一化 [n,wn]=ellipord(wp,ws,3,30);

[b,a]=ellip(n,3,30,wn); freqz(b,a,512,1000)

3》应用最小二乘拟合的方法设计一个8阶的低通滤波器，指定的频率点与对应的幅度分别由f和m给出，并画出它希望的频率响应和实际的频率响应。 f=[0 0.6 0.6 1]; m=[1 1 0 0];

[b,a]=yulewalk(8,f,m); [h,w]=freqz(b,a,128); plot(f,m,w/pi,abs(h),'--');

legend('ideal','yulewalk designed');

title('comparison of frequency response magnitudes')

2．实验内容5－2

确定通带截止频率为500HZ，阻带截止频率为600HZ以及通带的最大衰减量为1db，阻带的最小衰减量为50db，采样频率为2000HZ的数字低通巴特沃斯滤波器的阶数，并画出其幅频特性图。

实验内容5－3

采用直接设计法设计一个通带截止频率为500HZ～800HZ、阻带截止频率为600HZ～700HZ、通带内波纹的最大衰减为0.1db、阻带内波纹的最小衰减为50db、采样频率为2000HZ的带阻椭圆数字滤波器，并画出其幅频响应图。 实验内容5－4

模拟信号 xa（t）=sin(500*pi*t)+cos(50*pi*t)

有系统完成xa（t） A/D H(z) D/A ya(t)的处理，采样频率为1000HZ。

（1） 采用双线性变换法设计一个高通切比雪夫2型数字滤波器，要求所设计的滤波器阶数最小，以小于1db的衰减通过150hz的分量，以至少40db抑制100hz分量。滤波器应有单调的通带和等波纹的阻带。求此滤波器并画出幅频响应图（单位db）

（2） 产生上述信号xa（t）的150个样本点，用图形文件在时域和频域上分别

表示该信号。

（3） 然后通过上述所设计的滤波器对（2）所产生的信号滤波，并分别在时域和频域上结合滤波后的输出信号进行结果分析。

3．实验问题回答

（1）完成例5－1的程序并观察IIR滤波器的特点。

（2）在实验内容5-2、5-3和5－4中按要求完成程序，并回答提出的问题。

四．实验报告要求

1．按照实验内容要求完成相关实验程序，并得出相关的实验结果（包括图形结果）。

2．回答实验中提出的问题。

3．总结本次实验结果，按照实验报告格式要求，书写实验报告。

五．实验设备

PC机，MATLAB软件

实验六： FIR数字滤波器的设计

一．实验目的

1．加深对FIR滤波器的理解

2．学会用窗函数法设计满足响应要求的FIR数字滤波器

二．实验相关知识准备

1．基于窗函数设计具有标准频率响应的FIR滤波器

1）常用的窗函数

矩形序列：w=boxcar(n) 可产生一长度为n的矩形窗函数 三角窗函数：w=triang(n) 可产生n点的三角窗函数 bartlett窗：w=Bartlett(n) hamming窗：w=hamming(n) hanning窗：w=hanning(n) blackman窗：w=blackman(n)

chebyshev窗：w=chebwin(n,r) 可产生n点的chebyshev窗函数，其傅立叶变换后的旁瓣波纹低于主瓣rdb。注意，当n为偶数时，窗函数的长度为n+1 kaiser窗：w=Kaiser(n,beta) 可产生N点的kaiser窗函数，其中beta为影响窗函数旁瓣的β参数。若阻带最小衰减为As。

As2100.5842(As21)0.40.07886(As21)0.1102(As8.7)As5021As50

2）窗函数选择原则

 具有较低的旁瓣幅度，尤其是第一旁瓣幅度。  旁瓣幅度下降数率要大，以利增加阻带衰减。  主瓣的宽度要窄，以获得较陡的过渡带。

 通常上述几点很难同时满足，当选用主瓣宽度较窄时，虽然得到较陡的

过渡带，但通带和阻带的波纹明显增加；当选用最小的旁瓣幅度时，虽然得到匀滑的幅度相应和较小的阻带波动，但过渡带加宽。因此，实际选用的窗函数往往是它们的折中。在主瓣宽度达到一定要求的条件下，

适当牺牲主瓣宽度来换取旁瓣波动的减少。 3) 调用的函数fir1

fir1函数实现加窗线性相位FIR数字滤波器设计，可设计出标准的低通、带通、高通和带阻滤波器。

b=fir1(n,Wn) 可设计n阶低通FIR滤波器，这里Wn为规一化截止频

率，0≤Wn≤1,Wn=1相应于0.5fs。所设计的滤波器为加hamming窗的线性相位滤波器。当Wn=[W1 W2]时，fir1可得到带通滤波器，其通带为W1b=fir1(n,Wn,’ftype’) 当ftype=high时，设计高通FIR滤波器；

当ftype=stop时，设计带阻滤波器。其中n必须为偶数，当输入的阶次n为奇数时，fir1函数会自动将阶数加1。

b=fir1(n,Wn,window) window长度为n+1，如不指定window参数，

则fir1函数采用hamming窗。

b=fir1(n,Wn,‘ftype’,window)

2．基于频率采样法和窗函数设计具有任意频率响应的FIR滤波器

fir2函数用于设计具有任意频率响应的加窗数字FIR滤波器，其基本思

想是基于对滤波器的幅度响应进行频率采样。

b=fir(n,f,m) 可设计出一n阶的FIR滤波器，其滤波器的频率特性由矢量

f和m决定，f为频率采样点矢量，m为与之对应的幅度矢量，有关它们的约定可参见yulewalk函数。

b=fir2(n,f,m,window) 可设计具有任意频率响应形状的加窗线性相位FIR数字滤波器。其幅频特性由频率点向量f和幅度值向量m给出， 0≤f≤1，要求f为单增向量，而且从0开始，以1结束，1表示数字频率w=π。m与f等长度，m(k)表示频点f(k)的幅频响应值。Plot(f,m)可画出期望逼近的幅频响应曲线。如省略window，则自动取为hamming窗。

b=fir2(n,f,m,npt) 格式中，所希望的频率响应被等间隔地插值，长度变为

npt点，默认为512点。相应的b=fir2(n,m,npt,window)格式中，可指定窗函数。

b=fir2(n,f,m,npt,lab) 格式中，当所指定的f中，如果有2个相继的值相

同（例如在转折点的频率），可利用参数lab指定fir2在这2个频率点间插入必要的点数，以保证滤波器的频率响应能平滑过渡且具有必要的陡度，默认值为25点。相应的b=fir2(n,f,m,npt,lab,window)格式中，可指定窗函数。

三．实验内容 1． 例6－1

1》设计一个34阶的高通FIR滤波器，截止频率为0.48rad，并使用具有30db波纹的chebyshev窗。 window=chebwin(35,30); b=fir1(34,0.48,'high',window); freqz(b,1,512)

2》设计一个30阶的低通滤波器，使之与期望频率特性相近。 f=[0 0.6 0.6 1]; m=[1 1 0 0]; b=fir2(30,f,m); [h,w]=freqz(b,1,128); plot(f,m,w/pi,abs(h)) 2．实验内容6－2

设计一个阻带为0.4~0.7，阶数为38，窗函数为切比雪夫窗的带阻滤波器，并与窗函数为默认的hamming窗时的设计结果相比较，要求用图形文件分析比较结果并给出结论。

实验内容6－3

设计一个具有指定幅频响应的多带FIR滤波器，并与期望的幅频响应的结果进行比较，用图形文件说明比较结果。 实验内容6－4

设计一个最小阶次的低通FIR数字滤波器，性能指标为：通带0Hz～1500Hz,阻带截止频率2000Hz,通带波动1％，采样频率为8000Hz.用图形文件表示设计出的FIR低通滤波器的幅频响应曲线。 3．实验问题回答

（1）完成例6－1的程序并观察FIR滤波器的特点。

（2）在实验内容6-2、6-3和6－4中按要求完成程序，并回答提出的问题。

四．实验报告要求

1．按照实验内容要求完成相关实验程序，并得出相关的实验结果（包括图形结果）。

2. 回答实验中提出的问题。

3. 总结本次实验结果，按照实验报告格式要求，书写实验报告。

五．实验设备

PC机，MATLAB软件

附录A MATLAB系统的常用概念

1、命令窗口

在Windows 2000下启动MATLAB系统后，Windows 2000的工作平台上会弹出一个窗口，如下图所示，这个窗口称为MATLAB的命令窗口。MATLAB的命令窗口是用户与MATLAB解释器进行通信的工作环境，提示符‘>>’表示MATLAB解释器正等待用户输入命令。所有的MATLAB命令、MATLAB函数，以及MATLAB程序都要在这个窗口下运行。

在命令窗口，用户可以发出MATLAB命令。例如，为了生成一个3*3的矩阵，可以在提示符下，键入如下的命令：

A=[1 2 3;4 5 6; 7 8 9]

方括号命令表示矩阵，空格或逗号将每行的元素分开，而分号将矩阵的各行数值分开。再键入Enter后，MATLAB将回显如下的矩阵：

A=

1 2 3 4 5 6 7 8 9

为了求该矩阵的逆矩阵，则只要键入命令 ？B=inv(A);

MATLAB就将计算出相应的结果。如果不想在命令窗口中显示计算结果，只要如上所示，在该命令后多键入一个分号即可。此时，MATLAB系统只完成该命令所要求的计算任务，其计算结果不回显。这项功能在程序设计中是非常必要的。

MATLAB系统也可以说是一种新的语言，该语言十分容易掌握，其结构非常类似于数学式子的书写格式，用户花上很少的时间即可掌握MATLAB的大部分命令。

2、 图形窗口

MATLAB系统的强大功能之一是其优秀的图形功能。对于任何作图命令，MATLAB将打开另一个窗口来绘制输出图形，这样的窗口在MATLAB系统中被成为图形窗口。

在同一个图形窗口中，可以绘制多个图形，也可以生成多个图形窗口，并选择其中的一个图形窗口，在其中绘制图形。生成图形窗口的方法比较多，在没有图形窗口存在时，每个绘图函数都能自动生成一个图形窗口；也可以用figure 命令生成一个新的图形窗口；还可以用命令窗口的File菜单的New子菜单的Figure项来打开一个新的图形窗口。

3、搜索路径

MATLAB管理着一条搜索路径，它在搜索路径下寻找与命令相关的函数文件。例如，如果在MATLAB提示符下输入example, MATLAB解释器将按照下面的步骤来处理这条字符串：

（1）检查example是不是一个变量；

（2）如果不是，检查example是不是一个内部函数；

（3）如果不是，检查在当前文件夹下是否存在名为example.mex,example.dll或

example.m的文件。MEX文件是MATLAB 的执行文件，将优先执行；

（4）如果不存在，检查在MATLAB 的搜索路径的目录下是否存在名为

MEX,example.mex,example.dll,或example.m的文件。MEX文件优先执行。

用户可以打开路径浏览器（Path Browse）查看MATLAB系统的当前搜索路径，也可以在其中加入自己的路径。

4、 文件类型

在MATLAB系统中，根据功能可将MATLAB系统所使用的外部文件分为几类，并用不同的扩展名作为其标识，我们用的主要是M文件。M文件以字母m为其扩展名，例如startup.m。一般说来，M文件是以ASCII码文本文件，可以用任何文本编辑器进行编辑。在MATLAB系统中，有两类M文件。一类称为程序M文件，简称M文件；另一类称为函数M文件，或简称为函数，统称为 M文件。M文件的内容是由符合MATLAB语法的语句构成的，函数M文件的第一行必须是以关键字function开始的函数说明语句。两类M文件的共同特征是：在MATLAB命令窗口中的命令提示符下键入文件名，来执行M文件中的所有语句规定的计算任务或完成一定的功能。它们的区别在于以下两方面：第一，程序M文件中创建的变量都是MATLAB工作空间中的变量，工作空间中的其他程序或函数可以共享，而函数M文件中创建的所有变量除了全程变量外，均局限于函数运行空间内的局部变量；第二，函数M文件的调用式中可以有输入参数和输出参数，而程序M文件则没有这种功能。

5、语言语法要素

MATLAB只管理一种对象——矩阵。可以使用下列的任何一种方法在MATLAB环境下创建或输入一个矩阵：

1) 显示的输入一个元素序列；

2) 用MATLAB的内部函数创建一个矩阵；

3) 在M文件中用MATLAB语句创建一个矩阵； 4) 从一个外部数据文件中装载并创建一个矩阵。

在MATLAB中有两个基本概念：变量和表达式。变量由变量名表示，函数名作为特殊的变量名看待，每个变量名由一个字母后面跟随任意个字母或数字（包括下划线）组成，但

MATLAB只能分辨前19个字符。MATLAB能区分组成变量名的大小写字母。

MATLAB的语句则是两种形式之一：变量名=表达式或者表达式。 在前一种语句形式下，MATLAB将运算的结果赋给“变量名”；而在第二种语句形式下，将运算的结果赋给MATLAB的永久变量ans，每条语句以回车符结束。一般的，运算的结果在命令窗口中显示出来。如果语句的最后一个字符是分号“；”，那末，MATLAB仅执行赋值运算，不再显示运算的结果。

与C语言一样，MATLAB将字符串当作数组或矩阵处理。

6、MATLAB的基本运算符 矩阵运算符 A' A+B,A-B A*B A.*B < > <= >= == ~= 矩阵A的转置。 矩阵A与B的和与差。 矩阵A与B的乘法。 矩阵A与B的对应元素相乘 小于 小于或等于 大于 大于或等于 等于 不等于 关系运算符 7、特殊运算符

在MATLAB的M文件中，可以加入解释行。解释行的标识符为“%”，该标识符将被作为注解内容。程序执行时，注解被忽略。

方括号“[ ]”用于生成矩阵。特别的，语句A= [ ]生成空矩阵A。 行分隔符“；”用于MATLAB语句后时，表示该语句的执行结果不被回显出来，这可避免显示一些不感性趣的结果。

冒号“：“最主要的作用是生成向量，从下面的例子中可以看出它的使用方法： j:k 生成向量[j,j+1,j+2,…,k]

j:i:k 生成向量[j,j+i,j+2*i,…,k]，如果j>k，则生成空矩阵 A(:,j) 矩阵A的第j列 A(I,:) 矩阵A的第I行

A(j:k) 向量A(j),A(j+1),…,A(k)

A(:,j:k) 从第 j列到第k列的矩阵子块 换行连接符“…”，有时一条MATLAB语句会很长，在命令窗口的一行内很可能写不下，此时只要在该语句中加入三连点，再回车即可在下一行接着写该语句。

8、MATLAB的常用数学函数

sin Cos Tan Asin Atan Sinh 正弦函数 余弦函数 正切函数 反正弦函数 反正切函数 双曲正弦函数 三角函数 Cosh Tanh Asinh Acosh Atanh Acos Abs Angle Sqrt Real Imag Conj Round Ceil Rem Exp Log log10 双曲余弦函数 双曲正切函数 反双曲正弦函数 反双曲余弦函数 反双曲正切函数 反余弦函数 实数的绝对值、复数的模、字符串的ASIIC值 复数的幅角 方根函数 复数的实部 复数的虚部 复共轭运算 最邻近整数截断（四舍五入） 不大于自变量的最大整数 不小于自变量的最小整数 自然指数函数（以e为底） 自然对数函数（以e为底） 以10为底的对数函数 初等函数 9、程序流控制

与其他的程序设计语言一样，MATLAB语言也提供了条件语句。下面分别予以介绍。 1) for循环语句

MATLAB也有自己的for循环语句。如果要反复执行的一组语句的循环次数是已知或预定义的，就可以用for循环语句。例如：

for I=1:n x(I)=0; end

这条语句将向量x的前n个元素赋予零值，这里的变量n必须预先给定。 注意：每一个for必须与end配对使用。 2) while循环语句

MATLAB提供有while循环语句，它的作用是在一定的逻辑条件的控制下，不断的循环执行一条或一组语句，直到逻辑条件不满足为止。While 语句的一般形式是

While 表达式 语句组 end

3) if条件语句和bread语句

break 语句用于退出循环体，if条件语句有两种形式，分别是 if 表达式 语句组1 else

语句组2 end

和

if 表达式1 语句组1 elseif 表达式2 语句组2 else

语句组3 end

10、MATLAB的在线帮助

用户可以随时利用MATLAB的在线帮助查询自己不懂得用法的函数的具体用法，例如：在命令窗口键入help abs 后的显示如下：

ABS Absolute value ABS(X) is the absolute value of the elements of X. When X is complex, ABS(X) is the complex modulus (magnitude) of the elements of X.

See also SIGN,ANGLE,UNWRAP. Overloaded methods Help sym/abs.m

将abs函数的主要用法和用途都列了出来。

11 实验上机的具体过程如下：

1) 在windows 2000/XP的桌面上找到MATLAB的图标单击进入MATLAB的命

令窗口，或在开始菜单里选择程序再找到MATLAB单击也可进入MATLAB的命令窗口。

2) 进入命令窗口后，在菜单File中选择open，打开已存在的文件。如果是新文

件，在菜单File中选择New，再找到M-file即可打开MATLAB的编辑窗口，在编辑窗口内输入你的源程序后存盘。

3) 为了使程序能够在MATLAB中运行，需要在搜索路径中加入你的路径，加

入路径的过程如下，在MATLAB的命令窗口中选择File 菜单的Set path选

项，则打开一个Path Browser窗口，如下图所示：

在该窗口中单击Browse，选择你的路径加入，然后在File中选择Save Path后退出即可。这样，你可以在MATLAB的命令窗口中键入你的源程序文件名，或在编辑窗口中选择Tools中的Run即可编译运行。若编译无错，则可得出结果；若编译有错，则可根据命令窗口的提示进行修改后再编译运行，直至得出正确的结果。

附录B 信号处理工具箱函数

MATLAB包含了进行信号处理的许多工具箱函数，有关这些工具箱函数的使用可通过Help命令得到。为使用方便，在这里将给出几个常用到的函数的使用说明。 函数形式 X=sawtooth（t,width） 产生锯齿波width用于确定最大值的位置，即从0到 函数功能 关于函数参数的说明 或三角波。 2*width函数从-1上升到+1。 X=square(t,duty) Y=abs(x) 产生方波 求绝对值 Duty用于指定正半周期的比例 当x为复数时，得到的是复数模（幅值），若x为字符串，得到的是各个字符的ASCII码。 C=conv(a,b) 求卷积 求取矢量a和b的卷积,c的长度为a和b的长度和减去1。 [h,w]=freqs(b,a,n) 模拟滤波器.b,a为滤波器的冲击响应s变换的分子和分的频率响应 母多项式的系数，在n个频率点计算频率响应 [h,f]=freqz(b,a,n,Fs) 数字滤波器的频率响应。 Fs为采样频率，b,a为滤波器的冲击响应的Z变换的分子和分母多项式的系数，该函数的作用是在0~Fs/2频率范围内选取n个点（记在f中），并计算相应的频率响应。 [h,t]=impz(b,a,n) 数字滤波器b,a为滤波器的冲击响应s变换的分子和分的冲击响应 母多项式的系数，计算出冲击响应h，取样点树为n. [n,Wn]=buttord( Wp,Ws,Rp,Rs,’s’) Butterworth滤波器阶的选择。 Wp和Ws分别为通带和阻带的截止频率，皆大于0小于1。Rp和Rs分别是通带和阻带的波纹系数，’s’表示模拟域，也可不加‘s’，则为数字域。 [b,a]=butter(n,Wn,’ftype’,’s’) Butterworth模拟和数字滤波器设计。 [n,Wn]=cheblord( Wp,Ws,Rp,Rs,’s’) chebyshevI滤波器阶的选择。 设计阶数为n，截止频率为Wn的滤波器，ftype指滤波器的类型，‘high’是高通，’stop’是带阻，无此参数则是低通，’s’指模拟域，无则表示数字域，b,a是对应变换的分子分母多项式的系数。 Wp和Ws分别为通带和阻带的截止频率，皆大于0小于1。Rp和Rs分别是通带和阻带的波纹系数，’s’表示模拟域，也可不加‘s’，则为数字域。 [b,a]=cheby1(n,RpWn,Chebyshev(’ftype’,’s’) 切比雪夫)I型模拟和数字滤波器设计（通带等波纹）。 [b,a]=cheby2(n,RsWn,’ftype’,’s’) Chebyshev(设计阶数为n，截止频率为Wn的滤波器，通带内的波纹由 Rp（分贝）确定。Ftype指滤波器的类型，‘high’是高通，’stop’是带阻，无此参数则是低通，’s’指模拟域，无则表示数字域，b,a是对应变换的分子分母多项式的系数。 设计阶数为n，截止频率为Wn的滤波器，切比雪夫)II阻带内的波纹由 Rs（分贝）确定。ftype指型模拟和数滤波器的类型，‘high’是高通，’stop’是带阻，字滤波器设计（阻带等波纹）。 W=boxcar(n) W=hamming(n) W=hanning(n) 无此参数则是低通，’s’指模拟域，无则表示数字域，b,a是对应变换的分子分母多项式的系数。 产生矩形窗 产生一长度为n的矩形窗函数。 产生哈明窗 产生一长度为n的哈明窗函数。 产生汉宁窗。 产生一长度为n的汉宁窗。 W=blackman(n) 产生布莱克曼窗。 产生一长度为n的布莱克曼窗。 W=kaiser(n,beta) 产生凯泽而窗 Beta为影响窗函数旁瓣的参数，其与阻带衰减As的关系可参考书本。 z,p,k为s域传递函数的零点、极点和增益，Fs为取样频率,zd,pd,kd为经双线性变换后z域传递函数的零点、极点和增益。 [zd,pd,kd]=brilinear(z,p,k,Fs) 双线性变换。 [bz,az]=impinvar(b,a,F冲击响应不s) 变法实现模拟到数字的滤波器变换。 将模拟滤波器(b,a)变换成数字滤波器(bz,az)，两者的冲击响应不变。