Description of EMD3D3V

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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
