Crossing threshold for a certain duration

Share icon

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

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 (16.25 KB)

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 (88.06 KB)

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