繼續完善“類EMD”方法系列，本篇是繼EEMD、CEEMD、CEEMDAN、VMD、ICEEMDAN后的第6篇，想要看前幾種方法的點擊鏈接可以跳轉。

LMD（local mean decomposition，局部均值分解）方法是2005年由Smith等人[1]提出的，本質上是根據信號的包絡特征自適應地將一個非線性、非平穩信號按頻率遞減的順序逐級分離。LMD的提出也是用來解決EMD分解的端點效應和模態混疊問題，最開始是用來處理腦電數據的。

1. LMD的概念

與EMD衍生的一系列方法不同，經過LMD分解得到的分量被稱作“乘積函數(PF)”，即每個PF都是通過包絡函數乘以純調頻函數得到的。其中包絡函數是PF的瞬時幅值，純調頻函數的頻率是PF的瞬時頻率。

曲線為待分解曲線

圖解一下分解步驟：

（1）找到信號中的全部極值點??。

?

（2）求出相鄰極值點的局部均值點?，也就是相鄰兩個極值點的中點。公式為：?

?

藍色標記為局部均值點

（3）求出相鄰極值點的局部包絡??，也就是相鄰兩個極值點的幅值絕對值差值的一半。公式為?

?

綠色標記為局部包絡值

（4）將局部均值的折線連線進行平滑處理，得到局部均值函數?，如下圖中的橙色曲線。

?

橙色曲線為局部均值函數（示意）

（5）將局部包絡同樣進行平滑處理，得到局部包絡估計函數??，如下圖紫色曲線。

?

局部包絡估計函數（示意）

（6）將從原始信號的序列里面分割出來，得到零均值信號。

?

?

（7）對??進行解調，通過??除以??可以得到??。反復重復以上（1）~（6）步，直到包絡估計函數??，此時得到的??就是純調頻信號。

?

（8）將上述迭代獲得的全部局部包絡估計函數進行乘法運算就可以得到包絡信號??。

?

（9）第一個PF分量就可以寫為：

?（10）將原始信號減掉PF1，得到的剩余信號重復（1）~（9）步，直到剩余信號為單調函數為止，此時原始信號就被分解為k個PF和一個單調剩余信號。即：

分解完成！

2.?LMD的特點

對于LMD和EMD的區別，這里就直接引用論文了

（1）PF 分量和 IMF 分量的含義不相同。經過 EMD 分解后獲得的 IMF 屬于調頻信號，而利用 LMD 分解后獲得的 PF 分量屬于調幅調頻信號。要想得到 IMF 分量，必須保證原始信號極值點的數量絕對等于過零點的數量，或者兩個數值的數值差的結果小于等于 1，因此 IMF 分量并不會顯現出不過零點的局部波動；但 PF并不需要滿足這個條件。綜合上所述可以發現 PF 分量能夠更準確的反映原始信號的所有特征信息。
（2）針對局部均值函數 EMD 和 LMD 的求解方法存在明顯的差異。EMD 是分別對所有極值點采用三次樣條插值獲得原始信號的上包絡線和下包絡線，接下來使用平均值的方法就可以得到局部均值函數，采用這種方法更容易形成過包絡或欠包絡等缺陷；對于LMD 求解局部均值函數，求取相鄰兩個極值的平均值，并利用滑動平均算法對其進行平滑處理；LMD 能夠避免 EMD 中存在的過包絡和欠包絡的缺點。因此通過對比可以發現LMD 的分解結果更加準確。
（3）LMD 和 EMD 對瞬時頻率的求解思路不同。在 EMD 中，必須求解 Hilbert 才能獲得 IMF 的瞬時頻率，然后再利用對其瞬時相位的求取倒數，最終獲得瞬時頻率，當若干個 IMF 中的一個瞬時相位具有突變的時候，求解出來的瞬時頻率可能會出現工程實際中難以解釋的負值；但對于 LMD 則不會出現頻率為負值的情況，因為瞬時頻率的值可以直接通過分解后的 PF 分量直接計算得到，這種方法不但簡單而且求出的瞬時頻率值都屬于正值。因此求解瞬時頻率的時候 LMD 方法更占優勢。
（4）LMD 和 EMD 的整個分解過程計算量有所不同。針對 EMD 的求解過程主要存在兩個迭代過程，一個是獲取所需要的若干個 IMF 分量，另一個則是將所有的 IMF 從原始信號中分離出來，最終得到一個殘余分量；而 LMD 方法就相對于 EMD 復雜一些，其主要包含了三個迭代過程，第一個迭代是利用滑動平均算法對局部均值和局部包絡的折線進行平滑處理，得到局部均值函數以及包絡估計函數，第二個迭代是通過迭代過程獲得一個純調頻函數，第三個迭代就是將所有的 PF 分量全部計算出來。針對計算量來說，LMD的并不占優勢，其計算量相比于 EMD 稍大一些。

3.?LMD的MATLAB編程實現

在互聯網上沒有找到十分靠譜的LMD程序，于是筆者決定動手寫一個出來。

但是寫出來的代碼運行卻很有問題。

模態混疊端點效應都算小問題了。

PF1分量中存在模態混疊

某些信號分量會出現巨大的局部畸形，且其數量級遠原始信號。

比如對上圖中的原始信號加入微量噪聲后，分解結果就會變成這樣：

注：res還沒達到單調，如果要單調需要分解到十幾個分量，這里強行停止并將最后幾個分量重構了

前三個PF算是勉強還能看的話，從PF4分量就開始暴走了。

所以目前的代碼狀態是不可用。但是處于交流的目的，代碼還是貼上來，如果有同學解決了這個問題記得告訴我一聲哦。

function pf = kLMD(data,tol,MaxNum) % LMD分解函數 % 輸入： % data：待分解信號 % tol：迭代停止閾值 % MaxNum：最大PF分量數，設置此數是為了防止分解難以達到停止條件，陷入過多次循環。默認為10 % 輸出： % pf：乘積函數，即分解得到的分量，每一行為一個分量% 鑒于互聯網上沒找到特別靠譜的代碼，按照方法流程重頭寫了這個版本。 % 參考論文（即方法提出論文）為： % Smith, Jonathan S . The local mean decomposition and its application to EEG perception data[J]. % Journal of The Royal Society Interface, 2005, 2(5):443-454. % 筆者對于該方法的介紹文章（類EMD的“信號分解方法”及MATLAB實現（第六篇）——LMD）為： % https://zhuanlan.zhihu.com/p/444277130 % % Copyright (c) 2021 Mr.看海 All rights reserved.if ~exist('MaxNum') %如果未設置最大pf數，則設置為10MaxNum = 10; end k = 0; while 1dataTemp{1} = data;k = k+1;[pf(k,:)] = getpf(dataTemp{k},tol,'off');if k == 1dataTemp{k+1} = dataTemp{k}(:)-(pf(:)); %將PF從信號中分割出來elsedataTemp{k+1} = dataTemp{k}(:)'-(pf(k,:)); %將PF從信號中分割出來endif sum(diff(dataTemp{k+1})>0)==0 || sum(diff(dataTemp{k+1})<0)==0 %如果殘值單調，結束分解pf(k+1,:) = dataTemp{k+1}; %殘值breakendif k>=MaxNum %達到最大迭代次數，停止分解pf(k+1,:) = dataTemp{k+1}; %殘值breakend endendfunction [pf] = getpf(data,tol,figflag) % 通過迭代過程獲得一個純調頻函數 i = 0; while 1i = i + 1;if i == 1[m11,a11{i},h11,s11{i}] = getamhs(data,figflag); %調用函數，求得m11，a11，s11等else [m11,a11{i},h11,s11{i}] = getamhs(s11{i-1},figflag); %調用函數，求得m11，a11，s11等endif sum((a11{i} - 1).^2) < tol %達到閾值，迭代停止breakend enda = 1; for j = 1:ia = a.*a11{j}; %包絡信號 end pf = a.*s11{i}; %得到此階的乘積函數 pf = pf'; %轉變為行向量 endfunction [mm11,aa11,h11,s11] = getamhs(data,figflag) data = data(:)'; len = length(data); %% 1.求極值點 [pksP,locsP] = findpeaks(data); %找到最大值 [pksN,locsN] = findpeaks(-data); pksN = -pksN; %找到最小值 %% 2.求局部均值點m ns = [locsP,locsN;pksP,pksN]; %極值點 [B,I] = sort(ns(1,:)); %獲取排序索引 nsSort = ns(:,I); %對極值點m按順序排列 nsSortT = [[1;data(1)],nsSort,[len;data(len)]]; %將端點作為一個極值點加入 for i = 1:length(nsSortT(1,:))-1ms(i) = (nsSortT(2,i) + nsSortT(2,i+1))/2; msLocs(i) = (nsSortT(1,i) + nsSortT(1,i+1))/2; ai(i) = abs(nsSortT(2,i) - nsSortT(2,i+1))/2;msLocs(i) = (nsSortT(1,i) + nsSortT(1,i+1))/2; end %% 3.求局部均值函數 m11 = movmean([msLocs',ms'],2); %滑動平均 try[fm,~] = fit(m11(:,1),m11(:,2),'smoothingspline'); %曲線擬合 catch MEm11 = movmean([msLocs',ms'],2); %滑動平均[fm,~] = fit(m11(:,1),m11(:,2),'smoothingspline'); %曲線擬合 end mm11 = fm(1:len); %% 4.局部包絡 a11 = movmean([msLocs',ai'],2); %滑動平均 try[fa,~] = fit(a11(:,1),a11(:,2),'smoothingspline'); %曲線擬合 catch MEa11 = movmean([msLocs',ai'],2); %滑動平均[fa,~] = fit(a11(:,1),a11(:,2),'smoothingspline'); %曲線擬合 end %[fa,~] = fit(msLocs',ai','smoothingspline'); %曲線擬合 aa11 = fa(1:len); %% 5.求零均值信號 h11 = data-mm11'; %零均值信號 %% 6.調頻信號 s11 = h11(:)./aa11(:); end

測試用的仿真信號為：

fs = 1e3; t = 0:1/fs:1-1/fs; x = 8*t.^2 + cos(4*pi*t+10*pi*t.^2) + ...[cos(60*pi*(t(t<=0.5))) cos(200*pi*(t(t>0.5)-10*pi))]; n = rand(1,1000); sig = x+0.1*n;

其他

EMD、EEMD、CEEMD、CEEMDAN、ICEEMDAN、VMD以及HHT相關的程序筆者封裝了畫圖函數。關于EMD、EEMD、CEEMD、VMD和HHT的相關介紹可以看這里：

Mr.看海：這篇文章能讓你明白經驗模態分解（EMD）——EMD在MATLAB中的實現方法

Mr.看海：希爾伯特譜、邊際譜、包絡譜、瞬時頻率/幅值/相位——Hilbert分析衍生方法及MATLAB實現

Mr.看海：類EMD的“信號分解方法”及MATLAB實現（第一篇）——EEMD

Mr.看海：類EMD的“信號分解方法”及MATLAB實現（第二篇）——CEEMD

Mr.看海：類EMD的“信號分解方法”及MATLAB實現（第三篇）——CEEMDAN

Mr.看海：類EMD的“信號分解方法”及MATLAB實現（第四篇）——VMD

Mr.看海：類EMD的“信號分解方法”及MATLAB實現（第五篇）——ICEEMDAN

參考

^[1]Smith, Jonathan S . The local mean decomposition and its application to EEG perception data[J]. Journal of The Royal Society Interface, 2005, 2(5):443-454.

總結

以上是生活随笔為你收集整理的类EMD的“信号分解方法”及MATLAB实现（第六篇）——LMD的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。