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 scalingfindlevels/D=risingX/EDGE=1 target,5findlevels/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 waveVariable nRising = numpnts(risingX)Variable nFalling = numpnts(fallingX)if( nRising > nFalling )// difference should be at most 1InsertPoints nFalling, 1, fallingX
fallingX[nFalling]= rightx(target)endif// debugging for Graph0Duplicate/O risingX, risingY; risingY=target(risingX[p])// risingY should be approx thresholdDuplicate/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 & fallingXEnd
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 510FUNCTION 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 anotherExtract/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 eventsMatrixOp/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 eventsExtract/FREE/INDX/O wStart_Switch_temp, wStart_Dex_temp, wStart_Switch_temp!=1//Makes a wave of the start indicesMake/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 aboveMatrixOp/O/FREE wStop_Switch_temp=shiftVector(wIndices, -1,0)-wIndices
Extract/FREE/INDX/O wStop_Switch_temp, wStop_Dex_temp, wStop_Switch_temp!=1Make/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 anotherExtract/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 eventsMatrixOp/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 eventsExtract/FREE/INDX/O wStart_Switch_temp, wStart_Dex_temp, wStart_Switch_temp!=1//Makes a wave of the start indicesMake/O/D/N=(numpnts(wStart_Dex_temp)) wStart_dex_2=wIndices[wStart_Dex_temp]//Makes the wave for the stop indicesMake/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 pointsMake/O/N=(1e5) DataWave=gnoise(1)
DataWave[1,*]=0.999*DataWave[p-1]+DataWave[p]// The threshold value to seach forVariable Threshold=20// Calculates for each value in the DataWave if the value is above the threshold// This will create a wave looking like: 0001111001100101Duplicate/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: 0001234001200101Duplicate/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 valueVariable MaxDuration=WaveMax(ThresholdDurationWave)Make/O/N=(MaxDuration+1) HistogramWave
Histogram/B=2 ThresholdDurationWave, HistogramWave
// Displays all the wavesDoWindow/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 dominateSetAxis/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:
For simple cases, e.g., requiring 3 consecutive points above threshold:
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.
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
--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.
//****Method 2****
END
September 18, 2017 at 03:04 pm - Permalink
September 19, 2017 at 06:35 am - Permalink