数字信号处理实验报告汇总.doc

物理与电子电气工程学院

实验报告

课程名称：
数字信号处理
院 系：
物电学院
专 业：
电子信息科学与技术
班 级：

学 号：

姓 名：

实验报告（1）
实验名称
常见离散信号产生与实现
实验日期
2016年9月13日
指导教师
曹凤莲

一、实验目的
1. 掌握MATLAB语言的基本操作，学习基本的编程功能。
2. 掌握几种基本的离散时间信号（包括单位抽样序列，单位阶跃序列，单频正弦序列，单频复指数序列，实指数序列等）。
3.能够熟练利用MATLAB产生这些基本的离散时间信号。
二、预习要求
1．熟悉各种常用序列，掌握序列matlab实现。
2．预习MATLAB中序列产生的调用函数及绘制图形函数stem()。
三、实验内容
利用MATLAB编程产生和绘制下列有限长序列：
（1）单位脉冲序列
（2）单位阶跃序列
（3）矩形序列
（4）正弦型序列
（5）任意序列

四、实验报告
1．实现上述各序列。
(1) 单位冲激序列
function y=impDT(n) %定义impDT函数
y=(n==0);
n=-3:3;
x=impDT(n);
stem(n,x,'fill'); %离散时间信号的波形绘制在 MATLAB中一般用 stem 函数
xlabel('n');
grid on;
title('单位冲激序列 姓名 学号')
axis([-3 3 -0.1 1.1])

(2)单位阶跃序列
function y=uDT(n)
y=(n>=0);
n=-3:5;
x=uDT(n);
stem(n,x,'fill');
xlabel('n');
grid on
title('单位阶跃序列 姓名 学号')
axis([-3 5 -0.1 1.1])

（3）矩形序列
function y=uDT(n)
y=(n>=0);
n=-5:10;
x=uDT(n)-uDT(n-8); %此处调用了uDT.m 函数
stem(n,x,'fill');
xlabel('n');
grid on
axis([-5 10 -0.1 1.1])
title('矩形序列 姓名 学号')

（4）正弦型序列
n=-20:20;
x=sin(pi/5*n+pi/3);
stem(n,x,'fill');
xlabel('n');
grid on
title('正弦序列 姓名 学号')
axis([-20 20 -1.5 1.5])

（5）任意序列1：
function [x,n]=impseq (n0,n1,n2)
if((n0<n1)|(n0>n2)|(n1>n2))
error('参数必须满足n1<=n0<=n2')
end
n=[n1:n2];
x=[(n-n0)==0]; %以上定义impseq.m函数
n=-10:10;
x=impseq(0,-10,10)+2*impseq(1,-10,10)+3*impseq(2,-10,10)+4*impseq(3,-10,10)+5*impseq(4,-10,10);
stem(n,x,'fill');
xlabel('n');
grid on
title('任意序列1 姓名 学号')

任意序列2：
n=-10:10;
h=impseq(0,-10,10)+2*impseq(1,-10,10)+impseq(2,-10,10)+2*impseq(3,-10,10);
%同样的调用5.1中的impseq.m函数
stem(n,h,'fill');
xlabel('n');
ylabel('h');
grid on
title('任意序列2 姓名 学号')
axis([-10 10 -1 3])

2. 画出各序列的图形，并对结果进行分析。

分析总结：

Matlab的函数定义要放在M文件中，并且文件名与函数名一致才能运行。

此次实验大致巩固了Mtalab基本的语法与表达式。同时，我了解了单位阶跃序列、单位冲激序列、矩形序列、正弦序列的实现方法，并且通过他们的绘图了解这函数。从代码到图像，具体形象地理解了这些信号的特征。

实验报告（2）
实验名称
离散时间系统的时域分析
实验日期
2016年9月20日
指导教师
曹凤莲

一、实验目的
1．学会运用MATLAB求解离散时间系统的零状态响应。
2．学会运用MATLAB求解离散时间系统的单位取样响应。
3．学会运用MATLAB求解离散时间系统的卷积和。
二、预习要求
1. 预习实验中基础知识，熟悉MATLAB指令及y=filter(b,a,x)、impz(b,a,N)、y=conv(x,h)函数。
2. 结合实验内容，提前编制相应的程序。
3．思考改变差分方程的形式，单位抽样响应将如何变化。
三、实验内容
1. 试用MATLAB命令求解以下离散时间系统的单位取样响应。
（1）
（2）
2. 已知某系统的单位取样响应为，试用MATLAB求当激励信号为时，系统的零状态响应。
四、实验报告
1．根据内容要求，写出调试好的MATLAB语言程序及对应的MATLAB运算结果。
1. 试用MATLAB命令求解以下离散时间系统的单位取样响应。
%在m文件中定义两个函数以便调用
function y=impDT(n);
y=(n==0);
function y=uDT(n);
y=(n>=0
（1）
a=[3 4 1];
b=[1 1];
n=0:10;
x=impDT(n);
h=filter(b,a,x);
stem(n,h,'fill'),grid on;
xlabel('n');
title('系统单位取样响应 姓名 学号')

（2）
a=[2.5 6 10];
b=[1];
n=0:50;
x=impDT(n);
h=filter(b,a,x);
stem(n,h,'fill'),grid on;
xlabel('n');
title('系统单位取样响应 姓名 学号')

2. 已知某系统的单位取样响应为，试用MATLAB求当激励信号为时，系统的零状态响应。
nh=-10:20;
nx=-5:5;
x=uDT(nx)-uDT(nx-5);
h=(7./8).^nh.*(uDT(nh)-uDT(nh-10));
y=conv(x,h);
ny1=nx(1)+nh(1);
ny=ny1+(0:(length(nx)+length(nh)-2));
subplot(311);
stem(nx,x,'fill');
grid on;
xlabel('n');
title('x(n) 姓名 学号');
axis([-5 11 -0.1 1.1]);
subplot(312);
stem(nh,h,'fill');
grid on;
xlabel('n');
title('h(n) 姓名 学号');
axis([-5 11 -0.1 1.1]);
subplot(313);
stem(ny,y,'fill');
grid on;
xlabel('n');
title('y(n)=x(n)*h(n) 姓名 学号');
axis([-5 18 -1 5]);

分析总结：
函数的调用。编辑好M文件后调用单位冲激序列和单位阶跃序列,并利用这些简单序列的调用描述其他的序列，从而简化工作。

离散时间信号的波形绘制在 MATLAB 中一般用 stem 函数，由于 MATLAB 中矩阵元素的个数有限，所以MATLAB 只能表示一定时间范围内有限长度的序列；而对于无限序列，也只能在一定时间范围内表示出来。

实验报告（3）
实验名称
离散时间LTI系统的z域分析
实验日期
2016年9月27日
指导教师
曹凤莲

一、实验目的
1. 学会运用MATLAB分析离散时间系统的系统函数的零极点。
2. 学会运用MATLAB分析系统函数的零极点分布与其时域特性的关系。
3. 学会运用MATLAB进行离散时间系统的频率特性分析。
二、预习要求
1. 预习实验中基础知识，熟悉MATLAB指令及freqz（）、abs（）、angle（）函数。
2. 结合实验内容，提前编制相应的程序。
3．熟悉系统稳定的充要条件及学过的稳定判据。
三、实验内容
1. 试用MATLAB画出下列因果系统的系统函数零极点分布图，并判断系统的稳定性。
（1）
（2）
2. 试用MATLAB绘制系统的频率响应曲线。
四、实验报告
1．根据内容要求，写出调试好的MATLAB语言程序及对应的MATLAB运算结果。
1. 试用MATLAB画出下列因果系统的系统函数零极点分布图，并判断系统的稳定性。
（1）
>> B=[2,-1.6,-0.9];
>> A=[1,-2.5,1.96,-0.48];
>> zplane(B,A),grid on %通过调用 zplane(B,A) 绘出零极点
>> legend('零点','极点'); %标注名称与符号
>> title('零极点分布图 姓名 学号')

因果系统稳定条件是所有极点在单位圆内，系统有极点在单位圆外，故系统不稳定。

（2）
>> B=[1,-1];
>> A=[1,-0.9,-0.65,0.873,0];
>> zplane(B,A),grid on %通过调用 zplane(B,A) 绘出零极点
>> legend('零点','极点');
>> title('零极点分布图 姓名 学号');

因果系统稳定条件是所有极点在单位圆内，该因果系统所有极点全部在单位圆内，故系统稳定。
2. 试用MATLAB绘制系统的频率响应曲线。
A=[1];
B=[1,-3/4,1/8];
[H,w]=freqz(A,B,500,'whole'); %调用用于求离散系统频响特性的函数freqz()
Hm=abs(H);
Hp=angle(H); %angle用来求复数矩阵相位角的弧度值，其取值为-pi到pi
subplot(211);
plot(w,Hm);
grid on;
xlabel('\omega(rad)');
ylabel('Magnitude');
title('离散系统幅频特性曲线 学号 姓名');
subplot(212);
plot(w,Hp);
grid on;
xlabel('\omega(rad)');
ylabel('Phase');
title('离散系统相频特性曲线 学号 姓名')

2. 根据实验，思考freqz有哪些调用方式？频率响应有何特征？
MATLAB 提供了求离散时间系统频响特性的函数freqz，调用freqz 的格式主要有两种形式：
1.[H,w]=freqz(B,A,N)。其中，B 与A分别表示H (z)的分子和分母多项式的系数向量；N为正整数，默认值为512；返回值w包含[0,π]范围内N个频率等分点；返回值 H 则是离散时间系统频率响应在0-π范围内N个频率处的值。
2.[H,w]=freqz(B,A,N,’whole’)与第一种方式不同之处在于角频率的范围由[0,π]扩展到[0,2π]。
分析总结：
判断因果系统稳定性：如果极点全部在单位圆内，则系统是稳定的，反之则不稳定。
本次实验我掌握了调用roots或[z,p,k]=tf2zp(B,A)来求零极点的分布。在最后一个频率响应图像的横坐标单位应当为rad.我熟悉了系统函数的零极点分布、幅/相频特性曲线，同时也系统稳定判断有了更直观的理解。此次的教训是不能盲目地接受已有的实验例子，我们需要认识函数模型的物理意义，才能得出精确的结果。

实验报告（4）
实验名称
用FFT进行谱分析
实验日期
2016年10月10日
指导教师
曹凤莲

一、实验目的
掌握快速傅立叶变换的应用方法。
二、预习要求
1. 预习实验中的基础知识，运行编制好的MATLAB语句，熟悉fft()函数。
2. 根据实验中各的X(k)值以及频谱图，说明参数的变化对信号频谱产生哪些影响？
3. 思考如果周期信号的周期预先不知道，如何用FFT进行分析？
三、实验内容
1． 模拟信号，以进行采样，求：
（1）N＝40点FFT的幅度频谱，从图中能否观察出信号的2个频谱分量?
（2）提高采样点数，如N＝128，再求该信号的幅度频谱，此时幅度频谱发生了什么变化?信号的2个模拟频率和数字频率各多少？FFT频谱分析结果是否和理论一致？
N=40;
n=0:N-1;
t=0.01*n;
x=2*sin(4*pi*t)+5*cos(8*pi*t);
k=0:N/2;
w=2*pi/N*k;
X=fft(x,N); %快速傅里叶变换
magX=abs(X(1:N/2+1)); %绘制 x(n)的幅度谱
subplot(2,1,1);
stem(n,x,'.');
title('signal x(n) 姓名 学号');
subplot(2,1,2);
stem(w/pi,magX);
title('FFT N=40 姓名 学号');
xlabel('f(unit:pi)');
ylabel('|X|');
grid on;

N=128; %是2的幂数
n=0:N-1;
t=0.01*n;
x=2*sin(4*pi*t)+5*cos(8*pi*t);
k=0:N/2;
w=2*pi/N*k;
X=fft(x,N);
magX=abs(X(1:N/2+1));%绘制 x(n)的幅度谱
subplot(2,1,1);
stem(n,x,'.');
title('signal x(n) 姓名 学号');
subplot(2,1,2);
stem(w/pi,magX);
title('FFT N=128 姓名 学号');
xlabel('f(unit:pi)');ylabel('|X|');
grid on;

N=150;
n=0:N-1;
t=0.01*n;
x=2*sin(4*pi*t)+5*cos(8*pi*t);
k=0:N/2;
f=20/N*k;
X=fft(x,N);
magX=abs(X(1:N/2+1)); %绘制 x(n)的幅度谱
subplot(2,1,1);
stem(n,x,'.');
title('signal x(n) 姓名 学号');
subplot(2,1,2);
stem(f,magX);
title('FFT N=150 姓名 学号');
xlabel('f (unit:Hz)');ylabel('|X|');
grid on;

二. 利用MATLAB编程产生和绘制下列有限长序列并编写程序实现信号的谱分析。
1.
N=8;
n=0:N-1;
x=cos(n*pi/8);
k=0:N/2;
w=2*pi/N*k;
X=fft(x,N);
magX=abs(X(1:N/2+1));
subplot(2,1,1);stem(n,x,'.');
title('signal x(n) 姓名 学号');
subplot(2,1,2);stem(w/pi,magX);
title('FFT N=8 姓名 学号');
xlabel('f(unit:pi)');
ylabel('|X|');grid on;

N=16;
n=0:N-1;
x=cos(n*pi/8);
k=0:N/2;
w=2*pi/N*k;
X=fft(x,N);
magX=abs(X(1:N/2+1));
subplot(2,1,1);stem(n,x,'.');
title('signal x(n) 姓名 学号');
subplot(2,1,2);stem(w/pi,magX);
title('FFT N=16 姓名 学号');
xlabel('f(unit:pi)');
ylabel('|X|');
grid on;

N=32;
n=0:N-1;
x=cos(n*pi/8);
k=0:N/2;
w=2*pi/N*k;
X=fft(x,N);
magX=abs(X(1:N/2+1));%绘制 x(n)的幅度谱
subplot(2,1,1);
stem(n,x,'.');
title('signal x(n) 姓名 学号');
subplot(2,1,2);
stem(w/pi,magX);
title('FFT N=32 姓名 学号');
xlabel('f(unit:pi)');
ylabel('|X|');
grid

2.
N=10;
n=0:N-1;
x=cos(0.125*pi*n)+2.*cos(0.25*pi*n);
k=0:N/2;
w=2*pi/N*k;
X=fft(x,N);
magX=abs(X(1:N/2+1));%绘制 x(n)的幅度谱
subplot(2,1,1);
stem(n,x,'.');
title('signal x(n) 姓名 学号');
subplot(2,1,2);
stem(w/pi,magX);
title('FFT N=10 姓名 学号');
xlabel('f(unit:pi)');ylabel('|X|');
grid on;

N=16;
n=0:N-1;
x=cos(0.125*pi*n)+2.*cos(0.25*pi*n);
k=0:N/2;
w=2*pi/N*k;
X=fft(x,N);%绘制 x(n)的幅度谱
magX=abs(X(1:N/2+1));
subplot(2,1,1);stem(n,x,'.');
title('signal x(n) 姓名 学号');
subplot(2,1,2);stem(w/pi,magX);t
itle('FFT N=16 姓名 学号');
xlabel('f(unit:pi)');
ylabel('|X|');
grid on;

N=32;
n=0:N-1;
x=cos(0.125*pi*n)+2.*cos(0.25*pi*n);
k=0:N/2;
w=2*pi/N*k;
X=fft(x,N);
magX=abs(X(1:N/2+1));
subplot(2,1,1);stem(n,x,'.');
title('signal x(n) 姓名 学号');
subplot(2,1,2);stem(w/pi,magX);
title('FFT N=32 姓名 学号');
xlabel('f(unit:pi)');ylabel('|X|');
grid on;

N=64;
n=0:N-1;
x=cos(0.125*pi*n)+2.*cos(0.25*pi*n);
k=0:N/2;
w=2*pi/N*k;
X=fft(x,N);
magX=abs(X(1:N/2+1));
subplot(2,1,1);stem(n,x,'.');title('signal x(n) 姓名 学号');
subplot(2,1,2);stem(w/pi,magX);title('FFT N=64 姓名 学号');
xlabel('f(unit:pi)');ylabel('|X|');grid

四、实验报告
1．根据内容要求，写出调试好的MATLAB语言程序及对应的结果。
2. 结合实验中所得给定典型序列幅频特性曲线，与理论结合比较，并分析说明误差产生的原因以及用FFT作谱分析时有关参数的选择方法。
3. 总结实验所得主要结论。

分析总结：
周期信号的频谱是离散的，只有当整数倍周期的长度做FFT，得到的离散谱才能代表周期信号的频谱。
对于有限长序列做DFT时，首先要满足乃奎斯特抽样定理。其次，若截取不当则造成频谱泄露。抽样点N的选择要适当。如果太小，分辨率不够；对周期序列应当取其整数倍来抽样。因此N要适当得大一些。

误差产生的原因 ：
1.对周期序列的截取不当，造成频谱泄露；
2.抽样点数N太少，频率分辨率不够 用FFT做谱分析时参数的选择 ；
3.抽样频率要满足奈奎斯特准则，不小于信号最高频率的2倍；
4.在抽样频率一定的情况下，抽样点数N要适当。（太小会造成频率分辨力不够，太大会造成数据冗余。对周期序列，最好截取周期的整数倍进行谱分析。）

实验报告（5）
实验名称
实验五 数字滤波器结构的实现
实验日期
2016年10月17日
指导教师
曹凤莲

一、实验目的
（1） 加深对数字滤波器分类与结构的了解；
（2） 明确数字滤波器的基本结构及其相互间的转换方法；
（3） 掌握用MATLAB进行数字滤波器各种结构相互间转换的子函数及程序编写方法。
二、实验原理
一个离散LSI系统可用系统函数来表示；

也可用差分方程来表示：

当至少有一个不为0时，则在有限z平面上存在极点，表示一个IIR数字滤波器；当全都为0时，系统不存在极点，表示一个FIR系统。
IIR数字滤波器的基本结构分为直接Ⅰ型、直接Ⅱ型、级联型和并联型。
FIR数字滤波器的基本结构分为横截型、级联型、并联型、、线性相位型和频率抽样型。
三、实验仪器
微型计算机、MATLAB
四、实验内容
(1) 已知一个IIR系统的系统函数为

将其从直接型转换为级联型和并联型结构，并画出各种结构的流程图。
直接型:

级联型:
b=[0.1 -0.4 0.4 -0.1];
a=[1 0.3 0.55 0.2];
[sos,g]=tf2sos(b,a) %使用tf2sos将系统函数分解成一系列二阶子系统的级联形式 
[r,p,k]=residuez(b,a) %多项式分母若无重根的情况，应可进行因式分解，并求得其系数，此可以利用residue指令求解
sos =
1.0000 -2.6180 0 1.0000 0.3519 0
1.0000 -1.3820 0.3820 1.0000 -0.0519 0.5683

g =
0.1000

r =
-0.2893 + 0.0001i
-0.2893 - 0.0001i
1.1786

p =
0.0260 + 0.7534i
0.0260 - 0.7534i
-0.3519

k =
-0.5000

H（z）=0.1*(1-2.618Z-1)/(1+0.3519 Z-1 )* (1-1.382Z-1 + 0.382Z-2)/(1-0.0519 Z-1+ 0.5683Z-2 )

并联型：
b=[0.1 -0.4 0.4 -0.1];
a=[1 0.3 0.55 0.2];
[C,B,A]=tf2par(b,a)%使用tf2par将系统函数分解成并联形式
C =
-0.5000

B =
-0.5786 0.0148
1.1786 0

A =
1.0000 -0.0519 0.5683
1.0000 0.3519 0
H(z)=-0.5+(-0.5786 +0.0148Z-1)/(1-0.0519Z-1)+ 0.5683Z-2)+ 1.1786/(1+0.3519Z-1)

(2) 已知一个FIR系统的系统函数为

将其从横截型转换为级联型结构，并画出各种结构的流程图。
b=[0.2 0.88 0.212 0.212 0.885];
a=[1];
[sos,g]=tf2sos[b,a] %使用tf2sos将系统函数分解成一系列二阶子系统的级联形式 
sos =
1.0000 5.2595 4.6234 1.0000 0 0
1.0000 -0.8595 0.9571 1.0000 0 0

g =
0.2000

[b,a]=sos2tf(sos,g) %使用sos2tf反过来将二次分式转换为系统函数
b =
0.2000 0.8800 0.2120 0.2120 0.8850

a =
1 0 0 0

H(z)=0.2*(1+5.2595Z-1+4.6234Z-2)*(1-0.8595Z-1+0.9571Z-2)

分析总结：
(1) 直接I型。
缺点：
①需要2N个延迟器（z-1），太多；
②系数ai、bi对滤波器性能的控制不直接，对极、零点的控制难，一个ai、bi的改变会影响系统的零点或极点分布；
③对字长变化敏感（对ai、bi的准确度要求严格）；
④易不稳定，阶数高时，影响更大。
（2）直接II型(典范型)：
优点：延迟线减少一半，为N个，可节省寄存器或存储单元。
通常在实际中很少采用上述两种结构实现高阶系统，而是把高阶变成一系列不同组合的低阶系统（一、二阶）来实现。
（3）级联型
优点：
①简化实现，用一个二阶节，通过变换系数就可实现整个系统；
②极、零点可单独控制、调整，从而使频率响应调节更方便；
③各二阶节零、极点的搭配可互换位置，优化组合以减小运算误差；
④可流水线操作。
缺点：
阶节电平难控制，电平大易导致溢出，电平小则使信噪比减小。
（4） 并联型
优点：
①系统实现简单，只需一个二阶节，系统通过改变输入系数即可完成；
②极点位置可单独调整；
③运算速度快（可并行进行）；
④各二阶网络的误差互不影响，总的误差小，对字长要求低。
缺点：
不能直接调整零点，因多个二阶节的零点并不是整个系统函数的零点，当需要准确的传输零点时，级联型最合适。

实验报告（6）
实验名称
实验六 IIR数字滤波器的设计
实验日期
2016年10月25日
指导教师
曹凤莲

一、实验目的
1. 要求掌握IIR数字滤波器的设计原理、设计方法和设计步骤。
2. 能根据给定的滤波器指标进行滤波器设计。
3. 掌握数字巴特沃斯滤波器、数字切比雪夫滤波器设计原理和步骤。
二、预习要求
1. 预习实验中的基础知识，运行编制好的MATLAB语句，熟悉函数bilinear，impinvar，butter，buttord。
2. 思考双线性变换法中Ω和ω之间的关系是怎样的。
3. 能否利用公式完成脉冲响应不变法的数字滤波器设计？为什么？
三、实验内容
1、基于Butterworth型模拟滤波器原型使用冲激不变转换方法设计数字滤波器，要求具有下面的参数指标：
通带截止频率： 通带波动值：
阻带截止频率： 阻带波动值：
clear all;
Rp=1;
As=15;
wp=0.2*pi;
ws=0.3*pi;
T=1;
wap=wp/T;
was=ws/T; %得到模拟滤波器的频率—采用脉冲响应不变法的频率转换形式
[N,wn]=buttord(wap,was,Rp,As,'s'); %计算模拟滤波器的最小阶数
wn=wap./(10.^(0.1.*Rp)-1).^(1./(2.*N));
[bs,as]=butter(N,wn,'s'); %设计出所需的模拟低通滤波器
[bz,az]=impinvar(bs,as,1/T); %冲激不变法AF到DF
Rip=10^(-Rp/20);
Atn=10^(-As/20);
[H,w]=freqz(bz,az,512);
mag=abs(H);
db=20*log10(mag/max(mag));
figure;
subplot(2,1,1);
plot(w/pi,mag);
title('幅频特性 姓名 学号');
xlabel('w(/pi)');
ylabel('|H(jw)|');
axis([0,1,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 0.5 1]);
set(gca,'YTickMode','manual','YTick',[0 Atn Rip 1]);
grid on;
subplot(2,1,2);
plot(w/pi,db);
title('幅频特性(db) 姓名 学号');
xlabel('w(/pi)');
ylabel('dB');
axis([0,1,-40,5]);
set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 0.5 1]);
set(gca,'YTickMode','manual','YTick',[-40 -As -Rp 0]);
grid on;

2、基于Butterworth型模拟滤波器原型使用双线性不变法设计数字滤波器，要求具有下面的参数指标：
通带截止频率： 通带波动值：
阻带截止频率： 阻带波动值：
clear all;
Rp=1; %通带最大衰减
As=15; %阻带最小衰减
wp=0.2*pi; %通带截止频率
ws=0.3*pi; %阻带截止频率
T=1;
wap=2/T*tan(wp/2);
was=2/T*tan(ws/2);
[N,wn]=buttord(wap,was,Rp,As,'s'); %计算模拟滤波器的最小阶数
[bs,as]=butter(N,wn,'s'); %设计出所需的模拟低通滤波器
[bz,az]=bilinear(bs,as,1/T); %双线性变换法
Rip=10^(-Rp/20);
Atn=10^(-As/20);
[H,w]=freqz(bz,az,512);
mag=abs(H);
db=20*log10(mag/max(mag));
figure;
subplot(2,1,1);
plot(w/pi,mag);
title('幅频特性 姓名 学号');
xlabel('w(/pi)');
ylabel('|H(jw)|');
axis([0,1,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 0.5 1]);
set(gca,'YTickMode','manual','YTick',[0 Atn Rip 1]);
grid on;
subplot(2,1,2);
plot(w/pi,db);
title('幅频特性(db) 姓名 学号');
xlabel('w(/pi)');
ylabel('dB');
axis([0,1,-40,5]);
set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 0.5 1]);
set(gca,'YTickMode','manual','YTick',[-40 -As -Rp 0]);
grid on;

四、实验报告要求
1．按照实验步骤及要求，比较各种情况下的滤波性能。
2. 总结实验所得主要结论。

分析总结：

此次实验通过冲激响应法及双线性变换法设计出巴特沃斯低通滤波器，分别利用impinvar（）函数及bilinear（）函数进行设计。
主要过程为根据数字滤波器性能指标计算出相应的模拟滤波器的设计指标，然后再设计出相应的模拟滤波器的系统函数H(s)，最后根据变换得到相应的数字滤波器。 
其次是两种变换法优缺点的比较。脉冲响应法最大缺点是有频率响应的混叠效应，所以脉冲响应法只适用于限带的模拟滤波器，而且高频衰减越快，混叠效应越小。双线性变换法主要优点是避免频率响应的混叠现象，这是因为S平面与Z平面是单值一一对应关系。