基于時變零相位濾波的變轉速滾動軸承故障診斷

陳向民 張 亢 晉風華 李錄平

長沙理工大學能源與動力工程學院，長沙，410076

0 引言

滾動軸承是旋轉機械設備中最為常用的零部件之一，其工作狀態正常與否直接影響到整臺設備的性能。當滾動軸承出現點蝕、裂紋等局部故障時，故障元件的固有頻率會被故障元件的通過頻率所調制[1]。有效解調出滾動軸承故障信號中的調制頻率是準確診斷滾動軸承故障的關鍵。

變轉速下機械設備的振動信號中往往包含了豐富的狀態信息，一些在平穩運行時不易體現的故障特征可能會被充分表現出來[2-3]。當滾動軸承處于變轉速下運行時，其調制信號中的調制頻率與轉速相關，也為一時變量，直接對解調信號進行頻譜分析會產生“頻率混疊”現象[4]。

階次跟蹤是工程上常用的一種變轉速工況下的信號分析方法。它根據參考軸的轉速對信號進行等角度重采樣，能提取振動信號中與轉速有關的信息[5-6]。但由于機械設備中的零部件眾多，其振動信號中與轉速相關的調制信號成分也較多，且各種調制信號成分之間的能量強弱不同，特別對微弱故障而言，其故障特征容易被其他特征所淹沒，因此，需預先對調制信號中的各調制成分進行分離。

常用的自適應信號分離方法主要有EMD(empirical mode decomposition)[7-8]、EEMD(ensemble empirical mode decomposition)[9-10]、LMD(local mean decomposition)[11-12]，VMD(variational mode decomposition)[13]等方法，這些方法在滾動軸承故障特征提取中均取得了一定的效果。但當兩調制信號的頻率比較接近時，卻無法有效分離。對于多調制信號的分離問題，文獻[14]提出了基于線調頻小波路徑追蹤(chirplet path pursuit，CPP)的自適應時變濾波方法，但該方法濾取的信號存在時延。為此，文獻[15]提出采用零相位濾波，有效解決了時延問題，并將其應用于載波頻率(即齒輪嚙合頻率)和調制頻率(即軸轉頻)的分析中。但對于滾動軸承來說，其載波頻率(即共振頻率)為一定值，因此，不能直接使用該方法對變轉速下的滾動軸承故障振動信號進行分析，且當載波頻率被多個調制頻率同時調制時也無法有效分離。

綜上所述，針對變轉速下滾動軸承故障特征的提取與分離，本文將CPP算法[16-18]與零相位濾波相結合，并將其用于變轉速下的滾動軸承故障特征的提取與分離。對變轉速下的滾動軸承故障振動信號進行算法仿真和應用實例分析，結果表明，本文方法可在強噪聲環境下有效分離變轉速滾動軸承的故障調制特征。

1 線調頻小波路徑追蹤算法

線調頻小波路徑追蹤(CPP)算法采用的多尺度線調頻基元函數庫為

D(haμbμI)={haμbμI(t)}=
{KaμbμIe[-i(aμt+bμt2)]LI(t)}

(1)

其中，D表示基元函數庫；haμbμI(t)為多尺度線調頻基元函數；I表示動態分析時間段，I=[kN2-j，(k+1)N2-j]；N為分析信號長度；k為動態時間段的序號，k=0,1,…,2j-1；j為分析尺度系數，j=0，1，…，lb(N-1)；KaμbμI為歸一化系數，滿足‖haμbμI‖=1；aμ和bμ分別為頻率偏置系數和調頻率，且滿足aμ+2bμt

采用多尺度線性調頻基元函數對分析信號進行逐段投影，并計算每個時間分析段I內的投影系數和對應的線調頻基元函數。分析信號與多尺度線性調頻基元函數的相似性越高，其投影系數也就越大，此時，基元函數對應的能量也就越大。因此，需尋求一種動態連接算法Π，使得所連接基元函數對應的信號在整個分析時間內的總能量最大，且連接算法Π應覆蓋整個分析時間段，不重疊，即

(2)

式中，1I(t)為矩形濾波器，數字1表示矩形濾波器的幅值為常量1。

此時，對應的投影系數集合和基元函數集合分別為

(3)

CPP方法中連接算法Π的連接步驟如下：

(1)初始化。以i表示分析時間段序號，di表示第i個時間段之前的分解信號，ei表示第i個時間段的分解信號，Edi表示第i個分析時間段之前分解信號的總能量，li表示連接到第i個分析時間段的前一分析時間段序號，Eei表示第i個分析時間段投影系數對應的分解信號的能量。初始化時，Edi=0,li=0。

(2)對于動態分析時間段集合I={I1,I2…}中的每一個元素Ii，查找出與其相鄰的所有下一個動態分析時間段集合{Ij},如果

Edi+Eei>Edj

(4)

則有

(5)

連接算法Π可保證在整個分析時間段內基元函數組合所對應的信號與分析信號最為相似。基元函數在動態分析時間支持區Ii內的瞬時頻率fIi(t)=aμ+2bμt,ti∈Ii，將所有動態時間段集合I={I1,I2…}中所對應的頻率曲線集合fI={fI1,fI2…}按時間先后順序連接成線則為信號在整個分析時間段內的瞬時頻率估計。

2 時變零相位濾波器

零相位濾波是一種濾波思想，不涉及具體的濾波器，可以是FIR濾波器，也可以是IIR濾波器。零相位濾波中的FRR(forward-reverse filtering, reverse output)法是先將輸入信號x(t)按順序進行濾波，然后將所得結果逆轉后反向通過濾波器，再將結果逆轉后輸出，便可得到濾波信號y(t)[19]，其頻域表達式為

(6)

式中，X(ejω)為輸入信號x(t)的頻譜；Y(ejω)為輸出信號的頻譜；H(ejω)為濾波器的幅頻特性。

由式(6)可得

Y(ejω)=X(ejω)|H(ejω)|2

(7)

由式(7)可看出，輸出信號與輸入信號之間僅在幅值上相差|H(ejω)|2，相位上并無時延或畸變，即實現了零相位濾波。

時變零相位濾波器可由經典零相位濾波方法設計而來，即根據信號的頻率變化，在每一時刻均設計一個零相位濾波器[14]，即

H(ejω,ti)=Hi(ejω)i=1,2,…,N

(8)

假定分析信號為x(ti),i=1,2，…，N，CPP算法估計出的頻率為f(ti)，調制信號的階次為O，則時變零相位濾波器的設計步驟如下：

(9)

(10)

時變零相位濾波方法可根據信號的頻率變化特點自適應地進行濾波，因而非常適合頻率變化下的信號分析。

3 變轉速下的滾動軸承故障診斷原理

當齒輪箱中的滾動軸承出現點蝕、裂紋等局部故障時，故障元件的固有頻率會被故障元件的通過頻率所調制。滾動軸承內外圈局部故障的故障特征頻率計算公式分別如下[1]：

(11)

(12)

式中，Z為滾動體個數；fr為軸的轉頻；d0為滾動體直徑；D為節徑；α為接觸角。

齒輪箱的轉速發生改變，軸承的故障特征頻率也會隨之發生改變，此時，基于平穩分析的包絡譜分析不能很好地表征軸承的故障特征。但對于同一軸承，無論轉速如何變化，其故障特征階次卻是一定值(即只與軸承本身的參數相關)。內外圈的故障特征階次計算公式分別為

(13)

(14)

因此，本文針對變轉速下滾動軸承故障調制信息的提取與分離，將CPP算法和零相位濾波器設計相結合，提出了基于時變零相位濾波的變轉速滾動軸承故障診斷方法。該方法先采用CPP算法從故障軸承信號中估計齒輪箱中的嚙合頻率，并根據嚙合頻率與齒輪齒數求取轉速，同時，采用Hilbert包絡解調方法求取軸承故障信號的包絡信號；然后根據轉速信號分別設計各故障特征階次處的時變零相位濾波器，再采用各時變零相位濾波器對包絡信號進行分析，便可得到各階調制信號；最后根據轉速信號對各階調制信號進行階次分析，并根據階次譜診斷軸承故障。本文算法流程如圖1所示。

圖1 算法流程框圖Fig.1 Algorithm flow chart

4 算法仿真

為驗證本文方法分離軸承故障調制特征的有效性，設置仿真合成信號s為

s=sbe+sge+sno

(15)

(16)

sge=1.8[1+0.5cos(6πt3+10.8t2+50πt)]·
cos(108πt3+194.4πt2+900πt)

(17)

其中，sbe為軸承故障信號，軸承的固有頻率為1400 Hz，其三階故障特征階次分別為1.6、3.2和4.8，軸承故障信號分量如圖2a所示；sge為齒輪信號分量，其中信號的幅值為1.8，模擬齒數為18，嚙合頻率被1倍轉頻調制，調制信號幅值為0.5，齒輪信號分量如圖2b所示；sno為強度為-6 dB的高斯白噪聲。合成信號s的時域波形如圖2c所示。

(a)軸承故障信號分量

(b)齒輪信號分量

(c)合成信號圖2 仿真信號分量及合成信號Fig.2 Simulation signal components and its composite signal

采用Hilbert解調分析求取合成信號的包絡，得到的包絡信號如圖3所示。

圖3 仿真合成信號的包絡信號Fig.3 Envelope signal of the composite signal

設置頻率偏置系數的搜尋范圍為430～530 Hz，搜尋分辨率為1 Hz；調頻率的搜尋范圍為-200～200 Hz/s，搜尋分辨率為1 Hz/s(需要說明的是，對于頻率偏置系數和調頻率的取值，可采用短時傅里葉變換等時頻分析方法大致估計信號的頻率變化范圍，再適當將變化范圍擴大作為取值范圍，其取值范圍的大小只影響計算時間，對計算結果無影響)。根據上述參數采用CPP方法對合成信號進行分析，估計出的齒輪嚙合頻率如圖4中的虛線所示。圖4所示的實線為實際的嚙合頻率。

圖4 估計的嚙合頻率和實際嚙合頻率對比Fig.4 Comparison of estimated meshing frequency and actual meshing frequency

根據圖4中估計的嚙合頻率設計前3階故障特征階次處的時變零相位濾波器，得到的3個時變零相位濾波器的時頻特性如圖5所示，圖中淺色區域為時變零相位濾波器的通帶。圖5a～圖5c分別對應軸承3個調制頻率處時變濾波器的時頻特性。

(a)一階時變零相位濾波器

(b)二階時變零相位濾波器

(c)三階時變零相位濾波器圖5 各時變零相位濾波器的時頻特性Fig.5 Time-frequency characteristic of each time-varying zero-phase filter

根據圖5中的3個時變零相位濾波器分別對圖3所示包絡信號進行時頻濾波，得到的3個信號分量分別如圖6～圖8所示，其中，實線為原始信號，虛線為時變濾波獲取的信號。圖7b和圖8b分別為圖7a和圖8a在時間段0.1～0.3 s的局部放大圖。對比圖6、圖7b、圖8b中的實線和虛線可知，時變濾波信號與原始信號僅在幅值上存在一定的差異，而相位上無畸變。

圖6 一階濾波信號Fig.6 The first order filtered signal

(a)二階濾波信號

(b)二階濾波信號局部放大圖圖7 二階濾波信號及其局部放大圖Fig.7 The second order filtered signal and itspartial enlargement

(a)三階濾波信號

(b)三階濾波信號局部放大圖圖8 三階濾波信號及其局部放大圖Fig.8 The third order filtered signal and its partial enlargement

對提取的三個濾波信號分別進行階次譜分析，結果如圖9所示。圖9a～圖9c中，分別在階次1.591、3.218、4.809處出現了顯著峰值，與信號設置的階次1.6、3.2、4.8相符，驗證了本文方法提取和分離軸承故障調制特征的有效性。

(a)一階濾波信號的階次譜

(b)二階濾波信號的階次譜

(c)三階濾波信號的階次譜圖9 濾波信號的階次譜Fig.9 Order spectrums of the filtered signals

為增加對比，采用圖4中的估計嚙合頻率直接對圖3中的包絡信號進行階次分析，得到的階次譜如圖10所示。圖10中，在軸承故障階次1.591和3.218處出現了峰值，但在三階軸承故障階次4.8處未出現峰值，且在轉頻階次0.998處出現了明顯峰值。

圖10 仿真合成信號的包絡階次譜Fig.10 Envelope order spectrum of the simulation composite signal

同時，采用EEMD方法對圖3所示包絡信號進行分解，并對分解后的IMF進行階次分析，以獲取各IMF的階次譜。圖11a、圖11b分別為第4和第5個IMF的階次譜，圖中，在軸承故障階次1.591、3.181、4.809處分別出現了峰值，但在非軸承故障調制階次3.921處出現了未知峰值。

(a)IMF4的階次譜

(b)IMF5的階次譜圖11 IMF4與IMF5的階次譜Fig.11 Order spectrum of IMF4 and IMF5

5 應用實例

為驗證本文方法在實測變轉速軸承故障振動信號中分離軸承故障調制特征的有效性，在齒輪箱上進行變轉速軸承故障試驗。試驗齒輪箱主動齒與從動齒的齒數均為37。試驗軸承為6307E型滾動軸承，軸承參數如表1所示。將表1中參數代入式(13)、式(14)計算可得，軸承外圈局部故障和內圈局部故障的一階故障調制特征階次分別為3.06和4.93。

表1 滾動軸承參數

為模擬軸承故障，利用激光分別在軸承外圈和內圈上切割寬為0.15 mm、深為0.13 mm的槽，以模擬外圈和內圈局部故障。軸承故障樣本如圖12所示。

(a)軸承外圈故障 (b)軸承內圈故障圖12 軸承故障樣本Fig.12 Faulted bearings sample

試驗采用調速電機使軸承在變轉速下運行。試驗時，采用LMS采集儀同步采集振動加速度信號和轉速信號，以便進行對比分析。

5.1 變轉速軸承外圈局部故障振動信號分析

對外圈故障滾動軸承進行試驗，采樣頻率為8192 Hz，采樣時長為2 s。軸承外圈故障振動信號的時域波形如圖13所示。圖14 為圖13的包絡信號。

圖13 變轉速軸承外圈故障振動信號Fig.13 Vibration signal of bearing with outer-race failure under variable rotational speed

圖14 變轉速軸承外圈故障的包絡信號Fig.14 Envelope signal of bearing with outer-race failure under variable rotational speed

采用CPP方法對軸承外圈故障信號進行分析，估計齒輪的嚙合頻率如圖15中的虛線所示。圖15中的實線為實測的嚙合頻率。

圖15 估計的嚙合頻率與實測嚙合頻率對比Fig.15 Comparison of estimated meshing frequency and actual meshing frequency

采用本文方法對變轉速軸承外圈故障振動信號進行分析，得到的階次譜如圖16所示，圖16a～圖16c中，分別在階次3.086，6.171，9.257處出現了顯著的峰值，與軸承外圈故障的1～3階故障階次吻合，表明軸承外圈出現了局部故障。

根據圖15中估計的嚙合頻率信息對圖14所示包絡信號直接進行階次分析，得到包絡信號的階次譜如圖17所示。圖17中出現了較多的峰值，但這些峰值均與軸承外圈故障的故障階次有一定差異。

(a)一階濾波信號的階次譜

(b)二階濾波信號的階次譜

(c)三階濾波信號的階次譜圖16 軸承外圈故障振動信號的時變濾波分析結果Fig.16 Analysis result of the time-varying filtered vibration signal of bearing with outer-race failure

圖17 變轉速軸承外圈故障包絡信號的階次譜Fig.17 Order spectrum of envelope signal of bearing with outer-race failure under variable rotational speed

采用EEMD方法對圖14所示包絡信號進行分解，并根據圖15中估計的嚙合頻率曲線對分解后的IMF進行階次分析，得到的階次譜如圖18所示。圖18a～圖18c分別為IMF6～IMF8的階次譜，在階次3.086、6.171、9.171處出現峰值，與軸承外圈故障特征階次基本相符。但在轉頻階次1.029、1.971處也出現峰值，且二階和三階軸承故障特征階次6.171、9.171同處于分量IMF6中，未能有效分離。

5.2 變轉速軸承內圈局部故障振動信號分析

對內圈局部故障滾動軸承進行試驗。采樣頻率為20480 Hz，采樣時長為0.8 s。圖19為軸承內圈故障振動信號的時域波形圖。

(a)IMF6的階次譜

(b)IMF7的階次譜

(c)IMF8的階次譜圖18 變轉速軸承外圈故障IMF6～IMF8的階次譜Fig.18 Order spectrum of IMF6～IMF8 of bearing with outer-race failure under variable rotational speed

圖19 變轉速軸承內圈故障振動信號Fig.19 Vibration signal of bearing withinner-racefailure under variable rotational speed

(a)一階濾波信號的階次譜

(b)二階濾波信號的階次譜圖20 軸承內圈故障振動信號的時變濾波分析結果Fig.20 Analysis result of the time-varying filtered vibration signal of bearing with inner-race failure

采用本文的時變濾波方法對圖19信號進行分析，分析結果如圖20所示。圖中，在階次4.982和階次9.817處存在顯著峰值，與軸承內圈局部故障的一階和二階特征階次相符，故可判斷軸承內圈出現了局部故障。

同時，直接采用包絡階次分析方法對圖19所示信號進行分析，得到的階次譜如圖21所示。圖21中，在階次4.982處存在峰值，與軸承內圈局部故障的一階特征階次相符，但二階特征階次處無顯著峰值，且階次譜中存在眾多的干擾階次。

圖21 變轉速軸承內圈故障包絡信號的階次譜Fig.21 Order spectrum of envelope signal of bearing with outer-race failure under variable rotational speed

采用基于EEMD的包絡階次分析方法對圖19所示信號進行分析，分析結果如圖22所示。圖22中，在階次9.894和4.947處存在峰值，與軸承內圈故障的一、二階特征階次相符。但與圖20相比，圖22中存在較多的噪聲干擾，分析效果要遜色于圖20分析結果。

(a)IMF6的階次譜

(b)IMF7的階次譜圖22 變轉速軸承內圈故障IMF6～IMF7的階次譜Fig.22 Order spectrum of IMF6～IMF7 of bearing with outer-race failure under variable rotational speed

6 結論

(1)時變零相位濾波器根據調制信號的頻率變化特點，能自適應地濾取調制信號，無相位畸變，且具有較強的抗噪能力。

(2)與EEMD等方法進行了對比研究，結果表明，時變零相位濾波方法不僅能有效提取軸承故障調制特征，而且可實現其各調制特征的有效分離。

(3)由于濾波器在通帶內的幅頻特性會有一定的波動，因此，濾取的調制信號的幅值較原始信號會有一定的幅值畸變，如何減少濾波后信號幅值的畸變可進一步研究。

[1] 丁康，李巍華，朱小勇. 齒輪及齒輪箱故障診斷實用技術[M]. 北京：機械工業出版社，2005.

DING Kang，LI Weihua，ZHU Xiaoyong. The Gear and Gearbox Fault Diagnosis Practical Technology[M]. Beijing：China Machine Press，2005.

[2] 李志農，丁啟全，吳昭同，等. 旋轉機械升降速過程中的雙譜-FHMM識別方法[J]. 振動工程學報，2003，16(2): 171-174.

LI Zhinong, DING Qiquan, WU Zhaotong, et al. Study on Bispectrum-FHMM Recognition Method in Speed-up and Speed-down Process of Rotating Machinery[J]. Journal of Vibration Engineering, 2003,16(2):171-174.

[3] 張亢，程軍圣. 基于LMD和階次跟蹤分析的滾動軸承故障診斷[J]. 振動、測試與診斷，2016，36(3)：586-591.

ZHANG Kang, CHENG Junsheng. A Fault Diagnosis Method for Rolling Bearing Based on LMD and Order Tracking Analysis[J]. Journal of Vibration, Measurement & Diagnosis,2016, 36(3): 586-591.

[4] 康海英，欒軍英，鄭海起，等. 基于階次跟蹤和經驗模態分解的滾動軸承包絡解調分析[J]. 機械工程學報，2007，43(8)：119-122.

KANG Haiying, LUAN Junying, ZHENG Haiqi, et al. Envelope Demodulation Analysis of Bearing Based on Order Tracking and Empirical Mode Decomposition[J]. Chinese Journal of Mechanical Engineering, 2007, 43(8): 119-122.

[5] 郭瑜，秦樹人，湯寶平，等. 基于瞬時頻率估計的旋轉機械階比跟蹤[J]. 機械工程學報，2003，39(3): 32-36.

GUO Yu, QIN Shuren, TANG Baoping, et al. Order Tracking of Rotating Machinery Based on Instantaneous Frequency Estimation[J]. Chinese Journal of Mechanical Engineering, 2003, 39(3): 32-36.

[6] BORGHESANI P，PENNACCHI P，RANDALL R B，et al. Order Tracking for Discrete-random Separation in Variable Speed Conditions[J].Mechanical Systems and Signal Processing，2012，30：1-22.

[7] 孟宗，閆曉麗，王賽. 基于HMM校正與神經網絡延拓的EMD端點效應抑制方法[J]. 中國機械工程，2015，26(14): 1920-1925.

MENG Zong, YAN Xiaoli, WANG Sai. Restraining Method of End Effect for EMD Based on Error Calibration by HMM and Neural Network[J]. China Mechanical Engineering, 2015, 26(14): 1920-1925.

[8] CHENG Junsheng, YU Dejie, YANG Yu. A Fault Diagnosis Approach for Roller Bearings Based on EMD Method and AR Model[J]. Mechanical System and Signal Processing, 2006, 20: 350-362.

[9] 何青，褚東亮，毛新華. 基于EEMD和MFFOA-SVM的滾動軸承故障診斷[J]. 中國機械工程，2016，27(9): 1191-1197.

HE Qing, CHU Dongliang, MAO Xinhua. Study on Rolling Bearing Fault Diagnosis Based on EEMD and MFFOA-SVM[J]. China Mechanical Engineering, 2016, 27(9): 1191-1197.

[10] LEI Yaguo, HE Zhengjia, ZI Yanyang. Application of the EEMD Method to Rotor Fault Diagnosis of Rotating Machinery[J].Mechanical System and Signal Processing, 2009, 23(4): 1327-1338.

[11] 孟宗，閆曉麗，王亞超. 基于LMD和HMM的旋轉機械故障診斷[J]. 中國機械工程，2014，25(21): 2942-2946.

MENG Zong, YAN Xiaoli, WANG Yachao. Rotating Machinery Fault Diagnosis Based on Local Mean Decomposition and Hidden Markov Model[J]. China Mechanical Engineering, 2014, 25(21): 2942-2946.

[12] CHENG Junsheng， ZHANG Kang, YANG Yu. An Order Tracking Technique for the Gear Fault Diagnosis Using Local Mean Decomposition Method[J]. Mechanism and Machine Theory，2012，55: 67-76.

[13] 唐貴基，王曉龍. 變分模態分解方法及其在滾動軸承早期故障診斷中的應用[J]. 振動工程學報，2016，29(4): 638-648.

TANG Guiji, WANG Xiaolong. Variational Mode Decomposition Method and Its Application on Incipient Fault Diagnosis of Rolling Bearing[J]. Journal of Vibration Engineering, 2016, 29(4): 638-648.

[14] 彭富強，于德介，吳春燕. 基于自適應時變濾波階比跟蹤的齒輪箱故障診斷[J]. 機械工程學報，2012，48(7): 77-85.

PENG Fuqiang, YU Dejie, WU Chunyan. Self-adaptively Time-varying Filter Based Order Tracking Method and Its Application in Gearbox Fault Diagnosis[J]. Chinese Journal of Mechanical Engineering, 2012, 48(7): 77-85.

[15] WU Chunyan, LIU Jian, PENG Fuqiang, et al. Gearbox Fault Diagnosis Using Adaptive Zero Phase Time-varying Filter Based on Multi-scale Chirplet Sparse Signal Decomposition[J]. Chinese Journal of Mechanical Engineering, 2013, 26(4): 831-838.

[17] LUO Jiesi, YU Dejie, LIANG Ming. Gear Fault Detection Under Time-varying Rotating Speed via Joint Application of Multiscale Chirplet Path Pursuit and Multiscale Morphology Analysis[J]. Structural Health Monitoring, 2012, 11(5)：526-537.

[18] 陳向民, 于德介, 任凌志，等. 基于線調頻小波路徑追蹤階比能量解調的齒輪輪齒裂紋故障診斷[J]. 中國機械工程, 2011, 22(21): 2598-2603.

CHEN Xiangmin, YU Dejie, REN Lingzhi, et al. Order Energy Demodulating Approach Based on Chirplet Path Pursuit and Its Application to Gear Tooth Crack Fault Diagnosis[J]. Journal of Vibration Engineering, 2011, 22(21): 2598-2603.

[19] 紀躍波，秦樹人，湯寶平. 零相位數字濾波的方法與實現[J]. 振動、測試與診斷，2000，20(增刊): 167-172.

JI Yuebo, QIN Shuren, TANG Baoping.Method and Realization of Zero Phase Digital Filter[J]. Journal of Vibration, Measurement & Diagnosis, 2000, 20(S): 167-172.