說明：研究LMD局域均值分解有3個(gè)月左右，能找到的相關(guān)文章也基本上看了一遍，覺得是個(gè)很好的方法，號稱是EMD經(jīng)驗(yàn)?zāi)B(tài)分解的改進(jìn)版。但是網(wǎng)絡(luò)上一直沒有找到該算法的matlab程序，只見文章說的天花亂墜。后來自己寫了一個(gè)，但是使用有一個(gè)問題沒有解決，就是分解的時(shí)候怎么去除騎行波的問題。先把自己寫的程序貢獻(xiàn)出來，讓大家分析下，一起討論下，到底LMD的程序怎樣寫才能如文獻(xiàn)中說的那樣達(dá)到目的。歡迎大家熱烈討論！

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

文件夾包含，找出純調(diào)頻函數(shù)，計(jì)算瞬時(shí)幅值，瞬時(shí)頻率的函數(shù)等

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

%lmd1原始版

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

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

%局域包絡(luò)函數(shù)=|上包絡(luò)-下包絡(luò)|/2

%相關(guān)文章見《一種基于EMD的振動(dòng)信號時(shí)頻分析新方法研究》-胡勁松，楊世錫

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

%最后一個(gè)PF是殘余分量

%A是瞬時(shí)賦值

%SI是純調(diào)頻函數(shù)，求它的瞬時(shí)頻率就是需要的頻率

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);

%停止調(diào)節(jié)

%1.emd用的是三次樣條求包絡(luò)，要求至少3個(gè)極值點(diǎn)，所以這里c的極值點(diǎn)個(gè)數(shù)也應(yīng)該至少為3

%2.如果上一個(gè)PF的極值點(diǎn)數(shù)比下一個(gè)PF的極值點(diǎn)數(shù)少，說明結(jié)果也不正確(這個(gè)也可以作為停止條件考慮進(jìn)去)

%上面一句是否可以等價(jià)于當(dāng)前PF的極值點(diǎn)個(gè)數(shù)一定要大于等于殘量(c)的極值點(diǎn)個(gè)數(shù)(目前是用這個(gè)作為停止條件的一個(gè)參考寫入程序)

%3.當(dāng)前PF分量的能量應(yīng)該大于殘量c的能量(這個(gè)有待商榷)

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

if

n_pf

if

length(c_pos)<3 ||

n_c

n_pfn_l

PF(k+1,:)=c;

break

end

end

end

%%%%%%%%%%%%%%%%%%nengliang函數(shù)%%%%%%%%%%

function p=nengliang(y)

% my=mean(y);

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

% p=sum(p);

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

end

%%%%%%%%%%%%%找純函數(shù)%%%%%%%%%%%%%%%%%

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;?% 邊界延拓點(diǎn)數(shù)

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);

% 判斷極值點(diǎn)個(gè)數(shù)

if (length(indmin) + length(indmax)

< 3)

error('not enough

extrema')

end

% 插值的邊界條件

if indmax(1) < indmin(1)%

第一個(gè)極值點(diǎn)是極大值

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

以第一個(gè)極大值為對稱中心

lmax =

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

lmin =

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

lsym =

indmax(1);

else%

如果第一個(gè)采樣值小于第一個(gè)極小值，則將認(rèn)為該值是一個(gè)極小值，以該點(diǎn)為對稱中心

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))% 以第一個(gè)極小值為對稱中心

lmax =

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

lmin =

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

lsym =

indmin(1);

else%

如果第一個(gè)采樣值大于第一個(gè)極大值，則將認(rèn)為該值是一個(gè)極大值，以該點(diǎn)為對稱中心

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

%

將序列根據(jù)對稱中心，鏡像到兩邊

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%

如果對稱的部分沒有足夠的極值點(diǎn)

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%

這種情況不應(yīng)該出現(xiàn)，程序直接中止

error('bug')

end

lsym = 1;?%

直接關(guān)于第一采樣點(diǎn)取鏡像

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

% 延拓點(diǎn)上的取值

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

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

% 極值點(diǎn)和過零點(diǎn)位置提取

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

% 當(dāng)連續(xù)多個(gè)采樣值相同時(shí)，把最中間的一個(gè)值作為極值點(diǎn)，處理方式與連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)

%一個(gè)輸出

%功能：找序列極值點(diǎn)位置坐標(biāo)

%y:輸入序列

%pos:極值點(diǎn)的序列位置坐標(biāo)

[indmin,indmax]=position(y);

minmax=cat(2,indmin,indmax);

poss=sort(minmax);

end

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

%返還極大值和極小值位置下標(biāo)

%兩個(gè)輸出

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

%說明每個(gè)程序要單獨(dú)保存成m文件，放在同一個(gè)文件夾下調(diào)用

使用實(shí)例：

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('殘量信號');

從上圖可以看出，成功將第一個(gè)分量和第二個(gè)分量分離出來，但是殘余分量存在很大的問題，這是因?yàn)榉纸膺^程中的騎行波沒有去處導(dǎo)致的，至于怎樣完善LMD局域均值分解算法，目前個(gè)人沒有時(shí)間研究了，希望得到指點(diǎn)。

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。