Description of EMD1DNV

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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
