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:
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.
#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 ) DeletePoints0, 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
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]
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
// 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
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),
I hope this helps,
A.G.
WaveMetrics, Inc.
September 12, 2017 at 01:57 pm - Permalink
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
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.
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
--
J. J. Weimer
Chemistry / Chemical & Materials Engineering, UAH
September 13, 2017 at 03:58 pm - Permalink
(Experiment attached)
--Jim
#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.
September 13, 2017 at 03:59 pm - Permalink
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.
//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]
//****Method 2****
END
September 18, 2017 at 03:04 pm - Permalink
// 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