EMD3D3V_parallel

 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

 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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