Description of EMD2D2V

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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
