FrontPage

2015/10/03からのアクセス回数 3905

ここで紹介したSageワークシートは、以下のURLからダウンロードできます。

http://www15191ue.sakura.ne.jp:8000/home/pub/66/

また、私の公開しているSageのサーバ（http://www15191ue.sakura.ne.jp:8000/）にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。

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

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

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

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

fig-7.png

↑

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

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

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

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

     [ 式 for 変数 in リスト ]

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

sageへの入力:

# 要素単位への演算をするため、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])

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

sageへの入力:

# 直前の結果は、_で参照できる
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])

sageへの入力:

round(_)

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で調節しています。

sageへの入力:

list_plot(_, figsize=5)

sage0.png

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

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

sageへの入力:

# 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]

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

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

sageへの入力:

# 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)]
↑

データフィッティング †

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

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

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

sageへの入力:

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

グラフの重ね合わせ †

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

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

sageへの入力:

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

sage1.png

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

sageへの入力:

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

sage2.png

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

sageへの入力:

# 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]
↑

音声データ †

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

fig-8.png

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

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

↑

準備 †

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

sageへの入力:

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

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

sageへの入力:

# 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のリスト内包を使って簡単にできます。

sageへの入力:

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を使って、 生成したサウンドデータの波形とスペクトログラムを表示し、音声をブラウザー確認できるようにします。

sageへの入力:

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

th_wave.jpg

sageへの入力:

# waveをファイルに保存
recordSound(wave, 'swep441k20-20k.wav', smpRate)

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

sageへの入力:

# specgramの形が異なるため、以前Interfaceで紹介されていたvoice_howling.wavで検証
# ファイルからの読み込み
rec, Fs = sf.read(DATA+"voice_howling.wav")
sample =rec.flatten()

sageへの入力:

# 以前と同じ結果がでたので、問題なしとする
plotSpecgram(sample, Fs)
saveAndShowPlot('sample.png')

th_sample.jpg

↑

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

図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を定義します。

sageへの入力:

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にセットします。

sageへの入力:

# eq２をRsについて整理
Rs_eq = solve(eq2, Rs); Rs_eq

[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}{B + T_2 log \left [ \frac{R (v_e - v_0) }{R_2 v_0} \right ]} \end{eqnarray} $$

sageへの入力:

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

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

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

sageへの入力:

# 他の変数に値をセットして、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)

sageへの入力:

plot(tsFunc(v0)-273, [v0, 0.1, 0.9], figsize=5)

sage3.png

↑

コメント †

選択肢 投票
おもしろかった 1  
そうでもない 0  
わかりずらい 0  

皆様のご意見、ご希望をお待ちしております。


(Input image string)

添付ファイル: fileth_wave.jpg 1127件 [詳細] fileth_sample.jpg 1029件 [詳細] fileswep441k20-20k.wav 1027件 [詳細] filesage3.png 1069件 [詳細] filesage2.png 1251件 [詳細] filesage1.png 999件 [詳細] filesage0.png 1109件 [詳細] filefig-11.png 538件 [詳細] filefig-9.png 540件 [詳細] filefig-8.png 1124件 [詳細] filefig-7.png 1233件 [詳細]
トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2015-10-04 (日) 14:12:22 (3397d)
SmartDoc