python时间序列峰值检测_算法-实时时间序列中的峰值信号检测d

2024年1月29日发(作者：)

function [signals,avgFilter,stdFilter] = ThresholdingAlgo(y,lag,threshold,influence)% Initialise signal resultssignals = zeros(length(y),1);% Initialise filtered seriesfilteredY = y(1:lag+1);% Initialise filtersavgFilter(lag+1,1) = mean(y(1:lag+1));stdFilter(lag+1,1) = std(y(1:lag+1));% Loop over all datapoints y(lag+2),...,y(t)for i=lag+2:length(y)% If new value is a specified number of deviations awayif abs(y(i)-avgFilter(i-1)) > threshold*stdFilter(i-1)if y(i) > avgFilter(i-1)% Positive signalsignals(i) = 1;else% Negative signalsignals(i) = -1;end% Make influence lowerfilteredY(i) = influence*y(i)+(1-influence)*filteredY(i-1);else% No signalsignals(i) = 0;filteredY(i) = y(i);end% Adjust the filtersavgFilter(i) = mean(filteredY(i-lag:i));stdFilter(i) = std(filteredY(i-lag:i));end% Done, now return resultsend例：% Data

y = [1 1 1.1 1 0.9 1 1 1.1 1 0.9 1 1.1 1 1 0.9 1 1 1.1 1 1,...1 1 1.1 0.9 1 1.1 1 1 0.9 1 1.1 1 1 1.1 1 0.8 0.9 1 1.2 0.9 1,...1 1.1 1.2 1 1.5 1 3 2 5 3 2 1 1 1 0.9 1,...1 3 2.6 4 3 3.2 2 1 1 0.8 4 4 2 2.5 1 1 1];% Settingslag = 30;threshold = 5;influence = 0;% Get results[signals,avg,dev] = ThresholdingAlgo(y,lag,threshold,influence);figure; subplot(2,1,1); hold on;x = 1:length(y); ix = lag+1:length(y);area(x(ix),avg(ix)+threshold*dev(ix),'FaceColor',[0.9 0.9 0.9],'EdgeColor','none');area(x(ix),avg(ix)-threshold*dev(ix),'FaceColor',[1 1 1],'EdgeColor','none');plot(x(ix),avg(ix),'LineWidth',1,'Color','cyan','LineWidth',1.5);plot(x(ix),avg(ix)+threshold*dev(ix),'LineWidth',1,'Color','green','LineWidth',1.5);plot(x(ix),avg(ix)-threshold*dev(ix),'LineWidth',1,'Color','green','LineWidth',1.5);plot(1:length(y),y,'b');subplot(2,1,2);stairs(signals,'r','LineWidth',1.5); ylim([-1.5 1.5]);R代码ThresholdingAlgosignalsfilteredYavgFilterstdFilteravgFilter[lag]stdFilter[lag]for (i in (lag+1):length(y)){if (abs(y[i]-avgFilter[i-1]) > threshold*stdFilter[i-1]) {if (y[i] > avgFilter[i-1]) {signals[i]} else {signals[i]

filteredY(i) = y(i);endavgFilter(i) = mean(filteredY(i-lag:i));stdFilter(i) = std(filteredY(i-lag:i));endend% Demo screen functionfunction [] = DemoScreen(n,lag,threshold,influence)figure('Position',[200 100,1000,500]);subplot(2,1,1);title(sprintf(['Draw data points (%.0f max) [settings: lag = %.0f, '...'threshold = %.2f, influence = %.2f]'],n,lag,threshold,influence));ylim([0 5]); xlim([0 50]);H = gca; subplot(2,1,1);set(H, 'YLimMode', 'manual'); set(H, 'XLimMode', 'manual');set(H, 'YLim', get(H,'YLim')); set(H, 'XLim', get(H,'XLim'));xg = []; yg = [];for i=1:ntry[xi,yi] = ginput(1);catchreturn;endxg = [xg xi]; yg = [yg yi];if i == 1subplot(2,1,1); hold on;plot(H, xg(i),yg(i),'r.');text(xg(i),yg(i),num2str(i),'FontSize',7);endif length(xg) > lag[signals,avg,dev] = ...ThresholdingAlgo(yg,lag,threshold,influence);area(xg(lag+1:end),avg(lag+1:end)+threshold*dev(lag+1:end),...'FaceColor',[0.9 0.9 0.9],'EdgeColor','none');

std::vector smoothedZScore(std::vector input){//lag 5 for the smoothing functionsint lag = 5;//3.5 standard deviations for signalfloat threshold = 3.5;//between 0 and 1, where 1 is normal influence, 0.5 is halffloat influence = .5;if (() <= lag + 2){std::vector emptyVec;return emptyVec;}//Initialise variablesstd::vector signals((), 0.0);std::vector filteredY((), 0.0);std::vector avgFilter((), 0.0);std::vector stdFilter((), 0.0);std::vector subVecStart((), () + lag);avgFilter[lag] = mean(subVecStart);stdFilter[lag] = stdDev(subVecStart);for (size_t i = lag + 1; i < (); i++){if (std::abs(input[i] - avgFilter[i - 1]) > threshold * stdFilter[i - 1]){if (input[i] > avgFilter[i - 1]){signals[i] = 1; //# Positive signal}else{signals[i] = -1; //# Negative signal}//Make influence lower

filteredY[i] = influence* input[i] + (1 - influence) * filteredY[i - 1];}else{signals[i] = 0; //# No signalfilteredY[i] = input[i];}//Adjust the filtersstd::vector subVec(() + i - lag, () + i);avgFilter[i] = mean(subVec);stdFilter[i] = stdDev(subVec);}return signals;}brad answered 2019-05-20T02:30:31Z3 votesC ++实现#include#include#include#include#include#include#includeusing namespace std;typedef long double ld;typedef unsigned int uint;typedef std::vector::iterator vec_iter_ld;/*** Overriding the ostream operator for pretty printing vectors.*/templatestd::ostream &operator< vec) {os << "[";

if (() != 0) {std::copy((), () - 1, std::ostream_iterator(os, " "));os << ();}os << "]";return os;}/*** This class calculates mean and standard deviation of a subvector.* This is basically stats computation of a subvector of a window size qual to "lag".*/class VectorStats {public:/*** Constructor for VectorStats class.** @param start - This is the iterator position of the start of the window,* @param end - This is the iterator position of the end of the window,*/VectorStats(vec_iter_ld start, vec_iter_ld end) {this->start = start;this->end = end;this->compute();}/*** This method calculates the mean and standard deviation using STL function.* This is the Two-Pass implementation of the Mean & Variance calculation.*/void compute() {ld sum = std::accumulate(start, end, 0.0);uint slice_size = std::distance(start, end);ld mean = sum / slice_size;std::vector diff(slice_size);std::transform(start, end, (), [mean](ld x) { return x - mean; });

ld sq_sum = std::inner_product((), (), (), 0.0);ld std_dev = std::sqrt(sq_sum / slice_size);this->m1 = mean;this->m2 = std_dev;}ld mean() {return m1;}ld standard_deviation() {return m2;}private:vec_iter_ld start;vec_iter_ld end;ld m1;ld m2;};/*** This is the implementation of the Smoothed Z-Score Algorithm.** @param input - input signal* @param lag - the lag of the moving window* @param threshold - the z-score at which the algorithm signals* @param influence - the influence (between 0 and 1) of new signals on the mean and standard deviation* @return a hashmap containing the filtered signal and corresponding mean and standard deviation.*/unordered_map> z_score_thresholding(vector input, int lag, ld threshold, ld influence) {unordered_map> output;uint n = (uint) ();vector signals(());vector filtered_input((), ());vector filtered_mean(());vector filtered_stddev(());VectorStats lag_subvector_stats((), () + lag);

* Uses a rolling mean and a rolling deviation (separate) to identify peaks in a vector** @param y - The input vector to analyze* @param lag - The lag of the moving window (i.e. how big the window is)* @param threshold - The z-score at which the algorithm signals (i.e. how many standard deviations away from the movingmean a peak (or signal) is)* @param influence - The influence (between 0 and 1) of new signals on the mean and standard deviation (how much a peak(or signal) should affect other values near it)* @return - The calculated averages (avgFilter) and deviations (stdFilter), and the signals (signals)*/public HashMap> thresholdingAlgo(List y, Long lag, Double threshold, Double influence) {//init stats instanceSummaryStatistics stats = new SummaryStatistics()//the results (peaks, 1 or -1) of our algorithmList signals = new ArrayList(s((), 0))//filter out the signals (peaks) from our original list (using influence arg)List filteredY = new ArrayList(y)//the current average of the rolling windowList avgFilter = new ArrayList(s((), 0.0d))//the current standard deviation of the rolling windowList stdFilter = new ArrayList(s((), 0.0d))//init avgFilter and stdFilter(0..lag-1).each { ue(y[it as int]) }avgFilter[lag - 1 as int] = n()stdFilter[lag - 1 as int] = (ulationVariance()) //getStandardDeviation() uses sample variance (not whatwe want)()//loop input starting at end of rolling window(lag..()-1).each { i ->//if the distance between the current value and average is enough standard deviations (threshold) awayif (((y[i as int] - avgFilter[i - 1 as int]) as Double) > threshold * stdFilter[i - 1 as int]) {//this is a signal (i.e. peak), determine if it is a positive or negative signalsignals[i as int] = (y[i as int] > avgFilter[i - 1 as int]) ? 1 : -1//filter this signal out using influencefilteredY[i as int] = (influence * y[i as int]) + ((1-influence) * filteredY[i - 1 as int])

case (s: Double, i: Int) =>// if the distance between the current value and average is enough standard deviations (threshold) awayif ((s - avgFilter(i - 1)) > threshold * stdFilter(i - 1)) {// this is a signal (i.e. peak), determine if it is a positive or negative signalsignals(i) = if (s > avgFilter(i - 1)) 1 else -1// filter this signal out using influencefilteredY(i) = (influence * s) + ((1 - influence) * filteredY(i - 1))} else {// ensure this signal remains a zerosignals(i) = 0// ensure this value is not filteredfilteredY(i) = s}// update rolling average and ()(i - lag, i).foreach(s => ue(s))avgFilter(i) = nstdFilter(i) = (ulationVariance) // getStandardDeviation() uses sample variance (not what we want)}println()println()println(signals)h {case(x: Int, idx: Int) =>if (x == 1) {println(idx + " " + y(idx))}}val data = { case (s: Double, i: Int) => Map("x" -> i, "y" -> s, "name" -> "y", "row" -> "data") } ++ { case (s: Double, i: Int) => Map("x" -> i, "y" -> s, "name" -> "avgFilter", "row" -> "data") } ++ { case (s: Double, i: Int) => Map("x" -> i, "y" -> (s - threshold * stdFilter(i)), "name" -> "lower", "row" -> "data") } ++ { case (s: Double, i: Int) => Map("x" -> i, "y" -> (s + threshold * stdFilter(i)), "name" -> "upper", "row"-> "data") } ++

* @param influence - The influence (between 0 and 1) of new signals on the mean and standard deviation (how much a peak(or signal) should affect other values near it)* @return - The calculated averages (avgFilter) and deviations (stdFilter), and the signals (signals)*/fun smoothedZScore(y: List, lag: Int, threshold: Double, influence: Double): Triple, List, List> {val stats = SummaryStatistics()// the results (peaks, 1 or -1) of our algorithmval signals = MutableList(, { 0 })// filter out the signals (peaks) from our original list (using influence arg)val filteredY = ArrayList(y)// the current average of the rolling windowval avgFilter = MutableList(, { 0.0 })// the current standard deviation of the rolling windowval stdFilter = MutableList(, { 0.0 })// init avgFilter and (lag).forEach { s -> ue(s) }avgFilter[lag - 1] = dFilter[lag - 1] = (tionVariance) // getStandardDeviation() uses sample variance (not what we want)()//loop input starting at end of rolling window(lag.. - 1).forEach { i ->//if the distance between the current value and average is enough standard deviations (threshold) awayif ((y[i] - avgFilter[i - 1]) > threshold * stdFilter[i - 1]) {//this is a signal (i.e. peak), determine if it is a positive or negative signalsignals[i] = if (y[i] > avgFilter[i - 1]) 1 else -1//filter this signal out using influencefilteredY[i] = (influence * y[i]) + ((1 - influence) * filteredY[i - 1])} else {//ensure this signal remains a zerosignals[i] = 0//ensure this value is not filteredfilteredY[i] = y[i]}//update rolling average and deviation

integer, dimension(size(y)) :: PeakDetectreal, dimension(size(y)) :: filteredY, avgFilter, stdFilterinteger :: lag, iireal :: threshold, influence! Executing partPeakDetect = 0filteredY = 0.0filteredY(1:lag+1) = y(1:lag+1)avgFilter = 0.0avgFilter(lag+1) = mean(y(1:2*lag+1))stdFilter = 0.0stdFilter(lag+1) = std(y(1:2*lag+1))if (stdFilter(lag+1)/avgFilter(lag+1)>0.1) then ! If the coefficient of variation exceeds 10%, the signal is too uneven at thestart, possibly because of a (unit=*,fmt=1001)1001 format(1X,'Warning: Peak detection might have failed, as there may be a peak at the edge of the frequency range.',/)end ifdo ii = lag+2, size(y)if (abs(y(ii) - avgFilter(ii-1)) > threshold * stdFilter(ii-1)) then! Find only the largest outstanding value which is only the one greater than its predecessor and its successorif (y(ii) > avgFilter(ii-1) .AND. y(ii) > y(ii-1) .AND. y(ii) > y(ii+1)) thenPeakDetect(ii) = 1end iffilteredY(ii) = influence * y(ii) + (1 - influence) * filteredY(ii-1)elsefilteredY(ii) = y(ii)end if! Modified with respect to the original code. Mean and standard deviation are calculted symmetrically around the currentpointavgFilter(ii) = mean(filteredY(ii-lag:ii+lag))stdFilter(ii) = std(filteredY(ii-lag:ii+lag))end doend function PeakDetectreal function mean(y)!> @brief Calculates the mean of vector y

signals[i] += 1 # postive signalelsesignals[i] += -1 # negative signalend# Make influence lowerfilteredY[i] = influence*y[i] + (1-influence)*filteredY[i-1]elsesignals[i] = 0filteredY[i] = y[i]endavgFilter[i] = mean(filteredY[i-lag+1:i])stdFilter[i] = std(filteredY[i-lag+1:i])endreturn (signals = signals, avgFilter = avgFilter, stdFilter = stdFilter)end# Datay = [1,1,1.1,1,0.9,1,1,1.1,1,0.9,1,1.1,1,1,0.9,1,1,1.1,1,1,1,1,1.1,0.9,1,1.1,1,1,0.9,1,1.1,1,1,1.1,1,0.8,0.9,1,1.2,0.9,1,1,1.1,1.2,1,1.5,1,3,2,5,3,2,1,1,1,0.9,1,1,3,2.6,4,3,3.2,2,1,1,0.8,4,4,2,2.5,1,1,1]# Settings: lag = 30, threshold = 5, influence = 0lag = 30threshold = 5influence = 0results = SmoothedZscoreAlgo(y, lag, threshold, influence)upper_bound = results[:avgFilter] + threshold * results[:stdFilter]lower_bound = results[:avgFilter] - threshold * results[:stdFilter]x = 1:length(y)yplot = plot(x,y,color="blue", label="Y",legend=:topleft)yplot = plot!(x,upper_bound, color="green", label="Upper Bound",legend=:topleft)yplot = plot!(x,results[:avgFilter], color="cyan", label="Average Filter",legend=:topleft)yplot = plot!(x,lower_bound, color="green", label="Lower Bound",legend=:topleft)signalplot = plot(x,results[:signals],color="red",label="Signals",legend=:topleft)plot(yplot,signalplot,layout=(2,1),legend=:topleft)Matt Camp answered 2019-05-20T02:35:30Z