lmd介紹.doc

1 介紹 由于在現(xiàn)實中存在著大量的非平穩(wěn)信號，因此針對非平穩(wěn)信號的處理與分析一直是信號處理領域里的一個研究熱點，在目前的非平穩(wěn)信號處理方法中Wigner-Villy分布、小波變換和經(jīng)驗模式分解法（EMD）等典型方法已經(jīng)得到了很廣泛的應用，不過也存在一定的局限性，其中Wigner-Villy分布由于是二次型時頻表示，對于多分量信號，存在交叉項干擾。小波變換對時頻面是一種機械的格型分解，所以無自適應性。EMD是一種自適應的信號處理方法，它將復雜的多分量信號自適應地分解為若干個IMF（Intrinsic mode function，簡稱IMF）分量之和，進一步對每個IMF分量進行Hilbert變換求出瞬時頻率和瞬時幅值，從而得到原始信號完整的時頻分布，由于EMD的自適應分解性等優(yōu)點使他在機械故障診斷、地震信號分析、海洋信號分析等領域已經(jīng)得到了廣泛的應用，但是EMD方法在理論上還存在一些問題，如過包絡、欠包絡、模態(tài)混淆和端點效應問題，這些問題仍在研究中。 最近Jonathan S. Smith在前人的研究基礎上提出了一種新的自適應非平穩(wěn)信號的處理方法局部均值分解（LMD）。并將這種方法應用于腦電信號分析，取得了不錯的效果。LMD自適應的將任何一個復雜的非平穩(wěn)信號分解成若干個瞬時頻率具有物理意義的PF分量之和，其中每一個PF分量由一個包絡信號和一個純調頻信號相乘而得到，包絡信號就是該PF分量的瞬時幅值，而PF分量的瞬時頻率則可由純調頻信號直接求出，進一步將所有PF分量的瞬時頻率和瞬時幅值組合，便可以得到原始信號完整的時頻分布。2 局部均值分解方法 LMD方法的實質是從原始信號中分離出純調頻信號和包絡信號，將純調頻信號和包絡信號相乘便可以得到一個瞬時頻率具有物理意義的PF分量，迭代處理至所有的PF分量分離出來，便可以得到原始信號的時頻分布。對于任意信號()xt，其分解過程如下。 1） 找出原始信號所有的局部極值點，求出所有相鄰的局部極值點的平均值： 將所有相鄰的平均值點用直線連接起來，然后用滑動平均法進行平滑處理，得到局部均值函數(shù)。 2）求出包絡估計值 將所有相鄰兩個包絡估計值用直線連接，然后采用滑動平均方法進行平滑處理，得到包絡估計函數(shù)。 3）將局部均值函數(shù)從原始信號中分離出來，得到 4） 用除以包絡估計函數(shù)以對進行解調 對重復上述步驟便能得到的包絡估計函數(shù)，假如不等于1，說明不是一個純調頻信號，需要重復上述迭代過程n次，直至為一個純調頻信號，也即的包絡估計函數(shù)，所以，有 式中， 迭代終止的條件為 實際應用中，在不影響分解效果的前提下，為了減少迭代次數(shù)，降低運算時間，可以用 作為迭代終止的條件。 5） 把迭代過程中產生的所有包絡估計函數(shù)相乘便可以得到包絡信號（瞬時幅值函數(shù)） 6） 將包絡信號和純調頻信號相乘便可以得到原始信號的第一個PF分量 它包含了原始信號中最高的頻率成分，是一個單分量的調幅-調頻信號，其瞬時幅值就是包絡信號，其瞬時頻率則可由純調頻信號求出，即 7）將第一個PF分量從原始信號中分離出來，得到一個新的信號，將作為原始數(shù)據(jù)重復以上步驟，循環(huán)k次，直到為一個單調函數(shù)為止。 原始信號能夠被所有的PF分量和重構，即 說明LMD分解沒有造成原信號信息的丟失。 3. 端點效應 LMD分解會受到端點效應的影響，由于局域均值包絡函數(shù)在端點處均存在一段未知的信號，若對端點不進行處理，在程序運行時，會自動給這部分信號添加一些虛假信息，從而對LMD分解產生影響。LMD端點效應首先發(fā)生在端點附近，然后在迭代過程中不斷向內部擴散，迭代次數(shù)越多端點效應污染整個數(shù)據(jù)段的程度就越嚴重。端點效應會使分解得到的各分量在端點附近產生一些變形，從而使結果不容易滿足循環(huán)終止條件，增加了循環(huán)次數(shù)，嚴重的時候會使數(shù)據(jù)產生嚴重失真。為減小端點效應對算法的影響，在分解前要對端點經(jīng)行一定的處理，應用最多的處理方法是鏡像延拓算法，鏡像延拓是在端點以外延拓一段信號。實際處理的信號兩端點一般不是極值點，這時候可以采用鏡像延拓的方法進行拓展，在LMD算法中只延拓一個極值點就可以很好的消除端點效應的影響。延拓方法如下圖所示，以離端點最近的一個極值點為對稱軸（如圖中虛線所示），將離端點次近的極值點向外延拓。在求解局部均值函數(shù)和局部包絡函數(shù)的時候將這個延拓的極值點代入，即可求得完整的局部均值函數(shù)和局部包絡函數(shù)。極值點延拓示意圖實驗用的數(shù)據(jù)fs=1000;N=1000;ts=1/fs;n=0:N-1;t=n/fs;x=(t) (1+0.5*cos(9*pi*t).*cos(200*pi*t+2*cos(10*pi*t)+sin(pi*t).*sin(30*pi*t);以下是圖和解釋這個是沒有經(jīng)過鏡像延拓法后的結果 依次是pf1 pf2 pf3 可以發(fā)現(xiàn)端點附近的值失真嚴重對應的顯示程序subplot(311);plot(t,pf(1,:);subplot(312);plot(t,pf(2,:);subplot(313);plot(t,pf(3,:);這是對于pf1 pf2的純調頻函數(shù) 依次是si1 si2對于的顯示程序subplot(211);plot(t,si(1,:);subplot(212);plot(t,si(2,:);這個是對應pf1 pf2的包絡信號 也就是瞬時幅值函數(shù) 依次是a1 a2對應的顯示程序subplot(211);plot(t,a(1,:);subplot(212);plot(t,a(2,:);這個是瞬時頻率乘以采樣率 就是ff1 ff2程序用的直接法求取瞬時頻率對應的顯示程序subplot(211);plot(t,ff(1,:);subplot(212);plot(t,ff(2,:);這個是經(jīng)過鏡像延拓法后的pf1 pf2 pf3 端點效應明顯減弱對應的顯示程序subplot(311);plot(t,pf(1,:);subplot(312);plot(t,pf(2,:);subplot(313);plot(t,pf(3,:);這個是經(jīng)過鏡像延拓法后的 這是對于pf1 pf2的純調頻函數(shù) 依次是si1 si2對于的顯示程序subplot(211);plot(t,si(1,:);subplot(212);plot(t,si(2,:);上面一個是包絡信號 也就是瞬時幅值函數(shù)下面一個是對包絡信號進行線性調頻Z變換（CZT）的頻譜圖function hua_xihua(x,fs,f1,fn,jingdu,style)對應的顯示程序subplot(211);plot(t,a(1,:);subplot(212);hua_xihua(a1,fs,0,100,256,2);上面一個是瞬時頻率乘以采樣率 下面一個是對瞬時頻率進行線性調頻Z變換（CZT）的頻譜圖function hua_fft1(y,fs)對應的顯示程序subplot(211);plot(t,ff(1,:);subplot(212);hua_xihua(f1,fs,0,100,256,2