sequence/correlation.hで用いている計算式と理論
(Formula and theory used in sequence/correlation.h)
1. 定義域の異なる時系列データ間の相互相関関数の定義
(Definitions of cross correlation functions (CCFs)
between time series data of different definition ranges)
相互相関関数の計算に用いる2つの時系列データを\(f(t)\), \(g(t)\)とする。
\(f(t)\)の定義域を\(t_{min}^f \leq t \leq t_{max}^f\)、
\(g(t)\)の定義域を\(t_{min}^g \leq t \leq t_{max}^g\)とする。
このとき\(f(t)\)と\(g(t)\)の間の相互相関関数には2種類の定義が考えられる。
Let \(f(t)\) and \(g(t)\) be the two time series data
used for the calculation of a cross correlation function (CCF).
Assume that \(f(t)\) and \(g(t)\) are defined in
\(t_{min}^f \leq t \leq t_{max}^f\) and
\(t_{min}^g \leq t \leq t_{max}^g\), respectively.
Then CCF between \(f(t)\) and \(g(t)\) can be defined
in two different ways.
1つ目の定義では\(f(t)\), \(g(t)\)が定義域外で値0を取るものと見なして
\(-\infty<t<\infty\)の時系列データに拡張する。
この場合、相互相関関数\(C_1(f,g;\tau)\)は以下のように定義される。
\[\begin{equation}
C_1(f,g;\tau) \equiv
\frac{\int_{-\infty}^{\infty} f(t)g(t+\tau)dt}
{\sqrt{\int_{-\infty}^{\infty}
f(t)^2dt \int_{-\infty}^{\infty} g(t+\tau)^2dt}}
=\frac{\int_{t_{min}^c}^{t_{max}^c} f(t)g(t+\tau)dt}
{\sqrt{\int_{t_{min}^f}^{t_{max}^f} f(t)^2dt
\int_{t_{min}^g}^{t_{max}^g} g(t)^2dt}}
\label{eq.C1}
\end{equation}\]
In the first definition, \(f(t)\) and \(g(t)\) are extended to
time series data in \(-\infty<t<\infty\),
assuming that thay are zero outside of their original definition ranges.
The CCF under this assumption, denoted as \(C_1(f,g;\tau)\),
is defined by eq. (\ref{eq.C1}).
2つ目の定義では
全ての被積分関数が定義される範囲のみで積分を行う。
この場合、相互相関関数\(C_2(f,g;\tau)\)は以下のように定義される。
\[\begin{equation}
C_2(f,g;\tau) \equiv
\frac{\int_{t_{min}^c}^{t_{max}^c} f(t)g(t+\tau)dt}
{\sqrt{\int_{t_{min}^c}^{t_{max}^c} f(t)^2dt
\int_{t_{min}^c}^{t_{max}^c} g(t+\tau)^2dt}}
\label{eq.C2}
\end{equation}\]
In the second definition,
integrals are conducted only in the range where
all the integrands are defined.
The CCF for this case, denoted as \(C_2(f,g;\tau)\),
is thus defined by eq. (\ref{eq.C2}).
ここで\([t_{min}^c,t_{max}^c]\)は
\(f(t)\)と\(g(t+\tau)\)の定義域の共通部分を表す。
\(f(t)\)の定義域は\(t_{min}^f \leq t \leq t_{max}^f\)である。
\(g(t+\tau)\)の定義域は\(t_{min}^g \leq t+\tau \leq t_{max}^g\)
すなわち\(t_{min}^g-\tau \leq t \leq t_{max}^g-\tau\)である。
したがって
\[\begin{equation}
t_{min}^c=\max\{t_{min}^f,t_{min}^g-\tau\}
\label{eq.tminc}
\end{equation}\]
\[\begin{equation}
t_{max}^c=\max\{t_{min}^c,\min\{t_{max}^f,t_{max}^g-\tau\}\}
\label{eq.tmaxc}
\end{equation}\]
であることが分かる。
Here, \([t_{min}^c,t_{max}^c]\) indicates the range of \(t\)
in which both \(f(t)\) and \(g(t+\tau)\) are defined.
The function \(f(t)\) is defined in \(t_{min}^f \leq t \leq t_{max}^f\),
while the function \(g(t+\tau)\) is defined
in \(t_{min}^g \leq t+\tau \leq t_{max}^g\),
i.e., \(t_{min}^g-\tau \leq t \leq t_{max}^g-\tau\).
Therefore \(t_{min}^c\) and \(t_{max}^c\) are given by
eqs. (\ref{eq.tminc}) and (\ref{eq.tmaxc}).
2. 相互相関関数の定義域
(The definition range of CCF)
\(|\tau|\)が非常に大きくなると(\ref{eq.C1})(\ref{eq.C2})式の分子の
\(f(t)\)と\(g(t+\tau)\)の定義域が共通部分を持たなくなる。
このときデータの値によらず\(C_1=0\)となるので\(C_1\)を計算する意味が無く、
\(C_2\)は0/0の不定値で定義不能となる。
このことから相互相関関数は
\(f(t)\)と\(g(t+\tau)\)の定義域が共通部分を持つ\(\tau\)の範囲でのみ
計算すれば良いことが分かる。
When \(|\tau|\) is very large,
the definition ranges of \(f(t)\) and \(g(t+\tau)\),
which are in the numerators of eqs. (\ref{eq.C1}) and (\ref{eq.C2}),
do not overlap.
In this case, \(C_1\) needs not be calculated
as it always becomes zero regardless of the data values,
and \(C_2\) cannot be calculated because it is 0/0.
Therefore it is enough to calculate the CCF
only for the range of \(\tau\) in which
the definition ranges of \(f(t)\) and \(g(t+\tau)\) overlap.
この条件は\(t_{min}^c<t_{max}^c\)と表され、
(\ref{eq.tminc})(\ref{eq.tmaxc})式を代入すると
\[\begin{equation}
\max\{t_{min}^f,t_{min}^g-\tau\} < \min\{t_{max}^f,t_{max}^g-\tau\}
\label{eq.tau.range.original}
\end{equation}\]
となる。(\ref{eq.tau.range.original})式は
\[\begin{equation}
t_{min}^f < t_{max}^f
\label{eq.tau.range.1}
\end{equation}\]
\[\begin{equation}
t_{min}^f < t_{max}^g-\tau
\label{eq.tau.range.2}
\end{equation}\]
\[\begin{equation}
t_{min}^g-\tau < t_{max}^f
\label{eq.tau.range.3}
\end{equation}\]
\[\begin{equation}
t_{min}^g-\tau < t_{max}^g-\tau
\label{eq.tau.range.4}
\end{equation}\]
が全て満たされることと等価であるが、
(\ref{eq.tau.range.1})(\ref{eq.tau.range.4})は自明であるので
(\ref{eq.tau.range.2})(\ref{eq.tau.range.3})を連立させて
\[\begin{equation}
t_{min}^g-t_{max}^f < \tau < t_{max}^g-t_{min}^f
\label{eq.tau.range}
\end{equation}\]
を得る。すなわち相互相関関数は(\ref{eq.tau.range})式の\(\tau\)の範囲で
計算すれば良い。
This requirement is expressed as \(t_{min}^c<t_{max}^c\).
Inserting eqs. (\ref{eq.tminc}) and (\ref{eq.tmaxc}) into this requirement
results in (\ref{eq.tau.range.original}).
This condition is equivalent to requiring all of
eqs. (\ref{eq.tau.range.1})-(\ref{eq.tau.range.4}).
However, eqs. (\ref{eq.tau.range.1}) and (\ref{eq.tau.range.4}) are trivial.
Using the remaining equations
(\ref{eq.tau.range.2}) and (\ref{eq.tau.range.3}),
we obtain (\ref{eq.tau.range}).
Therefore the CCF needs to be calculated
for the range of \(\tau\) given by eq. (\ref{eq.tau.range}).
3. 時間領域における\(C_1\)の離散化
(Discretization of \(C_1\) in the time domain)
時系列データ\(f(t)\)が時刻
\[\begin{equation}
t=t_0^f + k\Delta t \hspace{1em} (k=0,1,2,\cdots,K-1)
\label{eq.t.f}
\end{equation}\]
で、\(g(t)\)が時刻
\[\begin{equation}
t=t_0^g + l\Delta t \hspace{1em} (l=0,1,2,\cdots,L-1)
\label{eq.t.g}
\end{equation}\]
で与えられているものとし
\[\begin{equation}
\tau=\tau_0+m\Delta t \hspace{1em} (m=0,1,2,\cdots,M-1)
\label{eq.tau}
\end{equation}\]
における相互相関関数の値を求めることを考える。
以下では
\[\begin{equation}
f_k \equiv f(t_0^f + k\Delta t)
\label{eq.fk}
\end{equation}\]
\[\begin{equation}
g_l \equiv f(t_0^g + l\Delta t)
\label{eq.gl}
\end{equation}\]
\[\begin{equation}
C_1^m \equiv C_1(f,g;\tau_0 + m\Delta t)
\label{eq.C1m}
\end{equation}\]
という表記を用いる。
Let the time series data \(f(t)\) and \(g(t)\) be defined
at times given by eqs. (\ref{eq.t.f}) and (\ref{eq.t.g}), respectively.
Consider to calculate the CCF between \(f(t)\) and \(g(t)\)
for \(\tau\) given at (\ref{eq.tau}).
Below, we use the notation given by
eqs. (\ref{eq.fk})-(\ref{eq.C1m}).
このとき
\[\begin{equation}
t_{min}^f=t_0^f
\label{eq.tminf.discrete}
\end{equation}\]
\[\begin{equation}
t_{max}^f=t_0^f+(K-1)\Delta t
\label{eq.tmaxf.discrete}
\end{equation}\]
\[\begin{equation}
t_{min}^g=t_0^g
\label{eq.tming.discrete}
\end{equation}\]
\[\begin{equation}
t_{max}^g=t_0^g+(L-1)\Delta t
\label{eq.tmaxg.discrete}
\end{equation}\]
であり、これらを(\ref{eq.tau.range})式に代入すると\(\tau\)の最小値は
\[\begin{equation}
\tau_0
=t_{min}^g-t_{max}^f
=t_0^g-t_0^f-(K-1)\Delta t
\label{eq.tau_min.discrete}
\end{equation}\]
最大値は
\[\begin{equation}
\tau_0+(M-1)\Delta t
=t_{max}^g-t_{min}^f
=t_0^g-t_0^f+(L-1)\Delta t
\label{eq.tau_max.discrete}
\end{equation}\]
となる。なお(\ref{eq.tau.range})式は等号を含まない不等式になっているが、
連続変数としての\(\tau\)に対する不等式であるので
1サンプル(\(\Delta t\))だけ範囲を縮めることはできない。
そのため両端を含むように\(\tau\)の範囲を取っている。
(\ref{eq.tau_max.discrete})\(-\)(\ref{eq.tau_min.discrete})より
\[\begin{equation}
(M-1)\Delta t=(L-1)\Delta t+(K-1)\Delta t
\label{eq.M.derive}
\end{equation}\]
となり、ここから
\[\begin{equation}
M=K+L-1
\label{eq.M}
\end{equation}\]
が得られる。
Then eqs. (\ref{eq.tminf.discrete})-(\ref{eq.tmaxg.discrete}) hold,
and inserting these equations into (\ref{eq.tau.range}) results in
eqs. (\ref{eq.tau_min.discrete}) and (\ref{eq.tau_max.discrete})
for the minimum and maximum values, respectively, of \(\tau\).
Here, although eq. (\ref{eq.tau.range}) does not include the equality,
it is a condition for a continuous variable \(\tau\),
meaning that the range of \(\tau\) cannot be shortened
by one sample \(\Delta t\).
We therefore used the range of \(\tau\) for the discrete case
to include the both ends.
Subtracting eq. (\ref{eq.tau_max.discrete}) by (\ref{eq.tau_min.discrete})
results in (\ref{eq.M.derive}),
which can be arranged as (\ref{eq.M}).
以上の準備のもとで(\ref{eq.C1})式の積分を離散化して計算する。
分母の離散化は容易である。台形公式を用いると
\[\begin{eqnarray}
\int_{-\infty}^{\infty} f(t)^2 dt
&\sim& \sum_{k=-\infty}^{\infty}
\frac{f(t_0^f+k\Delta t)^2+f(t_0^f+(k+1)\Delta t)^2}{2}\Delta t
\nonumber \\
&=& \Delta t \sum_{k=-\infty}^{\infty} f(t_0^f+k\Delta t)^2 \nonumber \\
&=& \Delta t \sum_{k=0}^{K-1} f(t_0^f+k\Delta t)^2 \nonumber \\
&=& \Delta t \sum_{k=0}^{K-1} f_k^2
\label{eq.C1.f2}
\end{eqnarray}\]
となり、同様に
\[\begin{eqnarray}
\int_{-\infty}^{\infty} g(t+\tau)^2 dt
&=& \int_{-\infty}^{\infty} g(t)^2 dt \nonumber \\
&\sim& \Delta t \sum_{l=0}^{L-1} g_l^2
\label{eq.C1.g2}
\end{eqnarray}\]
も成り立つ。
Based on these preparations, now the integrals in eq. (\ref{eq.C1})
are calculated in discretized forms.
The calculation of the denominator is easy.
Using the trapezoid approximations,
eqs. (\ref{eq.C1.f2}) and (\ref{eq.C1.g2}) are obtained.
(\ref{eq.C1})式の分子を
\[\begin{equation}
h(\tau) \equiv \int_{-\infty}^{\infty} f(t)g(t+\tau)dt
\label{eq.h}
\end{equation}\]
として
\[\begin{equation}
h_m \equiv h(\tau_0+m\Delta t)
\label{eq.hm}
\end{equation}\]
とおく。\(h_m\)は台形公式を用いて
\[\begin{eqnarray}
h_m
&=& \int_{-\infty}^{\infty} f(t)g(t+\tau_0+m\Delta t)dt \nonumber \\
&\sim& \sum_{k=-\infty}^{\infty}
f(t_0^f+k\Delta t)g(t_0^f+k\Delta t+\tau_0+m\Delta t)\Delta t
\nonumber \\
&=& \Delta t \sum_{k=-\infty}^{\infty}
f(t_0^f+k\Delta t)g(t_0^g-(K-1)\Delta t+k\Delta t+m\Delta t)
\nonumber \\
&=& \Delta t \sum_{k=-\infty}^{\infty} f_k g_{k+m-K+1}
\label{eq.hm.sum}
\end{eqnarray}\]
と書けるが、この和の各項は
\(f_k\neq 0\)かつ\(g_{k+m-K+1}\neq 0\)の場合にのみノンゼロの値を持つ。
したがって和を取る範囲としては
- \(f_k\neq 0\)となる条件から\(0\leq k\leq K-1\)
- \(g_{k+m-K+1}\neq 0\)となる条件から\(0\leq k+m-K+1\leq L-1\)、
すなわち\(K-m-1\leq k\leq K+L-m-2\)
とすれば良い。
この両方の条件が成り立つ\(k\)についてのみ和を取れば良いのであるから
\[\begin{equation}
h_m=\Delta t \sum_{k=\max\{0,K-m-1\}}^{\min\{K-1,K+L-m-2\}} f_k g_{k+m-K+1}
\label{eq.hm.sum.tdomain}
\end{equation}\]
が得られる。
Next, the numerator of eq. (\ref{eq.C1}) is represented by
eq. (\ref{eq.h}), and a notation of eq. (\ref{eq.hm}\) is introduced.
Using the trapezoid approximation, \(h_m\) is calculated
as (\ref{eq.hm.sum}).
Each term of this summation has a non-zero value
only when \(f_k\neq 0\) and \(g_{k+m-K+1}\neq 0\).
Therefore the summation range of \(k\) can be:
- \(0\leq k\leq K-1\) based on the condition of \(f_k\neq 0\), and
- \(0\leq k+m-K+1\leq L-1\), i.e., \(K-m-1 \leq k \leq K+L-m-2\),
based on the condition of \(g_{k+m-K+1}\neq 0\).
The summation needs to be taken only for \(k\) that meets
both of these conditions,
giving eq. (\ref{eq.hm.sum.tdomain}).
(\ref{eq.C1.f2})(\ref{eq.C1.g2})(\ref{eq.hm.sum.tdomain})を
(\ref{eq.C1})に代入すれば
\[\begin{equation}
C_1^m=\frac{\sum_{k=\max\{0,K-m-1\}}^{\min\{K-1,K+L-m-2\}} f_k g_{k+m-K+1}}
{\sqrt{\sum_{k=0}^{K-1} f_k^2 \cdot \sum_{l=0}^{L-1} g_l^2}}
\label{eq.C1m.discrete}
\end{equation}\]
が得られる。
Inserting eqs. (\ref{eq.C1.f2}), (\ref{eq.C1.g2}),
and (\ref{eq.hm.sum.tdomain}) into (\ref{eq.C1})
results in (\ref{eq.C1m.discrete}).
4. 周波数領域における\(C_1\)の離散化
(Discretization of \(C_1\) in the frequency domain)
\(h(\tau)\)の計算を周波数領域で行うことを考える。
\(f(t)\), \(g(t)\), \(h(\tau)\)のフーリエスペクトルを
それぞれ\(F(\omega)\), \(G(\omega)\), \(H(\omega)\)とするとき、
\(H(\omega)=F(\omega)^{∗} G(\omega)\)の関係が成り立つ。
ここで\(^{∗}\)は複素共役を表す。
この関係から、\(f(t)\), \(g(t)\)をフーリエ変換し、
そのクロススペクトルを計算することで\(H(\omega)\)を求め、
それをフーリエ逆変換すれば\(h(\tau)\)が求まると考えられる。
Let us now consider to calculate \(h(\tau)\) in the frequency domain.
Let \(F(\omega)\), \(G(\omega)\), and \(H(\omega)\) be
the Fourier spectra of \(f(t)\), \(g(t)\), and \(h(\tau)\), respectively.
Then a relation \(H(\omega)=F(\omega)^{∗} G(\omega)\) holds,
where \(^{∗}\) denotes the complex conjugate.
Based on this relation, \(h(\tau)\) is expected to be obtained by
calculating the Fourier spectra of \(f(t)\) and \(g(t)\),
then computing the cross spectrum of them
which is equal to \(H(\omega)\),
and finally applying the Fourier inverse transformation to it.
有限の長さの離散的な時系列データに対してこの考え方を適用するためには
\(f(t)\), \(g(t)\)の末尾にダミーの0を埋めて同じ長さにしなければならない。
さもないとフーリエスペクトルの周波数刻みが
\(F(\omega)\)と\(G(\omega)\)で異なってしまい、
クロススペクトルを計算できないからである。
以下ではこのダミーの0を埋めて作成する時系列データの
適切な長さについて検討する。
To apply this approach to discrete time series data of finite lengths,
dammy zeroes must be padded after \(f(t)\) and \(g(t)\)
to make them equal lengths.
By this, the frequency intervals of \(F(\omega)\) and \(G(\omega)\)
become equal, which enables the calculation of the cross spectrum.
Below, let us consider the appropriate length of the time series data
to create by padding dammy zeroes.
まずは\(h_m\)の離散フーリエ変換\(H_n\)の具体的な表記を求める。
離散フーリエ変換の定義式に
(\ref{eq.hm.sum})(\ref{eq.tau_min.discrete})式を代入すると
\[\begin{eqnarray}
H_n
&=& \exp\frac{2\pi i n \tau_0}{M \Delta t}\Delta t
\sum_{m=0}^{M-1} h_m\exp\frac{2\pi i n m}{M}
\nonumber \\
&=& \exp\frac{2\pi i n [t_0^g-t_0^f-(K-1)\Delta t]}{M \Delta t}
(\Delta t)^2
\nonumber \\
& & \sum_{m=0}^{M-1}\sum_{k=-\infty}^{\infty}
f_k g_{k+m-K+1}\exp\frac{2\pi i n m}{M}
\label{eq.Hn.derive1}
\end{eqnarray}\]
となる。
(\ref{eq.Hn.derive1})式において
\(k<0\), \(k\geq K\)の項は\(f_k=0\)なので和に対して寄与しない。
したがって\(k\)に関する和の範囲は\(k=0\)から\(K-1\)までとすれば十分である。
このことを用いて更に変形を進めると
\[\begin{eqnarray}
H_n
&=& \exp\frac{2\pi i n [t_0^g-t_0^f-(K-1)\Delta t]}{M \Delta t}
(\Delta t)^2
\nonumber \\
& & \sum_{m=0}^{M-1}\sum_{k=0}^{K-1}
f_k g_{k+m-K+1}\exp\frac{2\pi i n m}{M}
\nonumber \\
&=& \exp\frac{2\pi i n t_0^g}{M \Delta t}
\exp\left(-\frac{2\pi i n t_0^f}{M \Delta t}\right)
\exp\left[-\frac{2\pi i n (K-1)}{M}\right]
(\Delta t)^2
\nonumber \\
& & \sum_{k=0}^{K-1}\sum_{m=0}^{M-1}
f_k g_{k+m-K+1}\exp\frac{2\pi i n m}{M}
\nonumber \\
&=& \exp\frac{2\pi i n t_0^g}{M \Delta t}
\exp\left(-\frac{2\pi i n t_0^f}{M \Delta t}\right)
(\Delta t)^2
\nonumber \\
& & \sum_{k=0}^{K-1}\sum_{m=0}^{M-1}
f_k g_{k+m-K+1}\exp\frac{2\pi i n (m-K+1)}{M}
\nonumber \\
&=& \exp\frac{2\pi i n t_0^g}{M \Delta t}
\exp\left(-\frac{2\pi i n t_0^f}{M \Delta t}\right)
(\Delta t)^2
\nonumber \\
& & \sum_{k=0}^{K-1}\sum_{m=0}^{M-1}
f_k g_{k+m-K+1}
\exp\left(-\frac{2\pi ink}{M}\right)
\exp\frac{2\pi i n (k+m-K+1)}{M}
\nonumber \\
&=& \exp\frac{2\pi i n t_0^g}{M \Delta t}
\exp\left(-\frac{2\pi i n t_0^f}{M \Delta t}\right)
(\Delta t)^2
\nonumber \\
& & \sum_{k=0}^{K-1}\sum_{l=k-K+1}^{k+M-K}
f_k g_l
\exp\left(-\frac{2\pi ink}{M}\right)
\exp\frac{2\pi inl}{M}
\label{eq.Hn.derive2}
\end{eqnarray}\]
となる。
(\ref{eq.Hn.derive2})式において\(l<0\), \(l\geq L\)の項は
\(g_l=0\)なので和に対して寄与しない。
したがって\(l\)に関する和の範囲を狭めることができて
\[\begin{eqnarray}
H_n
&=& \exp\frac{2\pi i n t_0^g}{M \Delta t}
\exp\left(-\frac{2\pi i n t_0^f}{M \Delta t}\right)
(\Delta t)^2
\nonumber \\
& & \sum_{k=0}^{K-1}\sum_{l=\max\{k-K+1,0\}}^{\min\{k+M-K,L-1\}}
f_k g_l
\exp\left(-\frac{2\pi ink}{M}\right)
\exp\frac{2\pi inl}{M}
\label{eq.Hn.derive3}
\end{eqnarray}\]
と書き直せる。ところで\(k\leq K-1\)より\(k-K+1\leq 0\)である。
また(\ref{eq.M})式から\(k+M-K=k+L-1\geq L-1\)である。
したがって\(l\)に関する和の範囲は結局0から\(L-1\)までとなり、
\[\begin{eqnarray}
H_n
&=& \exp\frac{2\pi i n t_0^g}{M \Delta t}
\exp\left(-\frac{2\pi i n t_0^f}{M \Delta t}\right)
(\Delta t)^2
\nonumber \\
& & \sum_{k=0}^{K-1}\sum_{l=0}^{L-1}
f_k g_l
\exp\left(-\frac{2\pi ink}{M}\right)
\exp\frac{2\pi inl}{M}
\nonumber \\
&=& \exp\left(-\frac{2\pi i n t_0^f}{M \Delta t}\right)\Delta t
\sum_{k=0}^{K-1}f_k\exp\left(-\frac{2\pi ink}{M}\right)
\nonumber \\
& & \exp\frac{2\pi i n t_0^g}{M \Delta t}\Delta t
\sum_{l=0}^{L-1}g_l\exp\frac{2\pi inl}{M}
\label{eq.Hn}
\end{eqnarray}\]
という式が得られる。
First, let us derive an expression for
the discrete Fourier transformation, \(H_n\), of \(h_m\).
Inserting eqs. (\ref{eq.hm.sum}) and (\ref{eq.tau_min.discrete})
into the definition of the discrete Fourier transformation
results in eq. (\ref{eq.Hn.derive1}).
In this equation, terms for \(k<0\) and \(k\geq K\)
do not contribute to the summation because \(f_k=0\).
Therefore it is enough to take the summation of \(k\)
from \(k=0\) to \(K-1\).
Using it, the equation can further be arranged as (\ref{eq.Hn.derive2}).
In eq. (\ref{eq.Hn.derive2}), terms for \(l<0\) and \(l\geq L\)
do not contribute to the summation because \(g_l=0\).
Therefore the summation range for \(l\) can be narrowed,
giving eq. (\ref{eq.Hn.derive3}).
Here, \(k-K+1\leq 0\) because \(k\leq K-1\).
Also from eq. (\ref{eq.M}), it is derived that
\(k+M-K=k+L-1\geq L-1\).
Therefore the summation range for \(l\) is from 0 to \(L-1\),
finally yielding eq. (\ref{eq.Hn}).
(\ref{eq.Hn})式の右辺は
\(f_k\)の末尾にダミーの0を追加して
データ点数を\(M\)にした時系列データの離散フーリエ変換と、
\(g_l\)の末尾にダミーの0を追加して
データ点数を\(M\)にした時系列データの離散フーリエ変換との
クロススペクトルになっている。
このことから、周波数領域で相互相関関数を計算するためには
データ点数が\(M\)になるように
\(f_k\), \(g_l\)の末尾にダミーの0を追加すれば良いことが分かる。
The right hand side of eq. (\ref{eq.Hn}) is
the cross spectrum between
the discrete Fourier transformations of the time series data
obtained by adding dammy zeroes at the ends of \(f_k\) and \(g_l\)
to make their lengths \(M\).
This indicates that the CCF can be computed in the frequency domain
by using data of \(M\) samples,
created by adding dammy zeroes at the ends of \(f_k\) and \(g_l\).
5. 時間領域における\(C_2\)の離散化
(Discretization of \(C_2\) in the time domain)
\(C_1\)の計算式の導出にあたっては
\(f(t)\), \(g(t)\)が定義域外で値0を取るという仮定を繰り返し用いた。
例えば(\ref{eq.hm.sum})式では
無限区間の積分を\(k=-\infty\)から\(k=\infty\)までの和で近似し、
その後でノンゼロの値を持つ項のみを残すことで
有限項の和(\ref{eq.hm.sum.tdomain})を得た。
一方、\(C_2\)については最初の定義式自体が有限区間の積分になっている
(\ref{eq.C2}式)。
そのため離散化のアプローチも異なってくる。
In the derivation of the formula for \(C_1\),
an assumption, that \(f(t)\) and \(g(t)\) take the value of 0
outside of the definition range, was repeatedly used.
For example, eq. (\ref{eq.hm.sum}) is an approximation of
an integral in an infinite range by a summation
from \(k=-\infty\) to \(k=\infty\),
which is then modified to a summation of a finite number of terms
(eq. \ref{eq.hm.sum.tdomain}) by remaining only the terms
that have non-zero values.
In case of \(C_2\), however, the original definition (eq. \ref{eq.C2})
is already an integral in a finite range,
so that the discretization procedure is different from the case of \(C_1\).
\(C_2\)の計算ではまず積分範囲を表す\(t_{min}^c\), \(t_{max}^c\)の
離散化表記を求める。
(\ref{eq.tminc})に
(\ref{eq.tminf.discrete})
(\ref{eq.tming.discrete})
(\ref{eq.tau})
(\ref{eq.tau_min.discrete})
を代入すると
\[\begin{eqnarray}
t_{min}^c
&=& \max\{t_{min}^f,t_{min}^g-\tau\}
\nonumber \\
&=& \max\{t_0^f,t_0^g-(\tau_0+m\Delta t)\}
\nonumber \\
&=& \max\{t_0^f,t_0^g-\tau_0-m\Delta t\}
\nonumber \\
&=& \max\{t_0^f,t_0^g-[t_0^g-t_0^f-(K-1)\Delta t]-m\Delta t\}
\nonumber \\
&=& \max\{t_0^f,t_0^f+(K-m-1)\Delta t\}
\label{eq.tminc.discrete}
\end{eqnarray}\]
が得られる。
また(\ref{eq.tmaxc})に
(\ref{eq.tmaxf.discrete})
(\ref{eq.tmaxg.discrete})
(\ref{eq.tau})
(\ref{eq.tau_min.discrete})
を代入すると
\[\begin{eqnarray}
t_{max}^c
&=& \max\{t_{min}^c,\min\{t_{max}^f,t_{max}^g-\tau\}\}
\nonumber \\
&=& \max\{t_{min}^c,
\min\{t_0^f+(K-1)\Delta t,
t_0^g+(L-1)\Delta t-(\tau_0+m\Delta t)\}\}
\nonumber \\
&=& \max\{t_{min}^c,
\min\{t_0^f+(K-1)\Delta t,t_0^g-\tau_0+(L-m-1)\Delta t\}\}
\nonumber \\
&=& \max\left\{t_{min}^c,
\min\left\{t_0^f+(K-1)\Delta t,
\right.\right. \nonumber \\
& & \left.\left.
t_0^g-[t_0^g-t_0^f-(K-1)\Delta t]+(L-m-1)\Delta t
\right\}\right\}
\nonumber \\
&=& \max\{t_{min}^c,
\min\{t_0^f+(K-1)\Delta t,t_0^f+(K+L-m-2)\Delta t\}\}
\label{eq.tmaxc.discrete}
\end{eqnarray}\]
が得られる。
したがって積分を和で置き換える際の\(k\)の範囲は
\[\begin{equation}
k_{min} \equiv \max\{0,K-m-1\}
\label{eq.k_min}
\end{equation}\]
から
\[\begin{equation}
k_{max} \equiv \max\{k_{min},\min\{K-1,K+L-m-2\}\}
\label{eq.k_max}
\end{equation}\]
までとすれば良いことが分かる。
First, let us derive the discretized representation
for the integral range \([t_{min}^c,t_{max}^c]\).
Inserting eqs. (\ref{eq.tminf.discrete}), (\ref{eq.tming.discrete}),
(\ref{eq.tau}), and (\ref{eq.tau_min.discrete})
into (\ref{eq.tminc}) results in (\ref{eq.tminc.discrete}).
Also, inserting eqs. (\ref{eq.tmaxf.discrete}), (\ref{eq.tmaxg.discrete}),
(\ref{eq.tau}), and (\ref{eq.tau_min.discrete})
into (\ref{eq.tmaxc}) results in (\ref{eq.tmaxc.discrete}).
Therefore the range of \(k\) used to replace
an integral by a summation should be from \(k_{min}\) to \(k_{max}\),
which are defined by eqs. (\ref{eq.k_min}) and (\ref{eq.k_max}),
respectively.
この\(k_{min}\), \(k_{max}\)を用いれば
(\ref{eq.C2})式の右辺に登場する積分は
以下のように台形公式近似することができる。
\[\begin{eqnarray}
\int_{t_{min}^c}^{t_{max}^c}f(t)g(t+\tau)dt
&=& \int_{t_{min}^c}^{t_{max}^c} f(t)g(t+\tau_0+m\Delta t)dt
\nonumber \\
&=& \int_{t_{min}^c}^{t_{max}^c} f(t)g(t+t_0^g-t_0^f-(K-1)\Delta t+m\Delta t)dt
\nonumber \\
&\sim& \sum_{k=k_{min}}^{k_{max}}
f(t_0^f+k\Delta t)
\nonumber \\
& & g(t_0^f+k\Delta t+t_0^g-t_0^f-(K-1)\Delta t+m\Delta t)
\Delta t \nonumber \\
&=& \Delta t \sum_{k=k_{min}}^{k_{max}}
f(t_0^f+k\Delta t)g(t_0^g+(k+m-K+1)\Delta t)
\nonumber \\
&=& \Delta t \sum_{k=k_{min}}^{k_{max}}
f_k g_{k+m-K+1}
\label{eq.int_fg_C2.discrete}
\end{eqnarray}\]
同様にして
\[\begin{equation}
\int_{t_{min}^c}^{t_{max}^c} f(t)^2 dt
\sim \Delta t \sum_{k=k_{min}}^{k_{max}} f_k^2
\label{eq.int_f2_C2.discrete}
\end{equation}\]
\[\begin{equation}
\int_{t_{min}^c}^{t_{max}^c} g(t+\tau)^2 dt
\sim \Delta t \sum_{k=k_{min}}^{k_{max}} g_{k+m-K+1}^2
\label{eq.int_g2_C2.discrete}
\end{equation}\]
ここで(\ref{eq.tau})(\ref{eq.tau_min.discrete})を用いた。
これらを(\ref{eq.C2})に代入すると
\[\begin{equation}
C_2^m\equiv C_2(f,g;\tau_0+m\Delta t)
\sim\frac{\sum_{k=k_{min}}^{k_{max}} f_k g_{k+m-K+1}}
{\sqrt{\left(\sum_{k=k_{min}}^{k_{max}} f_k^2\right)
\left(\sum_{k=k_{min}}^{k_{max}} g_{k+m-K+1}^2\right)}}
\label{eq.C2m.discrete}
\end{equation}\]
が得られる。
Using these \(k_{min}\) and \(k_{max}\),
the integrals in the right hand side of eq. (\ref{eq.C2})
can be calculated, in the forms of trapezoid approximations, as
eqs. (\ref{eq.int_fg_C2.discrete})-(\ref{eq.int_g2_C2.discrete}),
where eqs. (\ref{eq.tau}) and (\ref{eq.tau_min.discrete}) were used.
Inserting these relations into (\ref{eq.C2})
results in (\ref{eq.C2m.discrete}).
6. \(C_1\)と\(C_2\)のどちらを用いるべきか
(Which of \(C_1\) or \(C_2\) be used?)
6-1. 短い時系列データ同士の相互相関関数を計算する場合
(In cases where the CCF between short time series data is calculated)
\(f(t)\), \(g(t)\)の長さが短く、
相互相関関数を最大のラグタイムまで計算する場合には
\(C_1\)を用いるべきである。
最大のラグタイムとは\(m=0\)または\(m=M-1=K+L-2\)のことであり、
\(f(t)\)と\(g(t+\tau)\)の定義域が
ちょうど1サンプルだけ重なる\(\tau\)に対応する。
この場合、\(C_2\)は1サンプルの相互相関なので\(\pm 1\)になってしまう。
\(C_1\)の方は\(f(t)\), \(g(t+\tau)\)の定義域外の値が0の部分を計算に含めるので
\(\pm 1\)にはならない。
\(C_1\) should be used in cases where
the lengths of \(f(t)\) and \(g(t)\) are short
and the CCF is computed over the maximum lag time.
Here the maximum lag time means \(m=0\) and \(m=M-1=K+L-2\),
which correspond to the cases where
only one sample overlaps between the definition ranges
of \(f(t)\) and \(g(t+\tau)\).
In this case, \(C_2=\pm 1\) because it is a cross correlation of one sample,
while \(C_1\) is not \(\pm 1\) because
\(f(t)\) and \(g(t+\tau)\) outside of their definition ranges
(with zero values) are used in the calculation.
例えば\(m=0\)のとき、\(C_1\), \(C_2\)は以下のようになる。
まず\(C_1\)については(\ref{eq.C1m.discrete})式に\(m=0\)を代入すると
\[\begin{eqnarray}
C_1^0
&=& \frac{\sum_{k=\max\{0,K-1\}}^{\min\{K-1,K+L-2\}} f_k g_{k-K+1}}
{\sqrt{\sum_{k=0}^{K-1} f_k^2 \cdot \sum_{l=0}^{L-1} g_l^2}}
\nonumber \\
&=& \frac{\sum_{k=K-1}^{K-1} f_k g_{k-K+1}}
{\sqrt{\sum_{k=0}^{K-1} f_k^2 \cdot \sum_{l=0}^{L-1} g_l^2}}
\nonumber \\
&=& \frac{f_{K-1}g_0}
{\sqrt{\sum_{k=0}^{K-1} f_k^2 \cdot \sum_{l=0}^{L-1} g_l^2}}
\label{eq.C1m.m0}
\end{eqnarray}\]
となり、\(\pm 1\)にはならない、一方、\(C_2\)の計算ではまず
(\ref{eq.k_min})(\ref{eq.k_max})式に\(m=0\)を代入すると
\[\begin{eqnarray}
k_{min}
&=& \max\{0,K-1\}
\nonumber \\
&=& K-1
\label{eq.k_min.m0}
\end{eqnarray}\]
\[\begin{eqnarray}
k_{max}
&=& \max\{k_{min},\min\{K-1,K+L-2\}\}
\nonumber \\
&=& \max\{K-1,K-1\}
\nonumber \\
&=& K-1
\label{eq.k_max.m0}
\end{eqnarray}\]
となり、これを(\ref{eq.C2m.discrete})に代入すると
\[\begin{eqnarray}
C_2^0
&=& \frac{\sum_{k=K-1}^{K-1}f_kg_{k-K+1}}
{\sqrt{\left(\sum_{k=K-1}^{K-1}f_k^2\right)
\left(\sum_{k=K-1}^{K-1}g_{k-K+1}^2\right)}}
\nonumber \\
&=& \frac{f_{K-1}g_0}{\sqrt{f_{K-1}^2 g_0^2}}
\nonumber \\
&=& \pm 1
\label{eq.C2m.m0}
\end{eqnarray}\]
となる。\(m=K+L-2\)のときも同様にして
\(C_1\neq \pm 1\), \(C_2=\pm 1\)となる。
For example, \(C_1\) and \(C_2\) for \(m=0\) are calculated as below.
Inserting \(m=0\) into eq. (\ref{eq.C1m.discrete})
results in the \(C_1\) value given by eq. (\ref{eq.C1m.m0}),
which is not \(\pm 1\).
To calculate \(C_2\), \(m=0\) is first inserted into
eqs. (\ref{eq.k_min}) and (\ref{eq.k_max}),
giving eqs. (\ref{eq.k_min.m0}) and (\ref{eq.k_max.m0}).
Inserting these relations into (\ref{eq.C2m.discrete})
reduces to (\ref{eq.C2m.m0}).
The case for \(m=K+L-2\) can be calculated in the same way,
and the result is that \(C_1\neq\pm 1\) and \(C_2=\pm 1\).
6-2. 長さが大きく異なる時系列データ間の相互相関関数を求める場合
(In cases where the CCF between time series data
of significantly different lengths is calculated
\(f(t)\), \(g(t)\)の長さが大きく異なる場合には
\(C_2\)を用いるべきである。
例えばマスターイベントの波形と連続波形との相互相関関数を用いて
イベントを検出する場合(マッチドフィルター解析)がこのケースに該当する。
もしもこの解析に\(C_1\)を用いると、
マスターイベントの波形において定義域外では値は全て0であると仮定することになり、
連続波形の方は0ではないので相互相関の値が小さくなってしまう。
例えば30秒間のマスターイベント波形を用いて1日分のイベントを検知する場合、
\(C_1\)を用いると
「30秒間だけ振幅がノンゼロで残りの23時間59分30秒間は振幅がゼロの
マスターイベント波形」
と実際の24時間分の連続波形との相互相関関数を計算することになる。
ノンゼロの部分同士がどれほどよく合っていたとしても
「23時間59分30秒間にわたって振幅がゼロ」という特徴が実際の観測波形と合わないので
相関値が非常に小さくなり、
特にマスターイベント自身との相互相関も1にならない。
したがってマッチドフィルター解析では\(C_2\)を用いる必要がある。
\(C_2\) should be used in cases where
the lengths of \(f(t)\) and \(g(t)\) are significantly different.
An example for such cases is detectings seismic events based on CCFs
between the waveform of a master event and a continuous waveform
(i.e., a matched filter analysis).
A use of \(C_1\) for this purpose is equivalent to assume that
all the values outside of the definition range are zero
in the master event waveform,
whereas the continuous waveform is not zero,
resulting in small values of cross correlations.
For example, consider to detect events on one day
using a master event waveform of 30 s long.
In this case, \(C_1\) is a CCF between
“a master event waveform that has zero amplitudes over the 24 hours
except for the 30 s of the original definition range”
and an actual continuous waveform of 24 hours.
Even though the non-zero parts between the master and continuous waveforms
are very similar, the characteristics of
“zero amplitudes over almost the 24 hours”
in the master event waveform is significantly different from
the actual continuous waveform, giving very small correlation values.
Even the correlation with the master event itself is not 1.
For this reason, \(C_2\) should be used for a matched filter analysis.
7. 相互相関関数におけるラグタイムの意味
(Meaning of the lag time in CCF)
\(\tau\)は相互相関関数における\(f(t)\)と\(g(t)\)のラグタイムを表すが、
その意味については注意深く考えて用いる必要がある。
\(\tau \geq 0\)の場合が一見して\(g(t)\)を後にずらした場合のように見えるが、
実は\(g(t)\)を前にずらした場合に対応する。
これは、\(x\)軸を正方向に伝搬する波が\(f(t+x/V)\)ではなく\((t-x/V)\)で
表される事情に似ている。
この節では相関係数を用いる代表的な問題について\(\tau\)の意味を検討する。
The meaning of \(\tau\), which is the lag time
between \(f(t)\) and \(g(t)\) in the CCF,
is confusing and needs a careful consideration.
A case for \(\tau \geq 0\) may look
as if \(g(t)\) is shifted relatively later than \(f(t)\).
Indeed, \(\tau \geq 0\) means that
\(g(t)\) is shifted relatively earlier than \(f(t)\).
The situation is similar to the fact that
the wave propagating in positive direction of \(x\)-axis
is denoted by \((t-x/V)\), not \(f(t+x/V)\).
In this section, meanings of \(\tau\) are described
for several typical examples for the CCF analyses.
7-1. マスターイベントと連続波形の相互相関関数を用いてイベントを検出する場合
(Detecting events using the CCF between a master event
and a continuous waveform)
マスターイベントの波形を\(f(t)\)すなわち関数の第1引数に割り当てると
\(\tau\)がイベントの発震時刻になる。
マスターイベントの時刻をずらしながら連続波形との相関を逐次計算していくことを考えると
\(g(t)\)の方にマスターイベントを割り当てたくなるが、逆なのである。
The event origin times are represented by \(\tau\)
if the master waveform is denoted by \(f(t)\).
As the time origin for the master event is shifted,
\(g(t)\) may look more appropriate for the master waveform,
but this is not correct.
簡単な例として、マスターイベントがインパルス\(\delta(t)\)で表され、
連続波形中の時刻\(t=t_0\)にインパルスがある場合を考える。
簡単のため、マスターイベントも連続波形もこのインパルス以外には
一切のシグナルを含まない理想的な場合を考える。
\(\tau=t_0\)で相互相関関数が1になれば成功である。
マスターイベントを\(f(t)\), 連続波形を\(g(t)\)に割り当てた場合、
\(f(t)=\delta(t)\), \(g(t)=\delta(t-t_0)\)であるので、
\(g(t+\tau)=\delta(t-t_0+\tau)\)である。
したがって\(g(t+\tau)\)は時刻\(t=t_0-\tau\)にインパルスを持ち、
相互相関関数が1となるのはそのインパルスと
マスターイベントのインパルスの時刻が一致する場合、
すなわち\(t_0-\tau=0\)の場合である。
このことからマスターイベントを\(f(t)\), 連続波形を\(g(t)\)に割り当てれば
\(\tau=t_0\)で相互相関関数が1となってこの時刻でイベントが検出されることになる。
もし間違って連続波形を\(f(t)\), マスターイベントを\(g(t)\)に割り当ててしまうと
\(\tau=-t_0\)でイベントが検出されてしまう。
Let us consider a simple example where
the master event is represented by an impulse \(\delta(t)\)
and an impluse lies at a time \(t=t_0\) in the continuous waveform.
For simplicity, consider an ideal case in which
both the master and continuous waveforms
do not contain signals other than the impulses mensioned above.
The detection is regarded to be success
if the CCF takes unity at \(\tau=t_0\).
First, consider a case where
the master and continuous waveforms are represented
by \(f(t)\) and \(g(t)\), respectively.
In this case, \(f(t)=\delta(t)\) and \(g(t)=\delta(t-t_0)\).
This means \(g(t+\tau)=\delta(t-t_0+\tau)\),
which indicates that \(g(t+\tau)\) contains an impulse at a time \(t=t_0-\tau\).
The CCF takes unity if
this time of impulse for \(g(t+\tau)\) equals
the time of impulse for the master event, i.e., \(t_0-\tau=0\).
Thus the CCF takes unity and the event is detected at \(\tau=t_0\)
if the master and the continuous waveforms are represented
by \(f(t)\) and \(g(t)\), respectively.
In opposite, if the master and the continuous waveforms are represented
by \(g(t)\) and \(f(t)\), respectively,
the event is detected at \(\tau=-t_0\).
7-2. 観測点間の走時差を求める場合
(Calculation of a travel time difference between two stations)
連続微動のように初動走時を使った震源決定が困難なイベントでは
相互相関関数を用いた観測点間の走時差の決定が
震源決定のための手段として有効である。
For events whose first arrival is difficult to use,
such as a continuous tremor,
travel time differences among the stations
obtained by the CCF analysis
is an approach to determine the hypocenter.
ここでも簡単のためインパルス波形を考える。
観測点1での波形を\(f(t)=\delta(t)\)とし、
このインパルスが観測点2には時間差\(T\)だけ遅れて到達するものとする。
このとき観測点2での波形は\(g(t)=\delta(t-T)\)で表される。
このとき相互相関関数は
\[\begin{equation}
C(f,g;\tau)=\int f(t)g(t+\tau)dt
=\int \delta(t)\delta(t-T+\tau)dt
=\delta(\tau-T)
\label{eq.correlation.traveltime}
\end{equation}\]
と書ける。すなわちラグタイム\(\tau=T\)で相互相関関数の値が1になる。
Here, waveforms composed of an impulse is considered for simplicity.
Let the waveform at station 1 be represented as \(f(t)=\delta(t)\),
and assume that the impulse arrives at station 2 with a delay time \(T\).
Then, the waveform at station 2 is represented as \(g(t)=\delta(t-T)\).
In this case, the CCF can be written as eq. (\ref{eq.correlation.traveltime}),
indicating that the CCF becomes unity for a lag time \(\tau=T\).
この結果の意味するところは、観測点1の波形を\(f(t)\)に、観測点2の波形を\(g(t)\)に
代入して相互相関関数を計算してみて正の時刻差での相関が大きくなったら
観測点1の方に波が先に到達したということである。
ラグタイム\(\tau>0\)とすることは\(g(t)\)を前にずらすことに相当するが、
もともと走時差によって後にずれていたものをその分だけ前にずらして
両観測点の走時を一致させた時に相関が大きくなると理解すれば良い。
This simple example indicates that
if the waveforms at stations 1 and 2 are used
for the functions \(f(t)\) and \(g(t)\) respectively,
and if the resultant CCF takes the maximum at a positive lag time \(\tau\),
then it means that the wave arrived at station 1 earlier than station 2.
This is because the CCF takes a large value
for a \(\tau\) for which the travel time difference of the two waveforms
is cancelled by a time shift corresponding to \(\tau\);
a positive \(\tau\) corresponds to shifting \(g(t)\) in earlier direction,
by which a travel time delay which originally existed in \(g(t)\)
is cancelled.