FrontPage

2012/03/31からのアクセス回数 3437

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

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

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

Sageで信号処理

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

sageの特徴は、以下の3点だと思います。

  • ブラウザーでどこからでも使える
  • 実行可能なドキュメントである
  • 既存ツールとの連携が可能

今回の例で、octaveとの連携および実行可能なドキュメントであることを余すとこなくご紹介します。

sageでoctaveを使う場合には、Sageのサーバマシンにoctaveと必要なパッケージがインストールされている 必要があります。 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を使って試してみましょう。

最初にデータの読み込み

最初にエコー混じりの音声データfilevoice_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')

fig1.png

周波数特性の表示

octave.evalでoctaveのコマンドを毎回括るのは大変ですし、ミスを多くなります。

そこで、octaveの処理をファイルにまとめて、sageのワークシートにアップロードして、 このファイルを利用すると便利です。

ファイルのアップロードは、

topmenu.png

の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')

fig2.png

フィルターの作成

エコーを除去するために、以下の様なBRFフィルタを作成します。

応答タイプBRF
設計手法IIR BRF 縦結合
サンプリング周波数Fs(Hz)22050
カット周波数2842.4, 8527.1
帯域幅(Hz)100

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

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);
% 5次のアナログ基準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')

fig3.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')

fig4.png

結果の保存

ノイズの取れた音声を、filedenoised.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')

fig5.png

デジタルフィルタの特性

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

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

fig6.png

アナログフィルタ

最初に100Hzのサンプリング周波数でカットオフ周波数を30Hzの1次のアナログLPF $$ G(s) = \frac{1}{1+ \frac{s}{6o \pi}} $$ の周波数特性をsageで表示してみます(図13の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])

fig7.png

双一次z変換

つぎにアナログフィルタをデジタルフィルタに変換するために、以下の双一次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)の置換を行い、周波数特性を表示します(図13の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)

fig8.png

残念ながら数式のままでは表示に失敗したので、代わりに数値に変換し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)

fig9.png

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

sageへの入力:

# 位相特性
f = var('f')
plot(angle(Fz(f)), [f, 0, 50])

fig10.png

sageは簡単

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

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

コメント

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

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

  • お世話になります。いま、Octave で信号処理に挑戦しようとしているのですが、共同研究者がプログラミングに慣れておらず、できるだけの工程を Sage で数式化できないものか、と考えております。いま、取り組むべきは相関関数の算出です。そこで音声ファイルをこのページのように読み込んでみたものの、Sage の vector にコンヴァートし内積を計算するところで、無限ループ化してしまいます。Octave では所定の操作を完遂できるのですが、それが数式で表現した場合どうなるのか提示したいのです。お手数ですが wavread() で読み込んだ『ステレオ音声』を L チャンネル vector、R チャンネル vector に切り分けて内積を取るには如何したら宜しいでしょうか? 以上、宜しくお願いします。 -- 蓮風? 2014-09-22 (月) 11:41:51
  • 蓮風さま、SageでWavデータ読み込み内積を求める例に置きました。参考にして下さい。 -- 竹本 浩? 2014-09-22 (月) 17:13:17
  • ありがとうございます。いま、scikits などパッケージのインストールから試みています。がんばります。 -- 蓮風? 2014-09-22 (月) 19:18:17
  • 蓮風さま、sageへのaudiolabのインストール方法は、以下を参考にしてください。 -- 竹本 浩? 2014-09-22 (月) 19:40:28
    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
    
  • パッケージもインストールできて、とりあえず相関係数が求められました。共同研究者に報告できます。ありがとうございます!! -- 蓮風? 2014-09-23 (火) 00:49:55

(Input image string)


添付ファイル: filedenoised.wav 316件 [詳細] filevoice_howling.wav 357件 [詳細] filetopmenu.png 514件 [詳細] filefig10.png 479件 [詳細] filefig9.png 460件 [詳細] filefig8.png 437件 [詳細] filefig7.png 472件 [詳細] filefig6.png 461件 [詳細] filefig5.png 458件 [詳細] filefig4.png 479件 [詳細] filefig3.png 467件 [詳細] filefig2.png 508件 [詳細] filefig1.png 472件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2014-09-23 (火) 00:49:55 (1010d)
SmartDoc