# Automatic indexing of individual cycles in continuous sinusoidal data

zephyr070

Mon, 05/14/2018 - 01:37 pm

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

#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

May 14, 2018 at 05:28 pm - Permalink

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

May 14, 2018 at 05:54 pm - Permalink

May 15, 2018 at 01:18 am - Permalink

#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.

May 15, 2018 at 12:09 pm - Permalink

indexcycles()

Andy

May 15, 2018 at 12:37 pm - Permalink

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

May 15, 2018 at 01:01 pm - Permalink

`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

May 16, 2018 at 02:37 am - Permalink

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 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

May 17, 2018 at 01:14 pm - Permalink

The procedure does throw a curve fit error at the end of execution: "Must contain at least as many points as there are parameters," or something to that effect, but the actual fit curve on the last cycle appears to be normal. Further, the row locations for the beginning and end of the cycle are correct.

#pragma rtGlobals=3 // Use modern global access method and strict wave access.

//Before executing program, perform FindLevels on data wave bounded by cursors, once for EDGE=1

//and again for EDGE=2, store these positions in Lvls_L and Lvls_R, respectively.

//Allocate #cycles

constant kNumCycle=1000000

Function IndexMax()

//Create wave to store results

Make/O/N=(kNumCycle,4) Limits

SetDimLabel 1,0,Cycle,Limits

SetDimLabel 1,1,StartRow,Limits

SetDimLabel 1,2,StopRow,Limits

SetDimLabel 1,3,Maxima,Limits

//Define source wave

Wave Pos,t_s

Wave Lvls_L,Lvls_R

//Define input variables, manually enter start and stop from cycle 0 of Lvls_L, Lvls_R

Variable cycle

Variable start=19194,stop=19210

Variable Xcross

Variable V_value

For(cycle=0;cycle<kNumCycle;cycle+=1)

If(cycle>0)

start=(Lvls_L[cycle])

stop=(Lvls_R[cycle])

endIf

CurveFit/Q/K={t_s[Lvls_L[cycle]]} poly_XOffset 3, Pos[Lvls_L[cycle],Lvls_R[cycle]] /X=t_s /D

Xcross = ((-K1+2*K2*t_s[Lvls_L[cycle]])/(2*K2))

FindValue/S=(Lvls_L[cycle]) /V=(Xcross) /T=0.015 t_s

Limits[cycle][%Cycle] = {cycle+1}

Limits[cycle][%StopRow] = {V_value+25}

Limits[cycle][%StartRow] = {V_value-25}

Limits[cycle][%Maxima] = {V_value}

endFor

Edit Limits.ld

InsertPoints 0,1,Limits

end

May 24, 2018 at 12:04 pm - Permalink