过渡带,加两个过渡带的低通滤波器,看其对滤波器性能的影响。设计与之相对应的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,追求最佳性能,并检验设计效果。
6
三、实验结果
(一)程序:
% 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)]; %偶对称
7
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)');
8
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);
9
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
10