EMD and BEMD algorithm implementations (MATLAB code)

EMD and BEMD algorithm implementations (MATLAB code)

Today combined with their own on the Internet to find some implementation code, slightly modified to test, no more experiments, may be in some of the problem of processing or relatively coarse.
Emd

function IMF = EMD (x)% empiricial Mode decomposition (Hilbert-huang Transform)% EMD decomposition or HHT transform% return value to cell type, one time IMF, two times IMF,
    ..., last residuals x = Transpose (x (:)); IMF = []; while ~ismonotonic (x) x1 = x;
    SD = INF; while (SD > 0.1) | |         ~ISIMF (x1) S1 = getspline (x1);       % maximum point spline curve s2 =-getspline (-x1);

        % minimum point spline h = x1-(S1+S2)/2;
        SD = SUM ((x1-h). ^2)/sum (x1.^2);
    x1 = h;
    End imf{end+1} = x1;
x = x-x1;
End imf{end+1} = x;
End% is monotonic function u = ismonotonic (x) u1 = Length (findpeaks1 (x)) *length (FINDPEAKS1 (-X));
If u1 > 0 u = 0;
else U = 1;
End end% is the IMF component function U = ISIMF (x) N = length (x);                     u1 = SUM (x (1:n-1). *x (2:n) < 0);    % over 0 points u2 = Length (Findpeaks1 (x)) +length (FINDPEAKS1 ());
The number of extreme points if ABS (U1-U2) > 1 u = 0;
else U = 1;
End end% constructs a spline function s = getspline (x) N = Length (x), according to the maximum value point;
p = findpeaks1 (x);
s = spline ([0 P n+1],[0 x (P) 0],1:n); End Function N = fiNDPEAKS1 (x)% Find peaks. To find the maximum point, return the coordinate of the corresponding maximum point n = Find (diff (diff (x) > 0) < 0);
% is equivalent to finding the second-order guide less than 0 points u = Find (x (n+1) > x (n));  
N (U) = N (u) +1;



 End

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59

BEMD

function [Imf_matrix]=bemd (IMG) percent input a grayscale image [ROW,COL,DEP] = Size (img);% row, col and depth of original image if DEP ~= 1
img = im2double (Rgb2gray (IMG));  
else img = im2double (IMG);
End%%%%% Main function% decomposition the number of the IMF is set to 3 (plus a residual amount of 4 decomposition amount) (can be modified according to the actual situation) m=4;
k=1;
input_img=img;  while (k<m) [Imf_de res_de]=decompsition (INPUT_IMG); The amount of the IMF component and remainder Imf_matrix (:,, K) =imf_de;%% preserve the IMF component input_img=res_de by decomposition;
Percent of the remainder as a new signal, decomposition k=k+1 again; End Imf_matrix (:,:, K) =res_de;%% save residual End Function [Imf_de res_de]=decompsition (input_img) [Width height]=size (input_i
MG);
X=1:width;
Y=1:height;

input_img_temple=input_img;  while (1) [Zmax IMAX zmin imin]=extrema2 (input_img_temple);
    %%%% image Surface Extreme Point [Xmax ymax]=ind2sub (Size (input_img_temple), IMAX);

    [Xmin ymin]=ind2sub (Size (input_img_temple), imin);  [Zmaxgrid,~,~]=gridfit (YMAX,XMAX,ZMAX,Y,X);


    %%%% surface fitting to find the extreme point of the envelope surface [Zmingrid,~,~]=gridfit (ymin,xmin,zmin,y,x);   Zavggrid= (Zmaxgrid+zmingrid)/2; %%%% envelope mean value%%%%%%IMF component judgment%%%%% imf_de=Input_img_temple-zavggrid;

    Sd=sum (sum (imf_de-input_img_temple). ^2)/sum (sum (imf_de). ^2);
    If sd<0.2 break else input_img_temple=imf_de;

End End Res_de=input_img-imf_de; End

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59

Note: Bemd uses EXTREMA2 to find the surface extremum and the gridfit surface fitting function to achieve the envelope surface acquisition.