Search code examples
pythonpeakutils

Why does peakutils.peak.indexes() seem to ignore the provided threshold value?

I'm retrieving the arrays holding the power levels and frequencies, respectively, of a signal from the plt.psd() method:

Pxx, freqs = plt.psd(signals[0], NFFT=2048, Fs=sdr.sample_rate/1e6, Fc=sdr.center_freq/1e6, scale_by_freq=True, color="green")

enter image description here Please ignore the green and red signals. Just the blue one is relevant for this question.

I'm able to have the peakutils.peak.indexes() method return the X and Y coordinates of a number of the most significant peaks (of the blue signal):

power_lvls = 10*log10(Pxx/(sdr.sample_rate/1e6))
indexes = peakutils.peak.indexes(np.array(power_lvls), thres=0.6/max(power_lvls), min_dist=120)
print("\nX: {}\n\nY: {}\n".format(freqs[indexes], np.array(power_lvls)[indexes]))

enter image description here As can be seen, the coordinates fit the blue peaks quite nicely.

What I'm not satisfied with is the number of peak coordinates I receive from the peak.indexes() method. I'd like to have only the coordinates of all peaks above a certain power level returned, e.g., -25 (which would then be exactly 5 peaks for the blue signal). According to the documentation of the peak.indexes() method this is done by providing the desired value as thres parameter.

But no matter what I try as thres, the method seems to entirely ignore my value and instead solely rely on the min_dist parameter to determine the number of returned peaks.

  • What is wrong with my threshold value (which I believe means "peaks above the lower 60% of the plot" in my code now) and how do I correctly specify a certain power level (instead of a percentage value)?

[EDIT]

I figured out that apparently the thres parameter can only take positive values between float 0. and 1. So, by changing my line slightly as follows I can now influence the number of returned peaks as desired:

indexes = peakutils.peak.indexes(np.array(power_lvls), thres=0.4, min_dist=1)

But that still leaves me with the question whether it's possible to somehow limit the result to the five highest peaks (provided num_of_peaks above thres >= 5).

I believe something like the following would return the five highest values:

print(power_lvls[np.argsort(power_lvls[indexes])[-5:]])

Unfortunately, though, negative values seem to be interpreted as the highest values in my power_lvls array. Can this line be changed such that (+)10 would be considered higher than, e.g., -40? Or is there another (better?) solution?

[EDIT 2]

These are the values I get as the six "highest" peaks: enter image description here

power_lvls = 10*log10(Pxx/(sdr.sample_rate/1e6))+10*log10(8/3)
indexes = peakutils.indexes(power_lvls, thres=0.35, min_dist=1)

power_lvls_max = power_lvls[np.argsort(power_lvls[indexes])[-6:]]
print("Highest Peaks in Signal:\nX: \n\nY: {}\n".format(power_lvls_max))

After trying various things for hours without any improvement I'm starting to think that these are neither valleys nor peaks, just some "random" values?! Which leads me to believe that there is a problem with my argsort line that I have to figure out first?!

[EDIT 3]

The bottleneck.partition() method seems to return the correct values (even if apparently it does so in random order, not from leftmost peak to rightmost peak):

import bottleneck as bn
power_lvls_max = -bn.partition(-power_lvls[indexes], 6)[:6]

enter image description here Luckily, the order of the peaks is not important for what I have planned to do with the coordinates. I do, however, have to figure out yet how to match the Y values I have now to their corresponding X values ...

Also, while I do have a solution now, for learning purposes it would still be interesting to know what was wrong with my argsort attempt.

Solution

  • I figured it out how to find the corresponding X values and get full coordinates of the six highest peaks:

    power_lvls = 10*log10(Pxx/(sdr.sample_rate/1e6))+10*log10(8/3)
indexes = peakutils.indexes(power_lvls, thres=0.35, min_dist=1)
print("Peaks in Signal 1\nX: {}\n\nY: {}\n".format(freqs[indexes], power_lvls[indexes]))
power_lvls_max = -bn.partition(-power_lvls[indexes], 6)[:6]
check = np.isin(power_lvls, power_lvls_max)
indexes_max = np.where(check)
print("Highest Peaks in Signal 1:\nX: {}\n\nY: {}\n".format(freqs[indexes_max], power_lvls[indexes_max]))

    Now I have my "peak filtering" (kind of), which I originally tried to achieve by messing around with the thres value of peakutils.peak.indexes(). The code above gives me just the desired result: enter image description here