EMD1D2V

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

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