How can find all peaks in float array?

zhmh

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]];

}

aha_1980

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

J.Hilk

You could use std::max_element

https://en.cppreference.com/w/cpp/algorithm/max_element

aha_1980

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

zhmh

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

J.Hilk

@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();
}

fcarney

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.

SimonSchroeder

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:

(Positive) peaks: 1 2 3 4 5 4 3 2 1 -> 5 is the peak in this sequence
Negative peaks, i.e. valleys: 5 4 3 2 1 2 3 4 5 -> Is 1 a peak in your definition?
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).
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...

Kent-Dorfman

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