MATLABで学ぶ信号処理 高速フーリエ変換（fast Fourier transform:fft）について 1
はじめに
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
MATLABプログラム以外にも機械エンジニアに有益な情報を発信していますので、ご興味のある方はぜひこちらもご確認いただけると幸いです。
Discussion