Home > FA-MVEMD > 1D > EMD1DNV.m

# EMD1DNV

## PURPOSE

Purpose:

## SYNOPSIS

function Results = EMD1DNV(u,t,param)

## DESCRIPTION

     Purpose:
     -To perform EMD on n(3-16) channels of 1 dimensional data

     Input:
     - u: Signal of size (length, number of channels)
     - 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
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

## CROSS-REFERENCE INFORMATION

This function calls:

- extrema EXTREMA Gets the global extrema points from a time series.
- filter_size1D

This function is called by:

- driver_1dnv_example cleanup

## SUBFUNCTIONS

- function Windows = filter_size1D(imax, imin, type)
- function H1 = Sift(H,w_sz)
- function [Max,Min] = OSF(H,w_sz)
- function f_signal = Ordfilt1(signal,order,window_size)
- function Env_med = Pad_smooth(Env_max,Env_min,w_sz)
- function [IO,Error] = Orth_index(Signal,IMF,Residue)
- function Plot_results(u,t,Results)
- function IMF_plot(signal,t,imf,name1,name2)

## SOURCE CODE

    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