![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
2012/03/31からのアクセス回数 6980
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
http://www.pwv.co.jp:8000/home/pub/22/
私のSageサーバ（http://www.pwv.co.jp:8000/）にユーザIDを作成することで、ダウンロードしたワークシートをアップロードし、実行したり、変更していろいろ動きを試すことができます。
インターフェース2011年12号にMATLABの特集があったので、Sageで1章の例題を動かしてみます。
sageの特徴は、以下の3点だと思います。
今回の例で、octaveとの連携および実行可能なドキュメントであることを余すとこなくご紹介します。
sageでoctaveを使う場合には、Sageのサーバマシンにoctaveと必要なパッケージがインストールされている 必要があります。 octaveのインストール
sageからoctaveを使う場合のグラフ処理を助ける関数を定義します。 （今回の改造で、グラフのアップデートに対応しました）
sageへの入力:
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))
インターフェース2011年12号の1章の信号処理では、エコー混じりの音声データからエコーノイズを除去する例を紹介しています。 同様の処理をoctaveを使って試してみましょう。
最初にエコー混じりの音声データ
voice_howling.wavを読み込みます。
sageへの入力:
# 音声データの読み込み
wav_file = DATA+"voice_howling.wav"
junk = octave.eval("[noised_sig,Fs] = wavread('%s');"%wav_file)
時間軸応答の表示は、sageのoctave.evalインタフェースをそのまま使う方法で実行してみます。
sageへの入力:
# 時間軸応答の表示
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')
octave.evalでoctaveのコマンドを毎回括るのは大変ですし、ミスを多くなります。
そこで、octaveの処理をファイルにまとめて、sageのワークシートにアップロードして、 このファイルを利用すると便利です。
ファイルのアップロードは、
のDataプルダウンメニューから"Upload or create file..."を選択して 行います。
今回の例題で周波数特性を表示するために、MyFreqs, MyFreqz, MySpectGramの関数を定義して sageからのコマンド入力を軽減します。
MyLib.mの内容は以下の通りです。（Dataプルダウンメニューからも内容を確認できます）
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への入力:
# Freqs, Freqz, Specgramの修正版MyLib.mを読み込む
script = DATA + "MyLib.m"
junk = octave.eval("source('%s')"%script)
最新のoctave（3.6.1）ではmusic法は使えないので、ピリオドグラムとスペクトグラムを表示します。
ピリオドグラムでは2800と8500Hzあたりにピンと飛び出たところが見え、スペクトグラムでは同じ周波数に 濃い赤の帯が見えます。これがエコーノイズの正体です。
sageへの入力:
# ノイズ信号の周波数応答の表示
octave.eval("MySpectGram(noised_sig, Fs)")
saveAndShowPlot('spectrogram.png')
エコーを除去するために、以下の様なBRFフィルタを作成します。
| 応答タイプ | BRF |
| 設計手法 | IIR BRF 縦結合 |
| サンプリング周波数Fs(Hz) | 22050 |
| カット周波数 | 2842.4, 8527.1 |
| 帯域幅(Hz) | 100 |
octaveでは、複数のフィルターを連結するための関数sysmultが用意されていますので、2842.4Hzと8527.1Hz の２つのバンド・リジェクト・フィルターを組み合わせてエコーノイズの除去のフィルターを作成します。
octaveの処理は、以下の様になります（createFilter.m）。
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への入力:
# フィルタの作成
script = DATA + "createFilter.m"
octave.eval("source('%s')"%script)
saveAndShowPlot('BRF.png')
出来上がったフィルターでエコーノイズを除去してみましょう。 信号処理に強いoctaveでは音声信号にフィルターを掛ける関数filterが用意されています。
エコーノイズを除去したdenoise_sigの周波数応答を表示すると、2842.4Hzと8527.1Hz の出っ張りがとれ、赤い帯も薄くなっています。
sageへの入力:
# エコーノイズの除去
octave.eval("denoise_sig = filter(Zb, Za, noised_sig);")
# スペクトグラムの表示
octave.eval("MySpectGram(denoise_sig, Fs)")
saveAndShowPlot('deNoisedSpectrogram.png')
ノイズの取れた音声を、denoised.wav に保存します。
再生すると、ピーという音が小さくなっているのが分かります。
sageへの入力:
# 結果のファイル出力
denoise_wav_file = DATA+"denoised.wav"
junk = octave.eval("wavwrite(denoise_sig, Fs, '%s')"%denoise_wav_file)
sageのインターラクティブ機能を使うと、フィルターの設計が楽にできます。
以下の様なフィルター作成関数MyFilterを作成します。
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; wcp1 = 2/T*tan(wc1*T/2); wcp2 = 2/T*tan(wc2*T/2); 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への入力:
# ユーザの入力でフィルターの形状を表示
# 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')
octave/IIRフィルタの伝達関数設計 では、octaveの関数を使ってアナログ・フィルタからデジタル・フィルタへの変換を行いました。
ここでは、インターフェース2006/09号の８章の図１３をsageを使って求めてみます。
最初に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への入力:
# アナログフィルタの特性
(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])
つぎにアナログフィルタをデジタルフィルタに変換するために、以下の双一次z変換をします。 $$ s = \frac{2}{T}\frac{1 - z^{-1}}{1 + z^{-1}} $$
sageへの入力:
# Gsをs-z変換
(T, z, s, f) = var('T z s f')
Gz = Gs.subs(s = 2/T * (1-z^(-1))/(1+z^(-1)))
デジタルフィルタに変換したGzの周波数特性を求めてみます。
zと\(\omega\)の関係式からGzにz=e^(I*2*pi*f*T)の置換を行い、周波数特性を表示します（図１３のb）。 $$ z = e^{i \omega T} $$
sageへの入力:
# z=exp(i w T)で変換し周波数特性を求める # サンプリング周波数を100Hzとして計算 Fz(f) = Gz.subs(z=e^(I*2*pi*f*T)).subs(T=1/100); view(Fz)
残念ながら数式のままでは表示に失敗したので、代わりに数値に変換しlist_plotで表示します。
sageへの入力:
# 振幅特性 # 数式のままではプロットできないので、数値に変換しlist_plotで表示する #plot(abs(Fz(f)), [f, 0, 50]) list_plot([ N(abs(Fz(f))) for f in range(1,50)],plotjoined=True)
位相特性は、そのまま計算できました。
sageへの入力:
# 位相特性
f = var('f')
plot(angle(Fz(f)), [f, 0, 50])
octaveでアナログ・デジタルの変換を行ったときには、関数を調べたり、表示を工夫したり しましたが、sageでは数式にそった形でアナログ・デジタルの変換をグラフにすることができました。
これが、数式処理システムsageの大きな特徴です。
皆様のご意見、ご希望をお待ちしております。
