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: This function is called by:

SUBFUNCTIONS ^

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
Generated on Thu 18-Apr-2019 12:22:00 by m2html © 2005