FIR数字滤波器的设计等波纹最佳逼近法.docx
《FIR数字滤波器的设计等波纹最佳逼近法.docx》由会员分享,可在线阅读,更多相关《FIR数字滤波器的设计等波纹最佳逼近法.docx(10页珍藏版)》请在冰豆网上搜索。
FIR数字滤波器的设计等波纹最佳逼近法
FIR数字滤波器的设计--等波纹最佳逼近法
一、等波最佳逼近的原理简介
等波纹最佳逼近法是一种优化设计法,即最大误差最小化准则,它克服了窗函数设计法和频率采样法的缺点,使最大误差(即波纹的峰值)最小化,并在整个逼近频段上均匀分布。
用等波纹最佳逼近法设计的FIR数字滤波器的幅频响应在通带和阻带都是等波纹的,而且可以分别控制通带和阻带波纹幅度,这就是等波纹的含义。
最佳逼近是指在滤波器长度给定的条件下,使加权误差波纹幅度最小化。
与窗函数设计法和频率采样法比较,由于这种设计法使滤波器的最大逼近误差均匀分布,所以设计的滤波器性能价格比最高。
阶数相同时,这种设计法使滤波器的最大逼近误差最小,即通带最大衰减最小,阻带最小衰减最大;指标相同时,这种设计法使滤波器阶数最低。
等波纹最佳逼近法的设计思想 。
用
表示希望逼近的幅度特性函数,要求设计线性相位FIR数字滤波器时,
必须满足线性相位约束条件。
用
表示实际设计的滤波器的幅度特性函数。
定义加权误差函数
为
式中,
为幅度误差加权函数,用来控制不同频带(一般指通带和阻带)的幅度逼近精度。
等波纹最佳逼近法的设计在于找到滤波器的系数向量
,使得在通带和阻带内的最大绝对值幅度误差
为最小,这也就是最大误差最小化问题。
二、等波纹逼近法设计滤波器的步骤和函数介绍
1.根据滤波器的设计指标的要求:
边界频率,通带最大衰减,阻带最大衰等估计滤波器阶数n,确定幅度误差加权函数
2.采用Parks-McClellan算法,获得所设计滤波器的单位脉冲响应
实现FIR数字滤波器的等波纹最佳逼近法的MATLAB信号处理工具函数为firpm和firpmord。
firpm函数采用数值分析中的多重交换迭代算法求解等波纹最佳逼近问题,求的满足等波纹最佳逼近准则的FIR数字滤波器的单位脉冲响应
。
firpmord根据逼近指标,计算采用Parks-McClellan算法等波纹最佳逼近滤波器的最低阶数,误差加权向量w,归一化边界频率向量f。
3对firpm和firpmord的说明
firpm
函数功能:
采用Parks-McClellan算法设计FIR滤波器
函数格式:
hn=firpm(n,f,m,w)
n是滤波器的阶
hn是数字滤波器的单位脉冲响应,其长度为n+1
f是希望滤波器的边界频率向量,要求f是单调增向量,并且从0开始,以1结束,1对于数字频率
m是与f对应的希望滤波器的幅度向量,m和f的长度相等,
表示希望滤波器在频率点
上的幅频响应,m和f给出了希望滤波器的幅度特性。
w是误差加权向量,其长度为f的一半。
表示对m中第i个频率段幅度逼近精度的加权值。
w缺省时,函数默认w全为1,即每个频率段的逼近误差加权值相同。
firpmord函数
函数功能:
根据逼近指标,计算采用Parks-McClellan算法等波纹最佳逼近滤波器的最低阶数n,误差加权向量w和归一化边界频率f。
其返回参数作为firpm函数的调用参数。
函数格式:
[n,f,m,w]=firpmord(f,m,rip,fs)
f可以是归一化边界数字频率向量,也可以是模拟边界频率向量,但必须以0开始,以1结束或fs/2结束,并且其中省略了0和fs/2两个频率点。
fs是时域采样频率,单位Hz。
fs缺省时,函数默认fs=2Hz.但这是f的长度(包括省略0和fs/2两个频率点)是m的两倍,即m中的每个元素表示f给定的一个逼近频段上的希望逼近的幅度值。
rip表示f和m描述的各逼近段允许的波纹振幅(幅频响应最大偏差)。
4.注意事项:
省略fs时,f必须是归一化的数字频率
有时计算的阶数n略小,使设计结果达不到指标要求,这时要取n+1或n+2
三、程序
低通滤波器设计
%Lowpassfilter:
wp=,ws=,N=26
%peakpassbandrippleis,peakstopbandrippleis
clc;clear;clearall;
wp_l=*pi;ws_l=*pi;N1=26;deltal_1=;deltal_2=;%初始化参数
f_l=[0wp_lws_lpi]/pi;%设定归一化的firpm函数参数
deltal=[deltal_1,deltal_2];
a_l=[1100];%幅值
[n,f,m,w_l]=firpmord(f_l(2:
3),[10],deltal);%计算权值
fil_l=firpm(N1,f_l,a_l,w_l);%生成滤波器
figure
(1);
stem(fil_l,'.');
title('impulseresponseoflowpassfilter');xlabel('timesequence');ylabel('amplitude');
fft_l=fft(fil_l,1024);
db_l=20*log10(abs(fft_l));%对数表示
figure
(2);
plot([0:
length(fft_l)-1]/length(fft_l),db_l);%将横坐归化到0-1,单位为1/1024
title('Log-magnituderesponseoflowpassfilter');xlabel('Normalizedfrequency');ylabel('Magnituderesponse(dB)');
axis([0-10010]);%设坐标轴的范围,只取0-2pi的一半
afl=abs(fft_l);%取傅变的幅值
afl_p=zeros(1,1024/2);afl_s=zeros(1,1024/2);%初始化两个一行512列的零序列
afl_p(1:
fix((wp_l/(2*pi))*1024))=afl(1:
fix((wp_l/(2*pi))*1024))-1;%将通带幅值赋给序列一
afl_s(fix((ws_l/(2*pi))*1024)+3:
end)=afl(fix((ws_l/(2*pi))*1024)+3:
end/2);%将阻带幅值赋给序列二
figure(3);
plot([0:
1024/2-1]/1024,afl_p+afl_s);%将通带和阻带的波形显示到一幅图像上
title('passbandandstopbandapproximationerroroflowpassfilter');xlabel('Normalizedfrequency');ylabel('Passbandandstopbandripple');
带通滤波器设计
%Bandpassfilter:
ws1=,wp1=,wp2=,ws2=,N=80
%peakpassbandrippleis,eachpeakstopbandrippleis
ws_b1=*pi;wp_b1=*pi;wp_b2=*pi;ws_b2=*pi;N2=80;deltab_1=;deltab_2=;deltab_3=;
deltab=[deltab_1,deltab_2,deltab_3];
f_b=[0ws_b1wp_b1wp_b2ws_b2pi]/pi;
a_b=[001100];
[n,f,m,w_b]=firpmord(f_b(2:
5),[010],deltab);
fil_b=firpm(N2,f_b,a_b,w_b);
figure(4);
stem(fil_b,'.');
title('impulseresponseofbandpassfilter');xlabel('timesequence');ylabel('amplitude');
fft_b=fft(fil_b,1024);
db_b=20*log10(abs(fft_b));
figure(5);
plot([0:
length(fft_b)-1]/length(fft_b),db_b);
title('Log-magnituderesponseofbandpassfilter');xlabel('Normalizedfrequency');ylabel('Magnituderesponse(dB)');
axis([0-10010]);
afb=abs(fft_b);
afb_p=zeros(1,1024/2);afb_s=zeros(1,1024/2);
afb_p(fix((wp_b1/(2*pi))*1024)+3:
fix((wp_b2/(2*pi))*1024))=afb(fix((wp_b1/(2*pi))*1024)+3:
fix((wp_b2/(2*pi))*1024))-1;
afb_s(fix((ws_b2/(2*pi))*1024)+2:
end)=afb(fix((ws_b2/(2*pi))*1024)+2:
end/2);
afb_s(1:
fix((ws_b1/(2*pi))*1024))=afb(1:
fix((ws_b1/(2*pi))*1024));
figure(6);plot([0:
1024/2-1]/(1024),afb_p+afb_s);
title('passbandandstopbandapproximationerrorofbandpassfilter');xlabel('Normalizedfrequency');ylabel('Passbandandstopbandripple');
高通滤波器设计
%Highpassfilter:
ws=,wp=,N=34
%peakpassbandrippleis,peakstopbandrippleis
ws_h=*pi;wp_h=*pi;N3=34;deltah_1=;deltah_2=;
deltah=[deltah_1,deltah_2];
f_h=[0ws_hwp_hpi]/pi;
a_h=[0011];
[n,f,m,w_h]=firpmord(f_h(2:
3),[01],deltah)
fil_h=firpm(N3,f_h,a_h,w_h);
figure(7);
stem(fil_h,'.');
title('impulseresponseofhighpassfilter');xlabel('timesequence');ylabel('amplitude');
fft_h=fft(fil_h,1024);
db_h=20*log10(abs(fft_h));
figure(8);
plot([0:
length(fft_h)-1]/length(fft_h),db_h);
title('Log-magnituderesponseofhighpassfilter');xlabel('Normalizedfrequency');ylabel('Magnituderesponse(dB)');
axis([0-10010]);
afh=abs(fft_h);
afh_p=zeros(1,1024/2);afh_s=zeros(1,1024/2);
afh_p(fix((wp_h/(2*pi))*1024)+5:
end)=afh(fix((wp_h/(2*pi))*1024)+5:
end/2)-1;
afh_s(1:
fix((ws_h/(2*pi))*1024))=afh(1:
fix((ws_h/(2*pi))*1024));
figure(9);
plot([0:
1024/2-1]/(1024),afh_p+afh_s);
title('passbandandstopbandapproximationerrorofhighpassfilter');xlabel('Normalizedfrequency');ylabel('Passbandandstopbandripple');
仿真图
1)以下三幅图依次为低通滤波器的冲激响应,对数幅频响应和通带阻带波纹
2)以下三幅图依次为带通滤波器的冲激响应,对数幅频响应和通带阻带波纹
3)以下三幅图依次为高通滤波器的冲激响应,对数幅频响应和通带阻带波纹