Interpolating a peak for two values of x - Python

Question

I want to find two values of x that intersect a certain value of y on a x-y plot of a resonance curve. However, as I have few data points I will need to interpolate to find these values of x.

The curve I am looking at can be seen below: enter image description here

How do I find the two values of x that equal a y value (shown in red)?

I have tried np.interpolate by splitting the data into two arrays: first with gradient(y)>0 and another with gradient(y)<0, however this yielded incorrect values. However, this approach is far from elegant and I seek a simple solution. Cheers in advance for any help.

Additional information: Code used thus far:

from numpy import *
import pylab as pl
import numpy as np
import scipy as scipy
from scipy import optimize

#Get data
fn = '4K_peak_hiresGhz.txt'
F_values, S_values, P_values = loadtxt(fn, unpack=True, usecols=[1, 2, 3])

#Select Frequency range of peak
lower = 4.995
upper = 5.002
F_values_2 = F_values[(F_values>lower) & (F_values<upper)]
S_values_2 = S_values[(F_values>lower) & (F_values<upper)]
P_values_2 = P_values[(F_values>lower) & (F_values<upper)]

#Calculate Q-value of selected peak
S_Peak = max(S_values_2)
print('S21 peak (dB):')
print(S_Peak)
print

F_Peak = F_values_2[S_values_2.argmax()]
print('Freq at S21 peak (GHz)')
print(F_Peak)
print

width_S = S_Peak - 3
print('S21 peak - 3dB (dB)')
print(width_S)
print

f, axarr = pl.subplots(2, sharex=False)
axarr[0].set_xlabel('Frequency (GHz)')
axarr[0].plot(F_values_2,S_values_2)
axarr[0].plot(F_values_2,S_values_2,'.')
axarr[0].set_ylabel('S21 (dB)')
axarr[0].axhline(y=width_S, linewidth='0.7', ls='dashed',color='red')
axarr[0].axhline(y=S_Peak, linewidth='1', ls='dashed', color='black')
axarr[0].axvline(x=F_Peak, linewidth='1', ls='dashed', color='black')
axarr[1].scatter(F_values_2, P_values_2)
axarr[1].set_ylabel('Phase (deg)')
axarr[1].set_xlabel('Frequency (GHz)')
pl.show()

Also, the data analysed in this program is found in this dpaste: http://dpaste.com/13AMJ92/

it would be great if you could post the `x, y` data that generated this plot... — Saullo G. P. Castro, May 26 '14 at 18:18
Good point. Data added in a dpaste and my initial bits of code added to my question. — user2992169, May 26 '14 at 18:29
Take a look at the function `find_transition_times` in my answer here: http://stackoverflow.com/questions/15112964/digitizing-an-analog-signal/15114952#15114952 You could use a variation of that function. It finds where the curve *increases* through a threshold, but getting any intersection should be an easy modification. — Warren Weckesser, May 26 '14 at 21:39
You mention in a comment down the page that it should be a lorenzian fit, is this correct? — will, Nov 08 '14 at 11:18
So to one of the answers which seemed completely legitimate to me you commented that you need it to be a Lorentzian fit. Perhaps you should add that to the question. Googling Lorentzian fit in python also brings you here: http://mesa.ac.nz/?page_id=1800 — Gullydwarf, Nov 13 '14 at 10:04

score 7 · Answer 1 · answered May 27 '14 at 02:41

Prepare curve data:

from numpy import *
import pylab as pl
import numpy as np
import scipy as scipy
from scipy import optimize

#Get data
fn = '13AMJ92.txt'
F_values, S_values, P_values = loadtxt(fn, unpack=True, usecols=[1, 2, 3])

#Select Frequency range of peak
lower = 4.995
upper = 5.002
F_values_2 = F_values[(F_values>lower) & (F_values<upper)]
S_values_2 = S_values[(F_values>lower) & (F_values<upper)]

S_Peak = max(S_values_2)
F_Peak = F_values_2[S_values_2.argmax()]
width_S = S_Peak - 3

split at peak, and use interp():

idx = S_values_2.argmax()
x1 = np.interp(width_S, S_values_2[:idx+1], F_values_2[:idx+1])
x2 = np.interp(width_S, S_values_2[-1:idx-1:-1], F_values_2[-1:idx-1:-1])
print x1, x2

the output:

4.99902583799 4.99954573393

You can also use Shapely:

from shapely import geometry
curve = geometry.asLineString(np.c_[F_values_2, S_values_2])
hline = geometry.LineString([(F_values_2[0], width_S), (F_values_2[-1], width_S)])
print np.asarray(curve.intersection(hline))

the output:

[[  4.99902584 -21.958     ]
 [  4.99954573 -21.958     ]]

If it's ok to use cubic splines, then you can use InterpolatedUnivariateSpline:

from scipy import interpolate
us = interpolate.InterpolatedUnivariateSpline(F_values_2, S_values_2 - width_S)
x1, x2 = us.roots()
print x1, x2
pl.plot(F_values_2, S_values_2 - width_S)
x = np.linspace(F_values_2[0], F_values_2[-1], 100)
pl.plot(x, us(x))
pl.axhline(0)
pl.plot([x1, x2], [0, 0], "o")

the output:

4.99891956723 4.99960633369

the plot:

enter image description here

This is nice but it shifts the resonance peak upwards and only finds routes at S_values_2 = 0, rather than S_Peak - 3... — user2992169, May 27 '14 at 08:45

Saullo G. P. Castro · Accepted Answer · 2014-05-26T19:30:23.547

4

You can use np.interp() in this case as:

np.interp(width_S, S_values_2, F_values_2)

In your plot the ticklabels are different from the thickvalues, you can check this fact comparing:

print(axarr[0].get_xticks())
print([s.get_text() for s in axarr[0].get_xticklabels()])

The interpolated value using this formula is correct, you have to set the right tick values to see that:

axarr[0].set_xticklabels(map(str, axarr[0].get_xticks()))

edited May 26 '14 at 19:30

answered May 26 '14 at 18:46

Saullo G. P. Castro

56,802
26
179
234

Plus, I seek _two_ x values that correspond to width_S. Can interp do this? – user2992169 May 26 '14 at 18:51
For example, when I do that I get the following plot - using `value = np.interp(width_S, S_values_2, F_values_2)` and `axarr[0].axvline(x=value,linewidth='1', ls='dashed', color='blue')`. Plot: http://i60.tinypic.com/14o4svl.png. Blue dashed line shows the incorrect interp! – user2992169 May 26 '14 at 18:53
@user2992169 check the update, the problem with your plot is that the `ticklabels` are different from the `tickvalues`... – Saullo G. P. Castro May 26 '14 at 19:32

Interpolating a peak for two values of x - Python

2 Answers2