octave/信号処理（エコーノイズの除去）のバックアップの現在との差分(No.11)

追加された行はこの色です。
削除された行はこの色です。
[[FrontPage]]

#contents

2012/03/18からのアクセス回数 &counter;

** エコーノイズの除去をoctaveにやってみる [#tb637dd3]
Interface 2011/12号のMATLABの特集で紹介されていた信号処理の例題（エコーノイズの除去）をoctaveでやってみました。

octaveのパッケージとして、以下のものが必要です。
- audio
- signal

*** 音声+エコーノイズファイルの読み込みと再生 [#s2bd2658]
&ref(voice_howling.wav); を読み込み、再生してみます。

#pre{{
% WAVファイルをロードして時間軸応答を表示します。
[noised_sig,Fs] = wavread('voice_howling.wav');
sound(noised_sig,Fs);
}}
結構すごいハウリングが聞こえてきます。

*** ノイズの時間軸応答の表示 [#y7c2a6b7]
特集ではMATLABを使っているのですが、octaveに存在しない関数があり、すべてを実行することはできませんでした。

今回確認できたのは、以下の3つのグラフです。
- 時間応答のプロット
- ピリオドグラム(periodogram)
- スペクトグラム(specgram)

#pre{{
% 時間軸応答の表示
[noised_sig,Fs] = wavread('voice_howling.wav');
t = (0:length(noised_sig)-1)/Fs;
figure
plot(t,noised_sig),grid on
xlim([0 t(end)]),ylim([-1 1])
title('voice+howling noize time domain response'),xlabel('Time(sec)'),ylabel('Magnitude')
}}

&ref(fig1.png);

*** 音声＋ノイズのピリオドグラムとスペクトグラム [#oc442d88]
次に音声＋ノイズのピリオドグラムとスペクトグラムを表示してみます。

#pre{{
% ノイズ信号の周波数応答の表示
figure
subplot(2,1,1)
[s1, f1] = periodogram(noised_sig, hamming(length(noised_sig)), length(noised_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(noised_sig, 512, Fs, hamming(512), 256)
title('Spectrogram'),xlabel('Frequency(Hz)'),ylabel('Time(sec)')
}}

&ref(fig2.png);

よっと分かりずらいですが、2850Hzと8500Hzにピリオドグラムの突出した信号が見られます。またスペクトグラムだと2850Hzと8500Hzに2本の濃い赤の帯が見て取れます。

** エコーノイズの除去 [#d1dd7869]

エコーを除去するために、以下の様なBRFフィルタを作成します。
| 応答タイプ | BRF |
| 設計手法 | IIR BRF 縦結合 |
| サンプリング周波数Fs(Hz) | 22050 |
| カット周波数 | 2842.4, 8527.1 |
| 帯域幅(Hz) | 100 |

このフィルターを作るために[[octave/IIRフィルタの伝達関数設計]]を勉強しました。

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

フィルタの周波数特性は以下の様になります。

&ref(fig3.png);

最後に、音声+ノイズのデータにフィルタを掛けてエコーノイズを除去します。

#pre{{
%% エコーノイズの除去
denoise_sig = filter(Zb, Za, noised_sig);

figure
subplot(2,1,1)
[s1, f1] = periodogram(denoise_sig, hamming(length(noised_sig)), length(noised_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(denoise_sig, 512, Fs, hamming(512), 256)
title('Spectrogram'),xlabel('Frequency(Hz)'),ylabel('Time(sec)')
}}

&ref(fig4.png);

濃い赤の帯が細く、少し薄くなっています。

再生してみると、完全にノイズが除去されたわけではありませんが、ハウリングのキーンという音は小さくなっています。

#pre{{
% 再生と結果の保存
sound(denoise_sig, Fs);
wavwrite(denoise_sig, Fs, 'denoised.wav');
}}

エコーノイズを除去した結果：&ref(denoised.wav);

** コメント [#y184975c]
#vote(おもしろかった[7],そうでもない[0],わかりずらい[1])
#vote(おもしろかった[29],そうでもない[1],わかりずらい[2])

皆様のご意見、ご希望をお待ちしております。
- いくつかの図が抜けていました。申し訳ありません。 -- [[竹本　浩]] &new{2012-03-18 (日) 21:00:25};

#comment_kcaptcha
octave/信号処理（エコーノイズの除去） のバックアップの現在との差分(No.11)

octave/信号処理（エコーノイズの除去）のバックアップの現在との差分(No.11)
