Home > FA-MVEMD > 3D > EMD3D3V_parallel_var.m

# EMD3D3V_parallel_var

## PURPOSE

Purpose:

## SYNOPSIS

function Results = EMD3D3V_parallel_var(u,v,w,param)

## DESCRIPTION

     Purpose:
     -To perform EMD on 3 channels of 3 dimensional data

     Input:
     - u: Signal 1
     - v: Signal 2
     - w: Signal 3
     - param.
       -nimfs: Number of IMFs to be extracted
       -tol: Sifting tolerance value
       -type: type of window size to be used
       -plot: 'on' to plot results, default hides IMF plots
       -nslice: number of slices in volume plot
       -xend: Domain length in x (meters)
       -yend: Domain length in y (meters)
       -zend: Domain length in z (meters)

     Output:
     - Results
       - IMF (structure containing IMFs of all three signals)
       - Residue (structure containing residue of all three signals)
       - Windows (Window sizes (5 types) for each IMF)
       - Sift_cnt (Number of sifting iterations for each signal)
       - IO (Index of orthogonality for each signal)
       - Error (Error of the decomposition for each signal)

     References:

     Written by Mruthun Thirumalaisamy
     Graduate Student
     Department of Aerospace Engineering
     University of Illinois at Urbana-Champaign
     May 16 2018
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

## CROSS-REFERENCE INFORMATION

This function calls:

- MinimaMaxima3D V 1.0 Dec 13, 07
- TIMF_plot default plot attributes

This function is called by:

## SUBFUNCTIONS

- function Windows = filter_size_var(maxima_pos, minima_pos, param, Dmnsize, type)
- function Windows = filter_size(maxima_pos, minima_pos,type)
- function H1 = Sift(H,w_sz)
- function [Max,Min] = OSF(H,w_sz)
- function Signal = Separable_ordfilt3(Signal, order, w_sz)
- function f_signal = Ordfilt1(signal,order,window_size)
- function Env_med = Pad_smooth(Env_max,Env_min,w_sz)
- function [IO,Error] = Orth_index(Signal,IMF,Residue)
- function Plot_results(u,v,w,Results,param)
- function TIMF_plot(signal,Colour,nslice,imf,name1,name2)

## SOURCE CODE

    0001 % Purpose:
    0002 % -To perform EMD on 3 channels of 3 dimensional data
    0003 %
    0004 % Input:
    0005 % - u: Signal 1
    0006 % - v: Signal 2
    0007 % - w: Signal 3
    0008 % - param.
    0009 %   -nimfs: Number of IMFs to be extracted
    0010 %   -tol: Sifting tolerance value
    0011 %   -type: type of window size to be used
    0012 %   -plot: 'on' to plot results, default hides IMF plots
    0013 %   -nslice: number of slices in volume plot
    0014 %   -xend: Domain length in x (meters)
    0015 %   -yend: Domain length in y (meters)
    0016 %   -zend: Domain length in z (meters)
    0017 %
    0018 % Output:
    0019 % - Results
    0020 %   - IMF (structure containing IMFs of all three signals)
    0021 %   - Residue (structure containing residue of all three signals)
    0022 %   - Windows (Window sizes (5 types) for each IMF)
    0023 %   - Sift_cnt (Number of sifting iterations for each signal)
    0024 %   - IO (Index of orthogonality for each signal)
    0025 %   - Error (Error of the decomposition for each signal)
    0026 %
    0027 % References:
    0028 %
    0029 %
    0030 % Written by Mruthun Thirumalaisamy
    0031 % Graduate Student
    0032 % Department of Aerospace Engineering
    0033 % University of Illinois at Urbana-Champaign
    0034 % May 16 2018
    0035 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0036
    0037 function Results = EMD3D3V_parallel_var(u,v,w,param)
    0038
    0039 %Reading signal characteristics
    0040 [Nx,Ny,Nz] = size(u); %Signal dimensions
    0041 B          = size(v); %Signal dimensions
    0042 C          = size(w); %Signal dimensions
    0043
    0044 Dmnsize = [Nx,Ny,Nz];
    0045
    0046 %Preliminary checks
    0047 if ~isfield(param,'nimfs')
    0048     param.nimfs = 10;
    0049 end
    0050
    0051 if ~isfield(param,'tol')
    0052     param.tol = 0.05; % 0.1% of the minimum signal amplitude
    0053 end
    0054
    0055 if ~isfield(param,'type')
    0056     param.type = 6;
    0057 end
    0058
    0059 if ~isfield(param,'plot')
    0060     param.plot = 'off';
    0061 end
    0062
    0063 if(~all(ismember(param.type,[1,2,3,4,5,6,7])))
    0064     error('Please enter a valid window size type')
    0065 end
    0066
    0067 if(~all([Nx,Ny,Nz]==B) ||  ~all([Nx,Ny,Nz]==C))
    0068     error('Inconsistent dimensions between channels. Please check input data');
    0069 end
    0070 clearvars B C
    0071
    0072 if(param.tol<=0.005)
    0073    warning('Low sifting tolerance may cause oversifting');
    0074    answer = questdlg('Would you like to continue?', ...
    0075     'User set low sifting tolerance', ...
    0076     'Yes','No','No');
    0077     % Handle response
    0078     switch answer
    0079         case 'Yes'
    0080
    0081         case 'No'
    0082             return;
    0083     end
    0084 end
    0085
    0086 if ~isfield(param,'plot')
    0087     param.plot = 'off';
    0088 end
    0089
    0090 %Initialisations
    0091 IMF.u = zeros(Nx, Ny, Nz, param.nimfs);
    0092 IMF.v = zeros(Nx, Ny, Nz, param.nimfs);
    0093 IMF.w = zeros(Nx, Ny, Nz, param.nimfs);
    0094 Residue.u = u; Residue.v = v; Residue.w = w;
    0095
    0096 H = zeros(Nx,Ny,Nz,3);
    0097 H1 = zeros(Nx,Ny,Nz,3);
    0098 mse = zeros(3,1);
    0099
    0100 Windows = zeros(7,3,param.nimfs);
    0101
    0102 sift_cnt = zeros(1,param.nimfs);
    0103 imf = 1;
    0104 stopflag = 1;
    0105
    0106     while(imf <= param.nimfs && stopflag)
    0107         %Initialising intermediary IMFs
    0108         H(:,:,:,1) = Residue.u; H(:,:,:,2) = Residue.v; H(:,:,:,3) = Residue.w;
    0109
    0110         sift_stop = 0; %flag to control sifting loop
    0111
    0112         Combined = H(:,:,:,1)/sqrt(3) + H(:,:,:,2)/sqrt(3) + H(:,:,:,3)/sqrt(3); %Combining three signals with equal weights
    0113         [Maxima,MaxPos,Minima,MinPos] = MinimaMaxima3D(Combined,1,1,[],[]);  %Obtaining extrema of combined signal
    0114
    0115         %Checking whether there are too few extrema in the IMF
    0116         if (nnz(Maxima) < 3 || nnz(Minima) < 3)
    0117             warning('Fewer than three extrema found in extrema map. Stopping now...');
    0118             break;
    0119         end
    0120
    0121         %Window size determination by delaunay triangulation
    0122 %          Windows(:,imf) = filter_size(MaxPos,MinPos,param.type);
    0123         Windows(:,:,imf) = filter_size_var(MaxPos,MinPos,param,Dmnsize,param.type);
    0124         w_sz = Windows(param.type,:,imf); %extracting window size chosen by input parameter
    0125
    0126
    0127         if~(any(w_sz))
    0128            warning('EMD3D3V has stopped because the Delaunay Triangulation could not be created (collinear points)');
    0129            stopflag = 0; %#ok<NASGU>
    0130            break;
    0131         end
    0132
    0133         %Begin sifting iteration
    0134         while~(sift_stop)
    0135             sift_cnt(imf) = sift_cnt(imf) + 1; %Incrementing sift counter
    0136
    0137             %Entering parallel sift calculations
    0138
    0139             parfor i=1:3
    0140                H1(:,:,:,i) = Sift(H(:,:,:,i),w_sz);
    0141
    0142                mse(i) = immse(H1(:,:,:,i),H(:,:,:,i));
    0143             end
    0144
    0145             %Stop condition checks
    0146             if (mse(1)<param.tol && mse(2)<param.tol && mse(3)<param.tol && sift_cnt(imf)~=1)
    0147                 sift_stop = 1;
    0148             end
    0149
    0150             H(:,:,:,1) = H1(:,:,:,1); H(:,:,:,2) = H1(:,:,:,2); H(:,:,:,3) = H1(:,:,:,3);
    0151         end
    0152
    0153         %Storing IMFs
    0154         IMF.u(:,:,:,imf) = H(:,:,:,1); IMF.v(:,:,:,imf) = H(:,:,:,2); IMF.w(:,:,:,imf) = H(:,:,:,3);
    0155
    0156         %Subtracting from Residual Signals
    0157         Residue.u = Residue.u - IMF.u(:,:,:,imf);
    0158         Residue.v = Residue.v - IMF.v(:,:,:,imf);
    0159         Residue.w = Residue.w - IMF.w(:,:,:,imf);
    0160
    0161         %Incrementing IMF counter
    0162         imf = imf + 1;
    0163
    0164     end
    0165
    0166     %Checking for oversifting
    0167     if(any(sift_cnt>=5*ones(size(sift_cnt))))
    0168         warning('Decomposition may be oversifted. Checking if window size increases monotonically...');
    0169
    0170         if( any (diff(Windows(param.type,:)) <= zeros(1,param.nimfs-1)) )
    0171         warning('Filter window size does not increase monotonically')
    0172         end
    0173     end
    0174
    0175     %Organising results
    0176     Results.IMF = IMF;
    0177     Results.Residue = Residue;
    0178     Results.Windows = Windows;
    0179     Results.Sifts = sift_cnt;
    0180
    0181     %Error and orthogonality
    0182     [Results.IO.u,Results.Error.u] = Orth_index(u,IMF.u,Residue.u);
    0183     [Results.IO.v,Results.Error.v] = Orth_index(v,IMF.v,Residue.v);
    0184     [Results.IO.w,Results.Error.w] = Orth_index(w,IMF.w,Residue.w);
    0185
    0186     switch(param.plot)
    0187         case 'on'
    0188             Plot_results(u,v,w,Results,param)
    0189     end
    0190 end
    0191
    0192 function Windows = filter_size_var(maxima_pos, minima_pos, param, Dmnsize, type)
    0193 % Purpose:
    0194 % -To determine the window size for order statistics filtering of a signal.
    0195 %
    0196 % Inputs:
    0197 % -Two matrices of extrema positions
    0198 % -Parameter array
    0199 % -Domain size array
    0200 %
    0201 % Outputs:
    0202 % -Calculated value of the window size
    0203 %
    0204 %
    0205 % -Our focus is to ensure physical significance of the resulting decomposition
    0206 % -We can take advantage of the fact that the sample rate in each direction is known \textit{a priori}
    0207 % -After obtaining an extrema map (index locations of the extrema), calculate the Delaunay Triangulation for this
    0208 % -The Euclidean distances between the nearest neighbours is calculated in the Cartesian co-ordinate space and
    0209 %  not the index location space
    0210 % -This results in a window size that is in metres and not indices (or sample length)
    0211 % -Now, we can multiply this window size with the sampling rate in each direction to get respective window sizes
    0212 %  in terms of sample lengths!
    0213 %
    0214 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0215
    0216 %generating the grid information
    0217 xcoord = linspace(0,param.xend,Dmnsize(1))';
    0218 ycoord = linspace(0,param.yend,Dmnsize(2))';
    0219 zcoord = linspace(0,param.zend,Dmnsize(3))';
    0220
    0221 %calculating the sampling rates
    0222 fsx = Dmnsize(1)/(xcoord(end)-xcoord(1));
    0223 fsy = Dmnsize(2)/(ycoord(end)-ycoord(1));
    0224 fsz = Dmnsize(3)/(zcoord(end)-zcoord(1));
    0225
    0226 %use delaunay triangulation to determine the nearest neighbours and hence
    0227 %the filter size
    0228
    0229 %processing d_max
    0230 max_nearest = zeros(length(maxima_pos),1);
    0231
    0232 try
    0233 TRI_max = delaunay(maxima_pos);
    0234 catch
    0235     warning('Maxima points are collinear. Exiting without further iterations');
    0236     Windows = [[0, 0, 0, 0, 0, 0, 0];[0, 0, 0, 0, 0, 0, 0];[0, 0, 0, 0, 0, 0, 0]];
    0237     return
    0238 end
    0239
    0240
    0241 maxima_pos_x = xcoord(maxima_pos(:,1));
    0242 maxima_pos_y = ycoord(maxima_pos(:,2));
    0243 maxima_pos_z = zcoord(maxima_pos(:,3));
    0244
    0245 %Calculating 6 edge lengths for each tetrahedron
    0246 e1 = sqrt( (maxima_pos_x(TRI_max(:,2))- maxima_pos_x(TRI_max(:,1))).^2 + (maxima_pos_y(TRI_max(:,2))- maxima_pos_y(TRI_max(:,1))).^2 + (maxima_pos_z(TRI_max(:,2))- maxima_pos_z(TRI_max(:,1))).^2 );
    0247 e2 = sqrt( (maxima_pos_x(TRI_max(:,3))- maxima_pos_x(TRI_max(:,1))).^2 + (maxima_pos_y(TRI_max(:,3))- maxima_pos_y(TRI_max(:,1))).^2 + (maxima_pos_z(TRI_max(:,3))- maxima_pos_z(TRI_max(:,1))).^2 );
    0248 e3 = sqrt( (maxima_pos_x(TRI_max(:,3))- maxima_pos_x(TRI_max(:,2))).^2 + (maxima_pos_y(TRI_max(:,3))- maxima_pos_y(TRI_max(:,2))).^2 + (maxima_pos_z(TRI_max(:,3))- maxima_pos_z(TRI_max(:,2))).^2 );
    0249 e4 = sqrt( (maxima_pos_x(TRI_max(:,4))- maxima_pos_x(TRI_max(:,1))).^2 + (maxima_pos_y(TRI_max(:,4))- maxima_pos_y(TRI_max(:,1))).^2 + (maxima_pos_z(TRI_max(:,4))- maxima_pos_z(TRI_max(:,1))).^2 );
    0250 e5 = sqrt( (maxima_pos_x(TRI_max(:,4))- maxima_pos_x(TRI_max(:,2))).^2 + (maxima_pos_y(TRI_max(:,4))- maxima_pos_y(TRI_max(:,2))).^2 + (maxima_pos_z(TRI_max(:,4))- maxima_pos_z(TRI_max(:,2))).^2 );
    0251 e6 = sqrt( (maxima_pos_x(TRI_max(:,4))- maxima_pos_x(TRI_max(:,3))).^2 + (maxima_pos_y(TRI_max(:,4))- maxima_pos_y(TRI_max(:,3))).^2 + (maxima_pos_z(TRI_max(:,4))- maxima_pos_z(TRI_max(:,3))).^2 );
    0252
    0253 %Calculating nearest neighbours for each maxima point
    0254 %Comparing tetrahedron edges associated with each vertex
    0255 em1 = min([e1, e2, e4],[],2); %Comparing edges 1, 2 and 4 (vertex 1)
    0256 em2 = min([e1, e3, e5],[],2); %Comparing edges 1, 3 and 5 (vertex 2)
    0257 em3 = min([e2, e3, e6],[],2); %Comparing edges 2, 3 and 6 (vertex 3)
    0258 em4 = min([e4, e5, e6],[],2); %Comparing edges 4, 5 and 6 (vertex 4)
    0259
    0260 e = [em1 ,em2, em3, em4];
    0261
    0262 %Making sure that the smallest edge associated with the each vertex is stored
    0263 %correctly
    0264 for i=1:length(em1)
    0265     for j=1:4
    0266         if max_nearest(TRI_max(i,j)) > e(i,j) || max_nearest(TRI_max(i,j)) == 0
    0267             max_nearest(TRI_max(i,j)) = e(i,j);
    0268         end
    0269     end
    0270 end
    0271
    0272 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0273 %processing d_min
    0274 min_nearest = zeros(length(minima_pos),1);
    0275
    0276 try
    0277 TRI_min = delaunay(minima_pos);
    0278 catch
    0279     warning('Minima points are collinear. Exiting without further iterations');
    0280     Windows = [[0, 0, 0, 0, 0, 0, 0];[0, 0, 0, 0, 0, 0, 0];[0, 0, 0, 0, 0, 0, 0]];
    0281     return
    0282 end
    0283 minima_pos_x = xcoord(minima_pos(:,1));
    0284 minima_pos_y = ycoord(minima_pos(:,2));
    0285 minima_pos_z = zcoord(minima_pos(:,3));
    0286
    0287 %Calculating 6 edge lengths for each tetrahedron
    0288 e1 = sqrt( (minima_pos_x(TRI_min(:,2))- minima_pos_x(TRI_min(:,1))).^2 + (minima_pos_y(TRI_min(:,2))- minima_pos_y(TRI_min(:,1))).^2 + (minima_pos_z(TRI_min(:,2))- minima_pos_z(TRI_min(:,1))).^2 );
    0289 e2 = sqrt( (minima_pos_x(TRI_min(:,3))- minima_pos_x(TRI_min(:,1))).^2 + (minima_pos_y(TRI_min(:,3))- minima_pos_y(TRI_min(:,1))).^2 + (minima_pos_z(TRI_min(:,3))- minima_pos_z(TRI_min(:,1))).^2 );
    0290 e3 = sqrt( (minima_pos_x(TRI_min(:,3))- minima_pos_x(TRI_min(:,2))).^2 + (minima_pos_y(TRI_min(:,3))- minima_pos_y(TRI_min(:,2))).^2 + (minima_pos_z(TRI_min(:,3))- minima_pos_z(TRI_min(:,2))).^2 );
    0291 e4 = sqrt( (minima_pos_x(TRI_min(:,4))- minima_pos_x(TRI_min(:,1))).^2 + (minima_pos_y(TRI_min(:,4))- minima_pos_y(TRI_min(:,1))).^2 + (minima_pos_z(TRI_min(:,4))- minima_pos_z(TRI_min(:,1))).^2 );
    0292 e5 = sqrt( (minima_pos_x(TRI_min(:,4))- minima_pos_x(TRI_min(:,2))).^2 + (minima_pos_y(TRI_min(:,4))- minima_pos_y(TRI_min(:,2))).^2 + (minima_pos_z(TRI_min(:,4))- minima_pos_z(TRI_min(:,2))).^2 );
    0293 e6 = sqrt( (minima_pos_x(TRI_min(:,4))- minima_pos_x(TRI_min(:,3))).^2 + (minima_pos_y(TRI_min(:,4))- minima_pos_y(TRI_min(:,3))).^2 + (minima_pos_z(TRI_min(:,4))- minima_pos_z(TRI_min(:,3))).^2 );
    0294
    0295 %Calculating nearest neighbours for each minima point
    0296 %Comparing tetrahedron edges associated with each vertex
    0297 emn1 = min([e1, e2, e4],[],2); %Comparing edges 1, 2 and 4 (vertex 1)
    0298 emn2 = min([e1, e3, e5],[],2); %Comparing edges 1, 3 and 5 (vertex 2)
    0299 emn3 = min([e2, e3, e6],[],2); %Comparing edges 2, 3 and 6 (vertex 3)
    0300 emn4 = min([e4, e5, e6],[],2); %Comparing edges 4, 5 and 6 (vertex 4)
    0301
    0302 e = [emn1 ,emn2, emn3, emn4];
    0303
    0304 %Making sure that the smallest edge associated with the each vertex is stored
    0305 %correctly
    0306 for i=1:length(emn1)
    0307     for j=1:4
    0308         if min_nearest(TRI_min(i,j)) > e(i,j) || min_nearest(TRI_min(i,j)) == 0
    0309             min_nearest(TRI_min(i,j)) = e(i,j);
    0310         end
    0311     end
    0312 end
    0313
    0314 %Window size calculations
    0315 d1 = min( min(max_nearest) , min(min_nearest) );
    0316 d2 = max( min(max_nearest) , min(min_nearest) );
    0317 d3 = min( max(max_nearest) , max(min_nearest) );
    0318 d4 = max( max(max_nearest) , max(min_nearest) );
    0319 d5 = (d1+d2+d3+d4)/4 ;
    0320 d6 = median([min_nearest; max_nearest]);
    0321 d7 = mode([min_nearest; max_nearest]);
    0322
    0323 Windows = [ fsx*[d1, d2, d3, d4, d5, d6, d7];...
    0324             fsy*[d1, d2, d3, d4, d5, d6, d7];...
    0325             fsz*[d1, d2, d3, d4, d5, d6, d7]    ]';
    0326
    0327 %making sure w_size is an odd integer
    0328 Windows = 2*(floor(Windows./2))+1;
    0329
    0330 check = Windows(type,:)<3;
    0331
    0332 if(any(check))
    0333     warning('WARNING: Calculated Window size less than 3');
    0334     warning('Overriding calculated value and setting window size = 3');
    0335     Windows(type,check) = 3;
    0336 end
    0337 end
    0338
    0339 function Windows = filter_size(maxima_pos, minima_pos,type)
    0340 % Purpose:
    0341 % -To determine the window size for order statistics filtering of a signal.
    0342 % The determination of the window size is based on the work of Bhuiyan et
    0343 % al.
    0344 %
    0345 % Inputs:
    0346 % -Two matrices of extrema positions
    0347 %
    0348 % Outputs:
    0349 % -Calculated value of the window size
    0350 %
    0351 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0352
    0353 %use delaunay triangulation to determine the nearest neighbours and hence
    0354 %the filter size
    0355
    0356 %processing d_max
    0357 max_nearest = zeros(length(maxima_pos),1);
    0358
    0359 try
    0360 TRI_max = delaunay(maxima_pos);
    0361 catch
    0362     warning('Maxima points are collinear. Exiting without further iterations');
    0363     Windows = [0, 0, 0, 0, 0, 0, 0];
    0364     return
    0365 end
    0366
    0367
    0368 maxima_pos_x = maxima_pos(:,1);
    0369 maxima_pos_y = maxima_pos(:,2);
    0370 maxima_pos_z = maxima_pos(:,3);
    0371
    0372 %Calculating 6 edge lengths for each tetrahedron
    0373 e1 = sqrt( (maxima_pos_x(TRI_max(:,2))- maxima_pos_x(TRI_max(:,1))).^2 + (maxima_pos_y(TRI_max(:,2))- maxima_pos_y(TRI_max(:,1))).^2 + (maxima_pos_z(TRI_max(:,2))- maxima_pos_z(TRI_max(:,1))).^2 );
    0374 e2 = sqrt( (maxima_pos_x(TRI_max(:,3))- maxima_pos_x(TRI_max(:,1))).^2 + (maxima_pos_y(TRI_max(:,3))- maxima_pos_y(TRI_max(:,1))).^2 + (maxima_pos_z(TRI_max(:,3))- maxima_pos_z(TRI_max(:,1))).^2 );
    0375 e3 = sqrt( (maxima_pos_x(TRI_max(:,3))- maxima_pos_x(TRI_max(:,2))).^2 + (maxima_pos_y(TRI_max(:,3))- maxima_pos_y(TRI_max(:,2))).^2 + (maxima_pos_z(TRI_max(:,3))- maxima_pos_z(TRI_max(:,2))).^2 );
    0376 e4 = sqrt( (maxima_pos_x(TRI_max(:,4))- maxima_pos_x(TRI_max(:,1))).^2 + (maxima_pos_y(TRI_max(:,4))- maxima_pos_y(TRI_max(:,1))).^2 + (maxima_pos_z(TRI_max(:,4))- maxima_pos_z(TRI_max(:,1))).^2 );
    0377 e5 = sqrt( (maxima_pos_x(TRI_max(:,4))- maxima_pos_x(TRI_max(:,2))).^2 + (maxima_pos_y(TRI_max(:,4))- maxima_pos_y(TRI_max(:,2))).^2 + (maxima_pos_z(TRI_max(:,4))- maxima_pos_z(TRI_max(:,2))).^2 );
    0378 e6 = sqrt( (maxima_pos_x(TRI_max(:,4))- maxima_pos_x(TRI_max(:,3))).^2 + (maxima_pos_y(TRI_max(:,4))- maxima_pos_y(TRI_max(:,3))).^2 + (maxima_pos_z(TRI_max(:,4))- maxima_pos_z(TRI_max(:,3))).^2 );
    0379
    0380 %Calculating nearest neighbours for each maxima point
    0381 %Comparing tetrahedron edges associated with each vertex
    0382 em1 = min([e1, e2, e4],[],2); %Comparing edges 1, 2 and 4 (vertex 1)
    0383 em2 = min([e1, e3, e5],[],2); %Comparing edges 1, 3 and 5 (vertex 2)
    0384 em3 = min([e2, e3, e6],[],2); %Comparing edges 2, 3 and 6 (vertex 3)
    0385 em4 = min([e4, e5, e6],[],2); %Comparing edges 4, 5 and 6 (vertex 4)
    0386
    0387 e = [em1 ,em2, em3, em4];
    0388
    0389 %Making sure that the smallest edge associated with the each vertex is stored
    0390 %correctly
    0391 for i=1:length(em1)
    0392     for j=1:4
    0393         if max_nearest(TRI_max(i,j)) > e(i,j) || max_nearest(TRI_max(i,j)) == 0
    0394             max_nearest(TRI_max(i,j)) = e(i,j);
    0395         end
    0396     end
    0397 end
    0398
    0399 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0400 %processing d_min
    0401 min_nearest = zeros(length(minima_pos),1);
    0402
    0403 try
    0404 TRI_min = delaunay(minima_pos);
    0405 catch
    0406     warning('Minima points are collinear. Exiting without further iterations');
    0407     Windows = [0, 0, 0, 0, 0, 0, 0];
    0408     return
    0409 end
    0410 minima_pos_x = minima_pos(:,1);
    0411 minima_pos_y = minima_pos(:,2);
    0412 minima_pos_z = minima_pos(:,3);
    0413
    0414 %Calculating 6 edge lengths for each tetrahedron
    0415 e1 = sqrt( (minima_pos_x(TRI_min(:,2))- minima_pos_x(TRI_min(:,1))).^2 + (minima_pos_y(TRI_min(:,2))- minima_pos_y(TRI_min(:,1))).^2 + (minima_pos_z(TRI_min(:,2))- minima_pos_z(TRI_min(:,1))).^2 );
    0416 e2 = sqrt( (minima_pos_x(TRI_min(:,3))- minima_pos_x(TRI_min(:,1))).^2 + (minima_pos_y(TRI_min(:,3))- minima_pos_y(TRI_min(:,1))).^2 + (minima_pos_z(TRI_min(:,3))- minima_pos_z(TRI_min(:,1))).^2 );
    0417 e3 = sqrt( (minima_pos_x(TRI_min(:,3))- minima_pos_x(TRI_min(:,2))).^2 + (minima_pos_y(TRI_min(:,3))- minima_pos_y(TRI_min(:,2))).^2 + (minima_pos_z(TRI_min(:,3))- minima_pos_z(TRI_min(:,2))).^2 );
    0418 e4 = sqrt( (minima_pos_x(TRI_min(:,4))- minima_pos_x(TRI_min(:,1))).^2 + (minima_pos_y(TRI_min(:,4))- minima_pos_y(TRI_min(:,1))).^2 + (minima_pos_z(TRI_min(:,4))- minima_pos_z(TRI_min(:,1))).^2 );
    0419 e5 = sqrt( (minima_pos_x(TRI_min(:,4))- minima_pos_x(TRI_min(:,2))).^2 + (minima_pos_y(TRI_min(:,4))- minima_pos_y(TRI_min(:,2))).^2 + (minima_pos_z(TRI_min(:,4))- minima_pos_z(TRI_min(:,2))).^2 );
    0420 e6 = sqrt( (minima_pos_x(TRI_min(:,4))- minima_pos_x(TRI_min(:,3))).^2 + (minima_pos_y(TRI_min(:,4))- minima_pos_y(TRI_min(:,3))).^2 + (minima_pos_z(TRI_min(:,4))- minima_pos_z(TRI_min(:,3))).^2 );
    0421
    0422 %Calculating nearest neighbours for each minima point
    0423 %Comparing tetrahedron edges associated with each vertex
    0424 emn1 = min([e1, e2, e4],[],2); %Comparing edges 1, 2 and 4 (vertex 1)
    0425 emn2 = min([e1, e3, e5],[],2); %Comparing edges 1, 3 and 5 (vertex 2)
    0426 emn3 = min([e2, e3, e6],[],2); %Comparing edges 2, 3 and 6 (vertex 3)
    0427 emn4 = min([e4, e5, e6],[],2); %Comparing edges 4, 5 and 6 (vertex 4)
    0428
    0429 e = [emn1 ,emn2, emn3, emn4];
    0430
    0431 %Making sure that the smallest edge associated with the each vertex is stored
    0432 %correctly
    0433 for i=1:length(emn1)
    0434     for j=1:4
    0435         if min_nearest(TRI_min(i,j)) > e(i,j) || min_nearest(TRI_min(i,j)) == 0
    0436             min_nearest(TRI_min(i,j)) = e(i,j);
    0437         end
    0438     end
    0439 end
    0440
    0441 %Window size calculations
    0442
    0443 d1 = min( min(max_nearest) , min(min_nearest) );
    0444 d2 = max( min(max_nearest) , min(min_nearest) );
    0445 d3 = min( max(max_nearest) , max(min_nearest) );
    0446 d4 = max( max(max_nearest) , max(min_nearest) );
    0447 d5 = (d1+d2+d3+d4)/4 ;
    0448 d6 = median([min_nearest; max_nearest]);
    0449 d7 = mode([min_nearest; max_nearest]);
    0450
    0451 Windows = [d1, d2, d3, d4, d5, d6, d7];
    0452
    0453 %making sure w_size is an odd integer
    0454 Windows = 2*(floor(Windows./2))+1;
    0455
    0456 if(Windows(type)<3)
    0457     warning('WARNING: Calculated Window size less than 3');
    0458     warning('Overriding calculated value and setting window size = 3');
    0459     Windows(type) = 3;
    0460 end
    0461 end
    0462
    0463 function H1 = Sift(H,w_sz)
    0464
    0465 %Envelope Generation
    0466 [Env_max,Env_min] = OSF(H,w_sz);
    0467
    0468 %padding
    0469 Env_med = Pad_smooth(Env_max,Env_min,w_sz);
    0470
    0471 %Subtracting from residue
    0472 H1 = H - Env_med;
    0473
    0474 end
    0475
    0476 function [Max,Min] = OSF(H,w_sz)
    0477 %Order statistics filtering to determine maximum and minmum envelopes
    0478             Max = Separable_ordfilt3(H, 'max', w_sz); %Max envelope
    0479             Min = Separable_ordfilt3(H, 'min', w_sz); %Min envelope
    0480
    0481             function Signal = Separable_ordfilt3(Signal, order, w_sz)
    0482                 % Purpose:
    0483                 % -To perform separable order statistics filtering of 3D
    0484                 % signals
    0485                 % -Boundary condition is always symmetric
    0486
    0487                 [X,Y,Z] = size(Signal);
    0488
    0489                 %Separable Filtering
    0490                 %First Dimension (X)
    0491                 for k = 1:Z
    0492                     for j = 1:Y
    0493                         Signal(:,j,k) = Ordfilt1(Signal(:,j,k),order,w_sz(1));
    0494                     end
    0495                 end
    0496
    0497                 %Second Dimension (Y)
    0498                 for k = 1:Z
    0499                     for i = 1:X
    0500                         Signal(i,:,k) = Ordfilt1(Signal(i,:,k),order,w_sz(2));
    0501                     end
    0502                 end
    0503
    0504                 %Third Dimension (Z)
    0505                 for j = 1:Y
    0506                     for i = 1:X
    0507                         Signal(i,j,:) = Ordfilt1(Signal(i,j,:),order,w_sz(3));
    0508                     end
    0509                 end
    0510
    0511                 function f_signal = Ordfilt1(signal,order,window_size)
    0512
    0513                     %1-D Rank order filter function
    0514
    0515                     %Pre-processing
    0516                     [a,b,c] = size(signal);           %Original signal size
    0517                     signal  = squeeze(signal);        %Removing the singleton dimensions
    0518                     L       = length(signal);         %Length of the signal
    0519                     signal  = reshape(signal, [L,1]); %Ensure that the processed signal is always a column vector
    0520
    0521                     r = (window_size-1)/2;
    0522
    0523                     %Padding boundaries
    0524                     x = [flip(signal(1:r)); signal ;flip(signal(end-(r-1):end))];
    0525
    0526                     [M,~] = size(x);
    0527                     y = zeros(size(x));
    0528
    0529                     switch order
    0530                         case 'max'
    0531                             for m = 1+r:M-r
    0532                                 % Extract a window of size (2r+1) around (m)
    0533                                 temp = x((m-r):(m+r));
    0534                                 w = sort(temp);
    0535                                 y(m) = w(end); % Select the greatest element
    0536                             end
    0537                         case 'min'
    0538                             for m = 1+r:M-r
    0539                                 % Extract a window of size (2r+1) around (m)
    0540                                 temp = x((m-r):(m+r));
    0541                                 w = sort(temp);
    0542                                 y(m) = w(1); % Select the smallest element
    0543                             end
    0544                         otherwise
    0545                             error('No such filering operation defined')
    0546                     end
    0547
    0548                     f_signal = y(1+r:end-r);
    0549
    0550                     f_signal = reshape(f_signal,[a,b,c]); %Restoring Signal size
    0551                 end
    0552             end
    0553 end
    0554
    0555 function Env_med = Pad_smooth(Env_max,Env_min,w_sz)
    0556 hx = floor(w_sz(1)/2);
    0557 hy = floor(w_sz(2)/2);
    0558 hz = floor(w_sz(3)/2);
    0559 %Padding
    0560 temp = padarray(Env_max,[hx hy],'replicate');
    0561 temp1 = permute(temp,[3 2 1]); %interchanging dimensions
    0562 temp = padarray(temp1,[hz 0],'replicate');
    0563 Env_maxp = permute(temp,[3 2 1]); %restoring dimensions
    0564
    0565 temp = padarray(Env_min,[hx hy],'replicate');
    0566 temp1 = permute(temp,[3 2 1]); %interchanging dimensions
    0567 temp = padarray(temp1,[hz 0],'replicate');
    0568 Env_minp = permute(temp,[3 2 1]); %restoring dimensions
    0569
    0570 %Smoothing
    0571
    0572 temp1 = movmean(Env_maxp,w_sz(3),3,'endpoints','discard');
    0573 temp2 = movmean(temp1,w_sz(2),2,'endpoints','discard');
    0574 Env_maxs = movmean(temp2,w_sz(1),1,'endpoints','discard');
    0575
    0576 temp1 = movmean(Env_minp,w_sz(3),3,'endpoints','discard');
    0577 temp2 = movmean(temp1,w_sz(2),2,'endpoints','discard');
    0578 Env_mins = movmean(temp2,w_sz(1),1,'endpoints','discard');
    0579
    0580 %Calculating mean envelope
    0581 Env_med = (Env_maxs + Env_mins)./2;
    0582
    0583 end
    0584
    0585 function [IO,Error] = Orth_index(Signal,IMF,Residue)
    0586 % Purpose:
    0587 % To calculate the index of orthogonality of a decomposition and its mean
    0588 % squared error
    0589
    0590 n_imf = size(IMF,4);
    0591 numerator = zeros(size(Signal));
    0592 I = sum(IMF,4) + Residue;
    0593
    0594 Error.map = (Signal-I)./Signal;
    0595 Error.global = immse(I,Signal);
    0596
    0597 for j = 1:n_imf
    0598     for k = 1:n_imf
    0599         if(j~=k)
    0600            numerator = numerator + IMF(:,:,:,j).*IMF(:,:,:,k);
    0601         end
    0602     end
    0603 end
    0604
    0605 IO.map = numerator/sum(sum(sum(Signal.^2))); %wrong
    0606 IO.global = sum(sum(sum(IO.map)));
    0607 end
    0608
    0609 function Plot_results(u,v,w,Results,param)
    0610 % default plot attributes
    0611 set(groot,'defaultaxesfontname','times');
    0612 set(groot,'defaultaxesfontsize',12);
    0613 set(groot,'defaulttextInterpreter','latex');
    0614 set(groot,'defaultLineLineWidth',2);
    0615
    0616 Colour = parula;
    0617 nslice  = param.nslice;
    0618
    0619 figure(1)
    0620         subplot(1,3,1)
    0621         TIMF_plot(u,Colour,nslice,0,'Signal','u');
    0622         subplot(1,3,2)
    0623         TIMF_plot(v,Colour,nslice,0,'Signal','v');
    0624         subplot(1,3,3)
    0625         TIMF_plot(w,Colour,nslice,0,'Signal','w');
    0626
    0627
    0628     for i=1:param.nimfs
    0629      figure(i+1)
    0630         subplot(1,3,1)
    0631         TIMF_plot(Results.IMF.u(:,:,:,i),Colour,nslice,i,'IMF','u');
    0632         subplot(1,3,2)
    0633         TIMF_plot(Results.IMF.v(:,:,:,i),Colour,nslice,i,'IMF','v');
    0634         subplot(1,3,3)
    0635         TIMF_plot(Results.IMF.w(:,:,:,i),Colour,nslice,i,'IMF','w');
    0636     end
    0637
    0638     figure(i+2)
    0639     subplot(1,3,1)
    0640         TIMF_plot(Results.Residue.u,Colour,nslice,0,'Residue','u');
    0641         subplot(1,3,2)
    0642         TIMF_plot(Results.Residue.v,Colour,nslice,0,'Residue','v');
    0643         subplot(1,3,3)
    0644         TIMF_plot(Results.Residue.w,Colour,nslice,0,'Residue','w');
    0645 end
    0646
    0647 function TIMF_plot(signal,Colour,nslice,imf,name1,name2)
    0648
    0649     [Nx, Ny, Nz] = size(signal);
    0650
    0651     xslice = linspace(1,Nx,nslice);
    0652     yslice = linspace(1,Ny,nslice);
    0653     zslice = linspace(1,Nz,nslice);
    0654     volume = slice(signal,xslice,yslice,zslice);
    0655     axis equal;
    0656     xlabel('x');
    0657     ylabel('y');
    0658     zlabel('z');
    0659     set(gca,'TickLabelInterpreter','latex')
    0660     switch(name1)
    0661         case 'IMF'
    0662             title(sprintf('%s %d %s',name1,imf,name2));
    0663         case 'Signal'
    0664             title(sprintf('%s %s',name1,name2));
    0665         case 'Residue'
    0666             title(sprintf('%s %s',name1,name2));
    0667     end
    0668     colorbar;
    0669     set(volume,'EdgeColor','none',...
    0670         'FaceColor','interp',...
    0671         'FaceAlpha','interp')
    0672     alpha('color')
    0673     view(30,30);
    0674     alphamap('rampup')
    0675     alphamap('decrease',.1)
    0676     colormap(Colour);
    0677 %     caxis([-3 3]);
    0678     hcb = colorbar;
    0679     colorTitleHandle = get(hcb,'Title');
    0680     titleString = '$\frac{u}{U_{\infty}}$';
    0681     set(colorTitleHandle ,'String',titleString,'Interpreter','latex','FontSize',14);
    0682     set(hcb,'TickLabelInterpreter','latex');
    0683
    0684 end