extract_peak: Extract Peak
In majpark21/TSexploreR: Explorative Tools for Time Series

Description Usage Arguments Details Value See Also Examples

A wrapper for different peak extraction methods. In addition, has some options to detect peak position to give a more precise estimate of peak position.

extract_peak(y, method = c("spike", "bump"), all.winlen = NULL,
  z.winlen = NULL, z.nsd = 3, z.influence = 0.5, z.robust = TRUE,
  sp.winlen = NULL, sp.nsd = 3, sp.roi = NULL, bp.winlen = NULL,
  bp.mind = NULL, bp.maxderiv = 0.04, bp.nsd = 0.5, bp.npos = 10L,
  wt.widths = NULL, wt.snr = 1, wt.centile_noise = 10,
  ext.winlen = NULL, ext.what = "maxi", recenter = 2,
  filter.thresh = NULL, filter.winlen = all.winlen * 2,
  detect.early = NULL, show.warning = TRUE)

`y`	numeric vector. The series from which to extract peak positions.
`method`	Peak extraction method(s). One or several of ("zscore", "spike", "bump", "wavelet", "extrema").
`all.winlen`	numeric. Set a common window length for all methods relying on rolling mean. Will overwrite any xxx.winlen argument.
`z.winlen`	numeric. Window length for rolling mean in "zscore" method.
`z.nsd`	numeric. Number of standard deviation away from rolling mean for "zscore" method.
`z.influence`	numeric. How much should peaks affect rolling mean? 0 means that they are not taken into account for computing rolling mean. 1 means that the rolling mean is computed without regards of points being peaks or not. Used in "zscore" method.
`z.robust`	logical. If TRUE rolling median and MAD replace rolling mean and standard deviation. Used in "zscore" method.
`sp.winlen`	numeric. Window length for rolling mean in "spike" method.
`sp.nsd`	numeric. Number of standard deviation away from rolling mean for "spike" method.
`sp.roi`	numeric vector of length 2. Where to look for spikes in "spike" method.
`bp.winlen`	numeric. Window length for rolling mean in "bump" method.
`bp.mind`	numeric. Minimum distance between peaks for "bump" method. See ?peakpick::peakpick - neighlim.
`bp.maxderiv`	numeric. Upper limit for the estimatied derivative for a point to be considered for a peak call. See ?peakpick::peakpick - deriv.lim.
`bp.nsd`	numeric. Number of standard deviation away from rolling mean for "bump" method.
`bp.npos`	integer. Peak standard deviations and means will be estimated plus/minus npos positions from peak. See ?peakpick::peakpick - peak.npos.
`wt.widths`	numeric vector. widths to use for calculating the CWT matrix. In general, this range should cover the expected width of peaks of interest. for "wavelet" method. See scipy.signal.find_peaks_cwt documentation.
`wt.snr`	Minimum SNR ratio. See scipy.signal.find_peaks_cwt documentation.
`wt.centile_noise`	When calculating the noise floor, percentile of data points examined below which to consider noise. See scipy.signal.find_peaks_cwt documentation.
`ext.winlen`	Window length for looking for local extrema. For "extrema" method.
`ext.what`	One of "mini" or "maxi". For "extrema" method.
`recenter`	numeric or NULL. If numeric, attempt to reposition peak position to a more precise location. Specifically, starting from peak estimated position, look for maximum in a window centered around it. The length of the window corresponds to the parameter value. Highly recommended for "spike" and "bump" method which tend to give an imprecise estimation of peak position.
`filter.thresh`	numeric or NULL. If numeric, attempt at filtering peaks on the base of a minimal amplitude. To do so, TS are detrended by subtracting a rolling mean and amplitude to baseline is computed. Peak is kept if this amplitude is above the "filter" value. If NULL, no filter is applied. It is highly recommended to recenter the peaks if this option is used.
`filter.winlen`	size of the window for moving average used to detrend prior to filtering. It is highly recommended to check the detrended signal (part of the output if filter is specified): in the detrended signal, area between peaks should oscillate around 0. If not, try to adjust the value.
`detect.early`	numeric or NULL. If numeric, extrapolate the time series by the corresponding number of points at its start. This aims at improving peak detection at the beginning of the series. Extrapolation is performed by ARIMA model (see ?forecast::auto.arima)
`show.warning`	Logical, disable warning for the current execution. Useful when providing all.winlen and applying the function multiple times.

Zscore method is based on https://stackoverflow.com/questions/22583391/peak-signal-detection-in-realtime-timeseries-data#22640362.

Spike and Bump methods are wrapper for the function of the peakPick package: https://cran.r-project.org/web/packages/peakPick/peakPick.pdf. For the bump method, a rolling mean is first applied to the signal as recommended by the package authors.

The wavelet method relies on scipy (python module) peak finder https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks_cwt.html.

A data table containing the signal along with results all of selected methods. Each method result is stored in a different column. If filter is specified, also returns the detrended signal.

isolate_peak

# Curve coming from real-data:
y <- c(1.114, 1.146, 1.106, 1.049, 1.031, 1.01, 1, 0.992, 0.987, 0.981, 1.009, 1.072, 1.053, 1.018, 0.988,
 0.979, 0.965, 0.958, 0.956, 0.942, 0.939, 0.932, 0.926, 0.922, 0.919, 0.919, 0.916, 0.915, 0.918, 0.916,
 0.913, 0.911, 0.907, 0.908, 0.904, 0.905, 0.9, 0.899, 0.899, 0.905, 0.905, 0.908, 0.913, 0.919, 0.919,
 0.911, 0.922, 0.918, 0.916, 0.915, 0.92, 0.917, 0.915, 0.926, 0.964, 1.099, 1.075, 1.018, 0.999, 0.989,
 0.977, 0.975, 0.964, 0.96, 0.958, 0.945, 0.946, 0.944, 0.944, 0.938, 0.934, 0.936, 0.933, 0.926, 0.921,
 0.923, 0.924, 0.926, 0.929, 0.932, 0.926, 0.932, 0.939, 0.941, 0.976, 1.116, 1.09, 1.053, 1.034, 1.019,
 1.038, 1.104, 1.089, 1.059, 1.027, 1.012, 1.009, 1, 0.991, 0.986, 0.973, 0.969, 0.975, 0.966, 0.961,
 0.96, 0.949, 0.946, 0.943, 0.938, 0.942, 0.937, 0.934, 0.936, 0.941, 0.939, 0.97, 0.932, 0.931, 0.941,
 0.944, 0.949, 0.969, 1.048, 1.052, 1.059, 1.103, 1.083, 1.05, 1.036, 1.015, 1.01, 0.997, 0.995, 0.988,
 0.977, 0.979, 0.977, 0.981, 0.982, 0.976, 0.971, 0.97, 0.965, 0.966, 0.967, 0.967, 0.972, 0.904, 0.966,
 0.978, 0.983, 0.979, 0.975, 0.98, 0.973, 0.968, 0.97, 0.962, 0.966, 0.962, 0.968, 0.964, 0.956, 0.963,
 0.954, 0.958, 0.956, 0.958, 0.952, 0.95, 0.952, 0.946, 0.953, 0.952, 0.952, 0.953, 0.952, 0.954, 0.956)

# Extract peaks with "vanilla" bump method. First smooth signal with a moving average of width 5, peaks must be
# separated with a minimal distance of 2 points.
peaks_vanilla <- extract_peak(y, method = "bump", bp.winlen = 5, bp.mind = 1, recenter = NULL, filter.thresh = NULL)
# Add: Recenter peaks to the maximum in a neigbourhood of 2 points.
peaks_center <- extract_peak(y, method = "bump", bp.winlen = 5, bp.mind = 1, recenter = 2, filter.thresh = NULL)
# Add: detrend the signal with a moving average of width 4, and filter peaks which amplitude is below 0.05 in the detrended signal.
peaks_center_filter <- extract_peak(y, method = "bump", bp.winlen = 5, bp.mind = 1, recenter = 2, filter.thresh = 0.05, filter.winlen = 10)
# Add: detection of early peak by extrapolating 10 points
peaks_final <- extract_peak(y, method = "bump", bp.winlen = 5, bp.mind = 1, recenter = 2, filter.thresh = 0.05, filter.winlen = 10, detect.early = 10)

par(mfrow=c(3,2))
plot(y, type = "b", main = "Vanilla bump detection")
abline(v=which(peaks_vanilla$bump), col = "red", lty = "dashed")
plot(y, type = "b", main = "Peak recentering")
abline(v=which(peaks_center$bump), col = "red", lty = "dashed")
plot(y, type = "b", main = "Peak center + filter - Raw data")
abline(v=which(peaks_center_filter$bump), col = "red", lty = "dashed")
plot(peaks_center_filter$detrended.signal, type = "b", main = "Peak center + filter - Detrended data\nCheck that baseline is around 0")
abline(v=which(peaks_center_filter$bump), col = "red", lty = "dashed")
plot(y, type = "b", main = "Peak center + filter + early detection")
abline(v=which(peaks_final$bump), col = "red", lty = "dashed")
# Plot arima prediction
y_extend <- rev(c(rev(y), forecast(auto.arima(rev(y)), h = 10)$mean))
plot(y_extend, main = "Extended Series for early peak detection.\nBlue dots are predicted by ARIMA.", type = "b")
points(1:10, y_extend[1:10], pch=20, col="blue")