Important: Please read the Qt Code of Conduct - https://forum.qt.io/topic/113070/qt-code-of-conduct

How can find all peaks in float array?



  • I want to find peak in the float array,In this code I have get wrong output when there is a peak on an endpoint.How can I fix it?

    // Util.cpp
    #include "Util.h"
    
    void diff(vector<float> in, vector<float>& out)
    {
    	out = vector<float>(in.size()-1);
    
        for(int i=1; i<(int)in.size(); ++i)
    		out[i-1] = in[i] - in[i-1];
    }
    
    void vectorProduct(vector<float> a, vector<float> b, vector<float>& out)
    {
    	out = vector<float>(a.size());
    
        for(int i=0; i<(int)a.size(); ++i)
    		out[i] = a[i] * b[i];
    }
    
    void findIndicesLessThan(vector<float> in, float threshold, vector<int>& indices)
    {
        for(int i=0; i<(int)in.size(); ++i)
    		if(in[i]<threshold)
    			indices.push_back(i+1);
    }
    
    void selectElements(vector<float> in, vector<int> indices, vector<float>& out)
    {
        for(int i=0; i<(int)indices.size(); ++i)
    		out.push_back(in[indices[i]]);
    }
    
    void selectElements(vector<int> in, vector<int> indices, vector<int>& out)
    {
        for(int i=0; i<(int)indices.size(); ++i)
    		out.push_back(in[indices[i]]);
    }
    
    void signVector(vector<float> in, vector<int>& out)
    {
    	out = vector<int>(in.size());
    
        for(int i=0; i<(int)in.size(); ++i)
    	{
    		if(in[i]>0)
    			out[i]=1;
    		else if(in[i]<0)
    			out[i]=-1;
    		else
    			out[i]=0;
    	}
    }
    
    
    void Peaks::findPeaks(vector<float> x0, vector<int>& peakInds)
    {
    	int minIdx = distance(x0.begin(), min_element(x0.begin(), x0.end()));
    	int maxIdx = distance(x0.begin(), max_element(x0.begin(), x0.end()));
    
    	float sel = (x0[maxIdx]-x0[minIdx])/4.0;
    
    	int len0 = x0.size();
    
    	vector<float> dx;
    	diff(x0, dx);
    	replace(dx.begin(), dx.end(), 0.0f, -Peaks::EPS);
    	vector<float> dx0(dx.begin(), dx.end()-1);
    	vector<float> dx1(dx.begin()+1, dx.end());
    	vector<float> dx2;
    
    	vectorProduct(dx0, dx1, dx2);
    
    	vector<int> ind;
    	findIndicesLessThan(dx2, 0, ind); // Find where the derivative changes sign
    	
    	vector<float> x;
    
    	vector<int> indAux(ind.begin(), ind.end());
    	selectElements(x0, indAux, x);
    	x.insert(x.begin(), x0[0]);
    	x.insert(x.end(), x0[x0.size()-1]);;
    
    
    	ind.insert(ind.begin(), 0);
    	ind.insert(ind.end(), len0);
    
    	int minMagIdx = distance(x.begin(), min_element(x.begin(), x.end()));
    	float minMag = x[minMagIdx];
    	float leftMin = minMag;
    	int len = x.size();
    
    	if(len>2)
    	{
    		float tempMag = minMag;
        	bool foundPeak = false;
        	int ii;
    
        	// Deal with first point a little differently since tacked it on
            // Calculate the sign of the derivative since we tacked the first
            //  point on it does not neccessarily alternate like the rest.
        	vector<float> xSub0(x.begin(), x.begin()+3);//tener cuidado subvector
        	vector<float> xDiff;//tener cuidado subvector
        	diff(xSub0, xDiff);
    
        	vector<int> signDx;
        	signVector(xDiff, signDx);
    
            if (signDx[0] <= 0) // The first point is larger or equal to the second
            {
                if (signDx[0] == signDx[1]) // Want alternating signs
                {
                    x.erase(x.begin()+1);
                    ind.erase(ind.begin()+1);
                    len = len-1;
                }
            }
            else // First point is smaller than the second
            {
                if (signDx[0] == signDx[1]) // Want alternating signs
                {
                	x.erase(x.begin());
                	ind.erase(ind.begin());
                    len = len-1;
                }
            }
    
        	if ( x[0] >= x[1] )
            	ii = 0;
        	else
            	ii = 1;
    
        	float maxPeaks = ceil((float)len/2.0);
        	vector<int> peakLoc(maxPeaks,0);
        	vector<float> peakMag(maxPeaks,0.0);
        	int cInd = 1;
        	int tempLoc;
        
        	while(ii < len)
        	{
            	ii = ii+1;//This is a peak
            	//Reset peak finding if we had a peak and the next peak is bigger
            	//than the last or the left min was small enough to reset.
            	if(foundPeak)
            	{
                	tempMag = minMag;
                	foundPeak = false;
                }
            
            	//Found new peak that was lager than temp mag and selectivity larger
            	//than the minimum to its left.
            
            	if( x[ii-1] > tempMag && x[ii-1] > leftMin + sel )
            	{
                	tempLoc = ii-1;
                	tempMag = x[ii-1];
            	}
    
            	//Make sure we don't iterate past the length of our vector
            	if(ii == len)
                	break; //We assign the last point differently out of the loop
    
            	ii = ii+1; // Move onto the valley
            	
            	//Come down at least sel from peak
            	if(!foundPeak && tempMag > sel + x[ii-1])
                {            	
    	            foundPeak = true; //We have found a peak
    	            leftMin = x[ii-1];
    	            peakLoc[cInd-1] = tempLoc; // Add peak to index
    	            peakMag[cInd-1] = tempMag;
    	            cInd = cInd+1;
    	        }
            	else if(x[ii-1] < leftMin) // New left minima
                	leftMin = x[ii-1];
                
            }
    
            // Check end point
            if ( x[x.size()-1] > tempMag && x[x.size()-1] > leftMin + sel )
            {
                peakLoc[cInd-1] = len-1;
                peakMag[cInd-1] = x[x.size()-1];
                cInd = cInd + 1;
            }
            else if( !foundPeak && tempMag > minMag )// Check if we still need to add the last point
            {
                peakLoc[cInd-1] = tempLoc;
                peakMag[cInd-1] = tempMag;
                cInd = cInd + 1;
            }
    
        	//Create output
        	if( cInd > 0 )
            {
    
    
            	vector<int> peakLocTmp(peakLoc.begin(), peakLoc.begin()+cInd-1);
                        selectElements(ind, peakLocTmp, peakInds);
    
            	//peakMags = vector<float>(peakLoc.begin(), peakLoc.begin()+cInd-1);
            }    	
    	}
    }
    
    

    in test code I set a peak at end of array and the output is :

    1.26083
    0.975976
    9.10844e-44 //wrong number print
    
    #include <iostream>
    #include "Util.h"
    
    using namespace std;
    
    
    int main()
    {
    	float signal[14] = {1.22933,1.07635,1.23343,1.17653,1.2322,1.26083,0.371609,0.975976,0.539859,0.039658,0.039658,1.039658,0.949407,1.64061};
    	vector<float> in(signal, signal + sizeof(signal) / sizeof(float));
    	vector<int> out;
    
        Peaks::findPeaks(in, out);
    
    	qDebug()<<"Maxima found:";
    
    	for(int i=0; i<out.size(); ++i)
    		qDebug()<<in[out[i]];
    
    }
    

  • Lifetime Qt Champion

    Hi @zhmh,

    at least for the diffing there is a std algorithm: https://en.cppreference.com/w/cpp/algorithm/adjacent_difference

    Afterwards, you would just need to find the maxima, for which an algorigthm probably exists too, I'm just not aware of.

    Regards


  • Moderators


  • Lifetime Qt Champion

    @J-Hilk Then I'd say he should use a combination of both algorithms :)



  • @J-Hilk the std::max_element just find one peak I want to find all peaks in array


  • Moderators

    @zhmh than you should adjust the title of the topic :P

    int main(int argc, char *argv[])
    {
        QApplication a(argc, argv);
        std::srand(std::time(0));
    
        std::vector<float> input;
        for(int i = 0; i < 100; i ++)
            input.push_back(rand()/static_cast<float>(RAND_MAX));
    
        std::vector<float> output;
        std::copy_if(input.cbegin(), input.cend(), std::back_inserter(output),[threshold = 0.6](const auto value) ->bool{ return value > threshold;});
    
        qDebug() << output;
    //    return a.exec();
    }
    


  • If you define a peak as:
    n is peak if abs(n-1) < abs(n) and abs(n+1) < abs(n)

    That means the end points are special cases. So either add zero to beginning and end of your array or handle the special cases separately.



  • I guess you need to first make clear what you understand to be peaks. At first glance your code seems to be overly complicated. From a quick look you seem to use a mathematical definition for peaks as extremas of a continuous function defined by points where the first derivative is 0. Let's have a look at a few examples to help us understand what you define as peaks. Here is my understanding:

    1. (Positive) peaks: 1 2 3 4 5 4 3 2 1 -> 5 is the peak in this sequence
    2. Negative peaks, i.e. valleys: 5 4 3 2 1 2 3 4 5 -> Is 1 a peak in your definition?
    3. Plateaus (i.e. not peaks): 1 2 3 4 4 4 5 6 -> None of the 4's is a peak because we go upwards after this again. Special case: 1 2 3 4 5 5 5 4 3 2 1 -> I assume there should be a peak here, but which 5? Take the middle one? What if we have an even number of 5's in this example? Maybe you are talking about a real function definition for your sequence so that plateaus will not occur (if you have strongly monotonic functions).
    4. When are the first and last number in a sequence defined to be peaks? Without plateaus they are either always or never peaks. If you allow for plateaus this is a lot more complicated... Mathematically speaking the ends of your sequence (not considering plateaus) are always local extrema.

    Number 1 is quite easy to check:

    void Peaks::findPeaks(vector<float> x0, vector<int>& peakInds)
    {
        for(int i = 1; i < x0.size()-1; ++i)
        {
            if(x0[i] > x0[i-1] && x0[i] > x0[i+1])
                peakInds.push_back(i);
            ...
        }
    }
    

    If you want to include number 2, just add this line inside the loop:

    if(x0[i] < x0[i-1] && x0[i] < x0[i+1])
        peakInds.push_back(i)
    

    Plateaus would be a lot more complicated to detect. Let us know if you want to include those as well. And just as I said: you can either always or never include the first and last element.

    BTW, there are two things about your code style I would change. First, there is no reason to put findPeaks into a class. If it is a namespace it would be okay, though. Second, with modern C++ there is no reason to use an out-parameter to your function. Instead of

    void Peaks::findPeaks(vector<float> x0, vector<int> &peakInds)
    {
        ...
    }
    

    you should write:

    vector<int> Peaks::findPeaks(const vector<float> &x0)
    {
        vector<int> peakInds;
        ...
        return peakInds;
    }
    

    There are now things like return-value-optimization and copy-ellision which will make returning a vector as fast as using out-parameters. Your approach only is advantageous when you want to add indices to an existing vector. Also note that I used const vector<float> &x0 to avoid a copy of this vector. This is just good modern style...



  • this problem would make an cool interview exercise for someone expected to write waveform analysis software. LOL


Log in to reply