Finding a threshold in a wave

Share icon

Tue, 11/26/2013 - 07:24 pm

Hi,

I have a lot of graphs (curves) that have some kind of "threshold" that I want to determine. The threshold is sometimes sharp; others gradual. Also, sometimes there might be "noise" at the beginning and end of the curves (could be ignored). See:

http://s23.postimg.org/c9bwqjil7/image.jpg
http://s29.postimg.org/m4cx54i07/image.jpg
http://s13.postimg.org/k6zhjb2wn/image.jpg

I have wrote a program called as: threshold("Y_wave_name","X_wave_name" ,start_after,end_before,num_segments,percent)

Here is the important part of the program (everything except plotting section):
#pragma rtGlobals=1
#include <WindowBrowser>
Function threshold(ywave_name,xwave_name,start_after, end_before,no_of_segments,percent)
    Variable start_after            // First [start_after]% of wave will be excluded.  to exclude noise.
    Variable end_before             // Last [end_before]% of wave will be excluded.  to exclude noise.
    Variable no_of_segments         // No. of segments the wave be divided into.
    Variable percent                // e.g. 25% difference in average btw two consecutive segments is considered a threshold.
    String ywave_name           // y_axis wave name
    String xwave_name           // x_axis wave name
    Wave DCI =$xwave_name           // DCI is my x axis wave
    String wave_list = WaveList(ywave_name, ";", "")    // get a list of all waves
    String wave_list_s = SortList (wave_list,";",0)     // sort the list
    Variable nwaves = ItemsInList(wave_list)            // get the # of waves
    Make/O/N=(nwaves) ywave = 0             // New wave for the Y axis
// Program Variables                            //
    Variable wave_length = numpnts($StringFromList(0, wave_list_s))     // all waves have the same length
    Variable i                              // Waves counter
    Variable c                              // Segments counter
    Variable e                              // Averaging counter
    Variable seg_sum                        //
    Variable myavg_old                      //
    Variable myavg_new                      //
    Variable start_point = wave_length * start_after/100    // convert to points
    Variable end_point = wave_length * end_before/100   // convert to points   
    Variable seg_length = wave_length / no_of_segments    // Segment length (points)
    Variable offset                             // Useful to start at different "offsets" after 1st point in a given segment
for (i=0; i lessthan nwaves; i+=1)                      // Loop all waves.
        Wave mywave = $StringFromList(i, wave_list_s)   // Get wave name from the list of available waves.
        for (offset=0; offset lessthan seg_length; offset+=1)   // For each offset, do the work below.
            for (c=start_point+offset; c lessthan wave_length-end_point; c=c+seg_length)    // Loop all segments
                seg_sum = 0                                     // reset seg_sum
                for (e=0; e lessthan seg_length; e+=1)      // Get the segment average
                    seg_sum += mywave[c+e]      //
                endfor                      //
                myavg_new = seg_sum / seg_length    //
                if (myavg_new greaterthan (1+percent/100)*myavg_old && myavg_old !=0)   // If yes, then we found the threshold.
                     ywave[i] = abs(DCI[c])     // Assign the corresponding xwave (DCI) point to the threshold wave.
                    break                       // And "break" to get out of for loop.
                endif                           //
                myavg_old = myavg_new           //
            endfor                          // For the next segment, "new" avg becomes "old".
        endfor                              //
    endfor                                  //
end


The code is not giving usefull thresholds, or lets say correct thresholds. I tried tuning it (e.g. 20,20,10,25: excluding first and last 20% of the wave, dividing it to 10 segments and defining a threshold as 25%). It looks like it works with some waves but not the others. I do not think the problem is in the code as a program, but in the algorithm used (I think it is too "naive").

P.S looks like the forum cut the code after either < or > so I replaced them by lessthan and greaterthan.
I haven't looked at your code too closely (btw, you can use the 'igor' html tag to make your code more readable in the forum). But it seems to me that your problem is calling for the 'FindLevel' (or FindLevels) function with the 'box' parameter applied for noise averaging. You want to find a specific threshold for each of your waves? Then just loop through them and call FindLevel for each, and then save the result somewhere.

November 27, 2013 at 02:33 am - Permalink

FindLEvel assumes I know the value at which the "crossing" will happen, but I do not. I only know it is a relative increase in amplitude by, say, 20%.

However, I found another function, EdgeStats...but even this do not give correct values.

November 27, 2013 at 08:43 pm - Permalink

ali8 wrote:
... relative increase in amplitude by, say, 20%.



Relative to what? The background?
What would be "correct values" for the three examples above?

November 27, 2013 at 10:53 pm - Permalink

Relative to the 'rest' of the signal.

For example, in the first plot, the threshold (x value) is ~ +/- 20, for the second plot it is hard to find it and it may even not be well defined, but probably around +/-50, third is +/- 10.

November 28, 2013 at 09:45 am - Permalink

This may not be the most efficient method, but I would implement a routine as follows for the data you posted.

1. Subtract baseline (say the mean of the first 10-20 data points)
2. Find max value and point number (use wavestats)
3. set reference value (a variable) to a % of the max value
4. Break the wave into two segments (right and left of max value location)
5. Evaluate each of the segments for crossing the reference value.

December 2, 2013 at 09:52 am - Permalink

I think you are assuming that the max value always appear in the middle. This actually not generally true.

December 3, 2013 at 02:57 pm - Permalink

I guess without stating clearly what you want to achieve it's difficult to understand your problem for everyone (or for me at least). Do you want to find the position where the signal 'switches'? Maybe first calculating a smoothed derivative of your data and looking for peaks instead will a more promising approach.

December 3, 2013 at 08:13 pm - Permalink

ali8 wrote:
I think you are assuming that the max value always appear in the middle. This actually not generally true.


Since my approach looks for the peak position, the max value does not have to be in the middle. I am, however, assuming that there is a rising and falling edge since your example data appeared this way. (A complete peak exists in the data).

Looking back on your graph with low signal amplitude, it may be best to smooth the data (since the low signal to noise ratio may make it difficult to actually determine the threshold crossing).

December 4, 2013 at 06:34 am - Permalink

@ Proland: I have already tried smoothing it is not very useful.

@Chozo: I will try the derivative and see how it works.

Thanks all, I will update you later on how things work.

December 4, 2013 at 09:34 am - Permalink