[[FrontPage]]

#contents

2012/04/03からのアクセス回数 &counter;

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

http://www.pwv.co.jp:8000/home/pub/23/

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

* Sageで信号処理 [#x436ae1b]

インターフェース2011年12号にMATLABの特集があったので、Sageで1章の例題を動かしてみます。

sageの特徴は、以下の3点だと思います。
* Sageでグラフを再現してみよう：微分方程式（オイラー法）編 [#kd0f2935]

- ブラウザーでどこからでも使える
- 実行可能なドキュメントである
- 既存ツールとの連携が可能
この企画は、雑誌や教科書にでているグラフをSageで再現し、
グラフの意味を理解すると共にSageの使い方をマスターすることを目的としています。

今回の例で、octaveとの連携および実行可能なドキュメントであることを余すとこなくご紹介します。
今回は、島根大学電子制御システム工学科の吉田和信氏の
[[Octaveによる動的システムシミュレーション入門>http://www.ecs.shimane-u.ac.jp/~kyoshida/octave(2002).pdf]]
の例題2.1の図2.1を題材にします。		

sageでoctaveを使う場合には、Sageのサーバマシンにoctaveと必要なパッケージがインストールされている
必要があります。
[[octaveのインストール]]
&ref(fig21.png);

** 例題2.1について [#y5692924]

*** octave連携の準備 [#xb5ce6aa]
以下のような連立一次常微分方程式をオイラー法を使って解く問題です。
$$
\dot x(t) = A x(t) + B u(t), x(0) = x_0
$$
オイラー法は、上記の式を差分近似を使って解く手法です。
$$
x(t +  \Delta t) \simeq x(t) + (A x(t) + B u(t)) \Delta t
$$

sageからoctaveを使う場合のグラフ処理を助ける関数を定義します。
（今回の改造で、グラフのアップデートに対応しました）
A, B, uは以下の様な行列とベクトルで与えられています。刻み幅\( \Delta t \) = 0.05, 終了時刻\( t_f \) = 20とします。

&ref(eq1.png);


*** octave連携の準備 [#tec3c795]

sageからoctaveを使う場合のグラフ処理を助ける関数を定義します。

sageへの入力:
#pre{{
import time
# octaveのグラフ出力ターミナルをpngにセット
#junk = octave.eval('putenv("GNUTERM", "png")')
# octave3.6.1の場合コメントアウトしてください。
# sageとoctaveがインストールされているマシンで、octaveのグラフを確認しながら
# 実行する場合には、上記の設定をコメントアウトしてください。 

# octaveのグラフをファイルに保存し、表示する関数を定義（画像のキャッシュをしない版）
def saveAndShowPlot(filename, fac=0.75):
    output = DATA + filename
    octave.eval("print -dpng '%s'"%output)
    width = int(640*fac)
    html('<img src="%s?%s" width="%spx">'%(filename, time.time(), width))
}}


** エコー混じりの音声データの処理 [#e0f8690d]

インターフェース2011年12号の1章の信号処理では、エコー混じりの音声データからエコーノイズを除去する例を紹介しています。
同様の処理をoctaveを使って試してみましょう。
** octaveを使ったオイラー法の解析例 [#da86c0f2]

最初は、octaveを使った解析例（解ex2_1.m）をそのままsageで実行します。ex2_1.mはワークシート上部のDataタグで内容を確認することができます。

*** 最初にデータの読み込み [#j6e413d8]
ベクトル\(x\)は、\( z, \dot z \)を要素に持ちます。

最初にエコー混じりの音声データ&ref(voice_howling.wav);を読み込みます。
$$
x(t) = \left[\begin{array}{c}
z(t)\\
\dot z (t)
\end{array}\right]	
$$		

図の青の線が、\(z(t)\)、緑の線が\(\dot z(t)\)に相当し、\(z(t)\)のピークで\(\dot z(t)\)が0になっているのが分かります。

sageへの入力:
#pre{{
# 音声データの読み込み
wav_file = DATA+"voice_howling.wav"
junk = octave.eval("[noised_sig,Fs] = wavread('%s');"%wav_file)
}}



*** 時間軸応答の表示 [#y3550e60]

時間軸応答の表示は、sageのoctave.evalインタフェースをそのまま使う方法で実行してみます。


sageへの入力:
#pre{{
# 時間軸応答の表示
octave.eval("t = (0:length(noised_sig)-1)/Fs;")
octave.eval("figure;")
octave.eval("plot(t,noised_sig),grid on;")
octave.eval("xlim([0 t(end)]),ylim([-1 1]);")
octave.eval("title('voice+howling noize time domain response');")
octave.eval("xlabel('Time(sec)'),ylabel('Magnitude');")
saveAndShowPlot('time_domain_response.png')
file = DATA + "ex2_1.m"
junk = octave.eval('source("%s")'%file);
saveAndShowPlot('octave_euler.png');
}}

&ref(fig1.png);

** 周波数特性の表示 [#ia9156c2]

octave.evalでoctaveのコマンドを毎回括るのは大変ですし、ミスを多くなります。
** Sageでのオイラー法の解析例 [#id7aaad8]

そこで、octaveの処理をファイルにまとめて、sageのワークシートにアップロードして、
このファイルを利用すると便利です。
次にSageを使って同じEuler法で解いてみます。ここではSageの特徴である、インタラクティブ機能を使って、
刻み幅dtを0.05, 01, 0.5, 1から選択できるようにしました。刻み幅によってグラフがどのように変化するかをみてください。

ファイルのアップロードは、
オイラー法は、連立一次常微分方程式を数値微分を使って解く方法であり、刻み幅によって解が正しく求まらないことがあることに注意が必要です。

&ref(topmenu.png);

のDataプルダウンメニューから"Upload or create file..."を選択して
行います。


*** 周波数特性を表示するためのライブラリ [#j7b0165f]

今回の例題で周波数特性を表示するために、MyFreqs, MyFreqz, MySpectGramの関数を定義して
sageからのコマンド入力を軽減します。


MyLib.mの内容は以下の通りです。（Dataプルダウンメニューからも内容を確認できます）

#pre{{
function H = MyFreqs(b, a, Ft)
    Fn = Ft/2;
	w = linspace(0, 2*pi*Fn, 512);
	f = w/(2*pi);
	% アナログフィルターの周波数特性を取得
	H = freqs(b, a, w);
	% プロット
	figure; 
	subplot(2, 1, 1);
	plot(f, abs(H));
	grid on;
	title('Filter response')
	xlabel('Frequency(Hz)')
	ylabel('Magnitude')
	
	subplot(2, 1, 2);
	plot(f, angle(H)*180);
	grid on
	xlabel('Frequency(Hz)')
	ylabel('Phase(degree)')
endfunction

function [H, W] = MyFreqz(b, a, Ft)
	[H, W] = freqz(b, a, 512, "half", Ft);
	Fn = Ft/2;
	w = linspace(0, 2*pi*Fn, 512);
	f = w/(2*pi);
	figure; 
	subplot(2, 1, 1);
	plot(f, abs(H));
	grid on;
	title('Filter response')
	xlabel('Frequency(Hz)')
	ylabel('Magnitude')
	
	subplot(2, 1, 2);
	plot(f, angle(H)*180);
	grid on
	xlabel('Frequency(Hz)')
	ylabel('Phase(degree)')
endfunction

function MySpectGram(sig, Fs)
	% ノイズ信号の周波数応答の表示
	figure
	subplot(2,1,1)
	[s1, f1] = periodogram(sig, hamming(length(sig)), length(sig), Fs);
	plot(f1, 20*log10(s1),'m-.'),grid
	xlim([0 f1(end)])
	title('Periodogram'),xlabel('Frequency(Hz)'),ylabel('Power/Frequency（dB/Hz）')
	subplot(2,1,2)
	specgram(sig, 512, Fs, hamming(512), 256)
	title('Spectrogram'),xlabel('Frequency(Hz)'),ylabel('Time(sec)')
endfunction
}}

MyLib.mの読み込みはoctaveのsourceコマンドを使います。
ワークシート内のファイルは、DATA+ファイル名で取り出すことができます。

sageへの入力:
#pre{{
# Freqs, Freqz, Specgramの修正版MyLib.mを読み込む
script = DATA + "MyLib.m"
junk = octave.eval("source('%s')"%script)
# ex2_1.mをsageで記述し、インタラクティブにdtを変更できるようにした例
from pylab import arange
A = matrix([[0, 1], [-1, -1]])
B = vector([0, 1])
u = 1
tf = 20
@interact
def _(dt=selector([0.05, 0.1, 0.5, 1])):
    # 微分方程式をEuler法で解く
    x = vector([0, 0])
    xx = []
    for t in arange(0, tf, dt):
        xx.append(x.list())
        dx = A*x + B*u
        x = x + dx*dt
    # 結果をプロット
    t = arange(0, tf, dt)
    x1_lst = [v[0] for v in xx]
    x2_lst = [v[1] for v in xx]
    lst_plt1 = list_plot(zip(t, x1_lst), plotjoined=True, linestyle ='-')
    lst_plt2 = list_plot(zip(t, x2_lst), plotjoined=True, linestyle =':')
    (lst_plt1 + lst_plt2).show()
}}

*** ノイズ信号の周波数応答の表示 [#x72db722]
&ref(fig2-1.png);

最新のoctave（3.6.1）ではmusic法は使えないので、ピリオドグラムとスペクトグラムを表示します。
&ref(fig2-2.png);

ピリオドグラムでは2800と8500Hzあたりにピンと飛び出たところが見え、スペクトグラムでは同じ周波数に
濃い赤の帯が見えます。これがエコーノイズの正体です。
** Sageの関数を使った例 [#wc6f2842]

Sageにも連立一次常微分方程式を解く関数がいくつか提供されています。ここでは、
desolve_system_rk4を使った例を示します。

sageへの入力:
#pre{{
# ノイズ信号の周波数応答の表示
octave.eval("MySpectGram(noised_sig, Fs)")
saveAndShowPlot('spectrogram.png')
}}
例題を\( z, \dot z \)で表すと、

&ref(fig2.png);
&ref(eq2.png);

*** フィルターの作成 [#k9982cfb]
これを展開すると、以下の様になります。

エコーを除去するために、以下の様なBRFフィルタを作成します。
$$
\left\{\begin{array}{l}
\dot z = \dot z\\
\ddot z = - z - \dot z + u(t)
\end{array} \right.
$$

| 応答タイプ | BRF |
| 設計手法 | IIR BRF 縦結合 |
| サンプリング周波数Fs(Hz) | 22050 |
| カット周波数 | 2842.4, 8527.1 |
| 帯域幅(Hz) | 100 |
desolve_system_rk4では、上記の右辺をと変数を与え、初期条件icsに\([ t_0, z(t_0), \dot z(t_0)]\)を与えて解きます。

octaveでは、複数のフィルターを連結するための関数sysmultが用意されていますので、2842.4Hzと8527.1Hz
の２つのバンド・リジェクト・フィルターを組み合わせてエコーノイズの除去のフィルターを作成します。

octaveの処理は、以下の様になります（createFilter.m）。
#pre{{
Ft = 22050;   % サンプリング周波数(Hz)
T = 1/Ft;   % サンプリング周期(sec)
Fn = Ft/2;  % ナイキスト周波数(Hz)
Fs1 = 2792.4; Ws1 = 2*pi*Fs1;
Fs2 = 2892.4; Ws2 = 2*pi*Fs2;
Fs3 = 8477.1; Ws3 = 2*pi*Fs3;
Fs4 = 8577.1; Ws4 = 2*pi*Fs4;
% プリワーピイング
Ws1p = 2/T*tan(Ws1*T/2);
Ws2p = 2/T*tan(Ws2*T/2);
Ws3p = 2/T*tan(Ws3*T/2);
Ws4p = 2/T*tan(Ws4*T/2);
% ５次のアナログ基準LPFを作成する
n = 5;
[z, p, g] = butter(n,1,'s');
% 基準LPFからバンドパスフィルタへの変換
[Z, P, G] = sftrans(z,p,g,[Ws1p Ws2p],true);
sys1 = zp2sys(Z, P, G, T);
[Z, P, G] = sftrans(z,p,g,[Ws3p Ws4p],true);
sys2 = zp2sys(Z, P, G, T);

[b, a] = sys2tf(sysmult(sys1, sys2));
[Zb, Za] = bilinear(b, a, 1/Ft);
MyFreqz(Zb, Za, Ft);    % プロット	
}}

完成したエコーノイズの除去のフィルターは以下の様な周波数特性となります。


sageへの入力:
#pre{{
# フィルタの作成
script = DATA + "createFilter.m"
octave.eval("source('%s')"%script)
saveAndShowPlot('BRF.png')
# Sageの連立一次常微分方程式を解く関数desolve_system_rk4を使った例
(z,z_dot,t) = var('z z_dot t')
P=desolve_system_rk4([z_dot,-z - z_dot + 1],[z,z_dot],ics=[0,0,0],ivar=t,end_points=20)
plt_f = list_plot([(i,j) for i, j, k in P], plotjoined=True, linestyle ='-')
plt_g = list_plot([(i,k) for i, j, k in P], plotjoined=True, linestyle =':')
(plt_f+plt_g).show()
}}

&ref(fig3.png);


*** エコーノイズの除去 [#c817a8e0]
** 伝達関数H(s)を使った解法 [#sec1544e]

出来上がったフィルターでエコーノイズを除去してみましょう。
信号処理に強いoctaveでは音声信号にフィルターを掛ける関数filterが用意されています。
\(\ddot z, \dot z, z, u\)をまとめると以下の２次の微分方程式になります。
$$
\ddot z + \dot z + z = u(t)
$$
上記の両辺をラプラス変換します。zのラブラス変換をZ, uのラブラス変換をUとすると以下の様になります。
$$
s^2 Z + s Z + Z = U
$$
ここで伝達関数$H(s) = \frac{Z}{U}$を計算します。
$$
H(s) = \frac{1}{s^2 + s + 1}
$$
単位ステップ応答のラプラス変換は、$\frac{1}{s}$なので、$Z(s)$は以下の様に求めまります。
$$
Z(s) = H(s) U(s) = H(s) \frac{1}{s} = \frac{1}{s^3 + s^2 + s}
$$
最後に$Z(s)$をラプラス逆変換すると$z(t)$の解析解が求まります。

エコーノイズを除去したdenoise_sigの周波数応答を表示すると、2842.4Hzと8527.1Hz
の出っ張りがとれ、赤い帯も薄くなっています。	
上記の処理をSageを使って行ってみます。


sageへの入力:
#pre{{
# エコーノイズの除去
octave.eval("denoise_sig = filter(Zb, Za, noised_sig);")
# スペクトグラムの表示
octave.eval("MySpectGram(denoise_sig, Fs)")
saveAndShowPlot('deNoisedSpectrogram.png')
(s, t) = var('s t')
P = 1/(s*s + s + 1)
# 単位ステップ応答は伝達関数/sの逆ラプラス変換で得られる
ip = inverse_laplace(P/s, s, t)
view(ip)
}}

&ref(fig4.png);


*** 結果の保存 [#xad02f09]
*** 結果のプロット [#r34495e0]

ノイズの取れた音声を、&ref(denoised.wav); に保存します。
再生すると、ピーという音が小さくなっているのが分かります。
求まった解をプロットすると図2.1と同じ結果となり、解析解と数値計算の結果が一致することが確認できました。

sageへの入力:
#pre{{
# 結果のファイル出力
denoise_wav_file = DATA+"denoised.wav"
junk = octave.eval("wavwrite(denoise_sig, Fs, '%s')"%denoise_wav_file)
plt_1 = plot(ip, [t, 0, 20], linestyle ='-')
plt_2 = plot(diff(ip,t), [t, 0, 20], linestyle =':')
(plt_1 + plt_2).show()
}}

&ref(fig5.png);

** フィルターの設計 [#pea8676c]

sageのインターラクティブ機能を使うと、フィルターの設計が楽にできます。
** 微分方程式の解 [#dc295304]

以下の様なフィルター作成関数MyFilterを作成します。
#pre{{
function [Zb, Za] = MyFilter(filter, type, fc1, fc2, Ft, n)
	T = 1/Ft;
	Fn = Ft/2;
	if type == "butter"
		[z, p, g] = butter(n,1,'s');
	elseif type == "cheby1"
		[z, p, g] = cheby1(n,3,1,'s');
	elseif type == "cheby2"
		[z, p, g] = cheby2(n,3,1,'s');
	else
		[z, p, g] = butter(n,1,'s');
	endif
	
	wc1 = 2*pi*fc1;
	wc2 = 2*pi*fc2;
最後に数式処理システムSageの微分方程式を解く関数desolveを使った例を示します。

	wcp1 = 2/T*tan(wc1*T/2);
	wcp2 = 2/T*tan(wc2*T/2);
desolveでは、変数tをvar関数で、変数tの関数ｚをfunction関数で定義し、
微分方程式deに\(\ddot z + \dot z + a = 1\)の関係式をセットし、
初期条件icsに\([ t_0, z(t_0), \dot z(t_0)]\)を与えて解きます。

	if filter == "LPF"
		[Z, P, G] = sftrans(z,p,g,wcp1,false);
	elseif filter == "HPF"
		[Z, P, G] = sftrans(z,p,g,wcp1,true);
	elseif filter == "BPF"
		[Z, P, G] = sftrans(z,p,g,[wcp1 wcp2],false);
	elseif filter == "BRF"
		[Z, P, G] = sftrans(z,p,g,[wcp1 wcp2],true);
	else
		[Z, P, G] = sftrans(z,p,g,wcp1,false);
	endif
	
	[b, a] = zp2tf(Z, P, G);
	[Zb, Za] = bilinear(b, a, 1/Ft);
endfunction
}}
解もグラフも先の伝達関数の例と同です。

@interactコマンドでブラウザーから入力した結果がその場で表示されます。type, filterの種類を変えて試してみて下さい。

sageへの入力:
#pre{{
# ユーザの入力でフィルターの形状を表示
# fc1:カットオフ周波数1, fc2:カットオフ周波数2, Ft:サンプリング周波数, n:フィルタの次数
@interact
def _(type=selector(["butter", "cheby1", "cheby2"]), filter = selector(["LPF", "HPF", "BPF", "BRF"]), fc1=30, fc2=40, Ft=100, n=5):
    #print filter, type, fc1, fc2, Ft, n
    octave.eval("[Zb, Za] = MyFilter('%s', '%s', %s, %s, %s, %s);"%(filter, type, fc1, fc2, Ft, n))
    octave.eval("MyFreqz(Zb, Za, %s);"%Ft)
    saveAndShowPlot('MyFilter.png')
t = var('t')
z = function('z', t)
de = diff(z, t, 2) + diff(z, t) + z == 1
sol = desolve(de, z, ics=[0,0,0]); view(sol)
}}

&ref(fig5.png);


** デジタルフィルタの特性 [#zf705861]

[[octave/IIRフィルタの伝達関数設計]]
では、octaveの関数を使ってアナログ・フィルタからデジタル・フィルタへの変換を行いました。

ここでは、インターフェース2006/09号の８章の図１３をsageを使って求めてみます。

&ref(fig6.png);

*** アナログフィルタ [#j4f6228c]

最初に100Hzのサンプリング周波数でカットオフ周波数を30Hzの１次のアナログLPF
$$
G(s) = \frac{1}{1+ \frac{s}{6o \pi}}
$$
の周波数特性をsageで表示してみます（図１３のaに相当）。

ラプラス演算子sと周波数の関係は
$$
s = i \omega = i (2 \pi f)
$$
なので、Gs.subsメソッドを使ってs=I*2*pi*fに置き換えて周波数f=0か50Hz
での特性をプロットします。

カットオフ周波数を30Hzで振幅が\(\frac{1}{\sqrt{2}}\)=(0.707)になっています。

sageへの入力:
#pre{{
# アナログフィルタの特性
(f, s) = var('f, s')
angle(x) = 180/pi*arctan(imag(x)/real(x))
Gs = 1/(1 + s/(60*pi))
# 振幅特性
plot(abs(Gs.subs(s=I*2*pi*f)), [f, 0, 50]).show()
# 周波数特性
g = Gs.subs(s=I*2*pi*f)
plot(angle(g), [f, 0, 50])
plt_sol1 = plot(sol, [t, 0, 20], linestyle ='-')
plt_sol2 = plot(diff(sol,t), [t, 0, 20], linestyle =':')
(plt_sol1 + plt_sol2).show()
}}

&ref(fig7.png);

*** 双一次z変換 [#u11e593e]

つぎにアナログフィルタをデジタルフィルタに変換するために、以下の双一次z変換をします。
$$
s = \frac{2}{T}\frac{1 - z^{-1}}{1 + z^{-1}}
$$
** まとめ [#p1f81fe3]

行列とベクトルで表される連立一次常微分方程式をオイラー法で数値微分することによって、
n次（今回は2次）の常微分方程式の解を求めることができる。

sageへの入力:
#pre{{
# Gsをs-z変換
(T, z, s, f) = var('T z s f')
Gz = Gs.subs(s = 2/T * (1-z^(-1))/(1+z^(-1)))
}}
さらにSageを使った解析解の求め方を紹介し、微分方程式をSageを使ってどのように解くかを図2.1
を例に説明しました。

是非、この機会にSageを使っていろいろな微分方程式を解いてみてください。

*** デジタルフィルタの周波数特性 [#p8a1cff1]

デジタルフィルタに変換したGzの周波数特性を求めてみます。		

zと\(\omega\)の関係式からGzにz=e^(I*2*pi*f*T)の置換を行い、周波数特性を表示します（図１３のb）。
$$ 
z = e^{i \omega T} 
$$


sageへの入力:
#pre{{
# z=exp(i w T)で変換し周波数特性を求める
# サンプリング周波数を100Hzとして計算
Fz(f) = Gz.subs(z=e^(I*2*pi*f*T)).subs(T=1/100); view(Fz)
}}

&ref(fig8.png);


残念ながら数式のままでは表示に失敗したので、代わりに数値に変換しlist_plotで表示します。


sageへの入力:
#pre{{
# 振幅特性
# 数式のままではプロットできないので、数値に変換しlist_plotで表示する
#plot(abs(Fz(f)), [f, 0, 50])
list_plot([ N(abs(Fz(f))) for f in range(1,50)],plotjoined=True)
}}

&ref(fig9.png);


位相特性は、そのまま計算できました。


sageへの入力:
#pre{{
# 位相特性
f = var('f')
plot(angle(Fz(f)), [f, 0, 50])
}}

&ref(fig10.png);

** sageは簡単 [#e004f63a]

octaveでアナログ・デジタルの変換を行ったときには、関数を調べたり、表示を工夫したり
しましたが、sageでは数式にそった形でアナログ・デジタルの変換をグラフにすることができました。

これが、数式処理システムsageの大きな特徴です。

** コメント [#na90a136]
#vote(おもしろかった,そうでもない,わかりずらい)

皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
トップ   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
SmartDoc