[algorithm] Peak signal detection in realtime timeseries data

Here is an altered Fortran version of the z-score algorithm. It is altered specifically for peak (resonance) detection in transfer functions in frequency space (Each change has a small comment in code).

The first modification gives a warning to the user if there is a resonance near the lower bound of the input vector, indicated by a standard deviation higher than a certain threshold (10% in this case). This simply means the signal is not flat enough for the detection initializing the filters properly.

The second modification is that only the highest value of a peak is added to the found peaks. This is reached by comparing each found peak value to the magnitude of its (lag) predecessors and its (lag) successors.

The third change is to respect that resonance peaks usually show some form of symmetry around the resonance frequency. So it is natural to calculate the mean and std symmetrically around the current data point (rather than just for the predecessors). This results in a better peak detection behavior.

The modifications have the effect that the whole signal has to be known to the function beforehand which is the usual case for resonance detection (something like the Matlab Example of Jean-Paul where the data points are generated on the fly won't work).

function PeakDetect(y,lag,threshold, influence)
    implicit none
    ! Declaring part
    real, dimension(:), intent(in) :: y
    integer, dimension(size(y)) :: PeakDetect
    real, dimension(size(y)) :: filteredY, avgFilter, stdFilter
    integer :: lag, ii
    real :: threshold, influence

    ! Executing part
    PeakDetect = 0
    filteredY = 0.0
    filteredY(1:lag+1) = y(1:lag+1)
    avgFilter = 0.0
    avgFilter(lag+1) = mean(y(1:2*lag+1))
    stdFilter = 0.0
    stdFilter(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 the start, possibly because of a peak.
        write(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 if
    do 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 successor
            if (y(ii) > avgFilter(ii-1) .AND. y(ii) > y(ii-1) .AND. y(ii) > y(ii+1)) then
                PeakDetect(ii) = 1
            end if
            filteredY(ii) = influence * y(ii) + (1 - influence) * filteredY(ii-1)
        else
            filteredY(ii) = y(ii)
        end if
        ! Modified with respect to the original code. Mean and standard deviation are calculted symmetrically around the current point
        avgFilter(ii) = mean(filteredY(ii-lag:ii+lag))
        stdFilter(ii) = std(filteredY(ii-lag:ii+lag))
    end do
end function PeakDetect

real function mean(y)
    !> @brief Calculates the mean of vector y
    implicit none
    ! Declaring part
    real, dimension(:), intent(in) :: y
    integer :: N
    ! Executing part
    N = max(1,size(y))
    mean = sum(y)/N
end function mean

real function std(y)
    !> @brief Calculates the standard deviation of vector y
    implicit none
    ! Declaring part
    real, dimension(:), intent(in) :: y
    integer :: N
    ! Executing part
    N = max(1,size(y))
    std = sqrt((N*dot_product(y,y) - sum(y)**2) / (N*(N-1)))
end function std

For my application the algorithm works like a charm! enter image description here

Examples related to algorithm

How can I tell if an algorithm is efficient? Find the smallest positive integer that does not occur in a given sequence Efficiently getting all divisors of a given number Peak signal detection in realtime timeseries data What is the optimal algorithm for the game 2048? How can I sort a std::map first by value, then by key? Finding square root without using sqrt function? Fastest way to flatten / un-flatten nested JSON objects Mergesort with Python Find common substring between two strings

Examples related to language-agnostic

IOException: The process cannot access the file 'file path' because it is being used by another process Peak signal detection in realtime timeseries data Match linebreaks - \n or \r\n? Simple way to understand Encapsulation and Abstraction How can I pair socks from a pile efficiently? How do I determine whether my calculation of pi is accurate? What is ADT? (Abstract Data Type) How to explain callbacks in plain english? How are they different from calling one function from another function? Ukkonen's suffix tree algorithm in plain English Private vs Protected - Visibility Good-Practice Concern

Examples related to time-series

How to convert dataframe into time series? Peak signal detection in realtime timeseries data Pandas: rolling mean by time interval Excel plot time series frequency with continuous xaxis How to calculate rolling / moving average using NumPy / SciPy? Plotting two variables as lines using ggplot2 on the same graph

Examples related to signal-processing

Creating lowpass filter in SciPy - understanding methods and units Peak signal detection in realtime timeseries data How to smooth a curve in the right way? Plotting power spectrum in python How to implement band-pass Butterworth filter with Scipy.signal.butter How to normalize a signal to zero mean and unit variance? How do I obtain the frequencies of each value in an FFT? How to apply a low-pass or high-pass filter to an array in Matlab? An implementation of the fast Fourier transform (FFT) in C#

Examples related to data-analysis

Python: pandas merge multiple dataframes How do I sum values in a column that match a given condition using pandas? Peak signal detection in realtime timeseries data How to sort a dataFrame in python pandas by two or more columns? Fitting polynomial model to data in R