function q = apitmemd(x, varargin)
% APIT-MEMD is an extension of the standard MEMD proposed by Rehman N. and
% Mandic D.P., 2010. This code has been modified based on the function MEMD (see below).
%
% function MEMD applies the "Multivariate Empirical Mode Decomposition" algorithm (Rehman and Mandic, Proc. Roy. Soc A, 2010)
% to multivariate inputs. We have verified this code by simulations for signals containing 3-16 channels.
%
% Syntax:
%
% imf = apitmemd(X)
% returns a 3D matrix 'imf(N,M,L)' containing M multivariate IMFs, one IMF per column, computed by applying
% the multivariate EMD algorithm on the N-variate signal (time-series) X of length L.
% - For instance, imf_k = IMF(k,:,:) returns the k-th component (1 <= k <= N) for all of the N-variate IMFs.
%
% For example, for hexavariate inputs (N=6), we obtain a 3D matrix IMF(6, M, L)
% where M is the number of IMFs extracted, and L is the data length.
%
% The default value of the alpha parameter is set 0.3
%
% imf = apitmemd(X,num_directions)
% where integer variable num_directions (>= 1) specifies the total number of projections of the signal
% - As a rule of thumb, the minimum value of num_directions should be twice the number of data channels,
% - for instance, num_directions = 6 for a 3-variate signal and num_directions= 16 for an 8-variate signal
% The default number of directions is chosen to be 64 - to extract meaningful IMFs, the number of directions
% should be considerably greater than the dimensionality of the signals
% The default value for variable alpha = 0.3 is used.
%
% imf = apitmemd(X, num_direction, 'alpha', alpha_var)
% Define the value of alpha variable which specifies the density of the relocated
% direction vectors; alpha=0, apitmemd operates similarly to standard memd
%
% imf = apitmemd(X,num_directions,'stopping criteria')
% uses the optional parameter 'stopping criteria' to control the sifting process.
% The available options are
% - 'stop' which uses the standard stopping criterion specified in [2]
% - 'fix_h' which uses the modified version of the stopping criteria specified in [3]
% The default value for the 'stopping criteria' is 'stop'.
%
% The settings num_directions=64 and 'stopping criteria' = 'stop' are defaults.
% Thus imf = MEMD(X) = MEMD(X,64) = MEMD(X,64,'stop') = MEMD(X,[],'stop'),
%
% imf = apitmemd(X, num_directions, 'stop', stop_vec)
% computes the IMFs based on the standard stopping criterion whose parameters are given in the 'stop_vec'
% - stop_vec has three elements specifying the threshold and tolerance values used, see [2].
% - the default value for the stopping vector is step_vec = [0.075 0.75 0.075].
% - the option 'stop_vec' is only valid if the parameter 'stopping criteria' is set to 'stop'.
%
% imf = apitmemd(X, num_directions, 'fix_h', n_iter)
% computes the IMFs with n_iter (integer variable) specifying the number of consecutive iterations when
% the number of extrema and the number of zero crossings differ at most by one [3].
% - the default value for the parameter n_iter is set to n_iter = 2.
% - the option n_iter is only valid if the parameter 'stopping criteria' = 'fix_h'
%
%
%
% This code allows to process unbalanced multivaraite signals having 3-16 channels, based on the multivariate EMD algorithm [1].
% - to perform MEMD on more than 16 channels, modify the variable 'Max_channels' on line 510 in the code accordingly.
% - to process 1- and 2-dimensional (univariate and bivariate) data using EMD, we recommend the toolbox from
% http://perso.ens-lyon.fr/patrick.flandrin/emd.html
%
% Acknowledgment: Part of this code is based on the bivariate EMD code, publicly available from
% http://perso.ens-lyon.fr/patrick.flandrin/emd.html. We would also like to thank
% Anh Huy Phan from RIKEN for helping us in optimizing the code and making it computationally efficient.
%
%
% Copyright: Naveed ur Rehman and Danilo P. Mandic, Oct-2009
%
%
% [1] Rehman and D. P. Mandic, "Multivariate Empirical Mode Decomposition", Proceedings of the Royal Society A, 2010
% [2] G. Rilling, P. Flandrin and P. Goncalves, "On Empirical Mode Decomposition and its Algorithms", Proc of the IEEE-EURASIP
% Workshop on Nonlinear Signal and Image Processing, NSIP-03, Grado (I), June 2003
% [3] N. E. Huang et al., "A confidence limit for the Empirical Mode Decomposition and Hilbert spectral analysis",
% Proceedings of the Royal Society A, Vol. 459, pp. 2317-2345, 2003
% [4] A. Hemakom, V. Goverdovsky, D. Looney and D. P. Mandic, "Adaptive-projection intrinsically-transformed multivariate
% empirical mode decomposition in cooperative brain-computer interface applications", Philosophical Transactions A, Vol. 374, No. 2065, pp. 20150199, 2016
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% Usage %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Case 1:
% inp = randn(1000,3);
% imf = memd(inp);
% imf_x = reshape(imf(1,:,:),size(imf,2),length(inp)); % imfs corresponding to 1st component
% imf_y = reshape(imf(2,:,:),size(imf,2),length(inp)); % imfs corresponding to 2nd component
% imf_z = reshape(imf(3,:,:),size(imf,2),length(inp)); % imfs corresponding to 3rd component
% Case 2:
% load syn_hex_inp.mat
% imf = memd(s6,256,'stop',[0.05 0.5 0.05])
global N N_dim;
[x, seq, t, ndir, N_dim, N, sd, sd2, tol, nbit, MAXITERATIONS, stop_crit, stp_cnt, alpha_var] = set_value(x, nargin, varargin{:});
r=x; n_imf=1;q = zeros(N_dim,1,N);
while ~stop_emd(r, seq, ndir,alpha_var)
% current mode
m = r;
% computation of mean and stopping criterion
if(strcmp(stop_crit,'stop'))
[stop_sift,env_mean] = stop_sifting(m,t,sd,sd2,tol,seq,ndir,alpha_var);
else
counter=0;
[stop_sift,env_mean,counter] = stop_sifting_fix(m,t,seq,ndir,stp_cnt,counter,alpha_var);
end
% In case the current mode is so small that machine precision can cause
% spurious extrema to appear
if (max(abs(m))) < (1e-10)*(max(abs(x)))
if ~stop_sift
warning('emd:warning','forced stop of EMD : too small amplitude')
else
disp('forced stop of EMD : too small amplitude')
end
break
end
% sifting loop
while ~stop_sift && nbit<MAXITERATIONS
%sifting
m = m - env_mean;
% computation of mean and stopping criterion
if(strcmp(stop_crit,'stop'))
[stop_sift,env_mean] = stop_sifting(m,t,sd,sd2,tol,seq,ndir,alpha_var);
else
[stop_sift,env_mean,counter] = stop_sifting_fix(m,t,seq,ndir,stp_cnt,counter,alpha_var);
end
nbit=nbit+1;
if(nbit==(MAXITERATIONS-1) && nbit > 100)
warning('emd:warning','forced stop of sifting : too many iterations');
end
end
q(:,n_imf,:)=m';
n_imf = n_imf+1;
r = r - m;
nbit = 0;
end
% Stores the residue
q(:,n_imf,:)=r';
%sprintf('Elapsed time: %f\n',toc);
end
%---------------------------------------------------------------------------------------------------
function stp = stop_emd(r, seq, ndir,alpha_var)
global N_dim;
ner = zeros(ndir,1);
for it=1:ndir
if (N_dim~=3) % Multivariate signal (for N_dim ~=3) with hammersley sequence
% Linear normalisation of hammersley sequence in the range of -1.00 - 1.00
b=2*seq(1:end,it)-1;
% Find angles corresponding to the normalised sequence
tht = atan2(sqrt(flipud(cumsum(b(N_dim:-1:2).^2))),b(1:N_dim-1)).';
% Find coordinates of unit direction vectors on n-sphere
dir_vec(1:N_dim) = [1 cumprod(sin(tht))];
dir_vec(1:N_dim-1) = cos(tht) .*dir_vec(1:N_dim-1);
else % Trivariate signal
