在Python中为两个x值插值峰值

Question

我想在一个共振曲线的 x-y 图上找到两个 x 值，这两个值会与某个特定的 y 值相交。不过，由于我手头的数据点不多，所以需要通过插值来找到这些 x 值。

我正在查看的曲线如下所示：

我该如何找到与红色标记的 y 值相等的两个 x 值呢？

我尝试过使用 np.interpolate，把数据分成两个数组：一个是 gradient(y)>0 的数据，另一个是 gradient(y)<0 的数据，但这样得到的值不正确。不过，这种方法看起来并不优雅，我希望能找到一个简单的解决方案。提前感谢任何帮助。

附加信息：

到目前为止使用的代码：

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

另外，这个程序分析的数据可以在这个链接找到： http://dpaste.com/13AMJ92/

Answer 1

准备曲线数据：

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

在峰值处进行分割，然后使用 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

输出结果：

4.99902583799 4.99954573393

你也可以使用 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))

输出结果：

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

如果可以使用三次样条插值，那么你可以使用 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")

输出结果：

4.99891956723 4.99960633369

绘图结果：

Answer 2

在这种情况下，你可以使用 np.interp()，用法如下：

np.interp(width_S, S_values_2, F_values_2)

在你的图表中，ticklabels（刻度标签）和 thickvalues（刻度值）是不一样的，你可以通过比较来确认这一点：

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

使用这个公式计算出的插值是正确的，你需要设置正确的刻度值才能看到这一点：

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