EMD3D2V_parallel

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

 Input: 
 - u: Signal 1
 - v: Signal 2
 - 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 (Modified: Dec 14 2018)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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