基于频率抽样设计法线性相位型FIR数字低通滤波器设计报告.docx

目录
任务书 2
摘要： 3
关键字： 3
FIR数字滤波器设计的原理 3
设计过程： 5
实验结果 5
1.程序： 5
2.截图 9
结束语

任务书
题目8
基于频率抽样设计法线性相位型FIR数字低通滤波器设计
主要
内容
1、根据设计指标，利用频率抽样法设计线性相位型FIR数字低通滤波器；
2、对比分析改变抽样点数及过渡点数对数字滤波器性能的影响。
设计
要求
1、根据设计指标，构建线性相位型FIR滤波器之频响抽样值序列；
2、确定数字滤波器的冲击响应及系统函数；
3、对比分析抽样点数及过渡点数的增减对系统频响的影响；
4、生成一个时间序列，利用已设计的FIR滤波器进行滤波，分析滤波效果。
主要
仪器
设备
1、计算机1台，安装MATLAB软件
主要
参考
文献
[美]维纳.K.恩格尔，约翰.G.普罗科斯著，刘树棠译.数字信号处理——使用MATLAB[M].西安：西安交通大学出版社，2002.
飞思科技产品研发中心编著.MATLAB7辅助信号处理技术与应用[M].北京：电子工业出版社，2005.
课程设计进度安排（起止时间、工作内容）
课程设计共设16个设计题目，每班3至4人为1组，1人1套设备，每组选作不同的题目，4个班共分4批。完整课程设计共20学时，为期1周，具体进度如下：
5学时 学习题目相关知识，掌握实现原理；
5学时 用MATLAB语言实现题目要求；
5学时 进一步完善功能，现场检查、答辩；
5学时 完成并提交课程设计报告。
课程设计开始日期
2012.12.31
课程设计完成日期
2013.1.6
课程设计实验室名称
信号处理实验室
地 点
实验楼501
资料下载地址
各班公共邮箱

摘要：
论文从基于频率抽样设计法线性相位型FIR数字低通滤波器的设计原理和设计过程方面来阐述我们的设计思路，并用matlab软件运行后得到了论文中的实验结果，我们还对课程设计中遇到的问题及解决方法进行了阐述，最后描述了自己的心得体会。

关键字：
频率抽样法 线性相位型FIR数字低通滤波器 matlab 有限长冲激响应
一、FIR数字滤波器设计的原理
频率抽样设计法
FIR低通滤波器的设计一般方法有两种，即频率抽样法和窗函数法，频率抽样法设计不同于窗函数法，窗函数是从时域出发，把理想的用一定形状得窗函数截取成有限长的，以此来近似理想的，这样得到的频率响应逼近于所要求的理想的频率响应。
频率抽样法则是从频域出发，把给定的理想频率响应加以等间隔抽样，即然后以此作为实际FIR数字滤波器的频率特性的抽样值，即令，知道后，由DFT定义，可以用频域的这N个抽样值来唯一确定有限长序列，而由的内插公式知道，利用这N个频域抽样值同样可求得FIR滤波器的系统函数及频率响应。这个或将逼近或，和的内插公式为

（2－2）

其中是内插函数
（2－3）
将式（2－3）代入（2－2）式，化简后可得

即

从内插公式（2－2）看到，在各频率抽样点上，滤波器的实际频率响应是严格地和理想频率响应数值相等，即。但是在抽样点之间的频率响应则是由各抽样点的加权内插函数的延伸叠加而形成的，因而有一定的逼近误差，误差大小取决于理想频率响应曲线形状，理想频率响应特性变化越平缓，则内插值越接近理想值，逼近误差越小，如下图梯形理想频率特性所示。

反之，如果抽样点之间的理想频率特性变化越陡，则内插值与理想值之误差就越大，因而在理想频率特性的不连续点附近，就会产生肩峰和波纹，

当为实数时，满足，由此得出，，也就是说，的模以为对称中心呈偶对称，的相角以为对称中心呈奇对称。再利用线性相位的条件，即可得到（N为偶数）：

二、设计过程：
（一）设计思想：
根据给定DLPF幅频特性要求（通带截止频率ωp=0.5π，通带最大衰减αp=0.5 dB，阻带截止频率ωs=0.6π，阻带最小衰减αs=50 dB）取得DLPF的H(k)。设冲激响应h(n)为偶对称，N为偶数。根据线性相位约束条件设计FIR线性相位数字滤波器的H(k)。根据H(k)生成DLPF的h(n)。利用频率抽样法设计加一个过渡带，加两个过渡带的低通滤波器，看其对滤波器性能的影响。设计与之相对应的DLPF，给出窗函数及所设计滤波器的幅度特性，对比分析DLPF幅频特性是否符合要求最后利用x=rand(1,sizex)函数随随机生成一个序列来验证设计的滤波器是否具有低通滤波器的特性。
（二）功能结构：
FIR数字滤波器具有严格的线性相位，低通滤波器只能让低频的通过而把高频的部分滤掉。
设计步骤：
（三）设计步骤
① 根据给定DLPF幅频特性要求（通带截止频率ωp=0.5π，通带最大衰减αp=0.5 dB，阻带截止
频率ωs=0.6π，阻带最小衰减αs=50 dB）取得DLPF的X(K)；
② 根据线性相位型数字滤波器条件，构建线性相位型DLPF的X(K)；
③ 根据X(K)生成DLPF的h(n)
④ 设计与之相对应的DLPF，给出窗函数及所设计滤波器的幅度特性，对比分析DLPF幅频特
性是否符合要求；
⑤ 试说明过渡点对所设计数字滤波器性能的影响；
⑥ 产生一个有干扰频率的时域序列（借助FFT分析说明其有干扰），使之通过所设计的DLPF，
对滤波输出结果作出分析，说明输出结果。
⑦ 扩展部分：自拟指标，设计一个DBPF，追求最佳性能，并检验设计效果。

三、实验结果
（一）程序：
% wp=0.5*pi;Rp=0.5dB; ws=0.6*pi;As=50dB;

频率抽样法

wp=0.5*pi;ws=0.6*pi;
tr_width=ws-wp
% 用频率抽样法设计FIR滤波器,过渡带内一个样本T1，N=40。
alpha=(N-1)/2; %N为偶数
l=[0:1:N-1]; wl=(2*pi/N)*l;
Hrs=[ones(1,11),T1,zeros(1,17),T1,ones(1,10)]; %偶对称
Hdr=[1 1 0 0 ];wdl=[0 0.5 0.6 1]; k1=0:(N/2-1); k2=(N/2+1):N-1;%依据公式7-107
angH=[-alpha*(2*pi)/N*k1,0,alpha*(2*pi)/N*(N-k2)];
H=Hrs.*exp(1i*angH);
h=real(ifft(H,N));[H,w]=freqz(h,1,1000,'whole');
%[db mag pha grd w]=freqz_m(h,1);
db=20*log10((abs(H)+eps)/max(abs(H)));%求FIR滤波器频响的dB值
delta_w=2*pi/1000; %将2pi等分1000份
%plot
figure(1);
subplot(221);plot(wl(1:21)/pi,Hrs(1:21),'o',wdl,Hdr,'linewidth',2);
title('理想滤波器频域波形');
axis[0 1 -0.1 ,1.2]); ylabel('Hr(k)');
set(gca,'XTickMode','manual','XTick',[0 0.5 0.6 1]);
set(gca,'YTickMode','manual','YTick',[0 T1 1]); grid;
subplot(222); stem(l,h,'m');
title('单位脉冲响应');axis([-1,N,-0.15,0.5]);ylabel('h(n)');
subplot(223); plot(w/pi,abs(H),wl(1:31)/pi,Hrs(1:31),'o','linewidth',2);axis([0 1 -0.2,1.2]);
title('频域抽样');xlabel('频率');ylabel('Hr(w)');grid;
set(gca,'XTickMode','manual','XTick',[0 0.5 0.6 1]);
set(gca,'YTickMode','manual','YTick',[0 T1 1]);
subplot(224); plot(w/pi,db,'r');axis([0,1,-100,10]);grid; ylabel('分贝');
title('幅频响应');xlabel('频率');
set(gca,'XTickMode','manual','XTick',[0,0.5,0.6,1]);
set(gca,'YTickMode','manual','YTick',[-50,0]);
% 用频率抽样法设计FIR滤波器,过渡带内无样本，N=20。
tr_width=ws-wp;
alpha=(N-1)/2; %N为偶数
l=[0:1:N-1]; wl=(2*pi/N)*l;
Hrs=[ones(1,6),zeros(1,9),ones(1,5)]; %偶对称
Hdr=[1 1 0 0 ];wdl=[0 0.5 0.6 1]; k1=0:(N/2-1); k2=(N/2+1):N-1;

angH=[-alpha*(2*pi)/N*k1,0,alpha*(2*pi)/N*(N-k2)]; H=Hrs.*exp(j*angH);
h=real(ifft(H,N));[H,w]=freqz(h,1,1000,'whole');
%[db mag pha grd w]=freqz_m(h,1);
db=20*log10((abs(H)+eps)/max(abs(H)));%求FIR滤波器频响的dB值
delta_w=2*pi/1000; %将2pi等分1000份
%plot
figure(2);clf;
subplot(221)plot(wl(1:11)/pi,Hrs(1:11),'o',wdl,Hdr,'linewidth',2);
title('理想滤波器频域波形');
axis([0 1 -0.1 ,1.2]); ylabel('Hr(k)');
set(gca,'XTickMode','manual','XTick',[0 0.5 0.6 1]);
%set(gca,'YTickMode','manual','YTick',[0 T1 1]);
grid;
subplot(222); stem(l,h,'m');
title('单位脉冲响应');axis([-1,N,-0.15,0.55]);ylabel('h(n)');
subplot(223); plot(w/pi,abs(H),wl(1:11)/pi,Hrs(1:11),'o','linewidth',2);axis([0 1 -0.2,1.2]);
title('频域抽样');xlabel('频率');ylabel('Hr(w)');grid;
%set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 1]);
%set(gca,'YTickMode','manual','YTick',[0 T1 1]);
subplot(224); plot(w/pi,db,'r');axis([0,1,-100,10]);grid; ylabel('分贝');
title('幅频响应');xlabel('频率');
set(gca,'XTickMode','manual','XTick',[0,0.5,0.6,1]);
set(gca,'YTickMode','manual','YTick',[-50,0]);
% 用频率抽样法设计FIR滤波器,过渡带内两个样本T1、T2，N=60。
tr_width=ws-wp;
alpha=(N-1)/2; %N为偶数
l=[0:1:N-1]; wl=(2*pi/N)*l;
Hrs=[ones(1,16),T1,T2,zeros(1,25),T2, T1,ones(1,15)]; %偶对称
Hdr=[1 1 0 0 ];wdl=[0 0.5 0.6 1]; k1=0:(N/2-1); k2=(N/2+1):N-1;
angH=[-alpha*(2*pi)/N*k1,0,alpha*(2*pi)/N*(N-k2)]; H=Hrs.*exp(j*angH);
h=real(ifft(H,N));[H,w]=freqz(h,1,1000,'whole');
%[db mag pha grd w]=freqz_m(h,1);
db=20*log10((abs(H)+eps)/max(abs(H)));%求FIR滤波器频响的dB值
delta_w=2*pi/1000; %将2pi等分1000份
%plot
figure(3);clf;
subplot(221);plot(wl(1:31)/pi,Hrs(1:31),'o',wdl,Hdr,'linewidth',2);
title('理想滤波器频域波形');
axis([0 1 -0.1 ,1.2]); ylabel('Hr(k)');
set(gca,'XTickMode','manual','XTick',[0 0.5 0.6 1]);
set(gca,'YTickMode','manual','YTick',[0 T2 T1 1]); grid;
subplot(222); stem(l,h,'m');
title('单位脉冲响应');axis([-1,N,-0.15,0.55]);ylabel('h(n)');
subplot(223); plot(w/pi,abs(H),wl(1:31)/pi,Hrs(1:31),'o','linewidth',2);axis([0 1 -0.2,1.2]);
title('频域抽样');xlabel('frequency in pi units');ylabel('Hr(w)');grid;
set(gca,'XTickMode','manual','XTick',[0 0.5 0.6 1]);
set(gca,'YTickMode','manual','YTick',[0 T2 T1 1]);
subplot(224); plot(w/pi,db,'r');axis([0,1,-100,10]);grid; ylabel('Decibels');
title('幅频响应');xlabel('频率');
set(gca,'XTickMode','manual','XTick',[0,0.5,0.6,1]);
set(gca,'YTickMode','manual','YTick'[-50,0]);
figure(4);
subplot(211); plot(w/pi,db,'r');axis([0,1,-100,10]);grid; ylabel('分贝');
title('幅频响应');xlabel('频率');
set(gca,'XTickMode','manual','XTick',[0,0.5,0.6,1]);
set(gca,'YTickMode','manual','YTick',[-50,0]);
subplot(212); plot(w/pi,angle(H),'b');axis([0,1,-3.5,3.5]);grid; ylabel('rad');
title('相频响应');
%生成一个时间序列，利用已设计的FIR滤波器进行滤波，分析滤波效果
alpha=(N-1)/2; %N为偶数
l=[0:1:N-1]; wl=(2*pi/N)*l;
Hrs=[ones(1,16),T1,T2,zeros(1,25),T2, T1,ones(1,15)]; %偶对称
Hdr=[1 1 0 0 ];wdl=[0 0.5 0.6 1]; k1=0:(N/2-1); k2=(N/2+1):N-1;
angH=[-alpha*(2*pi)/N*k1,0,alpha*(2*pi)/N*(N-k2)]; H=Hrs.exp(j*angH);
h=real(ifft(H,N));
x=rand(1,60);
LL=0:29;
X=fft(x,60);
XX=[X(1:30)];
y=conv(x,h); %x(n)和h(n)做卷积
yL=length(y);
disp('y(n)序列长度为N=');disp(yL);%求卷积后的序列长度
Y=fft(y,60);
YY=[Y(1:30)];
%plot
figure(5)
m=0:1:length(x)-1;
subplot(2,2,1);
title('序列时域波形');
stem(m,x,'.');
xlabel('时间序列 n');ylabel('x(n)');
subplot(2,2,2);
title('序列频域');
stem(LL,abs(XX),'.');
xlabel('k=0,1,2,...,29');ylabel('|X(k)|');axis([0,30,0,30]);
m=0:1:length(y)-1;
subplot(2,2,3);
title('卷积后波形');
stem(m,y,'.');
xlabel('时间序列 n');ylabel('y(n)');
subplot(2,2,4);
title('滤波后的频域抽样');
stem(LL,abs(YY),'.');
xlabel('k=0,1,2,...,29');ylabel('|Y(k)|');axis([0,30,0,30]);

窗函数法
% wp=0.5*pi;Rp=0.5dB; ws=0.6*pi;As=50dB;
% 查第237页表7.1可知用Hamming窗、Blackman窗均可(最小阻带衰减>=50dB)，
% 但Hamming窗具有较小的过渡带(6.6pi/M)，故选择Hamming窗。
wp=0.5*pi;ws=0.6*pi;
tr_width=ws-wp;
M=ceil(6.6*pi/tr_width)+1;
n=[0:1:M-1];
wc=(ws+wp)/2;
hd=ideal_lp(wc,M);
w_ham=(hamming(M))';
h=hd.*w_ham;
[H,w]=freqz(h,1,1000,'whole');
H=(H(1:1:501))';w=(w(1:1:501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
grd=grpdelay(h,1,w);
delta_w=2*pi/1000;
Rp=0.5;
As=50;
figure(1);
subplot(2,2,1);stem(n,hd);title('理想脉冲响应');
axis([0 M-1 -0.1 0.56]);xlabel('n');ylabel('hd(n)');
subplot(2,2,2);plot(w/pi,pha);
axis([0,1,-3.5,3.5]); ylabel('角度');xlabel('频率');
title('相位频率响应');
set(findobj(gcf,'Type','line','Color',[0 0 1]), 'Color','b', 'LineWidth',2);
grid;
subplot(2,2,3);stem(n,h);title('实际脉冲响应');
axis([0 M-1 -0.1 0.56]);xlabel('n');ylabel('h(n)');
subplot(2,2,4); plot(w/pi,db);
axis([0,1,-100,10]); ylabel('幅度');title('幅度响应');
set(gca,'XTickMode','manual','XTick',[0 0.2 0.3 1]);
set(gca,'YTickMode','manual','YTick',[-50 0]);grid;
pause

实验中用到的.m文件
hd=ideal_lp(wc,M)
function hd=ideallp(wc,M)
alpha=(M-1)/2;
n=[0:1:(M-1)];
m=n-alpha+eps;
hd=sin(wc*m)./(pi*m);
end

（二）截图

图1.添加一个过渡点所设计的滤波器、冲激响应h(n)、幅度特性（N=40，过渡点T1=0.37）

图2.不添加过渡点设计的低通滤波器、冲激响应h(n)、幅度特性（N=20，无过渡点）

图3.添加两个过渡点设计的低通滤波器、冲激响应h(n)、幅度特性（N=60，过渡点T1=0.59 T2=0.11）

图4.添加两个过渡点设计的低通滤波器幅频和相频响应

图5.添加两个过渡点设计的低通滤波器的冲激响应h(n)

图6.利用随机生成的序列x(n)验证添加两个过渡点的DLPF性能

图7.利用窗函数法设计的DLPF（冲激响应hdn)、h(n)、相频特性和幅频特性）

图8.对比频率采样法设计的DLPF和窗函数法设计的DLPF的幅度特性

（三）实验分析
1、设计与之相对应的DLPF，给出窗函数及所设计滤波器的幅度特性，对比分析DLPF幅频特性是否符合要求；
观察利用窗函数法设计数字低通滤波器和利用频率采样法设计的数字低通滤波器的幅度特性，通过对比分析，N=60，取两个过渡点T1=0.59 T2=0.11的设计参数下最大阻带最小衰减在50db以下，所设计的DLPF符合要求。

2、试说明过渡点对所设计数字滤波器性能的影响；
在两个频率取样点之间，频率响应值是由各个取样点间内插函数的加权确定，因此，存在着逼近误差，并且误差大小取决于理想频率响应的曲线形状和取样点数N的大小。为了减小逼近误差，即减小通带边缘处由于取样点的突然变化而引起的滤波器特性的振荡起伏，可以在理想滤波器频响的不连续点边缘增加一定的过度取样点，增加过渡带宽，来减小频带边缘的起伏振荡，进而增大阻带衰减。这些取样点上的取值不同，产生的效果也将不同。一般过渡带去1-3点取样值即可得到满意的结果。在本次数字低通滤波器的设计中，不加过渡取样点时，阻带最小衰减为-20dB；增加一个过渡取样点，阻带最小衰减将增为-44dB—54dB；增加两个过渡取样点，则阻带最小衰减增为-65dB—75dB。

3、产生一个有干扰频率的时域序列（借助FFT分析说明其有干扰），使之通过所设计的DLPF，对滤波输出结果作出分析，说明输出结果。
利用x=rand()产生了一个随机序列，通过FFT变换后可看到在频域上有干扰，使x(n)所设计的DLPF后，观察输出的y(n)和它的频谱特性，发现频率高于0.5π的信号被滤掉，滤波器性能良好。

四、结束语
通过本次数字信号课程设计，我熟练掌握了MATLAB工具软件在工程设计中的使用；熟练掌握了利用频率采样法设计FIR线性相位型数字低通滤波器的设计过程。通过利用窗函数法和频率采样法两种不同的方法设计数字低通滤波器，我理解了窗函数设计法是在时域内，以有限长冲激响应h(n)去逼近无限长冲激响应，从而实现FIR滤波器的设计。而本次实验采用的频率取样法则是在频域内，以有限个频率响应的取样，去逼近所期望的理想频率响应函数。频率采样法利用到第三章中DFT的应用和内插函数的定义。
本次实验要求设计线性相位的FIR数字滤波器，所以对H(k)有一定的约束条件，本次设计中采用h(n)为偶对称，N为偶数进行设计，即保证Hk=Hn-k, θk=-(N-1)/N*kπ。在频率采样法中，过渡采样点对滤波器的设计有着至关重要的作用，N值越大，过渡点越多，DLPF的频率响应越接近理想的频率响应，但应注意N值越大，计算量越大，过渡取样点的数量的增加也使过渡带宽增加，所以要设逼近理想的滤波器要合理设计参数。通过学习可知，频率抽样法的实现原理很简单，计算时可以借用DFT。通过改变阶数N和设置过渡点能够最终得到满意的结果，但是和窗函数设计法一样，其主要缺点是不易精确的确定滤波器通带和阻带的截止频率。
在本次课程设计中，我也出现了不少问题，虽然对频率采样法的原理有了深刻的理解，但是利用MATLAB实现起来还是要花费很大的功夫，再调用库函数和算法的实现过程中遇到了不少问题，但是通过查阅书籍和网上查找也一一解决。

;

20