gnuplot-info

[Gnuplot-info] Finding thresholds with gnuplot

[theozh <th...@gm...>]

Hi,
I would like to find the x-value of some data when its y-value passes a certain threshold.
The curves can be very different, so I think that fitting to a certain function is probably not a solution.
I know gnuplot is probably not made for that and other tools should be used. However, sometimes it is very convenient if there is a pure gnuplot solution.
The basic idea is: 
Subtracting the threshold from the data, take the square of this and look for the minimum.
Maybe, I am trying way too complicated? If there are better ways please let me know.
Thanks, Theo.

The following sample generates some simple dummy data to illustrate.


### finding thresholds with gnuplot
reset
# generate some dummy data
set print $Data
do for [i=1:100] {
    print sprintf("%g %.3f",i*2+rand(0), rand(0)*10+i)
}
set print
print $Data

THRESHOLD = 60

# $Data2 holds (difference to threshold) squared
set table $Dummy
plot $Data u 1:(($2-THRESHOLD)**2) with table
unset table

# find minimum with stats
stats $Data2 u 1:2 nooutput
MINIMUM = STATS_index_min_y

# find x and y value of $Data
set table $Dummy
plot $Data u (MIN_X = $1,1):(MIN_Y=$2,$2) every ::MINIMUM::MINIMUM with table
unset table

set key top left
set multiplot layout 2,1

set label 1 sprintf("Value %g at %g",MIN_Y, MIN_X) at graph 0.5, graph 0.9
set arrow 1 from first MIN_X, graph 0 to first MIN_X, THRESHOLD

unset logscale y
plot $Data u 1:2, THRESHOLD w l 

set logscale y
plot $Data2 u 1:2 w l

unset multiplot
### end gnuplot code

Re: [Gnuplot-info] Finding thresholds with gnuplot

[theozh <th...@gm...>]

sorry, there were errors in the code.
This should work now...


### finding thresholds with gnuplot
reset

# generate some dummy data
set print $Data
do for [i=1:100] {
    print sprintf("%g %.3f",i*2+rand(0), rand(0)*10+i)
}
set print
print $Data

THRESHOLD = 60

set multiplot layout 2,1
set key top left

# $Dummy holds (difference to threshold) squared
set table $Dummy
plot $Data u 1:(($2-THRESHOLD)**2) with table
unset table

set logscale y
plot $Dummy u 1:2 w l

# find minimum with stats
stats $Dummy u 1:2 nooutput
MINIMUM = STATS_index_min_y

# find x and y value of $Data
set table $Dummy
plot $Data u (MIN_X=$1,$1):(MIN_Y=$2,$2) every ::MINIMUM::MINIMUM with table
unset table

set label 1 sprintf("Value %g at %g",MIN_Y, MIN_X) at graph 0.5, graph 0.9
set arrow 1 from first MIN_X, graph 0 to first MIN_X, THRESHOLD

unset logscale y
plot $Data u 1:2, THRESHOLD w l 

unset multiplot
### end gnuplot code

Re: [Gnuplot-info] Finding thresholds with gnuplot

[Ethan A M. <eam...@gm...>]

On Thursday, 18 October 2018 20:45:18 theozh wrote:
> Hi,
> I would like to find the x-value of some data when its y-value passes a certain threshold.

Not clear.
Do you want one of these?

1) a list of all x values for which y(x) > threshold

2) an interpolated y value for each pair of x, x+1 where
   y(x) and y(x+1) are on opposite sides of the threshold

3) only the leading edge where y(x) < threshold and y(x+1) > threshold


> The curves can be very different, so I think that fitting to a certain function is probably not a solution.
> I know gnuplot is probably not made for that and other tools should be used. However, sometimes it is very convenient if there is a pure gnuplot solution.
> The basic idea is: 
> Subtracting the threshold from the data, take the square of this and look for the minimum.

I doubt that is what you really want.  It answers yet another question:  
  
4) for what x values is y(x) closer to threshold than both y(x-1) and y(x+1)?
Solutions to that do not necessarily involve crossing the threshold value at all.


> Maybe, I am trying way too complicated? If there are better ways please let me know.

It depends on what do you want to do with a value when you find it.
Do you need all points that satisfy your criterion or would any single
point suffice?

	Ethan

-- 

Re: [Gnuplot-info] Finding thresholds with gnuplot

[theozh <th...@gm...>]

you're right, I wasn't clear. There are several possibilities. I'll try again.
Actually, I had in mind either a point or range where you would say it "passes" the threshold.
Assuming, as in the example, increasing (noisy) data values.

a) range between two values: 
the first x1 value where y(x1) is above the threshold and the last value x2>x1 where y(x2) is below the threshold

b) one value:
where your "eye would say" this x value is the "crossing" of the threshold line.
This probably would be a linear interpolation around the threshold value. > 4) for what x values is y(x) closer to threshold than both y(x-1) and y(x+1)?> Solutions to that do not necessarily involve crossing the threshold value at all.you're right, I assumed that a crossing exists.

One value (version b) or a two values (version a) would be enough.

Re: [Gnuplot-info] Finding thresholds with gnuplot

[Ethan M. <eam...@gm...>]

On Thu, Oct 18, 2018 at 12:51 PM theozh <th...@gm...> wrote:

>
> you're right, I wasn't clear. There are several possibilities. I'll try
> again.
> Actually, I had in mind either a point or range where you would say it
> "passes" the threshold.
> Assuming, as in the example, increasing (noisy) data values.
>
> a) range between two values:
> the first x1 value where y(x1) is above the threshold and the last value
> x2>x1 where y(x2) is below the threshold
>
> b) one value:
> where your "eye would say" this x value is the "crossing" of the threshold
> line.
> This probably would be a linear interpolation around the threshold value.
> > 4) for what x values is y(x) closer to threshold than both y(x-1) and
> y(x+1)?> Solutions to that do not necessarily involve crossing the
> threshold value at all.you're right, I assumed that a crossing exists.
>
> One value (version b) or a two values (version a) would be enough.
>


Here is a proof-of-principle example based on running_avg.dem.


set key fixed right center vertical Left reverse autotitle box
set key invert samplen 4 spacing 1
set xtics 10
set ytics
set title "Demonstrate use of assignment and serial evaluation
operators\nto accumulate statistics as su
ccessive data lines are read in\n"
set xrange [ 0.00000 : 57.0000 ] noreverse nowriteback
filter(this) = (prev < threshold) && (this >= threshold)

prev = 999.
threshold = 4.4
datafile = "silver.dat"

plot datafile using 0:3 with lp, \
    '' using 0:(filter($3) ? ($3) : NaN) : (prev=$3)  with impulse, \
    threshold with lines

#
# Now save them to a table
#
prev = 999.
set table $HITS
plot datafile using 0:(filter($3) ? ($3) : NaN) : (prev=$3)  with impulse
unset table

print $HITS

Re: [Gnuplot-info] Finding thresholds with gnuplot

[theozh <th...@gm...>]

Thank you, Ethan, for your hint.
ah, once again the powerful ternary operator ;-)
Theo.

The first example gets the first x value where y(x) is above the threshold and the last x value where y(x) is still below the threshold.
Actually, this method could also be used to find FWHM of spectral data (see second example). However, if the spectrum is a Gaussian or Lorentzian or whatever peak, fitting with appropriate parameters might probably give better results.

First example:

### find x range when data passes a certain threshold 
# getting the first x value where y(x) is above the threshold
# and the last x value wher y(x) is below the threshold
reset session
# generate some dummy data
set print $Data
do for [i=1:100] {
    print sprintf("%g %.3f",i*2+rand(0), rand(0)*20+i)
}
set print

Threshold = 50
     
set table $Dummy
plot $Data u 1:($2<Threshold ? ($2) : NaN):($2>=Threshold ? ($2) : NaN) with table
unset table
# print $Dummy

stats $Dummy u 1:3 nooutput
Min_x = STATS_min_x
stats $Dummy u 1:2 nooutput
Max_x = STATS_max_x

set label sprintf("min: %.1f, max: %.1f", Min_x, Max_x) at Min_x, Threshold+20 center

set key top left
plot $Data using 1:2 with lp lt 6, \
    Threshold w l,\
    $Dummy u 1:2 w impulse lc rgb "green",\
    $Dummy u 1:3 w impulse lc rgb "red"
### end of code


Second example:
### find Full Width at Half Maximum (FWHM)
reset session 
# generate some dummy data
set print $Data
Intensity = rand(0)*100
Position = rand(0)*200+480
Width = rand(0)*50+50
do for [i=380:780] {
    print sprintf("%g %.3f",i, Intensity*exp(-((i-Position)/Width)**2)+rand(0)*Intensity*0.1)
}
set print

stats $Data nooutput
Int_max = STATS_max_y
set table $Dummy
plot $Data u (Peak=$1,$1) every ::STATS_index_max_y::STATS_index_max_y with table
unset table
Threshold = 0.5*Int_max

set table $FWHM
plot $Data u ($2>=Threshold ? ($1) : NaN):($2>=Threshold ? Threshold : NaN) with table
unset table
stats $FWHM nooutput
FWHM = STATS_max_x - STATS_min_x

reset
set label 1 sprintf("Peak %.1f", Peak) at Peak*1.1, Int_max*0.95 center boxed
set label 2 sprintf("FWHM %.1f", FWHM) at Peak, Threshold*1.1 center boxed
set arrow 1 nohead from Peak,graph 0 to Peak,Int_max
set yrange[*:*]

plot $Data using 1:2 with l lc rgb "red" t "Spectrum", \
    $FWHM using 1:2 w l lc rgb "blue" not
    
### end of code

Re: [Gnuplot-info] Finding thresholds with gnuplot

[Ethan A M. <eam...@gm...>]

On Saturday, 20 October 2018 10:18:54 theozh wrote:
> Thank you, Ethan, for your hint.
> ah, once again the powerful ternary operator ;-)
> Theo.
> 
> The first example gets the first x value where y(x) is above the threshold and the last x value where y(x) is still below the threshold.
> Actually, this method could also be used to find FWHM of spectral data (see second example). However, if the spectrum is a Gaussian or Lorentzian or whatever peak, fitting with appropriate parameters might probably give better results.


I have been thinking about the possibility of somehow providing a window
into the sequence of data points, always centered about the "current" point
when plotting.  I don't know how this would work exactly, but I am imagining
something like a reserved array name "window" that could only be accessed
only from inside a smoothing function.
For example to produce a running average over 5 points:

   Avg(N) = 1/N * sum [i=1:N] window[i]
   plot 'data' using 1:2 smooth window Avg(5) 

Your case is a bit more complicated, but having the local data
points in an array rather than in a mess of separate variables
maintained by the ternary operator seems like it would be a win.

Thoughts or suggestions welcome.

	Ethan



> First example:
> 
> ### find x range when data passes a certain threshold 
> # getting the first x value where y(x) is above the threshold
> # and the last x value wher y(x) is below the threshold
> reset session
> # generate some dummy data
> set print $Data
> do for [i=1:100] {
>     print sprintf("%g %.3f",i*2+rand(0), rand(0)*20+i)
> }
> set print
> 
> Threshold = 50
>      
> set table $Dummy
> plot $Data u 1:($2<Threshold ? ($2) : NaN):($2>=Threshold ? ($2) : NaN) with table
> unset table
> # print $Dummy
> 
> stats $Dummy u 1:3 nooutput
> Min_x = STATS_min_x
> stats $Dummy u 1:2 nooutput
> Max_x = STATS_max_x
> 
> set label sprintf("min: %.1f, max: %.1f", Min_x, Max_x) at Min_x, Threshold+20 center
> 
> set key top left
> plot $Data using 1:2 with lp lt 6, \
>     Threshold w l,\
>     $Dummy u 1:2 w impulse lc rgb "green",\
>     $Dummy u 1:3 w impulse lc rgb "red"
> ### end of code
> 
> 
> Second example:
> ### find Full Width at Half Maximum (FWHM)
> reset session 
> # generate some dummy data
> set print $Data
> Intensity = rand(0)*100
> Position = rand(0)*200+480
> Width = rand(0)*50+50
> do for [i=380:780] {
>     print sprintf("%g %.3f",i, Intensity*exp(-((i-Position)/Width)**2)+rand(0)*Intensity*0.1)
> }
> set print
> 
> stats $Data nooutput
> Int_max = STATS_max_y
> set table $Dummy
> plot $Data u (Peak=$1,$1) every ::STATS_index_max_y::STATS_index_max_y with table
> unset table
> Threshold = 0.5*Int_max
> 
> set table $FWHM
> plot $Data u ($2>=Threshold ? ($1) : NaN):($2>=Threshold ? Threshold : NaN) with table
> unset table
> stats $FWHM nooutput
> FWHM = STATS_max_x - STATS_min_x
> 
> reset
> set label 1 sprintf("Peak %.1f", Peak) at Peak*1.1, Int_max*0.95 center boxed
> set label 2 sprintf("FWHM %.1f", FWHM) at Peak, Threshold*1.1 center boxed
> set arrow 1 nohead from Peak,graph 0 to Peak,Int_max
> set yrange[*:*]
> 
> plot $Data using 1:2 with l lc rgb "red" t "Spectrum", \
>     $FWHM using 1:2 w l lc rgb "blue" not
>     
> ### end of code