橋梁に2つの加速度センサを取り付けて振動を測ったとします。片方は橋脚の根元、もう片方は橋の中央です。2つの波形を並べてみると、なんとなく似ている気もするし、似ていない気もする。ここで相関係数を計算すると、たとえば 0.42 という数字がひとつ出てきます。しかしこの 0.42 は、いったい何を教えてくれるのでしょうか。
実は、この単一の数字は多くの情報を潰してしまっています。橋の 2 Hz の固有振動については 2 点の動きがぴったり連動していて、20 Hz の高周波成分については風や交通のローカルな揺れが支配的でまったく連動していない、という状況を考えてください。相関係数はこの 2 つを混ぜて平均してしまうので、「まあまあ似ている」という曖昧な答えしか返しません。知りたいのは「似ているかどうか」ではなく、「どの周波数で似ているか」 なのです。
この問いに答えるのが コヒーレンス関数（coherence function） です。コヒーレンスは、2つの信号の線形な結びつきの強さを周波数ごとに 0 から 1 の数値で表します。ざっくり言えば「周波数ごとの相関係数の二乗」であり、機械振動の伝達関数測定（モーダル解析）では測定の信頼度を判定する必須の指標として、また脳波・脳磁図の解析では離れた脳部位が同じリズムで同期しているかを見る「機能的結合」の指標として、日常的に使われています。ほかにも、スピーカーとマイクの間のハウリング検出、レーダーの目標追尾における 2 チャネル間の相関判定、地震計アレイのノイズ評価など、応用範囲は非常に広い概念です。
そして、コヒーレンスには初学者が必ず一度は踏む有名な落とし穴があります。平均化をせずに 1 回の FFT だけでコヒーレンスを計算すると、どんな無関係な信号ペアでも答えが常に 1.0 になってしまう のです。この記事では、なぜそうなるのかを数式でごまかさずに示し、正しい計算手順と、コヒーレンスが 1 未満に落ちるときに何が起きているのかを徹底的に掘り下げます。
本記事の内容
- コヒーレンス関数の直感的な意味（周波数ごとの決定係数）
- クロススペクトル密度とマグニチュード二乗コヒーレンスの定義
- コーシー・シュワルツの不等式による $0 \le \gamma^2 \le 1$ の証明
- 単一セグメントでは $\hat{\gamma}^2 \equiv 1$ になることの証明と、平均化が必須である理由
- 平均化回数 $K$ に依存するバイアスと有意水準のしきい値
- 線形系 $y = h * x + n$ における $\gamma^2 = \mathrm{SNR}/(1+\mathrm{SNR})$ の導出
- コヒーレンスを下げる4つの原因（雑音・非線形・時間遅れ・未測定入力）
- 時間遅れによる低下が「窓関数の自己相関の二乗」で説明できることの導出
- Python（
scipy.signal.coherence）による検証と可視化
まず、この記事全体の見取り図を1枚の絵にしておきます。
左の2つの波形は、ゆっくりした 2 Hz の揺れ（両者で連動している）と細かい 20 Hz の揺れ（両者で無関係）が重なったものです。これを従来の相関係数に通すと $\rho_{xy}=0.42$ という数字がひとつ出てくるだけで、右上のように「2 Hz と 20 Hz が混ざって平均された」情報しか残りません。一方、右下のようにコヒーレンス関数を通せば、2 Hz では $\gamma^2 \approx 1$、20 Hz では $\gamma^2 \approx 0$ と、周波数軸に沿って結びつきの強さが分解されて見えてきます。この記事は、右下のグラフをどう計算し、どう読むかの話です。
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
コヒーレンス関数とは — 周波数ごとの相関係数
相関係数がこぼしてしまうもの
まず、普通の相関係数を思い出しましょう。2つの確率変数 $x, y$（平均 0 とします）に対して、ピアソンの相関係数は
$$ \rho_{xy} = \frac{E[xy]}{\sqrt{E[x^2]\,E[y^2]}} $$
でした。分子の $E[xy]$ が「連動の強さ」、分母が「それぞれの大きさ」で、割り算によってスケールに依存しない $-1 \le \rho_{xy} \le 1$ の量になります。この $\rho_{xy}^2$ は回帰分析における決定係数 $R^2$、すなわち「$y$ の分散のうち $x$ で線形に説明できる割合」に一致します。
ここで大事なのは、$E[xy]$ が 時間全体・周波数全体を一括で足し合わせた量 だという点です。冒頭の橋の例で言えば、2 Hz での強い連動と 20 Hz での無関係さが、足し算の中で混ざってしまいます。極端な場合、低周波では正に連動し高周波では負に連動していると、両者が打ち消し合って $\rho_{xy} \approx 0$ になり、「無関係」という完全に誤った結論に至ることさえあります。
実際にシミュレーションで確かめてみましょう。2 Hz の共通成分（両センサで連動）と、20 Hz 帯・広帯域のセンサ固有の揺れ（互いに無関係）を混ぜた 2 信号を作り、時間領域と周波数領域の両方で見比べます。
左の波形はぱっと見では「似ているような、似ていないような」状態で、中央の散布図から計算した相関係数も $\rho_{xy} = 0.43$（$\rho_{xy}^2 = 0.18$）という中途半端な1つの数字にしかなりません。ところが右のコヒーレンスを見ると、2 Hz 付近だけが $\gamma^2 \approx 1.0$ に跳ね上がり、それ以外の周波数は 0 に張り付いています。$\rho_{xy}^2$ の水平な破線（0.18）は、この鋭い山と広大な谷を平らにならしてしまった値にすぎない、ということが一目でわかります。
イコライザーのつまみごとに相関を測る
では、どうすればよいでしょうか。素直な発想は次のとおりです。2つの信号をそれぞれ狭い帯域のバンドパスフィルタに通し、その帯域だけを取り出した波形どうしで相関係数を計算する。 これを中心周波数を少しずつずらしながら繰り返せば、周波数ごとの相関の強さがグラフとして得られます。オーディオのグラフィックイコライザーのつまみを 1 本ずつ選んで、「このバンドだけで見たら 2 つの音は似ているか？」と尋ねていくイメージです。
先ほどと同じ 2 信号を、実際に 3 つのバンドパスフィルタに通して帯域ごとに相関係数を計算してみます。
1–3 Hz のバンドだけを取り出すと、2 本の波形はほとんど重なり、散布図は細い直線状になって相関係数は 0.97（$\rho^2 = 0.93$）に達します。一方 8–12 Hz と 18–22 Hz のバンドでは散布図が丸い雲になり、相関係数はそれぞれ $-0.02$、$-0.01$ とほぼゼロです。同じ 2 信号でも「どのつまみを選ぶか」で答えがまったく変わる — この帯域ごとの $\rho^2$ を周波数の連続関数にしたものが、次節で定義するコヒーレンスにほかなりません。
コヒーレンス関数がやっているのは、まさにこれです。ただし、実際にフィルタバンクを組むかわりに、フーリエ変換によって周波数成分を取り出し、スペクトル領域で相関係数の式をそのまま書き下します。バンドパスの中心周波数を $f$ とすれば、上の $\rho_{xy}$ の分子は「$x$ と $y$ の周波数 $f$ 成分どうしの積の期待値」、分母は「それぞれの周波数 $f$ 成分のパワー」に置き換わります。この置き換えを厳密に書いたものが、次節で定義するマグニチュード二乗コヒーレンスです。
もうひとつ、名前に含まれる「マグニチュード二乗」の理由を先に言っておきます。周波数領域では信号成分が複素数になるため、分子の「積の期待値」も複素数になります。複素数には位相があり、その位相は $x$ と $y$ の時間差を表します。時間差の情報を残したまま 0〜1 の指標にはできないので、絶対値をとって二乗し、大きさの情報だけを取り出します（位相の情報は捨てるのではなく、別途「クロススペクトルの位相」として活用します）。
ここまでで、コヒーレンスが「周波数ごとの $R^2$」であるという直感が得られました。次は、その計算に必要な部品であるクロススペクトル密度を用意します。
準備：パワースペクトル密度とクロススペクトル密度
パワースペクトル密度の復習
定常な確率過程 $x(t)$ のパワースペクトル密度（PSD）$S_{xx}(f)$ は、自己相関関数 $R_{xx}(\tau) = E[x(t)x(t+\tau)]$ のフーリエ変換として定義されます（ウィーナー・ヒンチンの定理）。
$$ S_{xx}(f) = \int_{-\infty}^{\infty} R_{xx}(\tau)\, e^{-j2\pi f \tau}\, d\tau $$
実用上は、長さ $T$ の観測区間のフーリエ変換 $X_T(f)$ を使った次の表現のほうが扱いやすいです。
$$ S_{xx}(f) = \lim_{T \to \infty} \frac{1}{T} E\!\left[\,|X_T(f)|^2\,\right] $$
ここで期待値 $E[\cdot]$ が付いていることが決定的に重要です。$|X_T(f)|^2$ は 1 回の観測から得られる乱雑な量（ピリオドグラム）であり、$T$ をいくら大きくしてもばらつきが減りません。真の PSD にたどり着くには、独立な観測を何度も行って平均をとる 必要があります。これがウェルチ法が「信号を複数のセグメントに分割して平均する」という手続きをとる理由でした。この「期待値＝平均化」の必要性は、後でコヒーレンスにおいてさらに劇的な形で現れます。
クロススペクトル密度
2つの信号を扱うために、自己相関を相互相関に置き換えます。相互相関関数 $R_{xy}(\tau) = E[x(t)\,y(t+\tau)]$ のフーリエ変換が クロススペクトル密度（cross spectral density, CSD） $S_{xy}(f)$ です。
$$ S_{xy}(f) = \int_{-\infty}^{\infty} R_{xy}(\tau)\, e^{-j2\pi f \tau}\, d\tau $$
フーリエ変換を使った表現では
$$ \begin{equation} S_{xy}(f) = \lim_{T \to \infty} \frac{1}{T} E\!\left[\,X_T^{*}(f)\, Y_T(f)\,\right] \end{equation} $$
となります。$X_T^{*}$ は複素共役です。PSD の $|X_T|^2 = X_T^{*} X_T$ の片方を $Y_T$ に差し替えただけ、と思えば覚えやすいでしょう。
CSD は PSD と違って 複素数 です。その意味を極形式で書くと見えてきます。
$$ S_{xy}(f) = |S_{xy}(f)|\, e^{j\theta_{xy}(f)} $$
- 絶対値 $|S_{xy}(f)|$ は、周波数 $f$ において 2 つの信号がどれだけ連動しているかの「大きさ」
- 偏角 $\theta_{xy}(f)$ は、周波数 $f$ の成分がどれだけ位相差を持っているかを表し、時間遅れ $\tau_d$ が一定なら $\theta_{xy}(f) = -2\pi f \tau_d$ という直線になります
この「複素数のまま平均する」という操作が何をしているのかを、複素平面上で見ておきましょう。
左は 2 信号が線形に結びついている場合で、各セグメントの $X_i^{*}Y_i$（灰色の矢印）が同じ方向を向くため、平均（太い矢印）は $|\hat{S}_{xy}| = 0.90$ と長いまま残ります。中央は無関係な場合で、矢印の向きがバラバラなので平均はほぼ原点に潰れ、$0.12$ まで小さくなりました。右の棒グラフが示すとおり、大きさだけを足した値（破線）は両者でほとんど変わらないのに、複素数のまま平均した値（実線）だけが激減します。この非対称性こそがコヒーレンスの心臓部です。
なお、共役をどちらに付けるかは文献によって流儀が分かれます。本記事では式 (1) の $E[X^{*}Y]$ を採用します。この流儀だと、後で見るように線形系 $Y = H X$ に対して $S_{xy} = H S_{xx}$ という余計な共役の付かない綺麗な形になります。SciPy の scipy.signal.csd(x, y) もこの流儀です。
これで分子と分母の材料が揃いました。あとは相関係数と同じ形に組み立てるだけです。
マグニチュード二乗コヒーレンスの定義
相関係数 $\rho_{xy} = E[xy] / \sqrt{E[x^2]E[y^2]}$ の各項を、周波数 $f$ における対応物に置き換えます。$E[xy] \to S_{xy}(f)$、$E[x^2] \to S_{xx}(f)$、$E[y^2] \to S_{yy}(f)$ とすれば、複素コヒーレンス
$$ \gamma_{xy}(f) = \frac{S_{xy}(f)}{\sqrt{S_{xx}(f)\, S_{yy}(f)}} $$
が得られます。これは複素数なので、絶対値の二乗をとって実数の指標にします。これが マグニチュード二乗コヒーレンス（magnitude-squared coherence, MSC） です。
$$ \begin{equation} \gamma^2_{xy}(f) = \frac{|S_{xy}(f)|^2}{S_{xx}(f)\, S_{yy}(f)} \end{equation} $$
単に「コヒーレンス」と言った場合、実務ではこの $\gamma^2$ を指すことがほとんどです。定義を眺めると、各項の物理的な役割がはっきりします。分子 $|S_{xy}|^2$ は「周波数 $f$ で 2 信号が連動している分のパワーの二乗」、分母 $S_{xx}S_{yy}$ は「それぞれが単独で持っているパワーの積」です。両者の比をとることで、$x$ や $y$ の絶対的な大きさ（センサの感度やアンプのゲイン）に依存しない、純粋な「結びつきの強さ」だけが残ります。
単位にも注意しておきましょう。$S_{xy}$ の単位は $[x][y]/\mathrm{Hz}$、$S_{xx}$ は $[x]^2/\mathrm{Hz}$、$S_{yy}$ は $[y]^2/\mathrm{Hz}$ なので、分子は $[x]^2[y]^2/\mathrm{Hz}^2$、分母も $[x]^2[y]^2/\mathrm{Hz}^2$ となり、$\gamma^2$ は 無次元 です。これも相関係数と同じ性質で、加速度と電圧のように単位の異なる 2 信号を比べられる理由になっています。
定義から $\gamma^2 \ge 0$ は明らかですが、上限が 1 であることは自明ではありません。次にこれを証明します。
$0 \le \gamma^2 \le 1$ の証明
複素確率変数の内積
証明の鍵はコーシー・シュワルツの不等式ですが、対象が複素数なので、まず「内積」を正しく定義しておきます。ある固定した周波数 $f$ における 2 つの複素確率変数 $X = X_T(f)$、$Y = Y_T(f)$ を考え、内積を
$$ \langle X, Y \rangle = E[X^{*} Y] $$
と定めます。この $\langle \cdot, \cdot \rangle$ が内積の公理を満たすことを確認しましょう。
第一に、$\langle X, X \rangle = E[|X|^2] \ge 0$ であり、等号は $X = 0$（確率 1 で）のときに限ります。これは正定値性です。第二に、$\langle Y, X \rangle = E[Y^{*}X] = \overline{E[X^{*}Y]} = \overline{\langle X, Y \rangle}$ で、エルミート対称性が成り立ちます。第三に、第 2 引数について線形性 $\langle X, \alpha Y_1 + \beta Y_2 \rangle = \alpha \langle X, Y_1 \rangle + \beta \langle X, Y_2 \rangle$ が明らかに成り立ちます。したがって、これは複素ヒルベルト空間の内積です。
内積が定義できたということは、コーシー・シュワルツの不等式がそのまま使えるということです。
コーシー・シュワルツの不等式の証明
不等式そのものを、内積の正定値性から導きます。任意の複素数 $\lambda$ に対して、ベクトル $Y – \lambda X$ のノルムの二乗は非負です。
$$ 0 \le \langle Y – \lambda X,\; Y – \lambda X \rangle $$
内積の線形性とエルミート対称性を使って展開します。$\langle Y-\lambda X, Y-\lambda X\rangle = \langle Y,Y\rangle – \lambda\langle Y,X\rangle – \lambda^{*}\langle X,Y\rangle + |\lambda|^2 \langle X,X\rangle$ なので、
$$ 0 \le E[|Y|^2] – \lambda\, \overline{E[X^{*}Y]} – \lambda^{*} E[X^{*}Y] + |\lambda|^2 E[|X|^2] $$
この不等式は 任意の $\lambda$ で成り立つので、右辺を最小にする $\lambda$ を代入すれば最も強い結論が得られます。$E[|X|^2] > 0$ のもとで右辺を $\lambda$ について最小化すると、$\lambda = E[X^{*}Y] / E[|X|^2]$ が最小点です（実際、右辺は $E[|X|^2]\left|\lambda – \frac{E[X^*Y]}{E[|X|^2]}\right|^2 + E[|Y|^2] – \frac{|E[X^*Y]|^2}{E[|X|^2]}$ と平方完成できます）。この $\lambda$ を代入すると、
$$ 0 \le E[|Y|^2] – \frac{|E[X^{*}Y]|^2}{E[|X|^2]} $$
両辺に $E[|X|^2]\;(>0)$ を掛けて移項すると、コーシー・シュワルツの不等式
$$ \begin{equation} |E[X^{*}Y]|^2 \le E[|X|^2]\, E[|Y|^2] \end{equation} $$
が得られます。
スペクトル密度への翻訳
式 (3) の各項を $1/T$ 倍して $T \to \infty$ の極限をとると、そのままスペクトル密度の言葉になります。
$$ |S_{xy}(f)|^2 \le S_{xx}(f)\, S_{yy}(f) $$
両辺を $S_{xx}S_{yy}\;(>0)$ で割れば、
$$ \gamma^2_{xy}(f) = \frac{|S_{xy}(f)|^2}{S_{xx}(f) S_{yy}(f)} \le 1 $$
分子が絶対値の二乗、分母が正なので $\gamma^2 \ge 0$ も自動的に成り立ちます。以上で $0 \le \gamma^2_{xy}(f) \le 1$ が示されました。$\square$
等号が成立するのはどんなときか
証明の議論をたどると、$\gamma^2 = 1$ となるのは平方完成した項がゼロ、すなわち
$$ Y = \lambda X \quad (\text{確率 } 1 \text{ で}) $$
が成り立つときだけだとわかります。$\lambda$ は複素定数です。周波数領域で $Y(f) = \lambda(f) X(f)$ ということは、時間領域では $y$ が $x$ を線形時不変システムに通しただけのもの であることを意味します（$\lambda(f)$ が周波数応答 $H(f)$ にあたります）。ゲインを掛けても、位相をずらしても、フィルタを通しても、$\gamma^2$ は 1 のままです。
逆に $\gamma^2 = 0$ となるのは $S_{xy}(f) = 0$、つまり周波数 $f$ において 2 信号が線形に無相関である場合です。ここで「線形に」という限定が付くのは相関係数と同じで、後で見るように 強い非線形の関係があっても $\gamma^2$ は 0 になり得ます。
理論上の性質はこれで押さえました。しかし、ここからが本題です。実際のデータで $\gamma^2$ を計算しようとすると、期待値 $E[\cdot]$ を推定量に置き換える必要があり、そこで最大の落とし穴が待ち構えています。
なぜ平均化が必須なのか — 単一セグメントでは $\hat{\gamma}^2 \equiv 1$
素直に推定量を作ってみる
手元にあるのは有限長のデータです。素直に考えれば、式 (2) の各スペクトル密度を、データから計算したピリオドグラムで置き換えればよさそうです。長さ $L$ のデータ 1 本から FFT で $X(f)$、$Y(f)$ を求め、
$$ \hat{S}_{xx}(f) = \frac{1}{L}|X(f)|^2, \quad \hat{S}_{yy}(f) = \frac{1}{L}|Y(f)|^2, \quad \hat{S}_{xy}(f) = \frac{1}{L}X^{*}(f) Y(f) $$
とします。これを定義式に代入しましょう。
$$ \hat{\gamma}^2(f) = \frac{|\hat{S}_{xy}(f)|^2}{\hat{S}_{xx}(f)\hat{S}_{yy}(f)} = \frac{\left|\frac{1}{L}X^{*}Y\right|^2}{\frac{1}{L}|X|^2 \cdot \frac{1}{L}|Y|^2} $$
分子を計算します。$|X^{*}Y|^2 = |X^{*}|^2|Y|^2 = |X|^2|Y|^2$ ですから、
$$ \hat{\gamma}^2(f) = \frac{\frac{1}{L^2}|X|^2 |Y|^2}{\frac{1}{L^2}|X|^2 |Y|^2} = 1 $$
分子と分母が完全に約分され、恒等的に 1 になってしまいました。 $X$ と $Y$ が何であっても、たとえまったく無関係な乱数どうしであっても、この推定量は全周波数で 1.0 を返します。
何が起きたのかを図形的に理解する
左は $K=1$ のときの複素平面で、$X(f)$ と $Y(f)$ がそれぞれ 1 本の矢印しかない状態です。2 本の間には位相差がありますが、それが「毎回同じか、それとも毎回バラバラか」は 1 回の観測からは判定しようがありません。中央の式のとおり分子と分母は完全に約分され、右のグラフのように無関係な乱数どうしでも全周波数で 1.0 の水平線になってしまいます。
なぜこうなるのか、複素平面で考えるとすっきりします。固定した周波数 $f$ において、$X$ と $Y$ は複素平面上の 2 つのベクトル（点）です。$X^{*}Y$ の絶対値は $|X||Y|$ であり、これは 2 つのベクトルの向きに一切依存しません。向きの情報はすべて偏角のほうに追いやられています。
つまり、単一の観測では 「$X$ と $Y$ の位相関係が毎回同じか、それともバラバラか」という情報がどこにも入っていない のです。コヒーレンスが本来測ろうとしているのはまさにこの「位相関係の再現性」であり、1 回の観測ではそれを測りようがありません。角度がひとつしかないとき、その角度が「安定している」かどうかは問えないのと同じです。
期待値 $E[\cdot]$ に立ち返れば、話は明快です。真の $S_{xy} = E[X^{*}Y]$ は複素数の 平均 であり、位相がバラバラなら足し合わせで打ち消し合って小さくなります。一方 $S_{xx} = E[|X|^2]$、$S_{yy} = E[|Y|^2]$ は非負の実数の平均なので打ち消し合いません。この「分子だけが打ち消し合う」という非対称性こそがコヒーレンスの本質であり、平均をとらなければその非対称性は決して現れません。
ウェルチ法によるセグメント平均
したがって、期待値をきちんと標本平均で近似する必要があります。ウェルチ法と同じ手続きです。長さ $N$ のデータを長さ $L$ のセグメントに（通常は 50% オーバーラップさせて）$K$ 個に分割し、各セグメントに窓関数 $w[n]$ を掛けて FFT します。$i$ 番目のセグメントのスペクトルを $X_i(f), Y_i(f)$ とすると、
$$ \hat{S}_{xy}(f) = \frac{1}{K}\sum_{i=1}^{K} X_i^{*}(f) Y_i(f), \quad \hat{S}_{xx}(f) = \frac{1}{K}\sum_{i=1}^{K} |X_i(f)|^2, \quad \hat{S}_{yy}(f) = \frac{1}{K}\sum_{i=1}^{K} |Y_i(f)|^2 $$
（正規化定数は共通なので $\gamma^2$ では約分されます）。これらを式 (2) に代入した推定量が
$$ \begin{equation} \hat{\gamma}^2(f) = \frac{\left|\sum_{i=1}^{K} X_i^{*}(f) Y_i(f)\right|^2}{\left(\sum_{i=1}^{K} |X_i(f)|^2\right)\left(\sum_{i=1}^{K} |Y_i(f)|^2\right)} \end{equation} $$
です。$K \ge 2$ であれば、分子の複素数の和が打ち消し合うため、もはや約分は起きません。$K = 1$ を代入すると先ほどの $\hat{\gamma}^2 = 1$ に戻ることも確認できます。
この事実の実務的な帰結は重大です。 「1周期分のデータで正確に測ろう」と考えて FFT 長を長くし、セグメント数を 1 にしてしまうと、コヒーレンスは 1.0 に張り付き、何も測れません。むしろ逆で、周波数分解能を多少犠牲にしてでもセグメント数 $K$ を確保するのが正しい設計です。
平均化回数 $K$ とバイアス
$K$ が有限である以上、$\hat{\gamma}^2$ にはバイアスが残ります。真のコヒーレンスが $\gamma^2 = 0$（完全に無関係）の場合、$K$ 個の独立なセグメントから計算した $\hat{\gamma}^2$ の確率密度は次のベータ分布型になることが知られています。
$$ p(u) = (K-1)(1-u)^{K-2}, \quad 0 \le u \le 1 $$
期待値を計算してみましょう。$\int_0^1 u\,(K-1)(1-u)^{K-2}\,du$ はベータ関数を使って $(K-1)B(2, K-1)$ と書けます。$B(2,K-1) = \frac{\Gamma(2)\Gamma(K-1)}{\Gamma(K+1)} = \frac{(K-2)!}{K!}$ なので、
$$ E[\hat{\gamma}^2] = (K-1)\cdot\frac{(K-2)!}{K!} = \frac{(K-1)!}{K!} = \frac{1}{K} $$
無関係な信号でも、平均的に $1/K$ のコヒーレンスが出てしまう ということです。$K=1$ なら $1/1 = 1$ で、先ほどの結果と一致します。$K=10$ なら 0.1、$K=100$ なら 0.01 です。一般の $\gamma^2$ に対しては、近似的に
$$ E[\hat{\gamma}^2] \approx \gamma^2 + \frac{(1-\gamma^2)^2}{K} $$
というバイアスが乗ります。真値が小さいほどバイアスの影響が相対的に大きい、という素直な形をしています。
左の確率密度 $p(u) = (K-1)(1-u)^{K-2}$ を見ると、$K=2$ では 0 から 1 まで一様（どんな値も等しく出る）で、$K$ を増やすにつれて 0 の近くに鋭く集中していく様子がわかります。右は期待値のグラフで、真値 0 の黒い曲線は $1/K$ に沿って下がり、$K=2$ では 0.5 もの見かけのコヒーレンスが出ることを示しています。真値が 0.9 のオレンジの曲線ではバイアス $(1-\gamma^2)^2/K$ が $K=2$ でも 0.005 しかなく、「真値が高いほどバイアスは効かない」という式の形がそのまま現れています。
有意水準のしきい値
バイアスがあるということは、「観測された $\hat{\gamma}^2 = 0.3$ は本当に意味のある結びつきなのか、それとも $K$ が小さいせいの見かけか」を判定する必要があるということです。上の確率密度を積分すると、$\gamma^2 = 0$ のもとで $\hat{\gamma}^2$ がしきい値 $c$ を超える確率は
$$ P(\hat{\gamma}^2 > c) = \int_c^1 (K-1)(1-u)^{K-2}\,du = (1-c)^{K-1} $$
です。これを有意水準 $\alpha$ と等しいと置いて $c$ について解くと、
$$ \begin{equation} \gamma^2_{\text{crit}} = 1 – \alpha^{1/(K-1)} \end{equation} $$
というしきい値が得られます。$\alpha = 0.05$、$K = 16$ なら $\gamma^2_{\text{crit}} = 1 – 0.05^{1/15} \approx 0.181$ です。つまり $K=16$ の平均化で $\hat{\gamma}^2 = 0.15$ が出ても、それは統計的には「無関係と区別がつかない」ということになります。
ひとつ注意点があります。ウェルチ法では通常 50% のオーバーラップを使いますが、オーバーラップしたセグメントは互いに独立ではありません。そのため実効的な自由度は名目上の $K$ より少なくなり、式 (5) のしきい値はやや甘め（低め）に出ます。後ほど Python で実測すると、$K=16$（ハン窓・50% オーバーラップ）での実際の 95 パーセンタイルは 0.190 程度と、理論値 0.181 よりわずかに高くなることが確認できます。
平均化の必要性がわかったところで、次はいよいよ「$\gamma^2$ が 1 未満になったとき、その数値は何を語っているのか」を定量化していきます。
線形系におけるコヒーレンス：$\gamma^2 = \mathrm{SNR}/(1+\mathrm{SNR})$
モデルの設定
コヒーレンスに具体的な意味を与える、最も基本的で最も有用なモデルを考えます。入力 $x[n]$ が線形時不変システム $h[n]$ を通り、出力側で独立な雑音 $n[n]$ が加わって $y[n]$ が観測される、という状況です。
$$ y[n] = (h * x)[n] + n[n] $$
たとえば、加振器で構造物を揺らして（$x$ が加振力）、加速度センサで応答を測る（$y$ が加速度）ような場面です。センサの電気的ノイズや、加振器とは無関係な環境振動が $n$ にあたります。仮定は 2 つだけです。(i) $h$ は線形時不変、(ii) $n$ は $x$ と無相関（$E[x \cdot n] = 0$）。
周波数領域では畳み込みが積になるので、
$$ Y(f) = H(f) X(f) + N(f) $$
と書けます。この $Y$ を式 (2) に代入していきます。
クロススペクトルの計算
まず分子側の $S_{xy}$ を求めます。定義 $S_{xy} = E[X^{*}Y]$ に $Y = HX + N$ を代入すると、
$$ S_{xy}(f) = E[X^{*}(HX + N)] = H\, E[|X|^2] + E[X^{*}N] $$
$H(f)$ は確率変数ではない（決定的なシステム特性）なので期待値の外に出せます。仮定 (ii) より $x$ と $n$ が無相関なので $E[X^{*}N] = 0$ となり、第 2 項が消えます。したがって
$$ \begin{equation} S_{xy}(f) = H(f)\, S_{xx}(f) \end{equation} $$
という、非常に見通しのよい関係が得られます。$E[X^{*}Y]$ という共役の付け方を選んだ理由がここにあります。この式は後で伝達関数の推定に直結します。
出力のパワースペクトル
次に分母側の $S_{yy}$ です。
$$ S_{yy}(f) = E[|HX+N|^2] = E\!\left[(HX+N)^{*}(HX+N)\right] $$
展開すると 4 項出てきます。
$$ S_{yy} = |H|^2 E[|X|^2] + H^{*}E[X^{*}N]^{*} + H E[X^{*}N] + E[|N|^2] $$
ここでも $E[X^{*}N] = 0$ により交差項が 2 つとも消えるので、
$$ \begin{equation} S_{yy}(f) = |H(f)|^2 S_{xx}(f) + S_{nn}(f) \end{equation} $$
「出力のパワー ＝ 信号がシステムを通ったパワー ＋ 雑音のパワー」という、直感どおりの分解です。
コヒーレンスの導出
式 (6) と式 (7) を式 (2) に代入します。
$$ \gamma^2(f) = \frac{|H S_{xx}|^2}{S_{xx}\left(|H|^2 S_{xx} + S_{nn}\right)} $$
分子は $|H|^2 S_{xx}^2$（$S_{xx}$ は実数の非負量なので絶対値の外に出ます）。分母と分子の $S_{xx}$ を 1 つ約分すると、
$$ \gamma^2(f) = \frac{|H|^2 S_{xx}}{|H|^2 S_{xx} + S_{nn}} $$
ここで、周波数 $f$ における信号対雑音比を
$$ \mathrm{SNR}(f) = \frac{|H(f)|^2 S_{xx}(f)}{S_{nn}(f)} $$
と定義します（分子は出力側に現れる「信号成分のパワー」、分母は「雑音のパワー」）。分子と分母を $S_{nn}$ で割ると、
$$ \begin{equation} \gamma^2(f) = \frac{\mathrm{SNR}(f)}{1 + \mathrm{SNR}(f)} \end{equation} $$
コヒーレンスは、出力側の信号対雑音比だけで決まる という、極めて明快な結論です。しかもシステム $H(f)$ の詳細（ゲイン・位相）は $\mathrm{SNR}$ の中に吸収されており、それ以外の形では現れません。
逆に解く：コヒーレンスから SNR を読む
式 (8) を $\mathrm{SNR}$ について解いてみましょう。$\gamma^2 (1 + \mathrm{SNR}) = \mathrm{SNR}$ より $\gamma^2 = \mathrm{SNR}(1 – \gamma^2)$、したがって
$$ \begin{equation} \mathrm{SNR}(f) = \frac{\gamma^2(f)}{1 – \gamma^2(f)} \end{equation} $$
デシベルなら $\mathrm{SNR}_{\mathrm{dB}} = 10\log_{10}\frac{\gamma^2}{1-\gamma^2}$ です。コヒーレンスを測れば、各周波数の SNR がそのまま読める ということです。これは測定現場で非常に強力な道具になります。「1200 Hz でコヒーレンスが 0.5 まで落ちている」という観測は、「その周波数では SNR が 0 dB しかない、つまり信号と雑音が同じ大きさだ」という診断に直訳できるからです。
具体的な対応を並べておきます。
| $\gamma^2$ | $\mathrm{SNR}$（真数） | $\mathrm{SNR}$ [dB] | 現場での解釈 |
|---|---|---|---|
| 0.10 | 0.111 | $-9.5$ | 雑音まみれ。この周波数の測定値は使えない |
| 0.50 | 1.0 | $0.0$ | 信号と雑音が同じ大きさ。かろうじてピークが見える程度 |
| 0.90 | 9.0 | $9.5$ | 実用の下限。伝達関数の推定に使える |
| 0.95 | 19.0 | $12.8$ | 良好 |
| 0.99 | 99.0 | $20.0$ | 非常に良好。よく管理された測定 |
この表をグラフにすると、コヒーレンスという 0〜1 の目盛りが SNR の物差しとしてどれだけ非線形なのかがよく見えます。
曲線は $\gamma^2 = 0.5$（SNR 0 dB）を境に、左半分では緩やかに、右半分では垂直に近く立ち上がります。つまり $\gamma^2$ が 0.9 から 0.99 へ動く「わずか 0.09」の変化が、SNR では 9.5 dB から 20 dB へ 10 dB 以上の改善に相当するということです。逆に言えば、0.99 と 0.995 の差を議論することにはほとんど意味がなく、注目すべきは緑の実用域（$\gamma^2 \ge 0.9$）から落ちてくる周波数がどこかという点になります。
振動試験の現場で「コヒーレンス 0.9 以上を確保せよ」という経験則がよく語られますが、その根拠がこの表です。$\gamma^2 = 0.9$ は SNR 約 10 dB に相当し、これを下回ると伝達関数の推定誤差が急速に悪化します。
さらに、式 (8) は $\gamma^2$ が SNR の単調増加関数であり、$\mathrm{SNR} \to \infty$ で $\gamma^2 \to 1$、$\mathrm{SNR} \to 0$ で $\gamma^2 \to 0$ となることも示しています。飽和が早い（SNR 20 dB でもう 0.99）ため、高 SNR 領域ではコヒーレンスの分解能が落ちる点は覚えておくとよいでしょう。
ここまでは雑音だけを考えました。しかし現実には、コヒーレンスが 1 を下回る原因は雑音だけではありません。次節でその全体像を整理します。
コヒーレンスが下がる4つの原因
$\gamma^2 < 1$ を見たとき、原因の切り分けができないと診断になりません。実務で遭遇する原因は、大きく次の 4 つに分類できます。
原因1：出力側の加法性雑音
すでに見たケースです。センサのノイズフロア、A/D 変換の量子化ノイズ、加振とは無関係な外乱振動などが該当します。特徴は $\gamma^2$ が周波数によって大きく変わる ことです。というのも、$\mathrm{SNR}(f) = |H(f)|^2 S_{xx}(f)/S_{nn}(f)$ は $|H(f)|^2$ に比例するため、システムの共振周波数（$|H|$ が大きい）ではコヒーレンスが高く、反共振周波数（$|H|$ が谷になる）では応答が雑音に埋もれてコヒーレンスがすとんと落ちます。実測の伝達関数とコヒーレンスを並べてプロットすると、この「反共振でコヒーレンスが落ちる」パターンが必ず見えます。
対処は素直です。加振レベルを上げて $S_{xx}$ を大きくする、平均化回数を増やす、あるいは入力信号のスペクトルを整形して弱い周波数帯にパワーを集中させます。
原因2：システムの非線形性
$y$ が $x$ の非線形関数である場合、コヒーレンスは 1 未満になります。定義の等号条件が $Y = \lambda X$（線形）だったので、これは当然です。興味深いのは、非線形性の「種類」によって $\gamma^2$ の落ち方がまったく違うことです。
入力 $x$ がガウス過程のとき、無記憶非線形 $y = g(x)$ に対しては バスガング（Bussgang）の定理 が使えます。この定理は、ガウス入力に対する非線形出力を
$$ g(x) = b\,x + d(x), \qquad b = E[g'(x)], \qquad E[x \cdot d(x)] = 0 $$
のように「最良の線形項 $bx$」と「$x$ と無相関な残差 $d$」に分解できることを述べています（$b = E[g'(x)]$ はスタインの補題から出ます）。$d$ は $x$ と無相関なので、コヒーレンスの観点からは原因1の雑音とまったく同じ役割を果たします。したがって
$$ \gamma^2 = \frac{b^2 \sigma_x^2}{\sigma_y^2} $$
となります。ここから 2 つの重要な帰結が出ます。
(a) 奇関数の非線形は部分的にコヒーレンスを下げる。 飽和特性 $g(x) = \tanh(ax)$ を考えます。単位分散のガウス入力に対して $a=3$ のとき、数値積分で $b = E[3\,\mathrm{sech}^2(3x)] = 0.7645$、$\sigma_y^2 = E[\tanh^2(3x)] = 0.7452$ となるので、$\gamma^2 = 0.7645^2/0.7452 = 0.784$ です。かなり強い飽和なのに、コヒーレンスは 0.78 と「そこそこ高い」値にとどまります。$a$ を変えると $a=1$ で 0.930、$a=5$ で 0.731 と、飽和が強くなるほど下がっていきます。
(b) 偶関数の非線形はコヒーレンスに一切映らない。 二乗特性 $g(x) = x^2$ では $b = E[2x] = 0$ です。したがって $\gamma^2 = 0$ — $y$ が $x$ から完全に決定論的に決まっているにもかかわらず、コヒーレンスはゼロ になります。これは相関係数の「非線形関係を見逃す」性質が周波数領域でもそのまま現れた例です。コヒーレンスが低いことは「関係がない」ことを意味しない、と肝に銘じるべき理由がここにあります。
この 2 つの帰結を図で並べておきます。
左の $\tanh(3x)$ では、データ点の雲の中心を貫く「最良の線形項」$bx$（$b = 0.764$）が確かに残っており、この傾きが $\gamma^2 = 0.784$ を生みます。中央の $y = x^2$ では放物線が左右対称なため、どんな傾きの直線を引いても平均的にはフィットせず、最良の線形項は $b = 0$ すなわち傾きゼロの直線になってしまいます（点線は直流成分 $E[y]=1$ で、これはコヒーレンスに寄与しません）。右のグラフは飽和の強さ $a$ を振ったときの $\gamma^2$ で、$a$ をいくら大きくしても 0.6 台までしか下がらず、奇関数の非線形は「そこそこ高いコヒーレンス」を保つことがわかります。
非線形性が疑われる場合は、入力振幅を変えて測定を繰り返すのが定石です。振幅を変えても $\gamma^2$ が変わらなければ雑音が主因、振幅を上げると $\gamma^2$ が下がるなら非線形が主因、と切り分けられます。
原因3：セグメント長に対して無視できない時間遅れ
これは最も見落とされやすく、かつ発見すると劇的に改善する原因です。$y[n] = x[n-D]$ という、雑音も非線形もない純粋な遅延だけの関係を考えます。理論上は $H(f) = e^{-j2\pi fD}$ という全通過系なので $\gamma^2 = 1$ のはずです。ところが実際に計算すると、$D$ が大きくなるにつれてコヒーレンスがどんどん下がります。
原因は、ウェルチ法が 有限長のセグメントを切り出す ことにあります。$x$ の $i$ 番目のセグメントには時刻 $[iR, iR+L)$ のデータが入りますが、$y$ の $i$ 番目のセグメントには時刻 $[iR-D, iR+L-D)$ のデータが入ります。$D$ がセグメント長 $L$ に近づくと、この 2 つの区間はほとんど重ならなくなり、同じセグメント番号どうしを比べているのに、実質的に無関係な区間を比べている 状態になります。
言葉だけではわかりにくいので、窓がずれる様子を絵にします。
左上で、青が $x$ 側のセグメント窓、赤が $y$ 側で実質的に使われている窓です。遅延 $D=64$ のぶんだけ横にずれており、紫に塗った重なり部分だけが「同じ時刻のデータを比べている」領域になります。左下はその重なりの積 $w[n]w[n+D]$ で、面積比が $r_w[64]/r_w[0] = 0.6592$。右のグラフのとおり、この比を二乗した $0.434$ がそのままコヒーレンスの値になります。雑音も非線形も一切ないのに、窓がずれただけで半分以下に落ちるのです。
これを定量化しましょう。$x$ を白色雑音（$E[x[n]x[m]] = \sigma^2\delta_{nm}$）とします。窓関数 $w[n]$（$n = 0,\dots,L-1$）を使ったセグメントスペクトルは
$$ X_i(f) = \sum_{n} w[n]\, x[n+iR]\, e^{-j2\pi fn}, \qquad Y_i(f) = \sum_{n} w[n]\, x[n+iR-D]\, e^{-j2\pi fn} $$
です。$Y_i$ の総和を $m = n – D$ で置換すると $Y_i(f) = e^{-j2\pi fD}\sum_m w[m+D]\,x[m+iR]\,e^{-j2\pi fm}$ となります。ここでクロススペクトルの期待値を計算します。
$$ E[X_i^{*} Y_i] = e^{-j2\pi fD}\sum_n \sum_m w[n]\,w[m+D]\; E[x[n+iR]x[m+iR]]\; e^{-j2\pi f(m-n)} $$
白色性より $E[x[n+iR]x[m+iR]] = \sigma^2\delta_{nm}$ なので、$m = n$ の項だけが残り、指数関数も $e^{0}=1$ になります。
$$ E[X_i^{*} Y_i] = \sigma^2 e^{-j2\pi fD}\sum_n w[n]\,w[n+D] = \sigma^2 e^{-j2\pi fD}\, r_w[D] $$
ここで $r_w[D] = \sum_n w[n]w[n+D]$ は 窓関数の自己相関 です。同様に $E[|X_i|^2] = E[|Y_i|^2] = \sigma^2 \sum_n w[n]^2 = \sigma^2 r_w[0]$ となります。これらをコヒーレンスの定義に入れると、指数部の絶対値は 1 なので消えて、
$$ \begin{equation} \gamma^2 = \frac{\sigma^4 \, r_w[D]^2}{\sigma^2 r_w[0]\cdot \sigma^2 r_w[0]} = \left(\frac{r_w[D]}{r_w[0]}\right)^2 = \rho_w[D]^2 \end{equation} $$
遅延によるコヒーレンス低下は、窓関数の正規化自己相関の二乗そのもの だとわかりました。矩形窓なら $\rho_w[D] = 1 – D/L$ なので $\gamma^2 = (1-D/L)^2$ という単純な放物線です。ハン窓では、$L=256$ に対して $D=64$（セグメント長の 1/4）でも $\rho_w[64] = 0.659$、$\gamma^2 = 0.434$ と半分以下まで落ちます。$D \ge L$ では自己相関がゼロになり $\gamma^2 = 0$ です。後ほど Python で $L=256$、$D=64$ の実測値 0.433 が理論値 0.4345 とぴたり一致することを確認します。
対処は 2 段階です。まず、クロススペクトルの位相 $\theta_{xy}(f) = -2\pi f D$ の傾きから、あるいは相互相関関数のピーク位置から遅延 $D$ を推定します。次に、片方の信号を $D$ だけシフトして時間軸を揃えてから改めてコヒーレンスを計算します。あるいは、セグメント長 $L$ を $D$ の 10 倍以上に取るのでも構いません（ただしその分 $K$ が減るので、データ長とのトレードオフになります）。音響測定でマイクとスピーカーの距離が数メートルある場合や、通信系で伝搬遅延が大きい場合は、この前処理を欠かすと診断を完全に誤ります。
なお、同じ理屈は「窓によるスペクトル漏れ」にも当てはまります。強い正弦波成分があってスペクトル漏れが大きいと、$x$ 側と $y$ 側で漏れの様子が違うためにコヒーレンスが不当に下がることがあります。窓関数の選択（ハン窓やブラックマン窓）はここでも効いてきます。
原因4：測っていない第2の入力
最後は、システム自体は完全に線形だが、$y$ に寄与する入力が $x$ 以外にもある場合です。
$$ y = h_1 * x_1 + h_2 * x_2 $$
で、$x_1$ しか測定していないとします。$x_2$ が $x_1$ と無相関なら、$h_2 * x_2$ は原因1の雑音と数学的に同一の役割を果たし、$\gamma^2_{x_1 y} = \mathrm{SNR}/(1+\mathrm{SNR})$（ここで $\mathrm{SNR} = |H_1|^2 S_{11}/(|H_2|^2 S_{22})$）となります。つまり コヒーレンスは「雑音」と「未測定の入力」を区別できません。
この状況を扱うには、$x_2$ も測って多入力の枠組みに移行します。$x_1$ の寄与を取り除いた上で $x_2$ と $y$ の結びつきを測る 部分コヒーレンス（partial coherence）、および全入力をまとめたときの説明力を表す 多重コヒーレンス（multiple coherence） が用意されています。考え方は重回帰分析における偏相関・重相関とまったく同じで、周波数ごとにそれをやっていると思えば理解しやすいでしょう。
原因の切り分けができるようになったところで、次はコヒーレンスの最も重要な応用である「伝達関数推定の信頼度評価」を見ていきます。
応用：伝達関数推定の信頼度とコヒーレント出力パワー
$H_1$ 推定量
式 (6) の $S_{xy} = H S_{xx}$ を $H$ について解けば、伝達関数の推定量が直ちに得られます。
$$ \hat{H}_1(f) = \frac{\hat{S}_{xy}(f)}{\hat{S}_{xx}(f)} $$
これを $H_1$ 推定量 と呼びます。出力側の雑音に強い（雑音が $x$ と無相関なら不偏）という性質があるため、加振試験では標準的に使われます。対になる $H_2$ 推定量 $\hat{H}_2 = \hat{S}_{yy}/\hat{S}_{yx}$ は入力側の雑音に強く、両者の比が $\hat{H}_1/\hat{H}_2 = \hat{\gamma}^2$ になるという美しい関係も成り立ちます。$\gamma^2$ が 1 に近いほど $H_1$ と $H_2$ が一致し、離れるほど 2 つの推定がずれる、というわけです。
推定誤差の式
$H_1$ 推定量のばらつきは、コヒーレンスと平均化回数で定量化できます。ゲイン $|\hat{H}_1|$ の正規化ランダム誤差（相対標準偏差）は
$$ \begin{equation} \varepsilon\left[|\hat{H}_1(f)|\right] = \frac{\sqrt{1-\gamma^2(f)}}{|\gamma(f)|\sqrt{2K}} \end{equation} $$
で与えられます（Bendat & Piersol の古典的な結果）。この式の読み方は次のとおりです。
- $\gamma^2 \to 1$ で誤差 $\to 0$。完璧な測定なら 1 回の平均でも正確
- $\gamma^2 \to 0$ で誤差 $\to \infty$。コヒーレンスが低い周波数の伝達関数は信用してはいけない
- 誤差は $1/\sqrt{K}$ で減る。精度を 2 倍にしたければ測定時間は 4 倍必要
具体例を出しましょう。$\gamma^2 = 0.5$（SNR 0 dB）で $K = 99$ のとき、$\varepsilon = \sqrt{0.5}/(\sqrt{0.5}\sqrt{198}) = 1/\sqrt{198} = 0.071$、つまりゲインの相対誤差は 7.1% です。$\gamma^2 = 0.909$（SNR 10 dB）なら $\varepsilon = \sqrt{0.091}/(\sqrt{0.909}\sqrt{198}) = 0.0225$ で 2.3% まで改善します。モンテカルロ実験（$L=256$、$K=99$、300 試行）で $\hat{H}_1$ のばらつきを実測すると、それぞれ 7.3%、2.3% となり、式 (11) とよく一致します。
左の対数軸のグラフを見ると、$\gamma^2$ が 0.2 を下回るあたりから誤差が急激に立ち上がり、$K=99$ でも 10% を超えることがわかります。$K$ を 10 から 300 へ 30 倍にしても曲線が下へ平行移動するだけで（$1/\sqrt{K}$ なので約 5.5 分の 1）、コヒーレンスの低さを平均化で埋め合わせるのは高くつく、という事実が読み取れます。右はモンテカルロ実測との比較で、$\gamma^2 = 0.5$ で 7.3%、$\gamma^2 = 0.909$ で 2.3% と、式 (11) の理論値（7.1%、2.2%）にほぼ乗っています。
コヒーレンス自体の推定誤差も同様に、$\varepsilon[\hat{\gamma}^2] = \sqrt{2}(1-\gamma^2)/(|\gamma|\sqrt{K})$ と書けます。$K$ が小さいとコヒーレンスの値そのものが信用できない、という当たり前の事実の定量版です。
コヒーレント出力パワー
もうひとつ、実務で便利な量を紹介します。式 (7) の $S_{yy} = |H|^2 S_{xx} + S_{nn}$ において、第1項（$x$ で説明できる部分）は $\gamma^2 S_{yy}$ に等しくなります。実際、
$$ \gamma^2 S_{yy} = \frac{|H|^2 S_{xx}}{|H|^2 S_{xx} + S_{nn}}\cdot\left(|H|^2 S_{xx} + S_{nn}\right) = |H|^2 S_{xx} $$
です。この $\gamma^2(f) S_{yy}(f)$ を コヒーレント出力パワー（coherent output power） と呼びます。対して $(1-\gamma^2)S_{yy} = S_{nn}$ が 非コヒーレント出力パワー、つまり $x$ で説明できない残差のスペクトルです。
これは強力な診断ツールになります。機械の振動スペクトルを測り、「このピークはモータ由来か、それとも別の要因か」を知りたいとします。モータの回転信号を $x$ として $\gamma^2 S_{yy}$ と $(1-\gamma^2)S_{yy}$ を並べてプロットすれば、各ピークがどちらの側に属するかが一目でわかります。回帰分析で「説明された変動」と「残差」に分解するのと同じことを、周波数ごとにやっているわけです。
理論と応用が揃いました。ここからは Python で実際に手を動かし、これまでの主張がすべて数値的に再現されることを確認していきます。
Pythonでの実装
実験1：単一セグメントでは常に $\gamma^2 = 1$
まず、この記事で最も重要な事実を確かめます。完全に無関係な 2 つの乱数列に対して、セグメント数を変えながらコヒーレンスを計算します。
import numpy as np
import matplotlib, matplotlib.pyplot as plt
from scipy import signal
# 日本語フォント設定
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
plt.rcParams["font.family"] = cand
break
plt.rcParams["axes.unicode_minus"] = False
rng = np.random.default_rng(0)
fs, L = 1000.0, 256
plt.figure(figsize=(11, 5))
for K_target, color in zip([1, 4, 16, 64], ["crimson", "darkorange", "seagreen", "royalblue"]):
# K個のセグメントが取れる長さのデータ（オーバーラップなし）
n = L * K_target
x = rng.standard_normal(n) # 互いに完全に無関係
y = rng.standard_normal(n)
f, C = signal.coherence(x, y, fs=fs, nperseg=L, noverlap=0, window="hann")
plt.plot(f, C, color=color, lw=1.2,
label=f"セグメント数 K={K_target}（平均 {C.mean():.3f}, 理論 1/K={1/K_target:.3f}）")
plt.xlabel("周波数 [Hz]"); plt.ylabel("マグニチュード二乗コヒーレンス $\\gamma^2$")
plt.title("無関係な2信号のコヒーレンス：K=1 では恒等的に 1 になる")
plt.ylim(-0.05, 1.15); plt.grid(alpha=0.3); plt.legend(loc="upper right", fontsize=9)
plt.tight_layout(); plt.show()
このグラフを見ると、$K=1$ の曲線が全周波数で 1.0 の水平線になっているはずです。まったく無関係な乱数どうしなのに「完全に線形な関係がある」と主張しているわけで、推定量として破綻しています。一方、$K$ を増やすと曲線は下がっていき、平均値が理論値 $1/K$（0.25, 0.0625, 0.0156）にきれいに近づきます。コヒーレンスは平均化して初めて意味を持つ、という本記事の中心的主張が、この 1 枚で確認できます。
実験2：$\gamma^2 = \mathrm{SNR}/(1+\mathrm{SNR})$ の検証
次に、線形系＋白色雑音のモデルで式 (8) を検証します。ここではシステムを全通過（$h = \delta$）にして、SNR を全周波数で一定にします。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
rng = np.random.default_rng(1)
fs, L, N = 1000.0, 1024, 200000
snr_db_list = np.arange(-15, 25.1, 2.5)
measured, theory = [], []
for snr_db in snr_db_list:
x = rng.standard_normal(N) # 入力（白色）
sigma_n = 10 ** (-snr_db / 20) # 雑音の標準偏差
y = x + sigma_n * rng.standard_normal(N) # y = x + n
f, C = signal.coherence(x, y, fs=fs, nperseg=L)
measured.append(C.mean())
snr = 10 ** (snr_db / 10)
theory.append(snr / (1 + snr))
plt.figure(figsize=(10, 5.5))
plt.plot(snr_db_list, theory, "k-", lw=2.5, label=r"理論 $\gamma^2 = \mathrm{SNR}/(1+\mathrm{SNR})$")
plt.plot(snr_db_list, measured, "o", ms=7, mfc="none", mec="crimson", mew=2,
label="実測（Welch法, K≈389）")
for g2, txt in [(0.5, "SNR=0dB"), (0.9, "SNR≈9.5dB")]:
plt.axhline(g2, color="gray", ls=":", lw=1)
plt.text(-14.5, g2 + 0.015, f"$\\gamma^2$={g2}（{txt}）", fontsize=9, color="gray")
plt.xlabel("出力側の SNR [dB]"); plt.ylabel(r"コヒーレンス $\gamma^2$")
plt.title("コヒーレンスと SNR の対応関係")
plt.grid(alpha=0.3); plt.legend(); plt.ylim(0, 1.05)
plt.tight_layout(); plt.show()
実測点が理論曲線の上にぴたりと乗ります。実際の数値では、SNR $=-10$ dB で実測 0.0929 に対し理論 0.0909、SNR $=0$ dB で実測 0.498 に対し理論 0.500、SNR $=+10$ dB で実測 0.9089 に対し理論 0.9091 と、小数第 2 位まで一致します。低 SNR 側で実測がわずかに高めなのは、$K$ が有限であることによる $(1-\gamma^2)^2/K$ のバイアスです（$\gamma^2$ が小さい領域では約 $1/K = 0.0026$ の下駄になります）。逆に $\gamma^2$ が 0.5 を超えるとバイアスは 0.0006 以下まで小さくなり、周波数方向のばらつきに埋もれて符号すら定まらなくなります。また曲線の形から、低 SNR 側では $\gamma^2 \approx \mathrm{SNR}$ とほぼ線形に、高 SNR 側では急速に 1 に飽和することも読み取れます。
実験3：周波数ごとに SNR が違うケース
現実のシステムでは $|H(f)|$ が周波数によって変わるので、コヒーレンスも周波数の関数になります。ローパス系に白色雑音を加えて、この様子を見てみましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
rng = np.random.default_rng(2)
fs, L, N = 1000.0, 1024, 400000
b = signal.firwin(63, 150, fs=fs) # 150 Hz のローパス系
x = rng.standard_normal(N)
y_clean = signal.lfilter(b, 1, x)
y = y_clean + 0.05 * rng.standard_normal(N) # 出力側に白色雑音
f, Cxy = signal.coherence(x, y, fs=fs, nperseg=L)
f, Syy = signal.welch(y, fs=fs, nperseg=L)
w, H = signal.freqz(b, worN=len(f), fs=fs)
snr_theory = np.abs(H) ** 2 * 1.0 / (0.05 ** 2) # |H|^2 Sxx / Snn（Sxx, Snn は同じ正規化）
g2_theory = snr_theory / (1 + snr_theory)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.8))
ax[0].semilogy(f, Syy, lw=1, color="royalblue", label=r"出力PSD $S_{yy}$")
ax[0].semilogy(f, Cxy * Syy, lw=1, color="seagreen", label=r"コヒーレント出力パワー $\gamma^2 S_{yy}$")
ax[0].semilogy(f, (1 - Cxy) * Syy, lw=1, color="crimson", label=r"非コヒーレント $(1-\gamma^2)S_{yy}$")
ax[0].set_xlabel("周波数 [Hz]"); ax[0].set_ylabel("PSD"); ax[0].legend(fontsize=9)
ax[0].set_title("出力パワーの分解"); ax[0].grid(alpha=0.3)
ax[1].plot(f, g2_theory, "k-", lw=2.5, label="理論 SNR/(1+SNR)")
ax[1].plot(f, Cxy, lw=1, color="crimson", alpha=0.8, label="実測コヒーレンス")
ax[1].set_xlabel("周波数 [Hz]"); ax[1].set_ylabel(r"$\gamma^2$"); ax[1].legend(fontsize=9)
ax[1].set_title("通過域では1に近く、阻止域では0に落ちる"); ax[1].grid(alpha=0.3)
plt.tight_layout(); plt.show()
右図では、150 Hz 以下の通過域でコヒーレンスがほぼ 1、遷移域で急落し、阻止域では 0 付近に張り付きます。阻止域では $|H(f)| \approx 0$ なので出力に信号成分がほとんど届かず、観測されるのは雑音だけ — つまり SNR がほぼゼロだからです。左図では、コヒーレント出力パワー（緑）が通過域で $S_{yy}$ とほぼ重なり、阻止域では非コヒーレント成分（赤）が支配的になっている様子が見えます。この 2 本の曲線の交差点が、「その周波数の測定が信用できるかどうか」の境界を教えてくれます。
実験4：3つの崩れ方の比較（雑音・非線形・遅延）
本記事の核心である 3 条件を、同じ図の上で比較します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
rng = np.random.default_rng(3)
fs, L, N = 1000.0, 256, 400000
x = rng.standard_normal(N)
x = x / x.std()
cases = {
"(基準) 完全な線形 y=x": x.copy(),
"(1) 白色雑音混入 SNR=6dB": x + 10 ** (-6 / 20) * rng.standard_normal(N),
"(2) 非線形歪み y=tanh(3x)": np.tanh(3 * x),
"(2') 偶関数の非線形 y=x^2": x ** 2,
"(3) 遅延 D=64 (=L/4)": np.roll(x, 64),
}
plt.figure(figsize=(11, 5.5))
for (name, y), c in zip(cases.items(), ["black", "royalblue", "seagreen", "darkorange", "crimson"]):
f, C = signal.coherence(x, y, fs=fs, nperseg=L)
plt.plot(f, C, lw=1.4, color=c, label=f"{name} 平均={C.mean():.3f}")
plt.xlabel("周波数 [Hz]"); plt.ylabel(r"コヒーレンス $\gamma^2$")
plt.title("コヒーレンスが崩れる3つの原因（セグメント長 L=256）")
plt.ylim(-0.05, 1.1); plt.grid(alpha=0.3); plt.legend(fontsize=9, loc="center right")
plt.tight_layout(); plt.show()
得られる平均値は、基準 1.000、雑音 SNR 6 dB で 0.799（理論 $3.98/4.98 = 0.799$）、$\tanh(3x)$ で 0.785（バスガング理論 0.784）、$y = x^2$ で 0.0004、遅延 $D=64$ で 0.433 です。特に注目すべきは 2 点あります。第一に、$y = x^2$ は $x$ から完全に決まる決定論的な関係なのに $\gamma^2 \approx 0.0004$ で、これは $K \approx 3100$ のときのバイアス $1/K \approx 0.0003$ とほぼ同じ値 — コヒーレンスから見ると「無関係な乱数」と区別がつきません。第二に、遅延だけのケースは雑音も非線形もないのに 0.433 まで落ちており、$\tanh$ の強い飽和より深刻な低下を招いています。遅延の見落としがいかに危険かがわかります。
実験5：遅延によるコヒーレンス低下は窓の自己相関で決まる
式 (10) の予測を確かめます。ハン窓の正規化自己相関の二乗と、実測のコヒーレンスを重ねます。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
rng = np.random.default_rng(4)
fs, L, N = 1000.0, 256, 400000
x = rng.standard_normal(N)
# ハン窓の正規化自己相関 rho_w[D]
w = signal.get_window("hann", L)
ac = np.correlate(w, w, "full")[L - 1:] # ラグ 0〜L-1
rho = np.append(ac / ac[0], 0.0) # ラグ L では自己相関 0
delays = np.arange(0, L + 1, 8)
meas = []
for D in delays:
f, C = signal.coherence(x, np.roll(x, int(D)), fs=fs, nperseg=L)
meas.append(C.mean())
plt.figure(figsize=(10, 5.5))
plt.plot(delays, rho[delays] ** 2, "k-", lw=2.5, label=r"理論 $\rho_w[D]^2$（ハン窓）")
plt.plot(delays, (1 - delays / L) ** 2, "--", color="gray", lw=1.5, label=r"参考：矩形窓 $(1-D/L)^2$")
plt.plot(delays, meas, "o", ms=6, mfc="none", mec="crimson", mew=2, label="実測コヒーレンス（平均）")
plt.axvline(L / 4, color="seagreen", ls=":", lw=1.5)
plt.text(L / 4 + 3, 0.85, "D = L/4 でも\n$\\gamma^2$ は 0.43 まで低下", fontsize=9, color="seagreen")
plt.xlabel("遅延 D [サンプル]（セグメント長 L=256）"); plt.ylabel(r"コヒーレンス $\gamma^2$")
plt.title("遅延によるコヒーレンス低下は窓関数の自己相関の二乗で決まる")
plt.grid(alpha=0.3); plt.legend(); plt.ylim(-0.03, 1.05)
plt.tight_layout(); plt.show()
実測点（赤丸）が理論曲線（黒）の上にほぼ完全に乗ります。数値で見ると $D=16$ で実測 0.9474 / 理論 0.9499、$D=64$ で 0.4328 / 0.4345、$D=128$ で 0.0310 / 0.0278 と、小数第 2 位まで一致します。式 (10) の導出が正しかったことの明快な確認です。灰色の破線（矩形窓）と比べると、$D$ が小さいうち（およそ $D < L/6$）はハン窓のほうが高い値を保ちますが、そこから先は逆転してハン窓のほうが速く落ち、$D \approx 0.6L$ 以降はほぼゼロに張り付きます。窓の端が滑らかなぶん、ラグが大きい領域では重なりが急速に失われるためです。実務では、まずクロススペクトルの位相の傾きから $D$ を推定して補正するのが正解です。
実験6：有意水準のしきい値
最後に、式 (5) のしきい値が実際に使えるかを確認します。無関係な信号ペアを大量に生成し、$\hat{\gamma}^2$ の分布を調べます。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
rng = np.random.default_rng(5)
L = 256
Ks = [8, 16, 32, 64]
fig, ax = plt.subplots(1, 2, figsize=(13, 4.8))
p95_meas, p95_theo, mean_meas = [], [], []
for K in Ks:
n = L * (K + 1) // 2 # 50%オーバーラップでK個のセグメント
vals = []
for _ in range(300):
a, b = rng.standard_normal(n), rng.standard_normal(n)
f, C = signal.coherence(a, b, fs=1.0, nperseg=L)
vals.append(C[1:-1])
vals = np.concatenate(vals)
p95_meas.append(np.percentile(vals, 95))
p95_theo.append(1 - 0.05 ** (1 / (K - 1)))
mean_meas.append(vals.mean())
ax[0].hist(vals, bins=60, density=True, histtype="step", lw=1.6, label=f"K={K}")
ax[0].set_xlabel(r"$\hat{\gamma}^2$"); ax[0].set_ylabel("確率密度")
ax[0].set_title(r"無関係な信号での $\hat{\gamma}^2$ の分布"); ax[0].legend(); ax[0].grid(alpha=0.3)
ax[1].plot(Ks, p95_theo, "k-o", lw=2, label=r"理論しきい値 $1-\alpha^{1/(K-1)}$ ($\alpha$=0.05)")
ax[1].plot(Ks, p95_meas, "s--", color="crimson", lw=1.5, label="実測 95 パーセンタイル")
ax[1].plot(Ks, mean_meas, "^:", color="seagreen", lw=1.5, label=r"実測平均（理論 $1/K$）")
ax[1].plot(Ks, [1 / k for k in Ks], "x", color="gray", ms=9, label=r"理論平均 $1/K$")
ax[1].set_xscale("log", base=2); ax[1].set_xlabel("平均化回数 K"); ax[1].set_ylabel(r"$\gamma^2$")
ax[1].set_title("有意判定のしきい値とバイアス"); ax[1].legend(fontsize=9); ax[1].grid(alpha=0.3)
plt.tight_layout(); plt.show()
左図のヒストグラムは、$K$ が大きくなるほど分布が 0 の近くに集中する様子を示しています。右図では、実測の平均値が理論値 $1/K$ にほぼ一致し、95 パーセンタイルも理論しきい値と近い値になります（$K=16$ で実測 0.190 / 理論 0.181）。実測がわずかに高いのは、50% オーバーラップしたセグメントが互いに独立でないためで、実効自由度が名目 $K$ より少し小さいことの現れです。実務では理論しきい値を目安にしつつ、境界付近の判定には余裕を見るのが安全です。
これらの実験によって、本記事で導出したすべての式が数値的に裏付けられました。最後に、現場で使うときの注意点をまとめておきます。
実務でコヒーレンスを使うときのチェックリスト
導出と実験を踏まえて、測定現場で押さえるべき点を整理します。
1. まず平均化回数 $K$ を確認する。 コヒーレンスのグラフを見せられたら、真っ先に「何回平均したのか」を聞くべきです。$K$ がわからないコヒーレンスは解釈不能です。$K$ が 10 未満なら、$1/K > 0.1$ のバイアスが乗っているので、0.2 程度の値には何の意味もありません。目安として $K \ge 30$、できれば $K \ge 100$ を確保したいところです。
2. 周波数分解能とのトレードオフを意識する。 総データ長 $N$ が決まっているとき、セグメント長 $L$ を大きくすれば周波数分解能 $f_s/L$ は上がりますが、$K \approx 2N/L$ が減ってコヒーレンスの信頼性が落ちます。「鋭い共振ピークを分離したい」という要求と「コヒーレンスを信用したい」という要求は正面から衝突します。データを長く取る以外に解決策はありません。
3. 遅延を先に補正する。 実験5 で見たとおり、セグメント長の 1/4 の遅延でもコヒーレンスは半分以下に落ちます。相互相関のピークやクロススペクトル位相の傾きから遅延を推定し、時間軸を揃えてから計算します。これを怠ると「システムに問題がある」と誤診しがちです。
4. 低コヒーレンス＝無関係、ではない。 $y = x^2$ の例が示すとおり、偶関数型の非線形関係はコヒーレンスにまったく映りません。$\gamma^2$ が測っているのは飽くまで 線形な 結びつきです。非線形が疑われるなら、入力振幅を振ってコヒーレンスの変化を見る、高次スペクトル（バイスペクトル）を使う、といった別の道具が必要です。
5. 反共振でコヒーレンスが落ちるのは正常。 構造物の伝達関数測定では、反共振周波数で応答が小さくなるため必然的に SNR が下がり、コヒーレンスも落ちます。これは測定の失敗ではなく物理です。共振ピーク付近で高いコヒーレンスが取れていれば、その測定は成功と判断してよいでしょう。
6. コヒーレンスは因果関係を語らない。 $\gamma^2_{xy}$ は対称（$\gamma^2_{xy} = \gamma^2_{yx}$）なので、$x$ が $y$ を引き起こしたのか、その逆か、あるいは共通の第三の要因があるのかを区別できません。方向性を知りたい場合は、クロススペクトル位相の符号や、グレンジャー因果性・有向コヒーレンスといった非対称な指標を検討します。
まとめ
本記事では、コヒーレンス関数について、定義から証明、崩れ方の定量化まで解説しました。
- 定義：マグニチュード二乗コヒーレンスは $\gamma^2_{xy}(f) = |S_{xy}(f)|^2/(S_{xx}(f)S_{yy}(f))$ で、相関係数の二乗（決定係数）を周波数ごとに計算したものにあたる無次元量
- 上下限：複素確率変数の内積 $\langle X,Y\rangle = E[X^{*}Y]$ にコーシー・シュワルツの不等式を適用すると $0 \le \gamma^2 \le 1$ が示され、$\gamma^2 = 1$ の等号条件は $Y = \lambda X$、すなわち $y$ が $x$ の線形時不変な応答であること
- 平均化の必然性：単一セグメントでは $|X^{*}Y|^2 = |X|^2|Y|^2$ により分子分母が完全に約分され、無関係な信号でも $\hat{\gamma}^2 \equiv 1$ になる。位相関係の再現性を測るには複数セグメントの平均が不可欠
- バイアスと有意性：真値 0 のとき $E[\hat{\gamma}^2] = 1/K$、有意水準 $\alpha$ のしきい値は $1 – \alpha^{1/(K-1)}$
- 雑音との関係：線形系＋出力雑音のモデルでは $\gamma^2 = \mathrm{SNR}/(1+\mathrm{SNR})$、逆に $\mathrm{SNR} = \gamma^2/(1-\gamma^2)$。$\gamma^2 = 0.9$ は SNR 約 9.5 dB に相当
- 非線形：バスガングの定理により $\gamma^2 = b^2\sigma_x^2/\sigma_y^2$（$b = E[g'(x)]$）。$\tanh(3x)$ で 0.784 だが、偶関数の $x^2$ では $b=0$ となり $\gamma^2 = 0$
- 遅延：セグメント長 $L$ に対する遅延 $D$ による低下は、窓関数の正規化自己相関を使って $\gamma^2 = \rho_w[D]^2$ と厳密に書ける。ハン窓・$D = L/4$ で 0.43 まで低下する
- 応用：$H_1$ 推定量のゲイン誤差は $\sqrt{1-\gamma^2}/(|\gamma|\sqrt{2K})$。コヒーレント出力パワー $\gamma^2 S_{yy}$ で、出力のどの成分が入力で説明できるかを分解できる
コヒーレンスは、スペクトル解析を「1本の信号のパワーを見る」段階から「2本以上の信号の関係を見る」段階へ引き上げる、最初の一歩にあたる概念です。ここから先には、部分コヒーレンス・多重コヒーレンスによる多入力系の同定、時間変化を追う時間周波数コヒーレンス（ウェーブレットコヒーレンス）、因果の方向を扱う有向コヒーレンスなど、豊かな世界が広がっています。
次のステップとして、以下の記事も参考にしてください。