こんなもんはどこの教科書にも書いてあるのだけれど、いつもの教科書置いてきてしまったので手を動かしてみる。PSD のコードも何度か書いた気がするがいつも正規化で迷う。だったらテストケースを持っといた方がいいだろということで。何も新しい事は書いてない。
FFT の定義は \[ Y(k) = \Sigma_{n=0}^{N-1}x(n) \exp{(-2πi n k / N)} \] という普通の。これで例題として \(x(n) = A \sin (2π n/4)\)という時系列を考える。単位は秒。 つまり 0.25 Hz の振動。これを一秒ごとに四回観測してそれぞれ \(x(0) = 0, x(1) = A, x(2) = 0, x(3) = -A\)という観測値を得る。
なんだか重箱だが、観測は「一秒という時間窓の代表」と考えるから観測時間 L = 4 秒であって 3-0 = 3 秒ではない。\(x(1)\)を 0.5 秒から 1.5 秒の観測の代表と考える。
因みに \(x'(n) = A \sin(2π n/2)\)という 0.5 Hz の振動を一秒ごとに評価すると\(x'(n) =0\)しか出てこないのでこれがナイキスト周波数というものかと分かる。
これを \(N=4\) で FFT すると \(Y(0) = 0, Y(1) = -2iA, Y(2) = 0, Y(3) = 2iA\)となる。パワースペクトルは片側とするから正規化していない振幅二乗和は PSD(1) \(= (-2iA)^2+(2iA)^2 = 8A^2\)、PSD(2) \(= 0\)となる。 これを\(N=4\) で割ると元の時系列の振幅二乗和\(2A^2\)に等しいというのが Parseval の公式。
パワースペクトル密度の横軸は観測時間\(L\)の逆数を単位とする周波数で PSD(1) とは PSD(0.25 Hz) で PSD(2) が PSD(0.5 Hz) だから短冊近似したパワーは正規化前は\(0.25 \times 8A^2 = 2A^2\)。
単位も考えておく。振幅が長さ m だと思って元々の一周期の振幅を二乗和して、パワーは力「率」だから一秒あたりとして [m2] 。値は 0.5A2。パワースペクトル密度の横軸が [1/s] だから縦軸の単位は [m2s] 。FFT を見えない\(dt=1\)秒のかかった積分だと思うと単位が [ms] になるからその二乗。となると正規化は分母に\(L\)があると次元が合う。4 で割るとよい、がこれだけでは\(L\)と\(N\)の関係が分からん。
試しに同じ\(x(n)\)を 0, 2/3, 4/3, 2, 8/3, 10/3 秒で観測して\(N=6\)としてみる。時系列は [0,√3A/2,√3A/2,0,-√3A/2,-√3A/2] で二乗和\(3A^2\)に短冊幅 2/3 秒をかけて\(L=4\)で割って前と同じくパワーは \(0.5A^2\)。これを FFT して[0, -3iA, 0, 0, 0, 3iA] で、PSD(0.25Hz)\(=18A^2\)で短冊幅 0.25 Hz をかけて \(4.5A^2\)。正規化するには 9 で割る。
\＃ もう FFT めんどいので Matlab でやった。
今度は同じ\(x(n)\)を 0, 1, 2, 3, 4, 5, 6, 7 で評価して\(L=8\)秒。時系列は [0,A,0,-A,0,A,0,-A] で FFT は [0, 0, -4Ai, 0, 0, 0, 4Ai, 0] で時系列のパワーが\(0.5A^2\)で、正規化前の二乗和が\(32A^2\)で短冊幅 0.125 Hz をかけて\(4A^2\)。正規化するには 8 で割る。
ここまでやると\((L/N^2)\)が見えてくる。あれ、\(L\)が分母だ。つまり FFT に「見えない\(dt=1\)」はかかっていないと解釈すべき。次元は正規化で合わす、と。
\(y(n) = B\sin(2π n / 17)\)を 0, 1/5, 2/5, 3/5, ..., (271-1/5) 秒後で観測して 1355 個の観測値を得る。これの二乗和\(733.536B^2\)に短冊幅 (1/5) をかけて観測時間 271 秒で割ってパワーは\(0.5016B^2\)でちょいと誤差入り。これを FFT したものの絶対値の二乗和は\(921005B^2\)で短冊幅 (1/271) をかけて\(3398.5B^2\)。果たして\(N^2/L\)= 6775 で割ると\(0.5016B^2\)でどんぴしゃ。
当然横軸を周波数\(f\)ではなくて角周波数\(2πf\)にすると\(L/(2πN^2)\)になる。これに関連して注意が必要なのは微分した量が周波数空間で周波数の掛け算になるとき。この時ばかりは\(f\)でなくて\(2πf\)を掛けないとおかしなことになる。