最小２乗法でフーリエ変換？：Lomb-Scargleスペクトルとコサイナー解析の一歩手前

　ピリオドグラム (Periodogram) の歴史について調べてみると、そもそもの目的は観測データに潜む周期性を見つけることのようです。これは、Arthur Schusterの論文 (1898)の導入に書いてありました。名前からして「周期 (period)」の「記録 (-gram)」なので、これは当然のことだと言えます。

　この点を踏まえると、確率過程のパワースペクトルをピリオドグラムと呼ぶことには少し違和感があります。なぜなら、確率過程のパワースペクトルでは、ピーク位置にのみ注目するのではなく、構造全体を読み解く必要があるからです。

　そんな私の感想はどうでもよいので、今回は、正弦波の最小二乗フィッティングを用いたピリオドグラム (パワースペクトル) の推定 について考えてみます。この方法は、Lomb-Scargleスペクトル、コサイナー解析と密接に関係していますので，これらの方法について理解するために、今回の内容は役に立つと思います。

最小２乗法でフーリエ変換をまねてみる
フーリエ変換もやってみる
フーリエ変換でも最小２乗法でも結果は一致
Rで検証
【補足】計算の説明

最小２乗法でフーリエ変換をまねてみる

　フーリエ変換を用いると、時系列データを異なる周波数の正弦波に分解できます。フーリエ変換の計算は、sin と cos の直交性を利用した内積に基づいています。

　ここではいったんフーリエ変換のことは忘れて、時系列データに対し sin と cos を最小二乗フィッティングすることで、「フーリエ変換っぽいこと」を試みます。

　時系列

$\displaystyle{ \{x[0], x[1], x[2], \cdots, x[N-1]\}}$

に対し、

$\displaystyle{ x[t] = a_f\, \cos (2 \pi f t) + b_f \, \sin (2 \pi f t) }$

を最小２乗法でフィットします．

　下の図はこの方法で振幅の周波数依存性を調べた結果です。実はこの最小２乗フィットの結果は、フーリエ変換を使った計算結果と一致します (周波数を同じ値で計算すれば完全一致)。

サンプル時系列 (上)にcosとsinをフィットして計算した振幅スペクトル (下)．A(f)は周波数fの成分の振幅 $\displaystyle{ A(f) = \sqrt{a _ f^ 2 + b _ f^ 2} }$ ．

\displaystyle{
A(f) = \sqrt{a _ f^ 2 + b _ f^ 2}
} — サンプル時系列 (上)にcosとsinをフィットして計算した振幅スペクトル (下)．A(f)は周波数fの成分の振幅 $\displaystyle{ A(f) = \sqrt{a _ f^ 2 + b _ f^ 2} }$ ．

　残差の２乗和を

$\displaystyle{ I_f (a_f, b_f) = \sum_{t=0}^{N-1} \left\{ x[t] - \left(a_f\, \cos (2 \pi f t) + b_f \, \sin (2 \pi f t) \right) \right\}^2 }$

として、これを最小化する $a _ f$ と $b _ f$ を求めます。この関数は、 $a _ f$ と $b _ f$ の２次関数なので、必ず最小値が存在します。ですので、

$\displaystyle{ \begin{align} \frac{\partial I_f (a_f, b_f)}{\partial a_f} &= 0 \\ \frac{\partial I_f (a_f, b_f)}{\partial b_f} &= 0 \\ \end{align} }$

となる、 $a _ f$ と $b _ f$ を求めれば良いということになります。

　ちょっと計算を省略して、行列の形で書けば、

になります。式を繰り返し書くのが面倒なので、以下の省略した表現を定義します。

$\displaystyle{ \begin{align} {\rm c}^2_{\rm sum} &= \sum_{t=0}^{N-1} \cos^2 ( 2 \pi f t )\\ {\rm s}^2_{\rm sum} &= \sum_{t=0}^{N-1} \sin^2 ( 2 \pi f t )\\ {\rm c}{\rm s}_{\rm sum} &= \sum_{t=0}^{N-1} \cos ( 2 \pi f t ) \sin ( 2 \pi f t ) \end{align} }$

この表記を使えば、

$\displaystyle{ \begin{align} \left[\begin{array}{cc} \displaystyle \sum_{t=0}^{N-1} \cos^2 ( 2 \pi f t ) & \displaystyle \sum_{t=0}^{N-1} \cos ( 2 \pi f t ) \sin ( 2 \pi f t ) \\ \displaystyle \sum_{t=0}^{N-1} \cos ( 2 \pi f t ) \sin ( 2 \pi f t ) & \displaystyle \sum_{t=0}^{N-1} \sin^2 ( 2 \pi f t ) \end{array}\right]\left[\begin{array}{c} a_f \\ b_f \end{array}\right] &= \left[\begin{array}{c} \displaystyle \sum_{t=0}^{N-1} x \left[t\right] \cos ( 2 \pi f t ) \\ \displaystyle \sum_{t=0}^{N-1} x \left[t\right] \sin ( 2 \pi f t ) \end{array}\right] \\ \left[\begin{array}{cc} \displaystyle {\rm c}^2_{\rm sum}& \displaystyle {\rm c}{\rm s}_{\rm sum} \\ \displaystyle {\rm c}{\rm s}_{\rm sum}& \displaystyle {\rm s}^2_{\rm sum} \end{array}\right]\left[\begin{array}{c} a_f \\ b_f \end{array}\right] &= \left[\begin{array}{c} \displaystyle \sum_{t=0}^{N-1} x \left[t\right] \cos ( 2 \pi f t ) \\ \displaystyle \sum_{t=0}^{N-1} x \left[t\right] \sin ( 2 \pi f t ) \end{array}\right] \\ \end{align} }$

上の式の左辺の行列の逆行列を両辺にかければ、 $a _ f$ と $b _ f$ が求まります。

$\displaystyle{ \begin{align} \left[\begin{array}{c} a_f \\ b_f \end{array}\right] = \frac{1}{{\rm c}^2_{\rm sum} \cdot {\rm s}^2_{\rm sum} - \left( {\rm c}{\rm s}_{\rm sum} \right)^2} \left[\begin{array}{cc} \displaystyle {\rm s}^2_{\rm sum} & \displaystyle - {\rm c}{\rm s}_{\rm sum} \\ \displaystyle - {\rm c}{\rm s}_{\rm sum} & \displaystyle {\rm c}^2_{\rm sum} \end{array}\right] & \left[\begin{array}{c} \displaystyle \sum_{t=0}^{N-1} x \left[t\right] \cos ( 2 \pi f t ) \\ \displaystyle \sum_{t=0}^{N-1} x \left[t\right] \sin ( 2 \pi f t ) \end{array}\right] \end{align} }$

この計算は、２次の行列の逆行列の公式を使っただけです。

　上の式は、めんどくさい形に見えるかもしれません。しかし、周波数がある条件を満たすときは、もっと簡単になります。

　周波数 $f$ を、離散フーリエ変換と同じく、 $\displaystyle{ k }$ を整数として、

$\displaystyle{ f_k=\frac{k}{N} }$

に設定すれば、

$\displaystyle{ \begin{align} {\rm c}^2_{\rm sum} &= \displaystyle \sum_{t=0}^{N-1} \cos^2 ( 2 \pi f_k t ) = \frac{N}{2} \\ {\rm s}^2_{\rm sum} &= \displaystyle \sum_{t=0}^{N-1} \sin^2 ( 2 \pi f_k t ) = \frac{N}{2} \\ {\rm c}{\rm s}_{\rm sum} &= \displaystyle \sum_{t=0}^{N-1} \cos ( 2 \pi f_k t ) \sin ( 2 \pi f_k t ) = 0 \\ \end{align} }$

となります (計算は【補足】を参照)。

　したがって、 $a _ {f _ k}$ と $b _ {f _ k}$ を、それぞれ、 $a _ k$ と $b _ k$ と書くことにして、

$\displaystyle{ \begin{align} \left[\begin{array}{c} a_k \\ b_k \end{array}\right] &= \frac{1}{{\rm c}^2_{\rm sum} \cdot {\rm s}^2_{\rm sum} - \left( {\rm c}{\rm s}_{\rm sum} \right)^2} \left[\begin{array}{cc} \displaystyle {\rm s}^2_{\rm sum} & \displaystyle - {\rm c}{\rm s}_{\rm sum} \\ \displaystyle - {\rm c}{\rm s}_{\rm sum} & \displaystyle {\rm c}^2_{\rm sum} \end{array}\right] \left[\begin{array}{c} \displaystyle \sum_{t=0}^{N-1} x \left[t\right] \cos ( 2 \pi f t ) \\ \displaystyle \sum_{t=0}^{N-1} x \left[t\right] \sin ( 2 \pi f t ) \end{array}\right]\\ &= \frac{2}{N}\left[\begin{array}{cc} \displaystyle 1 & 0 \\ \displaystyle 0 & 1 \end{array}\right] \left[\begin{array}{c} \displaystyle \sum_{t=0}^{N-1} x \left[t\right] \cos \left( 2 \pi \frac{k}{N} t \right) \\ \displaystyle \sum_{t=0}^{N-1} x \left[t\right] \sin \left( 2 \pi \frac{k}{N} t \right) \end{array}\right]\\ &= \left[\begin{array}{c} \displaystyle \frac{2}{N} \sum_{t=0}^{N-1} x \left[t\right] \cos \left( 2 \pi \frac{k}{N} t \right) \\ \displaystyle \frac{2}{N} \sum_{t=0}^{N-1} x \left[t\right] \sin \left( 2 \pi \frac{k}{N} t \right) \end{array}\right] \\ \end{align} }$

となります。つまり、

$\displaystyle{ \left\{\begin{array}{l} a_k = \displaystyle \frac{2}{N} \sum_{t=0}^{N-1} x \left[t\right] \cos \left( 2 \pi \frac{k}{N} t \right) \\ b_k = \displaystyle \frac{2}{N} \sum_{t=0}^{N-1} x \left[t\right] \sin \left( 2 \pi \frac{k}{N} t \right) \end{array}\right. }$

です。

　どこかで見たことのある形です。

$\displaystyle{ k }$ を整数として、周波数を $\displaystyle{ f _ k=\frac{k}{N} }$ で選んで、フィットした結果は、以下のようになります。

$\displaystyle{ x[t] = a_k \cos \left( 2 \pi \frac{k}{N} t \right)+ b_k \sin \left( 2 \pi \frac{k}{N} t \right) }$

ここで、 $a _ k$ と $b _ k$ は、

です。この式に登場する、和 $\sum$ の部分を

$\displaystyle{ \left\{\begin{array}{l} X_k^{(\cos)} = \displaystyle \sum_{t=0}^{N-1} x \left[t\right] \cos \left( 2 \pi \frac{k}{N} t \right) \\ X_k^{(\sin)} = \displaystyle \sum_{t=0}^{N-1} x \left[t\right] \sin \left( 2 \pi \frac{k}{N} t \right) \end{array}\right. }$

と表すことにします。そうすると、

$\displaystyle{ x[t] = \frac{2}{N} \left\{ X_k^{(\cos)} \cos\left( 2 \pi \frac{k}{N} t \right) + X_k^{(\sin)} \sin \left( 2 \pi \frac{k}{N} t \right) \right\} }$

となります。

フーリエ変換もやってみる

　上の最小２乗法の結果とフーリエ変換の結果を比較してみます。

　時系列

$\displaystyle{ \{x[0], x[1], x[2], \cdots, x[N-1]\}}$

のフーリエ変換は、 $i$ を虚数単位として、以下のようになりkます：

$\displaystyle{ \begin{align} X_k &= \sum_{t=0}^{N-1} x \left[t\right] \exp \left(- i 2 \pi \frac{k}{N} t \right) \\ &= \sum_{t=0}^{N-1} x \left[t\right] \left\{ \cos \left( - 2 \pi \frac{k}{N} t \right) + i \sin \left( - 2 \pi \frac{k}{N} t \right) \right\} \\ &= \sum_{t=0}^{N-1} x \left[t\right] \cos \left( 2 \pi \frac{k}{N} t \right) - i \sum_{t=0}^{N-1} x \left[t\right] \sin \left( 2 \pi \frac{k}{N} t \right) \end{align} }$

ここで、周波数は、

$\displaystyle{ f_k=\frac{k}{N} }$

なので、

$\displaystyle{ 0 \le f_k \le \frac{N-1}{N} < 1 }$

です。ただ、みなさんご存じのように、適切な周波数領域は、

$\displaystyle{ -\frac{1}{2} \le f_k \le \frac{1}{2} }$

ですので、注意してください (過去の記事参照)。

　先ほど定義した、

を使えば、フーリエ変換の結果は、

$\displaystyle{ \begin{align} X_k &= \sum_{t=0}^{N-1} x \left[t\right] \cos \left( 2 \pi \frac{k}{N} t \right) - i \sum_{t=0}^{N-1} x \left[t\right] \sin \left( 2 \pi \frac{k}{N} t \right) \\ &= X_k^{(\cos)} - i X_k^{(\sin)} \end{align} }$

と書けます。

　この $X_k$ が表す時間領域の周波数成分 $x_k [t$ ] は以下のようになります。

$\displaystyle{ x[t] = \frac{1}{N} \left\{ X_k^{(\cos)} \cos\left( 2 \pi \frac{k}{N} t \right) + X_k^{(\sin)} \sin \left( 2 \pi \frac{k}{N} t \right) \right\} }$

　この結果は、先ほど計算した最小２乗フィットの結果

$\displaystyle{ x[t] = \frac{2}{N} \left\{ X_k^{(\cos)} \cos\left( 2 \pi \frac{k}{N} t \right) + X_k^{(\sin)} \sin \left( 2 \pi \frac{k}{N} t \right) \right\} }$

と２倍だけ違います。

　フーリエ変換では、

$\displaystyle{ f_k = \frac{k}{N} }$

と、

$\displaystyle{ - f_k = -\frac{k}{N} }$

は、同じ周波数成分なので、両方の成分を足せば、フーリエ変換でも、

$\displaystyle{ x[t] = \frac{2}{N} \left\{ X_k^{(\cos)} \cos\left( 2 \pi \frac{k}{N} t \right) + X_k^{(\sin)} \sin \left( 2 \pi \frac{k}{N} t \right) \right\} }$

となり、最小２乗フィットの結果と一致します。

フーリエ変換でも最小２乗法でも結果は一致

　今回のポイントは、最小２乗法で正弦波をフィットしても、フーリエ変換してその結果から特定の周波数成分（正弦波の成分）を抽出しても結果は同じということです．

　この事実を使って、不等間隔にサンプリングされたデータのペリオドグラム（パワースペクトル）を推定できないか、と考えたのがLomb-Scargleスペクトルです。また、24時間周期のように、事前に特定の周期だけを仮定して分析するのがコサイナー解析です。

Rで検証

　最小２乗フィットでもフーリエ変換でも結果が一致することを確かめるRスクリプトの例は以下です。

# 振動成分をもつ2次自己回帰過程の時系列を生成する
#set.seed(1)  # 乱数の指定
n <- 256  # 時系列の長さ
a1 <- 0.5  # AR(2)の係数
a2 <- -0.75  # 
sig <- 0.5  # ノイズの標準偏差

# ホワイトノイズ
eps <- rnorm(n,mean=0,sd=sig)

# 時系列データの初期値を設定
x <- numeric(n)
x[1] <- eps[1]
x[2] <- a1*y[1]+eps[2]

# AR(2)プロセスの生成
for (t in 3:n) {
  x[t] <- a1 * x[t-1] + a2 * x[t-2] + eps[t]
}
#
time <- 1:n-1
#
FT <- fft(x)
A.FT <- abs(FT)/n
f.FT <- time/n
freq <- f.FT[f.FT > 0 & f.FT < 0.5]
n.freq <- length(freq)
Amp <- c()
###
for(i in 1:n.freq){
 f <- freq[i]
 # 最小二乗フィット
 xx <- cos(2*pi*f*time)
 zz <- sin(2*pi*f*time)
 #
 Y <- matrix(c(sum(x),sum(x*xx),sum(x*zz)),ncol=1) 
 X <- matrix(c(n,sx<-sum(xx),sz <- sum(zz),sx,sum(xx^2),sxz <- sum(xx*zz),sz,sxz,sum(zz^2)),ncol=3)
 X.inv <- solve(X)
 A <- X.inv %*% Y
 #　モデルパラメタの推定結果
 mesor.est <- A[1,1]
 phase.est <- atan2(A[3,1],A[2,1])
 Amp[i] <- sqrt(A[2,1]^2+A[3,1]^2)
}
##
# グラフのプロット
par(mfrow=c(3,1),mar=c(5,5,3,2),cex.main=1.8)
plot(time, x,"l", xlab = "t", ylab = "x[t]", pch = 16, col=4,lwd=2,xaxs="i",main="Sample time series of AR(2)")
plot(freq,Amp,"l",xaxs="i",col=2,lwd=2,xlim=c(0,0.5),main="Amplitude spectrum |A(f)| of least-square sinusoidal fitting",xlab="Frequency",ylab="")
plot(f.FT,A.FT*2,"l",xaxs="i",col=3,lwd=2,xlim=c(0,0.5),main="Amplitude spectrum |A(f)| of Fourier transform",xlab="Frequency",ylab="")

【補足】計算の説明

　以下の計算について説明しておきます。

$\displaystyle{ \sum_{t=0}^{N-1} \cos^2 \left( \frac{2\pi k}{N} t \right) }$

まずは、高校で習う公式

$\displaystyle{ \cos^2 x = \frac{1 + \cos 2x}{2} }$

を

$\displaystyle{ \cos^2 \left( \frac{2\pi k}{N} t \right) }$

に適用すると、

$\displaystyle{ \begin{align} \sum_{t=0}^{N-1} \cos^2 \left( \frac{2\pi k}{N} t \right) &= \sum_{t=0}^{N-1} \frac{1 + \cos \left( \frac{4\pi k}{N} t \right)}{2} \\ &= \frac{1}{2} \sum_{t=0}^{N-1} 1 + \frac{1}{2} \sum_{t=0}^{N-1} \cos \left( \frac{4\pi k}{N} t \right) \end{align} }$