はじめに
自己相関関数とフーリエ変換には密接な関係があり、パワースペクトルのフーリエ逆変換は自己相関と一致することが知られています。
そのため、データが多いときは定義通り自己相関関数を計算をするよりも、フーリエ変換を使って自己相関関数を計算する方が早い場合も多いです。
ですが、我々が数値計算できるのは離散化されたフーリエ変換です（データも離散化されているし…)。
そこで、離散化されたデータ、フーリエ変換でも自己相関関数とパワースペクトルの関係は成り立つのか？ということを数式を用いて考えます。
先に結果だけ述べておくと離散フーリエ変換から自己相関を計算するときには、自己相関の定義に気を付ける必要があるということがわかります。
普通の連続的なフーリエ変換の場合
普通のフーリエ変換とは以下の通り
\hat{f}(\omega) = \int ^\infty_{-\infty} f(t)e^{-i\omega t} dt \\
f(t) = \frac{1}{2\pi} \int ^\infty_{-\infty} \hat{f} (\omega) e^{i\omega t} dt \ .
定義の仕方はいろいろあると思うが今回は上記のものを使う。
このようなフーリエ変換の場合、自己相関関数を$\xi (t)$として
\begin{align}
\xi (t) &\equiv f^*(-t)*f(t) \\
&= \int^\infty_{-\infty}f^*(t'-t)f(t')dt' = \ \int^\infty_{-\infty}f^*(t')f(t'+t)dt' \\
&= \frac{1}{2\pi}\int^\infty_{-\infty} P(\omega)e^{i\omega t} d\omega
\end{align}
が成り立つ。ここで$P(\omega)$はパワースペクトル。
これはウィーナー-ヒンチン(Wiener-Khinchin)の定理と呼ばれ、有名な定理なので証明は割愛させてもらいます。
この式を見るとパワースペクトルのフーリエ逆変換が自己相関関数になることをがわかります。
これを離散フーリエ変換でやろうというのがこの記事の主題です。
離散フーリエ変換の場合
まず離散フーリエ変換の定義ですが、scipy.fftで使われる定義に従うことにする。UserGuideを見てみるとフーリエ変換とその逆変換は
y[k]=\sum^{N-1}_{n=0}e^{-2\pi j \frac{kn}{N}}x[n] \\
x[n]=\frac{1}{N}\sum^{N-1}_{k=0} e^{2\pi j \frac{kn}{N}} y[k]
という定義で書かれている。そのためこれをもとに計算していこうと思うのだが、いささか見ずらい。少しわかりやすく文字を変えて、
\hat{f}[k]=\sum^{N-1}_{n=0}e^{-2\pi i \frac{kn}{N}}f[n] \\
f[n]=\frac{1}{N}\sum^{N-1}_{k=0} e^{2\pi i \frac{kn}{N}} \hat{f}[k]
の形にして計算をしていく。
虚数を$j$から$i$にしたため、流儀の違う方に怒られそうだがご勘弁を。
さて計算していこう…と思うのだが、計算に必要なものがある。
それはデルタ関数だ。先ほどのウィーナー-ヒンチンの定理を証明するには、実はデルタ関数のフーリエ表示が必要である。
今回は離散データを考えてディラックのデルタ関数ではなく、クロネッカーのデルタを離散フーリエ表示することによって計算を進める。
クロネッカーのデルタの離散フーリエ表示
まずクロネッカーのデルタは以下のような定義ができるはずだ、
g[x] = \sum^{N-1}_{x'=0} g[x']\delta_{x,x'} \ .
ここで$g[x]$は有限の長さ$N$を持つ任意の関数である。
$\delta_{x,x'}$が、$x=x'$の場合に$\delta_{x,x'}=1$で、それ以外の場合に$\delta_{x,x'}=0$となるような関数であれば、上記は成り立つ。
つまり
\delta_{i,j}= \left\{
\begin{array}{ll}
1 & (i=j) \\
0 & (i \neq j)
\end{array}
\right. \
となる。
おそらく、この式は間違ってないだろう。
一方、$g[x]$について、フーリエ逆変換とフーリエ変換を用いると
\begin{align}
g[x]&=\frac{1}{N}\sum^{N-1}_{k=0} e^{2\pi i \frac{kx}{N}} \hat{g}[k] \\
&=\frac{1}{N}\sum^{N-1}_{k=0} e^{2\pi i \frac{kx}{N}}
\left(\sum^{N-1}_{x'=0}e^{-2\pi i \frac{kx'}{N}}g[x']\right) \\
&= \frac{1}{N}\sum^{N-1}_{k=0}\sum^{N-1}_{x'=0} g[x']
e^{2\pi i \frac{k}{N}(x-x')} \\
&= \sum^{N-1}_{x'=0} g[x'] \left(\frac{1}{N}\sum^{N-1}_{k=0} e^{2\pi i \frac{k}{N}(x-x')} \right)
\end{align}
この最後の式と最初に出てきた
g[x] = \sum^{N-1}_{x'=0} g[x']\delta_{x,x'}
見比べてもらえば
\delta_{x,x'} =\frac{1}{N}\sum^{N-1}_{k=0} e^{2\pi i \frac{k}{N}(x-x')}
であることが、わかってもらえると思う。
これで、クロネッカーのデルタのフーリエ表示ができた。
本題 パワースペクトルのフーリエ逆変換
さて、クロネッカーのデルタのフーリエ表示ができたところで本題に戻りたい。
本題は、離散化されたデータ、フーリエ変換でも自己相関関数とパワースペクトルの関係は成り立つのか？という話だった。
これを確かめるためパワースペクトルのフーリエ逆変換を計算してみる。
\begin{align}
\hat{P}[n]&=\frac{1}{N}\sum^{N-1}_{k=0}e^{2\pi i \frac{kn}{N}}P(k) \\
&=\frac{1}{N}\sum^{N-1}_{k=0}e^{2\pi i \frac{kn}{N}} \hat{f}^*[k]\hat{f}[k] \\
&=\frac{1}{N}\sum^{N-1}_{k=0}\sum^{N-1}_{k'=0}e^{2\pi i \frac{kn}{N}} \hat{f}^*[k]\hat{f}[k] \delta_{k',k} \\
&=\frac{1}{N}\sum^{N-1}_{k=0}\sum^{N-1}_{k'=0}e^{2\pi i \frac{kn}{N}} \hat{f}^*[k]\hat{f}[k'] \frac{1}{N}\sum^{N-1}_{n'=0} e^{2\pi i \frac{n'}{N}(k'-k)}\\
&=\frac{1}{N^2} \sum^{N-1}_{n'=0} \sum^{N-1}_{k=0} \sum^{N-1}_{k'=0} \hat{f}^*[k]\hat{f}[k'] e^{2\pi i \frac{1}{N}\left\{ k(n-n')+n'k' \right\}}\\
&=\sum^{N-1}_{n'=0} \left(\frac{1}{N} \sum^{N-1}_{k=0} \hat{f}^*[k] e^{2\pi i \frac{k}{N}(n-n') }\right)\left(\frac{1}{N} \sum^{N-1}_{k'=0} \hat{f}[k'] e^{2\pi i \frac{k'n'}{N}} \right) \\
&=\sum^{N-1}_{n'=0} \left(\frac{1}{N} \sum^{N-1}_{k=0} \hat{f}[k] e^{2\pi i \frac{k}{N}(n'-n) }\right)^*\left(\frac{1}{N} \sum^{N-1}_{k'=0} \hat{f}[k'] e^{2\pi i \frac{k'n'}{N}} \right) \\
&= \sum^{N-1}_{n'=0} f^*[n'-n] f[n']
\end{align}
となる。これをわかりやすいように書き直すと
\frac{1}{N}\sum^{N-1}_{k=0}e^{2\pi i \frac{kt}{N}}P(k) = \sum^{N-1}_{t'=0} f^*[t'-t] f[t'] \equiv \xi[t]
となる。$\xi[t]$は離散化での自己相関関数。
まさに見事に離散化でも自己相関とパワースペクトルの関係は成り立つではないか。
と思ったら、それは要注意である。
式の解釈
\xi[t] \equiv \sum^{N-1}_{t'=0} f^*[t'-t] f[t']
よく式を見てみよう。特に$f^*[t'-t]$のところだ。
そもそも我々がやってたことはなんだろうか？
それは離散的で有限範囲のデータのパワースペクトルと自己相関の関係を調べていた。そう、データは有限範囲である。
上の式ではデータ$f[t]$は$0\leq t \leq N-1$までしか定義されていない。
それなのに$f[t'-t]\ (t'=0,1,\cdots,N-1)$を考えると定義されていない配列を参照することになる。さて、どういうことだろうか。
そもそも、離散フーリエ変換とはどのようなものであっただろうか？
Qiita上でもいくつか解説記事があげられている。
詳しくはそれらを参照してほしいが、簡単にいってしまうと、
離散フーリエ変換は関数を周期関数として扱っている。
そのため、$f[N+a]=f[a]$となるはずだ。
つまり、上の式は周期関数の自己相関と同じになっているというだけの話である。
まとめ
わかったことは次の一点
離散フーリエ変換から計算する自己相関関数は、もとの関数を周期関数とみなして計算する自己相関関数と一致する
ということだ。
一般的な自己相関関数は周期関数とみなさない場合の方が多いかもしれない。
まあ、自己相関関数が自分の意図している定義と違うかもしれないということだ。
あとがき
乱雑に書いたので、随時修正していきたい。
内容間違ってたら教えてください。
2022/10/01 ウィーナー-ヒンチンの定理で一部間違いがあったため修正