为了从hd(n)得到一个FIR滤波器,必须同时在两边截取hd(n)。
而要得到一个因果的线性相位滤波器,它的h(n)长度为N,必须有:
这种操作叫做加窗,h(n)可以看做是hd(n)与窗函数w(n)的乘积:
h(n)=hd(n)w(n)
其中
根据w(n)的不同定义,可以得到不同的窗结构。
在频域中,因果FIR滤波器响应H(e^jw)由Hd(e^jw)和窗响应W(e^jw)的周期卷积得到,即
常用的窗函数有矩形窗、巴特利特(BARTLETT)窗、汉宁(HANNING)窗、海明(HAMMING)窗、布莱克曼(BLACKMAN)窗、凯泽(KAISER)窗等。
4.4FIR滤波器的窗函数设计法
FIR滤波器的设计方法有许多种,如窗函数设计法、频率采样设计法和最优化设计法等。
窗函数设计法的基本原理是用一定宽度窗函数截取无限脉冲响应序列获得有限长的脉冲响应序列,主要设计步骤为:
(1)通过傅里叶逆变换获得理想滤波器的单位脉冲响应hd(n)。
(2)由性能指标确定窗函数W(n)和窗口长度N。
(3)求得实际滤波器的单位脉冲响应h(n),h(n)即为所设计FIR滤波器系数向量b(n)。
五、设计内容
5.1设计题目:
1-1.试用MATLAB设计一巴特沃斯低通数字滤波器,要求通带截至频率Wp=30HZ,主带截至频率为Ws=35HZ,通带衰减不大于0.5DB,主带衰减不小于40DB,抽样频Fs=100HZ。
1-2.基于Butterworth模拟滤波器原型,使用双线性状换设计数字滤波器:
各参数值为:
通带截止频率Omega=0.2*pi,阻带截止频率Omega=0.3*pi,通带波动值Rp=1dB,阻带波动值Rs=15dB,设Fs=20000Hz。
1-3设计一巴特沃斯高通数字滤波器,要求通带截止频率0.6*pi,通带衰减不大于1dB,阻带衰减15DB,抽样T=1。
1-4.设计一巴特沃斯带阻数字滤波器,要求通带上下截至频率为0.8*PI、0.2*PI,通带衰减不大于1DB,阻带上下截至频率0.7*PI、0.4*PI阻带衰减不小于30DB,
2-1.用窗函数法设计一个线性相位FIR低通滤波器,并满足性能指标:
通带边界频率
Wp=0.5*pi,阻带边界频率Ws=0.66*pi,阻带衰减不小于40dB,通带波纹不大于3dB。
选择汉宁窗。
2-4.用海明窗设计一个FIR滤波器,其中Wp=0.2*pi,Ws=0.3*pi,通带衰减不大于0.25dB,阻带衰减不小于50dB。
5.2设计程序代码及结果:
1-1一.试用MATLAB设计一巴特沃斯低通数字滤波器,要求通带截至频率Wp=30HZ,阻带截至频率为Ws=35HZ,通带衰减不大于0.5DB,阻带衰减不小于40DB,抽样频Fs=100HZ。
代码为:
fp=30;
fs=35;
Fs=100;
wp=2*pi*fp/Fs;
ws=2*pi*fs/Fs;
wp=tan(wp/2);
ws=tan(ws/2);%通带最大衰减为0.5dB,阻带最小衰减为40dB
[N,wn]=buttord(wp,ws,0.5,40,'s');%模拟低通滤波器极零点
[z,p,k]=buttap(N);%由极零点获得转移函数参数
[b,a]=zp2tf(z,p,k);%由原型滤波器获得实际低通滤波器
[B,A]=lp2lp(b,a,wp);
[bz,az]=bilinear(B,A,.5);
[h,w]=freqz(bz,az,256,Fs);
figure
plot(w,abs(h))
gridon
图1巴特沃斯数字低通滤波器
1-2基于Butterworth模拟滤波器原型,使用双线性状换设计数字滤波器:
各参数值为:
通带截止频率Omega=0.2*pi,阻带截止频率Omega=0.3*pi,通带波动值Rp=1dB,阻带波动值Rs=15dB,设Fs=4000Hz。
代码:
wp=0.2*pi;ws=0.3*pi;
Fs=4000;T=1/Fs;
OmegaP=(2/T)*tan(wp/2);
OmegaS=(2/T)*tan(ws/2);
rp=1;rs=15;as=15;
ripple=10^(-rp/20);attn=10^(-rs/20);
[n,wn]=buttord(OmegaP,OmegaS,rp,rs,'s');
[z,p,k]=Buttap(n);
[b,a]=zp2tf(z,p,k);
[bt,at]=lp2lp(b,a,wn);
[b,a]=bilinear(bt,at,Fs);
[db,mag,pha,grd,w]=freqz_m(b,a);
%
%下面绘出各条曲线
subplot(2,2,1);plot(w/pi,mag);title('MagnitudeFrequency幅频特性');
xlabel('w(/pi)');ylabel('|H(jw)|');
axis([0,1,0,1.1]);
set(gca,'XTickMode','manual','XTick',[00.20.31]);
set(gca,'YTickMode','manual','YTick',[0attnripple1]);grid
subplot(2,2,2);plot(w/pi,db);title('MagnitudeFrequency幅频特性(db)');
xlabel('w(/pi)');ylabel('dB');
axis([0,1,-30,5]);
set(gca,'XTickMode','manual','XTick',[00.20.31]);
set(gca,'YTickMode','manual','YTick',[-60-as-rp0]);grid
subplot(2,2,3);plot(w/pi,pha/pi);title('PhaseFrequency相频特性');
xlabel('w(/pi)');ylabel('pha(/pi)');
axis([0,1,-1,1]);
subplot(2,2,4);plot(w/pi,grd);title('GroupDelay群延时');
xlabel('w(/pi)');ylabel('Sample');
axis([0,1,0,15]);
set(gca,'XTickMode','manual','XTick',[00.20.31]);grid
运行结果:
图2巴特沃思数字低通滤波器幅频-相频特性
1-3设计一巴特沃斯高通数字滤波器,要求通带截止频率0.6*pi,通带衰减不大于1dB,阻带衰减15DB,抽样T=1。
Wp=0.6*pi;
Ws=0.4*pi;
Ap=1;
As=15;
[N,wn]=buttord(Wp/pi,Ws/pi,Ap,As);%计算巴特沃斯滤波器阶次和截止频率
%频率变换法设计巴特沃斯高通滤波器
[db,mag,pha,grd,w]=freqz_m(b,a);%数字滤波器响应
plot(w,mag);
title('数字滤波器幅频响应|H(ej\Omega)|')
图3巴特沃斯数字高通滤波器
2-1用窗函数法设计一个线性相位FIR低通滤波器,并满足性能指标:
通带边界频率
Wp=0.5*pi,阻带边界频率Ws=0.66*pi,阻带衰减不小于40dB,通带波纹不大于3dB。
选择汉宁窗。
代码:
wp=0.5*pi;
ws=0.66*pi;
wdelta=ws-wp;
N=ceil(8*pi/wdelta)
ifrem(N,2)==0
N=N+1;
end
);
1数字滤波器及传统设计方法
数字滤波器可以理解为是一个计算程序或算法,将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列,并在转化过程中,使信号按预定的形式变化。
数字滤波器有多种分类,根据数字滤波器冲激响应的时域特征,可将数字滤波器分为两种,即无限长冲激响应(IIR)滤波器和有限长冲激响应(FIR)滤波器。
IIR数字滤波器具有无限宽的冲激响应,与模拟滤波器相匹配。
所以IIR滤波器的设计可以采取在模拟滤波器设计的基础上进一步变换的方法。
FIR数字滤波器的单位脉冲响应是有限长序列。
它的设计问题实质上是确定能满足所要求的转移序列或脉冲响应的常数问题,设计方法主要有窗函数法、频率采样法和等波纹最佳逼近法等。
在对滤波器实际设计时,整个过程的运算量是很大的。
例如利用窗函数法【2】设计M阶FIR低通滤波器时,首先要根据
(1)式计算出理想低通滤波器的单位冲激响应序列,然后根据
(2)式计算出M个滤波器系数。
当滤波器阶数比较高时,计算量比较大,设计过程中改变参数或滤波器类型时都要重新计算。
(1)
(2)
设计完成后对已设计的滤波器的频率响应要进行校核,要得到幅频相频响应特性,运算量也是很大的。
我们平时所要设计的数字滤波器,阶数和类型并不一定是完全给定的,很多时候都是要根据设计要求和滤波效果不断的调整,以达到设计的最优化。
在这种情况下,滤波器的设计就要进行大量复杂的运算,单纯的靠公式计算和编制简单的程序很难在短时间内完成设计。
利用MATLAB强大的计算功能进行计算机辅助设计,可以快速有效的设计数字滤波器,大大的简化了计算量,直观简便。
2数字滤波器的MATLAB设计
2.1FDATool界面设计
2.1.1FDATool的介绍
FDATool(FilterDesign&AnalysisTool)是MATLAB信号处理工具箱里专用的滤波器设计分析工具,MATLAB6.0以上的版本还专门增加了滤波器设计工具箱(FilterDesignToolbox)。
FDATool可以设计几乎所有的基本的常规滤波器,包括FIR和IIR的各种设计方法。
它操作简单,方便灵活。
FDATool界面总共分两大部分,一部分是DesignFilter,在界面的下半部,用来设置滤波器的设计参数,另一部分则是特性区,在界面的上半部分,用来显示滤波器的各种特性。
DesignFilter部分主要分为:
FilterType(滤波器类型)选项,包括Lowpass(低通)、Highpass(高通)、Bandpass(带通)、Bandstop(带阻)和特殊的FIR滤波器。
DesignMethod(设计方法)选项,包括IIR滤波器的Butterworth(巴特沃思)法、ChebyshevTypeI(切比雪夫I型)法、ChebyshevTypeII(切比雪夫II型)法、Elliptic(椭圆滤波器)法和FIR滤波器的Equiripple法、Least-Squares(最小乘方)法、Window(窗函数)法。
FilterOrder(滤波器阶数)选项,定义滤波器的阶数,包括SpecifyOrder(指定阶数)和MinimumOrder(最小阶数)。
在SpecifyOrder中填入所要设计的滤波器的阶数(N阶滤波器,SpecifyOrder=N-1),如果选择MinimumOrder则MATLAB根据所选择的滤波器类型自动使用最小阶数。
FrenquencySpecifications选项,可以详细定义频带的各参数,包括采样频率Fs和频带的截止频率。
它的具体选项由FilterType选项和DesignMethod选项决定,例如Bandpass(带通)滤波器需要定义Fstop1(下阻带截止频率)、Fpass1(通带下限截止频率)、Fpass2(通带上限截止频率)、Fstop2(上阻带截止频率),而Lowpass(低通)滤波器只需要定义Fstop1、Fpass1。
采用窗函数设计滤波器时,由于过渡带是由窗函数的类型和阶数所决定的,所以只需要定义通带截止频率,而不必定义阻带参数。
MagnitudeSpecifications选项,可以定义幅值衰减的情况。
例如设计带通滤波器时,可以定义Wstop1(频率Fstop1处的幅值衰减)、Wpass(通带范围内的幅值衰减)、Wstop2(频率Fstop2处的幅值衰减)。
当采用窗函数设计时,通带截止频率处的幅值衰减固定为6db,所以不必定义。
WindowSpecifications选项,当选取采用窗函数设计时,该选项可定义,它包含了各种窗函数。
2.1.2带通滤波器设计实例
本文将以一个FIR滤波器的设计为例来说明如何使用MATLAB设计数字滤波器:
在小电流接地系统中注入83.3Hz的正弦信号,对其进行跟踪分析,要求设计一带通数字滤波器,滤除工频及整次谐波,以便在非常复杂的信号中分离出该注入信号。
参数要求:
96阶FIR数字滤波器,采样频率1000Hz,采用Hamming窗函数设计。
本例中,首先在FilterType中选择Bandpass(带通滤波器);在DesignMethod选项中选择FIRWindow(FIR滤波器窗函数法),接着在WindowSpecifications选项中选取Hamming;指定FilterOrder项中的SpecifyOrder=95;由于采用窗函数法设计,只要给出通带下限截止频率Fc1和通带上限截止频率Fc2,选取Fc1=70Hz,Fc2=84Hz。
设置完以后点击DesignFilter即可得到所设计的FIR滤波器。
通过菜单选项Analysis可以在特性区看到所设计滤波器的幅频响应、相频响应、零极点配置和滤波器系数等各种特性。
设计完成后将结果保存为1.fda文件。
在设计过程中,可以对比滤波器幅频相频特性和设计要求,随时调整参数和滤波器类型,
以便得到最佳效果。
其它类型的FIR滤波器和IIR滤波器也都可以使用FDATool来设计。
图1滤波器幅频和相频响应(特性区)
Fig.1MagnitudeResponseandPhaseResponseofthefilter
2.2程序设计法
在MATLAB中,对各种滤波器的设计都有相应的计算振幅响应的函数【3】,可以用来做滤波器的程序设计。
上例的带通滤波器可以用程序设计:
c=95;%定义滤波器阶数96阶
w1=2*pi*fc1/fs;
w2=2*pi*fc2/fs;%参数转换,将模拟滤波器的技术指标转换为数字滤波器的技术指标
window=hamming(c+1);%使用hamming窗函数
h=fir1(c,[w1/piw2/pi],window);%使用标准响应的加窗设计函数fir1
freqz(h,1,512);%数字滤波器频率响应
在MATLAB环境下运行该程序即可得到滤波器幅频相频响应曲线和滤波器系数h。
篇幅所限,这里不再将源程序详细列出。
3Simulink仿真
本文通过调用Simulink中的功能模块构成数字滤波器的仿真框图,在仿真过程中,可以双击各功能模块,随时改变参数,获得不同状态下的仿真结果。
例如构造以基波为主的原始信号,,通过Simulink环境下的DigitalFilterDesign(数字滤波器设计)模块导入2.1.2中FDATool所设计的滤波器文件1.fda。
仿真图和滤波效果图如图2所示。
图2Simulink仿真图及滤波效果图
Fig.2Simulatedconnectionsandwaveform
可以看到经过离散采样、数字滤波后分离出了83.3Hz的频率分量(scope1)。
之所以选取上面的叠加信号作为原始信号,是由于在实际工作中是要对已经经过差分滤波的信号进一步做带通滤波,信号的各分量基本同一致,可以反映实际的情况。
本例设计的滤波器已在实际工作中应用,取得了不错的效果。
4结论
利用MATLAB的强大运算功能,基于MATLAB信号处理工具箱(SignalProcessingToolbox)的数字滤波器设计法可以快速有效的设计由软件组成的常规数字滤波器,设计方便、快捷,极大的减轻了工作量。
在设计过程中可以对比滤波器特性,随时更改参数,以达到滤波器设计的最优化。
利用MATLAB设计数字滤波器在电力系统二次信号处理软件和微机保护中,有着广泛的应用前景。
通带截止频率fp=1KHz,阻带截止频率fs=1.5KHz,通带衰减Rp=1.5dB,阻带衰减Rs=25dB,采样频率Fs=10kHz。
wp=1k/10k
ws=1.5/10
rp=1.5db
rs=25
[N,Wn]=BUTTORD(wp,ws,rp,rs)
[b,a]=butter(n,wn)
freqz(b,a)