Searching for specific maxima

Hello Everyone,

I have a specific question. Let us say I have some random data set that looks like this:

1 2 1 2 3 2 3 4 122 3 4 5 3 51 4 2

I am expecting two "maxima" or peaks (not necessarily 122 or 51). I want to be able to autonomously detect the two largest maxima, save the two maxima values and the index where they occur. Is there an easy way to do this?

Thanks for your help.

[Walter Roberson]

What if more than one location has the same maximum (or second highest) value?

Answer 1

[the cyclist]

x = [1 2 1 2 3 2 3 4 122 3 4 5 3 51 4 2];
[sorted_x indexToSortedValues] = sort(x,'descend');

This will sort the values from high-to-low, and tell you the indices of the sorted values. Be wary of the potential uncertainty that Walter mentions in his comment.

[Sarah]

Hey thanks a lot! This is so simple, I really should have thought of that! I made an algorithm, but it was waaay too complex. So referring to the problem that Walter raised, what can I do to make it more robust?

The way I have set up my code, the parameters will be constantly changing. I am detecting peaks in a signal. So when the parameters change, I will have new maxima for every iteration...though I still expect two peaks for each iteration. But after a certain point in my analysis, the data will converge to only one peak. What can I do to make the code automatically sense that only one prominent "peak" exists, and stop at that iteration?

[Walter Roberson]

There are a number of peak-picking routines in the MATLAB File Exchange; they should work better than looking at the maxima.

[Image Analyst]

I agree with Walter. The cyclist's code will do what he says but it won't necessarily find peaks or help you locate them, unless your second peak is the second highest value. For example if you still have two peaks but your data looks like x = [1 2 1 2 3 2 3 4 122 120 4 5 3 51 4 2] your second value will be 120, which was part of the first peak, and your second peak is now located as the third value of the sorted array - probably not where you expected to find it. (This is a different watchout than Walter's warning of multiple elements with the same value.)

[the cyclist]

I have to admit I totally missed the fact that you seem to be looking for localized peaks, as opposed to just the two largest values (i.e. maxima). There was an FEX "Pick of the Week" for finding local extrema, as discussed here: http://blogs.mathworks.com/pick/2008/05/09/finding-local-extrema/

[Sarah]

Sorry guys, I should have been more clear, it is my fault. So I looked at two peak detection functions: findpeaks (built-in) and extrema (what the cyclist specified). Both do a good job of detecting the peaks. The problem is, are they adaptive? My signal is changing with every iteration, and my expectation is that eventually the two significant peaks that I care about will eventually become one. I want the peak detection function to realize this. Is that possible?

Answer 2

[Andrei Bobrov]

A = [1 2 1 2 3 2 3 4 122 120 3 4 5 3 51 4 2]
[pks,locs] = findpeaks(A)
[ix,ix] = sort(pks,'descend')
outpks = pks(ix(1:2))
outidx = locs(ix(1:2))

[the cyclist]

Note that the findpeaks() function is part of the Signal Processing Toolbox, so you may not have access to it.

[Sarah]

I have access to it, and it works pretty well for me, but I don't think it is robust enough to detect the transition from two peaks to one peak.

[Andrei Bobrov]

Hi cyclist! Variant without 'findpeaks'.

i1 = diff(A(:).')>0;

ix = strfind(i1,[1 0]);

[c i2] = sort(A(ix),'descend')

outpks = c(1:2)

outidx = ix(i2(1:2))

Answer 3

[Dr. Seis]

Similar to andrei:

A = [1 2 1 2 3 2 3 4 122 120 3 4 5 3 51 4 2]
[pks,locs] = findpeaks(A)
maxpks = max(pks);
threshold = 0.4; % 40% threshold
outpks = pks(pks>threshold*maxpks);
outidx = locs(pks>threshold*maxpks);

Set the threshold so that it will only pick up either: 1) the biggest two picks or 2) one peak

[Sarah]

Thanks! But how did you choose 40% for the threshold, or was that just a random value? How do you know that your threshold will be robust enough to truly differentiate the true number of peaks?

[Dr. Seis]

The threshold value is data dependent. It was just a random value that I picked that would yield the largest of the two peaks and would yield only one peak should those two transition together (since if they joined together there would be no other individual peak greater than 40% of 122+51). If, in reality, the peaks are going to be a lot closer in amplitude than the difference between 122 and 51, then you can make that threshold larger (like 80%). However, if you are expecting a larger range of differences between the largest peak and the second largest peak (before they transition into one peak) then you will need to set that threshold lower (like 20-30%).

[Sarah]

This works, but it still does not robustly detect the peaks. :( I have a more clear explanation below, if you don't mind taking a look at that...

Answer 4

[Sarah]

Hey guys,

I am afraid that I still haven't found a solution to my problem :( Let me just explain what I am doing first to make everything clear.

I have two signals: a sinusoidal signal with unity frequency and amplitude, and a "rectangular pulse" signal. The pulse has a constant baseline for everywhere except a specific time range, where it becomes one rectangular pulse with a user specified amplitude. I combine these two signals together. Here is what that looks like:

http://i.imgur.com/4rc0c.jpg

I process this data and transform it into the frequency domain. Here is a link which shows what the frequency domain looks like with the combined signal.

http://i.imgur.com/IdvKe.jpg

So, there are clearly two peaks. The first peak represents the exact start of the rectangular pulse, and the second peak represents the exact end of the rectangular pulse. Everywhere else is generally constant because of the sinusoidal function.

Now, this is my study. At a specific sample time, I am gradually decreasing the width of the pulse using a for loop. At some point, it will not be clear to see those two peaks like you did before. At some specified width, you might see something like this, for example:

http://i.imgur.com/RzcOW.jpg

Clearly, though the point representing the second peak is still there, it is difficult to determine and can be construed as one peak. I would like to detect the value of the WIDTH at which this phenomenon occurs. How can I make a robust peak detection algorithm to accurately do this? Simply setting thresholds for findpeaks does not work...any help would be greatly appreciated.

[Dr. Seis]

You don't happen to have a short example (like the size you posted before) where the threshold method breaks down? Is it that in certain circumstances the second peak falls way below the defined threshold... or does findpeaks not actually see a double peak when those two peaks begin to merge?

[Sarah]

I don't think I can come up with a short example. The peak detection algorithms work like they are supposed to. But if you look in this image here:

http://i.imgur.com/RzcOW.jpg

It is much harder to discern two separate peaks, because the area around the base appears to be "mixed". Does that make sense?

I am measuring the limits of temporal resolution of a particular algorithm that I am using.

[Image Analyst]

It looks like you're plotting only the markers and maybe they oscillate on an element by element basis giving the appearance of two separate plots. Can you set xlim([9.475 9.525]) and plot with lines, like plot(data, 'bo-') so we can see what's going on?

[Sarah]

I did what you suggested, and here is the result:

http://i.imgur.com/8K2ZZ.jpg

So clearly, there is one prominent peak. Here is the code that I have to detect peaks:

[Peaks,Index] = findpeaks(Freq(1,:));

MaxPeak = max(Freq(1,:));

OutPks = Peaks(Peaks>Threshold*MaxPeak);

OutIdx = Index(Peaks>Threshold*MaxPeak);

[OutRow,OutColumn] = size(OutPks);

I can set the threshold to some value so that only peaks that are greater than Threshold*MaxPeak is detected. The problem with that is, the peaks fluctuate with each iteration of width, so it will confuse the current algorithm that I have. I need something more robust than that.

[Sarah]

If I follow the threshold approach, I need to set my threshold to some "robust" value so that I can detect two peaks for when there really are two peaks, and then just one peak if the case is like what I provided above in my most recent previous post. Does that make any sense?

[Dr. Seis]

Sorry, one more plot. Can you plot <http://i.imgur.com/RzcOW.jpg> using the same limits as <http://i.imgur.com/8K2ZZ.jpg> (i.e., [9.475 9.525])? So, just to be sure, in the earlier picture where there are two peaks, the robust peak finder should in fact return two peaks?

[Sarah]

Hey, I'm sorry, I forgot the values I inputted to get that exact plot :( But to answer your second question, yes, if there are two peaks like in that image, the robust peak finder should return both peaks and the indexes. If there is only one peak like in the second image link, the peak finder should detect and return one peak and one index corresponding to the peak. I feel like my code is close, but I just don't really know how to code that algorithm...

[Sarah]

In the following case, I have set the width of the rectangular pulse to Wid = 0.8935. Here is the plot of the frequency content:

http://i.imgur.com/PpGQj.jpg

There are two distinct peaks that correspond to the starting and ending points of the rectangular pulse. However, my current peak detection algorithm only returns one peak:

Peak: 999.6653

Index: 9488

It should have returned two peaks there.

[Dr. Seis]

So... I think I see the problem. There are going to be cases where you might need a threshold of 10% so that (as in the case immediately above) it will detect both the primary and the secondary peak. But (maybe) in other cases a threshold of 10% will return peaks that it should not have.

[Dr. Seis]

Is there, perhaps, a window of indices around the main peak that it should check? For example, it should only check for extra peaks +/- 500 indices from the main peak?

[Sarah]

Yes, you are right. That is my problem, and why I say that the detection algorithm isnt quite robust enough to tackle this type of problem.

And as for your second question, I am not sure, as the values fluctuate a lot. But now that you bring that up, maybe I can try? I don't really know how I could define a "robust" range though. Do you have any idea?

[Dr. Seis]

Try a range of 1.1 times the number of indices associated with the width of your pulse.

Answer 5

[Dr. Seis]

Next idea, instead of envelope of Freq we replace any negative (or really small) frequencies with linear interpolation of positive ones:

tempFreq = Freq(1,:);
maxFreq = max(tempFreq);
for i = 2 : length(tempFreq)-1
  if tempFreq(i) <= 0.01*maxFreq
      tempFreq(i) = mean(tempFreq([i-1,i+1]));
  end
end
[Peaks,Index] = findpeaks(tempFreq); 
[MaxPeak,IndPeak] = max(Peaks); 
IndRange = abs(Index(IndPeak)-Index);
IndThres = 500; % Set this value to the number of indices 
              % associated with your pulse width * 1.1         
OutPks = Peaks( IndRange <= IndThres );
OutPks = Index( IndRange <= IndThres ); 
[OutRow,OutColumn] = size(OutPks);

You will still need to define IndThres to be the about 1.1 times your pulse width above.