Dawn@4: function [Maxima, countMaxima] = findMaxima(f, step)
Dawn@4: 
Dawn@4: %
Dawn@4: % MAXIMA ESTIMATION
Dawn@4: %
Dawn@4: % function [Maxima, countMaxima] = findMaxima(f, step);
Dawn@4: %
Dawn@4: % This function estimates the local maxima of a sequence
Dawn@4: %
Dawn@4: % ARGUMENTS:
Dawn@4: % f: the input sequence
Dawn@4: % step: the size of the "search" window
Dawn@4: %
Dawn@4: % RETURN:
Dawn@4: % Maxima: [2xcountMaxima] matrix containing: 
Dawn@4: %         1. the maxima's indeces
Dawn@4: %         2. tha maxima's values
Dawn@4: % countMaxima: the number of maxima
Dawn@4: %
Dawn@4: 
Dawn@4: 
Dawn@4: 
Dawn@4: %
Dawn@4: % STEP 1: find maxima:
Dawn@4: % 
Dawn@4: 
Dawn@4: countMaxima = 0;
Dawn@4: for (i=1:length(f)-step-1) % for each element of the sequence:
Dawn@4:     if (i>step)
Dawn@4:         if (( mean(f(i-step:i-1))< f(i)) && ( mean(f(i+1:i+step))< f(i)))  
Dawn@4:             % IF the current element is larger than its neighbors (2*step window)
Dawn@4:             % --> keep maximum:
Dawn@4:             countMaxima = countMaxima + 1;
Dawn@4:             Maxima(1,countMaxima) = i;
Dawn@4:             Maxima(2,countMaxima) = f(i);
Dawn@4:         end
Dawn@4:     else
Dawn@4:         if (( mean(f(1:i))<= f(i)) && ( mean(f(i+1:i+step))< f(i)))  
Dawn@4:             % IF the current element is larger than its neighbors (2*step window)
Dawn@4:             % --> keep maximum:
Dawn@4:             countMaxima = countMaxima + 1;
Dawn@4:             Maxima(1,countMaxima) = i;
Dawn@4:             Maxima(2,countMaxima) = f(i);
Dawn@4:         end
Dawn@4:         
Dawn@4:     end
Dawn@4: end
Dawn@4: 
Dawn@4: %
Dawn@4: % STEP 2: post process maxima:
Dawn@4: %
Dawn@4: 
Dawn@4: MaximaNew = [];
Dawn@4: countNewMaxima = 0;
Dawn@4: i = 0;
Dawn@4: while (i<countMaxima)
Dawn@4:     % get current maximum:
Dawn@4:     i = i + 1;
Dawn@4:     curMaxima = Maxima(1,i);
Dawn@4:     curMavVal = Maxima(2,i);
Dawn@4:     
Dawn@4:     tempMax = Maxima(1,i);
Dawn@4:     tempVals = Maxima(2,i);
Dawn@4:     
Dawn@4:     % search for "neighbourh maxima":
Dawn@4:     while ((i<countMaxima) && ( Maxima(1,i+1) - tempMax(end) < step / 2))
Dawn@4:         i = i + 1;
Dawn@4:         tempMax(end+1) = Maxima(1,i);
Dawn@4:         tempVals(end+1) = Maxima(2,i);
Dawn@4:     end
Dawn@4:     
Dawn@4:    
Dawn@4:     % find the maximum value and index from the tempVals array:
Dawn@4:     %MI = findCentroid(tempMax, tempVals); MM = tempVals(MI);
Dawn@4:     
Dawn@4:     [MM, MI] = max(tempVals);
Dawn@4:         
Dawn@4:     if (MM>0.02*mean(f)) % if the current maximum is "large" enough:
Dawn@4:         countNewMaxima = countNewMaxima + 1;   % add maxima
Dawn@4:         % keep the maximum of all maxima in the region:
Dawn@4:         MaximaNew(1,countNewMaxima) = tempMax(MI); 
Dawn@4:         MaximaNew(2,countNewMaxima) = f(MaximaNew(1,countNewMaxima));
Dawn@4:     end        
Dawn@4:     tempMax = [];
Dawn@4:     tempVals = [];
Dawn@4: end
Dawn@4: 
Dawn@4: Maxima = MaximaNew;
Dawn@4: countMaxima = countNewMaxima;
Dawn@4: 
Dawn@4: