Sageで信号処理

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

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

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

sageでoctaveを使う場合には、Sageのサーバマシンにoctaveと必要なパッケージがインストールされている 必要があります。 「octaveのインストール」参照

octave連携の準備

sageからoctaveを使う場合のグラフ処理を助ける関数を定義します。 （今回の改造で、グラフのアップデートに対応しました）

エコー混じりの音声データの処理

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

最初にデータの読み込み

最初にエコー混じりの音声データvoice_howling.wavを読み込みます。

時間軸応答の表示

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

周波数特性の表示

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+ファイル名で取り出すことができます。

フィルターの作成

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

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);    % プロット	

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

エコーノイズの除去

出来上がったフィルターでエコーノイズを除去してみましょう。 信号処理に強いoctaveでは音声信号にフィルターを掛ける関数filterが用意されています。

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

結果の保存

ノイズの取れた音声を、denoised.wavに保存します。 再生すると、ピーという音が小さくなっているのが分かります。

フィルターの設計

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の種類を変えて試してみて下さい。

双一次z変換

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

デジタルフィルタの周波数特性

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

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

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

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

sageは簡単

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

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