片対数グラフのスプライン補間

14718 ワード
numpy scipy グラフ matplotlib Python3 scipy テキストリンク

励みになるので、気軽にコメント、いいね、ぜひ。

片対数グラフでスプライン補間がうまくいかない

例として、ある実験のデータを片対数でプロットしてみる。

data
import numpy as np
np.array([[ 1.00e+03,  2.01e+01],
          [ 1.00e+04,  2.00e+01],
          [ 2.00e+04,  1.97e+01],
          [ 3.00e+04,  1.90e+01],
          [ 4.30e+04,  1.76e+01],
          [ 6.00e+04,  1.55e+01],
          [ 1.00e+05,  1.17e+01],
          [ 3.00e+05,  1.80e+00],
          [ 1.00e+06, -1.12e+01]])

#グラフの設定
plt.xscale('log')
plt.xlabel("Frequency[Hz]")
plt.ylabel("Gain[dB]")
plt.grid(True)

#散布図をプロット
plt.scatter(x, y, marker="x",s=50)

#スプライン補間した関数をプロット
from scipy import interpolate
f = interpolate.interp1d(x, y, kind='quadratic')
xnew = np.logspace(3, 6, num=90)
plt.plot(xnew,f(xnew))

▼出力されるグラフ

10^5～10^6あたりの補完がおかしい。
これは横軸をlogでとっているにもかかわらず、スプライン補間を元の軸のまま行っていることが原因である。

正しいスプライン補間をする方法

logをとってからスプライン補間した関数を求めれば良い。
以下に具体的なコードを示す。

#　～　省略　～
#スプライン補間した関数をプロット
from scipy import interpolate
f = interpolate.interp1d(np.log10(x), y, kind='quadratic') #logをとってスプライン補間
xnew = np.logspace(3, 6, num=90)
plt.plot(xnew,f(np.log10(xnew))) #plotするときもlogをとってから渡す

▼出力されるグラフ

今度はうまく補間されている！めでたい。

実際に試したい人向けコード

jupyter notebookを使っている人向け。


# coding: utf-8

# In[1]:


import numpy as np
import matplotlib.pyplot as plt

#プロットするデータ。
data = np.array([[ 1.00e+03,  2.01e+01],
                 [ 1.00e+04,  2.00e+01],
                 [ 2.00e+04,  1.97e+01],
                 [ 3.00e+04,  1.90e+01],
                 [ 4.30e+04,  1.76e+01],
                 [ 6.00e+04,  1.55e+01],
                 [ 1.00e+05,  1.17e+01],
                 [ 3.00e+05,  1.80e+00],
                 [ 1.00e+06, -1.12e+01]])

#logのグラフを扱うような人はグラフはグレースケールだろう
plt.style.use('grayscale')
plt.rcParams["axes.facecolor"] = (1,1,1,0)
plt.rcParams["figure.facecolor"] = (1,1,1,0)


# In[2]:


x = data[:,0]
y = data[:,1]

plt.xscale('log')
plt.xlabel("Frequency[Hz]")
plt.ylabel("Gain[dB]")
plt.grid(True)

plt.scatter(x, y, marker="x",s=50)

from scipy import interpolate
f = interpolate.interp1d(x, y, kind='quadratic')
xnew = np.logspace(3, 6, num=90)
plt.plot(xnew,f(xnew))


# In[3]:


x = data[:,0]
y = data[:,1]

plt.xscale('log')
plt.xlabel("Frequency[Hz]")
plt.ylabel("Gain[dB]")
plt.grid(True)

plt.scatter(x, y, marker="x",s=50)

from scipy import interpolate
f = interpolate.interp1d(np.log10(x), y, kind='quadratic')
xnew = np.logspace(3, 6, num=90)
plt.plot(xnew,f(np.log10(xnew)))