Crossing threshold for a certain duration

Share icon

Tue, 09/12/2017 - 01:00 pm

I have waves in which I'd like to find out if a threshold is crossed continuously for a minimum of 3 seconds and where the x-location of the final threshold crossing is (ie. determining the time during which the values above threshold are maintained). I am having some difficulty implementing this using Findlevels operation with this and would appreciate any advice. Am I perhaps using the wrong function?
I'm sure you could simply write a procedure to perform this task but I find it entertaining to think how this could be done using MatrixOP. My approach would be the following:
MatrixOP/O ot=greater(srcWave,thresh).

So far you are just generating a wave that gives you 1 where data are over the specified thresh. Next you need to look at your sampling and figure out how many consecutive points you need to maintain (we will call that N). Here you need to accept that you will not be able to make any judgement about the beginning and the end of the wave (N-1 points). Next consider the product:
 
MatrixOP/O ot=ot*rotateRows(ot,1).

For simple cases, e.g., requiring 3 consecutive points above threshold:
 
MatrixOP/O ot=ot*rotateRows(ot,1)*rotateRows(ot,2),
etc.

I hope this helps,

A.G.
WaveMetrics, Inc.

September 12, 2017 at 01:57 pm - Permalink

I would expect that FindLevels would be pretty useful for this application.

Maybe if you posted your code and an example of the data as an experiment file we could help you better.

--Jim Prouty
Software Engineer, WaveMetrics, Inc.

September 12, 2017 at 08:45 pm - Permalink

This is what I have so far. I'm checking if two points are less than 3 seconds apart as a test but it's probably not foolproof.
Any suggestions for a better way to continuously test if the threshold is maintained for 3 seconds after the first point above threshold is detected.

function tryit()
wave tempstorage
wave target
variable duration
variable i=0
variable x1 = 100
variable x2 = 1000

Findlevels /D=tempstorage /Edge=1 /Q /R= (x1,x2) target,5


    for(i=0;i<numpnts(tempstorage)-1;i+=1) 
        if ((tempstorage[i+1] - tempstorage[i] < 3) && (tempstorage[i+2] - tempstorage[i+1] < 3))
            duration =  tempstorage[i+2] - tempstorage[0]  
       
           break
           
        endif          
    endfor                     




print duration

end
   

September 13, 2017 at 12:00 pm - Permalink

You have your starting point. For the duration criteria, integrate the wave and look for the crossing of Intensity*Time > (Threshold*3).

--
J. J. Weimer
Chemistry / Chemical & Materials Engineering, UAH

September 13, 2017 at 03:58 pm - Permalink

Something like this?

(Experiment attached)

--Jim

#pragma TextEncoding = "Windows-1252"
#pragma rtGlobals=3     // Use modern global access method and strict wave access.

Macro tryit()

    Make/O/N=1000 target= sin(x*x/2500)*10 // point scaling
    findlevels/D=risingX/EDGE=1 target,5
    findlevels/D=fallingX/EDGE=2 target,5
    // There are two cases: rising edge is earliest, or falling edge is earliest.
    // Since we want above the threshold, we discard an earlier falling edge.
    // This code assumes X is increasing with point number.
    Variable firstRisingX= risingX[0]
    Variable firstFallingX= fallingX[0]
    if( firstFallingX < firstRisingX )
        DeletePoints 0, 1, fallingX
    endif
   
    // if a rising edge has no corresponding falling edge, add one at the end of the wave
    Variable nRising = numpnts(risingX)
    Variable nFalling = numpnts(fallingX)
    if( nRising > nFalling ) // difference should be at most 1
        InsertPoints nFalling, 1, fallingX
        fallingX[nFalling]= rightx(target)
    endif
   
    // debugging for Graph0
    Duplicate/O risingX, risingY; risingY=target(risingX[p]) // risingY should be approx threshold
    Duplicate/O fallingX, fallingY; fallingY=target(fallingX[p]) // ditto
   
    // compute the duration above the threshold as falling-rising times.
    Duplicate/O risingX, duration
    duration = fallingX - risingX // duration above threshold

    // create list of durations (by index) that are >=3.
    Extract/INDX duration, indexOfDurationsLongEnough, duration>=3
   
    // the indexOfDurationsLongEnough wave contains only the indexes of the durations that are >= 3.0.
    // durations shorter than that are not included, so it may well be shorter than risingX & fallingX
End

--Jim Prouty
Software Engineer, WaveMetrics, Inc.
DurationAboveThreshold.pxp

September 13, 2017 at 03:59 pm - Permalink

I had a similar situation where I needed to find the start and stop indices of events based on the values in a state wave. I came up with a method similar to the MatrixOp version above to get the start and stop indices. Subtracting the stop and start to get the length and pulling out the events that are greater than a certain duration would be easy to add. I haven't had any issues so far, and it handles events that are only one point long or are at either the start or stop of the state wave. I came up with two ways of finding the stop indices. The first is a little slower than the second, but seems a little more robust. I tried a couple of methods using FindLevels, and the Extract/MatrixOP version ended up being faster.

EDIT: Using shiftVector(wIndices, 1,0) to get the start indices causes the code to miss events that start at point p=1 (second point in the wave). Changing the fill value to NaN seems to solve this.

#pragma rtGlobals=3     // Use modern global access method and strict wave access.

//Gets the start and stop indices for all events where the state is equal to 510
FUNCTION Get_Start_Stop_Indices()

    Wave/Z wState=root:wState
       
    //****Method 1****
    //Grabs the indices where wState==510.  The values in this wave will increment by 1 during events, and increment by a larger number when going from one event to another
    Extract/FREE/INDX/O wState, wIndices, wState==510                                                  
   
    //Get the start indices
    //Makes a shifted wave and subtracts that from the unshifted wave.  The resulting wave will be 1 during an event and greater than one when switching between events
    MatrixOp/O/FREE wStart_Switch_temp=wIndices-shiftVector(wIndices, 1,NaN)
    //Gets the indices where the differences are not equal to 1.  These are the start indices of the events
    Extract/FREE/INDX/O wStart_Switch_temp, wStart_Dex_temp, wStart_Switch_temp!=1
    //Makes a wave of the start indices
    Make/O/D/N=(numpnts(wStart_Dex_temp)) wStart_dex_1=wIndices[wStart_Dex_temp]       

    //Get the stop indices
    //Notice it's shifted-unshifted here, and the reverse above
    MatrixOp/O/FREE wStop_Switch_temp=shiftVector(wIndices, -1,0)-wIndices                     
    Extract/FREE/INDX/O wStop_Switch_temp, wStop_Dex_temp, wStop_Switch_temp!=1
    Make/O/D/N=(numpnts(wStop_Dex_temp)) wStop_dex_1=wIndices[wStop_Dex_temp]
     //****Method 1****
     
    //***********************************
           
    //****Method 2****
    //Grabs the indices where wState==510.  The values in this wave will increment by 1 during events, and increment by a larger number when going from one event to another
    Extract/FREE/INDX/O wState, wIndices, wState==510
   
    //Get the start indices
    //Makes a shifted wave and subtracts that from the unshifted wave.  The resulting wave will be 1 during an event and greater than one when switching between events
    MatrixOp/O/FREE wStart_Switch_temp=wIndices-shiftVector(wIndices, 1,NaN)
    //Finds the indices where the differences are not equal to 1.  These are the start indices of the events
    Extract/FREE/INDX/O wStart_Switch_temp, wStart_Dex_temp, wStart_Switch_temp!=1
    //Makes a wave of the start indices
    Make/O/D/N=(numpnts(wStart_Dex_temp)) wStart_dex_2=wIndices[wStart_Dex_temp]
   
    //Makes the wave for the stop indices
    Make/O/D/N=(numpnts(wStart_Dex_temp)) wStop_dex_2
    //The last stop index is the last point in wIndices
    wStop_dex_2[numpnts(wStop_dex_2)-1]=wIndices[numpnts(wIndices)-1]
    //The rest of the points are the points right before the next start index          
    IF(numpnts(wStop_dex_2)>1)
        wStop_dex_2[0,numpnts(wStart_Dex_temp)-2]=wIndices[wStart_Dex_temp[p+1]-1]
   
IF
//****Method 2****

END
Find_State_Start_Stop.pxp

September 18, 2017 at 03:04 pm - Permalink

This would be my approach to the problem. It creates a histogram with the number of values above a threshold in sequence. The number of times a value was greater than the threshold for any given number of points (I used points instead of seconds to make it easier) can then be read from the resulting graph.

Function Test()
//  Test Function

    //  Creates a realistic-looking wave with fake data points
    Make/O/N=(1e5) DataWave=gnoise(1)
    DataWave[1,*]=0.999*DataWave[p-1]+DataWave[p]

    //  The threshold value to seach for
    Variable Threshold=20
   
    //  Calculates for each value in the DataWave if the value is above the threshold
    //  This will create a wave looking like: 0001111001100101
    Duplicate/O DataWave ThresholdWave
    ThresholdWave[]=DataWave[p]>Threshold
   
    //  Creates a threshold duration wave which counts up by one for each consecutive data point above the threshold
    //  For the above example this would create a wave looking like: 0001234001200101
    Duplicate/O ThresholdWave ThresholdDurationWave
    ThresholdDurationWave[0]=ThresholdWave[0]
    ThresholdDurationWave[1,*]=(ThresholdDurationWave[p-1]+ThresholdWave[p])*ThresholdWave[p]
   
    //  Creates a histogram wave with the number of events in which a given number or more of consecutive data points are above the threshold value
    Variable MaxDuration=WaveMax(ThresholdDurationWave)
    Make/O/N=(MaxDuration+1) HistogramWave
    Histogram /B=2 ThresholdDurationWave, HistogramWave
   
    //  Displays all the waves
    DoWindow/K DataWindow
    Display /W=(100,100,500,300) /N=DataWindow /K=1 DataWave

    DoWindow/K ThresholdWindow
    Display /W=(100,330,900,430) /N=ThresholdWindow /K=1 ThresholdWave

    DoWindow/K ThresholdDurationWindow
    Display /W=(100,460,900,560) /N=ThresholdDurationWindow /K=1 ThresholdDurationWave

    DoWindow/K HistogramWindow
    Display /W=(510,100,910,300) /N=HistogramWindow /K=1 HistogramWave
   
    //  Cuts off the event with no data points in sequence above the threshold which would otherwise dominate
    SetAxis /W=HistogramWindow bottom 1, MaxDuration
    SetAxis/A=2 left
end

September 19, 2017 at 06:35 am - Permalink