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

目录

任务书 (2)

摘要： (3)

关键字： (3)

FIR数字滤波器设计的原理 (3)

设计过程： (5)

实验结果 (5)

1.程序： (5)

2.截图 (9)

结束语

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

关键字：

频率抽样法 线性相位型FIR 数字低通滤波器 matlab 有限长冲激响应

一、FIR 数字滤波器设计的原理 频率抽样设计法

FIR 低通滤波器的设计一般方法有两种，即频率抽样法和窗函数法，频率抽样法设计不同于窗函数法，窗函数是从时域出发，把理想的()d h n 用一定形状得窗函数截取成有限长的()h n ，以此()h n 来近似理想的()d h n ，这样得到的频率响应()jw H e 逼近于所要求的理想的频率响应()jw d H e 。

频率抽样法则是从频域出发，把给定的理想频率响应()jw d H e 加以等间隔抽样，即2()|

()jw d d w k N

H e H k π

==然后以此()d H k 作为实际FIR 数字滤波器的频率特性

的抽样值()H k ，即令2()()()|

0,1,,1jw d d w k N

H k H k H e k N π====-，知道()

H k 后，由DFT 定义，可以用频域的这N 个抽样值()H k 来唯一确定有限长序列()h n ，而由()X z 的内插公式知道，利用这N 个频域抽样值()H k 同样可求得FIR 滤波器的系统函数()H z 及频率响应()jw H e 。这个()H z 或()jw H e 将逼近()d H z 或

()jw d H e ，()H z 和()jw H e 的内插公式为

1

1

01()

()1N

N k k N

z H k H z N

W z ----=-=

-∑ 10

2()()()N jw

k H e H k w k N

π

-==Φ-

∑ （2－2）

其中()w Φ是内插函数

1()

2

sin(

)12

()sin()2

N jw wN

w e

w N --Φ= （2－3） 将式（2－3）代入（2－2）式，化简后可得

11

()20

sin()12()()sin()

2k N N j j w jw N k wN H e e H k e w k N N

ππ----==-∑ 即 11

(1)()2

sin[()]12()()sin()

2k N N j N j w jw N k w k N N H e e H k e w k N N

πππ----=-=??-∑

21

0()()N j nk N

n H k h n e

π--==∑

从内插公式（2－2）看到，在各频率抽样点上，滤波器的实际频率响应是严格地和理想频率响应数值相等，即22()()()()j k j

k N

N

d d H e

H k H k H e

ππ===。但是在

抽样点之间的频率响应则是由各抽样点的加权内插函数的延伸叠加而形成的，因而有一定的逼近误差，误差大小取决于理想频率响应曲线形状，理想频率响应特性变化越平缓，则内插值越接近理想值，逼近误差越小，如下图梯形理想频率特性所示。

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

当()h n 为实数时，满足**()(())()()N N H k H N k R k H N k =-=-，由此得出

()()H k H N k =-，()()Q k Q N k =--，也就是说，()H k 的模()H k 以2k N =为对称中心呈偶对称，()H k 的相角()Q k 以2k N =为对称中心呈奇对称。再利用线性相位的条件1

()2

jw N Q e w -=-

，即可得到（N 为偶数）

： 210,,122()0221()1,,122N N k k N N

Q k k N N

N k k N N

ππ?-??

??

-=- ?

???????

?

==

???-????

-=+- ? ??????

?

21()2

21

()()2()0,,12()02()1,,1

2N j k N N j N k N N H k e k N

H k k N H N k e k N ππ----???=-? ?

???

?==

??

???-=+-? ????

()12

(0)sin 2sin 2N j j N H H e e N ωωωω--???

?????=?

??? ?

????

()121

sin sin 22sin sin 22N

k k k N N H k N N k k N N 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=ideal_lp(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)