当計測コラムでは、当社お客様相談室によくお問い合わせいただくご質問をとりあげ、回答内容をご紹介しています。今回はFFT解析により得られたパワースペクトルからパワースペクトル密度(PSD)を計算する方法をご紹介します。
ランダム信号などの周期的ではない信号(連続スペクトルをもつ信号)をFFT解析すると、得られたパワースペクトルの各成分の振幅の値はFFT 計算の周波数分解能Δf に依存し、解析条件を変更すると値が変わってしまいます。そのような場合は、解析結果の値がΔfに依存しないように単位周波数幅(1 Hz)あたりのパワー値で表現するパワースペクトル密度(PSD)が使われます。
パワースペクトルと周波数分解能 Δfの関係
FFT解析をおこなう際の周波数レンジ、サンプル点数と周波数分解能の間には次の関係があります。
- ライン数 [点] = サンプル点数 [点] ÷ 2.56
- 周波数分解能Δf [Hz] = 周波数レンジ[Hz] ÷ ライン数 [点]
図1に連続スペクトルをもつ信号をFFT解析して得られたパワースペクトルを示します。周波数レンジはいずれも2 kHzで、サンプル点数は512点、1024点、2048点の三通りです。周波数分解能Δfはそれぞれ10 Hz、5 Hz、2.5 Hzです。
1 kHz成分の振幅値はそれぞれ0.515 m/s2、0.337 m/s2、0.234 m/s2でした。理想的なランダム信号であれば、周波数分解能Δfが半分になると、振幅値は√2 分の1になります。
-
図1 パワースペクトルのサンプル点数による差(上段: 512点/Δf = 10 Hz、中段: 1024点/Δf = 5 Hz、下段: 2048点/Δf = 2.5 Hz)
周波数レンジ: 2 kHz、サンプル点数: 512点でFFT解析をおこなった場合、周波数分解能は10 Hzになり、10 Hz刻みのパワースペクトルが得られます。得られたパワースペクトルの1000 Hz成分は995 Hzから1005 Hzまでの成分の合計値で、1010 Hz成分は1005 Hzから1015 Hzまでの合計値といったように、それぞれの値は周波数分解能の幅を持った帯域の合計値だと考える事ができます。
周波数レンジ: 2 kHz、サンプル点数: 1024点でFFT解析をおこなった場合、周波数分解能は5 Hzになり、合計を求める帯域の幅が半分になりますので、512点で解析した場合と比べるとパワー値(振幅の2乗値)は半分になるため、振幅値は√2分の1になります。
厳密にいうと、512点での解析結果の1000 Hzの値には、995 Hzから1005 Hzまでの成分だけが含まれるのではなく、ウィンドウ関数(ハニング窓、フラットトップ窓)の影響をうけてそれより若干広い範囲の成分が含まれますが、サンプル点数・周波数分解能がかわってもウィンドウ関数の影響は同じであるために、周波数分解能Δfが半分になると振幅値は√2 分の1になるという関係はそのまま成り立ちます。
パワースペクトル密度(PSD)の計算方法
FFT解析で得られたパワースペクトルですが、パワー値(振幅2乗値)が周波数分解能Δfに比例するため、単位周波数幅(1 Hz)あたりのパワー値で表示するには、パワースペクトルのパワー値(振幅2乗値)をΔfで割ればよいことになります。
実際にはΔfで割っただけではウィンドウ関数(ハニング窓、フラットトップ窓)の影響を受けて正しい値にならないため、パワースペクトル密度(PSD)は次の計算式で求めます。
PSD=パワー値÷(Δf×Wf)
=振幅値2÷(Δf×Wf)
ここでWfはウィンドウ(窓関数)ごとの補正値です。計測コラム179号(2016年8月)で、オーバーオール(OA)・パーシャルオーバーオール(POA)を計算する際の補正値Hfをご紹介しましたが、Hfの逆数がパワースペクトル密度を計算する際の補正値Wfになります。
表1 PSDを求める際のウィンドウ(窓関数)の補正値
| ウィンドウ(窓関数) | OA計算時の補正値(Hf) | PSD計算時の補正値(Wf) |
| レクタンギュラ(矩形窓) | 1 | 1 |
| ハニング | 2/3 = 約0.6667 | 3/2 = 1.5 |
| フラットトップ | 1/3.6714416356 = 約0.2724 | 3.6714416356 |
| フォース | 1 | 1 |
| 指数 | 1 | 1 |
パワースペクトル密度(PSD)はパワー値(振幅2乗値)を周波数分解能Δfで割ったもので、その単位は[振幅値の単位2 / Hz]という単位になります。振幅値の単位がm/s2 であれば[(m/s2)2 / Hz]、振幅値の単位が Vであれば[V2 / Hz]という単位になります。
パワースペクトル密度(PSD)の値は[振幅値の単位2 / Hz]という単位で表示する事が多いのですが、その平方根をとった値で表示する事もあります。この場合、PSDの平方根の値の単位は[m/s2 / √Hz]、[V / √Hz]などと表現されます。過去に測定した結果や、測定対象製品の仕様の比較する場合、それらのPSDの値の単位が[振幅値の単位2 / Hz]なのか[振幅値の単位/ √Hz]なのかを確認し、単位をあわせて測定する必要があります。
当社FFTアナライザには、パワースペクトル密度(PSD)を表示する機能があります。その際は表示値をV2かVから選ぶことができ、V2を選ぶとPSDの値が[振幅値の単位2 / Hz]という単位で表示されます。Vを選ぶとPSDの平方根の値が[振幅値の単位/ √Hz]という単位で表示されます。
加速度パワースペクトルからのパワースペクトル密度の算出
FFTアナライザで解析した加速度パワースペクトルをエクセルに読み込んだ例を表2に示します。周波数レンジ(C6セル)は2000 Hz、サンプル点数(B7セル)は1024点ですので、周波数分解能は5 Hzになります。A17 ~ A417セルが周波数の値で0 Hzから2 kHzまで5 Hzおきの値が並んでいます。
Y軸スケール(B14セル)がLinになっていますので、これはY軸スケールをLinに設定して計測したデータで、B17 ~ B417セルの値は各周波数成分の物理値(加速度値)です。Y軸スケールがMagLogになっているデータはY軸スケールをLog/MagLogに設定して計測したデータで、この場合もB17 ~ B417セルの値は各周波数成分の振幅値(加速度値)です。
D17 ~ D417セルに表2のような数式を入力すると同セルにパワースペクトル密度(PSD)の値が表示されます。この値の単位は[(m/s2)2 / Hz]です。
Y軸Magnitude(B16セル)がrmsになっていますのでB17 ~ B417セルの値は実効値(RMS値)です。そのため、上記の方法で求めたパワースペクトル密度(PSD)の値も実効値(RMS値)です。
なお、表2の内容は以下のリンク からダウンロードいただけます。
表2 加速度パワースペクトルからのパワースペクトル密度の算出例
表2 加速度パワースペクトルからのパワースペクトル密度の算出例
| A | B | C | D | E | |
| 1 | Label: | CH2: パワースペクトル | |||
| 2 | DateTime: | Tue Oct 18 20:35:50 2016 | |||
| 3 | DataKind: | CH2 | PowerSpec | Mag | |
| 4 | DataPoints: | 402 | Filter: | FLAT | : |
| 5 | DataCalc: | ||||
| 6 | Frequency: | 0 | 2000 | Hz | |
| 7 | Sample: | 1024 | Internal | ||
| 8 | Average: | 681 | Power/Sum | ||
| 9 | Voltage(CH2): | -10 | dBVrms | ||
| 10 | EU/V(CH2): | 1.00E+03 | 0dBRef.(CH2): | 1.00E+00 | |
| 11 | Window(CH2): | Hann | |||
| 12 | X-AxisScale: | Lin | 周波数分解能[Hz] | ||
| 13 | X-AxisUnit: | Hz | 5 | =C6/(B7/2.56) | |
| 14 | Y-AxisScale: | Lin | |||
| 15 | Y-AxisUnit: | m/s2 | PSD値 | セルの数式 | |
| 16 | Y-AxisMagnitude: | rms | (m/s2)^2/Hz | ||
| 17 | 0.0 | 0.211 | 0.005917 | =(B17*B17)/($D$13*1.5) | |
| 18 | 5.0 | 0.151 | 0.003042 | =(B18*B18)/($D$13*1.5) | |
| 19 | 10.0 | 0.077 | 0.000800 | =(B19*B19)/($D$13*1.5) | |
| 20 | 15.0 | 0.173 | 0.003980 | =(B20*B20)/($D$13*1.5) | |
| 21 | 20.0 | 0.201 | 0.005400 | =(B21*B21)/($D$13*1.5) | |
| 22 | 25.0 | 0.222 | 0.006581 | =(B22*B22)/($D$13*1.5) | |
| 23 | 30.0 | 0.245 | 0.008021 | =(B23*B23)/($D$13*1.5) | |
| … | |||||
電圧信号パワースペクトルからのパワースペクトル密度の算出
FFTアナライザで解析してデシベル値で表示した電圧信号のパワースペクトルをエクセルに読み込んだ例を表3に示します。周波数レンジ(C6セル)は2000 Hz、サンプル点数(B7セル)は1024点ですので、周波数分解能は5 Hzになります。A17 ~ A417セルが周波数の値で0 Hzから2 kHzまで5 Hzおきの値が並んでいます。
Y軸スケール(B14セル)がLogになっていますので、これはY軸スケールをLogに設定して計測したデータで、B17 ~ B417セルの値は各周波数成分のデシベル値です。
表3 電圧信号スペクトルからのオーバーオールの算出例
| A | B | C | D | E | |
| 1 | Label: | CH2: パワースペクトル | |||
| 2 | DateTime: | Tue Oct 18 20:37:00 2016 | |||
| 3 | DataKind: | CH2 | PowerSpec | Mag | |
| 4 | DataPoints: | 402 | Filter: | FLAT | AnalogFilter(CH2): |
| 5 | DataCalc: | ||||
| 6 | Frequency: | 0 | 2000 | Hz | |
| 7 | Sample: | 1024 | Internal | ||
| 8 | Average: | 681 | Power/Sum | ||
| 9 | Voltage(CH2): | -10 | dBVrms | ||
| 10 | EU/V(CH2): | 1.00E+00 | 0dBRef.(CH2): | 1.00E+00 | ADOverHold(CH2): |
| 11 | Window(CH2): | Hann | |||
| 12 | X-AxisScale: | Lin | 周波数分解能[Hz] | ||
| 13 | X-AxisUnit: | Hz | 5 | =C6/(B7/2.56) | |
| 14 | Y-AxisScale: | Log | |||
| 15 | Y-AxisUnit: | V | PSD値 | セルの数式 | |
| 16 | Y-AxisMagnitude: | rms | V2/Hz | ||
| 17 | 0.0 | -73.53 | 5.9171E-09 | =(10^(B17/10))/($D$13*1.5) | |
| 18 | 5.0 | -76.42 | 3.0421E-09 | =(10^(B18/10))/($D$13*1.5) | |
| 19 | 10.0 | -82.22 | 7.9964E-10 | =(10^(B19/10))/($D$13*1.5) | |
| 20 | 15.0 | -75.25 | 3.9802E-09 | =(10^(B20/10))/($D$13*1.5) | |
| 21 | 20.0 | -73.93 | 5.3996E-09 | =(10^(B21/10))/($D$13*1.5) | |
| 22 | 25.0 | -73.07 | 6.5810E-09 | =(10^(B22/10))/($D$13*1.5) | |
| 23 | 30.0 | -72.21 | 8.0211E-09 | =(10^(B23/10))/($D$13*1.5) | |
| … | |||||
D17 ~ D417セルに表3のような数式を入力すると同セルにパワースペクトル密度(PSD)の値が表示されます。この値の単位は[V2 / Hz]です。デシベル値を振幅値に変換する式は10^(デシベル値/20) ですが、今回はデシベル値をパワー値(振幅2乗値)に変換したいので、10^(デシベル値/10) という式を使いました。
Y軸Magnitude(B16セル)がrmsになっていますのでB17 ~ B417セルの値は実効値(RMS値)です。そのため、上記の方法で求めたパワースペクトル密度(PSD)の値も実効値(RMS値)です。
なお、表3の内容は以下のリンクからダウンロードいただけます。
まとめ
今回はFFT解析により得られたパワースペクトルからパワースペクトル密度(PSD)を計算する方法をご紹介しました。
当社FFTアナライザにはパワースペクトル密度(PSD)を計算・表示する機能がありますが、その機能を使わずに保存した場合などは今回紹介した方法で計算する事ができます。
パワースペクトル密度(PSD)は1 Hzあたりのパワー値(振幅値の2乗値)で表示され、その単位は[振幅値の単位2 / Hz] すなわち[(m/s2)2 / Hz]、[V2 / Hz]などになります。
FFT解析をおこない、過去に測定した結果や測定対象製品の仕様の比較する場合、過去の結果や仕様がPSDなのか単なるパワースペクトルなのか、PSDであれば値の単位が[振幅値の単位2 / Hz]なのか[振幅値の単位/ √Hz]なのかを確認したうえで、単位や測定条件をあわせて測定をおこなってください。
(2016年10月20日発行メールマガジンより抜粋)