說明：研究LMD局域均值分解有3個月左右，能找到的相關文章也基本上看了一遍，覺得是個很好的方法，號稱是EMD經驗模態分解的改進版。但是網絡上一直沒有找到該算法的matlab程序，只見文章說的天花亂墜。后來自己寫了一個，但是使用有一個問題沒有解決，就是分解的時候怎么去除騎行波的問題。先把自己寫的程序貢獻出來，讓大家分析下，一起討論下，到底LMD的程序怎樣寫才能如文獻中說的那樣達到目的。歡迎大家熱烈討論！

程序仍存在不收斂的問題，拿出來分享只是希望高手指點一二，寫的不好歡迎拍磚！

文件夾包含，找出純調頻函數，計算瞬時幅值，瞬時頻率的函數等

%%%%%%%%%%%主程序%%%%%%%%%%

%lmd1原始版

%通過emd.m的三次樣條+鏡像延拓求出上下包絡及均值

%局域均值函數=(上包絡+下包絡)/2

%局域包絡函數=|上包絡-下包絡|/2

%相關文章見《一種基于EMD的振動信號時頻分析新方法研究》-胡勁松，楊世錫

function[PF,A,SI]=lmd1(m)

%最后一個PF是殘余分量

%A是瞬時賦值

%SI是純調頻函數，求它的瞬時頻率就是需要的頻率

c=m;

k=0;

wucha1=0.001;

n_l=nengliang(m);

while 1

k=k+1;

a=1;

h=c;

[pf,a,si]=zhaochun1(a,h,wucha1);

c=c-pf;

PF(k,:)=pf;

A(k,:)=a;

SI(k,:)=si;

c_pos=pos(c);

n_c=nengliang(c);

n_pf=nengliang(pf);

%停止調節

%1.emd用的是三次樣條求包絡，要求至少3個極值點，所以這里c的極值點個數也應該至少為3

%2.如果上一個PF的極值點數比下一個PF的極值點數少，說明結果也不正確(這個也可以作為停止條件考慮進去)

%上面一句是否可以等價于當前PF的極值點個數一定要大于等于殘量(c)的極值點個數(目前是用這個作為停止條件的一個參考寫入程序)

%3.當前PF分量的能量應該大于殘量c的能量(這個有待商榷)

%4.殘余能量不能大于信號能量

if

n_pf

if

length(c_pos)<3 ||

n_c

n_pfn_l

PF(k+1,:)=c;

break

end

end

end

%%%%%%%%%%%%%%%%%%nengliang函數%%%%%%%%%%

function p=nengliang(y)

% my=mean(y);

% p=(y-my).*(y-my);

% p=sum(p);

p=sum(abs(y).^2);

end

%%%%%%%%%%%%%找純函數%%%%%%%%%%%%%%%%%

function [pf,a,si]=zhaochun1(a,h,wucha1)

chun_num=0;

while 1

chun_num=chun_num+1;

t=1:length(h);

[envmin,envmax,envmoy,indmin,indmax,indzer] =

envelope(t,h,'spline');

mi=(envmax+envmin)./2;

ai=abs(envmax-envmin)./2;

a=a.*ai;

si=(h-mi)./ai;

h=si;

ai_funmax=max(ai);

ai_funmin=min(ai);

if

(ai_funmax<=1+wucha1&&ai_funmin>=1-wucha1)

break

end

end

pf=a.*si;

chun_num

function [envmin, envmax,envmoy,indmin,indmax,indzer] =

envelope(t,x,INTERP)

%computes envelopes and mean with various interpolations

NBSYM = 2;?% 邊界延拓點數

DEF_INTERP = 'spline';

if nargin < 2

x = t;

t = 1:length(x);

INTERP = DEF_INTERP;

end

if nargin == 2

if ischar(x)

INTERP = x;

x = t;

t = 1:length(x);

end

end

if ~ischar(INTERP)

error('interp parameter must be ''linear'''',

''cubic'' or ''spline''')

end

if ~any(strcmpi(INTERP,{'linear','cubic','spline'}))

error('interp parameter must be ''linear'''',

''cubic'' or ''spline''')

end

if min([size(x),size(t)]) > 1

error('x and t must be vectors')

end

s = size(x);

if s(1) > 1

x = x';

end

s = size(t);

if s(1) > 1

t = t';

end

if length(t) ~= length(x)

error('x and t must have the same length')

end

lx = length(x);

[indmin,indmax,indzer] = extr(x,t);

%boundary conditions for interpolation

[tmin,tmax,xmin,xmax] =

boundary_conditions(indmin,indmax,t,x,NBSYM);

% definition of envelopes from interpolation

envmax = interp1(tmax,xmax,t,INTERP);?envmin = interp1(tmin,xmin,t,INTERP);

if nargout > 2

envmoy =

(envmax + envmin)/2;

end

end

function [tmin,tmax,xmin,xmax] =

boundary_conditions(indmin,indmax,t,x,nbsym)

% computes the boundary conditions for interpolation (mainly mirror

symmetry)

lx = length(x);

% 判斷極值點個數

if (length(indmin) + length(indmax)

< 3)

error('not enough

extrema')

end

% 插值的邊界條件

if indmax(1) < indmin(1)%

第一個極值點是極大值

if x(1) > x(indmin(1))%

以第一個極大值為對稱中心

lmax =

fliplr(indmax(2:min(end,nbsym+1)));

lmin =

fliplr(indmin(1:min(end,nbsym)));

lsym =

indmax(1);

else%

如果第一個采樣值小于第一個極小值，則將認為該值是一個極小值，以該點為對稱中心

lmax =

fliplr(indmax(1:min(end,nbsym)));

lmin =

[fliplr(indmin(1:min(end,nbsym-1))),1];

lsym =

1;

end

else

if x(1) <

x(indmax(1))% 以第一個極小值為對稱中心

lmax =

fliplr(indmax(1:min(end,nbsym)));

lmin =

fliplr(indmin(2:min(end,nbsym+1)));

lsym =

indmin(1);

else%

如果第一個采樣值大于第一個極大值，則將認為該值是一個極大值，以該點為對稱中心

lmax =

[fliplr(indmax(1:min(end,nbsym-1))),1];

lmin =

fliplr(indmin(1:min(end,nbsym)));

lsym =

1;

end

end

%

序列末尾情況與序列開頭類似

if indmax(end) < indmin(end)

if x(end) <

x(indmax(end))

rmax =

fliplr(indmax(max(end-nbsym+1,1):end));

rmin =

fliplr(indmin(max(end-nbsym,1):end-1));

rsym =

indmin(end);

else

rmax =

[lx,fliplr(indmax(max(end-nbsym+2,1):end))];

rmin =

fliplr(indmin(max(end-nbsym+1,1):end));

rsym =

lx;

end

else

if x(end) >

x(indmin(end))

rmax =

fliplr(indmax(max(end-nbsym,1):end-1));

rmin =

fliplr(indmin(max(end-nbsym+1,1):end));

rsym =

indmax(end);

else

rmax =

fliplr(indmax(max(end-nbsym+1,1):end));

rmin =

[lx,fliplr(indmin(max(end-nbsym+2,1):end))];

rsym =

lx;

end

end

%

將序列根據對稱中心，鏡像到兩邊

tlmin = 2*t(lsym)-t(lmin);

tlmax = 2*t(lsym)-t(lmax);

trmin = 2*t(rsym)-t(rmin);

trmax = 2*t(rsym)-t(rmax);

% in case symmetrized parts do not extend enough%

如果對稱的部分沒有足夠的極值點

if tlmin(1) > t(1) | tlmax(1)

> t(1)% 對折后的序列沒有超出原序列的范圍

if lsym == indmax(1)

lmax =

fliplr(indmax(1:min(end,nbsym)));

else

lmin =

fliplr(indmin(1:min(end,nbsym)));

end

if lsym == 1%

這種情況不應該出現，程序直接中止

error('bug')

end

lsym = 1;?%

直接關于第一采樣點取鏡像

tlmin =

2*t(lsym)-t(lmin);

tlmax =

2*t(lsym)-t(lmax);

end?%

序列末尾情況與序列開頭類似

if trmin(end) < t(lx) | trmax(end)

< t(lx)

if rsym == indmax(end)

rmax =

fliplr(indmax(max(end-nbsym+1,1):end));

else

rmin =

fliplr(indmin(max(end-nbsym+1,1):end));

end

if rsym == lx

error('bug')

end

rsym = lx;

trmin =

2*t(rsym)-t(rmin);

trmax =

2*t(rsym)-t(rmax);

end

% 延拓點上的取值

xlmax =x(lmax);

xlmin =x(lmin);

xrmax =x(rmax);

xrmin =x(rmin);

% 完成延拓

tmin = [tlmin t(indmin) trmin];

tmax = [tlmax t(indmax) trmax];

xmin = [xlmin x(indmin) xrmin];

xmax = [xlmax x(indmax) xrmax];

end

%---------------------------------------------------------------------------------------------------

% 極值點和過零點位置提取

function [indmin, indmax, indzer] = extr(x,t);

%extracts the indices corresponding to extrema

if(nargin==1)

t=1:length(x);

end

m = length(x);

if nargout > 2

x1=x(1:m-1);

x2=x(2:m);

indzer = find(x1.*x2<0);

if any(x == 0)

iz = find( x==0 );

indz = [];

if any(diff(iz)==1)

zer = x == 0;

dz = diff([0 zer 0]);

debz = find(dz == 1);

finz = find(dz == -1)-1;

indz = round((debz+finz)/2);

else

indz = iz;

end

indzer = sort([indzer

indz]);

end

end

d = diff(x);

n = length(d);

d1 = d(1:n-1);

d2 = d(2:n);

indmin = find(d1.*d2<0 &

d1<0)+1;

indmax = find(d1.*d2<0 &

d1>0)+1;

% when two or more consecutive points have the same value we

consider only one extremum in the middle of the constant area

% 當連續多個采樣值相同時，把最中間的一個值作為極值點，處理方式與連0類似

if any(d==0)

imax = [];

imin = [];

bad = (d==0);

dd = diff([0 bad 0]);

debs = find(dd == 1);

fins = find(dd == -1);

if debs(1) == 1

if

length(debs) > 1

debs = debs(2:end);

fins = fins(2:end);

else

debs = [];

fins = [];

end

end

if length(debs) > 0

if fins(end)

== m

if length(debs) > 1

debs = debs(1:(end-1));

fins = fins(1:(end-1));

else

debs = [];

fins = [];

end?end

end

lc = length(debs);

if lc > 0

for k =

1:lc

if d(debs(k)-1) > 0

if d(fins(k)) < 0

imax = [imax round((fins(k)+debs(k))/2)];

end

else

if d(fins(k)) > 0

imin = [imin round((fins(k)+debs(k))/2)];

end

end

end

end

if length(imax) > 0

indmax =

sort([indmax imax]);

end

if length(imin) > 0

indmin =

sort([indmin imin]);

end

end?end

end

%%%%%%%%%%%%%%%%%%pos%%%%%%%%%%%%

function poss=pos(y)

%一個輸出

%功能：找序列極值點位置坐標

%y:輸入序列

%pos:極值點的序列位置坐標

[indmin,indmax]=position(y);

minmax=cat(2,indmin,indmax);

poss=sort(minmax);

end

%%%%%%%%%%%%%%position%%%%%%%

%返還極大值和極小值位置下標

%兩個輸出

function [indmin,indmax]=position(y)

m = length(y);

d = diff(y);

n = length(d);

d1 = d(1:n-1);

d2 = d(2:n);

indmin = find(d1.*d2<0 &

d1<0)+1;

indmax = find(d1.*d2<0 &

d1>0)+1;

if any(d==0)

imax = [];

imin = [];

bad = (d==0);

dd = diff([0 bad 0]);

debs = find(dd == 1);

fins = find(dd == -1);

if debs(1) == 1

if

length(debs) > 1

debs = debs(2:end);

fins = fins(2:end);

else

debs = [];

fins = [];

end

end

if length(debs) > 0

if fins(end)

== m

if length(debs) > 1

debs = debs(1:(end-1));

fins = fins(1:(end-1));

else

debs = [];

fins = [];

end?end

end

lc = length(debs);

if lc > 0

for k =

1:lc

if d(debs(k)-1) > 0

if d(fins(k)) < 0

imax = [imax round((fins(k)+debs(k))/2)];

end

else

if d(fins(k)) > 0

imin = [imin round((fins(k)+debs(k))/2)];

end

end

end

end

if length(imax) > 0

indmax =

sort([indmax imax]);

end

if length(imin) > 0

indmin =

sort([indmin imin]);

end

end

end

%說明每個程序要單獨保存成m文件，放在同一個文件夾下調用

使用實例：

x=@(t)

(2+cos(90*t).*cos(500*t+1800.*t.*t));

fs=5000;

t=0:1/fs:0.341;

y=x(t);

subplot(5,1,1);plot(t,y);xlabel('原始信號');

[pf,a,si]=lmd1(y);

subplot(5,1,2);plot(t,pf(1,:));xlabel('PF1');

subplot(5,1,3);plot(t,pf(2,:));xlabel('PF2');

subplot(5,1,4);plot(t,pf(3,:));xlabel('PF3');

subplot(5,1,5);plot(t,pf(4,:));xlabel('殘量信號');

從上圖可以看出，成功將第一個分量和第二個分量分離出來，但是殘余分量存在很大的問題，這是因為分解過程中的騎行波沒有去處導致的，至于怎樣完善LMD局域均值分解算法，目前個人沒有時間研究了，希望得到指點。

總結

以上是生活随笔為你收集整理的线谱法 时钟分量的提取 matlab,LMD局域均值分解的matlab程序及示例的全部內容，希望文章能夠幫你解決所遇到的問題。

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