sage/Sageで信号処理
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2012/03/31からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www.pwv.co.jp:8000/home/pub/22/
私のSageサーバ（http://www15191ue.sakura.ne.jp:8000/）に...
* Sageで信号処理 [#p4db677b]
インターフェース2011年12号にMATLABの特集があったので、Sag...
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...
}}
** エコー混じりの音声データの処理 [#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...
}}
*** 時間軸応答の表示 [#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 respo...
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,...
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...
plot(f1, 20*log10(s1),'m-.'),grid
xlim([0 f1(end)])
title('Periodogram'),xlabel('Frequency(Hz)'),ylabel('Pow...
subplot(2,1,2)
specgram(sig, 512, Fs, hamming(512), 256)
title('Spectrogram'),xlabel('Frequency(Hz)'),ylabel('Tim...
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が...
の２つのバンド・リジェクト・フィルターを組み合わせてエコ...
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では音声信号にフィルターを掛ける関数f...
エコーノイズを除去したdenoise_sigの周波数応答を表示すると...
の出っ張りがとれ、赤い帯も薄くなっています。
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')"%deno...
}}
** フィルターの設計 [#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コマンドでブラウザーから入力した結果がその場で表...
sageへの入力:
#pre{{
# ユーザの入力でフィルターの形状を表示
# fc1:カットオフ周波数1, fc2:カットオフ周波数2, Ft:サンプ...
@interact
def _(type=selector(["butter", "cheby1", "cheby2"]), filt...
#print filter, type, fc1, fc2, Ft, n
octave.eval("[Zb, Za] = MyFilter('%s', '%s', %s, %s, ...
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の...
$$
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に置き換えて周波...
での特性をプロットします。
カットオフ周波数を30Hzで振幅が\(\frac{1}{\sqrt{2}}\)=(0.7...
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]
つぎにアナログフィルタをデジタルフィルタに変換するために...
$$
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)の置換を行い...
$$
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);
残念ながら数式のままでは表示に失敗したので、代わりに数値...
sageへの入力:
#pre{{
# 振幅特性
# 数式のままではプロットできないので、数値に変換しlist_pl...
#plot(abs(Fz(f)), [f, 0, 50])
list_plot([ N(abs(Fz(f))) for f in range(1,50)],plotjoine...
}}
&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(おもしろかった[5],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
- お世話になります。いま、Octave で信号処理に挑戦しようと...
- 蓮風さま、SageでWavデータ読み込み[[内積を求める例>http:...
- ありがとうございます。いま、scikits などパッケージのイ...
- 蓮風さま、sageへのaudiolabのインストール方法は、以下を...
#pre{{
Sageをシェルモードで起動してpythonのパッケージをインスト...
$ /usr/local/sage-6.0/sage -sh
(sage-sh) $
MacOSXの場合
新たにソースから作成した。
ソースのダウンロード
(sage-sh) $ tar xzvf ~/Downloaded/scikits.audiolab-0.11.0...
(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...
(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
}}
- パッケージもインストールできて、とりあえず相関係数が求...
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2012/03/31からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www.pwv.co.jp:8000/home/pub/22/
私のSageサーバ（http://www15191ue.sakura.ne.jp:8000/）に...
* Sageで信号処理 [#p4db677b]
インターフェース2011年12号にMATLABの特集があったので、Sag...
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...
}}
** エコー混じりの音声データの処理 [#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...
}}
*** 時間軸応答の表示 [#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 respo...
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,...
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...
plot(f1, 20*log10(s1),'m-.'),grid
xlim([0 f1(end)])
title('Periodogram'),xlabel('Frequency(Hz)'),ylabel('Pow...
subplot(2,1,2)
specgram(sig, 512, Fs, hamming(512), 256)
title('Spectrogram'),xlabel('Frequency(Hz)'),ylabel('Tim...
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が...
の２つのバンド・リジェクト・フィルターを組み合わせてエコ...
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では音声信号にフィルターを掛ける関数f...
エコーノイズを除去したdenoise_sigの周波数応答を表示すると...
の出っ張りがとれ、赤い帯も薄くなっています。
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')"%deno...
}}
** フィルターの設計 [#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コマンドでブラウザーから入力した結果がその場で表...
sageへの入力:
#pre{{
# ユーザの入力でフィルターの形状を表示
# fc1:カットオフ周波数1, fc2:カットオフ周波数2, Ft:サンプ...
@interact
def _(type=selector(["butter", "cheby1", "cheby2"]), filt...
#print filter, type, fc1, fc2, Ft, n
octave.eval("[Zb, Za] = MyFilter('%s', '%s', %s, %s, ...
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の...
$$
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に置き換えて周波...
での特性をプロットします。
カットオフ周波数を30Hzで振幅が\(\frac{1}{\sqrt{2}}\)=(0.7...
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]
つぎにアナログフィルタをデジタルフィルタに変換するために...
$$
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)の置換を行い...
$$
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);
残念ながら数式のままでは表示に失敗したので、代わりに数値...
sageへの入力:
#pre{{
# 振幅特性
# 数式のままではプロットできないので、数値に変換しlist_pl...
#plot(abs(Fz(f)), [f, 0, 50])
list_plot([ N(abs(Fz(f))) for f in range(1,50)],plotjoine...
}}
&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(おもしろかった[5],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
- お世話になります。いま、Octave で信号処理に挑戦しようと...
- 蓮風さま、SageでWavデータ読み込み[[内積を求める例>http:...
- ありがとうございます。いま、scikits などパッケージのイ...
- 蓮風さま、sageへのaudiolabのインストール方法は、以下を...
#pre{{
Sageをシェルモードで起動してpythonのパッケージをインスト...
$ /usr/local/sage-6.0/sage -sh
(sage-sh) $
MacOSXの場合
新たにソースから作成した。
ソースのダウンロード
(sage-sh) $ tar xzvf ~/Downloaded/scikits.audiolab-0.11.0...
(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...
(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
}}
- パッケージもインストールできて、とりあえず相関係数が求...
#comment_kcaptcha
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
