对图像进行emd分解程序得到各个imf分量

emdm文件

function imf = emd(x)

% Empiricial Mode Decomposition (Hilbert-Huang Transform)

% EMD分解或HHT变换

% 返回值为cell类型，依次为一次IMF、二次IMF、、最后残差

x = transpose(x(:));

imf = [];

while ~ismonotonic(x)

x1 = x;

sd = Inf;

while (sd > 01) || ~isimf(x1)

s1 = getspline(x1); % 极大值点样条曲线

s2 = -getspline(-x1); % 极小值点样条曲线

x2 = x1-(s1+s2)/2;

sd = sum((x1-x2)^2)/sum(x1^2);

x1 = x2;

end

imf{end+1} = x1;

x = x-x1;

end

imf{end+1} = x;

% 是否单调

function u = ismonotonic(x)

u1 = length(findpeaks(x))length(findpeaks(-x));

if u1 > 0

u = 0;

else

u = 1;

end

% 是否IMF分量

function u = isimf(x)

N = length(x);

u1 = sum(x(1:N-1)x(2:N) < 0); % 过零点的个数

u2 = length(findpeaks(x))+length(findpeaks(-x)); % 极值点的个数

if abs(u1-u2) > 1

u = 0;

else

u = 1;

end

% 据极大值点构造样条曲线

function s = getspline(x)

N = length(x);

p = findpeaks(x);

s = spline([0 p N+1],[0 x(p) 0],1:N);

这是对信号进行分解的程序，看看对你有没有帮助

用findpeaks函数

可以用后面的选项限制返回峰的大小和多少，除去一些因为噪声而产生的小峰

[] = findpeaks(x,'minpeakheight',mph) 峰值大于mph才返回[] = findpeaks(x,'minpeakdistance',mpd) 某峰前mpd个点和后mpd个点之间的峰忽略[] = findpeaks(x,'threshold',th) 与相邻值的差值大于th才返回

[] = findpeaks(x,'npeaks',np) 总共返回峰的个数[] = findpeaks(x,'sortstr',str) 按峰高排序

这些条件你可以自己根据数据选择，以滤除你不想要的峰

我看的你图，你应该只想要x轴范围在100~150里面的那个大的峰

所以，大概可以加个条件

[pks,locs] = findpeaks(xd,'minpeakheight',200,'sortstr','descend');

就会返回大于200的所有峰，而且峰高从大到小排列

plot(1:length(xd),xd);hold on;plot(ind(1),pks(1),'k');hold off; %返回第一个就是最高的

或者

[pks,locs] = findpeaks(xd,'minpeakdistance',30,'sortstr','descend');

plot(1:length(xd),xd);hold on;plot(ind(1),pks(1),'k');hold off; %返回第一个就是最高的

具体参数你可以自己调一下，你可数据是不够平滑

如果是找很大，很宽的峰，可以再适当平滑一下数据

你这样的数据多半是找到那个大峰上面偏右边的那个小突起

function imf = emd(x,n);%%最好把函数名改为emd1之类的，以免和Grilling的emd冲突

%%n为你想得到的IMF的个数

c = x('; % copy of the input signal (as a row vector)

N = length(x);-

% loop to decompose the input signal into n successive IMFs

imf = []; % Matrix which will contain the successive IMF, and the residuefor t=1:n

% loop on successive IMFs

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

% inner loop to find each imf

h = c; % at the beginning of the sifting process, h is the signal

SD = 1; % Standard deviation which will be used to stop the sifting process

while SD > 03 % while the standard deviation is higher than 03 (typical value) %%筛选停止准则

% find local max/min points

d = diff(h); % approximate derivative %%求各点导数

maxmin = []; % to store the optima (min and max without distinction so far)

for i=1:N-2

if d(i)==0 % we are on a zero %%导数为0的点，即”驻点“，但驻点不一定都是极值点，如y=x^3的x=0处

if sign(d(i-1))~=sign(d(i+1)) % it is a maximum %%如果驻点两侧的导数异号（如一边正，一边负），那么该点为极值点

maxmin = [maxmin, i]; %%找到极值点在信号中的坐标（不分极大值和极小值点）

end

elseif sign(d(i))~=sign(d(i+1)) % we are straddling a zero so%%如y=|x|在x=0处是极值点，但该点倒数不存在，所以不能用上面的判

断方法

maxmin = [maxmin, i+1]; % define zero as at i+1 (not i) %%这里提供了另一类极值点的判断方法

end

end

if size(maxmin,2) < 2 % then it is the residue %%判断信号是不是已经符合残余分量定义

break

end

% divide maxmin into maxes and mins %% 分离极大值点和极小值点

if maxmin(1)>maxmin(2) % first one is a max not a min

maxes = maxmin(1:2:length(maxmin));

mins = maxmin(2:2:length(maxmin));

else % is the other way around

maxes = maxmin(2:2:length(maxmin));

mins = maxmin(1:2:length(maxmin));

end % make endpoints both maxes and mins

maxes = [1 maxes N];

mins = [1 mins N];

%------------------------------------------------------------------------- % spline interpolate to get max and min envelopes; form imf

maxenv = spline(maxes,h(maxes),1:N); %%用样条函数插值拟合所有的极大值点

minenv = spline(mins, h(mins),1:N); %%用样条函数插值拟合所有的极小值点

m = (maxenv + minenv)/2; % mean of max and min enveloppes %%求上下包络的均值

prevh = h; % copy of the previous value of h before modifying it %%h为分解前的信号

h = h - m; % substract mean to h %% 减去包络均值

% calculate standard deviation

eps = 00000001; % to avoid zero values

SD = sum ( ((prevh - h)^2) / (prevh^2 + eps) ); %% 计算停止准则

end

imf = [imf; h]; % store the extracted IMF in the matrix imf

% if size(maxmin,2)<2, then h is the residue

% stop criterion of the algo if we reach the end before n

if size(maxmin,2) < 2

break

end

c = c - h; % substract the extracted IMF from the signal

end

return

%此版本为ALAN 版本的整合注释版

function imf = emd(x)

% Empiricial Mode Decomposition (Hilbert-Huang Transform)

% imf = emd(x)

% Func : findpeaks

x= transpose(x(:));%转置为行矩阵

imf = [];

while ~ismonotonic(x) %当x不是单调函数，分解终止条件

x1 = x;

sd = Inf;%均值

%直到x1满足IMF条件，得c1

while (sd > 01) | ~isimf(x1) %当标准偏差系数sd大于01或x1不是固有模态函数时，分量终止条件

s1 = getspline(x1);%上包络线

s2 = -getspline(-x1);%下包络线

x2 = x1-(s1+s2)/2;%此处的x2为文章中的h

sd = sum((x1-x2)^2)/sum(x1^2);

x1 = x2;

end

以上就是关于对图像进行emd分解程序得到各个imf分量全部的内容，包括:对图像进行emd分解程序得到各个imf分量、眼电信号用matlab求峰值问题，求大神、matlab编的有关EMD去噪的程序，处理的是核磁共振测井信号，急求等相关内容解答，如果想了解更多相关内容，可以关注我们，你们的支持是我们更新的动力！