Commit b3a2cfbc authored by Neil Vaytet's avatar Neil Vaytet
Browse files

fixed old peak finding method

parent cd9fa9af
......@@ -6,6 +6,12 @@
# __author__ = "Marcos Duarte, https://github.com/demotu/BMC"
# __version__ = "1.0.5"
# __license__ = "MIT"
#
# Modified: 2019/01 - neil.vaytet@esss.se - European Spallation Source
import numpy as np
def detect_peaks(x, mph=None, mpd=1, threshold=0, edge="rising",
kpsh=False, valley=False):
......
......@@ -8,8 +8,9 @@ try:
except ImportError:
from .detect_peaks import detect_peaks
peak_finding = "old"
print("Warning: scipy.signal.find_peaks could not be imported. "
print("Deprecation warning: scipy.signal.find_peaks could not be imported. "
"Reverting to the old peak finding method. "
"This will soon be deprecated. "
"Install scipy>=1.1.0 to fix this problem.")
################################################################################
......@@ -87,18 +88,18 @@ def check_peaks(peaks, nwindows):
"""
This function check if the number of peaks found is the correct one.
"""
nvalleys = len(peaks)-2
if nvalleys > 2*(nwindows - 1):
print("Way too many valleys: should be {}, found {}.".format(nwindows-1, nvalleys))
npeaks = len(peaks)-2
if npeaks > 2*(nwindows - 1):
print("Way too many peaks: should be {}, found {}.".format(nwindows-1, npeaks))
return 3 # Continue looping
elif nvalleys > (nwindows - 1):
print("Too many valleys: should be {}, found {}.".format(nwindows-1, nvalleys))
elif npeaks > (nwindows - 1):
print("Too many peaks: should be {}, found {}.".format(nwindows-1, npeaks))
return 2 # Continue looping
elif nvalleys < (nwindows - 1):
print("Too few valleys: should be {}, found {}.".format(nwindows-1, nvalleys))
elif npeaks < (nwindows - 1):
print("Too few peaks: should be {}, found {}.".format(nwindows-1, npeaks))
return 1 # Continue looping
else:
print("OK: Found {} valleys.".format(nvalleys))
print("OK: Found {} peaks.".format(npeaks))
return 0 # Exit the loop
################################################################################
......@@ -251,15 +252,18 @@ def get_wfm_windows(workspace=None, data=None, filename=None, nwindows=6,
if peak_finding == "scipy":
print("The peak prominence is:", prominence)
# Find valleys using `detect_peaks` function from scipy.signal
# Find valleys using `find_peaks` function from scipy.signal
peaks, _ = find_peaks(-y[lbound:rbound+1], prominence=prominence*(ymax-ymin))
else:
# Find valleys using `detect_peaks` function from Marcos Duarte
print("The minimum peak distance (mpd) is:", mpd)
peaks = detect_peaks(y[lbound:rbound+1], mpd=mpd, valley=True)
# Added start and end to peaks array
peaks = np.concatenate(([i_start], peaks+lbound, [i_end]))
# Make sure we only have peaks between left and right edges
peaks += lbound
peaks = peaks[np.where(np.logical_and(peaks > i_start, peaks < i_end))]
# Add start and end to peaks array
peaks = np.concatenate(([i_start], peaks, [i_end]))
# Check how many peaks were found
loop = check_peaks(peaks, nwindows)
# If loop is 3, it means way too many peaks were found, so we re-run
......@@ -274,7 +278,7 @@ def get_wfm_windows(workspace=None, data=None, filename=None, nwindows=6,
# Try to remove peaks if there are too many
if loop == 2:
print("Too many peaks found. Attempting automatic dismissal")
print("Attempting automatic dismissal")
# if the value at the bottom of the valley is over the lowest window
# average, then this is not a real valley
removed_a_peak = True
......@@ -283,7 +287,7 @@ def get_wfm_windows(workspace=None, data=None, filename=None, nwindows=6,
for p in range(1, len(peaks)-1):
rmean = np.average(y[peaks[p]:peaks[p+1]])
lmean = np.average(y[peaks[p-1]:peaks[p]])
mean = min(lmean,rmean)
mean = min(lmean, rmean)
if y[peaks[p]] > mean:
print("Removed peak number {} at x,y = {},{}".format(p,x[peaks[p]],y[peaks[p]]))
peaks = np.delete(peaks, p)
......@@ -374,7 +378,7 @@ def get_wfm_windows(workspace=None, data=None, filename=None, nwindows=6,
ax.plot(x[peaks[-1]], y[peaks[-1]], "o", color="yellow", label="Trailing edge")
ax.set_xlabel("Arrival time at detector")
ax.set_ylabel("Amplitude")
ax.set_ylim([0.0,1.05*np.amax(data[:,1])])
ax.set_ylim([ymin,1.05*np.amax(data[:,1])])
ax.plot(x[peaks[0]], -ymax, "o", color="r",label="Valleys")
ax.legend(loc=(0,1.02),ncol=4, fontsize=10)
ax2 = fig.add_subplot(122)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment