[[FrontPage]]
#contents
2012/03/31からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
http://www.pwv.co.jp:8000/home/pub/22/
私のSageサーバ（http://www15191ue.sakura.ne.jp:8000/）にユーザIDを作成することで、ダウンロードしたワークシートをアップロードし、実行したり、変更していろいろ動きを試すことができます。
* Sageで信号処理 [#p4db677b]
インターフェース2011年12号にMATLABの特集があったので、Sageで1章の例題を動かしてみます。
sageの特徴は、以下の3点だと思います。
- ブラウザーでどこからでも使える
- 実行可能なドキュメントである
- 既存ツールとの連携が可能
今回の例で、octaveとの連携および実行可能なドキュメントであることを余すとこなくご紹介します。
sageでoctaveを使う場合には、Sageのサーバマシンにoctaveと必要なパッケージがインストールされている
必要があります。
[[octaveのインストール]]
*** octave連携の準備 [#k90c932a]
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))
}}
** エコー混じりの音声データの処理 [#nb2fc3ed]
インターフェース2011年12号の1章の信号処理では、エコー混じりの音声データからエコーノイズを除去する例を紹介しています。
同様の処理をoctaveを使って試してみましょう。
*** 最初にデータの読み込み [#x233c881]
最初にエコー混じりの音声データ&ref(voice_howling.wav);を読み込みます。
sageへの入力:
#pre{{
# 音声データの読み込み
wav_file = DATA+"voice_howling.wav"
junk = octave.eval("[noised_sig,Fs] = wavread('%s');"%wav_file)
}}
*** 時間軸応答の表示 [#k2103cb1]
時間軸応答の表示は、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')
}}
&ref(fig1.png);
** 周波数特性の表示 [#me4b137d]
octave.evalでoctaveのコマンドを毎回括るのは大変ですし、ミスを多くなります。
そこで、octaveの処理をファイルにまとめて、sageのワークシートにアップロードして、
このファイルを利用すると便利です。
ファイルのアップロードは、
&ref(topmenu.png);
のDataプルダウンメニューから"Upload or create file..."を選択して
行います。
*** 周波数特性を表示するためのライブラリ [#c12c47e8]
今回の例題で周波数特性を表示するために、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)
}}
*** ノイズ信号の周波数応答の表示 [#ff7cbb51]
最新のoctave（3.6.1）ではmusic法は使えないので、ピリオドグラムとスペクトグラムを表示します。
ピリオドグラムでは2800と8500Hzあたりにピンと飛び出たところが見え、スペクトグラムでは同じ周波数に
濃い赤の帯が見えます。これがエコーノイズの正体です。
sageへの入力:
#pre{{
# ノイズ信号の周波数応答の表示
octave.eval("MySpectGram(noised_sig, Fs)")
saveAndShowPlot('spectrogram.png')
}}
&ref(fig2.png);
*** フィルターの作成 [#udd26ac7]
エコーを除去するために、以下の様なBRFフィルタを作成します。
| 応答タイプ | BRF |
| 設計手法 | IIR BRF 縦結合 |
| サンプリング周波数Fs(Hz) | 22050 |
| カット周波数 | 2842.4, 8527.1 |
| 帯域幅(Hz) | 100 |
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')
}}
&ref(fig3.png);
*** エコーノイズの除去 [#h37bd0a5]
出来上がったフィルターでエコーノイズを除去してみましょう。
信号処理に強いoctaveでは音声信号にフィルターを掛ける関数filterが用意されています。
エコーノイズを除去したdenoise_sigの周波数応答を表示すると、2842.4Hzと8527.1Hz
の出っ張りがとれ、赤い帯も薄くなっています。
sageへの入力:
#pre{{
# エコーノイズの除去
octave.eval("denoise_sig = filter(Zb, Za, noised_sig);")
# スペクトグラムの表示
octave.eval("MySpectGram(denoise_sig, Fs)")
saveAndShowPlot('deNoisedSpectrogram.png')
}}
&ref(fig4.png);
*** 結果の保存 [#y9191755]
ノイズの取れた音声を、&ref(denoised.wav); に保存します。
再生すると、ピーという音が小さくなっているのが分かります。
sageへの入力:
#pre{{
# 結果のファイル出力
denoise_wav_file = DATA+"denoised.wav"
junk = octave.eval("wavwrite(denoise_sig, Fs, '%s')"%denoise_wav_file)
}}
** フィルターの設計 [#w58374f1]
sageのインターラクティブ機能を使うと、フィルターの設計が楽にできます。
以下の様なフィルター作成関数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;
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への入力:
#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')
}}
&ref(fig5.png);
** デジタルフィルタの特性 [#n8f31f5c]
[[octave/IIRフィルタの伝達関数設計]]
では、octaveの関数を使ってアナログ・フィルタからデジタル・フィルタへの変換を行いました。
ここでは、インターフェース2006/09号の８章の図１３をsageを使って求めてみます。
&ref(fig6.png);
*** アナログフィルタ [#q17ee4d3]
最初に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])
}}
&ref(fig7.png);
*** 双一次z変換 [#vfe864bc]
つぎにアナログフィルタをデジタルフィルタに変換するために、以下の双一次z変換をします。
$$
s = \frac{2}{T}\frac{1 - z^{-1}}{1 + z^{-1}}
$$
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)))
}}
*** デジタルフィルタの周波数特性 [#c7254717]
デジタルフィルタに変換した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は簡単 [#m8c380b8]
octaveでアナログ・デジタルの変換を行ったときには、関数を調べたり、表示を工夫したり
しましたが、sageでは数式にそった形でアナログ・デジタルの変換をグラフにすることができました。
これが、数式処理システムsageの大きな特徴です。
** コメント [#yf5a77ed]
#vote(おもしろかった[4],そうでもない[0],わかりずらい[0])
#vote(おもしろかった[5],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
- お世話になります。いま、Octave で信号処理に挑戦しようとしているのですが、共同研究者がプログラミングに慣れておらず、できるだけの工程を Sage で数式化できないものか、と考えております。いま、取り組むべきは相関関数の算出です。そこで音声ファイルをこのページのように読み込んでみたものの、Sage の vector にコンヴァートし内積を計算するところで、無限ループ化してしまいます。Octave では所定の操作を完遂できるのですが、それが数式で表現した場合どうなるのか提示したいのです。お手数ですが wavread() で読み込んだ『ステレオ音声』を L チャンネル vector、R チャンネル vector に切り分けて内積を取るには如何したら宜しいでしょうか？ 以上、宜しくお願いします。 -- [[蓮風]] &new{2014-09-22 (月) 11:41:51};
- 蓮風さま、SageでWavデータ読み込み[[内積を求める例>http://www15191ue.sakura.ne.jp:8000/home/pub/46/]]に置きました。参考にして下さい。 -- [[竹本 浩]] &new{2014-09-22 (月) 17:13:17};
- ありがとうございます。いま、scikits などパッケージのインストールから試みています。がんばります。 -- [[蓮風]] &new{2014-09-22 (月) 19:18:17};
- 蓮風さま、sageへのaudiolabのインストール方法は、以下を参考にしてください。 -- [[竹本 浩]] &new{2014-09-22 (月) 19:40:28};
#pre{{
Sageをシェルモードで起動してpythonのパッケージをインストールします。
$ /usr/local/sage-6.0/sage -sh
(sage-sh) $
MacOSXの場合
新たにソースから作成した。
ソースのダウンロード
(sage-sh) $ tar xzvf ~/Downloaded/scikits.audiolab-0.11.0.tar.gz
(sage-sh) $ cd scikits.audiolab-0.11.0/
(sage-sh) $ cat >site.cfg <<EOF
# for MacOSX
[sndfile]
include_dirs = /opt/local/include
library_dirs = /opt/local/lib
sndfile_libs = sndfile
EOF
(sage-sh) $ python setup.py build
(sage-sh) $ python setup.py install
CentOSの場合
(sage-sh) $ wget https://pypi.python.org/packages/source/s/scikits.audiolab/scikits.audiolab-0.11.0.tar.gz#md5=f93f17211c7763d8631e0d10f37471b0
(sage-sh) $ cd ~/local/
(sage-sh) $ tar xvf ~/src/scikits.audiolab-0.11.0.tar.gz
(sage-sh) $ cd scikits.audiolab-0.11.0/
(sage-sh) $ cat >site.cfg <<EOF
[sndfile]
[sndfile]
include_dirs=/usr/include
library_dirs=/usr/lib64/
libraries=sndfile
EOF
(sage-sh) $ python setup.py build
(sage-sh) $ python setup.py install
}}
- パッケージもインストールできて、とりあえず相関係数が求められました。共同研究者に報告できます。ありがとうございます！！ -- [[蓮風]] &new{2014-09-23 (火) 00:49:55};
#comment_kcaptcha
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
