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

信号の周波数変換による位相限定相関関数の挙動の変化
Variability of Behavior of Phase-Only Correlation Functions
under Frequency Transformation of Signals○佐藤有弥，八巻俊輔，吉澤誠○ Yuya Sato, Shunsuke Yamaki, Makoto Yoshizawa
東北大学

Tohoku University
キーワード : 連続時間信号 (continuous-time signals)，位相限定相関関数 (Phase-Only Correlationfunction)，周波数変換 (frequency transformation)

連絡先 : 〒 980-8579 仙台市青葉区荒巻字青葉 6-6-05 東北大学 青葉山キャンパス 電子情報システム・応物系 1 号館 521 号室 吉澤・杉田研究室
佐藤有弥，Tel.: (022)795-7130，Fax.: (022)795-7125，E-mail: yuya.sato.q4@dc.tohoku.ac.jp

1.

まえがき

ない場合の POC 関数の挙動について明らかに

位相限定相関 (POC:Phase-Only Correlation)関数は２つの信号の類似度を評価する関数として知られている．また，POC 関数は２信号の幾何学的な関係，例えば信号の位置ずれなど，を求めることができる．そのため画像マッチン
グ [1]，生体認証アルゴリズム [2]，古いフィル
ムの位置ずれ補正 [3]，周期性を持つ DNA 配列の探索 [4] などに応用されてきた．
２つの信号の位相スペクトルが等しい場合，つまり位相スペクトル差が 0 であるとき POC関数はデルタ関数になる．また，２つの信号の位相スペクトルが異なる場合，つまり位相スペクトル差が 0 ではないときには POC 関数はデルタ関数とは異なる．２つの信号の位相スペクトルが完全に一致するということはまれである
ため，２つの信号の位相スペクトル差が 0 ではする必要があった．文献 [5] では，信号の位相スペクトル差が確率変数であると仮定し，POC関数の期待値と分散を導出することで POC 関数の挙動を統計的に解析してきた．
信号は，連続な時間軸上で定義される連続時間信号と離散的な時間軸上で定義される離散時間信号の二つに分類される．連続時間信号は信号処理を施す場合，任意の周期でサンプリングすることで離散時間信号に変換される．連続時間信号を一定のサンプリング周期でサンプリングして離散時間信号を得るとき，連続時間信号
では (−∞, ∞) で定義されていた位相が離散時間信号においては (−π, π] に制限されてしまう．
離散時間信号における POC 関数の性質については先行研究 [6] などによって解析がなされているが，連続時間信号における POC 関数について–1–は解析がなされていないため，これを考える必

ここで，α(k) = θ(k)−ϕ(k) であり信号の位相ス要がある．また，サンプリングに伴う位相スペペクトル差を表す．つぎに，位相スペクトル差クトルの変化によって，位相スペクトルを確率α(k) を確率変数とすることで，POC 関数の期変数と仮定した場合に従う確率分布が，連続時待値と分散を求める．位相スペクトル差 α(k) が間信号においてはどういった分布であったのか確率分布に従う確率変数であるとすると，POCを離散時間信号の位相スペクトルの分布から復関数の期待値と分散は位相因子の期待値 A を用元することができなくなる．連続時間における

いて次のように表せる．

位相スペクトルの情報を保存するためには，信変化を抑える必要がある．その方法として，本E[r(m)] = Aδ(m)1(1 − |A|2 )Var[r(m)] =N

稿では連続時間における投影法を用いた位相ス

よって，位相因子の期待値 A を求めることがでペクトルの非線形な変換を提案する．

きれば，POC 関数の期待値と分散を計算するこ号のサンプリング時に起こる位相スペクトルの(4)(5)とができる．

2.

位相スペクトルが確率分布に従う場合，確率離散時間における POC 関数

密度関数の特性関数から位相因子の期待値を求

まず，先行研究 [5] で明らかになっている事実めることができる．位相スペクトルの従う確率として，離散時間信号における POC 関数の定義

分布の確率密度関数を p(α(k)) とすると，特性とその統計的性質を述べる．信号長が N である

関数 Ψ(t) は次の式で表される．
Z ∞ejα(k)t p(α(k))dα(k)Ψ(t) =２つの離散時間複素信号をそれぞれ x(n), y(n)とする．これらの離散時間フーリエ変換（DFT）−∞
jα(k)t= E[eX(k), Y (k) はそれぞれ以下の式で表される．

=

相因子の期待値 A = E[ejα(k) ] を計算すること

x(n)WNkn = |X(k)|ejθ(k) (1)n=0=ができる．以上より，離散時間信号の POC 関数の期待値と分散は，位相スペクトル差が確率Y (k) = DFT[y(n)]N−1X(6)

よって，式 (6) において t = 1 とすることで，位X(k) = DFT[x(n)]N−1X]分布に従う場合，その特性関数を計算することy(n)WNkn = |Y (k)|ejϕ(k) (2)で求めることが可能である．

n=0

ここで，WN = exp(−j2π/N ) であり，DFT の回転因子を表す．また，θ(k), ϕ(k) はそれぞれX(k), Y (k) の位相スペクトルである．離散時間

3.

連続時間における POC 関数

3.1

POC 関数の定義

信号の POC 関数は，正規化クロスパワースペ
クトル R(k) を離散フーリエ逆変換したものとして定義される．よって，離散時間信号の POC関数は以下のようになる．

r(m) ==1
N1
N

N−1X
k=0
N−1X
k=0析するため，まず連続時間信号における POC関数を定義する．信号長が T である２つの連続時間複素信号を，f (t), g(t) とする．これらの信R(k)WN−mkejα(k) WN−mk

連続時間信号における POC 関数の挙動を解(3)号のフーリエ級数展開は以下の式で表される．
P∞
jω0 kt (0 < t < T )= =f (t) = Rk=−∞ F (k)e1 Tjω
kt
0F (k) = T 0 f (t)edt= |F (k)|ejθ(k) (−∞ <=k<= ∞)

–2–
P∞
jω0 kt (0 < t < T )
g(t)
=

= =k=−∞ G(k)eR
1 Tjωkt
0G(k) = T 0 g(t)edt= |G(k)|ejϕ(k) (−∞ <=k<= ∞)POC 関数の期待値は次の式で表される．
"#NX1jα(k) jω0 kt
E[r(t)] = Eee
2N + 1k=−N1E[ejα(k) ]ejω0 kt
2N + 1k=−N
(2N +1)ω0 tsin21(t ̸= 0, T )Aω
0t(10)=2N + 1sin
2A (t = 0, T )ここで，ω0 = 2π/T であり，θ(k), ϕ(k) はそれ=ぞれ F (k), G(k) の位相を表している．２つの信号 f (t), g(t) の正規化クロスパワースペクトルR(k) は次のように表される．
R(k) = F̂ (k)Ĝ∗ (k)F (k)G∗ (k)=|F (k)||G∗ (k)|= ejα(k)NX

また，POC 関数の分散は次の式で表される．

Var[r(t)] = E[r(t)r∗ (t)] − E[r(t)]E[r∗ (t)]1=(1 − AA∗ )(11)
2N + 1(7)ここで，α(k) = θ(k) − ϕ(k) であり，位相スペクトル差を表している．

ここで，位相因子の期待値を A としている．位

連続時間信号の POC 関数は正規化クロスパ

相スペクトル差 α(k) を確率密度関数 p(k) に従ワースペクトルをフーリエ係数とするフーリエう確率変数と仮定すると，特性関数 Ψ(s) は次級数展開の形で表される．また，正規化のためのように表せる．
Z ∞Ψ(s) =ejα(k)t p(α(k))dα(k)以下の条件を付加する．
Z ∞|r(t)|2 dt = Const.

−∞jα(k)t(8)= E[e−∞この条件の下，連続時間信号における POC 関数 r(t) を定義すると，以下のようになる．

r(t) ==

12N + 1NX

これは，離散時間信号の場合で考えた特性関数
時間信号の場合と同じように，特性関数を考えることで位相因子の期待値を求めることができる．

k=−N(9)k=−N

3.2(12)と一致するため，連続時間信号においても離散R(k)ejω0 kt

NX1ejα(k) ejω0 kt
2N + 1]

3.3

期待値と分散計算例：位相スペクトル差 α(k) が一様分布に従う場合

位相スペクトル差 α(k) が平均 µ の一様分布

式 (9) で定義される連続時間信号の POC 関数の統計的性質を見ていく．位相スペクトル差α(k) が確率分布に従うと仮定する．このとき，

U (−w + µ, w + µ) に従うものと仮定し，このときの POC 関数の期待値と分散を求める．位相因子の期待値は次のように計算される．

A=sin(w) −jµew(13)

これを POC 関数の期待値と分散の一般式に代入することで，位相スペクトル差が平均 µ の一様分布 U (−w + µ, w + µ) に従う場合の POC 関–3–

(a) 実測値(a) 実測値(b) 理論式(b) 理論式

Fig. 1: 連続時間信号のPOC 関数の期待値

Fig. 2: 連続時間信号のPOC 関数の分散数の期待値と分散の理論式が以下のように求め4.

られる．

位相スペクトルの周波数変換

 4.1 位相スペクトル(2N +1)ω0 t
sin

2sin(w)1
位相スペクトル差は，２つの信号の位相スペ 2N + 1 w
ω0 tsin 2
E[r(t)] =(14) クトルがどれだけ異なっているのかというデー(t ̸= 0, T )タである．連続時間信号において位相スペクト
 sin(w) (t = 0, T )w
ルは (−∞, ∞) の値をとりうる線形確率分布に
1sin(w)Var[r(t)] =1−(15) 従うと考えられる．それに対して離散時間信号
2N + 1wの場合，位相スペクトルは (−π, π] の範囲に制シミュレーションによって得られた実測値と，連限されており，先行研究 [6] では離散時間信号
続時間信号の POC 関数の期待値と分散の理論式の位相スペクトルは巻き込み分布に従うものでがそれぞれ次の図 1(a), 図 1(b), 図 2(a), 図 2(b)あると考えられている．連続時間信号はサンプになる．図 1(a) と図 1(b)，図 2(a) と図 2(b) かリングすることで離散時間信号へと変換できる．
ら，連続時間信号における POC 関数の期待値
よって，この位相スペクトルの従う分布の変化と分散の理論式が実測値に合致しているといえはサンプリングによって生じるものと考えられる．また，期待値のピークが t = 0, T であるこる．離散時間信号の位相スペクトル αd としたと，w が大きくなった場合，期待値は減少する
とき，巻き込み法による位相スペクトルの変換傾向にあり，分散は増加する傾向にあることが分かった．

–4–

式 (17),(19) から，p′ (α̃) の特性関数 Ψ(s) は以下は式 (16) によって与えられる．

αd = (α − π)(mod 2π) − π(16)巻き込み法による位相スペクトルの変換の特徴として，変換が不可逆であることが挙げられの式で表される．
Z ∞Ψ(s) =p′ (α̃)ej α̃s dα̃−∞Z ∞=p(α)ej2 arctan(α)s dα (20)−∞

る．これは，位相スペクトルが従っていた確率分布についても同様に言える. そのため，離散以上のことから，投影法を用いた変換を位相時間信号の位相スペクトルの情報から連続時間スペクトル差に施した場合の位相因子の期待値信号の位相スペクトルの情報を完全には復元でA についても，確率密度関数の特性関数を考えるきない．巻き込み分布に従う位相スペクトルのことによって求めることができる．よって，POC取りうる値の範囲は (−π, π] であり，これはサ関数の期待値と分散の一般式にここで求めた位ンプリング周波数 Fs で正規化されている．し相因子の期待値 A を代入することで，位相スペたがって，位相の変換によって位相スペクトルクトル差を変換した後の POC 関数の期待値との値の取りうる範囲を (−πFs , πFs ) に変換でき分散を求めることができる．

れば，位相スペクトルは連続時間信号の位相スペクトル情報を保存したまま離散時間信号とし計算例：位相スペクトル差が一様分布に従う場合4.3て扱うことができると考えられる．位相スペク
トルの値の取りうる範囲を変換する方法として，以下の式による周波数変換を提案する．

位相スペクトル差に周波数変換を施さない場
合と同様に，位相スペクトル差 α(k) が平均 µα̃(k) = 2 arctan(α(k)) (−π < α̃(k) < π) (17)の一様分布 U (−w + µ, w + µ) に従うものと仮定し，このときの POC 関数の期待値と分散を4.2求める．位相因子の期待値は次のように計算さ投影法を用いた変換れる．

式 (17) で与えられる変換を行った場合の POC

A=関数の挙動について見ていく．変換を行った後の POC 関数は以下の式になる．

r(t) =NX
1ej α̃(k) ejω0 kt
2N + 1(18)k=−Narctan(µ + w) + arctan(µ − w)−1w

1
1 + (µ + w)2+jlog(21)2w
1 + (µ − w)2

このとき，µ = 0 とすると，虚部が 0 となるので，次のようになる．

変換前後で POC 関数の変化は位相スペクトル

A=差 α(k) が α̃(k) となったことのみであるため，この周波数変換による変化を調べるためには位相因子の期待値 A の変化を考えればよい．位相2 arctan(w)−1w(22)

これを，POC 関数の一般式に代入すると，期待値と分散の理論式は以下のように求められる．

因子の期待値を考えるために，位相スペクトルの変化について考える．α̃(k) の確率密度関数をp′ (α̃) とするとこれは α, p(α) を用いて以下のように表される．

p′ (α̃) = p(α)

1 + α22(19)–5–E[r(t)] = sin (2N +1)ω0 t
22 arctan(w)1−1 2N + 1ωw
sin 20 t(t ̸= 0, T )
2arctan(w)− 1 (t = 0, T )w

(a) 実測値(a) 実測値(b) 理論式(b) 理論式

Fig. 3: 連続時間信号のPOC 関数の期待値
1Var[r(t)] =2N + 1
1−2 arctan(w)−1w

Fig. 4: 周波数変換後のPOC 関数の分散
2 !

は w の値の増加に伴い単調に増加する．w > w1
では，w の値の増加に伴い，期待値は単調減少

シミュレーションによって得られた実測値と，し分散は単調増加することがわかる．

連続時間信号の POC 関数の期待値と分散の理論式がそれぞれ次の図 3(a), 図 3(b), 図 4(a), 図5.

4(b) になる．図 3(a) と図 3(b)，図 4(a) と図 4(b)から，周波数変換後の POC 関数の期待値と分散の理論式が実測値に合致しているといえる．またこの図から，周波数変換を施した場合の POC関数の期待値と分散の挙動は周波数変換を施していない場合と異なり，w の値の増加に伴って増減しないことがわかる．図 3(b),4(b) の t = 0 のとき，POC 関数が w の値に伴って変化する様子を表した図が図 5(a),5(b) である．この図から，

POC 関数の期待値の最小値を与える w の値と分散の最大値を与える w の値は等しいことがわかる．この w の値を w1 とすると，w <= w1 では

おわりに

本研究では，連続時間信号をサンプリングして得られる離散時間信号の POC 関数の挙動について調べる前段階として連続時間信号における POC 関数の挙動の解析と，位相スペクトル差の巻き込み変換による影響を回避する手段としての位相スペクトル差の非線形変換を考えた．
今後の研究として，位相スペクトルの周波数変換を行った連続時間信号と周波数変換を行わない連続時間信号のそれぞれをサンプリングして離散時間信号にした場合の POC 関数の挙動について調べる必要がある．

w の値の増加に伴い，期待値は単調減少し分散–6–

spectrum diﬀerences for frequency indices,”2012 3rd IEEE International Conference on
Network Infrastructure and Digital Content,Beijing, 2012, pp. 360-364.
6) Shunsuke
YAMAKI,
Masahide
ABE,
Masayuki KAWAMATA, Statistical Analysis of Phase-Only Correlation Functions
with Phase-Spectrum Diﬀerences Following
Wrapped Distributions, IEICE Transactions
on Fundamentals of Electronics, Communications and Computer Sciences, 2016, Volume
E99.A, Issue 10, Pages 1790-1798, ReleasedOctober 01, 2016.

(a) 期待値 (t = 0)(b) 分散 (t = 0)

Fig. 5: 周波数変換後の POC 関数の w に伴う変化

参考文献
1) C. D. KUGLIN. The phase correlation imagealignment mathod. Proc. Int. Conf. on Cibernetics and Society, 1975, pages 163-165, 1975.
2) K. Miyazawa, K. Ito, T. Aoki, K. Kobayashi
and H. Nakajima, ”An Eﬀective Approachfor Iris Recognition Using Phase-Based Image
Matching,” in IEEE Transactions on PatternAnalysis and Machine Intelligence, vol. 30, no.
10, pp. 1741-1756, Oct. 2008.
3) Mizuki HAGIWARA, Masahide ABE,Masayuki KAWAMATA, Estimation Methodof Frame Displacement for Old Films Using
Phase-Only Correlation, Journal of Signal
Processing, 2004, 8 巻, 5 号, p. 421-429.
4) A. K. Brodzik, ”Phase-only ﬁltering for themasses (of DNA Data): a new approach to sequence alignment,” in IEEE Transactions onSignal Processing, vol. 54, no. 6, pp. 24562466, June 2006.
5) S. Yamaki, J. Odagiri, M. Abe and M. Kawamata, ”Eﬀects of stochastic phase spectrum
diﬀerences on phase-only correlation functions: Part I: Statistically constant phase–7–