EMD1D3V

 Purpose: 
 -To perform EMD on 3 channels of 1 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

 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 1 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 %
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
0031 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0032 
0033 function Results = EMD1D3V(u,v,w,t,param) 
0034 
0035 %Reading signal characteristics
0036 [Nx] = length(u); %Signal dimensions
0037 B    = length(v); %Signal dimensions
0038 C    = length(w); %Signal dimensions
0039 
0040 %Preliminary checks
0041 if ~isfield(param,'nimfs')
0042     param.nimfs = 10;
0043 end
0044 
0045 if ~isfield(param,'tol')
0046     param.tol = 0.05; % 0.1% of the minimum signal amplitude
0047 end
0048 
0049 if ~isfield(param,'type')
0050     param.type = 6;
0051 end
0052 
0053 if ~isfield(param,'plot')
0054     param.plot = 'off';
0055 end
0056 
0057 if(~all(ismember(param.type,[1,2,3,4,5,6,7])))
0058     error('Please enter a valid window size type')
0059 end
0060 
0061 if(~all(Nx==B) ||  ~all(Nx==C))
0062     error('Inconsistent dimensions between channels. Please check input data');
0063 end
0064 clearvars B C
0065 
0066 if(param.tol<=0.005)
0067    warning('Low sifting tolerance may cause oversifting');
0068    answer = questdlg('Would you like to continue?', ...
0069     'User set low sifting tolerance', ...
0070     'Yes','No','No');
0071     % Handle response
0072     switch answer
0073         case 'Yes'
0074             
0075         case 'No'
0076             return;
0077     end
0078 end
0079 
0080 if ~isfield(param,'plot')
0081     param.plot = 'off';
0082 end
0083 
0084 %Initialisations
0085 IMF.u = zeros(Nx, param.nimfs); 
0086 IMF.v = zeros(Nx, param.nimfs);
0087 IMF.w = zeros(Nx, param.nimfs);
0088 Residue.u = u; Residue.v = v; Residue.w = w;
0089 
0090 Windows = zeros(7,param.nimfs);
0091 
0092 sift_cnt = zeros(1,param.nimfs);
0093 imf = 1;
0094 
0095     while(imf <= param.nimfs)
0096         %Initialising intermediary IMFs
0097         H.u = Residue.u; H.v = Residue.v; H.w = Residue.w;
0098 
0099         sift_stop = 0; %flag to control sifting loop
0100         
0101         Combined = H.u/sqrt(3) + H.v/sqrt(3) + H.w/sqrt(3) ; %Combining two signals with equal weights
0102         [Maxima,MaxPos,Minima,MinPos] = extrema(Combined);  %Obtaining extrema of combined signal
0103         
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 extrema map. Stopping now...');
0108             break;
0109         end
0110         
0111         %Window size determination by delaunay triangulation
0112         Windows(:,imf) = filter_size1D(MaxPos,MinPos,param.type);        
0113         w_sz = Windows(param.type,imf); %extracting window size chosen by input parameter
0114         
0115         %Begin sifting iteration
0116         while~(sift_stop)            
0117             sift_cnt(imf) = sift_cnt(imf) + 1; %Incrementing sift counter
0118             %Envelope Generation
0119             Env = OSF(H,w_sz);
0120             
0121             %padding
0122             Env = Pad_smooth(Env,w_sz);
0123            
0124             %Calculating mean envelope
0125             Env.u.med = (Env.u.maxs + Env.u.mins)./2;
0126             Env.v.med = (Env.v.maxs + Env.v.mins)./2;
0127             Env.w.med = (Env.w.maxs + Env.w.mins)./2;
0128             
0129 %             osfplot(t, H, Env);
0130             
0131             %Subtracting from residue
0132             H1.u = H.u - Env.u.med; H1.v = H.v - Env.v.med; H1.w = H.w - Env.w.med;     
0133             
0134 %             projenv_IMF_plot(t, Combined, MaxPos,MinPos)
0135              
0136             %Stop condition checks
0137             mse_u = immse(H1.u,H.u); mse_v = immse(H1.v,H.v); mse_w = immse(H1.w,H.w);       
0138             if (mse_u<param.tol && mse_v<param.tol && mse_w<param.tol && sift_cnt(imf)~=1)
0139                 sift_stop = 1;
0140             end
0141             
0142             H.u = H1.u; H.v = H1.v; H.w = H1.w;           
0143         end
0144         
0145         %Storing IMFs
0146         IMF.u(:,imf) = H.u; IMF.v(:,imf) = H.v; IMF.w(:,imf) = H.w;
0147 
0148         %Subtracting from Residual Signals
0149         Residue.u = Residue.u - IMF.u(:,imf);
0150         Residue.v = Residue.v - IMF.v(:,imf);
0151         Residue.w = Residue.w - IMF.w(:,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(1,param.nimfs)))
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     [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,t,Results,param)
0181     end
0182 end
0183 
0184 function Windows = filter_size1D(imax, imin, type)
0185 %
0186 % Purpose:
0187 % -To determine the window size for order statistics filtering of a signal.
0188 % The determination of the window size is based on the work of Bhuiyan et
0189 % al.
0190 %
0191 % Inputs:
0192 % -Two 1D extrema maps
0193 %
0194 % Outputs:
0195 % -Calculated value of the window size
0196 %
0197 % Written by Mruthun Thirumalaisamy
0198 % Graduate Student
0199 % Department of Aerospace Engineering
0200 % University of Illinois at Urbana-Champaign
0201 % March 30 2018
0202 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0203 
0204 edge_len_max = diff(sort(imax));
0205 edge_len_min = diff(sort(imin));
0206     
0207     
0208         %Window size calculations
0209         
0210         d1 = min( min(edge_len_max) , min(edge_len_min) );
0211         d2 = max( min(edge_len_max) , min(edge_len_min) );
0212         d3 = min( max(edge_len_max) , max(edge_len_min) );
0213         d4 = max( max(edge_len_max) , max(edge_len_min) );
0214         d5 = (d1+d2+d3+d4)/4 ;
0215         d6 = median([edge_len_max; edge_len_min]);
0216         d7 = mode([edge_len_max; edge_len_min]);
0217         
0218         Windows = [d1, d2, d3, d4, d5, d6, d7];
0219 
0220 %making sure w_size is an odd integer
0221 Windows = 2*(floor(Windows./2))+1;
0222          
0223 if(Windows(type)<3)
0224     warning('WARNING: Calculated Window size less than 3');
0225     warning('Overriding calculated value and setting window size = 3');
0226     Windows(type) = 3;
0227 end
0228 
0229 end
0230 
0231 function Env = OSF(H,w_sz)
0232 %Order statistics filtering to determine maximum and minmum envelopes
0233             Env.u.max = Ordfilt1(H.u, 'max', w_sz); %Max envelope u
0234             Env.u.min = Ordfilt1(H.u, 'min', w_sz); %Min envelope u
0235              
0236             Env.v.max = Ordfilt1(H.v, 'max', w_sz); %Max envelope v
0237             Env.v.min = Ordfilt1(H.v, 'min', w_sz); %Min envelope v
0238             
0239             Env.w.max = Ordfilt1(H.w, 'max', w_sz); %Max envelope w
0240             Env.w.min = Ordfilt1(H.w, 'min', w_sz); %Min envelope w
0241                 
0242         function f_signal = Ordfilt1(signal,order,window_size)
0243 
0244             %1-D Rank order filter function
0245 
0246             %Pre-processing
0247             [a,b,c] = size(signal);           %Original signal size
0248             signal  = squeeze(signal);        %Removing the singleton dimensions
0249             L       = length(signal);         %Length of the signal
0250             signal  = reshape(signal, [L,1]); %Ensure that the processed signal is always a column vector
0251 
0252             r = (window_size-1)/2;
0253 
0254             %Padding boundaries
0255             x = [flip(signal(1:r)); signal ;flip(signal(end-(r-1):end))];
0256 
0257             [M,~] = size(x);
0258             y = zeros(size(x));
0259 
0260             switch order
0261                 case 'max'
0262                     for m = 1+r:M-r
0263                         % Extract a window of size (2r+1) around (m)
0264                         temp = x((m-r):(m+r));
0265                         w = sort(temp);
0266                         y(m) = w(end); % Select the greatest element
0267                     end
0268                 case 'min'
0269                     for m = 1+r:M-r
0270                         % Extract a window of size (2r+1) around (m)
0271                         temp = x((m-r):(m+r));
0272                         w = sort(temp);
0273                         y(m) = w(1); % Select the smallest element
0274                     end
0275                 otherwise
0276                     error('No such filering operation defined')
0277             end
0278 
0279             f_signal = y(1+r:end-r);
0280 
0281             f_signal = reshape(f_signal,[a,b,c]); %Restoring Signal size
0282         end
0283       
0284 end
0285 
0286 function Env = Pad_smooth(Env,w_sz)
0287 h = floor(w_sz/2);
0288 
0289 %Padding
0290 %u
0291 Env.u.maxp = padarray(Env.u.max,[h 0],'symmetric');
0292 Env.u.minp = padarray(Env.u.min,[h 0],'symmetric');
0293 %v
0294 Env.v.maxp = padarray(Env.v.max,[h 0],'symmetric');
0295 Env.v.minp = padarray(Env.v.min,[h 0],'symmetric');
0296 %w
0297 Env.w.maxp = padarray(Env.w.max,[h 0],'symmetric');
0298 Env.w.minp = padarray(Env.w.min,[h 0],'symmetric');
0299 
0300 %Smoothing
0301 %u
0302 Env.u.maxs = movmean(Env.u.maxp,w_sz,1,'endpoints','discard');
0303 Env.u.mins = movmean(Env.u.minp,w_sz,1,'endpoints','discard');
0304 %v
0305 Env.v.maxs = movmean(Env.v.maxp,w_sz,1,'endpoints','discard');
0306 Env.v.mins = movmean(Env.v.minp,w_sz,1,'endpoints','discard');
0307 %w
0308 Env.w.maxs = movmean(Env.w.maxp,w_sz,1,'endpoints','discard');
0309 Env.w.mins = movmean(Env.w.minp,w_sz,1,'endpoints','discard');
0310 
0311 end
0312 
0313 function [IO,Error] = Orth_index(Signal,IMF,Residue)
0314 % Purpose:
0315 % To calculate the index of orthogonality of a decomposition and its mean
0316 % squared error
0317 
0318 n_imf = size(IMF,2);
0319 numerator = zeros(size(Signal));
0320 I = sum(IMF,2) + Residue;
0321 
0322 Error.map = (Signal-I)./Signal;
0323 Error.global = immse(I,Signal);
0324 
0325 for j = 1:n_imf
0326     for k = 1:n_imf
0327         if(j~=k)
0328            numerator = numerator + IMF(:,j).*IMF(:,k);
0329         end
0330     end
0331 end
0332 
0333 IO.map = numerator/(sum(Signal.^2)); %wrong
0334 IO.global = sum(IO.map);
0335 end
0336 
0337 function Plot_results(u,v,w,t,Results,param)
0338 % default plot attributes
0339 set(groot,'defaultaxesfontname','times');
0340 set(groot,'defaultaxesfontsize',12);
0341 set(groot,'defaulttextInterpreter','tex');
0342 set(groot,'defaultLineLineWidth',2);
0343 
0344 figure(1)   
0345         subplot(3,1,1)
0346         IMF_plot(u,t,0,'Channel','1');
0347         subplot(3,1,2)
0348         IMF_plot(v,t,0,'Channel','2');
0349         subplot(3,1,3)
0350         IMF_plot(w,t,0,'Channel','2');
0351 
0352 
0353      for i=1:param.nimfs   
0354      figure(i+1)    
0355         subplot(3,1,1)
0356         IMF_plot(Results.IMF.u(:,i),t,i,'IMF','Channel 1');
0357         subplot(3,1,2)
0358         IMF_plot(Results.IMF.v(:,i),t,i,'IMF','Channel 2');
0359         subplot(3,1,3)
0360         IMF_plot(Results.IMF.w(:,i),t,i,'IMF','Channel 3');
0361      end
0362     
0363     figure(i+2)  
0364     subplot(3,1,1)
0365     IMF_plot(Results.Residue.u,t,0,'Residue','Channel 1');
0366     subplot(3,1,2)
0367     IMF_plot(Results.Residue.v,t,0,'Residue','Channel 2');
0368     subplot(3,1,3)
0369     IMF_plot(Results.Residue.w,t,0,'Residue','Channel 3');
0370 end
0371 
0372 function IMF_plot(signal,t,imf,name1,name2)    
0373 
0374     plot(t,signal,'-k');
0375     axis ([0 6*pi -5 5]);
0376     xlabel('t');
0377     set(gca,'TickLabelInterpreter','tex')
0378     switch(name1)
0379         case 'IMF'
0380             title(sprintf('%s %d %s',name1,imf,name2));
0381         case 'Channel'
0382             title(sprintf('%s %s',name1,name2));
0383         case 'Residue'
0384             title(sprintf('%s %s',name1,name2));
0385     end
0386     
0387     
0388 end
0389