FrontPage

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

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

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

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

  • audio
  • signal

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

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

fig1.png

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

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

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

fig2.png

よっと分かりずらいですが、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);
% 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);    % プロット

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

fig3.png

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

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

fig4.png

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

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

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

エコーノイズを除去した結果:filedenoised.wav

コメント

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

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

  • いくつかの図が抜けていました。申し訳ありません。 -- 竹本 浩? 2012-03-18 (日) 21:00:25

(Input image string)


添付ファイル: filefig3.png 528件 [詳細] filefig2.png 595件 [詳細] filedenoised.wav 578件 [詳細] filefig4.png 531件 [詳細] filevoice_howling.wav 936件 [詳細] filefig1.png 543件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2016-12-14 (水) 21:44:41 (160d)
SmartDoc