octave/信号処理（エコーノイズの除去）のバックアップ(No.2) - PukiWiki

[ トップ ] [ 新規 | 一覧 | 単語検索 | 最終更新 | ヘルプ ]

エコーノイズの除去をoctaveにやってみる
エコーノイズの除去
コメント

2012/03/18からのアクセス回数 15270

エコーノイズの除去をoctaveにやってみる †

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

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

audio
signal

音声+エコーノイズファイルの読み込みと再生 †

voice_howling.wav を読み込み、再生してみます。

% WAVファイルをロードして時間軸応答を表示します。
[noised_sig,Fs] = wavread('voice_howling.wav');
sound(noised_sig,Fs);

結構すごいハウリングが聞こえてきます。

ノイズの時間軸応答の表示 †

特集ではMATLABを使っているのですが、octaveに存在しない関数があり、すべてを実行することはできませんでした。

今回確認できたのは、以下の3つのグラフです。

時間応答のプロット
ピリオドグラム(periodogram)
スペクトグラム(specgram)

% 時間軸応答の表示
[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')

音声＋ノイズのピリオドグラムとスペクトグラム †

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

% ノイズ信号の周波数応答の表示
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)')

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

エコーノイズの除去 †

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

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

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

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を作成する
[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);    % プロット

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

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

%% エコーノイズの除去
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)')

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

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

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

エコーノイズを除去した結果：denoised.wav

コメント †

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

SmartDoc