## Crossing threshold for a certain duration

epiphenom
Posts: 77
Joined: 2012-08-21

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?

Posts: 756
Joined: 2007-06-29
Location: United States

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.

Posts: 619
Joined: 2007-07-19
Location: United States

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.

epiphenom
Posts: 77
Joined: 2012-08-21

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

jjweimer
Posts: 1262
Joined: 2007-08-14
Location: United States

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

Posts: 619
Joined: 2007-07-19
Location: United States

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.

AttachmentSize
DurationAboveThreshold.pxp16.25 KB

[ last edited September 13, 2017 - 15:02 ]
KZarzana
Posts: 1
Joined: 2017-09-13
Location: United States

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]
ENDIF
//****Method 2****

END

AttachmentSize
Find_State_Start_Stop.pxp88.06 KB

[ last edited September 18, 2017 - 14:04 ]
olelytken
Posts: 54
Joined: 2015-10-12
Location: Germany

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