フーリエ変換した後のグラフについて

参照信号にディーゼルエンジン音を与えLMSアルゴリズムを適応した後,求められた出力信号と所望信号との誤差を誤差関数としてプロットしました.

下の図の黄色が入力,青が誤差信号,オレンジが所望信号をそれぞれFFTしたグラフです.

このグラフを見ても低減できているか確信を得られないため,なにがどうなっていれば低減できているかおしえていただきたくことはできますでしょうか.

よろしくお願いいたします.

[Kaede]

FFTをしたときは必ず縦軸は0以上ですか？0以上にしたらこのようなグラフになりました.

私は大学4年ですが,まだMATLABを使い始めたばかりなので,上手くいかないことばかりで苦戦しています.みなさんはどのようにMATLABを学習されましたか.

[Shoumei]

周波数特性を見ても、DC（0Hz）成分が大きく、どれがディーゼルエンジンの信号なのかよくわかりません。

元の信号の取得方法に問題ありませんか？どこの周波数帯の信号を低減させてどのような結果を得たいのでしょうか？

MATLABの勉強はひたすらドキュメントを読んで、例を実行・改造してみました。

[Kaede]

コメントありがとうございます.元の信号はwavファイルであるためaudioreadで取得いたしました.

使用しているブロック図はまずxにaをかけた所望信号dとxにLMSアルゴリズムを適応したyの誤差eを出力しています.

そのxに今回はディーゼルエンジン音を入れANCをするという目的です.

プログラムはこちらを使用しておりますが今回はaを単純遅延にさせて行っています.

教えていただきありがとうございます.わからないことばかりで課題が全く進まなくて困っているため一つ一つ理解していこうと思います.頑張ります.

[Shoumei]

wavファイルは何で録音しましたか？

通常、ノイズ除去とか特定周波数の除去などの目的を持って適応フィルタを利用しますが、何のためにディーゼルエンジン音に対して適応フィルタをかけていますか？

もう少し研究室の同僚、先輩方、指導教官にご指導を仰いだほうが良いと思います。ここではMATLABの使い方を教えることはできますが、勉強を教えることはできません。

[Kaede]

大変申し訳ないのですが,何で録音したかまでは知らなくてただデータを渡されてこれを入力としてLMSアルゴリズムを適応してみてくださいと言われただけなのでなんとも言えないです.

研究ではディーゼルエンジン音のノック音の低減をやるための準備としてアクティブノイズコントロールを行っているのだと思います.

はい,かしこまりました.今はオンラインで行っているため聞きづらいところがありますが少しづつやってみます。

[Shoumei]

人間の聴覚で認識するノイズ削減をしたいのであれば、人間の聴覚に近いログスケールで見たほうが良いと思います。表示するときはplotではなく、semilogxを使ってみると良いと思います。

また、プログラムcomplete.mを拝見したところ、元の時間軸データ点数287925点に対して、fft点数が1024点になっていました。そうすると1025点目以降のデータはFFTに含まれないことになりますので、FFT点数は指定せずに使ったほうが良いと思います。

また、適応フィルタというのは、最初は効いておらず、時々刻々とノイズが減ってくるような応答をすると思いますので、FFTで静的な解析をするより、スペクトログラムなどで時間応答を確認したほうが良いと思います。

例があるのでご覧ください。

音響ノイズキャンセリング

正規化された LMS 適応フィルターを使用した、Simulink での低周波ノイズの除去

MATLABのspectrogram関数

プログラムの一部を書き直すと、以下のような特性が得られました。

31, 35, 50Hzの周期ノイズ、残念ながら消えてませんね。

win = hamming(length(d1));
d1win = win.*d1';
e1win = win.*e1';
y1win = win.*y1';

% データは287925点ありますが、FFTを1024点で行うと
% 1025点目以降のデータは切り捨てられます。
% FFT点数を指定せずにデータ点数でDFTを行うようにして使いましょう。
fft_e1=fft(y1);
dB_e1=mag2db(abs(fft_e1));
fft_d1=fft(d1);
dB_d1=mag2db(abs(fft_d1));
fft_y1=fft(y1);
dB_y1=mag2db(abs(fft_y1));

f = (0:length(fft_e1)-1)*fs/length(fft_e1);

figure
semilogx(f,dB_e1)
xlim([0 fs/2])
xlabel('frequency[Hz]')
ylabel('Magnitude[dB]')
hold on
semilogx(f,dB_d1)
semilogx(f,dB_y1)
grid on, hold off
legend('Error signal','Desired signal','Output signal');

あと、目的と理由をちゃんと考えてやりましょうね。