計測自動制御学会東北支部 第 294 回研究集会 (2015.5.29)
資料番号 294-5

巻き込み 2 変量正規分布に従う位相スペクトルをもつ
2 信号間の位相限定相関関数の統計的性質
Statistical Property of Phase-Only Correlation Functions
Between Two Signals With Phase Spectrum Following Wrapped
Bivariate Normal Distributions
鈴木亮 ∗ , 八巻俊輔 ∗ , 川又政征 ∗ , 吉澤誠 ∗Ryo Suzuki∗ , Shunsuke Yamaki∗ , Masayuki Kawamata∗ , Makoto Yoshizawa∗
∗∗

東北大学

Tohoku University

キーワード : 方向統計学 (directional statistics), 位相限定相関関数 (phase-only correlation functions), 巻
き込み 2 変量正規分布 (wrapped bivariate normal distributions)

連絡先 : 仙台市青葉区荒巻字青葉 6-6-05 東北大学大学院工学研究科 電気エネルギーシステム専攻 吉澤・杉田研究室
鈴木亮, Tel.:022-795-7130, E-mail: r-suzuki@yoshizawa.ecei.tohoku.ac.jp

1.

はじめにの位相スペクトル差は 0 である．位相スペクト

2 つの信号の類似度を評価する方法の一つとして，信号のもつ位相情報を用いる位相限定相関関数（POC 関数:Phase-only-Correlation Func-tions）を用いる方法がある．POC 関数は振幅スペクトルを 1 に正規化した信号（位相限定信号）に対して相関を計算することで求められる．
また POC 関数は，2 信号間の幾何学的な関係を求めることができる．例えば，信号の位置ずれを POC 関数のピークの出現する位置のずれに変換して表現することができる．これらの特徴から，POC 関数は指紋認証技術 1) や画像マッチング技術 2) ，周期性のある DNA 配列の探索
3) などに応用されてきた．

2 つの評価したい信号が等しいときは，それぞれの信号の持つ位相情報は等しく，2 つの信号ル差が 0 となる時の POC 関数はデルタ関数になる．また，2 つの信号の位相スペクトル差が 0ではない時，POC 関数はデルタ関数と異なる．
しかし，これまでの研究では，2 信号間の位相スペクトル差が 0 ではない時に，POC 関数がどれだけデルタ関数と異なるかについて数学的な根拠が与えられていなかった．
著者らのグループでは，これまでに 2 信号間の位相スペクトル差が 0 でない場合について，信号の位相スペクトル差を 1 変量確率変数と仮定して，POC 関数の統計的解析を行った 4) 5) ．文
献 4) では，位相スペクトル差が直線上で与えられる確率分布に従うと仮定して，POC 関数の期待値と分散を導出した．また文献 5) では，位相スペクトル差のもつ周期性を考慮し，位相ス–1–

ペクトル差が円周上で与えられる確率分布に従て以下の式で与えられる．

うと仮定し，POC 関数の期待値と分散を導出しN −1r(m) =た．
これらに対して本稿では，2 信号の位相スペクトルを 2 変量確率変数と仮定して，POC 関=数の統計的解析を行う．2 信号の位相スペクト

1 ∑ X(k)Y ∗ (k)W −mkN|X(k)||Y (k)| N1
Nk=0N−1∑ejαk WN−mk(3)k=0

ルが平面上で与えられる 2 変量正規分布に従う

ここで ∗ は複素共役を表し，αk = θk − ϕk はと仮定して POC 関数の期待値と分散を導出す

2 つの信号の位相スペクトル差を表す．つまりる．次に 2 信号の位相スペクトルが平面上で与位相スペクトル差 αk がわかると POC 関数を求えられる確率密度関数に従うとした場合の問題めることができる．

点を指摘し，巻き込み法を用いて 2 変量正規密度関数を球面上の確率密度関数に変換する．そ2.2の結果，2 信号の位相スペクトルが巻き込み 2変量正規分布に従うと仮定たときの POC 関数の期待値と分散は，平面上で与えられる 2 変量正規分布を仮定したときの結果と一致すること

POC 関数の期待値と分散の一般式

著者らのグループでは，文献 4) で，位相スペクトル差 αk を確率変数と仮定して，位相スペクトル差 αk が 0 でない場合の POC 関数の期待値と分散を導出した．POC 関数の期待値 E[r(m)]を示す．

および分散 Var[r(m)] の一般式は以下の式で表される．

2.

POC 関数の統計的解析

2.1

位相限定相関関数の定義

E[r(m)] = Aδ(m)1
Var[r(m)] =(1 − |A|2 )N

まず最初に位相限定相関関数 (POC 関数) を定義する．信号長が N の 2 つの複素信号を x(n)(4)(5)

ここで位相因子 ejαk の期待値をと y(n) とする. これらの信号 x(n) と y(n) の 1

A = E[ejαk ](6)次元離散フーリエ変換は以下の式で表される．

X(k) =Y (k) =N−1∑
n=0
N−1∑

とおいている．式 (4) と (5)，(6) より位相スペクx(n)WNnk = |X(k)|ejθk(1)トル差 αk の確率密度関数を与えることで，POC関数の統計的性質を表すことができる．

y(n)WNnk = |Y (k)|ejϕk(2)n=02.3

ここで k = 0, …, N − 1 は離散周波数インデックスであり，WN = exp(−j2π/N ) は離散フーリエ変換の回転因子を表す．また |X(k)| は信号x(n) の振幅スペクトルであり，θk は x(n) の位相スペクトルである．同様に |Y (k)| は信号 y(n)の振幅スペクトルであり，ϕk は y(n) の位相スペクトルである．POC 関数は 2 つの信号のクロ特性関数を用いた POC 関数の期待値と分散の導出

位相スペクトル差 αk の確率密度関数が与えられた時，式 (6) の計算にその特性関数を用いることができる．確率変数 α の確率密度関数 p(α)の特性関数 ψα (t) は以下で定義される．
∫ ∞jαtψα (t) = E[e ] =ejαt p(α)dα

スパワースペクトルの離散フーリエ逆変換とし–2–−∞(7)

ここで t = 1 とすることで，位相因子の期待値を求めることができる．つまり，上式は以下の

ψα (1) = E[ejα ] =1f (θ, ϕ) =ように表せる．

∫ ∞−∞(8)従って，位相スペクトル差 α が確率変数に従う場合は、上式を用いることで POC 関数の期待値と分散を導出することができる．

eΘΣ

(9)

2π|Σ|()
Θ= θ ϕ()σθ2
ρσθ σϕΣ=ρσθ σϕ
σϕ2ejα p(α)dα= A−1 Θt12

このとき，Σ は共分散行列を表し，σθ2 , σϕ2 はそれぞれ 2 信号の位相スペクトル (θk , ϕk ) の分散を表している．このとき，σθ , σϕ は σθ > 0 とσϕ > 0 の関係を満たす．また ρ は 2 変量 θ, ϕ の3.

2 変量正規分布に従う位相スペクトルをもつ 2 信号間の POC関数の統計的性質関係を表すもので相関係数と呼ばれる．相関係数 ρ は −1 <=ρ<= 1 を満たす．相関係数 ρ = 1の時は，θ = ϕ であり，ρ = −1 の時は，θ = −ϕとなる．相関係数が 0 < ρ < 1 の時は，θ と ϕ著者らのグループがこれまでに行ってきた POCがともに増加する関係を満たし，−1 < ρ < 0 の関数の統計的解析では，2 信号の位相スペクト時は，θ と ϕ がともに減少する関係を満たす．

ル θk と ϕk について次のように仮定していた．

次に位相スペクトル差 αk と 2 変量 (θk , ϕk ) の信号 x(n) の位相スペクトル θk を確定信号とし，関係について考える．位相スペクトル差 αk は 2

一方，信号 y(n) の位相スペクトル ϕk を確率信信号の位相スペクトル (θk , ϕk ) の差であるため，号であると仮定をおいていた．この時，ϕk は周

位相スペクトル差 αk と (θk , ϕk ) の間には以下の波数領域で確率的に変動する特性をもつと仮定

ような関係を満たす行列 T が与えられる．
()
θkαk = T= θk − ϕk(10)ϕk

し，POC 関数の統計的解析をしていた．実際の信号処理においては，2 信号 x(n) と y(n) の位相スペクトル θk と ϕk がともに確率信号である(場合も考えられる．本研究では，2 信号の位相T =スペクトル θk と ϕk に周波数領域で確率的に変動する雑音が重畳すると仮定して，POC 関数の統計的解析を行う．この時，2 信号の位相スペ位相スペクトル差 αk の確率密度関数 p(αk ) は以下の式で表される．

クトル (θk , ϕk ) は，2 変量の確率変数であると

2αk1e− 2σ2p(αk ) = √2πσ2σ = T ΣT t仮定している．

3.12 変量正規分布に従う場合= σθ2 − 2ρσθ σϕ + σϕ2

2 信号の位相スペクトル (θk , ϕk ) が 2 変量 (θ, ϕ)の確率変数であると仮定する．(θ, ϕ) が平均 0 =(0, 0)，共分散行列 Σ の 2 変量正規分布 N (0, Σ)に従う時，その確率密度関数は以下の式で定義される．

1 −1)(11)(12)つまり 2 信号の位相スペクトル (θk , ϕk ) が平均0 = (0, 0)，共分散行列 Σ の 2 変量正規分布N (0, Σ) で与えられる場合，位相差 αk は平均0，分散 σ 2 の 1 変量正規分布 N (0, σ 2 ) に従うことがわかる．この時，位相因子の期待値は式 (8)–3–

より

A = ψα (1) = E[ejα ] =∫ ∞ejα p(α)dα−∞2=e− σ2と導出できる．この時，式 (4) から (6) と (13)より
σ2E[r(m)] = e− 2 δ(m)12(1 − e−σ )Var[r(m)] =N(13)(14)Fig. 1 2 信号の位相スペクトルの相関係数が変化した時の POC 関数の期待値の挙動を得る．式 (12) を式 (14) と (15) に代入すると，
以下のように POC 関数の期待値と分散の一般式を導出できる．

E[r(m)] = e− 2 (σθ −2ρσθ σϕ +σϕ ) δ(m)
(15)
{}
221
Var[r(m)] =1 − e−(σθ −2ρσθ σϕ +σϕ ) (16)N
13.222Fig. 2 2 信号の位相スペクトルの相関係数が変化した時の POC 関数の分散の挙動

POC 関数の挙動

信号 x(n) の位相スペクトル θk の分散 σθ2 を 1とし，信号 y(n) の位相スペクトル ϕk の分散 σϕ24.

位相スペクトルが巻き込み 2 変を 1 として，相関係数 ρ の変化に対する POC 関量正規分布に従う場合の POC

数 r(m) の期待値 E[r(0)] および分散 Var[r(m)]関数の統計的性質のそれぞれの値を Fig. 1 と 2 に図示する．σθ2 =1 かつ σϕ2 = 1 である時，式 (15) と (16) に代入
すると以下の式を得る．

3 章で定義した 2 変量正規分布は確率変数 (θ, ϕ)が平面上で定義されているものであった．確率変数 (θ, ϕ) は 2 信号それぞれの位相スペクトE[r(m)] = e−(1−ρ) δ(m)(17){}
1Var[r(m)] =1 − e−2(1−ρ)(18)
N
式 (17) より m = 0 の時の POC 関数の期待値は，以下のようにできる．

E[r(0)] = e−(1−ρ)

ルを表すものであり，定義域 −π <= θ < π，−π <= ϕ < π を満たさなければいけない．周期的な性質をもつデータは，角度データとしてみなすことができ，方向統計学に基づいて解析することができる．位相スペクトル (θ, ϕ) は 2π[rad](19)の周期性をもつため，2 変量の角度データとして取り扱うことができ，球面上の点と対応付け

式 (18) と (19) より相関係数の値が −1 から 1ることができる．今回は，方向統計学の概念をに増加するに従い，それぞれ POC 関数の分散導入し確率変数 (θ, ϕ) を球面上の分布を用いてVar[r(m)] は減少し，期待値 E[r(0)] は増加する

解析する．

ことがわかる．

–4–

Fig. 5線形ヒストグラム

Fig. 6円周ヒストグラム

Fig. 3 2 変量位相スペクトル (θk , ϕk ) を平面上にプロットした例

Fig. 4 2 変量位相スペクトル (θk , ϕk ) を球面上にプロットした例

4.2

方向統計学の基礎

方向統計学とは，気象学における風向データや渡り鳥の移動方向のデータなど，方向や角度4.1

角度データの統計的処理データを統計的に解析する学問である．角度デー

3 章では，2 変量 (θk , ϕk ) は平面上で表される
1 点と対応づけて，統計的に POC 関数の解析をタは円周上の点と対応させ，円周上の分布を用いて解析することができる 6) ．例としてある地してきた．例えば，2 変量 (θk , ϕk ) = ( π4 , π4 ) の点で観測した風向きの方向を線形ヒストグラム

時，平面上では Fig. 3 のように表すことができと円周ヒストグラムで表した場合について，Fig.

る．位相スペクトル (θk , ϕk ) は，本来 Fig. 3 の5 と 6 にそれぞれ示す． Fig. 5 では，角度デー点線で示される範囲をとりうる．この時，Fig. 3

タを線形データのように扱い，図示した．Fig.の斜線で示される範囲で表される点は，(θk , ϕk )5 を見ると，北西と北東の 2 つの方向に風向がの定義域外での値をとっている．角度データの集中しているように見える．一方で，Fig. 6 のもつ周期性を考慮するために，角度データを球ように円周上にプロットすると，北の方向に風面上の 1 点に対応付けしてこの問題を解決する．

向が集中していることが確認できる．Fig. 5 と2 変量 (θk , ϕk ) = ( π4 , π4 ) の時，球面上では Fig.

6 の 2 つは，同じデータを用いているが線形ヒ4 のように表すことができる．

ストグラムと円周ヒストグラムの 2 つを用いてデータの表示方法を変えている．風向きの方向の例では，円周ヒストグラムで表示すると風向きの集中方向を理解しやすい．このように，角–5–関数を表している．

f (θ, ϕ) =Q=2πσθ σϕ∞∑∞∑

1 − ρ2 k1 =−∞ k2 =−∞(θ − 2πk1 )2 (ϕ − 2πk2 )2+σθ2
σϕ2−Fig. 71√Qe2(21)2ρ(θ − 2πk1 )(ϕ − 2πk2 )σθ σϕ1 変量正規分布の巻き込み分布

4.5
度データを取り扱う際には，線形データを統計的に解析する方法とは違う統計的手法がある．

位相因子の期待値の導出

位相因子の期待値の導出を行う．式 (4) と (5)，(6) により位相因子 ejαk の期待値 A がわかれば，POC 関数の期待値および分散を導出すること4.3

巻き込み分布ができる．式 (6) で与えられている α は 2 信号

球面上で与えられる確率密度関数は周期性を間の位相スペクトル差 αk を表しているため，2持つ必要がある．平面上で与えられた確率密度変量の位相スペクトル (θk , ϕk ) を用いて以下の関数を球面上に巻き込んだ分布を巻き込み分布ように表すことができる．という．2 変量確率密度関数の巻き込み分布は

A = E[ejαk ] = E[ej(θk −ϕk ) ]以下の式で与えられる．

f (θ, ϕ) =∞
∑∞
∑

(22)

オイラーの公式を使い，式 (22) の右辺を以下のg(θ + 2πk1 , ϕ + 2πk2 )(20) ように展開する．

k1 =−∞ k2 =−∞この時，g(θ, ϕ) と f (θ, ϕ) はそれぞれ，平面上E[ej(θk −ϕk ) ] = E[cos(θk − ϕk ) + j sin(θk − ϕk )]と球面上の確率密度関数を表す．例として，平= E[cos(θk − ϕk )] + jE[sin(θk − ϕk )]均 0，分散 3 の正規分布 N (0, 3) を円周上に巻き= E[cos θk cos ϕk ] + E[sin θk sin ϕk ]込んだ分布を Fig. 7 に示す．ここでは簡単のた+ jE[sin θk cos ϕk ] − jE[cos θk sin ϕk ](23)

め，円周上に巻き込む分布を 1 変量正規分布と定義している．この時，g(θ) と f (θ) は，それぞれ 1 変量正規分布と巻き込み 1 変量正規分布を表している．

4.6

POC 関数の期待値の導出

文献 7) により，2 変量 (θ, ϕ) が巻き込み 2 変

4.4量正規分布に従う場合に，2 変量の角度データ

巻き込み 2 変量正規分布

式 (9) で表される平均 0 = (0, 0)，共分散行列

Σ の 2 変量正規密度関数 N (0, Σ) を球面上に巻き込んだ分布を以下の式で表す．この時，f (θ, ϕ)は，平均方向 (0, 0) の巻き込み 2 変量正規密度を表す統計量として以下の式が与えられる．
}{
1 22E[sin θ sin ϕ] = exp − (σθ + σϕ ) sinh(ρσθ σϕ )2(24){}1E[cos θ cos ϕ] = exp − (σθ2 + σϕ2 ) cosh(ρσθ σϕ )2(25)E[sin θ cos ϕ] = E[cos θ sin ϕ] = 0–6–(26)式 (24) と (25)，(26) を式 (23) に代入すると，位参考文献

相因子の期待値は以下のように導出できる．
{}
1 2j(θk −ϕk )2E[e] = exp − (σθ − 2ρσθ σϕ + σϕ ) (27)21) H. Nakajima, K. Kobayashi, T. Aoki, andT. Higuchi,

Pattern collation apparatus

この時，POC 関数の期待値と分散は式 (23) を

based on spatial frequency characteristics

式 (4) と (5) に代入することにより，以下の式で(USP 5915034), US patent, May 1995.

導出できる．

2) C. D. Kuglin and D. C. Hines,

The Phase

Correlation Image Alignment Method,E[r(m)] = Aδ(m){}1
Proc. IEEE, Int. Conf. on Cybernetics and= exp − (σθ2 − 2ρσθ σϕ + σϕ2 ) δ(m)2Society, pp. 163-165, 1975.
(28)3) Andrzej K. Brodzik,Phase only filter1Var[r(m)] = (1 − |A|2 )N(ing for the masses(of DNA data): A new)
{ ()}
12
2=
1 − exp − σθ − 2ρσθ σϕ + σϕapproach to sequence alignment, IEEEN
Trans. Signal Processing, VOL. 54, NO. 6,(29)pp. 2456-2466, June 2014.
ここで式 (28) と (29) は 2 変量 (θ , ϕ ) に 2 変kk

量正規分布を仮定して導出された式 (15) と (16)4) S. Yamaki, J. Odagiri, M. Abe and M.とそれぞれ一致している．すなわち，2 変量正Kawamata, ”Eﬀects of Stochastic Phase

規分布に巻き込み分布を仮定した場合の結果は，Spectrum Diﬀerences on Phase-Only Cor-

もとの平面上で表される確率分布を仮定した場

relation Functions Part I: Statistically Con-合の結果と一致することがわかった．

stant Phase Spectrum Diﬀerences for Frequency Indices,” Proceedings of IEEE 3rd5.

まとめInternational Conference on Network In-

2 信号の位相スペクトルを 2 変量正規分布に従うと仮定して，位相スペクトルの変化に対す

frastructure and Digital Content, pp. 360364, September 2012.

る POC 関数の変化を調べた．この時，位相スペ5) S. Yamaki, M. Abe and M. Kawamata,クトル差は 1 変量正規分布に従い，分散は 2 変

Statistical Analysis of Phase-Only Cor-量正規分布の共分散行列を線形変換することでrelation Functions Based on Directional求められる．また，2 信号の位相スペクトルをStatistics,巻き込み 2 変量正規分布に従うと仮定して POC

VOL. E97-A, NO. 12, pp. 2601-2610, DE-関数の期待値と分散を導出した．その結果，2 信

CEMBER 2014.

号の位相スペクトルが巻き込み 2 変量正規分布に従うと仮定したときの POC 関数の期待値と分散は，平面上で与えられる 2 変量正規分布を仮定したときの結果と一致することを証明した．
これより，2 信号の位相スペクトルが巻き込み 2IEICE Trans. Fundamentals,6) KANTIV. MARDIA and PETER E,
JUPP,

Directional Statistics,

John Wi-

ley and Sons Ltd, 2000.
7) N. I. Fisher and A. J. Lee,A correlation

変量正規分布に従う場合の POC 関数の統計的

coeﬃcient for circular data, ” Biometrical.

性質は，方向統計学の概念を用いて説明できる

70. 2, pp. 327-332, 1983.

ことが示された．

–7–