emd,eemd,vmd,频谱图，分解图对比matlab代码

作为自己学习的一个记录吧。

对这个信号进行实验，其中公式是截图截的，懒得改了，f1就是s1。

对这个s信号进行分解。下面开始代码操作：

原始信号生成：运行该段代码，生成一个s.mat数据，并作图。

clear
clc
close all
t = 0:0.001:2;
s1 = cos(4*pi.*t);
figure
plot(t,s1)
%%
s2 = 1/4*cos(48*pi.*t);
figure
plot(t,s2)
%%
s3 = 1/16*cos(576*pi.*t);
figure
plot(t,s3)
%%
s4 = 3*wgn(1,length(t),10*log10(0.01));
figure
plot(t,s4)
figure
s = s1+s2+s3+s4;
plot(t,s)save s.mat s

先上main_emd的主函数

这段代码就是对原始信号s进行emd分解，并出图。

clear
clc
close all
load s.mat
%% EEMD分解
[u its]=emd(s);
t = 0:0.001:2;
% [a b]=size(u);
a = 8;figure(1);
imfn=u;
subplot(a+1,1,1);
plot(t,s); %故障信号
ylabel('s','fontsize',12,'fontname','宋体');for n1=1:asubplot(a+1,1,n1+1);plot(t,u(n1,:));%输出IMF分量，a(:,n)则表示矩阵a的第n列元素，u(n1,:)表示矩阵u的n1行元素ylabel(['IMF' int2str(n1)]);%int2str(i)是将数值i四舍五入后转变成字符，y轴命名endxlabel('时间\itt/s','fontsize',12,'fontname','宋体');%%
figure('Name','频谱图','Color','white');K = a;
L = length(t);
fs = 2001;
for i = 1:Kp=abs(fft(u(i,:))); %并fft，得到p，就是包络线的fft---包络谱 subplot(K,1,i);plot((0:L-1)*fs/L,p)   %绘制包络谱xlim([0 fs/2]) %展示包络谱低频段，这句代码可以自己根据情况选择是否注释if i ==1title('频谱图'); ylabel(['IMF' int2str(i)]);%int2str(i)是将数值i四舍五入后转变成字符，y轴命名elseif i==Kxlabel('频率');  ylabel(['IMF' int2str(i)]);%int2str(i)是将数值i四舍五入后转变成字符，y轴命名elseylabel(['IMF' int2str(i)]);%int2str(i)是将数值i四舍五入后转变成字符，y轴命名end
end
set(gcf,'color','w');

再上main_eemd的主函数

这段代码就是对原始信号s进行eemd分解，并出图。其实和emd的主函数差不多。

clear
clc
% close all
load s.mat
%% EEMD分解
Nstd = 0.5;
NR = 500;
MaxIter = 5000;
[u its]=eemd(s,Nstd ,NR,MaxIter);
t = 0:0.001:2;
% [a b]=size(u);
a = 8;figure(1);
imfn=u;
subplot(a+1,1,1);
plot(t,s); %故障信号
ylabel('s','fontsize',12,'fontname','宋体');
for n1=1:asubplot(a+1,1,n1+1);plot(t,u(n1,:));%输出IMF分量，a(:,n)则表示矩阵a的第n列元素，u(n1,:)表示矩阵u的n1行元素ylabel(['IMF' int2str(n1)]);%int2str(i)是将数值i四舍五入后转变成字符，y轴命名
endxlabel('时间\itt/s','fontsize',12,'fontname','宋体');%%
figure('Name','频谱图','Color','white');K = a;
L = length(t);
fs = 2001;
for i = 1:Kp=abs(fft(u(i,:))); %并fft，得到p，就是包络线的fft---包络谱 subplot(K,1,i);plot((0:L-1)*fs/L,p)   %绘制包络谱xlim([0 fs/2]) %展示包络谱低频段，这句代码可以自己根据情况选择是否注释if i ==1title('频谱图'); ylabel(['IMF' int2str(i)]);%int2str(i)是将数值i四舍五入后转变成字符，y轴命名elseif i==Kxlabel('频率');  ylabel(['IMF' int2str(i)]);%int2str(i)是将数值i四舍五入后转变成字符，y轴命名elseylabel(['IMF' int2str(i)]);%int2str(i)是将数值i四舍五入后转变成字符，y轴命名end
end
set(gcf,'color','w');

然后是main_vmd的主函数

clear
clc
close all
load s.mat
X = s;
t = 0:0.0005:1;
%--------- some sample parameters forVMD：对于VMD样品参数进行设置---------------
alpha = 2000;       % moderate bandwidth constraint：适度的带宽约束/惩罚因子
tau = 0;          % noise-tolerance (no strict fidelity enforcement)：噪声容限（没有严格的保真度执行）
K = 8;             % modes：分解的模态数，可以自行设置，这里以8为例。
DC = 0;             % no DC part imposed：无直流部分
init = 1;           % initialize omegas uniformly  ：omegas的均匀初始化
tol = 1e-10;
%--------------- Run actual VMD code:数据进行vmd分解---------------------------
% [u,omega] = pVMD(s,fs, alpha, K, tol);
[u, u_hat, omega] = VMD(X, alpha, tau, K, DC, init, tol); %其中u为分解得到的IMF分量
a = K;
figure(1);
imfn=u;
subplot(a+1,1,1);
plot(t,s); %故障信号
ylabel('s','fontsize',12,'fontname','宋体');
for n1=1:asubplot(a+1,1,n1+1);plot(t,u(n1,:));%输出IMF分量，a(:,n)则表示矩阵a的第n列元素，u(n1,:)表示矩阵u的n1行元素ylabel(['IMF' int2str(n1)]);%int2str(i)是将数值i四舍五入后转变成字符，y轴命名
endxlabel('时间\itt/s','fontsize',12,'fontname','宋体');%%
figure('Name','频谱图','Color','white');K = a;
L = length(t);
fs = 2001;
for i = 1:Kp=abs(fft(u(i,:))); %并fft，得到p，就是包络线的fft---包络谱 subplot(K,1,i);plot((0:L-1)*fs/L,p)   %绘制包络谱xlim([0 fs/2]) %展示包络谱低频段，这句代码可以自己根据情况选择是否注释if i ==1title('频谱图'); ylabel(['IMF' int2str(i)]);%int2str(i)是将数值i四舍五入后转变成字符，y轴命名elseif i==Kxlabel('频率');  ylabel(['IMF' int2str(i)]);%int2str(i)是将数值i四舍五入后转变成字符，y轴命名elseylabel(['IMF' int2str(i)]);%int2str(i)是将数值i四舍五入后转变成字符，y轴命名end
end
set(gcf,'color','w');

最后再上emd、eemd、vmd的函数代码

首先是emd的：

function [imf,ort,nbits] = emd(varargin)[x,t,sd,sd2,tol,MODE_COMPLEX,ndirs,display_sifting,sdt,sd2t,r,imf,k,nbit,NbIt,MAXITERATIONS,FIXE,FIXE_H,MAXMODES,INTERP,mask] = init(varargin{:});if display_siftingfig_h = figure;
end%main loop : requires at least 3 extrema to proceed
while (~stop_EMD(r,MODE_COMPLEX,ndirs) && (k < MAXMODES+1 || MAXMODES == 0) && ~any(mask))% current modem = r;% mode at previous iterationmp = m;%computation of mean and stopping criterionif FIXE[stop_sift,moyenne] = stop_sifting_fixe(t,m,INTERP,MODE_COMPLEX,ndirs);elseif FIXE_Hstop_count = 0;[stop_sift,moyenne] = stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H,MODE_COMPLEX,ndirs);else[stop_sift,moyenne] = stop_sifting(m,t,sd,sd2,tol,INTERP,MODE_COMPLEX,ndirs);end% in case the current mode is so small that machine precision can cause% spurious extrema to appearif (max(abs(m))) < (1e-10)*(max(abs(x)))if ~stop_siftwarning('emd:warning','forced stop of EMD : too small amplitude')elsedisp('forced stop of EMD : too small amplitude')endbreakend% sifting loopwhile ~stop_sift && nbit<MAXITERATIONSif(~MODE_COMPLEX && nbit>MAXITERATIONS/5 && mod(nbit,floor(MAXITERATIONS/10))==0 && ~FIXE && nbit > 100)disp(['mode ',int2str(k),', iteration ',int2str(nbit)])if exist('s','var')disp(['stop parameter mean value : ',num2str(s)])end[im,iM] = extr(m);disp([int2str(sum(m(im) > 0)),' minima > 0; ',int2str(sum(m(iM) < 0)),' maxima < 0.'])end%siftingm = m - moyenne;%computation of mean and stopping criterionif FIXE[stop_sift,moyenne] = stop_sifting_fixe(t,m,INTERP,MODE_COMPLEX,ndirs);elseif FIXE_H[stop_sift,moyenne,stop_count] = stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H,MODE_COMPLEX,ndirs);else[stop_sift,moyenne,s] = stop_sifting(m,t,sd,sd2,tol,INTERP,MODE_COMPLEX,ndirs);end% displayif display_sifting && ~MODE_COMPLEXNBSYM = 2;[indmin,indmax] = extr(mp);[tmin,tmax,mmin,mmax] = boundary_conditions(indmin,indmax,t,mp,mp,NBSYM);envminp = interp1(tmin,mmin,t,INTERP);envmaxp = interp1(tmax,mmax,t,INTERP);envmoyp = (envminp+envmaxp)/2;if FIXE || FIXE_Hdisplay_emd_fixe(t,m,mp,r,envminp,envmaxp,envmoyp,nbit,k,display_sifting)elsesxp=2*(abs(envmoyp))./(abs(envmaxp-envminp));sp = mean(sxp);display_emd(t,m,mp,r,envminp,envmaxp,envmoyp,s,sp,sxp,sdt,sd2t,nbit,k,display_sifting,stop_sift)endendmp = m;nbit=nbit+1;NbIt=NbIt+1;if(nbit==(MAXITERATIONS-1) && ~FIXE && nbit > 100)if exist('s','var')warning('emd:warning',['forced stop of sifting : too many iterations... mode ',int2str(k),'. stop parameter mean value : ',num2str(s)])elsewarning('emd:warning',['forced stop of sifting : too many iterations... mode ',int2str(k),'.'])endendend % sifting loopimf(k,:) = m;if display_siftingdisp(['mode ',int2str(k),' stored'])endnbits(k) = nbit;k = k+1;r = r - m;nbit=0;end %main loopif any(r) && ~any(mask)imf(k,:) = r;
endort = io(x,imf);if display_siftingclose
end
end%---------------------------------------------------------------------------------------------------
% tests if there are enough (3) extrema to continue the decomposition
function stop = stop_EMD(r,MODE_COMPLEX,ndirs)
if MODE_COMPLEXfor k = 1:ndirsphi = (k-1)*pi/ndirs;[indmin,indmax] = extr(real(exp(i*phi)*r));ner(k) = length(indmin) + length(indmax);endstop = any(ner < 3);
else[indmin,indmax] = extr(r);ner = length(indmin) + length(indmax);stop = ner < 3;
end
end%---------------------------------------------------------------------------------------------------
% computes the mean of the envelopes and the mode amplitude estimate
function [envmoy,nem,nzm,amp] = mean_and_amplitude(m,t,INTERP,MODE_COMPLEX,ndirs)
NBSYM = 2;
if MODE_COMPLEXswitch MODE_COMPLEXcase 1for k = 1:ndirsphi = (k-1)*pi/ndirs;y = real(exp(-i*phi)*m);[indmin,indmax,indzer] = extr(y);nem(k) = length(indmin)+length(indmax);nzm(k) = length(indzer);[tmin,tmax,zmin,zmax] = boundary_conditions(indmin,indmax,t,y,m,NBSYM);envmin(k,:) = interp1(tmin,zmin,t,INTERP);envmax(k,:) = interp1(tmax,zmax,t,INTERP);endenvmoy = mean((envmin+envmax)/2,1);if nargout > 3amp = mean(abs(envmax-envmin),1)/2;endcase 2for k = 1:ndirsphi = (k-1)*pi/ndirs;y = real(exp(-i*phi)*m);[indmin,indmax,indzer] = extr(y);nem(k) = length(indmin)+length(indmax);nzm(k) = length(indzer);[tmin,tmax,zmin,zmax] = boundary_conditions(indmin,indmax,t,y,y,NBSYM);envmin(k,:) = exp(i*phi)*interp1(tmin,zmin,t,INTERP);envmax(k,:) = exp(i*phi)*interp1(tmax,zmax,t,INTERP);endenvmoy = mean((envmin+envmax),1);if nargout > 3amp = mean(abs(envmax-envmin),1)/2;endend
else[indmin,indmax,indzer] = extr(m);nem = length(indmin)+length(indmax);nzm = length(indzer);[tmin,tmax,mmin,mmax] = boundary_conditions(indmin,indmax,t,m,m,NBSYM);envmin = interp1(tmin,mmin,t,INTERP);envmax = interp1(tmax,mmax,t,INTERP);envmoy = (envmin+envmax)/2;if nargout > 3amp = mean(abs(envmax-envmin),1)/2;end
end
end%-------------------------------------------------------------------------------
% default stopping criterion
function [stop,envmoy,s] = stop_sifting(m,t,sd,sd2,tol,INTERP,MODE_COMPLEX,ndirs)
try[envmoy,nem,nzm,amp] = mean_and_amplitude(m,t,INTERP,MODE_COMPLEX,ndirs);sx = abs(envmoy)./amp;s = mean(sx);stop = ~((mean(sx > sd) > tol | any(sx > sd2)) & (all(nem > 2)));if ~MODE_COMPLEXstop = stop && ~(abs(nzm-nem)>1);end
catchstop = 1;envmoy = zeros(1,length(m));s = NaN;
end
end%-------------------------------------------------------------------------------
% stopping criterion corresponding to option FIX
function [stop,moyenne]= stop_sifting_fixe(t,m,INTERP,MODE_COMPLEX,ndirs)
trymoyenne = mean_and_amplitude(m,t,INTERP,MODE_COMPLEX,ndirs);stop = 0;
catchmoyenne = zeros(1,length(m));stop = 1;
end
end%-------------------------------------------------------------------------------
% stopping criterion corresponding to option FIX_H
function [stop,moyenne,stop_count]= stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H,MODE_COMPLEX,ndirs)
try[moyenne,nem,nzm] = mean_and_amplitude(m,t,INTERP,MODE_COMPLEX,ndirs);if (all(abs(nzm-nem)>1))stop = 0;stop_count = 0;elsestop_count = stop_count+1;stop = (stop_count == FIXE_H);end
catchmoyenne = zeros(1,length(m));stop = 1;
end
end%-------------------------------------------------------------------------------
% displays the progression of the decomposition with the default stopping criterion
function display_emd(t,m,mp,r,envmin,envmax,envmoy,s,sb,sx,sdt,sd2t,nbit,k,display_sifting,stop_sift)
subplot(4,1,1)
plot(t,mp);hold on;
plot(t,envmax,'--k');plot(t,envmin,'--k');plot(t,envmoy,'r');
title(['IMF ',int2str(k),';   iteration ',int2str(nbit),' before sifting']);
set(gca,'XTick',[])
hold  off
subplot(4,1,2)
plot(t,sx)
hold on
plot(t,sdt,'--r')
plot(t,sd2t,':k')
title('stop parameter')
set(gca,'XTick',[])
hold off
subplot(4,1,3)
plot(t,m)
title(['IMF ',int2str(k),';   iteration ',int2str(nbit),' after sifting']);
set(gca,'XTick',[])
subplot(4,1,4);
plot(t,r-m)
title('residue');
disp(['stop parameter mean value : ',num2str(sb),' before sifting and ',num2str(s),' after'])
if stop_siftdisp('last iteration for this mode')
end
if display_sifting == 2pause(0.01)
elsepause
end
end%---------------------------------------------------------------------------------------------------
% displays the progression of the decomposition with the FIX and FIX_H stopping criteria
function display_emd_fixe(t,m,mp,r,envmin,envmax,envmoy,nbit,k,display_sifting)
subplot(3,1,1)
plot(t,mp);hold on;
plot(t,envmax,'--k');plot(t,envmin,'--k');plot(t,envmoy,'r');
title(['IMF ',int2str(k),';   iteration ',int2str(nbit),' before sifting']);
set(gca,'XTick',[])
hold  off
subplot(3,1,2)
plot(t,m)
title(['IMF ',int2str(k),';   iteration ',int2str(nbit),' after sifting']);
set(gca,'XTick',[])
subplot(3,1,3);
plot(t,r-m)
title('residue');
if display_sifting == 2pause(0.01)
elsepause
end
end%---------------------------------------------------------------------------------------
% defines new extrema points to extend the interpolations at the edges of the
% signal (mainly mirror symmetry)
function [tmin,tmax,zmin,zmax] = boundary_conditions(indmin,indmax,t,x,z,nbsym)lx = length(x);if (length(indmin) + length(indmax) < 3)error('not enough extrema')end% boundary conditions for interpolations :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);elselmax = fliplr(indmax(1:min(end,nbsym)));lmin = [fliplr(indmin(1:min(end,nbsym-1))),1];lsym = 1;endelseif x(1) < x(indmax(1))lmax = fliplr(indmax(1:min(end,nbsym)));lmin = fliplr(indmin(2:min(end,nbsym+1)));lsym = indmin(1);elselmax = [fliplr(indmax(1:min(end,nbsym-1))),1];lmin = fliplr(indmin(1:min(end,nbsym)));lsym = 1;endendif 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);elsermax = [lx,fliplr(indmax(max(end-nbsym+2,1):end))];rmin = fliplr(indmin(max(end-nbsym+1,1):end));rsym = lx;endelseif 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);elsermax = fliplr(indmax(max(end-nbsym+1,1):end));rmin = [lx,fliplr(indmin(max(end-nbsym+2,1):end))];rsym = lx;endendtlmin = 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 enoughif tlmin(1) > t(1) || tlmax(1) > t(1)if lsym == indmax(1)lmax = fliplr(indmax(1:min(end,nbsym)));elselmin = fliplr(indmin(1:min(end,nbsym)));endif lsym == 1error('bug')endlsym = 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));elsermin = fliplr(indmin(max(end-nbsym+1,1):end));endif rsym == lxerror('bug')endrsym = lx;trmin = 2*t(rsym)-t(rmin);trmax = 2*t(rsym)-t(rmax);end zlmax =z(lmax); zlmin =z(lmin);zrmax =z(rmax); zrmin =z(rmin);tmin = [tlmin t(indmin) trmin];tmax = [tlmax t(indmax) trmax];zmin = [zlmin z(indmin) zrmin];zmax = [zlmax z(indmax) zrmax];
end%---------------------------------------------------------------------------------------------------
%extracts the indices of extrema
function [indmin, indmax, indzer] = extr(x,t)if(nargin==1)t=1:length(x);
endm = length(x);if nargout > 2x1=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);elseindz = iz;endindzer = sort([indzer indz]);end
endd = 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 successive points have the same value we consider only one extremum in the middle of the constant area
% (only works if the signal is uniformly sampled)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) == 1if length(debs) > 1debs = debs(2:end);fins = fins(2:end);elsedebs = [];fins = [];endendif length(debs) > 0if fins(end) == mif length(debs) > 1debs = debs(1:(end-1));fins = fins(1:(end-1));elsedebs = [];fins = [];endendendlc = length(debs);if lc > 0for k = 1:lcif d(debs(k)-1) > 0if d(fins(k)) < 0imax = [imax round((fins(k)+debs(k))/2)];endelseif d(fins(k)) > 0imin = [imin round((fins(k)+debs(k))/2)];endendendendif length(imax) > 0indmax = sort([indmax imax]);endif length(imin) > 0indmin = sort([indmin imin]);endend
end%---------------------------------------------------------------------------------------------------function ort = io(x,imf)
% ort = IO(x,imf) computes the index of orthogonality
%
% inputs : - x    : analyzed signal
%          - imf  : empirical mode decompositionn = size(imf,1);s = 0;for i = 1:nfor j =1:nif i~=js = s + abs(sum(imf(i,:).*conj(imf(j,:)))/sum(x.^2));endend
endort = 0.5*s;
end
%---------------------------------------------------------------------------------------------------function [x,t,sd,sd2,tol,MODE_COMPLEX,ndirs,display_sifting,sdt,sd2t,r,imf,k,nbit,NbIt,MAXITERATIONS,FIXE,FIXE_H,MAXMODES,INTERP,mask] = init(varargin)x = varargin{1};
if nargin == 2if isstruct(varargin{2})inopts = varargin{2};elseerror('when using 2 arguments the first one is the analyzed signal X and the second one is a struct object describing the options')end
elseif nargin > 2tryinopts = struct(varargin{2:end});catcherror('bad argument syntax')end
end% default for stopping
defstop = [0.05,0.5,0.05];opt_fields = {'t','stop','display','maxiterations','fix','maxmodes','interp','fix_h','mask','ndirs','complex_version'};defopts.stop = defstop;
defopts.display = 0;
defopts.t = 1:max(size(x));
defopts.maxiterations = 2000;
defopts.fix = 0;
defopts.maxmodes = 0;
defopts.interp = 'spline';
defopts.fix_h = 0;
defopts.mask = 0;
defopts.ndirs = 4;
defopts.complex_version = 2;opts = defopts;if(nargin==1)inopts = defopts;
elseif nargin == 0error('not enough arguments')
endnames = fieldnames(inopts);
for nom = names'if ~any(strcmpi(char(nom), opt_fields))error(['bad option field name: ',char(nom)])endif ~isempty(eval(['inopts.',char(nom)])) % empty values are discardedeval(['opts.',lower(char(nom)),' = inopts.',char(nom),';'])end
endt = opts.t;
stop = opts.stop;
display_sifting = opts.display;
MAXITERATIONS = opts.maxiterations;
FIXE = opts.fix;
MAXMODES = opts.maxmodes;
INTERP = opts.interp;
FIXE_H = opts.fix_h;
mask = opts.mask;
ndirs = opts.ndirs;
complex_version = opts.complex_version;if ~isvector(x)error('X must have only one row or one column')
endif size(x,1) > 1x = x.';
endif ~isvector(t)error('option field T must have only one row or one column')
endif ~isreal(t)error('time instants T must be a real vector')
endif size(t,1) > 1t = t';
endif (length(t)~=length(x))error('X and option field T must have the same length')
endif ~isvector(stop) || length(stop) > 3error('option field STOP must have only one row or one column of max three elements')
endif ~all(isfinite(x))error('data elements must be finite')
endif size(stop,1) > 1stop = stop';
endL = length(stop);
if L < 3stop(3)=defstop(3);
endif L < 2stop(2)=defstop(2);
endif ~ischar(INTERP) || ~any(strcmpi(INTERP,{'linear','cubic','spline'}))error('INTERP field must be ''linear'', ''cubic'', ''pchip'' or ''spline''')
end%special procedure when a masking signal is specified
if any(mask)if ~isvector(mask) || length(mask) ~= length(x)error('masking signal must have the same dimension as the analyzed signal X')endif size(mask,1) > 1mask = mask.';endopts.mask = 0;imf1 = emd(x+mask,opts);imf2 = emd(x-mask,opts);if size(imf1,1) ~= size(imf2,1)warning('emd:warning',['the two sets of IMFs have different sizes: ',int2str(size(imf1,1)),' and ',int2str(size(imf2,1)),' IMFs.'])endS1 = size(imf1,1);S2 = size(imf2,1);if S1 ~= S2if S1 < S2tmp = imf1;imf1 = imf2;imf2 = tmp;endimf2(max(S1,S2),1) = 0;endimf = (imf1+imf2)/2;endsd = stop(1);
sd2 = stop(2);
tol = stop(3);lx = length(x);sdt = sd*ones(1,lx);
sd2t = sd2*ones(1,lx);if FIXEMAXITERATIONS = FIXE;if FIXE_Herror('cannot use both ''FIX'' and ''FIX_H'' modes')end
endMODE_COMPLEX = ~isreal(x)*complex_version;
if MODE_COMPLEX && complex_version ~= 1 && complex_version ~= 2error('COMPLEX_VERSION parameter must equal 1 or 2')
end% number of extrema and zero-crossings in residual
ner = lx;
nzr = lx;r = x;if ~any(mask) % if a masking signal is specified "imf" already exists at this stageimf = [];
end
k = 1;% iterations counter for extraction of 1 mode
nbit=0;% total iterations counter
NbIt=0;
end

然后是eemd的,注意eemd的调用了emd的函数，所以eemd不要单独使用哦。

function [modos its]=eemd(x,Nstd,NR,MaxIter)
%--------------------------------------------------------------------------
%WARNING: this code needs to include in the same
%directoy the file emd.m developed by Rilling and Flandrin.
% -------------------------------------------------------------------------
%  OUTPUT
%   modos: contain the obtained modes in a matrix with the rows being the modes
%   its: contain the iterations needed for each mode for each realization
%
%  INPUT
%  x: signal to decompose
%  Nstd: noise standard deviation
%  NR: number of realizations
%  MaxIter: maximum number of sifting iterations allowed.
% -------------------------------------------------------------------------
%   Syntax
%
%   modos=eemd(x,Nstd,NR,MaxIter)
%  [modos its]=eemd(x,Nstd,NR,MaxIter)
% -------------------------------------------------------------------------
%  NOTE:   if Nstd=0 and NR=1, the EMD decomposition is obtained.
% -------------------------------------------------------------------------desvio_estandar=std(x);
x=x/desvio_estandar;
xconruido=x+Nstd*randn(size(x));
[modos, o, it]=emd(xconruido,'MAXITERATIONS',MaxIter);
modos=modos/NR;
iter=it;
if NR>=2for i=2:NRxconruido=x+Nstd*randn(size(x));[temp, ort, it]=emd(xconruido,'MAXITERATIONS',MaxIter);temp=temp/NR;lit=length(it);[p liter]=size(iter);if lit<literit=[it zeros(1,liter-lit)];end;if liter<lititer=[iter zeros(p,lit-liter)];end;iter=[iter;it];[filas columnas]=size(temp);[alto ancho]=size(modos);diferencia=alto-filas;if filas>altomodos=[modos; zeros(abs(diferencia),ancho)];end;if alto>filastemp=[temp;zeros(abs(diferencia),ancho)];end;modos=modos+temp;end;
end;
its=iter;
modos=modos*desvio_estandar;
end

VMD的脚本函数之前文章有说过，有需要的来这里粘贴吧，只粘贴vmd.m即可哦。

(52条消息) VMD分解，matlab代码，包络线，包络谱，中心频率，峭度值，能量熵，近似熵，包络熵，频谱图，希尔伯特变换，包含所有程序MATLAB代码，-西储大学数据集为例_包络谱绘制_今天吃饺子的博客-CSDN博客

也不能光上代码，结果图全部都在这里哦

首先是原始信号的图，原谅我太懒，没整坐标轴。这里简单说一下，从上往下依次是s1,s2,s3,s4,s，这里就是对s信号进行的分解哦。

然后是emd信号的分解图和频谱图

eemd的分解图和频谱图

VMD信号的分解图和频谱图

话说emd和eemd的差别看似不大哦，但是多少还是有点区别的，稍微降低了一些模态重叠。反观VMD在模态重叠这方面，确实无敌哦，基本无重叠。

okkkkkkkkkkkkkkkkkkkkkkk，祝大家代码永不出错！

emd,eemd,vmd,频谱图，分解图对比matlab代码相关推荐

Python实现“EMD\EEMD\VMD+Hilbert时频图”与“CWT小波时频图”
Python实现"EMD\EEMD\VMD+Hilbert时频图"与"CWT小波时频图" 信号处理中常需要分析时域统计量.频率成分,但不平稳信号的时域波形往 ...
matlab 频谱图例子_利用matlab怎样进行频谱分析
图像的频率是表征图像中灰度变化剧烈程度的指标,是灰度在平面空间上的梯度.如:大面积的沙漠在图像中是一片灰度变化缓慢的区域,对应的频率值很低:而对于地表属性变换剧烈的边缘区域在图像中是一片灰度变化剧烈的 ...
分酒问题matlab代码,matlab葡萄酒分类数据归一化问题
该楼层疑似违规已被系统折叠隐藏此楼查看此楼 matlab葡萄酒分类数据归一化问题% 选定训练集和测试集 % 将第一类的1-30,第二类的60-95,第三类的131-153做为训练集 train_wi ...
5种样式实现div容器中三图摆放实例对比说明
代码地址如下: http://www.demodashi.com/demo/11593.html 效果演示: demo点查看效果需求说明: 如下图所示为设计图,希望在图片上传无规则无规律的情况下实现 ...
VMD分解，matlab代码，包络线，包络谱，中心频率，峭度值，能量熵，近似熵，包络熵，频谱图，希尔伯特变换，包含所有程序MATLAB代码，-西储大学数据集为例
目录 1.选取数据 2.VMD函数-matlab代码 3.采用matlab脚本导入数据并做VMD分解 4. VMD分解图 5.计算中心频率 6.画包络线 7. 画包络谱 8. 计算峭度值 9.计算能量 ...
机器学习之MATLAB代码--CEEMDAN+EEMD+EMD+VMD+IMF重构络（十八）
机器学习之MATLAB代码--CEEMDAN+EEMD+EMD+VMD+IMF重构络(十八) 压缩分量的EEMD代码压缩分量的EEMD数据压缩分量的EEMD结果 CEEMDAN代码 CEEMDAN ...
如何绘制业务架构图 — 3.分解图
分解图:是对研究对象的有序分离.或是对细粒度要素的有序归集. 分解图是业务架构三视图的第二张图,其目的有两个:一是自上而下的"分解",二是自下而上的"汇集".但 ...
炫舞机器人内部零件图_汽车零件分解图高清图,求汽车内部零部件结构分
摩托车前减震结构图是什么样的? 摩托车前减震器结构图如图.小弹簧是在减震器的下端. 摩托车减震器分类: 1.根据安装位置分,有前减震器和后减震器; 2.按结构形式分,有(a)伸缩... 摩托车的化油 ...
信源编码作业（1）——绘制并分析清浊音频谱图
浊音以a.o为例绘制相应频谱图浊音时域频域 a o 结论: 1.时域波形呈现明显的周期特征. 2.频域在3kHz一下有多处明显波峰. 清音以p.t为例绘制相应频谱图清音时域频域 p t ...

emd,eemd,vmd,频谱图，分解图对比matlab代码

emd,eemd,vmd,频谱图，分解图对比matlab代码相关推荐

最新文章

热门文章