% function [allmode,newYsM, newYsm]=MA_EMD(Y,spikeIDMx,spikeIDMn,fs)

%
% This is an Minimum Arclength - EMD program
% Write and tested on Matlab 2016a and 2017a
%
% INPUT:
%   Y: Input data;1-d data only
%   spikeIDMx: index of the local maxima that are regarded as spike, 1xNsp vector
%   spikeIDmx: index of the local minima that are regarded as spike, 1xNsm vector
%   TNM : total number of modes

% OUTPUT:
%   allmode: IMFs in each column. The first column is the raw data. The second column (IMF1) is by MAMA-EMD
%   newYsM: the new points derived by minimum-arclength criterion (local max)
%         2*Nspike*10 matrix. 
%           (1, ~, ~) is the position (y-value)
%           (2, ~, ~) is to mark of the optimum point can be find (1: yes; 0: no)
%           if the optimum point can't be found by minimum-arclenth criteria,
%           the new osition is replaced by average of the adjacent local
%           maxima
%   newYsm: the new points derived by minimum-arclength criterion (local min)

% Code writer: Hui-Wen Yang 2019 @NTU and NCU
% Please cite:
% Hui-Wen Yang ; Shyh-Kang Jeng ; Hsu-Wen V. Young ; Chen Lin ; Yung-Hung Wang ; 
% Kun Hu ; Men-Tzung Lo. (2019). A Minimum Arclength Method for Removing Spikes 
% in Empirical Mode Decomposition. IEEE Access PP(99):1-1 ?January 2019

% This code is modified from the EEMD of Zhaohua Wu 2009
% References:   
%  Wu, Z., and N. E Huang (2008), 
%  Ensemble Empirical Mode Decomposition: a noise-assisted data analysis method. 
%   Advances in Adaptive Data Analysis. Vol.1, No.1. 1-41. 
function [allmode,newYsM, newYsm]=MA_EMD(Y,spikeIDMx,spikeIDMn,TNM)

if size(Y,1) > 1, Y = Y'; end

%part1.read data
xsize=length(Y);
dd=1:1:xsize;
%part2.evaluate TNM as total IMF number, assign 0 to TNM2 matrix
if nargin < 4
    TNM = 6;
end
TNM2=TNM+2;

allmode = zeros(xsize,TNM2);    %matrix to save the result IMFs
newYsM = zeros(2,size(spikeIDMx,1),10); % matrix to save the adjusted spike point
newYsm = zeros(2,size(spikeIDMn,1),10); 
% note that 10 represents number of sifting
% however, in Minimum-Arclength, we only do 4 siftings (sifgint#1,4,7,10) for efficiency

%part3 Do MA-EMD  -----EEMD loop start
NE = 1; % We keep the original EEMD structure but only do 1 ensemble (i.e. EMD)
for iii=1:1:NE,   
  
    %part4 --Add noise to original data,we have X1

    X1=Y;
    
    %part4 --assign original data in the first column  
    mode = zeros(xsize,TNM2);
    mode(:,1) = Y;
    
    %part5--give initial 0 to xorigin and xend
    xorigin = X1;
    xend = xorigin;
    
    %part6--start to find an IMF-----IMF loop start
    for nmode = 1:TNM,
        xstart = xend; %last loop value assign to new iteration loop 
                       %xstart -loop start data
        %part7--sift 10 times to get IMF---sift loop  start 
        iter = 1;
        while iter <= 10,
            [spmax, spmin, flag]=extrema(xstart);  %call function extrema 
%             display(['iteration' num2str(iter)])
        %% If nmode = 1 (the first IMF), we perfrom Minimum arclength quiteria
        if nmode == 1
            
            for ispike = 1:length(spikeIDMx)
                % find the real spike position
                [v,skipID2] = max(xstart(spikeIDMx(ispike)-5:spikeIDMx(ispike)+5)); 
                spikeIDMx(ispike)=spikeIDMx(ispike)-6+skipID2;
                spikeI = find(spmax(:,1)==spikeIDMx(ispike));
                
                if ~isempty(spikeI) % if the spike realy exists
                    % the left/right bound for calculating arclength
                    leftbnd = max(1,spikeI-5); rightbnd = min(spikeI+5,size(spmax,1));
                    newY = MinArcLength(spmax(leftbnd:rightbnd,1),spmax(leftbnd:rightbnd,2),6);
                    if isempty(newY)    % if the minimum arclength point does not exists
                        newY = mean(spmax([spikeI-1,spikeI+1],2));
                        newYsM(2,ispike,iter) = 0;
                    else
                        newYsM(2,ispike,iter) = 1;
                    end
                    spmax(spikeI,2) = newY;
                    newYsM(1,ispike,iter) = newY;
                end
            end
            for ispike = 1:length(spikeIDMn)
                % find the real spike position
                [v,skipID2] = min(xstart(spikeIDMn(ispike)-5:spikeIDMn(ispike)+5));
                spikeIDMn(ispike)=spikeIDMn(ispike)-6+skipID2;
                spikeI = find(spmin(:,1)==spikeIDMn(ispike));
                
                if ~isempty(spikeI)
                    % the left/right bound for calculating arclength
                    leftbnd = max(1,spikeI-5); rightbnd = min(spikeI+5,size(spmin,1));
                    newY = MinArcLength(spmin(leftbnd:rightbnd,1),spmin(leftbnd:rightbnd,2),6);
                    if isempty(newY)
                        newY = mean(spmin([spikeI-1,spikeI+1],2));
                        newYsm(2,ispike,iter) = 0;
                    else
                        newYsm(2,ispike,iter) = 1;
                    end
                    spmin(spikeI,2) = newY;
                    newYsm(1,ispike,iter) = newY;
                end
            end
            iter = iter+3;  % for the MA-sifting, we do only 1,4,7,10th sifting
        end
      
            %the usage of  spline ,please see part11.  
        %% Change original spline into adjusted spline
            upper= spline(spmax(:,1),spmax(:,2),dd); %upper spline bound of this sift 
            lower= spline(spmin(:,1),spmin(:,2),dd); %lower spline bound of this sift 

%%
            
            mean_ul = (upper + lower)/2;%spline mean of upper and lower  
            
            xstart = xstart - mean_ul;%extract spline mean from Xstart
            iter = iter + 1;
        end
        %part7--sift 10 times to get IMF---sift loop  end      
        
        %part8--subtract IMF from data ,then let the residual xend to start to find next IMF 
         xend = xend - xstart;
        
        %part9--after sift 10 times,that xstart is this time IMF 
        
        mode(:,nmode+1) = xstart;
        
    end
    %part6--start to find an IMF-----IMF loop end

    %part 10--after gotten  all(TNM) IMFs ,the residual xend is over all trend
    %                        put them in the last column  
    mode(:,nmode+2)=xend;
    
    %after part 10 ,original +TNM-IMF+overall trend  ---those are all in mode    
     allmode=allmode+mode;
    
end
%part3 Do EEMD  -----EEMD loop end

%part10--devide EEMD summation by NE,std be multiply back to data
allmode=allmode/NE;


%part11--the syntax of the matlab function spline
%yy= spline(x,y,xx); this means
%x and y are matrixs of n1 points ,use n1 set (x,y) to form the cubic spline
%xx and yy are matrixs of n2 points,we want know the spline value yy(y-axis) in the xx (x-axis)position
%after the spline is formed by n1 points ,find coordinate value on the spline for [xx,yy] --n2 position. 

