高速旋轉(zhuǎn)機(jī)械用消除頻率混疊現(xiàn)象的小波變換算法研究

楊 鵬

(1.核工業(yè)理化工程研究院，天津 300180；2.國(guó)防科技工業(yè)核材料技術(shù)創(chuàng)新中心，天津 300180)

小波理論被認(rèn)為是對(duì)傅里葉分析的重大突破，已成為當(dāng)今從應(yīng)用數(shù)學(xué)到信號(hào)與圖像處理等眾多領(lǐng)域的研究熱點(diǎn)。小波變換是由法國(guó)科學(xué)家Morlet于1980年在進(jìn)行地震數(shù)據(jù)分析時(shí)提出的。由于小波分析方法良好的時(shí)頻分辨能力，近年來在工程上獲得了廣泛應(yīng)用。但在實(shí)際信號(hào)分析過程中發(fā)現(xiàn)，小波變換可能引起嚴(yán)重的頻率混疊現(xiàn)象。針對(duì)以上問題，本文試圖從原理上分析小波分析產(chǎn)生頻率混疊的原因，并介紹小波分析消除混疊改進(jìn)算法的Matlab實(shí)現(xiàn)過程[1]。

1 小波算法及其頻率混疊

Mallat算法在小波分析中的地位相當(dāng)于快速傅里葉變換(FFT)在傅里葉變換中的地位。Mallat算法又稱塔式算法，它由小波濾波器H、G和h、g對(duì)信號(hào)進(jìn)行分解和重構(gòu)，算法如下。

分解算法為：

(1)

式中：t為離散時(shí)間序列號(hào)，t=1,2,3,…,N；f(t)為原始信號(hào)；j為層數(shù)，j=1,2,3,…,J，J=log2N；H和G為時(shí)域中小波分解濾波器的系數(shù)；Aj為信號(hào)f(t)在第j層的近似部分(即低頻部分)的小波系數(shù)；Dj為信號(hào)f(t)在第j層的細(xì)節(jié)部分(即高頻部分)的小波系數(shù)[2]。

重構(gòu)算法為：

(2)

式中：j為分解的層數(shù)，若分解的最高層即分解的深度為J，則j=J-1,J-2,…,1,0；h和g為時(shí)域中小波重構(gòu)濾波器的系數(shù)。

Mallat算法中包括3個(gè)關(guān)鍵環(huán)節(jié)：與正交鏡像濾波器卷積、隔點(diǎn)采樣及隔點(diǎn)插零。這些環(huán)節(jié)要求正交鏡像濾波器必須具有理想的截止特性，即h和H為理想低通濾波器，g和G為理想高通濾波器。而Mallat算法所用的正交鏡像濾波器實(shí)際上并非理想濾波器[3]，以Daubechies(dbN)小波為例，若信號(hào)的采樣頻率為800 Hz，N=5和N=20時(shí)的dbN小波的頻域特性如圖1所示[4]。由圖1可看出，低通部分及帶通部分在交界處有明顯重疊，重疊的長(zhǎng)度隨N的增大而減少。

2 消除頻率混疊的小波分解與重構(gòu)算法及其Matlab實(shí)現(xiàn)

小波變換的頻域定義相當(dāng)于用1族帶通濾波器對(duì)信號(hào)進(jìn)行濾波，隨著小波變換尺度的減小，濾波器的中心頻率向高頻率移動(dòng)的同時(shí)，其帶通寬度也隨之增加。小波理論提供了包括傅里葉分析所采用的三角基函數(shù)以外的多種小波基函數(shù)，使小波分析具備多種情況下的分析能力，作為一種信號(hào)的時(shí)頻分析方法，具有多分辨 分析的特點(diǎn)，它將信號(hào)分解成低頻部分和高頻部分，分解過程中，低頻部分失去的信息由高頻部分捕獲[5]。在第j+1層的分解中，只分解第j層的低頻部分，保留第j層的高頻部分，逐層分解下去，可進(jìn)行很深層次的分解。在機(jī)械設(shè)備故障診斷中，當(dāng)機(jī)器發(fā)生故障時(shí)，由于機(jī)器各零部件的結(jié)構(gòu)不同，輸出信號(hào)在各頻率波段中的表現(xiàn)也不同。在故障診斷的研究中，信號(hào)采集后，要對(duì)信號(hào)進(jìn)行分析處理，提取并考察其特征參量的變化。

圖1 N=5和N=20時(shí)dbN小波的頻域特性Fig.1 Frequency of dbN with N=5 and 20

基于上述原因，亟需一種消除小波分解和重構(gòu)過程中產(chǎn)生的頻率混疊的算法[6]。改進(jìn)算法的主要思想為：利用傅里葉變換和傅里葉逆變換來去掉多余的頻率成分。基于這一思想的改進(jìn)算法如圖2所示。此方法的Matlab實(shí)現(xiàn)過程[7]如下。

圖2 改進(jìn)的小波分解與重構(gòu)快速算法Fig.2 Improvement algorithm for resolvement and reconstitution of wavelet

2.1 消除取得近似部分Aj過程中頻率混疊的方法

近似部分是能量的低階頻率部分，其混疊會(huì)導(dǎo)致低階信號(hào)混亂，影響信號(hào)基礎(chǔ)頻率的輸出，其消除方法如下：

2) 將FFT結(jié)果中頻率大于1+fs/2j(fs為信號(hào)采樣頻率)部分譜值置零；

3) 對(duì)置零后的結(jié)果進(jìn)行快速傅里葉逆變換；

4) 對(duì)快速傅里葉逆變換的結(jié)果進(jìn)行隔點(diǎn)采樣，采樣后的結(jié)果作為Aj再進(jìn)行下一步分解。

在Matlab中，此過程的實(shí)現(xiàn)通過定義一個(gè)函數(shù)m文件來實(shí)現(xiàn)，設(shè)函數(shù)名為r_appcoef。首先用[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters(wname)語句獲得小波分解所要用到的正交鏡像濾波器，其中Lo_D,Hi_D,Lo_R,Hi_R分別對(duì)應(yīng)G、H、g、h；wfilters(wname)表示所選擇的小波基函數(shù)名[8-9]。

利用fft(convn(A,Lo_D))命令即可實(shí)現(xiàn)步驟1。

值得探討的是步驟2中的頻率問題。從濾波的角度看，Mallat算法是將信號(hào)f(t)分解到一系列子帶的濾波過程。各子帶的頻率范圍與信號(hào)f(t)的采樣頻率有關(guān)，設(shè)原始信號(hào)的采樣頻率為f(t)，Mallat算法中各子帶的頻率范圍列于表1。

表1 Mallat算法的頻率范圍Table 1 Limit’s frequency of Mallat algorithm

步驟2中，在得到FFT的結(jié)果后，首先應(yīng)構(gòu)建連續(xù)頻率軸，然后再通過隔點(diǎn)插零的形式進(jìn)行分解。其流程圖示于圖3。

圖3 消除取得近似部分頻率混疊的流程圖Fig.3 Flow chart of reducing frequency mixed of approximate part

2.2 消除取得細(xì)節(jié)部分Dj過程中頻率混疊的方法

細(xì)節(jié)部分是能量的高階頻率部分，其混疊會(huì)導(dǎo)致高頻信號(hào)混亂，影響信號(hào)細(xì)節(jié)頻率的輸出，其消除方法如下：

2) 將FFT結(jié)果中f<1+fs/2j部分譜值置零；

3) 對(duì)置零后的結(jié)果進(jìn)行快速傅里葉逆變換，結(jié)果作為DJ，不進(jìn)行隔點(diǎn)采樣。

消除流程示于圖4。

圖4 消除取得細(xì)節(jié)部分頻率混疊的流程圖Fig.4 Flow chart of reducing frequency mixed of detail part

2.3 消除近似部分Aj重構(gòu)過程中頻率混疊的方法

近似部分是能量的一個(gè)低階重要部分，其重構(gòu)過程中的混疊會(huì)導(dǎo)致低階信號(hào)混亂，影響信號(hào)的輸出，其消除方法如下：在第2j尺度上，由未被隔點(diǎn)采樣的信號(hào)開始重構(gòu)，這樣在2j尺度上不進(jìn)行隔點(diǎn)插零，而是從下1步，即第2j-1尺度上開始隔點(diǎn)插零，由于第1步重構(gòu)不產(chǎn)生虛假頻率成分，而后續(xù)的虛假成分一般可被h濾掉，因此就消除了Aj重構(gòu)過程中的頻率混疊[10]。

2.4 消除細(xì)節(jié)部分Aj重構(gòu)過程中頻率混疊的方法

細(xì)節(jié)部分是能量的一個(gè)高頻重要部分，其混疊會(huì)導(dǎo)致高頻信號(hào)混亂，影響信號(hào)細(xì)節(jié)頻率的輸出，其消除方法如下：

1) 設(shè)dj是由Dj重構(gòu)的結(jié)果，對(duì)dj進(jìn)行FFT；

2) 將FFT結(jié)果中非[1+fs/2j，fs/2j]部分譜值置零；

3) 對(duì)置零后的結(jié)果進(jìn)行快速傅里葉逆變換，其結(jié)果即為真實(shí)的dj。

3 改進(jìn)后的算法分析實(shí)例

在實(shí)現(xiàn)過程中，由理論上的推導(dǎo)轉(zhuǎn)化成實(shí)際的編程，還應(yīng)注意以下問題。

1) 信號(hào)在與濾波器函數(shù)卷積中的長(zhǎng)度問題。在將信號(hào)與濾波器函數(shù)進(jìn)行卷積時(shí)，使用的是Matlab中的convn函數(shù)，x=convn(a,b)計(jì)算的長(zhǎng)度為size(a)+size(b)-1。在分解與重構(gòu)過程中都有信號(hào)與濾波器卷積的環(huán)節(jié)，這樣最后得到的結(jié)果其長(zhǎng)度必然與原信號(hào)的長(zhǎng)度不同，因此需利用Matlab中的wkeep函數(shù)來控制信號(hào)的長(zhǎng)度。設(shè)N為原信號(hào)長(zhǎng)度，A1為第1層近似部分重構(gòu)后的結(jié)果，則運(yùn)行A1=wkeep(A1,N)命令即可取得A1的長(zhǎng)度為N的中心部分。這樣就保證了信號(hào)長(zhǎng)度的一致性。

2) 信號(hào)隔點(diǎn)插零和隔點(diǎn)采樣的方式。Matlab提供了dyaddown和dyadup兩個(gè)小波分析中用于隔點(diǎn)采樣和隔點(diǎn)插零的函數(shù)。dyaddown的用法為Y=dyaddown(X,EVENODD)，當(dāng)EVENODD取偶數(shù)時(shí)，Y(k)=X(2k)；當(dāng)EVENODD取奇數(shù)時(shí)，Y(k)=X(2k+1)。dyadup的用法為Y=dyadup(X,EVENODD)，當(dāng)EVENODD取偶數(shù)時(shí)，Y(2k-1)=X(k)，Y(2k)=0；當(dāng)EVENODD取奇數(shù)時(shí)，Y(2k-1)=0，Y(2k)=X(k)。為保證信號(hào)長(zhǎng)度在經(jīng)隔點(diǎn)采樣和隔點(diǎn)插零后不發(fā)生長(zhǎng)度上的變化，應(yīng)取適當(dāng)?shù)腅VENODD值以保持一致。在使用dyaddown時(shí)EVENODD取0，在使用dyadup時(shí)EVENODD取1。

為更清楚地體現(xiàn)改進(jìn)算法的效果，對(duì)某實(shí)驗(yàn)過程中的測(cè)量信號(hào)進(jìn)行分析。待分析信號(hào)為該實(shí)驗(yàn)采集到的微弱突變信號(hào)，采樣頻率為500 Hz，采樣點(diǎn)數(shù)取5 000，取dbN小波，N取40。利用默認(rèn)小波分析算法和改進(jìn)小波分析算法對(duì)信號(hào)進(jìn)行3層分解，結(jié)果分別示于圖5、6。由圖5可看出，信號(hào)D2和D1、D3和D2的頻譜存在明顯混疊。而從圖6可看出，各層的頻率都在理論上的頻段范圍內(nèi)，未出現(xiàn)頻率混疊現(xiàn)象。

4 結(jié)論

本文設(shè)計(jì)了Mallat小波變換改進(jìn)算法，利用傅里葉變換和傅里葉逆變換構(gòu)成Mallat小波變換在分解和重構(gòu)過程中所需的嚴(yán)格正交鏡像濾波器。改進(jìn)前后的小波算法對(duì)實(shí)際信號(hào)的分析結(jié)果表明，本文設(shè)計(jì)的Mallat小波變換改進(jìn)算法可消除信號(hào)分析中出現(xiàn)的頻率混疊現(xiàn)象。

圖5 默認(rèn)小波分析算法分析結(jié)果Fig.5 Analysis result of wavelet algorithm

圖6 改進(jìn)小波分析算法分析結(jié)果Fig.6 Analysis result of improvement wavelet algorithm