Description of EMD2D3V

0001 % Purpose:
0002 % -To perform EMD on 3 channels of 2 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, 'off' to hide IMF plots
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 %   [1] Bhuiyan et. al, 'Fast and Adaptive Bidimensional EmpiricalMode
0025 %       Decomposition Using Order-Statistics Filter Based
0026 %       Envelope Estimation',2008
0027 %
0028 %   [2] FABEEMD (Matthew Koll, Dept. of Aerospace Engineering, University
0029 %                of Illinois at Urbana-Champaign)
0030 %
0031 %
0032 % Written by Mruthun Thirumalaisamy
0033 % Graduate Student
0034 % Department of Aerospace Engineering
0035 % University of Illinois at Urbana-Champaign
0036 % May 11 2018
0037 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0038 
0039 function Results = EMD2D3V(u,v,w,param) 
0040 
0041 %Reading signal characteristics
0042 [Nx,Ny] = size(u); %Signal dimensions
0043 B       = size(v); %Signal dimensions
0044 C       = size(w); %Signal dimensions
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]==B) ||  ~all([Nx,Ny]==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 %Initialisations
0087 IMF.u = zeros(Nx, Ny ,param.nimfs); 
0088 IMF.v = zeros(Nx, Ny ,param.nimfs);
0089 IMF.w = zeros(Nx, Ny ,param.nimfs);
0090 Residue.u = u; Residue.v = v; Residue.w = w;
0091 
0092 Windows = zeros(7,param.nimfs);
0093 
0094 sift_cnt = zeros(1,param.nimfs);
0095 imf = 1;
0096 stopflag = 1;
0097 
0098     while(imf <= param.nimfs && stopflag)
0099         %Initialising intermediary IMFs
0100         H.u = Residue.u; H.v = Residue.v; H.w = Residue.w;
0101 
0102         sift_stop = 0; %flag to control sifting loop
0103         
0104         Combined = H.u/sqrt(3) + H.v/sqrt(3) + H.w/sqrt(3); %Combining three signals with equal weights
0105         [maxima,minima] = Identify_max_min(Combined);  %Obtaining extrema of combined signal
0106         
0107         %Checking whether there are too few extrema in the IMF
0108         if (nnz(maxima) < 3 || nnz(minima) < 3)
0109             warning('Fewer than three extrema found in maxima map. Stopping now...');
0110             break;
0111         end
0112         
0113         %Window size determination by delaunay triangulation
0114         Windows(:,imf) = filter_size(maxima,minima,param.type);        
0115         w_sz = Windows(param.type,imf); %extracting window size chosen by input parameter
0116         
0117         if~(w_sz)
0118            warning('EMD2D3V has stopped because the Delaunay Triangulation could not be created (collinear points)'); 
0119            stopflag = 0; %#ok<NASGU>
0120            break;
0121         end
0122         
0123         %Begin sifting iteration
0124         while~(sift_stop)            
0125             sift_cnt(imf) = sift_cnt(imf) + 1; %Incrementing sift counter
0126             %Envelope Generation
0127             Env = OSF(H,w_sz);
0128             
0129             %padding
0130             Env = Pad_smooth(Env,w_sz);
0131            
0132             %Calculating mean envelope
0133             Env.u.med = (Env.u.maxs + Env.u.mins)./2;
0134             Env.v.med = (Env.v.maxs + Env.v.mins)./2;
0135             Env.w.med = (Env.w.maxs + Env.w.mins)./2;
0136             
0137             %Subtracting from residue
0138             H1.u = H.u - Env.u.med; H1.v = H.v - Env.v.med; H1.w = H.w - Env.w.med;         
0139              
0140             %Stop condition checks
0141             mse_u = immse(H1.u,H.u); mse_v = immse(H1.v,H.v); mse_w = immse(H1.w,H.w);         
0142             if (mse_u<param.tol && mse_v<param.tol && mse_w<param.tol && sift_cnt(imf)~=1)
0143                 sift_stop = 1;
0144             end
0145             
0146             H.u = H1.u; H.v = H1.v; H.w = H1.w;                
0147         end
0148         
0149         %Storing IMFs
0150         IMF.u(:,:,imf) = H.u; IMF.v(:,:,imf) = H.v; IMF.w(:,:,imf) = H.w;
0151 
0152         %Subtracting from Residual Signals
0153         Residue.u = Residue.u - IMF.u(:,:,imf);
0154         Residue.v = Residue.v - IMF.v(:,:,imf);
0155         Residue.w = Residue.w - IMF.w(:,:,imf);  
0156         
0157         %Incrementing IMF counter
0158         imf = imf + 1;
0159         
0160     end
0161     
0162     %Checking for oversifting
0163     if(any(sift_cnt>=5*ones(size(sift_cnt))))
0164         warning('Decomposition may be oversifted. Checking if window size increases monotonically...');
0165         
0166         if( any (diff(Windows(param.type,:)) <= zeros(1,param.nimfs-1)) )
0167         warning('Filter window size does not increase monotonically')
0168         end
0169     end
0170     
0171     %Organising results
0172     Results.IMF = IMF;
0173     Results.windowtype = param.type;
0174     Results.Residue = Residue;
0175     Results.Windows = Windows;
0176     Results.Sifts = sift_cnt;
0177     
0178     %Error and orthogonality
0179     [Results.IO.u,Results.Error.u] = Orth_index(u,IMF.u,Residue.u);
0180     [Results.IO.v,Results.Error.v] = Orth_index(v,IMF.v,Residue.v);
0181     [Results.IO.w,Results.Error.w] = Orth_index(w,IMF.w,Residue.w);
0182     
0183     switch(param.plot)
0184         case 'on'
0185             Plot_results(u,v,w,Results,param)
0186     end
0187 end
0188 
0189 function [maxima,minima] = Identify_max_min(signal)
0190 % Purpose:
0191 % To identify the maxima and minima locations in a two or three dimensional array
0192 mask = ones(3); mask(5) = 0; %Window size for the extrema detection fixed at 3x3 (Bhuiyan et.al)
0193 
0194 B = ordfilt2(signal,8,mask);
0195 C = ordfilt2(signal,1,mask);
0196 maxima = signal >= B;
0197 minima = signal <= C;
0198 end
0199 
0200 function Windows = filter_size(maxima_map, minima_map,type)
0201 % Purpose:
0202 % -To determine the window size for order statistics filtering of a signal.
0203 % The determination of the window size is based on the work of Bhuiyan et
0204 % al.
0205 %
0206 % Inputs:
0207 % -Two 2D extrema maps
0208 %
0209 % Outputs:
0210 % -Calculated value of the window size
0211 %
0212 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0213 
0214 %use delaunay triangulation to determine the nearest neighbours and hence
0215 %the filter size
0216 
0217 %processing d_max
0218 [maxima_pos_y,maxima_pos_x] = find(maxima_map);
0219 
0220 max_nearest = zeros(length(maxima_pos_y),1);
0221 
0222 try
0223 TRI_max = delaunay([maxima_pos_x maxima_pos_y]);
0224 catch
0225     warning('Maxima points are collinear. Exiting without further iterations');
0226     Windows = [0, 0, 0, 0, 0, 0, 0];
0227     return
0228 end
0229 
0230 %Calculating 3 edge lengths for each triangle
0231 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 );
0232 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 );
0233 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 );
0234 
0235 %Calculating nearest neighbours for each maxima point
0236 %Comparing edges associated with each vertex
0237 em1 = min([e2, e1],[],2); %Comparing edges 2 and 1 (vertex 1)
0238 em2 = min([e3, e1],[],2); %Comparing edges 3 and 1 (vertex 2)
0239 em3 = min([e3, e2],[],2); %Comparing edges 3 and 2 (vertex 3)
0240 
0241 e = [em1 ,em2, em3]; %Smaller edge for each vertex in each triangle (since one vertex is associated with two edges in a triangle)
0242 
0243 %Making sure that the smallest edge associated with the each vertex is stored
0244 %correctly
0245 for i=1:length(em1)
0246     for j=1:3
0247         if max_nearest(TRI_max(i,j)) > e(i,j) || max_nearest(TRI_max(i,j)) == 0
0248             max_nearest(TRI_max(i,j)) = e(i,j);
0249         end
0250     end
0251 end
0252 
0253 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0254 %processing d_min
0255 
0256 [minima_pos_y,minima_pos_x] = find(minima_map);
0257 min_nearest = zeros(length(minima_pos_y),1);
0258 
0259 try
0260 TRI_min = delaunay([minima_pos_x minima_pos_y]);
0261 catch
0262     warning('Minima points are collinear. Exiting without further iterations');
0263     Windows = [0, 0, 0, 0, 0, 0, 0];
0264     return
0265 end
0266 
0267 %Calculating 3 neighbour distances for each minima point
0268 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 );
0269 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 );
0270 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 );
0271 
0272 %Calculating nearest neighbours for each maxima point
0273 
0274 %Comparing triangle edges associated with each vertex
0275 emn1 = min([e2, e1],[],2); %Comparing edges 2 and 1 (vertex 1)
0276 emn2 = min([e3, e1],[],2); %Comparing edges 3 and 1 (vertex 2)
0277 emn3 = min([e3, e2],[],2); %Comparing edges 3 and 2 (vertex 3)
0278 
0279 e = [emn1 ,emn2, emn3]; %Smaller edge for each vertex in each triangle (since one vertex is associated with two edges in a triangle)
0280 
0281 %Making sure that the smallest edge associated with the each vertex is stored
0282 %correctly
0283 for i=1:length(emn1)
0284     for j=1:3
0285         if min_nearest(TRI_min(i,j)) > e(i,j) || min_nearest(TRI_min(i,j)) == 0
0286             min_nearest(TRI_min(i,j)) = e(i,j);
0287         end
0288     end
0289 end
0290 
0291 %Window size calculations
0292 
0293 d1 = min( min(max_nearest) , min(min_nearest) );
0294 d2 = max( min(max_nearest) , min(min_nearest) );
0295 d3 = min( max(max_nearest) , max(min_nearest) );
0296 d4 = max( max(max_nearest) , max(min_nearest) );
0297 d5 = (d1+d2+d3+d4)/4 ;
0298 d6 = median([min_nearest; max_nearest]);
0299 d7 = mode([min_nearest; max_nearest]);
0300 
0301 Windows = [d1, d2, d3, d4, d5, d6, d7];
0302 
0303 %making sure w_size is an odd integer
0304 Windows = 2*(floor(Windows./2))+1;
0305          
0306 if(Windows(type)<3)
0307     warning('WARNING: Calculated Window size less than 3');
0308     warning('Overriding calculated value and setting window size = 3');
0309     Windows(type) = 3;
0310 end
0311 end
0312 
0313 function Env = OSF(H,w_sz)
0314 %Order statistics filtering to determine maximum and minmum envelopes
0315             Env.u.max = ordfilt2(H.u ,w_sz.^2, true(w_sz),'symmetric'); %Max envelope u
0316             Env.u.min = ordfilt2(H.u ,1, true(w_sz),'symmetric');       %Min envelope u
0317             
0318             Env.v.max = ordfilt2(H.v ,w_sz.^2, true(w_sz),'symmetric'); %Max envelope v
0319             Env.v.min = ordfilt2(H.v ,1, true(w_sz),'symmetric');       %Min envelope v
0320             
0321             Env.w.max = ordfilt2(H.w ,w_sz.^2, true(w_sz),'symmetric'); %Max envelope w
0322             Env.w.min = ordfilt2(H.w ,1, true(w_sz),'symmetric');       %Min envelope w
0323 end
0324 
0325 function Env = Pad_smooth(Env,w_sz)
0326 h = floor(w_sz/2);
0327 
0328 %Padding
0329 %u
0330 Env.u.maxp = padarray(Env.u.max,[h h],'replicate');
0331 Env.u.minp = padarray(Env.u.min,[h h],'replicate');
0332 %v
0333 Env.v.maxp = padarray(Env.v.max,[h h],'replicate');
0334 Env.v.minp = padarray(Env.v.min,[h h],'replicate');
0335 %w
0336 Env.w.maxp = padarray(Env.w.max,[h h],'replicate');
0337 Env.w.minp = padarray(Env.w.min,[h h],'replicate');
0338 
0339 %Smoothing
0340 %u
0341 temp = movmean(Env.u.maxp,w_sz,2,'endpoints','discard');
0342 Env.u.maxs = movmean(temp,w_sz,1,'endpoints','discard');
0343 temp = movmean(Env.u.minp,w_sz,2,'endpoints','discard');
0344 Env.u.mins = movmean(temp,w_sz,1,'endpoints','discard');
0345 %v
0346 temp = movmean(Env.v.maxp,w_sz,2,'endpoints','discard');
0347 Env.v.maxs = movmean(temp,w_sz,1,'endpoints','discard');
0348 temp = movmean(Env.v.minp,w_sz,2,'endpoints','discard');
0349 Env.v.mins = movmean(temp,w_sz,1,'endpoints','discard');
0350 %w
0351 temp = movmean(Env.w.maxp,w_sz,2,'endpoints','discard');
0352 Env.w.maxs = movmean(temp,w_sz,1,'endpoints','discard');
0353 temp = movmean(Env.w.minp,w_sz,2,'endpoints','discard');
0354 Env.w.mins = movmean(temp,w_sz,1,'endpoints','discard');
0355 
0356 end
0357 
0358 function [IO,Error] = Orth_index(Signal,IMF,Residue)
0359 % Purpose:
0360 % To calculate the index of orthogonality of a decomposition and its mean
0361 % squared error
0362 
0363 n_imf = size(IMF,3);
0364 numerator = zeros(size(Signal));
0365 I = sum(IMF,3) + Residue;
0366 
0367 Error.map = (Signal-I)./Signal;
0368 Error.global = immse(I,Signal);
0369 
0370 for j = 1:n_imf
0371     for k = 1:n_imf
0372         if(j~=k)
0373            numerator = numerator + IMF(:,:,j).*IMF(:,:,k);
0374         end
0375     end
0376 end
0377 
0378 IO.map = numerator/sum(sum(Signal.^2));
0379 IO.global = sum(sum(IO.map));
0380 end
0381 
0382 function Plot_results(u,v,w,Results,param)
0383 % default plot attributes
0384 set(groot,'defaultaxesfontname','times');
0385 set(groot,'defaultaxesfontsize',12);
0386 set(groot,'defaulttextInterpreter','latex');
0387 set(groot,'defaultLineLineWidth',2);
0388 
0389 Colour = redblue;
0390 
0391 figure(1)   
0392         subplot(3,1,1)
0393         BIMF_plot(u,Colour,0,'Signal','u');
0394         subplot(3,1,2)
0395         BIMF_plot(v,Colour,0,'Signal','v');
0396         subplot(3,1,3)
0397         BIMF_plot(w,Colour,0,'Signal','w');
0398 
0399 
0400     for i=1:param.nimfs
0401      figure(i+1)   
0402         subplot(3,1,1)
0403         BIMF_plot(Results.IMF.u(:,:,i),Colour,i,'IMF','u');
0404         subplot(3,1,2)
0405         BIMF_plot(Results.IMF.v(:,:,i),Colour,i,'IMF','v');
0406         subplot(3,1,3)
0407         BIMF_plot(Results.IMF.w(:,:,i),Colour,i,'IMF','w');
0408     end
0409     
0410     figure(i+2)
0411     subplot(3,1,1)
0412         BIMF_plot(Results.Residue.u,Colour,0,'Residue','u');
0413         subplot(3,1,2)
0414         BIMF_plot(Results.Residue.v,Colour,0,'Residue','v');
0415         subplot(3,1,3)
0416         BIMF_plot(Results.Residue.w,Colour,0,'Residue','w');
0417 end
0418 
0419 function BIMF_plot(signal,Colour,imf,name1,name2)
0420 %Masking wall data
0421 % load('MASK_file','MASK');
0422 % signal = MASK.*signal;
0423 
0424     imAlpha=ones(size(signal));
0425     imAlpha(isnan(signal))=0;    
0426     imagesc(signal,'AlphaData',imAlpha);
0427     set(gca,'color',0*[1 1 1]);
0428     xlabel('$x$')
0429     ylabel('$y$')
0430     axis equal;
0431     axis tight;
0432 
0433     switch(name1)
0434         case 'IMF'
0435             title(sprintf('%s %d %s ',name1,imf,name2));
0436         case 'Residue'
0437             title(sprintf('%s %s ',name1,name2));
0438         case 'Signal'
0439             title(sprintf('%s %s ',name1,name2));
0440     end
0441     set(gca,'TickLabelInterpreter','latex')
0442     colormap(Colour);
0443     hcb = colorbar;
0444     colorTitleHandle = get(hcb,'Title');
0445     titleString = 'u';
0446     set(colorTitleHandle ,'String',titleString,'Interpreter','latex','FontSize',14);
0447     set(hcb,'TickLabelInterpreter','latex');
0448 end
EMD2D3V

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
