零相位数字滤波器在非平稳信号处理中的应用

零相位数字滤波器在非平稳信号处理中的应用

数字滤波是数字数据处理的一般方法。普通数字滤波在滤波过程中会发生一定的相位位移。为了解决这个问题，引入了零相位数字滤波。两个主要的零相位数字滤波算法（之前的序列、欺诈和结果序列）和rr（不安全、不安全、结果序列和结果序列）。首先，我们研究了采用合适的滤波初期模型、信号延伸度、frr和rr的方法，改善了零相位数字滤波的边界问题，并对稳定信号滤波产生了良好的效果。在对非平稳信号进行滤波时,与小波分解重构和经验模态分解等无相移的非平稳信号处理方法相比,零相位数字滤波器具有截止频率明确、计算量小的优点.尤其在计算能力较弱的嵌入式应用系统中,零相位数字滤波器的这一优点更加突出.调幅信号、调频信号是两类典型的非平稳信号,很多工程实际信号可以借助这两种信号进行模拟.例如,电力系统的次同步振荡、电压暂态振荡、开关设备动作时的振动、发电机组轴系的扭振等就可用调幅信号模拟;旋转机械调速过程中的转速则多为调频信号.文献利用线性调频信号对零相位数字滤波器进行了分析.因此,本文作者将主要针对调幅信号来讨论零相位数字滤波器在非平稳信号处理中的特点.利用调幅信号对一种典型的带通零相位数字滤波器进行了仿真分析;对小波包分解重构、经验模态分解及零相位数字滤波器方法的滤波效果进行了比较;探讨了在调幅信号处理中存在的过渡过程,及其对零相位数字滤波的幅值特性以及频率特性的影响;论述了影响滤波误差的因素;提出了利用分段零相位滤波的方法改善零相位数字滤波器过渡过程的方法;用实测的次同步振荡信号对两种方法进行了验证.1零相位滤波、ir数字滤波器对于一维时间序列,零相位数字滤波器的基本原理是,利用正向时间序列和翻转时间序列通过滤波器时的相移相互抵消,从而实现滤波结果的零相移.利用傅里叶变换的时域位移性质和时间翻转性质,即可从理论上论证其正确性.以FRR方法为例,将长度为N+1的单边时间序列x[n]通过冲激响应为h[n]的滤波器,得到一次滤波序列y1[n],再将一次滤波序列进行时序反转、平移得到数据序列y2[n],y2[n]通过滤波器h[n]得到二次滤波序列y3[n],最后将二次滤波序列进行时序反转、平移得到最终滤波输出y[n].从而,有y1[n]=x[n]*h[n](1)y2[n]=y1[Ν-n](2)y3[n]=y2[n]*h[n](3)y[n]=y3[Ν-n](4)设x[n]↔X(ω)、y1[n]↔Y1(ω)、y2[n]↔Y2(ω)、y3[n]↔Y3(ω)、y[n]↔Y(ω)、h[n]↔H(ω)构成傅里叶变换对,且考虑到对实序列h[n]有H(-ω)=H*(ω).依式(1～4),FRR滤波输出的频域表达如下式推导Y(ω)=e-jωΝY3(-ω)=e-jωΝY2(-ω)Η(-ω)=e-jωΝejωΝY1(ω)Η(-ω)=Y1(ω)Η(-ω)=X(ω)Η(ω)Η(-ω)=X(ω)|Η(ω)|2(5)因而,FRR方法实现了零相移滤波,其传递函数频域表示如下ΗFRR(ω)=Η(ω)Η(-ω)(6)同理可得,RRF方法实现零相移滤波时,其传递函数频域表示如下ΗRRF(ω)=Η(-ω)Η(ω)(7)选定冲激响应序列为h[n]的普通数字滤波器,实际使用时可以通过两个分立的滤波过程来实现零相位滤波,也可以直接构造零相位数字滤波器,其冲激响应序列分别为hFRR[n]和hRRF[n].hFRR[n]=h[n]*h[-n](8)hRRF[n]=h[-n]*h[n](9)其中h[-n]为h[n]的时间翻转序列.h[n]的选择可有FIR滤波器和IIR滤波器两种类型.在零相位滤波器的应用中,IIR滤波器在相频特性上的劣势可以忽略,且在相同阶数下IIR滤波器比FIR滤波器具有更好的截止特性,因而h[n]一般选用IIR数字滤波器.23khs.自适应变化为比较小波包分解重构、经验模态分解和零相位数字滤波器3种方法的差异.选取含有3个衰减振荡分量的输入信号f(t),信号采样率为1kHz.设定3个分量的振荡频率分别为10Hz、23Hz、39Hz,其计算公式如下f1(t)=e-3(t-1)u(t-1)sin(20π(t-1))(10)f2(t)=e-2(t-1)u(t-1)sin(46π(t-1))(11)f3(t)=e-(t-1)u(t-1)sin(78π(t-1))(12)输入信号f(t)的计算公式如下f(t)=f1(t)+f2(t)+f3(t)(13)1小波包的滤波、并存在重要差异利用小波包可将信号分解到不同节点,各节点对应了不同频带的子波,对目标节点系数进行重构即可得到所需的滤波输出.考虑到重构信号的不失真和光滑,采用了bior、sym和db小波基对多分量信号进行处理,最终选定db8小波基对f(t)进行5层分解重构,由于信号采样率为1kHz,小波包5层分解的各子波频带宽度为15.625Hz,使10Hz、23Hz、39Hz3个频率成分位于、、节点对应的频带内,滤波结果如图1所示(图中虚线为滤波输出,实线为输入分量).小波包分解重构的滤波输出,突变点附近存在过渡过程,其余部分信号不光滑,在幅值和相位上都存在抖动,其中23Hz和39Hz分量抖动比较严重.2mt经验模态分解mt模式,m经验模态分解方法是自适应的将信号分解为不同时间尺度的本征模态函数(IMF)分量和趋势项,本征模态函数是一类零均值的调幅调频信号.对含随机宽带噪声的信号,经验模态分解是类小波的二进滤波器组,第一个本征模态函数的频带最高.由于m(t)由3个调幅正弦信号组成,经验模态分解得到的前3个本征模态函数分量即对应其滤波输出.其中,第一个本征模态函数对应39Hz分量、第二个本征模态函数对应23Hz分量、第三个本征模态函数对应10Hz分量,滤波结果如图2所示(图中虚线为滤波输出,实线为输入分量).经验模态分解方法的滤波输出,突变点附近存在过渡过程,但其余部分信号很光滑,幅值和相位的稳定性强于小波包分解重构方法,但是仍然存在抖动.我们又采用标准差为0.005、聚合次数为50次的聚合经验模态分解对f(t)进行滤波,发现能改善过渡过程.但是,同时造成了39Hz分量混叠到23Hz分量中,为减少混叠效应需要增加聚合次数,极大的增加了计算量.3过渡过程中的误差检验数字滤波器在对有限长信号进行滤波时,相当于给信号加了一个矩形窗,在滤波输出信号上造成衰减振荡.零相位数字滤波器进行了两次滤波,在输出信号的首尾两端形成失真.FRR方法第二次滤波的首端是输入信号的尾端,滤波输出的首端失真较小;RRF方法第二次滤波的首端是输入信号的首端,滤波输出的尾端失真较小.对于本文考虑的调幅正弦信号,通过在输入信号首尾两端进行信号延拓的方法,来消除端点效应.依据对端点效应误差的要求,延拓的数据长度应比拟于过渡过程时间,本文取1s长度的数据点.消除端点效应的FRR方法和RRF,在使用中的差异可忽略不计.本文的后续讨论中,零相位数字滤波器均采用FRR方法实现,且采用了首尾端的端点延拓,不再特别说明.综合考虑滤波器阶跃响应的过渡过程时间和振荡次数等因素,选定用4阶巴特沃斯带通滤波器(直接型系数见表1)构造零相位数字滤波器,中心频率分别为10Hz、23Hz、39Hz,带宽为3Hz,记为hbp10、hbp23、hbp39.滤波结果见图3(图中虚线为滤波输出,实线为输入分量).滤波输出很光滑,除了突变点附近,幅值和相位与原始信号重合,效果好于小波包分解重构和经验模态分解方法.不考虑信号的突变,则零相位数字滤波器实现了较严格的零相移滤波.3种方法都存在过渡过程.以10Hz分量的误差为例来考察过渡过程,将3种方法所得结果前3s的误差示于图4中.从上至下依次为小波包分解重构、经验模态分解、零相位滤波器误差,实线为原始误差,虚线为幅值放大100倍后的误差.其中,经验模态分解误差极值最大,零相位滤波器与小波包分解重构的误差极值较小;零相位滤波与经验模态分解误差稳定时间较短,小波包分解重构的误差稳定时间较长;但零相位滤波器反时序的特性造成突变点前有信号的提前输出.下面我们将对零相位滤波器的过渡过程详细探讨,并研究其幅值和频率特性,从而对其过渡过程性能进行改善.3过渡过程3.1同上转台滤波过滤算法零相位数字滤波器解决了滤波中存在的相移问题,但是在对非平稳信号进行滤波时,仍需注意滤波中可能存在的过渡过程问题.尤其是在信号发生突变的位置影响较大,这是因为突变点的存在类似于将信号在该点截断,引入了新的端点.考察一个在第2秒时刻幅值由0跃变为1在第8秒再跳变为0的10Hz正弦信号x(t),其计算公式如下x(t)=sin(20π(t-2))[u(t-2)-u(t-8)](14)其中u(t)为单位阶跃信号.将x(t)通过hbp10进行滤波,滤波结果如图5所示,存在过渡过程.为了更清楚地表达过渡过程,将其幅值响应曲线示于图6.输入信号的突变在滤波输出的对应位置上存在过渡过程,稳定到±2%的时间约0.48s,最大超调量为3.4%;进入平稳阶段后,滤波输出与输入信号能很好的吻合,误差均值在10-5量级,而误差的方差在10-7量级,平稳阶段时间越长,误差越小.对这类幅值发生阶跃的正弦信号,可以采用滤波输出信号最大幅值的中值时刻来近似定位突变点,从而去除突变点前的虚假信号.采用该方法估计x(t)的突变点的估计误差只有0.005s,实际应用中受波形形态及噪声的影响,误差变大.事实上,不同特点的信号通过同一零相位数字滤波器,滤波信号幅值和频率表现出的特性也会不同.为进行具体的仿真分析,选取衰减振荡、发散振荡、发散振荡经采取措施后又衰减振荡这3种情况进行研究,利用幅值衰减、幅值发散、幅值先发散后衰减这3种不同调幅正弦信号进行模拟.3.23滤波器hbp10滤波特性先构造调幅正弦信号的幅值包络a1(t)、a2(t)、a3(t),其中a1(t)为衰减幅值、a2(t)为发散幅值、a3(t)为先发散后衰减幅值,其计算公式如下a1(t)=e-(t-1)u(t-1)(15)a2(t)=e(t-9)[1-u(t-9)](16)a3(t)=e(t-5)[1-u(t-5)]+e-(t-5)u(t-5)(17)对应的调幅正弦信号分别为m1(t)、m2(t)、m3(t),其计算公式如下m1(t)=a1(t)sin(20πt)(18)m2(t)=a2(t)sin(20πt)(19)m3(t)=a3(t)sin(20πt)(20)采用一个高斯白噪声信号wn(t)作为噪声干扰信号,wn(t)的统计分布为N(0,0.01).对应的含噪声的输入信号分别为s1(t)、s2(t)、s3(t),其表达式分别为下列公式所示s1(t)=m1(t)+wn(t)(21)s2(t)=m2(t)+wn(t)(22)s3(t)=m3(t)+wn(t)(23)s1(t)、s2(t)、s3(t)的波形如图7所示,经过零相位带通数字滤波器hbp10滤波后,滤波输出信号如图8所示.可明显看出,这3种信号的滤波输出的最大幅值小于的x(t)的滤波输出的最大幅值;过渡过程时间也小于x(t)的过渡过程时间;其中s3(t)的滤波结果最好,误差均值和方差都在10-5量级.可见输入信号的形态特点对滤波结果有影响.4振幅值和频率特性4.1滤波结果的频率方差滤波输出的幅值和频率可利用希尔伯特变换求得.其基本原理如下.对信号ν(t),其希尔伯特变换由下式定义ν˜(t)=1π∫-∞+∞ν(τ)t-τdτ(24)对应的幅值包络a(t)可按下式计算a(t)=ν2(t)+ν˜2(t)(25)对应的相位θ(t)可按下式计算θ(t)=arctanν˜(t)ν(t)(26)对应的瞬时频率定义为f(t)=12πdθ(t)d(t)(27)下面利用希尔伯特变换获取滤波结果的幅值包络和频率,对零相位数字滤波器的幅值响应和频率响应特点进行详细分析.4.2放大电路一般特性s1(t)、s2(t)、s3(t)经过零相位带通数字滤波器hbp10滤波后,理想的幅值应为a1(t)、a2(t)、a3(t).实际按式(25)计算出的滤波输出信号幅值分别如图9～图11所示.可明显看出,s1(t)、s2(t)、s3(t)的滤波输出的最大幅值分别只有0.81、0.80和0.93,明显小于的x(t)的滤波输出的最大幅值1.034;s3(t)的最大值没有明显延迟,s1(t)、s2(t)过渡过程时间为0.22s、0.22s,也明显小于对x(t)滤波时的过渡过程时间0.48s;信号幅值较小的部分噪声影响比较大.由于不像x(t)有稳定的信号阶段可以使滤波输出完成过渡达到稳定,上述3类信号衰减系数越大幅值损失越大.将不同衰减系数的10Hz分量滤波输出的幅值示于图12中,其结果印证了该结论.对于零相位带通数字滤波器而言,与普通带通数字滤波器相似,其幅值响应具有如下特点:1)滤波器类型一定时,较宽的频带可以获得较大的幅值和较小的时延;提高滤波器阶数可加速阻频带衰减,减少临近频带的干扰,但也会带来通带幅值的振荡.2)输入信号幅值包络变化时,滤波输出的时延从幅值包络变化率大的部分偏向变化率小的部分;幅值包络变化越平缓,滤波后的时延也越小;滤波器阶数大于等于4时,信号幅值包络越接近阶跃,超调量越大.3)因为白噪声无法完全滤除,幅值较小的部分噪声影响比较明显.4.3滤波、滤波效率及信号背景分析按式(27)式定义对m1(t)、m2(t)、m3(t)的滤波输出信号进行计算,得到的瞬时频率分别示于图13～图15中.其中,首尾两端的振荡是由于希尔伯特变换的边界效应造成的.另外需要注意的地方是,在信号突变点附近(约为一个过渡过程时间),滤波输出存在较大的频率突变,说明滤波输出信号在该位置存在畸变,在图13～图14中分别位于0.755s、9.246s处.取中间频率稳定的3～8s时间段进行分析,可得频率均值为10.000Hz,方差小于10-5量级.其中,m3(t)在中间突变点处的频率处于9.99Hz与10.01Hz之间,明显好于m1(t)、m2(t)在突变点处的滤波结果.为考察噪声对频率响应的影响,取s1(t)信号为例.图16给出了s1(t)经滤波后的瞬时频率曲线,可见滤波结果变差,信号幅值较小的部位上噪声造成相位的波动,使波形相对发生了较大的畸变.仍取3～8s时间段进行分析,可得频率均值为9.962Hz,方差在10-1量级.零相位带通数字滤波器的频率响应具有如下特点:1)中心频率和频带宽度变化对稳定的频率段影响很小,主要对信号突变点附近的频率有影响.2)当输入信号的幅值包络变化时,滤波输出的频率畸变点位置与幅值时延的方向相反,位于突变点偏幅值包络变化率较大的一侧.3)噪声的存在,造成信号的畸变,滤波输出频率存在波动,幅值较小的部分噪声影响比较明显.5分段零相位滤波由上述仿真结果可知,m3(t)信号的滤波结果明显好于m1(t)与m2(t).m3(t)与m1(t)、m2(t)的不同之处在于,其突变点两边的幅值是对称的.基于此,作者提出分段零相位滤波方法,用以改善零相位滤波器在处理存在突变的调幅信号时的过渡过程性能,其具体步骤如下:1)对信号进行零相位滤波.2)在滤波幅值包络的中值点附近,搜索滤波输出中绝对值最小的点,进行突变点的近似定位.3)将输入信号按定位结果分段,对分段后的信号进行端点延拓,再分别进行零相位滤波.4)将分段滤波结果组合,形成滤波输出.零相位滤波器采用hbp10,利用零相位滤波和分段零相位滤波方法分别对s1(t)进行滤波,滤波结果示于图17中.从图中可以明显观察到,分段零相位滤波的结果不但比原始的零相位滤波结果过渡过程缩短,而且也解决了反时序提前输出的问题.可见,分段零相位滤波能缩短过渡过程时间,消除反时序提前输出,极大改善零相位滤波器的滤波效果.其缺点是,增加了滤波次数和端点延拓次数,增加了计算量.6零相位数字滤波器和分段零相位数字滤波器为从次同步振荡中分离出主要的振荡分量,将该方法应用到次同步振荡信号的滤波中.在某次试验中测得振荡信号,采样率为1kHz,数据长度为95s,如图18所示.通过时频分析知该振荡信号中主要包含有13Hz、50Hz、25Hz等几个较大的频率成分,其中在次同步频域内13Hz振幅最大.以13Hz分量为例,仍采用4阶巴特沃斯带通滤波器(直接型的中心频率为13Hz,带宽3Hz),构造零相位