数字信号分析实验.doc

桂林工学院电子与计算机科学系
目录

目录 1
前言 信号处理实验简介 2
实验一 离散信号的时域运算与变换 4
实验二 因果离散线性系统的时域分析 9
实验三 DFT变换的性质及应用 13
实验四 数字低通巴特沃斯滤波器的设计 18
实验五 窗函数设计FIR滤波器 23
实验六 基于FFT的图象压缩与放大 26
附 录 31
附录1 ：MATLAB编程及绘图方法 31
附录2 ：MATLAB矩阵及矩阵的运算 38
前言 信号处理实验简介

对于一个信号处理系统来说，可以将流程分为三个阶段，首先是定义输入序列，第二是对输入的序列进行运算，第三是保存输出序列，并显示结果。

一、实验课程任务与要求
1. 实验目的

信号处理实验教学是为了将学生的计算机操作能力、分析能力、软件设计能力与应用实践结合起来，引导学生由浅入深地掌握信号处理理论与开发工具，具备实际应用的信号处理软件开发与制作基础。

2. 实验基本要求（以软件设计为主要表现形式）

（1）上机前应准备好实验的程序设计算法描述与关键分析内容；
（2）准备好程序测试数据和设备操作步骤，上机调试、运行；
（3）完成每个实验后进行数据与程序对比分析；
（4）递交实验结果的可执行程序、源程序并演示实验结果；
（5）写出实验报告，
二、实验学时安排
1．实验一：离散信号的时域运算与变换 （2学时）
2．实验二：因果离散线性系统的时域分析 （2学时）
3．实验三：DFT变换的性质及应用 （2学时）
4．实验四：数字低通巴特沃斯滤波器的设计 （2学时）
5. 实验五： 窗函数设计FIR滤波器 （2学时）
6. 实验六：基于FFT的图象压缩与放大 （2学时）
三、实验报告格式
实验报告

姓名：________________学号：_______________实验日期：_______________
实验题目：
实验目的：
实验内容：
实验地点：
实验结果：（包括列出实验编写的所有文件及各项实验结果的曲线，并加注必要的说明）
结果分析：
总结：
四、实验考核
（1）实验预习报告；
（2）实验签到；
（3）上机实际操作；
（4）实验设计报告；
五、实验仪器设备要求
（1）有快速的较高性能微机和较大内存与硬盘的设备；
（2）设备数量能适应学生人数；
（3）有Matlab程序设计环境；
六、教材及参考书
1．张志勇等. 精通MATLAB6.5. 北京航空航天大学出版社. 2003年3月
2. 邹鲲等.MATLAB6.x信号处理. 清华大学出版社. 2002年5月
3．陈怀琛等. MATLAB及在电子信息课程中的应用.电子工业出版社. 2002
4．程佩清.数字信号处理教程 [M].清华大学出版社.2003
实验一 离散信号的时域运算与变换
实验目的

1. 熟悉MATLAB编程特点
2. 了解离散序列的延迟、相加、相乘及平移、反折、及倒相变换
实验内容

1.设计一个实现序列移位的函数
将序列x(n)＝｛1,2,3,4,0,7｝，其中n＝0：5的每一个样本都移动3个周期，移位后的序列y(n)=x(n-3)

2.序列的奇偶分解
将序列x(n)＝｛0,1,2,3,4,3,2,1,0｝，其中n＝-3：5进行奇偶分解

用函数stem显示其奇偶序列
3.序列的加法运算

设x1(n)＝｛1,0.5,0.3,0.4｝其中n=-1：2 ；x2(n)＝｛0.2,0.3,0.4,0.5,0.8,1｝其中n=-2：3
4.实现2中序列的翻转

实验报告

1.列出本次实验编写的所有文件及各项实验结果的曲线，加注必要的说明
2.对本实验结果做理论计算，解释实验结果
3.总结实验体会及实验中存在的问题
思考：

1. 在matlab中应如何实现序列的相乘运算
2. 在matlab中应如何实现序列的倒相运算？
1. 实验原理：

熟悉MATLAB

首先，在用MATLAB表示离散序列并将其可视化时，我们还应注意以下几点：
第一、离散时间序列无法用符号运算来表示，要用矩阵的形式；
第二、由于在MATLAB 中矩阵的元素个数有限，因此，MATLAB无法表示无限长的序列；
第三、在绘制离散信号波形的函数stem命令，而不是plot命令。
下面是MATLAB绘制单位抽样序列及其移位序列的函数，n0是单位抽样序列的位移量，n1,n2是序列的起止时刻， ，调用该函数就可以绘出单位抽样序列及其移位序列的波形图。
function [x,n] = impuls (n0,n1,n2)
% Generates x(n) = delta(n-n0); n=n0 处建立一个单位抽样序列
% [x,n] = impuls (n0,n1,n2)
if ((n0 < n1) | (n0 > n2) | (n1 > n2))
error('arguments must satisfy n1 <= n0 <= n2')
end
n = [n1:n2];
x = [zeros(1,(n0-n1)), 1, zeros(1,(n2-n0))];

将上述文件存为：impuls.m，在命令窗口输入
>>n0=1,n1=-20,n2=30;
[x,n]=impuls (n0,n1,n2);figure（1），stem(n,x)

上图为右移了一个单位的单位抽样序列
其次，对于离散序列来说，序列相加、相乘是将两序列对应时间序号的值逐项相加或相乘，平移、反折、及倒相变换与连续信号的定义完全相同，但需要注意，与连续信号不同的是，在MATLAB中，离散序列的时域运算和变换不能用符号运算来实现，必须用向量表示的方法，即在MATLAB中离散序列的相加、相乘需表示成两个向量的相加、相乘，因而参加运算的两序列向量必须具有相同的维数。
• 在MATLAB中，用一个向量即可表示一个有限长度的序列。但是，这样的向量并没有包含对应的时间序号的信息。所以，要完整地表示离散信号要用两个向量

• 如： 序列: f(k)={1,2,-1,3,4,5}
• k=0
• 在MATLAB中应表示为：
• f=[1,2,-1,3,4,5]
• k=[-3,-2,-1,0,1,2]或是k＝-3：2
2. MATLAB常用信号生成函数：

熟悉MATLAB

• ★ ZEROS
• 功能：产生全零阵列
• 调用格式：X=ZEROS (N) %产生N行N列的全零矩阵
• X=ZEROS(M,N) %产生M行N列的全零矩阵
• ★    ONES
• 功能：产生全1阵列
• 调用格式：X=ONES(N) %产生N行N列的全1矩阵
• X=ONES(M,N) %产生M行N列的全1矩阵
• ★    SINC
• 功能：辛格函数
• 调用格式：Y= SINC(X) %
• ★    RECTPULS
• 功能：产生矩形脉冲信号
• 调用格式：Y= RECTPULS (T) %产生高度为1、宽度为1、关于T＝0对称的矩形脉冲
• Y= RECTPULS (T，W) %产生高度为1、宽度为W、关于T＝0对称的矩形脉冲
• ★    RAND
• 功能：产生伪随机序列
• 调用格式：Y= RAND (1,N) %产生[0,1]上均匀分布的随机序列
• Y= RANDN (1,N) %产生均值为0，方差为1的白噪声序列
• ★    SAWTOOTH
• 功能：产生周期锯齿波或三角波
• 调用格式：Y= SAWTOOTH (T) %产生幅值为＋1，－1，以2为周期的方波
• Y=SAWTOOTH(T,WIDTH)% 产生幅值为＋1，－1，以WIDTH *2为周期的方波
• ★    SQUARE
• 功能：产生方波
• 调用格式：Y= SQUARE (T) %产生幅值为＋1，－1，以2为周期的锯齿波
• Y= SQUARE(T,DUTY) % 产生幅值为＋1，－1，以占空比为DUTY的方波
• 例： t = 0:.0001:.0625;
• y = SQUARE(2*pi*30*t，80); plot(t,y)％产生一个占空比为80%的方波
• ★    FLIPLR
• 功能： 序列左右翻转
• 调用格式：Y = FLIPLR(X)
• % X = 1 2 3 翻转后 3 2 1
• 4 5 6 6 5 4
• ★    CUMSUM 、SUM
• 功能： 计算序列累加
• 调用格式：Y = CUMSUM(X) % 向量X元素累加，记录每一次的累加结果,而SUM只记录最后的结果
3. 参考程序

（1）实现序列移位的函数
• function [y,n] = segshift(x,n,n0)
• % 功能 y(n) = x(n-n0)
• % [y,n] = segshift(x,n,n0)
• n = n+n0; y = x;

（2）实现两序列相加
• function [y,n] = sigadd(x1,n1,x2,n2)
• % 实现 y(n) = x1(n)+x2(n)
• % -----------------------------
• % [y,n] = sigadd(x1,n1,x2,n2)
• % y = 在包含n1 和 n2 的n点上求序列和,
• % x1 = 在 n1上的第一序列
• % x2 = 在 n2上的第二序列(n2可与 n1不等)
• n = min(min(n1),min(n2)):max(max(n1),max(n2)); % y(n)的长度
• y1 = zeros(1,length(n)); y2 = y1; % 初始化
• y1(find((n>=min(n1))&(n<=max(n1))==1))=x1; % 具有y的长度的 x1
• y2(find((n>=min(n2))&(n<=max(n2))==1))=x2; % 具有y的长度的x2
• y = y1+y2; % 序列相加.

（3）实现序列奇偶分解
• m = -fliplr(n);
• m1 = min([m,n]); m2 = max([m,n]); m = m1:m2;
• x1 =[ zeros(1,(length(m)-length(n))),x];
• xe = 0.5*(x1 + fliplr(x1))
• xo = 0.5*(x1 - fliplr(x1))
• subplot(1,2,1),stem(m,xe),ylabel('xe'),
• subplot(1,2,2),stem(m,xo),ylabel('xo'),

奇偶分解结果入下图所示：

（4）正弦序列翻转
n=-1:10;
x=sin(0.4*pi*n);
y=fliplr(x);
n1=-fliplr(n);
subplot(2,1,1),stem(n,x)
subplot(2,1,2),stem(n1,y)
效果如右图所示：
42 / 43
实验二 因果离散线性系统的时域分析
实验目的

实现由差分方程构成的数字滤波器：

实验内容

1.用两种不同的程序计算数字滤波器的单位抽样响应，
给定差分方程：y(n)-y(n-1)+0.9y(n-2)=x(n)
a.编制文件himpz.m，实现数字滤波器的单位抽样响应
b. 编制函数文件hfilter.m，实现数字滤波器的单位抽样响应
2.给定计算此数字滤波器的单位阶跃响应，并判断系统的稳定性。用stem(n,y)画出相应的图形。参考程序gfilter.m
实验报告

1.列出本次实验编写的所有文件及各项实验结果的曲线，加注必要的说明
2.对本实验结果做理论计算，解释实验结果
3.总结实验体会及实验中存在的问题
思考：

1.在编写程序gfilter.m时，调用himpz.m和hfilter.m哪一个比较合理？
2.实验中给定的滤波器是什么性质的滤波器？
1.实验原理：

熟悉MATLAB

LTI离散系统的单位抽样响应定义为：当输入为单位抽样序列是系统产生的零状态响应，用h(n)表示。
LTI离散系统

(1)MATLAB提供了函数impz( )求离散系统的单位抽样响应，并绘制时域波形。
★ IMPZ
功能：函数将绘出由矢量a和b表示的数字滤波器在指定取样时刻范围内的抽样响应h(n)的时域波形，并求其数值解。
调用格式：H= IMPZ (B,A)%计算离散系统冲激响应，取样点数n由函数自动选取
或H=IMPZ(B,A,N)％计算指定时刻范围内（0:N-1）的离散系统抽样响应
说明：由矢量a和b构成的数字滤波器的差分方程为

其中： B=[b0,b1,…,bM,bM-1],A=[a0,a1,…,aN,aN-1]，要求A的首项不为0
(2)由于系统的单位冲激响应h(n)也就是系统输入为 时系统的零状态响应，除了用上述的impz求解外。还可以调用filter函数求h(n)，此时系统的输入为单位抽样序列。
离散系统差分方程求解的 filter函数
调用格式（1）y=filter（b,a,x）
说明：b,a是式

的系数组成的向量b=[b0,b1,b2,…bM-1]及a=[a0,a1a2,…aN-1]，x是输入向量数组，y是输出向量数组和x的长度相同，而且系数a0要保证不为零。
调用格式（2）y=filter（b,a,x,xic）
说明：xic是初始条件等效的输入序列，MATLAB提供了函数filtic（）来确定xic
xic调用格式 xic=filtic（b,a,x,y）
说明：b,a是差分方程的系数数组，y,x是初始条件数组，分别由y(n),x(n)的初始条件确定，即y=[y(-1),y(-2),….]及x=[x(-1),x(-2),…]
(3)要判断系统稳定性，一个方法是必须对h(n)求累加，如果结果不为inf(无穷大)，则系统稳定；另一个方法是求系统的极点（即差分方程等号右端为零时方程的解），判断极点是否都在单位园内，即模值是否小于1。MATLAB提供了一个函数ROOTS（C）求齐次方程的根
★ 函数ROOTS
功能：求多项式的根
调用格式：X＝ROOTS(C) %计算形如＝0的根，
％其中多项式系数矩阵为C， C有N+1个元素，
2.参考程序：

熟悉MATLAB

(1)himpz.m
a=[1,-1,0.9];b=1;
impz(b,a)
title('单位抽样响应')

(2)hfilter.m
a=[1,-1,0.9];b=1;
% 求h(n)
x=impuls (0,-10,100);n=[-10:100];
h=filter(b,a,x);
subplot(2,1,1);stem(n,h)
axis([-10,100,-1.1,1.1])
title(‘单位冲激响应’);xlabel(‘n’);ylabel(‘h(n)’)

（3）gfilter.m
% 求单位阶跃响应
a=[1,-1,0.9];b=1;
x=[zeros(1,10),ones(1,101)];
n=-10:100;
s=filter(b,a,x);
subplot(2,1,2);stem(n,s)
axis([-10,100,-.5,2.5])
title(‘单位阶跃响应’);xlabel(‘n’);ylabel(‘s(n)’)
%判断系统稳定
sum(abs(h))
z=roots(a);
magz=abs(z)
程序运行结果：
ans = 14.8337
magz =
0.9487
h(n)累加的结果是一个实数，说明系统稳定；magz说明极点的模小于1，所以系统稳定。

实验三 DFT变换的性质及应用
实验目的

1.实现信号的DFT变换
2. 了解DFT 应用：
(1)用DFT计算卷积
(2)用DFT 对序列进行谱分析
实验内容

1.用三种不同的DFT程序计算的傅立叶变换X(k)，并比较三种程序的计算机运行时间
步骤：
a.用for循环语句编制函数文件dft1.m，实现循环计算X(k)；
b.编写矩阵运算的函数文件dft2.m，实现矩阵计算X(k)；
根据定义： 为序列x(n)的DFT
则

c.调用FFT函数直接计算X(K)，如程序函数dft3.m
d.分别利用上述三种不同的方式编写的DFT程序计算序列x(n)的DFT变换X(k)，并画出幅频和相频特性，并比较3个程序的运行时间。
2.给定，利用DFT实现两序列的线性卷积运算，并研究DFT的点数与混叠的关系，并用stem(n,y)画出相应的图形。参考程序dft4.m

3. 讨论序列补零及增加数据长度对信号频谱的影响
(1)求出序列x(n)=cos(0.48 n)+cos(0.52 n)基于有限个样点n=10的频谱
(2) 求n=100时，取x(n)的前10个，后90个设为零，得到x(n)的频谱
(3) 增加x(n)有效的样点数，取100个样点得到x(n)的频谱
实验报告

参考程序dft5.m

1.列出本次实验编写的所有文件及各项实验结果的曲线，加注必要的说明
2.对本实验结果做理论计算，解释实验结果
3.总结实验体会及实验中存在的问题
思考：

1.直接计算N=256电DFT、FFT两者在理论上的计算量有何差别？
2.补零对x(n)计算出的谱线有什么影响，而增加有效数据长度对x(n)计算出的谱线有什么影响，这两种方式说明了什么问题？

参考程序：

熟悉MATLAB

1.dft.m
function [am,pha]=dft1(x)
N=length(x);
w=exp(-j*2*pi/N);
for k=1:N
sum=0;
for n=1:N
sum=sum+x(n)*w^((k-1)*(n-1));
end
am(k)=abs(sum);
pha(k)=angle(sum);
end

2.dft2.m
function [am,pha]=dft2(x)
N=length(x);
n=[0:N-1];
k=[0:N-1];
w=exp(-j*2*pi/N);

nk=n’*k;
wnk=w.^(nk);
Xk=x*wnk;
am= abs(Xk);
pha=angle(Xk);

3.dft3.m
function [amfft,phafft]=dft3(x)
N=length(x);
Xk=fft(x);
amfft= abs(Xk);
phafft=angle(Xk);
4. dft4.m
%N1+N2-1=23<32
N=32;
x=[0:15];
xx=[x,zeros(1,16) ];
h=[ones(1,8),zeros(1,24)];
Xk=fft(xx,N);
Hk=fft(h,N);
Yk=Xk.* Hk;
y=ifft(Yk,N);
n=0:N-1;
stem(n,y);
hold on
%N=N1=16
N1=16;
x1=[0:15];
h1=[ones(1,8),zeros(1,8)];
Xk1=fft(x1,N1);
Hk1=fft(h1,N1);
Yk1=Xk1.* Hk1;
y1=ifft(Yk1,N1);
n1=0:N1-1;
stem(n1,y1,’.’,’m’);
5.dft5.m
％x(n)基于10个样点的频谱
figure(1)
n=[0:1:99];
x=cos(0.48*pi*n)+cos(0.52*pi*n);
n1=[0:1:9];y1=x(1:1:10);
subplot(2,1,1);stem(n1,y1);title('signal x(n), 0 <= n <= 9');xlabel('n')
axis([0,10,-2.5,2.5])
Y1=fft(y1);magY1=abs(Y1(1:1:6));
k1=0:1:5;w1=2*pi/10*k1;
subplot(2,1,2);stem(w1/pi,magY1);title(‘10点DFT');
xlabel('w/pi')，axis([0,1,0,10])
%在10个样点的基础上添90个零，得到密度高的频谱
figure(2)
n3=[0:1:99];y3=[x(1:1:10) zeros(1,90)]; ％添90个零。得到100个数据
subplot(2,1,1);stem(n3,y3);title('signal x(n), 0 <= n <= 9 + 90 zeros');xlabel('n')
axis([0,100,-2.5,2.5])
Y3=fft(y3);magY3=abs(Y3(1:1:51));
k3=0:1:50;w3=2*pi/100*k3;
subplot(2,1,2);stem(w3/pi,magY3);title('100点DFT');xlabel('w/pi')
axis([0,1,0,10])
%增加x(n)有效的样点数，取100个样点
figure(3)
n=[0:1:99];
x=cos(0.48*pi*n)+cos(0.52*pi*n);
subplot(2,1,1);stem(n,x);title('signal x(n), 0 <= n <= 99');xlabel('n')
axis([0,100,-2.5,2.5])
X=fft(x);magX=abs(X(1:1:51));
k=0:1:50;w=2*pi/100*k;
subplot(2,1,2);stem(w/pi,magX);title('100点DFT);xlabel('w/pi')
axis([0,1,0,60])

figure（1）

figure（2）

figure（3）

实验四 数字低通巴特沃斯滤波器的设计
实验目的

1.掌握IIR数字滤波器的设计方法
实验原理

IIR数字滤波器的设计，主要采用间接法，即：首先设计出低通模拟滤波器H(S)；进行频率变换，将其转换为高通、带通、带阻滤波器；再用脉冲响应不变法或双线性变换法从模拟滤波器转换为数字滤波器。
对单极点的N阶H(S)用部分展开式：

冲激响应不变法取H(S)的单阶极点Sk的指数函数作为H（Z）的极点

双线性变换法是用代换H（S）中的S得到H(Z),双线性变换法可完全消除频率混叠失真但存在非线性频率失真，而冲激响应不变法存在混叠失真。
在不同的设计阶段MATLAB的信号处理工具箱都给出了相应的滤波器设计函数，这些函数代表了不同类型的逼近函数的滤波器，常用的有巴特沃斯滤波器、切比雪夫滤波器、椭圆滤波器类型。
设计流程图

下图所示为IIR滤波器设计过程及所用到的MATLAB工具箱中的函数。
模拟原型低通滤波器设计（buttap，cheb1ap，cheb2ap，bess1lap，ellipap函数）
求最小阶数N（buttord，cheb1ord，cheb2ord，ellipord函数）
频率转换（低通转低通lp2lp，低通转高通lp2hp，低通转带通lp2bp， 低通转带阻函数lp2bs）
模拟滤波器转数字滤波器（bilinear，impinvar函数）
合成一步的设计函数
（butter，cheb1，cheb2，ellip，besself）

滤波器系数
A,B

实验内容

1.设计一个低通巴特沃斯模拟滤波器，绘制滤波器的的幅频响应及零极点图。指标如下：
通带截止频率：WP＝1000HZ, 通带最大衰减：RP=3dB
阻带截止频率：Ws＝2000HZ, 阻带最小衰减：Rs=40 dB
参考程序butter1.m
2. 用冲激响应不变法和双线性变换法将一模拟低通滤波器转换为数字低通滤波器
并图释H(S)和H(Z)，采样频率Fs＝1000Hz
参考程序butter2.m
实验报告

1.列出本次实验编写的所有文件及各项实验结果的曲线，加注必要的说明
频率特性应包括幅频、相频特性
2.理论计算模拟低通原型滤波器的阶数和极点分布，说明实验中所用的冲激响应不变法和双线性变换法的原理及其编程根据
3.总结实验体会及实验中存在的问题
思考：

1.双线性变换法和冲激响应不变法比较有何优点？

IIR设计函数：

熟悉MATLAB

巴特沃斯滤波器为例介绍设计函数的功能和用法，其它设计函数的用法类似。
（1）求最小阶数N的函数buttord
调用格式1：[N, Wn] = BUTTORD(Wp, Ws, Rp, Rs, 's')
² 功能：求出巴特沃斯模拟滤波器的最小阶数N和频率参数Wn
² 说明：Wp, Ws,为通带截止频率和阻带截止频率如式所示；
，
Rp, Rs为通带最大衰减 Apass(dB)和阻带最小衰减Astop(dB)如式

's'对应模拟滤波器。
（2）模拟低通滤波器设计函数buttap
调用格式：[Z,P,K] = BUTTAP(N)
² 功能：按给定的阶数N设计出巴特沃斯模拟低通滤波器
² 说明：Z,P,K为返回的N阶模拟滤波器的零点、极点和增益系数。如：
[式3.5.4]
其中k为增益常数，zj和pi为系统函数H(Z)的N个零点和N个极点。
(3)模拟数字变换函数－双线性变换函数bilinear或脉冲响应不变法函数impinvar
调用格式：[NUMd,DENd] = BILINEAR(NUM,DEN,Fs)
² 功能： 把模拟滤波器系数为NUM,DEN变为近似等价的数字滤波器系数为NUMd,DENd
调用格式：[BZ,AZ] = IMPINVAR(B,A,Fs)
² 功能：设计出数字滤波器其单位抽样响应h(n)为模拟滤波器的单位冲激响应h(t)的采样值，采样频率为Fs
² 说明：B,A为模拟滤波器的系统函数H(S)的分子和分母多项式的系数，如下所示：

B=[b1,b2,…,bm,bm+1],A=[a1,a2,…,an,an+1]
(4) 合为一步的数字滤波器设计函数butter
调用格式1：[B,A] = BUTTER(N,Wn)
² 功能：设计N阶低通数字滤波器，系数矢量B，A按Z的降幂排列，当Wn为二元矢量时，Wn = [W1 W2]，返回一个2N阶带通滤波器，通带范围：W1 < W < W2，如式3.5.6
调用格式2：[B,A] = BUTTER(N, Wn, 'high') 设计高通数字滤波器.
调用格式3：[B,A] = BUTTER(N, Wn, 'stop') 设计带阻滤波器 Wn = [W1,W2].
调用格式4：[Z,P,K] = BUTTER(N,Wn) 返回低通数字滤波器的N阶零、极点矢量.
★ FREQS
² 调用格式： [H,W] = FREQS(B,A,M)
² 功能：此函数可以求出系统频率响应的数值解，并列出此系统的幅频及相频响应曲线。

参考程序：

熟悉MATLAB

1.butter1.m
％巴特沃兹滤波器的幅频响应图
subplot(1,2,1)； ％分两个窗口，幅频图在第一个窗口
wp=1000;ws=2000;rp=3;rs=40; ％设置指标
[N,wn]=buttord(wp,ws,rp,rs,'s') ％计算巴特沃斯低通滤波器的阶数和3dB截止频率
[B,A]=butter(N,wn,'s'); ％代入N和Wn设计巴特沃斯模拟低通滤波器
[Z,P,K]=buttap(N); ％计算滤波器的零、极点
[h,w]=freqs(B,A,1024); ％计算1024点模拟滤波器频率响应h，和对应的频率点w
％画频率响应幅度图
plot(w,20*log10(abs(h)/abs(h(1))))
grid;
xlabel('频率Hz');ylabel('幅度（dB）');％给x轴和y轴加标注
title('巴特沃斯幅频响应') ％给图形加标题
axis([0,3000, -40,3]);
line([0,2000],[-3,-3]);
line([1000,1000],[-40,3]);
%绘制巴特沃斯滤波器的极点图
subplot(1,2,2) ％在第二个窗口画极点图
p=P';q=Z';
x=max(abs([p,q]));
x=x+0.1;y=x;
axis([-x,x,-y,y]);
axis('square')
plot([-x,x],[0,0]);hold on
plot([0,0],[-y,y]);hold on
plot(real(p),imag(p),'x')
程序运行结果：
N = 7
wn = 1.0359e+003
Z = []
P = -0.2225 + 0.9749i -0.2225 - 0.9749i -0.6235 + 0.7818i -0.6235 - 0.7818i
-0.9010 + 0.4339i -0.9010 - 0.4339i -1.0000
K = 1.0000

图例1
2.butter2.m
设模拟低通滤波器的系统函数为：
b=1;a=[1,1000];
w=[0:1000]*2*pi;％模拟频率为2f，其中f取0~1000Hz
[h,w]=freqs(b,a,w);％计算模拟滤波器的频率响应
subplot(2,2,1)
plot(w/2/pi,abs(h)/abs(h(1)));grid;％画模拟滤波器幅频特性
title('模拟频率响应');
xlabel('f(Hz)');ylabel('幅度');
Fs=1000;
[bz,az]=impinvar(b,a,Fs); ％冲激响应不变法设计数字滤波器
[bzl,azl]=bilinear(b,a,Fs); ％双线性变换法设计数字滤波器
wz=[0:pi/512:pi];
hz1=freqz(bz,az,wz);
hz2=freqz(bzl,azl,wz);
subplot(2,2,2);plot(wz/pi,abs(hz1)/hz1(1));grid; ％画出冲激响应不变法滤波器的幅频图，
axis([0,1,0,1]) ％数字频率wz归一化为0~1
title('冲激响应不变法数字频率响应')
subplot(2,2,3);plot(wz/pi,abs(hz2)/hz2(1));grid; ％画出双线性变换法滤波器的幅频图，
axis([0,1,0,1]) ％数字频率wz归一化为0~1
title('双线性变换法数字频率响应');

图例2

实验五 窗函数设计FIR滤波器
实验目的

1.掌握利用窗函数法设计线性相位FIR数字滤波器的方法
实验原理

用窗函数法设计FIR滤波器时，先根据给定的Wp和N通过傅立叶反变换求得相应的理想滤波器单位脉冲响应，由于是一个无限长序列，所以第二步要选择一个合适的窗函数w(n)来截取到合适的长度（阶数）以保证实现要求的阻带衰减，最后得到FIR滤波器单位冲激响应h(n)＝w(n)。
即：
其中

实验内容

1. 利用窗生成函数boxcar，hanning，bartlett设计FIR低通数字滤波器，指标如下：采样频率Fs＝1000Hz，截止频率ws＝200Hz，滤波器阶数N＝65，画出滤波器的幅度频率特性和窗的形状，参考程序chfir1.m
2. 用基于窗函数的滤波器设计函数fir1，重新设计上述滤波器
实验报告

1.列出本次实验编写的所有文件及各项实验结果的曲线，加注必要的说明
2.写出理论计算低通滤波器的h(n)的方法
3.总结实验体会及实验中存在的问题
思考：

1.比较矩形窗、汉宁窗、巴特里特窗的滤波器特性

FIR设计函数：

熟悉MATLAB

（1）矩形窗生成函数hanning
调用格式：W =HANNING(N)
² 功能：生成N点的汉宁窗
² 说明：N为矩形窗时域的长度，W为行矢量，其元素为对应矩形窗时域每一采样点的幅度值。

（2）基于窗函数的滤波器设计函数fir1
调用格式1：B = FIR1(N,Wn,'ftype'，window)
² 功能：以给定的窗函数（window例如：bartlett(N+1)），和滤波器类型（high、low、stop、band）设计出预期理想频率特性的滤波器。
² 说明：N为滤波器阶数，Wn为通带截止频率如式3.5.3的Wp，B为滤波器系数，即单位抽样响应h(n)＝b(n+1) n＝0,1,2,…N；如果'ftype'和window缺省，该函数得到截止频率为Wn且满足线性相位条件的N阶FIR低通滤波器，window默认为hamming窗。
★