实验-三--用凯塞窗设计线性相位带通FIR滤波器---副本.doc

中国地质大学(北京)
实 验 报 告

课程名称： 数字信号处理
实验名称：用凯塞窗设计线性相位带通FIR滤波器
姓 名： 张淑坤
学 号: 1004133105
班 级： 10041331
指导教师： 陈玉东
评 分：
实验时间： 2015-12-31
用凯塞窗设计线性相位带通FIR滤波器
实验目的：
基于MATLAB环境，熟悉利用窗函数法设计线性相位FIR滤波器的原理和方法。
实验环境：
硬件环境:计算机,软件环境：MATLAB平台。
实验环境原理：
凯塞窗是对于给定的阻带衰减，使主瓣具有最大能量意义下的最佳窗函数，因此具有最陡的过渡带。
窗函数的形式为

其中b为形状参数，I0[x]是第一类修正零阶贝塞尔函数，其幂函数展开式为

给定滤波器的过渡带宽度Dw (rad)和阻带衰减As(dB)，则滤波器的长度和形状参数b可由下列经验公式给出：

设计流程图
用stem(wk) 显示wk(n)
显示N及b
输入带通指标：wp1、wp2、ws1、ws2、As、
编写函数dbpfilter_FIR.m用以计算理想带通滤波器单位抽样响应hd(n)
计算N(N为奇数)及b
编写函数kaiser_WF.m计算凯塞窗函数wk(n)
计算h (n)=hd(n)·wk(n)
计算滤波器的幅频及对数幅频特性
并打印结果，验证指标要求

实验内容要求：
1. 编写计算理想带通滤波器单位抽样响应hd(n)的M函数文件dbpfilter_FIR.m，各变量定义如下：

dbpfilter_FIR(d_omega, N, beta)
% d_omega ---输入数字频率数组（向量）；
% d_omega(1)---阻带下边缘截止频率
% d_omega(2)---通带下边缘截止频率
% d_omega(3)---通带上边缘截止频率
% d_omega(4)---阻带上边缘截止频率
% N ---数字带通滤波器的长度
% hd ---理想数字带通滤波器单位冲激响应
% h ---实际数字带通滤波器单位冲激响应（所设计的）
%wk ---凯泽窗序列
% beta ---凯泽窗参数

2. 编写计算凯泽窗函数wk(n)的M函数文件kaiser_WF.m，各变量定义如下：
function wk= kaiser_WF(N, beta)
% N ---凯泽窗序列的长度
% beta---凯泽窗参数
% wk ---凯泽窗序列

计算贝塞尔函数的参考程序bessel_IM.m如下：

function s=bessel_IM (x)
eps= 10^(-12);
n=1; s=1; D= 1;
while D>(eps*s)
T= x/(2*n);
D= D *T^2;
s= s+D;
n= n+l;
end

3. 编写 .m程序文件，通过调用dbpfr.m和kaiser_WF.m文件，设计下列带通FIR滤波器：
通带允许起伏 ≤l dB，wp1=0.3p，wp2=0.5p
阻带衰减 ≤40 dB，ws1=0.15p，ws2=0.65
实验内容及结果
1、MATLAB程序
function [hd]=ideal_LP(wc,N)
%Ideal Lowpass filter computation
%------------------------------------------
%[hd]=ideal_LP(wc,M)
%hd ---理想低通滤波器单位冲激响应（0<=n<=(n-1)）
%wc ---截至频率（单位弧度/秒）
%N ---理想低通滤波器的长度
%
%

alpha=(N-1)/2;
n=[0:(N-1)];
m=n-alpha+eps;
hd=sin(wc*m)./(pi*m);

function wk= kaiser_WF(N, beta)
for n=1:1:N
wk(n)=bessel_IM (beta*sqrt(1-(1-2*(n-1)/(N-1))^2))/bessel_IM(beta);
end

function s=bessel_IM (x)
N=49;
n=[0:1:N-1];
eps= 10^(-12);
n=1; s=1; D= 1;
while D>(eps*s)
T= x/(2*n);
D= D *T^2;
s= s+D;
n= n+1;
End

function dbpfilter_FIR(d_omega, N, beta)

As=60;
d_omega=[0.15,0.3,0.5,0.65]*pi
delta_w=d_omega(2)-d_omega(1);
N=(As-7.95)/(2.286*delta_w)
N=ceil(N)
M=N;
n=[0:1:N-1];
beta=0.1102*(As-8.7);
wk=kaiser_WF(N,beta);
wc_lower=(d_omega(1)+d_omega(2))/2

wc_upper=(d_omega(3)+d_omega(4))/2
hd=ideal_LP(wc_upper,N)-ideal_LP(wc_lower,N);
h=hd.*wk;

[H,w]=freqz(h,1,1000,'whole');
H=(H(1:501))';
w=(w(1:501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));

figure(1);clf;
subplot(221);
stem(n,hd);
title('理想数字带通滤波器单位冲激响应');
axis([-1,N,-0.4,0.6]);
ylabel('h_d(n)');

subplot(222);
stem(n,wk);
title('凯泽窗');
axis([-1,N,0,1.3]);
ylabel('w_k(n)');

subplot(223);
stem(n,h);
axis([-1,N,-0.4,0.6]);
ylabel('h(n)')
title('设计出的滤波器单位冲激响应');

subplot(224);
plot(w/pi,db);
axis([0,1,-100,10]);
ylabel('|H(w)|')
title('凯泽窗的累加幅度函数');

xlabel('w');
axis([0 1 -80 10]);
set(gca,'XTickMode','manual','XTick',d_omega/pi);
set(gca,'YTickMode','manual','YTick',[-60 0]);
grid;
2、 波形图

理想数字带通滤波器单位冲激响应

凯泽窗函数

实际带通滤波器单位冲激响应

设计出的带通滤波器

实验总结与体会
在老师和同学的热心帮助下，通过本次实验使我熟练掌握了运用MATLAB进行编程，掌握了凯泽窗的设计和FIR滤波器的基本原理，观察实现的曲线图形并比较与理想的误差。培养了我的实际应用的能力，提高分析问题、解决问题的能力。而且提高了我对这门课的了解深度，把理论知识和实践仿真相结合，增强了个人的动手和独立思考能力，也为以后后续课程的学习以及从事实际工作打下良好的基础。

1． 窗函数有哪些指标要求?
窗函数有截短和平滑的作用，窗函数选择的好，可以在相同阶次的情况下，提高滤波器的性能，或是在满足设计要求的情况下，减少滤波器阶数。
选窗标准：
1. 较低的旁瓣幅度，尤其是第一旁瓣；
2. 旁瓣幅度要下降得快，以利于增加阻带衰减；
3. 主瓣宽度要窄，这样滤波器过渡带较窄。
但这三点难以同时满足，当选用主瓣宽度较窄时，虽然得到的幅频特性较陡峭，但通带、阻带波动会明显增加；当选用较低的旁瓣幅度时，虽然得到的幅频特性较平缓匀滑，但过渡带变宽。因此，实际的选择往往是取折衷。
一般选这几个窗之一：矩形窗、三角窗、汉宁窗、海明窗、布拉克曼窗、凯塞窗，可以查查资料比较他们的旁瓣幅度，过渡带宽度和阻带最小衰减后再进行选择。

2. 用窗函数法设计FIR滤波器时，滤波器的过渡带宽度和阻带衰减各与哪些因素有关?
过渡带宽度与窗函数的形状和窗的宽度有关；阻带衰减只有窗函数的 形状决定，不受N的影响。

3. 凯塞窗函数的b参数一般选取范围是多少? b 的大小对窗函数形状以及频谱有何影响?
b是一个可自由选择的参数，它可以同时调整主瓣宽度与旁瓣幅值，b越大，则wk(n)窗宽变得越窄，频谱的旁瓣就越小，但主瓣宽度也相应增加。因而，改变b值就可以对主瓣宽度与旁瓣衰减进行选择，b=0时相当于矩形窗，b=5.44时相当于汉明窗，b=8.5时相当于布莱克曼窗。一般选择4<b<9，这相当于旁瓣幅度与主瓣幅度的比值由3.1%变到0.047%。
报告要求：
1. 列出本实验编写的所有文件及各项结果(包括数据、曲线)，并加注必要的说明。
2. 写出计算理想带通滤波器单位抽样响应hd(n)的方法。
3. 对给定指标要求的带通滤波器，理论计算用凯塞窗设计所需的滤波器长度N和形状参数b
4. 分析实验结果，写出实验体会及实验中存在的问题。

注：实验报告的写作格式见《中国地质大学（北京）小学期实习报告撰写规范》。

评分标准：满分100分

满分
评分
报告整体结构
有实验内容简介，格式规范，图表整洁。
10

能正确设计线性相位带通FIR滤波器
结果正确，附有程序。
60

问题回答
思考问题回答正确。
30

总分