前回の記事
前回の記事でパワースペクトルの内容について少しだけ記載しました。
スポンサーリンク
質問を頂きました
この内容について、大手企業に勤めているエンジニアの方から質問を頂きました。
振動音響用の市販のFFTアナライザでオートパワーもしくはスペクトルを計測して、その結果を利用してオクターブ分析やO.A.値の計算をすると結果が一致しないのですが、なぜでしょうか？
私の推測も交えながら、この質問を意図をもう少し詳細に記載したいと思います。
このエンジニアの方は「市販のFFTアナライザ」を使用しているとのことなので、小野測器のDSシリーズ、シーメンスのScadas、Brüel & KjærのPulseなどの計測器を使用していると考えられます。
その計測器でオートパワーもしくはスペクトルを計測している。また、オクターブ結果やO.A.結果も測定・計算している、ということだと思います。
そして、計測器で測定・計算したオクターブ結果やO.A.結果と、自分でオートパワーから計算したオクターブ結果やO.A.結果が一致しない、ということで質問してきたのだと思われます。
本記事だけでは不十分な場合は、下記文献が参考になると思います。
実効値と振幅値について
FFTで何をみたいの？（実効値なのか、振幅値なのか）
オートパワーの説明の前に、そもそもFFT（スペクトル）で何を求めたいのでしょうか？
自分自身でそのことをちゃんと理解しているのでしょうか？
例えば、図1のような信号があったときに、振幅値Aを求めたいのか、実効値（RMS）A／√2を求めたいのか、どっちでしょうか？（振動や騒音であれば、実効値で求める必要があります。）
また、実測で図1の信号を取得した場合は、センサーの校正値は実効値を測定するための値なのか、振幅値を測定するための値なのか、を理解しておく必要があります。（ほとんどの校正値は実効値ではなく振幅値を求めるための値です。）
図1 信号
オートパワーは実効値なのか、振幅値なのか
スペクトルの説明と同様に、実効値を求めたいのか、振幅値を求めたいのかを明確にしておく必要があります。
重要：窓関数を使用した場合はオートパワーやスペクトルは2種類ある！
窓関数を使用した場合、オートパワーやスペクトルは2種類あります。こんな説明をすると、先生方や研究者の方に「何言ってんの？」と言われてしまいそうですね。順を追って説明したいと思います。
先ほどの繰り返しになりますが、振動や騒音ではオートパワーやスペクトルでは実効値を扱います。ここは覚えておいてください。なので、スペクトル、オートパワーおよびPSDには下表のような関係があります。
| 振幅 | 実効値（の振幅） | |
| スペクトル | A | A/√2 |
| オートパワー | A^2 | A^2/2 |
| PSD | A^2/Δf | A^2/(Δf×2) |
ここまでは、みなさん問題ないと思います。
では、もう一度同じ質問をします。
そもそもFFT（スペクトル）やオートパワーで何を求めたいのでしょうか？
ほとんど下記の①か②に該当すると思います。
①オートパワーやスペクトルで実効値の振幅を正確に求めたい
②オートパワーやスペクトルからオクターブ分析やO.A.値を求めたい
例えばハニング窓を使用してスペクトル＆オートパワーを求めると、下記のようになります。ちなみに、青線が窓関数なし、赤線がハニング窓を適用しています。窓により振幅が小さくなる影響は補正していません。
この窓関数適用の影響を補正することになるのですが、補正して何をみたいのでしょうか？
スポンサーリンク
close all
fs=2^10*32;
t=0:1/fs:1;
hannWin = 0.5*ones(1,length(t)) - 0.5*cos(2.0*pi*t);
sig=sqrt(2)*sin(1000*pi*2*t);
% % fft
fft_sig=fft(sig,fs)/length(sig)/sqrt(2);
win_fft_sig=fft(sig.*hannWin,fs)/length(sig)/sqrt(2);
% % 虚像（負の周波数の補正）
fft_sig=fft_sig(1:fs/2)*2;
win_fft_sig=win_fft_sig(1:fs/2)*2;
figure
subplot(211)
plot(0:length(fft_sig)-1,abs(fft_sig))
hold on
plot(0:length(win_fft_sig)-1,abs(win_fft_sig),'r')
xlim([980 1020])
grid on
title('スペクトル')
xlabel('周波数 Hz')
ylabel('振幅 Pa')
subplot(212)
plot(0:length(fft_sig)-1,abs(fft_sig).^2)
hold on
plot(0:length(win_fft_sig)-1,abs(win_fft_sig).^2,'r')
xlim([980 1020])
grid on
title('オートパワー')
xlabel('周波数 Hz')
ylabel('パワー Pa^2')
①オートパワーやスペクトルで実効値の振幅を正確に求めたい
オートパワーやスペクトルで”実効値の振幅”を正確に求めたい場合、ハニング窓だとfft結果を「2倍」する必要があります。すると、”実効値の振幅”が一致します。
close all
fs=2^10*32;
t=0:1/fs:1;
hannWin = 0.5*ones(1,length(t)) - 0.5*cos(2.0*pi*t);
sig=sqrt(2)*sin(1000*pi*2*t);
% % fft
fft_sig=fft(sig,fs)/length(sig)/sqrt(2);
win_fft_sig=fft(sig.*hannWin,fs)/length(sig)/sqrt(2);
% % 虚像（負の周波数の補正）
fft_sig=fft_sig(1:fs/2)*2;
win_fft_sig=win_fft_sig(1:fs/2)*2;
% % ①の補正
win_fft_sig=win_fft_sig*2;
figure
subplot(211)
plot(0:length(fft_sig)-1,abs(fft_sig))
hold on
plot(0:length(win_fft_sig)-1,abs(win_fft_sig),'r')
xlim([980 1020])
grid on
title('スペクトル')
xlabel('周波数 Hz')
ylabel('振幅 Pa')
subplot(212)
plot(0:length(fft_sig)-1,abs(fft_sig).^2)
hold on
plot(0:length(win_fft_sig)-1,abs(win_fft_sig).^2,'r')
xlim([980 1020])
grid on
title('オートパワー')
xlabel('周波数 Hz')
ylabel('パワー Pa^2')
②オートパワーやスペクトルからオクターブ分析やO.A.値を求めたい
オートパワーやスペクトル結果からオクターブ分析結果やO.A.値を求めたい場合、周波数帯域内のエネルギーが一致している必要があります。そのため、ハニング窓だとfft結果を「√(8/3)倍」する必要があります。（オートパワーを8/3倍）
スポンサーリンク
close all
fs=2^10*32;
t=0:1/fs:1;
hannWin = 0.5*ones(1,length(t)) - 0.5*cos(2.0*pi*t);
sig=sqrt(2)*sin(1000*pi*2*t);
% % fft
fft_sig=fft(sig,fs)/length(sig)/sqrt(2);
win_fft_sig=fft(sig.*hannWin,fs)/length(sig)/sqrt(2);
% % 虚像（負の周波数の補正）
fft_sig=fft_sig(1:fs/2)*2;
win_fft_sig=win_fft_sig(1:fs/2)*2;
% % ②の補正
win_fft_sig=win_fft_sig*sqrt(8/3);
figure
subplot(211)
plot(0:length(fft_sig)-1,abs(fft_sig))
hold on
plot(0:length(win_fft_sig)-1,abs(win_fft_sig),'r')
xlim([980 1020])
grid on
title('スペクトル')
xlabel('周波数 Hz')
ylabel('振幅 Pa')
subplot(212)
plot(0:length(fft_sig)-1,abs(fft_sig).^2)
hold on
plot(0:length(win_fft_sig)-1,abs(win_fft_sig).^2,'r')
xlim([980 1020])
grid on
title('オートパワー')
xlabel('周波数 Hz')
ylabel('パワー Pa^2')
つまり、①は振幅を正確求めるための補正値で、②は帯域のエネルギーを求めるための補正値です。
一般的に、市販の測定器でスペクトルやオートパワーで測定した場合、補正値は①を使用されます。そのため、そのスペクトルやオートパワーを用いて、ご自身でO.A.値やオクターブ分析結果を求めても、①の補正値を使用したままなので、計算結果は間違ってしまいます。
市販の測定器で直接測定したO.A.値やオクターブ分析結果は②の補正値を使用して計算されています。
この理解がないままにご自身でプログラムを組んでしまうとお客さんに報告した結果が実は間違っていたということになりかねません。みなさん注意しましょう！特にJISやISOの試験で、分析をMATLABで行う場合は要注意です。計算が間違っていて、要求仕様を満たせない場合……….考えたくないですね。
教科書や論文ではここらへんのことは詳しく説明してくれません。計測器メーカーの人でも知らない人が多いですし、質問してもわかってないので、トンチンカンな回答をしてくることもあります。困ったものです。
ちなみに、最近流行りのデータサイエンスさん達はこの辺の知識が皆無です。
なので、機械の故障などの判定に振動・騒音を使用していて、閾値で判定している場合、その閾値の計算が間違っていたりします。まあ故障を判定できているなら問題ないですが。
今後、振動・騒音のfftやFRFを算出するfunctionファイルを提供していく予定なので、気になる方はこのサイトをチェックして頂けると幸いです。
今日はこのへんでGood luck!
コメント