【數字信号處理】基于matlab LMD算法和ELMD算法管道洩漏信号處理【含Matlab源碼 1985期】

一、局部分解理論研究

局部均值分解算法(Local Mean Decomposition， LMD) 作為處理非平穩随機信号的一種手段，得到了廣泛應用，并成熟地應用于機械故障診斷、信号特征提取與分析等方面。LMD算法最大的特點就在于其對信号的自适應分解能力， 這種自适應能力主要展現在該方法能夠通過資料自身特點，通過特定手段，将原始信号分為不同模态函數， 針對不同模态函數進一步處理。與此同時， 局部均值分解算法(LMD) 相較于模态分解的創始算法經驗模态分解算法(Empirical Mode Decomposition， EMD) 而言，其具備端點效應小、疊代次數少等優勢。本章将詳述LMD算法基本原理及分解流程，并針對LMD存在的基本問題， 進行改進研究。

1 局部均值分解算法理論研究

2005年， JonathanS.Smith提出的一種新的非線性和非平穩信号分析方法一局部均值分解算法(LMD) ， 并應用于腦電信号的分析中， 取得不錯的效果。此外， 局部均值分解在機械故障診斷中也得到良好的應用。局部均值分解可以依據信号本身的特征進行自适應分解的，産生具有真實實體意義的PF分量，并由此得到能夠清晰準确反映出信号能量在空間各尺度上分布規律的時頻分布，有利于更加細緻的對信号特征進行分析。

1.1局部均值算法(LMD) 分解流程

利用LMD分解， 可以将原始信号分解并産生若幹有效PF分量， 将所有PF分量相加即可重構出原始信号。其中每個PF分量都是一個純調頻信号和包絡信号的乘積，且每個PF分量的瞬時頻率具有實際實體意義。LMD分解的具體流程如下：(1)原始信号x(t)，找出x(t)上的所有局部極值點n，由相鄰的兩個極值點n，n.計算出的一個均值m，即

随後将得到的所有平均值m，用折線連接配接，連接配接過程采用滑動平均方法進行平滑處理，進而得到局部均值函數m，(t)。同時，利用相鄰極值點計算包絡估計值a，，即

通過公式(2-12)可以發現，信号重構過程中，不會出現資訊的丢失，保證了信号的完整性。具體流程圖如下圖1.1所示

(1)總體局部均值分解算法研究

模态混疊現象的發生一方面來自于算法本身，另一方面則會受到原始信号頻率特征的影響。當待處理信号确定之後，信号特征頻率變不會發生改變，是以，為抑制模态混疊現象的出現，從算法本身出發，引入噪聲輔助信号處理方法，即在信号中加入白噪聲來平滑脈沖幹擾。因為在LMD分解過程中， 需要對極值點進行處理， 得到局部均值函數及包絡估計函數，極值點的分布就會影響到包絡的拟合情況。如果信号的極值點分布不均勻， 則極易産生模态混疊。是以， 在LMD中， 借助高斯白噪聲輔助法， 對模态混疊進行抑制， 得到總體局部均值分解(Ensemble Local Mean Decomposition ELMD) 算法。ELMD分解， 是指在LMD分解前， 将不同有限幅值的白噪聲信号加入待分解信号，利用白噪聲均值為零，頻譜能量分布均勻的特性，使得白噪聲可以均勻的分布在整個時頻空間中，并且不同時間尺度的信号會自動分布到與背景噪聲相關的适當尺度上去。對于單次試驗，由于噪聲的添加，使得每次結果都會受到噪聲的影響産生偏差，這是由于在分解過程中，信号包括原始信号及附加的白噪聲。但随着試驗的次數增加，由于零均值的特性，噪聲将會互相抵消，進而得到消除，唯一持久穩固的部分便是信号本身， 是以可以認定內建均值的結果就是最終分解結果， 即EL MD算法。

總體局部均值EL MD分解算法流程圖如圖2.13所示， 即對原始信号分别加入n組不同高斯白噪聲， 分别進行LMD分解， 随後将得到的n組PF分量進行平均處理，得到最終的LMD分解結果。

二、部分源代碼

clear;clc;close all;
x1=xlsread('5.xlsx');
VarName2=x1(:,2);
x=VarName2';
Nstd =0.1;
NR =100;
%[modes its]=eemd(ecg,0.1,100,1000);
modes = mlmd(x,Nstd,NR);
[a, b]=size(modes);
N=length(x);
fs=1000;
Ts=1/fs;
t=0:Ts:N*Ts;
t=t(1:N);
% 繪制仿真信号和其頻譜圖
figure(1)
subplot(211)
plot(t,x)
subplot(212)
y2=x;
L=length(y2);
NFFT = 2^nextpow2(L);
Y = fft(y2,NFFT)/L;
f = fs/2*linspace(0,1,NFFT/2);
plot(f,2*abs(Y(1:NFFT/2)))
PF=modes;
line=size(PF,1); 
NN = length(PF(1,:));
n = linspace(0,1,NN);
for k1=0:4:line-1
     figure('Color',[1 1 1]);
     for k2=1:min(4,line-k1)
        subplot(4,2,2*k2-1);
        plot(t,PF(k1+k2,:));
        title(sprintf('第%d個PF', k1+k2))
        xlabel('Time/s')
        ylabel(sprintf('PF%d',k1+k2));
        subplot(4,2,2*k2)
        [yf, f] = FFTAnalysis(PF(k1+k2,:), Ts);
        plot(f, yf)
        title(sprintf('第%d個PF的頻譜', k1+k2))
        xlabel('f/Hz')
        ylabel('|PF(f)|');
     end
end;
% 頻譜分析
function [Y, f] = FFTAnalysis(y, Ts)
Fs = 1/Ts;
L = length(y);
NFFT = 2^nextpow2(L);
% y = y - mean(y);
Y = fft(y, NFFT)/L;
Y = 2*abs(Y(1:NFFT/2+1));
f = Fs/2*linspace(0, 1, NFFT/2+1);
end
function [pf,a,si,u] = lmd(x)
x=x';
c = x';
N = length(x);
A = ones(1,N);
PF = [];
AA=[];
SI=[];
U=[];
aii = 2*A;

while(1)

  si = c;
  a = 1;
  
   while(1)
    h = si;
    
      maxVec = [];
      minVec = [];
      
   % look for max and min point
      for i = 2: N - 1
         if h (i - 1) < h (i) && h (i) > h (i + 1)
            maxVec = [maxVec i];    
         end
         if h (i - 1) > h (i) && h (i) < h (i + 1)
            minVec = [minVec i];    
         end         
      end
      
   % check if it is residual
      if (length (maxVec) + length (minVec)) < 2
         break;
      end
           
  % handle end point 
      lenmax=length(maxVec);
      lenmin=length(minVec);
      %left end point
      if h(1)>0
          if(maxVec(1)<minVec(1))
              yleft_max=h(maxVec(1));
              yleft_min=-h(1);
          else
              yleft_max=h(1);
              yleft_min=h(minVec(1));
          end
      else
          if (maxVec(1)<minVec(1))
              yleft_max=h(maxVec(1));
              yleft_min=h(1);
          else
              yleft_max=-h(1);
              yleft_min=h(minVec(1));
          end
      end
      %right end point
      if h(N)>0
          if(maxVec(lenmax)<minVec(lenmin))
             yright_max=h(N);
             yright_min=h(minVec(lenmin));
          else
              yright_max=h(maxVec(lenmax));
              yright_min=-h(N);
          end
      else
          if(maxVec(lenmax)<minVec(lenmin))
              yright_max=-h(N);
              yright_min=h(minVec(lenmin));
          else
              yright_max=h(maxVec(lenmax));
              yright_min=h(N);
          end
      end
      %get envelop of maxVec and minVec using
      %spline interpolate
      maxEnv=spline([1 maxVec N],[yleft_max h(maxVec) yright_max],1:N);
      minEnv=spline([1 minVec N],[yleft_min h(minVec) yright_min],1:N);
      
    mm = (maxEnv + minEnv)/2;
    aa = abs(maxEnv - minEnv)/2;
    
    mmm = mm;
    aaa = aa;
    preh = h;
    h = h-mmm;
    si = h./aaa;
    a = a.*aaa;    
    aii = aaa;
    B = length(aii);
    C = ones(1,B);
    bb = norm(aii-C);
    if(bb < 1000)
        break;      

三、運作結果

四、matlab版本及參考文獻

1 matlab版本

2014a或2019b