Automatic indexing of individual cycles in continuous sinusoidal data

zephyr070
Posts: 7
Joined: 2017-05-23
Location: United States

I am revisiting a problem I had to shelve last July. Our lab has an instrument that acquires sinusoidal position data from an actuator at a rate of 50 samples/s, with an oscillation frequency of 1 Hz. I would like to index the minima on either side of each cycle (hundreds of thousands) by row location, essentially giving me the "start" and "stop" rows of the cycle. This is time-prohibitive to do manually using GetInfo cursors, but it's how I've been doing it. The intent is to use those start and stop row values to generate graphical analysis of an arbitrary cycle or subset of cycles out of the whole 100k+ set.

To complicate matters, with much staring at previous data sets, I am aware there is a ca. 11.5 microsecond/cycle error in the position signal coming from our instrument (permanent firmware issues), so the start and stop times of the nth cycle are not where you'd expect them to be based on those of the first cycle. I need to create a loop routine that accounts for this indeterminacy, and I have come up with the following strategy, but am having trouble coding it:

1) Make wave to store (Cycle#, StartRow, StopRow)
2) Set dimension labels according to step 1
3) Define data source wave
4) Define input range variables to use in WaveStats as [Start, Stop], these will need to be incremented by 50 rows after each loop iteration
5) Loop:
-Find V_minRowLoc between first and second cycles using WaveStats, with initial range values determined by manually placed cursors, call this point P
-Assign P as first StopRow value
-First StartRow value = (First StopRow - 50)
-Increment WaveStats range values to encompass next minimum, using (P-25) for Start and (P+25) for Stop

This seems like it might be overly complicated, and perhaps a more elegant solution exists. I do know I eventually require the precision of WaveStats, rather than FindPeak. Below is what I have so far, based on a gracious response from Andy Hegedus when I attempted to tackle this problem last time:

#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3		// Use modern global access method and strict wave access.
 
Function AutoIndex()
 
	//Create Wave to store results
	Make/O /N=(0,3) Limits
 
	setdimlabel 1,0,Cycle,Limits
	setdimlabel 1,1,StartRow,Limits
	setdimlabel 1,2,StopRow,Limits
 
	//Define sources of waves
	Wave Position
 
	//define the range of inputs
	SVAR CycleList
	Variable cycle, maxcycle
	Variable start, stop
	maxcycle = itemsinlist(CycleList)// allow for indeterminent number of cycles
 
	For(cycle=0;cycle<maxcycle;cycle+=1)
		WaveStats/R=[start,stop] Position
		Limits[cycle][%Cycle]= {cycle+1}
		Limits[cycle][%StopRow] = {V_minRowLoc}
		Limits[cycle][%StartRow] = {(StopRow-50)}
        //more code required
 
	Endfor
end


tony
Posts: 168
Joined: 2007-08-28
Location: Canada

If the signal is sinusoidal, at least locally, maybe you could fit a sine (perhaps sequentially over suitable subranges centered about the peak) to get the start/end of cycle positions? Each fit should give you good starting guesses for coefficients for the next.


hegedus
Posts: 327
Joined: 2009-03-21
Location: United States

Another scheme if the data has fairly good precision (i.e. not very much noise).

Take a derivative and then Findlevels /edge=2 derivativeWave 0

The resulting `W_FindLevels (or your own wave name if you use the /Dest flag) will give you all the maxima positions.

Andy


[ last edited May 14, 2018 - 16:54 ]
thomas_braun
Posts: 586
Joined: 2009-10-07
Location: Germany

@zephyr070 Could you attach some small sample data?


zephyr070
Posts: 7
Joined: 2017-05-23
Location: United States

Thank you all for your helpful comments. I have revised my procedure as follows, and it at least compiles without any errors:

#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3		// Use modern global access method and strict wave access.
 
Function IndexCycles()
 
	//Create Wave to store results
	Make/O /N=(0,3) Limits
 
	setdimlabel 1,0,Cycle,Limits
	setdimlabel 1,1,StartRow,Limits
	setdimlabel 1,2,StopRow,Limits
 
	//Define sources of waves
	Wave Position
 
	//define the range of inputs
	Variable cycle
	Variable start, stop
	Variable a=pcsr(A),b=pcsr(B)
	Variable V_minRowLoc
 
	For(cycle=0;cycle<10;cycle+=1)//max cycle count has to be entered manually
		If(cycle>0)
			a=V_minRowLoc+35
			b=V_minRowLoc+65
		endIf
			WaveStats/R=[start,stop] Position
			Limits[cycle][%Cycle]= {cycle+1} //make it 1 based
			Limits[cycle][%StopRow] = {V_minRowLoc} // include all increase more than one
			Limits[cycle][%StartRow] = {V_minRowLoc-50} // include all increase more than one
 
 
	Endfor
	Edit Limits
end

I'm sure I'll show my newbie colors with this question, but how do I execute the .ipf I've just complied from the command window, without copying and pasting all the code?

Per Thomas Braun's request, I have attached a sample experiment with 10 cycles.

AttachmentSize
Sinusoidal_data_10cycles.pxp18.42 KB

hegedus
Posts: 327
Joined: 2009-03-21
Location: United States

Just input

indexcycles()

Andy


hegedus
Posts: 327
Joined: 2009-03-21
Location: United States

Here is what I did with your problem:

1. Calculate a derivative: Differentiate/METH=1 Position/X=t_sec/D=Position_DIF//History label from menu item differentiate

2. Find where it crosses zero findLevels/P /Edge=2 Position_dif 0 //the edge flag says only consider crossings that are decreasing (maxima)

Igor replies: V_Flag= 1; V_LevelsFound= 10; I use point scaling so now you have a wave with 10 values on the maxima.

If you want the were it crosses zero as the marker in the cycle, you can take the second derivative and do the find levels at 0.

In the picture I created some additional waves to hold the Y position.

Andy

AttachmentSize
Example Method.png119.69 KB

_sk
Posts: 115
Joined: 2014-10-28
Location: Switzerland

You can get all the indices where a certain threshold value is crossed. Using some heuristics one can write something like:
make/o/n=881 w_max = (position[p] > 1197) ? p : nan

Now, w_max stores all indices where the sin wave is greater than the heuristic threshold value, 1197.

If you want you can zapnans the wave like so:
wavetransform/o zapnans w_max

Now, w_max will store only the indices without the nan.

Since the number of points of position and t_sec are the same, you can write the time at which these extrema occur in another wave, like so:
make/o/n=881 w_max_t = (position[p] > 1197) ? t_sec[p] : nan; wavetransform/o zapnans w_max_t;

And there really are many variations on the theme.

best,
_sk


[ last edited May 16, 2018 - 01:37 ]
zephyr070
Posts: 7
Joined: 2017-05-23
Location: United States

Thank you all for the input. Andy, I tried your method of differentiating and then using FindLevels on my most recent data set, which comprises 1,000,000 cycles. I Set /EDGE=1 so I could obtain the minima zero crossings. Unexpectedly, the operation returned V_LevelsFound = 1,000,010. I decided to check this by performing FindLevels on the position data (not the derivative) with /EDGE=2 and a level equal to half the amplitude and obtained V_LevelsFound = 1,000,000 as expected. I'm not sure why there is a discrepancy there.

The following updated function works well for #cycles up to several hundred thousand, but analysis of the derived point values for StartRow and StopRow show the data appear to be off by 2 whole cycles at the end of the execution: I counted backwards from the last cycle visually, then placed cursors according the StartRow and StopRow for that cycle. I added the move cursor commands so I could watch the cursors step through the entire data set as sort of a pseudo progress bar, to make sure they weren't "hopping over" a cycle by mistake, and they are always somewhere around half amplitude for each cycle step.

Two additional observations:
1) Often in the latter cycles, the difference between StartRow and StopRow is not 50, as is supposed to be according to my code. What could be causing that discrepancy? I thought maybe it was due to the type of numerical data being output, which according to WaveType is 64-bit floating point, but I am stumped.
2) A visual inspection of the graph at latter cycles shows V_minRowLoc is incorrect, and off by a couple points. How is WaveStats missing this?

#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3		// Use modern global access method and strict wave access.
 
Function IndexCycles()
 
	//Create Wave to store results
	Make/O /N=(0,5) Limits
 
	setdimlabel 1,0,Cycle,Limits
	setdimlabel 1,1,StartRow,Limits
	setdimlabel 1,2,StopRow,Limits
	setdimlabel 1,3,CsrAPos,Limits
	setdimlabel 1,4,CsrBPos,Limits
 
	//Define sources of waves
	Wave Pos
 
	//define the range of inputs
	Variable cycle
	Variable start=pcsr(A),stop=pcsr(B)
	Variable V_minRowLoc
 
 	//max cycle count has to be entered manually, = (#cycles in data set - 1)
	For(cycle=0;cycle<999999;cycle+=1)
		If(cycle>0)
			start=V_minRowLoc+37
			stop=V_minRowLoc+63
			Cursor A Pos start
			Cursor B Pos stop
		endIf
		WaveStats/Q/R=[start,stop] Pos
		Limits[cycle][%Cycle] = {cycle+1}
		Limits[cycle][%StopRow] = {V_minRowLoc}
		Limits[cycle][%StartRow] = {V_minRowLoc-50}
		Limits[cycle][%CsrAPos] = {pcsr(A)}
		Limits[cycle][%CsrBPos] = {pcsr(B)}
	Endfor
 
	//last cycle does not obey If-endif conditions, enter limits manually
	WaveStats/R=[V_minRowLoc,(V_minRowLoc+50)] Pos
	Edit Limits.ld
	InsertPoints 0,1,Limits
 
end


Back to top