Interface201511

3125 days ago by takepwave

Hiroshi TAKEMOTO (take@pwv.co.jp)

Sageでグラフを再現してみよう：Interface 2015/11

この企画は、雑誌や教科書にでているグラフをSageで再現し、 グラフの意味を理解すると共にSageの使い方をマスターすることを目的としています。

Interface 2015/11号にMathematicaを使った例題が紹介されていたので、 これをSageを使って試してみます。

最初に図7（Interface 2015/11から引用）に紹介されている例題をSageでやってみます。

10ビット三角関数テーブルを作る

図7の最初の部分は、sinのテーブルを作成することです。

Sageのリストでは、リストの要素単位に演算をできません。 このような場合には、numpyのarrayを使います。

pythonでリストを作るときには、以下の書式のリスト内包を使います。 

	[ 式 for 変数 in リスト ]
forで変数に値をセットするためのリストには、通常rangeやarangeが使われますが、 今回は最後の2πを含めたリストを作るために、0から2πを均等分割するlinspaceを使いました。

# 要素単位への演算をするため、numpyのarrayを使用 # arangeだと最後の0が含まれないのでlinspaceで21個出力するように変更 import numpy as np np.array([sin(N(x)) for x in np.linspace(0, 2*pi, 21)]) 
        
array([  0.00000000e+00,   3.09016994e-01,   5.87785252e-01,
         8.09016994e-01,   9.51056516e-01,   1.00000000e+00,
         9.51056516e-01,   8.09016994e-01,   5.87785252e-01,
         3.09016994e-01,   1.22464680e-16,  -3.09016994e-01,
        -5.87785252e-01,  -8.09016994e-01,  -9.51056516e-01,
        -1.00000000e+00,  -9.51056516e-01,  -8.09016994e-01,
        -5.87785252e-01,  -3.09016994e-01,  -2.44929360e-16])
array([  0.00000000e+00,   3.09016994e-01,   5.87785252e-01,
         8.09016994e-01,   9.51056516e-01,   1.00000000e+00,
         9.51056516e-01,   8.09016994e-01,   5.87785252e-01,
         3.09016994e-01,   1.22464680e-16,  -3.09016994e-01,
        -5.87785252e-01,  -8.09016994e-01,  -9.51056516e-01,
        -1.00000000e+00,  -9.51056516e-01,  -8.09016994e-01,
        -5.87785252e-01,  -3.09016994e-01,  -2.44929360e-16])

直前の結果を参照するには、_(アンダーバー)を使います。 以下の様に直前の結果を1000倍し、round関数で小数点以下を取り除き、 10ビット三角関数テーブルを作ります。

# 直前の結果は、_で参照できる 1023*_ 
        
array([  0.00000000e+00,   3.16124385e+02,   6.01304313e+02,
         8.27624385e+02,   9.72930816e+02,   1.02300000e+03,
         9.72930816e+02,   8.27624385e+02,   6.01304313e+02,
         3.16124385e+02,   1.25281368e-13,  -3.16124385e+02,
        -6.01304313e+02,  -8.27624385e+02,  -9.72930816e+02,
        -1.02300000e+03,  -9.72930816e+02,  -8.27624385e+02,
        -6.01304313e+02,  -3.16124385e+02,  -2.50562735e-13])
array([  0.00000000e+00,   3.16124385e+02,   6.01304313e+02,
         8.27624385e+02,   9.72930816e+02,   1.02300000e+03,
         9.72930816e+02,   8.27624385e+02,   6.01304313e+02,
         3.16124385e+02,   1.25281368e-13,  -3.16124385e+02,
        -6.01304313e+02,  -8.27624385e+02,  -9.72930816e+02,
        -1.02300000e+03,  -9.72930816e+02,  -8.27624385e+02,
        -6.01304313e+02,  -3.16124385e+02,  -2.50562735e-13])
round(_) 
        
array([    0.,   316.,   601.,   828.,   973.,  1023.,   973.,   828.,
         601.,   316.,     0.,  -316.,  -601.,  -828.,  -973., -1023.,
        -973.,  -828.,  -601.,  -316.,    -0.])
array([    0.,   316.,   601.,   828.,   973.,  1023.,   973.,   828.,
         601.,   316.,     0.,  -316.,  -601.,  -828.,  -973., -1023.,
        -973.,  -828.,  -601.,  -316.,    -0.])

結果をグラフにプロットする

リストをグラフにプロットするには、list_plot関数を使います。 そのままだとグラフが大きく表示されるので、figsize=5で調節しています。

list_plot(_, figsize=5) 
       

最初のリスト内包の式で、sin(N(x))としていましたが、sin(x)としたらどうなるのでしょう？

Sageでは、数式がそのまま使われるので、分数やsin等の関数がそのまま出力されます。

# Nを使わないと、数式のまま表示される [sin(x) for x in np.arange(0, 2*pi, pi/10)] 
        
[0,
 1/4*sqrt(5) - 1/4,
 sin(1/5*pi),
 1/4*sqrt(5) + 1/4,
 sin(2/5*pi),
 1,
 sin(3/5*pi),
 1/4*sqrt(5) + 1/4,
 sin(4/5*pi),
 1/4*sqrt(5) - 1/4,
 0,
 -1/4*sqrt(5) + 1/4,
 sin(6/5*pi),
 -1/4*sqrt(5) - 1/4,
 sin(7/5*pi),
 -1,
 sin(8/5*pi),
 -1/4*sqrt(5) - 1/4,
 sin(9/5*pi),
 -1/4*sqrt(5) + 1/4]
[0,
 1/4*sqrt(5) - 1/4,
 sin(1/5*pi),
 1/4*sqrt(5) + 1/4,
 sin(2/5*pi),
 1,
 sin(3/5*pi),
 1/4*sqrt(5) + 1/4,
 sin(4/5*pi),
 1/4*sqrt(5) - 1/4,
 0,
 -1/4*sqrt(5) + 1/4,
 sin(6/5*pi),
 -1/4*sqrt(5) - 1/4,
 sin(7/5*pi),
 -1,
 sin(8/5*pi),
 -1/4*sqrt(5) - 1/4,
 sin(9/5*pi),
 -1/4*sqrt(5) + 1/4]

リストの部分列を取り出すには、[始まりのインデクス:終わりのインデックス+1]の形式で指定します。

以下は、aから最初の0, 1, 2番目の要素を取り出すために、a[0:3]とした結果です。 丸括弧()で括られている部分は、タプルと呼ばれ、座標点などを表す時によく使われます。

# pythonは0始まりで、終わりのインデックスの前まで出力するので0:3の指定で、0, 1,2の要素が出力される a = [(N(x), sin(N(x))) for x in np.linspace(0, 2*pi, 21)]; a[0:3] 
        
[(0.000000000000000, 0.000000000000000),
 (0.314159265358979, 0.309016994374947),
 (0.628318530717959, 0.587785252292473)]
[(0.000000000000000, 0.000000000000000),
 (0.314159265358979, 0.309016994374947),
 (0.628318530717959, 0.587785252292473)]

データフィッティング

Sageでデータに最も合ったモデルを求める関数がfind_fit関数です。 find_fitには、データとモデルを入力としますが、モデルに使用する変数は予めvar関数で宣言する必要があります。

以下では5次関数をモデルとする、model(x)関数を変数（係数）を使って定義します。

solution_dict=Trueは、結果を辞書型で返すように指定します。このオプションを使うと、 find_fitの結果を使った関数がとても簡単に求まります。

(x, a0, a1, a2, a3, a4, a5) = var('x a0 a1 a2 a3 a4 a5') model(x) = a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4 + a5*x^5 fit = find_fit(a, model, solution_dict=True) #show(model) #show(fit) print model print fit 
        
x |--> a5*x^5 + a4*x^4 + a3*x^3 + a2*x^2 + a1*x + a0
{a4: 0.08628153266267449, a1: 0.8964456756065791, a5: -0.005492852966760598, a3: -0.38820422530842147, a2: 0.25248024823120263, a0: 0.004935764502076758}
x |--> a5*x^5 + a4*x^4 + a3*x^3 + a2*x^2 + a1*x + a0
{a4: 0.08628153266267449, a1: 0.8964456756065791, a5: -0.005492852966760598, a3: -0.38820422530842147, a2: 0.25248024823120263, a0: 0.004935764502076758}

グラフの重ね合わせ

データのリストプロットとフィッティングで求めた関数b(x)を合わせてプロットするには、 プロット関数の結果を+記号で足し合わせます。

関数のプロットには、plot関数を使い、[x, 0, 2*pi]でプロット範囲を指定し、 rgbcolor='green'で線色を緑に指定します。 図7ではSin関数とフィッティング曲線の重ね合わせですが、ここでは元のデータとフィッティング曲線の重ね合わを表示しています。

b(x) = model.subs(fit) list_plot(a, figsize=5)+plot(b(x), [x, 0, 2*pi], rgbcolor='green', figsize=5) 
       

次に、Sin曲線とフィッティング曲線の差をプロットします。plot関数の第1引数には式を指定するので、 b(x)-sin(x)をそのまま書くとSin曲線とフィッティング曲線の差がプロットできます。

plot(b(x)-sin(x), [x, 0, 2*pi], figsize=5) 
       

最初の三角関数テーブルを1行のリスト内包で書くと以下の様になります。

# pythonのリスト内包の表現は強力です [N(round(1023*sin(x))) for x in np.linspace(0, 2*pi, 21)] 
        
[0.000000000000000,
 316.000000000000,
 601.000000000000,
 828.000000000000,
 973.000000000000,
 1023.00000000000,
 973.000000000000,
 828.000000000000,
 601.000000000000,
 316.000000000000,
 0.000000000000000,
 -316.000000000000,
 -601.000000000000,
 -828.000000000000,
 -973.000000000000,
 -1023.00000000000,
 -973.000000000000,
 -828.000000000000,
 -601.000000000000,
 -316.000000000000,
 0.000000000000000]
[0.000000000000000,
 316.000000000000,
 601.000000000000,
 828.000000000000,
 973.000000000000,
 1023.00000000000,
 973.000000000000,
 828.000000000000,
 601.000000000000,
 316.000000000000,
 0.000000000000000,
 -316.000000000000,
 -601.000000000000,
 -828.000000000000,
 -973.000000000000,
 -1023.00000000000,
 -973.000000000000,
 -828.000000000000,
 -601.000000000000,
 -316.000000000000,
 0.000000000000000]

音声データ

図8（Interface 2015/11から引用）では、Mathematicaを使った音声データの処理が紹介されています。

Sageには音声データを扱う関数がないので、Pythonのライブラリを組み合わせて使います。

今回は、音声データの取り込みにSoundFileライブラリを使用し、音声データの書き込みと読み込み、 波形とスペクトログラムの表示を紹介するのに留めます。

準備

最初に、必要なライブラリをインポートします。

# 音声データの処理用パッケージ import soundfile as sf import pylab as pl import scipy.fftpack as Fourier from scipy.signal import butter, lfilter, freqz import time 
       

波形とスペクトログラムをプロットする関数と音声データの書き込み関数を定義します。

# pylabを使ったスペックグラムと波形の処理用に関数を定義 # 画像ファイルの保存と表示関数 def saveAndShowPlot(filename, fac=0.75): output = DATA + filename pl.savefig(output) width = int(640*fac) html('<img src="%s?%s" width="%spx">'%(filename, time.time(), width)) # specgramと波形のプロット関数 def plotSpecgram(wave, Fs): # FFTのサンプル数 N = 512 # FFTで用いるハミング窓 hammingWindow = np.hamming(N) # specgramのプロット t = range(0, len(wave)) pl.figure() pl.subplot(2, 1, 1) pl.specgram(wave, NFFT=N, Fs=Fs, noverlap=N/2, window=hammingWindow) # waveのプロット pl.subplot(2, 1, 2) pl.plot(t, wave, 'b') pl.grid() # ファイル保存 def recordSound(sources, name, Hz): filename = DATA + name sf.write(sources, filename, samplerate=Hz) # 結果を確認する方法 html('<a href="%s">録音結果（右クリックでファイルをダウンロードして再生してください）</a>' % name) 
       

サウンドデータの作成

サウンドデータとして1秒間で20Hz〜20kHzに上昇するスイープ音を生成します。

これもPythonのリスト内包を使って簡単にできます。

smpRate = 44100 # サンプリングレート44.1kHz fLo = 20 fHi = 20000 wave = [N(sin(2*pi*(fLo+(fHi/2-fLo)*t)*t)) for t in np.linspace(0, 1, smpRate+1)] 
       

音声データのプロットと保存用に定義したplotSpecgram関数とsaveAndShowPlot関数、recordSoundを使って、 生成したサウンドデータの波形とスペクトログラムを表示し、音声をブラウザー確認できるようにします。

plotSpecgram(wave, smpRate) saveAndShowPlot('wave.png') 
       

スペクトログラムの形が図８（ソノグラフ）と異なるので、以前InterfaceのMatlabの例題で使われていたエコーデータ のスペクトログラムをプロットしてみました。正しく表示されているので、Mathematicaのソノグラフの定義とスペクトログラム の定義が異なるのかも知れませんが、ここではこれ以上深入りはしないことにします。

# specgramの形が異なるため、以前Interfaceで紹介されていたvoice_howling.wavで検証 # ファイルからの読み込み rec, Fs = sf.read(DATA+"voice_howling.wav") sample =rec.flatten() 
       
# 以前と同じ結果がでたので、問題なしとする plotSpecgram(sample, Fs) saveAndShowPlot('sample.png') 
       

連立方程式を使った応用例

図11（Interface 2015/11から引用）では、サーミスタの抵抗と図12のサーミスタ温度測定回路の$R_s、R$の分圧の関係式 から$v_0$と$T_s$の関係式を求めています。 $$ \begin{eqnarray} eq_1 & = & R_s == R_2 e^{B \left ( \frac{1}{T_s} - \frac{1}{T_2} \right )} \\ eq_2 & = & v_0 == \frac{R}{R_s + R} v_e \end{eqnarray} $$

最初に使用する変数R, Rs, R2, B, Ts, T2, ve, v0を変数定義し、eq1, eq2を定義します。

R, Rs, R2, B, Ts, T2, ve, v0 = var('R Rs R2 B Ts T2 ve v0') # 以下の関係式が与えられた場合のTsとv0の関係を求める eq1 = Rs == R2*exp(B*(1/Ts - 1/T2)) eq2 = v0 == R/(Rs + R)*ve 
       

solve関数の使い方

solve関数を使って連立方程式を解きます。 本来なら、以下の様にすれば解けるのですが、 

	solve([eq1, eq2], [Ts,Rs])
eq2にTsが含まれていないため、Sageでは上手く解くことができませんでした。

ここでは、連立方程式の解法を順番にSageで行って解を求めていきます。

最初にeq2からRsを求め、その結果をRs_eqにセットします。

# eq２をRsについて整理 Rs_eq = solve(eq2, Rs); Rs_eq 
        
[Rs == -R*(v0 - ve)/v0]
[Rs == -R*(v0 - ve)/v0]

つぎにRs_eqをeq1に代入し、Tsを求め、右辺の式を取りだします。 solve関数の結果はリスト形式で返るので、[0]で最初の1個を取りだしています。

求まった解は、図11の表現とは異なりますが、logに-1を掛けることが逆数になることから同じ結果であることが確認できます。 $$ \begin{eqnarray} T_s & = & -\frac{B T_2}{T_2 log \left [ -\frac{R_2 v_0}{R v_0 - R v_e} \right ] - B} \\ & = & -\frac{B T_2}{T_2 log \left [ \frac{R_2 v_0}{R v_e - R v_0 } \right ] - B} \\ & = & \frac{-B T_2}{-T_2 log \left [ \frac{R v_e - R v_0 }{R_2 v_0} \right ] - B} \\ & = & \frac{B T_2}{T_2 log \left [ B + \frac{R v_e - R v_0 }{R_2 v_0} \right ]} \end{eqnarray} $$

# Rs_eqの結果で、eq1を解く sol = solve(eq1.subs(Rs_eq[0]), Ts)[0].rhs() #show(sol) sol 
        
-B*T2/(T2*log(-R2*v0/(R*v0 - R*ve)) - B)
-B*T2/(T2*log(-R2*v0/(R*v0 - R*ve)) - B)

次に変数R, ve, R2, T2, Bに値をセットし、その時のTsとv0の関係式を求めます。 これには、変数名:値の対をカンマで区切り、{}で囲んだ辞書型を使い、subsメソッドで変数の値と置換します。

これで目的とするTsを求める関数tsFunc(v0)が求まりました。 plot関数を使ってtsFuncをプロットすると図11のTsとv0の関係が表示されます。

# 他の変数に値をセットして、Tsを変数v0で定義した関数にする tsFunc(v0) = sol.subs({R:10000, ve:1, R2: 30000, T2: 273,B: 3420}) #show(tsFunc) tsFunc 
        
v0 |--> -311220/(91*log(-3*v0/(v0 - 1)) - 1140)
v0 |--> -311220/(91*log(-3*v0/(v0 - 1)) - 1140)
plot(tsFunc(v0)-273, [v0, 0.1, 0.9], figsize=5)