信号解析 第３回講義録
日時：２０１４年４月２８日
講義内容：スペクトルとピリオドグラム
担当者：情報知能工学科 小島史男
はじめに
1
前回の講義では時系列信号の特徴が平均値関数、自己共分散関数、自己相関関数とよばれる３つの統計
量で特徴づけられること、またその推定値は有限のサンプルデータから推定でき、とくに弱定常過程の
時系列信号では標本データをとればとるほど良い近似になっていることを説明しました。今週の講義で
は、時系列信号の統計量の特徴を周波数領域、すなわちスペクトルで記述することを学習します。今回
学ぶことは２年後期のスペクトル解析の復習になります。教科書では第３章の内容です。
前回の復習
2
宿題の宿題 (1) での統計ソフト R を用いて英国のガス使用量の時系列データ UKgas をここにまとめて
みましょう。以下のようにデータを呼び出し、結果を表示しますと次のようになります。一番左のコラ
ムは西暦で、次の４つのコラムで 1960 年から 1986 年までを４つの季節に分割しガスの使用量が記録
されています。またこのデータを並べるとデータ長 N = 108 の時系列データが得られます。コマンド
ts.plot(UKgas) を使うとこの時系列データを図示 (図 1) することができます。
英国のガス消費量
> y <- UKgas
> y
Qtr1
Qtr2
Qtr3
Qtr4
1960 160.1 129.7
84.8 120.1
1961 160.1 124.9
84.8 116.9
1962 169.7 140.9
89.7 123.3
1963 187.3 144.1
92.9 120.1
1964 176.1 147.3
89.7 123.3
------------------------------------------------------------1985 1087.0 534.7 281.8 787.6
1986 1163.9 613.1 347.4 782.8
>
この時系列データの標本自己相関関数を求めてみましょう。
UKgas
> acf(y, type="correlation",plot=TRUE)
> dev.copy2eps(file="UKgasCorrelation.eps")
この解析結果から見えてくる信号の性質について考えましょう。自己相関関数がある一定の周期をもっ
て関係が持続していることがわかります。つまりガスの使用量は季節の移り変わりと対応していること
がわかります。一方前回の宿題の (2) の結果はどうだったでしょうか。図 3 はその標本（サンプルごと
異なります）の自己相関関数をプロットしたものです。UKgas の場合は周期性を持っていますが、この
場合は最初だけ相関が強いことがわかります。つまり時系列が相関を持たないことがわかります。時系
列データの性質をさらに定量的に把握する方法について、それでは本日の学習に入ります。
1200
1000
800
600
200
400
UKgas
1960
1965
1970
1975
1980
1985
Time
Figure 1: UKgas の時系列信号
0.4
0.2
0.0
−0.2
ACF
0.6
0.8
1.0
Series y
0
1
2
3
4
Lag
Figure 2: UKgas の時系列信号
5
0.4
−0.2
0.0
0.2
ACF
0.6
0.8
1.0
Series z
0
5
10
15
20
Lag
Figure 3: 白色雑音の自己相関関数
時系列データと定常性
3
宿題 (2) の白色雑音の標本平均は 0 分散は 1 のあたりで一定値をとり、また自己相関関数からは各時間
で独立（無相関）となることから、時間差で一定値をとることはがわかります。つまり白色雑音は定常
性をもっているとみなしてよいことになります。
しかし宿題 (1) の UKgas はどうみても定常性の条件を満たしているようにはみえません。実データ
では各種の前処理を施すことによって定常性を検証します。UKGass のような長期にわたるデータの変
化をみるためにしばしば対数をとることがあります。図 4 はもとの時系列データの対数をとった時系列
をプロットしたものです。このデータからは季節変動が若干ですがならされて平均化されていることが
わかります。しかし依然として標本平均は一定値とならず上昇傾向となります。このような傾向をトレ
ンドと呼びます。ただこの上昇は一定値のようにも見て取れます。このような場合数列の階差をとるこ
とで上昇傾向を取り除くことができます。つまり UKgas の時系列データを y[n] として
z[n] = log(y[n]) − log(y[n − 1])
としてこれをプロットし、３つの時系列データを比較してみましょう（図 5）。対数をとり階差をとるこ
とで、この時系列データはほぼ定常性を満足していることがわかります。
以下はを R で実行するときの命令文の例です。各自確かめてみてください。
コマンドライン
>
>
>
>
>
>
>
>
par(mfrow=c(3,1))
y <- UKgas
lgy <- log(UKgas)
z <- diff(lgy, lag = 1)
ts.plot(y)
ts.plot(lgy)
ts.plot(z)
dev.copy2eps(file="plot3.eps")
1200
1000
800
200
400
600
UKgas
1960
1965
1970
1975
1980
1985
Time
y
200
600
1000
Figure 4: UKGass の対数データ
1960
1965
1970
1975
1980
1985
1975
1980
1985
1975
1980
1985
lgy
4.5
5.5
6.5
Time
1960
1965
1970
−0.5
z
0.5 1.0
Time
1960
1965
1970
Time
Figure 5: UKgas の元データ、対数データおよび階差データの比較
4
自己共分散関数の離散時間フーリエ変換 (DTFT)
時系列信号の自己共分散関数 Ck （もしくは自己相関関数 Rk ）は時系列信号の時間方向の特徴を統計的
に記述する有効な指標であることは前回学習しました。たとえば不規則性が強いと時間方向への相関が
急激に減少しますが、周期信号の場合は周期的に減衰していくことが確かめられます。これらの特徴を
周波数領域で捉えることは重要です。自己共分散関数 Ck の離散時間フーリエ変換 (DTFT) の１周期分
∞
∑
p(f ) = DTFT{CK } =
(
Ck exp(−2πjkf ) = hCk , exp(2πjkf )i
k=−∞
1
1
− ≤f ≤
2
2
)
を Ck のパワースペクトル密度関数 (power spectral density function) あるいはスペクトル (Spectrum)
と呼びます。なおここでは後の議論の便宜上 DTFT で用いた Ω = 2πf T の代わりに周波数 f の関数と
しています。DTFT では時間領域では離散でも、周波数領域では f に関する連続関数となることは２年
次に学習しています。さらに自己共分散関数は時間に関して偶対象、すなわち偶関数 (Ck = C−k ) とな
りますので、自己共分散関数が実数値関数であればスペクトルも実部のみの偶関数となります。具体的
には
p(f ) = C0 + 2
∞
∑
Ck cos(2πkf )
k=1
となります。この逆変換からスペクトルが与えられれば、自己相関関数は次式によって求めることがで
きます。
∫
Ck =
1
2
− 12
∫
p(f ) exp(2πjkf )df =
1
2
− 12
p(f ) cos(2πkf )df
一例として、時系列が理想的な雑音であるガウス型白色雑音（ホワイトノイズ）のスペクトルを計算し
てみましょう。ホワイトノイズとは自己共分散関数が C0 = σ 2 , Ck = 0 (k 6= 0) で与えられる時系列信号
です。公式を適用すれば、スペクトルは
p(f ) =
∞
∑
Ck cos(2πkf ) = C0 = σ 2
k=−∞
となります。つまり周波数全帯域で一定値 σ 2 をとりますので、その周波数上での積分量は無限大になっ
てしまいます。実際はこのような雑音をつくりだすことは不可能ですが、シミュレーションにおいては、
一様乱数から擬似的にこのような雑音を生成することができます。
5
標本共分散関数の離散フーリエ変換 (DFT)
ˆk を導入しました。これを利用すれ
前回の講義で時系列信号の統計量を推定する標本自己共分散関数 C
N −1
ば、同様にスペクトルの推定値をもとめることができます。時系列 {y[n]}n=0
があたえられたとき、標
本自己相関関数の離散フーリエ変換 (DFT)
N
−1
∑
pˆ[l] =
Cˆk exp(−2πjkfl ) = Cˆ0 + 2
k=−N +1
N
−1
∑
Cˆk cos(2πkfl )
k=1
をピリオドグラム (periodgram) と呼びます。ここで周波数は離散値 fl = l/N (l = 0, 1, 2, · · · , [N/2]) を
とるものとします。つまり、標本自己共分散関数が周期的に無限領域に拡張されています。またピリオ
ドグラムを連続的に拡張した
N
−1
∑
pˆ(f ) =
(
Cˆk exp(−2πjkf )
k=−N +1
1
1
− ≤f ≤
2
2
)
を標本スペクトルともいいます。また逆に標本スペクトルから標本自己共分散関数は
Cˆk =
∫
1
2
− 12
pˆ(f ) exp(2πjkf )df
(k = 0, 1, 2, · · · , N − 1)
で求めることができます。このようにして、スペクトルの推定値は標本自己共分散関数から容易にもと
めることができます。以下では計算のテクニックについて学習していきましょう。
6
ピリオドグラムの良い近似法
ˆk は一般的に時間差
標本自己共分散関数の計算式をみればわかるように、一定のデータ長のもとでは、C
k が大きくなり、データ長 N に近づくと平均をとるサンプル数が少なくなり、あまり意味をなさなくな
ります。実際にピリオドグラムを計算するには対象となるデータ長 N の時系列信号を適当な長さ L で再
分割し、これによって計算された L 個のピリオドグラムの平均をとる手法があります。すなわち、時系
列データを長さ N/L 個のデータに再分割し、ピリオドグラム
(i)
pˆ [l] = Cˆ0 + 2
(i)
L−1
∑
(i)
Cˆk cos(2πkfl )
(
k=1
l
fl =
2L
(
[ ] )
L
l = 0, 1, · · · ,
,
2
N
, i = 1, 2, · · · ,
L
)
を再分割した区間ごとに求めていき、この L 個のピリオドグラムの平均
pˆ[l] =
N/L
L ∑ (i)
pˆ [l]
N i=1
によって平均化すると良好な近似ができます。ただこの方法だけでは、区間 L の大きさによっては十分
な近似ができない場合が多いので、スペクトルウインドウと呼ばれる畳み込み演算（フィルタ）をかま
せて平滑化する方法が一般的に用いられます。具体的には
m
∑
p˜[l] =
Wi pˆ[l − i]
i=−m
によってスペクトルの推定値を求めます。Wi はハニング窓、ハミング窓などが用いられます。窓関数に
ついてはスペクトル解析で説明しました。このようにピリオドグラムの実際の計算に関しては、平均化
と平滑化により近似精度の向上を図ることができます。
7
ピリオドグラムの計算の高速化
では本日の講義の最終部分に移ります。ピリオドグラムの計算は基本的には自己共分散関数の離散フーリ
エ変換です。したがって高速フーリエ変換 (FFT) を使って計算の高速化が可能です。たとえばデータ長 N
の標本自己共分散関数のピリオドグラムを FFT で計算すれば計算量は o(N 2 ) から o(N/2 log N ) まで節
ˆk (k = 0, 1, 2, · · · , N − 1)
減できることはスペクトル解析で学習しました。ところで標本自己共分散関数 C
2
の計算には o(N /2) の演算が必要ですので、このままではあまり高速化できていません。時系列信号の
−1
データから直接 FFT の適用を考えたほうが高速化を達成することができます。いま時系列信号 {y[n]}N
n=0
の離散フーリエ変換は
)
(
N
−1
∑
2πjnl
X[l] =
y[n] exp −
N
n=0
で計算できます。このとき以下の関係が導出できます。
1
1
|X[l]|2 =
X[l]X[l]
N
N
{ −1
)} {N∑
)}
(
(
−1
1 N∑
2πjnl
2πjml
=
y[n] exp −
y[m] exp
N n=0
N
N
m=0
=
)
(
−1 N
−1
∑
1 N∑
2πj(n − m)l
y[n]y[m] exp −
N n=0 m=0
N
さらに m = n + |k| とおくと
1
|X[l]|2 =
N
=
N
−1
∑
n−|k|
∑
(
1
2πjkl
y[n]y[n + |k|] exp −
N
N
n=1
k=−N +1
N
+1
∑
(
2πjkl
Cˆk exp −
N
k=−N +1
= p[l]
)
)
−1
の関係式が得られます。このことを利用すれば、ピリオドグラム pˆ[l] は時系列信号時系列信号 {y[n]}N
n=0
の離散フーリエ変換 X[l] を FFT で高速計算するだけで直接計算することができます。スペクトル解析
の講義で説明しましたが、FFT の計算にあたってはデータの再分割は２のべき乗で行うべきです。ゼロ
拡張 (zero-padding) を行うと近似精度が極端に落ちる場合がありますので注意してください。
R を使ったおさらい
8
前回と同じ時系列信号（教科書 p41 の式 (3.18) の信号生成モデル）を使ってピリオドグラムを求めてみ
ましょう。時系列信号が
(
)
(
)
2πn
2πn
y[n] = cos
+ cos
+ w[n]
10
4
で与えられ、データ長は N = 400、また w[n] は平均値 0、分散が正定値の正規型白色雑音とします。R
でのコマンドラインは次のようになります。
コマンドライン
>
>
>
>
source("model20120423.R")
z <- model20120423(400, 0.01)
spec.pgram(z, type="h")
dev.copy2eps(file="periodgram.eps")
以下はこれを実行した結果の例です。図 6 は白色雑音の分散が 0.01 を設定していますが図??は白色雑音
1e−02
1e−04
spectrum
1e+00
1e+02
Series: z
Raw Periodogram
0.0
0.1
0.2
0.3
0.4
0.5
frequency
bandwidth = 0.000713
Figure 6: 時系列信号のピリオドグラム
の分散が 1 の場合の結果です。ノイズレベルが異なっても（もちろん雑音系列が異なっても）ピリオド
グラムをみるとノイズに埋もれた信号から２つの周波数成分の信号が確認できます。また白色雑音はそ
のまわりで周波数で全域にわたってほぼ一定になるのが見て取れます。
下記は宿題ではありませんが、先週の結果をやり直しながら、本日の復習をしてください。
練習： 先週の宿題 (1)(2) のピリオドグラムを求めてみましょう。またそれらの結果と自己相関関数を
比較してみてください。
© Copyright 2024 ExpyDoc
