MATLABで学ぶ信号処理　 高速フーリエ変換（fast Fourier transform:fft）について 1

2021.03.20

MATLABにはfft関数が用意されていて、簡単に誰でもfftができます。
しかし、振動・騒音関係でMATLABのfftを使用する場合は注意が必要です。

スポンサーリンク

どういうことかというと、騒音器や振動測定器で測定した時間波形をMATLABやPythonなどのプログラムで単純にfftをすると、騒音器や振動測定器で測定した”周波数分析結果”と一致しません。

今回は測定器の周波数分析結果と自作プログラムのFFT結果が異なる原因を説明したいと思います。
また、騒音器や振動測定器と同様の結果が得られるようにFFTのプログラムを作成するにはどうすればいいかを紹介いたします。

学生の方も読んでいるかも知れないので、順を追って説明いたします。

騒音や振動は基本的に実効値として扱うことになっています。
例えば、騒音の場合、騒音[dB]は下のような式で表されます。

このとき、マイクロフォン校正器と呼ばれる、マイクロフォンでの測定が正しくなるような調整機器のようなものでは、1.0k[Hz]の1[Pa]の正弦波が出力されます。これを測定して、時間波形を観測すると、下記のようになります。

スポンサーリンク

「校正器からは1.0k[Hz]の1[Pa]の正弦波が出力される」と言いましたが、実際には1[Pa]以上であることがわかります。これは、実効値で1[Pa]の正弦波が出力されていることが原因です。

fs=2^13;
fr=1000;
t=0:1/fs:1;
sig=sin(2*pi*fr*t)*sqrt(2);

  • MATLABのfftについて

では、上記の時間波形をMATLABで単純にfftした場合、どのようになるのでしょうか？

clear all;close all

fs=2^13;
fr=1000;
t=0:1/fs:1;
sig=sin(2*pi*fr*t)*sqrt(2);

p=fft(sig,length(sig));

freq=0:length(p)-1;
figure
subplot(211)
plot(freq,abs(p))
xlabel('Freq. [Hz]');ylabel('Ampli. [Pa]')
subplot(212)
plot(freq,angle(p))
xlabel('Freq. [Hz]');ylabel('Phase [rad]')

結果はこちらです。

振幅が1.0 [Pa]ではないこと、7193Hz(厳密には fs-1000+1 Hz)にもピークがでることがわかります。
つまり、何かがおかしいってことです。
では、何がおかしいのかについて説明していきます。

  • 原因①　虚像の理解不足


 

  • 原因②　fft窓の理解不足


  • 正しいfftプログラム

では、正しいプログラムを下記に示します。

スポンサーリンク 

clear all;close all

fs=2^13;
fr=1000;
t=0:1/fs:1;
sig=sin(2*pi*fr*t)*sqrt(2);


figure
subplot(211)
plot(t,sig)
xlabel('Time [s]')
ylabel('Pressure [Pa]')
subplot(212)
plot(t,sig)
xlim([0 t(20)])
xlabel('Time [s]')
ylabel('Pressure [Pa]')

p=fft(sig,length(sig))/length(sig)/sqrt(2);
p=[p(1) p(2:floor((length(p)/2)))*2 p(floor((length(p)/2))+1)];
freq=0:fs/2;
figure
subplot(211)
plot(freq,abs(p(1:length(freq))))
xlabel('Freq. [Hz]')
ylabel('Pressure [Pa]')

subplot(212)
plot(freq,angle(p))
xlabel('Freq. [Hz]')
ylabel('Phase [rad]')

結果はこちらです。1k[Hz]で1[Pa]となっていることがわかります。

 

今回はこのへんでGood luck

コメント