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

位相限定相関関数の確率分布の数理表現
Mathematical Expressions of Probability Distributions ofPhase-Only Correlation Functions○山口 孝志，八巻 俊輔，吉澤 誠○ Takashi Yamaguchi, Shunsuke Yamaki, Makoto Yoshizawa
東北大学

Tohoku University
キーワード : 位相限定相関関数 (phase-only correlation functions)，確率分布 (probability distributions)
連絡先 : 〒 980-8579 仙台市青葉区荒巻字青葉 6-6-05
東北大学大学院 医工学研究科 吉澤・杉田研究室
山口孝志，Tel.:022-795-7130，E-mail: takashi.yamaguchi.s8@dc.tohoku.ac.jp

1.

はじめに

別するための閾値を設ける必要がある．ピーク

位相限定相関 (POC:Phase-Only Correlation)
関数は，2 つの信号の類似度を評価する関数として知られている．そのため，画像マッチング
1, 2, 3, 4, 5, 6) など，幅広い分野で応用されてきた．これらの応用は，2 つの信号の位相スペク
トルが完全に等しいとき，すなわち位相スペク
トル差が 0 であるとき，POC 関数がデルタ関数になる性質を根拠としている．
しかしながら，実際の信号において 2 つの信号の位相スペクトルが完全に等しくなることはほとんどありえない．例えば，実際の画像マッ
チングにおいて，基準画像と観測画像が同一でとサイドローブはとりうる値の範囲が同じであるため，ピークとサイドローブを完全に区別するような閾値を決めることは不可能である．そのため，ピークをサイドローブと判断する，またはサイドローブをピークと判断するような誤りを最小にするような閾値の導出が求められる．
このような閾値を導出するためには POC 関数の確率分布を導出する必要がある．本稿では位相スペクトル差を線形の確率分布に従う確率変数であると仮定し，信号の点数を変化させたと
きの POC 関数の実部と虚部の確率分布を数値計算によって求めた．

あることはめったにない．このような場合 2 つの信号の位相スペクトルが完全に等しいときに

2.

POC 関数がデルタ関数になるという性質をマッ
チングに用いることはできない．

位相限定相関 (POC) 関数

本章では POC 関数の定義を示す．長さ N の2 つの複素離散信号 a(n) および b(n) を考える．

実際の信号の類似度を評価する関数として POC関数を用いる場合，ピークとサイドローブを区–1–

これらの信号の離散フーリエ変換はそれぞれ，

3.

A(k) = DFT[a(n)]
=

N
−1
∑

本章では，信号マッチングで求められている

a(n)WNnk

POC 関数のピークとサイドローブを分ける閾n=0= |A(k)|e

jθk

(1)

値を求めるためには，POC 関数の確率分布の導出が必要であることを説明する．まず，ピーク(k = 0, 1, · · · N − 1)とサイドローブがとりうる値の傾向について実B(k) = DFT[b(n)]=POC 関数の確率分布を導出す
る意義

験的に調べる．図 1 は N = 8 の POC 関数に関N
−1∑b(n)WNnk

して，POC 関数の定義式 (3) の位相スペクトルn=0= |B(k)|ejϕk(2)(k = 0, 1, · · · N − 1)

差 αk に正規分布 N (0, σ 2 ) に従う乱数を代入し，
ピーク r(0) とサイドローブ r(m ̸= 0) をそれぞれ複素平面にプロットする試行を 1000 回行っで与えられる．ここで j は虚数単位で，WN =た散布図である．図 1(a) と図 1(b)，図 1(c)，図exp(−j2π/N ) は離散フーリエ変換の回転因子で1(d) はそれぞれ σ 2 = 0, 0.25, 0.5, 1 としたとき

ある．|A(k)| と |B(k)| はそれぞれ信号 a(n) との散布図である．図 1 より位相スペクトル差のb(n) の振幅スペクトルで，θk と ϕk はそれぞれ信分散 σ 2 が増加するにつれて，POC 関数のピー

号 a(n) と b(n) の位相スペクトルである．このとクとサイドローブの分布の重なりが大きくなる

き 2 つの信号 a(n) および b(n) の POC 関数 r(m)ことがわかった．とりうる値の範囲が重なってはクロスパワースペクトル A(k)B ∗ (k)(“*” は複いるピークとサイドローブを分けることを考え

素共役を表す) を振幅スペクトル |A(k)||B(k)| でるときには，ピークをサイドローブと判断して割ることにより得られる正規化クロスパワース

しまう，もしくはサイドローブをピーク判断しペクトルの離散フーリエ逆変換 (IDFT) としててしまうという誤りを最小とするような閾値の以下で与えられる．

設定が求められる．

[r(m) = IDFT
=

=

1
N
1
N

N
−1
∑
k=0
N
−1
∑

A(k)B ∗ (k)|A(k)||B(k)|]同様の閾値の設定をディジタル通信の基底帯
域信号伝送で行われているため紹介する．基底帯域信号伝送では時間幅 Ts の矩形パルスが Ts

A(k)B ∗ (k)W −mk|A(k)||B(k)| N
ejαk WN−mk

ごとに送信されている．その振幅は，送信シン
ボルが 1 のときには +1，送信シンボルが 0 の

(3)

k=0

ときには −1 である．このとき送信系統は振幅0 を使わない NRZ 矩形パルス系列である．受信(m = 0, 1, · · · , N − 1)

波形は矩形 NRZ 信号に雑音が重畳されており，
ここで，αk = θk − ϕk は 2 つの信号の位相

受信信号は判定スレッショルドレベル Vref と比スペクトル差であり，ejαk は位相因子と呼ばれる．式 (3) より，位相スペクトル差 αk が決定す

較され，Vref より大きければ 1，小さければ 0と判定され，受信シンボルが生成される 7) ．本

れば，POC 関数の値は具体的に定まる．

来は +1 と −1 で完全に分かれていた矩形パル
スが雑音が重畳することで同じ範囲に値を持つようになっている．これは POC 関数のピークとサイドローブと同じであり，判定スレッショ–2–

ルドレベル Vref がピークとサイドローブを分け
1

るための閾値に相当している．このときディジ

peak r(0)sidelobe r(m≠ 0)

タル通信では誤り率を最小にするような Vref の0.5

Im[r(m)]決定の際に，確率密度関数を用いている．図 2
は，NRZ 信号伝送における受信シンボルの確率0密度関数である．送信シンボルは +1 と −1 し−0.5

かとりえないが，雑音が重畳することでそれぞれ，+1 と −1 を中心とする広がりを持つように

−1
−1

−0.5

なっていることがわかる．図 2 では雑音として

0
Re[r(m)]0.51(a) σ 2 = 0

正規分布 N (0, 1) を仮定している．図 2 におけ

1peak r(0)sidelobe r(m≠ 0)

る斜線部はシンボル誤りの確率を示しているた
め，誤り率を最小にするような Vref は斜線部の0.5

Im[r(m)]

面積を最小にするように決められる．以上のことからピークとサイドローブを分けるときに誤り率を最小にするような閾値を設定するために，

0

−0.5

POC 関数のピークとサイドローブそれぞれの確率分布を考える．

−1
−1

−0.5

0
Re[r(m)]0.51(b) σ 2 = 0.25

4.

POC 関数の確率分布の導出

4.1

まえがき

1
peak r(0)sidelobe r(m≠ 0)

0.5

Im[r(m)]

本章では POC 関数の確率分布の導出を行う．
まず確率分布の導出に必要となる基礎として，連続型の確率分布の定義と連続型の確率変数の

0

−0.5

変換，確率変数の和の分布について示す．そし−1−1て POC 関数の位相スペクトル差 αk が確率変数−0.5

0
Re[r(m)]0.51(c) σ 2 = 0.5

であると仮定したときの N = 1 のときの POC
1関数の確率分布を導出し，ヒストグラムとの比

peak r(0)sidelobe r(m≠ 0)較から導出した確率分布が妥当であることを示

0.5

Im[r(m)]す．次に N > 1 のときの POC 関数の確率分
布を導出する．具体例として N = 2 のときの0

POC 関数の確率分布を導出し，ヒストグラムと−0.5の比較から確率分布が妥当であることを示す．

−1
−1

4.2

確率分布について−0.5

0Re[r(m)]0.51(d) σ 2 = 1

本節では，POC 関数の確率分布の導出に必要となる基礎として，連続型の確率分布の定義と連続型の確率変数の変換，確率変数の和の分布

Fig. 1: 位相スペクトル差の分散 σ 2 を変化させたときの POC 関数のピークとサイドローブの散布図について示す．

–3–である．∆x → 0 とすると，

fY (y) = fX (x)

(
) dφ
dx
= fX φ−1 (y)
dy
dy

(8)

となる．同様に考えて，φ は単調減少で微分可能の場合，y − ∆y = φ(x + ∆X) とすると，x <X<= x + ∆x となる確率と y − ∆y < Y <
=yと

ܸ௥௘௙

Fig. 2: NRZ 信号伝送の受信シンボルの確率密度関数なる確率は同一で，

fY (y)(−∆y) = fX (x)∆x(9)

である．∆x → 0 とすると，

4.2.1

連続型の確率分布の定義 8)

fY (y) = −fX (x)確率変数 X が連続する変数 (連続型の確率変

(
) dφ
dx
= −fX φ−1 (y)
(10)
dy
dy

数) のとき，X が x から x + ∆x の間に入る確率

よって，φ は微分可能の場合，式 (8) と式 (10)P (x < X <= x + ∆x) を考える．この値は ∆x ををまとめて以下の式が成り立つ．

小さくすると 0 へ収束してしまうので，∆x で

fY (y) = fX (x)割って ∆x → 0 とした極限を fX (x)，すなわち，P (x < X <= x + ∆x)
∆x→0∆x

fX (x) = lim

(4)

4.2.3

(
) dφ
dx
= fX φ−1 (y)
(11)
dy
dy

確率変数の和の確率分布を考える．fX (x) は確率密度関数とよばれ，X がX を 2 つの変数 X1 と X2 の和 X = X1 +X2 とxmin と xmax (xmin と xmax は xmin < xmax をすると，この変数の確率密度関数 fX (x) は，X1

満足する任意の定数) の間に入る確率は，fX (x)と X2 が独立の場合，X1 の確率密度関数 fX1 (x)を定積分してと X2 の確率密度関数 fX2 (x) を用いて，∫ ∞fX (x) =fX1 (τ )fX2 (x − τ )dτ(12)∫ xmaxP (xmin < x <= xmax ) =fX (x)dx

(5)

xmin

−∞

で与えられる．

4.2.2となる 8) ．
ここで Xi (i = 1, 2, · · · ) が従う分布の特性関

連続型の確率変数の変換

数 CXi (t) は確率密度関数 fXi (x) を用いて
∫ ∞ejtx fXi (x)dx(13)
CXi (t) =確率変数 Y が X の関数，すなわち，

Y = φ(X)の場合を考える．φ は単調増加で微分可能とする．この場合，関数 φ の逆関数 φ−1 が存在し，

X = φ−1 (Y ) となる．X ，Y の確率密度関数をそれぞれ，fX (x)，fY (y) とする．y = φ(x)，y+∆y = φ(x+∆X) とすると，x < X <= x+∆xとなる確率と y < Y <= y + ∆y となる確率は同一で，fY (y)∆y = fX (x)∆x

−∞

(6)

(7)

で与えられる．ここで X = X1 + X2 が従う分布の特性関数 CX (t) を求めることを考える．X の確率密度関数 fX (z) を用いて式 (13) と同様に，∫ ∞ejtx fX (x)dx(14)
CX (t) =−∞とかける．式 (14) の右辺に式 (12) を代入して以下の式が得られる．
)
(∫ ∞
∫ ∞
jtx
fX1 (τ )fX2 (x − τ )dτ dx
e

–4–

−∞

−∞

∫ ∞∫ ∞=−∞−∞

ejtτ fX1 (τ )ejt(x−τ ) fX2 (x − τ )dτ dx(15)ここで υ = x − τ と変数変換すると，式 (15) は，∫ ∞∫ ∞
ejtτ fX1 (τ )ejtυ fX2 (υ)dτ dυ−∞ −∞) (∫ ∞
)
(∫ ∞
jtυ
jtτ
e fX2 (υ)dυ=e fX1 (τ )dτ−∞fα (α) で表す．N = 1 のときの POC 関数の実部と虚部の確率分布を導出するためには確率変数 α0 に対して，実部は cos(α0 )，虚部は sin(α0 )にそれぞれ確率変数の変換を行う必要がある．

4.3.2

−∞

= CX1 (t)CX2 (t)POC 関数の実部の確率分布の導出

確率変数 α0 を cos(α0 ) に変換する．確率変(16)

数 α0 の範囲は −∞ <= α0 <= ∞ である．α0 とと変形される．式 (16) より X が従う分布の特性

cos(α0 ) の関係を図 3 に示す．図 3 より，ひとつ関数 CX (t) は X1 と X2 が従う分布の特性関数の cos(θ) の値と cos(θ +2πq) と cos(−θ −2πq)(qの積で表されることがわかった．つまり，X =は任意の整数) は等しいことがわかる．このこと

X1 + X2 の確率密度関数 fX (x) を得るためにはから，cos(α0 ) と α0 は 1 対 1 対応でないことがわ

CX1 (t) と CX2 (t) の積を逆変換すればよい．よっ

かる．1 対 1 対応でないと逆関数が存在できないて以下の式が得られる．
∫ ∞1e−jtx CX1 (t)CX2 (t)dt
fX (x) =2π −∞ため，cos(α0 ) が単調増加する範囲 (−π+2πq) <

(17)

α0 <
= 2πq(q は任意の整数) と単調減少する範囲
2πq < α0 <= (π + 2πq)(q は任意の整数) で範囲ごとに変換を行う必要がある．そこで新しい確

4.3

N = 1 のときの POC 関数の確率分
布の導出

本節では最も基本的な問題設定として N = 1

率変数として α̃0 (0 <= α̃0 <= π) を定義し，α0 を

α̃0 に変換することを考える．α̃0 と α0 の関係は以下の式で表される．

のときの POC 関数の確率分布の導出を行う．はじめに，問題提起として N = 1 のときの POC関数を与える．次に，位相スペクトル差が確率変数と仮定した時の POC 関数の実部と虚部の確率分布の導出をする．そして，導出した確率分布が妥当かどうかを確率分布とヒストグラムの比較によって示す．

4.3.1

α̃0 = |arg(ejα0 )|jα0
ここで −π <= arg(e ) < π である．図 4 は α̃0と α0 の関係を表している．図 4 からわかるよ
うに cos(α̃0 ) と同じ値になる cos(α0 ) の α0 が α̃0にまとめられていることがわかる．α0 から α̃0への変換によって α̃0 と cos(α̃0 ) は 1 対 1 対応となるので確率変数の変換を行うことができる．

α̃0 の確率密度関数である fα̃ (α̃) は，α0 の確問題提起

率密度関数 fα (α) を用いて，以下の式で表すこN = 1 のときの POC 関数 r(m) は式 (3) よとができる．

り，m = 0 および k = 0 のみにおいて定義され

fα̃ (α̃) =る．N = 1 のときの POC 関数 r(0) は位相スペ

(18)

∞
∑
(

fα (α̃ + 2πq)q=−∞)+fα (−α̃ − 2πq)

クトル差 α0 を用いて以下の式で表される．

r(0) = cos(α0 ) + j sin(α0 )

(19)

(20)

次に α̃0 を cos(α̃0 ) に変換する．確率変数 Y =cos(α̃0 ) とし，確率変数 Y の確率密度関数を fY (y)位相スペクトル差 α0 を確率変数と仮定した時とおくと，式 (11) より以下の式が成り立つ．

の N = 1 のときの POC 関数 r(0) の実部と虚部の確率分布を導出する．α0 の確率密度関数は–5–fY (y) = fα̃ (cos−1 (y))

d cos−1 (y)
dy

(21)

ቑຍ

ቑຍ

ῶᑡ

ቑຍ

ῶᑡ

図 6 は α0 と α̂0 の関係である．図 6 より sin(α̂0 )ῶᑡと同じ値になる sin(α0 ) の α0 が α̂0 にまとめら ߠ

れていることがわかる．この変換によって α̂0 とെߠ െ ʹߨ ߠ െ ʹߨߠെߠെߠ ൅ ʹߨsin(α̂0 ) は 1 対 1 対応なので確率変数の変換を行ߠ ൅ ʹߨ

うことができる．
ῶᑡ

ቑຍ

ῶᑡ

ቑຍ

ῶᑡ

ߠ െ ʹߨ െߠ െ ߨߠെߠ ൅ ߨ

ቑຍ ῶᑡ

Fig. 3: α0 と cos(α0 ) の関係
α̃0πߠ ൅ ʹߨ െߠ ൅ ͵ߨ

π/2
0
−3π

−2π

−π

0

π

2π

α0

Fig. 4: α0 と α̃0 の関係
Fig. 5: α0 と sin(α0 ) の関係

このとき，

d cos−1 (y)1=√dy
1 − y2

αˆ0

π/2
0
−π/2
−3π

(22)

−2π

−π

0

π

2π

α0

であるため，式 (21) に式 (20) と式 (22) を代入

Fig. 6: α0 と α̂0 の関係

して以下の式が導出できる．

α̂0 の確率密度関数である fα̂ (α̂) は，α0 の確
fY (y)=

1

率密度関数を fα (α) を用いて，以下の式で表す

∞
∑
(

√
fα (cos−1 (y) + 2πq)
1 − y 2 q=−∞)+fα (− cos−1 (y) − 2πq)(23)ことができる．

fα̂ (α̂) =∞∑(fα (α̂ + 2πq)q=−∞+fα (−α̂ − π − 2πq)以上より，N = 1 のときの POC 関数の実部の確率分布を式 (23) で表すことができた．

4.3.3

)

(25)

次に確率変数 Z = sin(α̂0 ) の確率密度関数 fZ (z)を式 (11) を用いて導出すると以下のようになる．

POC 関数の虚部の確率分布の導出

fZ (z) = fα̂ (sin−1 (z))

本小節では 4.3.2 と同じ方法で確率変数 α0 を

sin(α0 ) に変換することで，N = 1 での POC 関

d sin−1 (z)
dz

(26)

上式において，1d sin−1 (z)=√dz
1 − z2

数の虚部の確率分布を導出することを考える．
図 5 に α0 と sin(α0 ) の関係を示す．図 5 より，(27)cos(α0 ) と同様に sin(α0 ) に単調増加の範囲と単であり，式 (26) に式 (25) と式 (27) を代入して調減少の範囲が周期的に現れることから，α0 と以下の式が導出できる．

sin(α0 ) は 1 対 1 対応ではないためそのまま確率

fZ (z)

変数の変換を行うことはできない．そのため新=π
しい確率変数として α̂0 (− π2 <= α̂0 <
= 2 ) を定義

し，α0 を α̂0 に変換することを考える．α0 と α̂0∞∑(1fα (sin−1 (z) + 2πq)21 − z q=−∞)+fα (− sin−1 (z) − π − 2πq)(28)√の関係は次式で与えられる．

α̂0 = |arg(ej(α0 − π2 )

π
)| −
2

以上より，N = 1 のときの POC 関数の虚部の(24)確率分布を式 (28) として表すことができた．

–6–

4.3.4

導出した確率分布とヒストグラムの比較

10frequency

位相スペクトル差の確率変数 α0 の確率密度関数 fα (α) を用いて，N = 1 のときの POC 関数の実部と虚部の確率密度関数はそれぞれ式 (23)64と式 (28) であると前小節までに導出した．本

2

小節では導出した確率分布とヒストグラムの比0−1較を行う．位相スペクトル差の確率変数 α0 が確率密度関数 fα (α) は20Re[r(0)]

0.5

1

2.5

]

histogram
pdf

2

(29)
frequency1fα (α) = √ exp
2π

−α2

−0.5

Fig. 7: N = 1 のときの POC 関数の実部の確率分布の理論式とヒストグラムの比較

N (0, 1) の正規分布に従うと仮定すると，α0 の[histogram
pdf

8となる．このとき POC 関数の実部と虚部の確

1.5
1率密度関数の導出式である式 (23) に式 (29) を0.5代入した式に対し数値解析を行った．確率密度

0
−1

−0.5

0
Im[r(0)]0.51関数と比較検討するためにヒストグラムを用い

Fig. 8: N = 1 のときの POC 関数の虚部の確率分布の理論式とヒストグラムの比較る．図 7 と図 8 はヒストグラムとそれぞれ POC
関数の実部と虚部の確率密度関数の導出式の数値解析結果である．ヒストグラムは面積が 1 になるように正規化している．図 7 と図 8 より，

4.4.1

問題提起

導出した N = 1 のときの POC 関数の実部と虚N = 1 のときと同様に位相スペクトル差 αk部の確率分布がそれぞれのヒストグラムをよくを確率変数と仮定した時の POC 関数の実部と近似していることがわかる．以上より導出した虚部の確率分布を導出する．N > 1 のときのN = 1 のときの POC 関数の実部の確率密度関

POC 関数 r(m) は式 (3) で表される．たとえば数である式 (23) と虚部の確率密度関数である式

N = 2 のときの POC 関数 r(m)(m = 0, 1) は式(28) は妥当であると考えられる．

(3) より以下の式で表すことができる．
1 ∑ jαke
21

4.4r(0) =

N > 1 のときの POC 関数の確率分
布の数値計算による導出

k=0

1{
=
(cos(α0 ) + cos(α1 ))2}+j (sin(α0 ) + sin(α1 ))本節では N > 1 のときの POC 関数の確率分(30)

布を導出する．まず，N = 1 のときの POC 関数の実部と虚部の確率分布の特性関数を求める．その特性関数を応用して，サイドローブ r(m ̸= 0)

k=0

1{
=
(cos(α0 ) − cos(α1 ))2}+j (sin(α0 ) − sin(α1 ))を持つ最小の偶数として N = 2 における POC
関数の確率分布を数値計算によって導出し，ヒストグラムと比較して確率分布が妥当であることを確かめる．

1 ∑ jαk −jπke e21r(1) =(31)

式 (30) と式 (31) からもわかるように，N > 1のときの POC 関数の実部と虚部の確率分布を求めるためには確率変数の和の確率分布を求め，–7–

その確率変数を N で割る必要がある．つまり，

y = cos(α̃)(0 < α̃ <= π) に変数変換をすると，

式 (30) の cos(α0 ) や cos(α1 ) のような，POC

式 (33) の右辺は
∫ π
eit cos(α̃) fY (cos(α̃)) sin(α̃)dα̃

関数 r(m) のそれぞれの項の確率変数 Xi (i =0, 1, · · · , N − 1) を用いて，一般に POC 関数

0r(m) の確率変数は (X0 + X1 + · · · + XN −1 )/N=で表され，この確率変数の確率分布を求める必る．確率変数を N で割ったときの確率分布は

4.2.2 の確率変数の変換を用いて求めることができる．確率変数 X ，Y が Y = X/N (N は任意(
eit cos(α̃) fα (α̃ + 2πq)q=−∞ 0)+fα (−α̃ − 2πq) dα̃

要がある．確率変数の和の確率分布は式 (17) より特性関数の積を逆変換することにより得られ

∞ ∫ π
∑

と変形できる．次に式 (34)に関してα̃ = α −2πq(−∞ <= α <= ∞) に変数変換をすることによりの自然数) の関係があり，X ，Y の確率密度関数∞ ∫ 2πq+π∑をそれぞれ fX (x)，fY (y) で表せるとき式 (11)q=−∞ 2πq

ejt cos(α−2πq)(より以下の関係が成り立つ．

fY (y) = N × fX (N y)

(34)

=

(32)

)
fα (α) + fα (−α) dα∞ ∫ 2πq+π∑ejt cos(α) fα (α)dαq=−∞ 2πq−π∫ ∞N > 1 のときの POC 関数の確率分布は，まず=−∞それぞれの項の特性関数を求め，次に特性関数

ejt cos(α) fα (α)dα(35)の積の逆変換で確率分布を求め，最後に確率変を得る．以上より式 (23) の特性関数は位相スペ数を N で割る確率変数の変換を行うという手順

クトル差の確率密度関数 fα (α) を用いて式 (35)でこれから求める．

のように表すことができた．
同様に式 (28) の特性関数 CZ (t)(−∞ <
=
= t <

4.4.2

∞) を導出する．式 (13) と式 (28) より特性関数

特性関数の導出N = 1 のときの POC 関数の実部の確率密度関数である式 (23) と虚部の確率密度関数である
式 (28) の特性関数をこれから導出する．式 (23)と式 (28) の特性関数を導出するのは，N > 1のとき POC 関数の項は cos(αk + θ) と sin(αk +は以下の式で表される．
∫ 1
CZ (t) =ejtz fZ (z)dzπy = sin(α̂)(− π2 < α̂ <= 2 ) に変数変換をすると式 (36) は∫ πθ)(k = 0, 1, · · · , N −1)(θ は任意の角度 [rad]) で2構成されており，cos(α0 ) の確率密度関数である

π2

式 (23) と sin(α0 ) の確率密度関数である式 (28)=の特性関数を応用して POC 関数の他の項の確導出する．式 (13) と式 (23) より特性関数は以下の式で表される．
∫ 1
CY (t) =ejty fY (y)dy

eit sin(α̂) fZ (sin(α̂)) cos(α̂)dα̂

∞ ∫ π
∑
2
π
q=−∞ − 2(
eit sin(α̂) fα (α̂ + 2πq))+fα (−α̂ − π − 2πq) dα̂

率密度関数の特性関数も導出できるからである．
まず式 (23) の特性関数 CY (t)(−∞ <=t<= ∞) を

(36)

−1

(37)

と変形できる．次に式 (37)に関してα̂ = α − π −

(33)

−1

–8–

2πq(−∞ <=α<= ∞) に変数変換をすると，グラムの比較を行っている．ヒストグラムは面

ejt sin(α− 2 −2πq)π

積が 1 になるように正規化している．図 9 と図

q=−∞ 2πq

(

=
=

π
π )
fα (α − ) + fα (−α − ) dα22
∞ ∫ 2πq+ π∑2ejt sin(α) fα (α)dα3πq=−∞ 2πq− 2∫ ∞jt sin(α)

e−∞fα (α)dα10，図 11，図 12 より N = 2 のときの POC 関
数の確率分布がそれぞれのヒストグラムをよく近似できていることがわかる．
2.5

(38)

2
frequency∞ ∫ 2πq+π∑率分布を示しており，いずれも計算式とヒストを得る．以上より式 (28) の特性関数は位相スペクトル差の確率密度関数 fα (α) を用いて式 (38)

1.5
1
0.5
0

のように表すことができた．

4.4.3

histogram
pdf

−0.5
−1

N = 2 のときの POC 関数の確率分布の数値計算による導出−0.50Re[r(0)]0.5

1Fig. 9: N = 2 のときの POC 関数 r(0) の実部の確率分布の計算式とヒストグラムの比較
1.4

4.4.2 で導出した式 (35) と式 (38) は式変形に1.2

よってこれ以上簡単な形に変形するのが困難で

histogram
pdf

frequency

1あるため，本小節では式 (35) と式 (38) の積分
を区分求積を用いて数値計算することにより確

0.8
0.6
0.4

率分布を導出することを考える．具体的には，

0.2
0
−1

サイドローブ r(m ̸= 0) が現れる最小の点数として N = 2 のときの POC 関数の確率分布を−0.50
Im[r(0)]0.5

1Fig. 10: N = 2 のときの POC 関数 r(0) の虚部の確率分布の計算式とヒストグラムの比較

数値計算によって導出する．N = 2 のときのPOC 関数 r(m)(m = 0, 1) はそれぞれ式 (30) と式 (31) で表される．αk (k = 0, 1) が同一かつ独立な分布に従うと仮定すると，それぞれの項の特性関数は式 (35) と式 (38) のどちらかで求め

5.

むすび

ることができる．この特性関数の積を式 (17) に

本稿では，位相スペクトル差が線形の確率分

基づいて逆変換を行い，得られた確率分布を式

布であると仮定したときの，POC 関数の実部と(32) に基づいて確率変数の変換を行えば，POC

虚部の確率分布を求めた．まず，N = 1 のとき

関数 r(m)(m = 0, 1) の実部と虚部の確率分布の POC 関数の実部と虚部の確率分布を位相スを求めることができる．具体例として確率変数

ペクトル差の確率密度関数を用いて定式化した．

αk (k = 0, 1) が N (0, 1) に従うと仮定した時の数

次に N = 1 のときの POC 関数の実部と虚部そ値計算を行う．図 9 と図 10 がそれぞれ N = 2 のれぞれの確率密度関数から求めた特性関数を用ときの POC 関数のピーク r(0) の実部 Re[r(0)]いて，N = 2 のときの POC 関数の実部と虚部と虚部 Im[r(0)] の確率分布を示しており，いずの確率分布を数値計算によって求めた．求めた

れも計算式とヒストグラムの比較を行っている．

確率分布についてはそれぞれに対応するヒスト

同様に，図 11 と図 12 がそれぞれサイドローブ

グラムから妥当性を示した．

r(1) における実部 Re[r(1)] と虚部 Im[r(1)] の確–9–

今後の研究として，N > 2 のときの POC 関数

IEEE Trans. Image Process., vol.5, no.8,

5
histogram
pdf

frequency

4

pp.1266-1271, Aug. 1996.

34) G. Wolberg and S. Zokai, “Robust im-2

age registration using logpolar trans-

1form,” Proc. IEEE Int. Conf. Image Pro-

0
−1

−0.5

0
Re[r(1)]

0.5

cess. (IEEE ICIP), pp.493-496, Vancou-

1ver, Canade, Sept. 2000.
Fig. 11: N = 2 のときの POC 関数 r(1) の実部の確率分布の計算式とヒストグラムの比較“Extention of phase correlation to sub-

1.4
1.2

5) H. Foroosh, J. Zerubia, and M. Berthod,

histogram
pdf

pixel registration,” IEEE Trans. Image

frequency

1

Process., vol.11, no.3, pp.188-200, March

0.8
0.6

2002.

0.4
0.2
0
−1

6) L.Chen and K.H. Yap, “An eﬀective tech−0.50Im[r(1)]0.5

1nique for subpixel image registration un-

Fig. 12: N = 2 のときの POC 関数 r(1) の虚部の確率分布の計算式とヒストグラムの比較

der noisy condition,” IEEE Trans. Syst.
Man. Cybern. A, Syst. Humans, vol.38,
no.4, pp.881-887, July 2008.

の確率分布を数値計算によって導出することで，POC 関数の確率分布の傾向を確認する．そして，一般的な点数 N の POC 関数について確率分布
を数値計算ではなく一般式として導出する．

7) 鈴木博, “ビット誤り率,” 電子 · 通信工学
=EKR-10 ディジタル通信の基礎 -ディジ
タル変復調による信号伝送-, pp.76-78, 株
式会社数理工学社, 2012.

8) 縄田和満, “連続型の確率分布,” 東京大学

参考文献

工学教程 基礎系数学 確率 · 統計 I, 東京大

1) C.D. Kuglin and D.C. Hines, “The phasecorrelation image alignment method,”
Proc. Int. Conf. Cybernetics and Society,pp.163-165, 1975.
2) Q. Chen, M. Defrise, and F. Deconinck,“Symmetric Phase-only matched filteringof fourier-mellin transforms for image registration and recognition,” IEEE Trans.
Pattern. Anal. Mach. Intell., vol.16, no.12,
pp.1156-1168, Dec. 1994.
3) B.S. Reddy and B.N. Chatterji, “An ﬀtbased technique for translation, rotation,and scale-invariant image registration,”– 10 –

学工学教程編纂委員会, 丸善出版株式会社,
平成 25 年.