加振器で構造物を揺すり、力センサと加速度センサの信号をFFTして、出力を入力で割る。教科書どおりにやったのに、共振ピークの高さが測るたびに違う。しかも同じデータから「別の割り算」をすると、ピークの高さが数dB変わってしまう。実験で周波数応答を測ったことのある人なら、一度は経験する不気味な現象です。
原因は、割り算の仕方が1通りではないことにあります。入力にも出力にもセンサ雑音が乗っているとき、「どちらの雑音を前提に最小二乗を解くか」で答えが変わります。出力側の雑音を前提にすれば $H_1 = S_{xy}/S_{xx}$ が、入力側の雑音を前提にすれば $H_2 = S_{yy}/S_{yx}$ が出てきます。この2つは同じデータから計算しても一致せず、しかも真値からのずれる向きが逆です。
この使い分けは、実験モーダル解析（構造物の固有振動数と減衰比を実測から同定し、有限要素モデルを更新する）と、制御系設計のためのプラント同定（ゲイン余裕・位相余裕を実測周波数応答から読む）という2つの現場で、そのまま結果の信頼性を左右します。共振ピークを $H_1$ で読めば減衰比を過大に見積もり、反共振の谷を $H_2$ で読めば谷が埋まって見える。どちらも「測定器のせい」ではなく、推定量の性質として説明できる現象です。本記事では、この2つの推定量を最小二乗の正規方程式から導き、バイアスの向きを式で示し、コヒーレンス $\gamma^2$ が実は $H_1/H_2$ そのものであることを証明します。
本記事の内容
- 「出力を入力で割る」素朴な推定がなぜ破綻するのか
- 出力雑音を前提とした最小二乗から $H_1 = S_{xy}/S_{xx}$ を導出する
- 入力雑音を前提とした最小二乗から $H_2 = S_{yy}/S_{yx}$ を導出する
- $H_1$ は過小評価、$H_2$ は過大評価というバイアスの向きを式で示す
- 共振で $H_1$ が崩れ、反共振で $H_2$ が崩れる物理的な理由
- コヒーレンス $\gamma^2 = H_1/H_2$ の導出と、その実務的な読み方
- 総和最小二乗（TLS）による $H_v$ 推定の導出と、不偏になる条件
- 2自由度振動系を使ったPythonでの数値実験（入力SNR・出力SNRを振って比較）
全体像を1枚の図にまとめておきます。
上段の「真の世界」では $Y = HX$ が厳密に成り立っていますが、私たちが手にできるのは下段の測定入力 $x_m$ と測定出力 $y_m$ の2本だけです。真の信号と雑音を分離する手段はないので、推定量を作るには「どちら側の雑音を前提にするか」という仮定を必ず1つ置くことになります。その仮定の選び方が右端の2つの分岐であり、$H_1$ は必ず真値より小さく、$H_2$ は必ず真値より大きくなる——本記事はこの図を式で埋めていく作業だと考えてください。
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
「出力を入力で割る」だけでは伝達関数は測れない
線形時不変システムの周波数応答 $H(\omega)$ は、入力のフーリエ変換 $X(\omega)$ と出力のフーリエ変換 $Y(\omega)$ の比として定義されます。理屈のうえでは、測定した $X$ と $Y$ を割り算すれば $H$ が出るはずです。
$$ \hat{H}(\omega) = \frac{Y(\omega)}{X(\omega)} $$
ところが、これを実データでやると3つの意味で困ります。
1つ目は、ゼロ割りです。 加振信号が広帯域雑音なら、$X(\omega)$ はどの周波数線でもランダムな複素数で、たまたまほぼ0になる周波数線が必ず出てきます。そこで割ると $\hat H$ が発散します。
2つ目は、平均化ができないことです。 1回の測定だけでは雑音の影響を減らせません。ところが $Y/X$ を何回も計算して平均しても、比の平均は平均の比ではないので、素直には改善しません。
この2つが実際にどう表れるかを、後で使う2自由度系のデータで見ておきます。
上段の灰色の線が1ブロックだけの $Y/X$ で、黒の真値のまわりを大きく暴れています。とくに下段で入力の線スペクトル $|X_k|^2$ がたまたま小さくなった 12 Hz 付近では、推定値が真値から19 dBも跳ね上がっています。下段の破線は下位2%の水準で、広帯域雑音で加振するかぎり「たまたま小さい周波数線」は必ず一定割合で現れることが分かります。一方、青の $H_1$（120ブロック平均）は同じデータから作ったのに黒線にぴたりと重なっており、比を平均するのではなくスペクトルを平均してから割るという順序の違いが効いていることが読み取れます。
3つ目が本質的な問題で、「雑音がどちらに乗っているか」を無視していることです。 実測では、力センサが読む値も加速度センサが読む値も、真の物理量そのものではありません。センサ雑音、アンプの熱雑音、AD変換の量子化誤差、外乱（隣の装置の振動）が必ず混ざります。$Y/X$ という書き方は、暗黙のうちに「$X$ も $Y$ も完璧に測れている」と仮定しているのです。
ここで、統計学でおなじみの状況を思い出してください。散布図に直線を当てはめるとき、$y$ を $x$ で回帰した直線と、$x$ を $y$ で回帰して傾きを逆数にした直線は一致しません。前者は「$y$ 方向の誤差」を、後者は「$x$ 方向の誤差」を最小化しているからです。相関が完全でない限り、この2本の直線は必ず開きます。周波数応答の推定でまったく同じことが起きます。$H_1$ が前者、$H_2$ が後者に対応し、2本の直線の傾きの比が決定係数 $r^2$ になるのと同じように、$H_1/H_2$ がコヒーレンス $\gamma^2$ になります。
同じ点の集まりに対して、青（縦の破線＝出力方向の誤差を最小化）と赤（横の破線＝入力方向の誤差を最小化）の2本の直線が引けています。青は傾き0.88と真の1より小さく、赤は1.35と大きい——これが後で式で示す「$H_1$ は過小、$H_2$ は過大」の原型です。2本の傾きの比が0.65で、これがそのままコヒーレンス $\gamma^2$ に対応します。相関が完全なら2本は重なり、相関が弱いほど開いていきます。
つまり、周波数応答の推定は「割り算」ではなく「回帰」です。回帰である以上、誤差をどこに置くかというモデルの仮定が必ず必要で、その仮定が結果を変えます。この視点を持てば、$H_1$ と $H_2$ の違いは謎ではなくなります。
まずは雑音の入り方をきちんとモデル化し、記号を整理するところから始めましょう。
記号と雑音モデルの整理
対象は線形時不変（LTI）システムとします。真の入力 $x(t)$、真の出力 $y(t)$ の間には、周波数領域で厳密に
$$ \begin{equation} Y(\omega) = H(\omega) \, X(\omega) \end{equation} $$
が成り立つとします。$H(\omega)$ が求めたい未知の複素数です。ここに測定雑音を入れます。私たちが実際に手にできるのは、真の信号ではなく
$$ \begin{align} x_m(t) &= x(t) + u(t) \quad \text{（入力側の雑音 } u \text{ が加わった測定入力）} \\ y_m(t) &= y(t) + v(t) \quad \text{（出力側の雑音 } v \text{ が加わった測定出力）} \end{align} $$
です。$u$ は力センサ・電流センサ側の雑音、$v$ は加速度センサ・電圧プローブ側の雑音だと思ってください。以下では次の3つを仮定します。
- $u$ と $v$ は互いに無相関
- $u$、$v$ はいずれも真の入力 $x$ と無相関
- 各信号は定常で、平均値は0
この仮定は「雑音がシステムの動作と関係なく乗ってくる」という、素直な状況を表しています。
次にスペクトルの記号を決めます。信号を長さ $T$ のブロックに切り、$k$ 番目のブロックのDFTを $X_k(\omega)$、$Y_k(\omega)$ と書きます。パワースペクトル密度（PSD）とクロススペクトル密度（CSD）は、ブロック平均として
$$ \begin{align} \hat S_{xx}(\omega) &= \frac{1}{n_d}\sum_{k=1}^{n_d} \frac{2}{T}\,|X_k(\omega)|^2 \\ \hat S_{xy}(\omega) &= \frac{1}{n_d}\sum_{k=1}^{n_d} \frac{2}{T}\,X_k^*(\omega)\, Y_k(\omega) \end{align} $$
と推定します。$n_d$ はブロック数（平均回数）、$*$ は複素共役です。この推定こそがWelch法そのもので、窓関数の適用と重複区間の扱いはWelch法によるパワースペクトル密度推定の理論と実装で扱ったとおりです。
ここで規約に注意してください。本記事では $S_{xy} = E[X^* Y]$ とします（入力側に共役をつける）。この規約だと $S_{xy}$ の位相がそのまま $H$ の位相になり、式が最も素直になります。SciPyの scipy.signal.csd(x, y) もこの規約です。共役の定義から、ただちに
$$ S_{yx} = E[Y^* X] = \left(E[X^* Y]\right)^* = S_{xy}^* $$
という関係が成り立ちます。この「$S_{yx}$ は $S_{xy}$ の共役」という何気ない事実が、後でコヒーレンスの導出に効いてきます。
最後に、仮定1〜3から測定量のスペクトルを書き下しておきます。無相関な項の期待値が消えることを使うと
$$ \begin{align} S_{x_m x_m} &= S_{xx} + S_{uu} \\ S_{y_m y_m} &= |H|^2 S_{xx} + S_{vv} \\ S_{x_m y_m} &= E[(X+U)^*(HX+V)] = H\,S_{xx} \end{align} $$
となります。3本目が重要です。クロススペクトルだけは、入力雑音も出力雑音も影響しません。 $U^* H X$、$X^* V$、$U^* V$ の3項がすべて無相関ゆえに期待値0で落ちるからです。雑音に汚されないのはクロススペクトルの分子だけ、という非対称性が、この後の議論すべての出発点になります。
左と中央の棒には、真の信号による成分（青）の上に雑音（斜線）が必ず積み上がります。ところが右のクロススペクトルだけは斜線部分がゼロで、真の $H S_{xx}$ がそのまま残ります。$H_1$ の分子と $H_2$ の分母がともにこのクロススペクトルであることを思い出すと、「$H_1$ は分母だけが汚れる」「$H_2$ は分子だけが汚れる」という後の結論が、この1枚から先読みできます。
役者がそろいました。ここから、最小二乗という1本の原理から $H_1$ と $H_2$ の両方が出てくる様子を見ていきます。
H1推定 — 出力雑音を前提にした最小二乗
まず「入力は正確に測れているが、出力にだけ雑音が乗る」という状況を考えます。加振力を測る力センサは高品位で、応答を測る加速度センサ側に電気雑音や外乱が乗っている、という設定です。このとき、ある1本の周波数線 $\omega$ に注目すると、$k$ 番目のブロックについて
$$ Y_{m,k} = H X_k + V_k $$
という関係が成り立ちます。$X_k$ はブロックごとに違う値をとる既知量、$V_k$ は未知のランダムな誤差、$H$ は全ブロック共通の未知の複素定数です。これはまさに「切片なしの単回帰」を複素数でやる問題です。$n_d$ 個の観測点 $(X_k, Y_{m,k})$ に、原点を通る直線 $Y = HX$ を当てはめる、と言い換えてもかまいません。
誤差 $V_k$ の二乗和を最小にする $H$ を探します。目的関数は
$$ J(H) = \sum_{k=1}^{n_d} \left| Y_{m,k} – H X_k \right|^2 $$
です。複素数の絶対値の二乗は $|z|^2 = z z^*$ なので、展開すると
$$ J(H) = \sum_{k=1}^{n_d} \left( Y_{m,k} – H X_k \right)\left( Y_{m,k}^* – H^* X_k^* \right) $$
各項を掛け合わせて整理します。
$$ J(H) = \sum_{k=1}^{n_d} \left( |Y_{m,k}|^2 – H X_k Y_{m,k}^* – H^* X_k^* Y_{m,k} + |H|^2 |X_k|^2 \right) $$
ここで複素変数に関する微分の作法を使います。$H$ と $H^*$ を形式的に独立な変数とみなして $H^*$ で偏微分する、というWirtinger微分です（実部・虚部それぞれで微分して連立させても同じ結果になりますが、こちらのほうが圧倒的に速いです）。$|H|^2 = H H^*$ を $H^*$ で微分すると $H$ が残ることに注意して
$$ \frac{\partial J}{\partial H^*} = \sum_{k=1}^{n_d} \left( – X_k^* Y_{m,k} + H |X_k|^2 \right) = 0 $$
これが正規方程式です。$H$ について整理すると
$$ H \sum_{k=1}^{n_d} |X_k|^2 = \sum_{k=1}^{n_d} X_k^* Y_{m,k} $$
両辺を $n_d$ で割り、さらにPSDの正規化定数 $2/T$ を掛ければ、左辺の和は $\hat S_{xx}$、右辺の和は $\hat S_{xy}$ そのものです。したがって
$$ \begin{equation} \boxed{\;H_1(\omega) = \frac{S_{xy}(\omega)}{S_{xx}(\omega)}\;} \end{equation} $$
を得ます。クロススペクトルを入力のパワースペクトルで割る、という見慣れた式が、出力雑音を前提とした最小二乗解として導かれました。
この形は、ベクトル表記にすると見通しがよくなります。$\bm{x} = (X_1, \dots, X_{n_d})^T$、$\bm{y} = (Y_{m,1}, \dots, Y_{m,n_d})^T$ と置けば
$$ H_1 = (\bm{x}^H \bm{x})^{-1} \bm{x}^H \bm{y} $$
これは1変数の最小二乗解 $(\bm{A}^H\bm{A})^{-1}\bm{A}^H\bm{y}$ にほかなりません。幾何学的には、$n_d$ 次元の複素ベクトル空間で、観測ベクトル $\bm{y}$ を $\bm{x}$ の張る1次元部分空間へ直交射影していることになります。この見方は射影行列と最小二乗法の理論で扱った構造と同じで、「残差が $\bm{x}$ と直交する」という条件がそのまま正規方程式です。
薄い青の帯が「$H$ を動かして作れる $Hx$ の全体」で、赤い測定出力ベクトル $\bm{y}$ はその帯の外にあります。$\bm{y}$ に最も近い帯の中の点が緑の $H_1\bm{x}$ で、そのとき残差（灰色の破線）は帯と直角に交わります。「残差が $\bm{x}$ と直交する」という条件を式にしたものが正規方程式そのもので、$H_1$ は難しい最適化を解かなくても射影1回で求まることが分かります。
$H_1$ の重要な性質を確認しておきます。仮定より $S_{x y_m} = H S_{xx}$ でしたから
$$ H_1 = \frac{H S_{xx}}{S_{xx}} = H $$
出力にだけ雑音があるとき、$H_1$ は不偏推定量です。 出力雑音がどれだけ大きくても、平均回数を増やせば $H_1$ は真値に収束します。これが $H_1$ が実験モーダル解析でも制御系同定でも既定値として使われる理由です。
しかし、この結論は「入力は正確に測れている」という仮定に全面的に依存しています。仮定が破れたときに何が起きるかを見るために、対称な立場の推定量を作ってみましょう。
H2推定 — 入力雑音を前提にした最小二乗
今度は立場を入れ替えます。「出力は正確に測れているが、入力にだけ雑音が乗る」という状況です。実験では珍しくありません。加振器と構造物の間に挟んだ力センサは、取り付け治具の質量の慣性力を余計に拾いますし（質量キャンセルが必要になる理由です）、微小な力を測る場面ではセンサのノイズフロアが効いてきます。一方、加速度センサの信号は大きく、SN比に余裕があることがしばしばあります。
このとき、周波数線 $\omega$ でのモデルは
$$ X_{m,k} = \frac{1}{H} Y_k + U_k $$
となります。誤差 $U_k$ が乗っているのは $x$ の側なので、最小化すべきは $x$ 方向の残差です。未知数を $G = 1/H$ と置けば、これは前節とまったく同じ形の問題で、$X$ と $Y$ の役割だけが入れ替わっています。したがって正規方程式も、添字を入れ替えるだけで書けます。
$$ G \sum_{k=1}^{n_d} |Y_k|^2 = \sum_{k=1}^{n_d} Y_k^* X_{m,k} $$
スペクトルの記号に直すと
$$ G = \frac{S_{yx}}{S_{yy}} $$
求めたいのは $H = 1/G$ ですから、逆数をとって
$$ \begin{equation} \boxed{\;H_2(\omega) = \frac{S_{yy}(\omega)}{S_{yx}(\omega)}\;} \end{equation} $$
を得ます。ここで大事なのは、$H_2$ は「$y$ を $x$ で最小二乗した結果」ではないという点です。$x$ を $y$ で最小二乗してから逆数をとったものです。回帰の向きが逆なので、$H_1$ とは別物になります。散布図に2本の回帰直線が引ける話とまったく同じ構図です。
同じ60点のデータに対して、左は各点から直線まで縦に、右は横に線を下ろしており、最小化している量が90度違うことがひと目で分かります。その結果、当てはめた直線の傾きは左が0.74、右が1.34と大きく食い違います。データも真の関係も同一なのに答えが変わるのは、推定量の優劣ではなく「誤差をどちらの軸に置いたか」という仮定の違いだけによるものです。
$H_2$ が入力雑音に対して不偏であることも確認しておきましょう。入力にだけ雑音がある場合、$S_{y y} = |H|^2 S_{xx}$、$S_{y x_m} = E[Y^*(X+U)] = H^* S_{xx}$ なので
$$ H_2 = \frac{|H|^2 S_{xx}}{H^* S_{xx}} = \frac{H H^*}{H^*} = H $$
きれいに真値になります。入力にだけ雑音があるとき、$H_2$ は不偏推定量です。
ここまでで、対称な2つの結論が出そろいました。
- 出力雑音のみ → $H_1$ が不偏
- 入力雑音のみ → $H_2$ が不偏
では、それぞれ「仮定が外れている」場合、推定値は真値からどちらへずれるのでしょうか。次節でその向きを式で確定させます。
バイアスを式で比べる
両方の雑音が同時に存在する一般の場合を計算します。前節までに用意した3本の式
$$ S_{x_m x_m} = S_{xx} + S_{uu}, \quad S_{y_m y_m} = |H|^2 S_{xx} + S_{vv}, \quad S_{x_m y_m} = H S_{xx} $$
をそのまま定義に代入するだけです。まず $H_1$ について、分子には雑音が入らず、分母だけが入力雑音の分だけ膨らむことに注目します。
$$ H_1 = \frac{S_{x_m y_m}}{S_{x_m x_m}} = \frac{H S_{xx}}{S_{xx} + S_{uu}} $$
分子・分母を $S_{xx}$ で割って整理すると
$$ \begin{equation} H_1 = H \cdot \frac{1}{1 + S_{uu}/S_{xx}} \end{equation} $$
$S_{uu}/S_{xx} > 0$ なので、右側の係数は必ず1未満です。つまり $H_1$ は必ず真値より小さくなります（過小評価）。しかも係数は入力の雑音対信号比だけで決まり、出力雑音には一切依存しません。
次に $H_2$ です。今度は分母のクロススペクトルが雑音に汚れず、分子の $S_{y_m y_m}$ だけが出力雑音の分だけ膨らみます。
$$ H_2 = \frac{S_{y_m y_m}}{S_{y_m x_m}} = \frac{|H|^2 S_{xx} + S_{vv}}{H^* S_{xx}} $$
分子の第1項を $H^* S_{xx}$ で割ると $|H|^2/H^* = H$ になるので、$H$ でくくり出せます。
$$ \begin{equation} H_2 = H \left( 1 + \frac{S_{vv}}{|H|^2 S_{xx}} \right) \end{equation} $$
括弧の中は必ず1より大きいので、$H_2$ は必ず真値より大きくなります（過大評価）。そしてこちらは出力の雑音対信号比だけで決まり、入力雑音には依存しません。
結果を表にまとめます。
| 推定量 | 入力雑音 $S_{uu}$ の影響 | 出力雑音 $S_{vv}$ の影響 | バイアスの向き |
|---|---|---|---|
| $H_1 = S_{xy}/S_{xx}$ | $\times\ \dfrac{1}{1+S_{uu}/S_{xx}}$ | 影響なし（不偏） | 常に過小 |
| $H_2 = S_{yy}/S_{yx}$ | 影響なし（不偏） | $\times\ \left(1+\dfrac{S_{vv}}{\|H\|^2 S_{xx}}\right)$ | 常に過大 |
この表は、実験で使えるかなり強力な道具です。真値は必ず $H_1$ と $H_2$ の間にある、と言い切れるからです。
$$ |H_1| \le |H| \le |H_2| $$
さらに、$H_1$ と $H_2$ の位相は完全に一致します。$H_1$ の位相は $\arg S_{xy}$、$H_2$ の位相は $-\arg S_{yx} = -\arg S_{xy}^* = \arg S_{xy}$ だからです。つまり2つの推定量の違いは振幅だけに現れ、位相特性を読むぶんにはどちらを使っても同じ答えになります。共振周波数や位相余裕の読み取りには影響しないが、振幅（したがって減衰比やゲイン余裕）には効く、というのが実務上のインパクトです。
横軸はその周波数でのSN比で、青が $H_1$、赤が $H_2$ のずれです。青は0 dBの線を決して上に越えず、赤は決して下に潜らないので、真値は必ず2本の間に挟まれます。SN比が20 dBもあればどちらも0.1 dB以内で実用上区別できませんが、SN比が0 dBまで落ちると $H_1$ は $-6$ dB、$H_2$ は $+6$ dBと、合わせて12 dBもの開きが生じます。青い点は後の数値実験のケースB（入力SNR 10 dB）に対応し、理論値 $-0.83$ dBが実測 $-0.72 \sim -0.80$ dBとよく合っています。
ここまでは「雑音のパワーが与えられたときの話」でした。しかし実験では、入力雑音・出力雑音のパワーは周波数によらずほぼ一定（白色）であることが多いのに対し、信号のパワーは周波数によって何桁も変わります。この食い違いが、$H_1$ と $H_2$ の崩れ方に強い周波数依存性を生みます。
なぜ共振でH1、反共振でH2が崩れるのか
$H_1$ のバイアス係数 $1/(1 + S_{uu}/S_{xx})$ は、入力の信号対雑音比だけで決まりました。ここで問うべきは「入力のパワー $S_{xx}$ が小さくなるのはどの周波数か」です。
加振器（シェーカ）で構造物を揺する場面を考えます。加振器を「内部インピーダンス $Z_s$ を持つ力源」とみなすと、実際に構造物へ入る力は
$$ F(\omega) = \frac{F_0(\omega)}{1 + Z_s \, Y_{\mathrm{mob}}(\omega)}, \qquad Y_{\mathrm{mob}} = \frac{\text{速度}}{\text{力}} = j\omega H(\omega) $$
と書けます。$Y_{\mathrm{mob}}$ はモビリティ（駆動点アドミッタンス）で、共振周波数でピークを持ちます。構造物が共振していると、わずかな力で大きく動いてしまう＝構造物側のインピーダンスが極端に低い状態になるので、加振器はそこに力を伝えられません。結果として、加振力のスペクトル $S_{xx}$ は共振周波数でストンと落ちます。これが実験の世界でよく知られた 力のドロップアウト（force drop-out） です。後の数値実験では、共振点で加振力PSDが約 $-18$ dB も落ち込むのを確認します。
$S_{xx}$ が落ちれば $S_{uu}/S_{xx}$ が跳ね上がり、$H_1$ のバイアス係数が1から大きく外れます。だから $H_1$ は共振ピークを低く見積もるのです。ピークが低く出るということは、半値幅法で読む減衰比が過大になり、「実際より減衰の効いた構造」に見えてしまうことを意味します。モーダル試験で $H_1$ をそのまま信じるのが危険な理由がここにあります。
上段が測ろうとしている真の周波数応答で、19.67 Hzと51.50 Hzに共振、その間の31.83 Hzに反共振があります。下段を見ると、加振器の内部インピーダンスがない灰色のケースでは加振力PSDが平坦なのに対し、$Z_s = 300$ を入れたオレンジのケースでは共振周波数とぴったり同じ位置で18.1 dBも落ち込んでいます。落ち込む場所が「ピークを最も正確に測りたい場所」と一致してしまうのが、力のドロップアウトのいちばん厄介なところです。
なお、共振ピークが低く出るもうひとつの原因に漏れ誤差（leakage）があります。分解能帯域幅より共振の半値幅が狭いと、窓関数のメインローブがピークを平滑化して低く潰します。こちらは雑音とは無関係の系統誤差で、$H_1$ か $H_2$ かを選んでも消えません。ブロック長を伸ばして分解能を上げるか、正弦掃引やバースト加振で周期性を確保して漏れ自体を断つのが正攻法です。本記事の数値実験では、分解能を十分に細かくして漏れの影響を無視できる状態にしたうえで、雑音の効果だけを取り出しています。
一方の $H_2$ のバイアス係数 $1 + S_{vv}/(|H|^2 S_{xx})$ を見ると、分母に $|H|^2 S_{xx}$、つまり真の出力パワーが入っています。出力雑音 $S_{vv}$ が周波数によらないノイズフロアだとすれば、この係数が大きくなるのは出力パワーが最も小さい周波数です。
駆動点の周波数応答は、2つの共振の間に反共振（antiresonance）を持ちます。反共振では応答が数十dB落ち込み、出力信号がノイズフロアに埋まります。したがって $S_{vv}/(|H|^2 S_{xx})$ が爆発的に大きくなり、$H_2$ は反共振の谷を埋めて、大幅に過大評価します。谷が埋まった周波数応答から零点（伝達関数の分子の根）を同定しようとすれば、まったく違う答えが出ます。反共振は構造物のどこに振動節ができるかを表す重要な情報なので、これも実害があります。
同じ論理は、応答が小さい高周波域全般にも当てはまります。質量制御域では $|H| \propto 1/\omega^2$ で応答が落ちていくので、どこかで必ずノイズフロアと交差し、そこから先は $H_2$ が使い物にならなくなります。
まとめると、崩れ方には次のような役割分担があります。
- $H_1$ が崩れる場所：入力パワーが落ちるところ＝共振点（力のドロップアウト）
- $H_2$ が崩れる場所：出力パワーが落ちるところ＝反共振点と高周波域（ノイズフロア）
「$H_1$ は共振で過小、$H_2$ は反共振で過大」という経験則は、こうしてバイアスの式とスペクトルの形から演繹できます。
ここで自然な疑問が湧きます。測定中に「今この周波数では入力雑音と出力雑音のどちらが効いているのか」を知る手がかりはないのでしょうか。実は $H_1$ と $H_2$ の比そのものが、その手がかりを与えてくれます。
コヒーレンス γ² = H1/H2 の導出
2つの信号がどれだけ線形に結びついているかを表す指標が、コヒーレンス関数です。定義は
$$ \begin{equation} \gamma^2(\omega) = \frac{|S_{xy}(\omega)|^2}{S_{xx}(\omega)\, S_{yy}(\omega)} \end{equation} $$
で、相関係数の二乗（決定係数）を周波数ごとに計算したものだと思ってください。$0 \le \gamma^2 \le 1$ であることは、コーシー・シュワルツの不等式 $|E[X^*Y]|^2 \le E[|X|^2]E[|Y|^2]$ から直ちに従います。
さて、$H_1$ と $H_2$ の比を計算してみましょう。定義をそのまま代入します。
$$ \frac{H_1}{H_2} = \frac{S_{xy}/S_{xx}}{S_{yy}/S_{yx}} $$
分数の割り算なので、分母の逆数を掛けます。
$$ \frac{H_1}{H_2} = \frac{S_{xy}}{S_{xx}} \cdot \frac{S_{yx}}{S_{yy}} = \frac{S_{xy} S_{yx}}{S_{xx} S_{yy}} $$
ここで前に確認した $S_{yx} = S_{xy}^*$ を使うと、分子は $S_{xy} S_{xy}^* = |S_{xy}|^2$ になります。
$$ \begin{equation} \boxed{\;\frac{H_1}{H_2} = \frac{|S_{xy}|^2}{S_{xx} S_{yy}} = \gamma^2\;} \end{equation} $$
コヒーレンスは、2つの推定量の比そのものでした。 導出はたった3行ですが、含んでいる情報は濃いです。
第一に、$\gamma^2$ は実数で $0$ 以上 $1$ 以下なので、この式から $|H_1| \le |H_2|$ が自動的に従います。前節で表からたどり着いた不等式が、恒等式として再導出されたことになります。
第二に、$H_1$ と $H_2$ の位相が一致することも、この式から即座に言えます。左辺の比が正の実数なら、両者の偏角は同じでなければなりません。
第三に、コヒーレンスは推定の信頼幅そのものを教えてくれます。$|H_1|$ と $|H_2|$ の比が $\gamma^2$ なので、dB表記では
$$ 20\log_{10}|H_2| – 20\log_{10}|H_1| = -10\log_{10}\gamma^2 $$
例えば $\gamma^2 = 0.9$ なら両者の差は $0.46$ dB しかなく、どちらを使ってもほぼ同じです。$\gamma^2 = 0.5$ なら $3.0$ dB、$\gamma^2 = 0.1$ なら $10$ dB も開きます。真値はこの区間のどこかにあるわけですから、コヒーレンスが低い周波数では「$H_1$ でも $H_2$ でも、どのみち10dBの不確かさがある」と読めます。
コヒーレンスの分解 — なぜ向きは分からないのか
さらに深い理解のために、$\gamma^2$ を雑音モデルで書き下してみます。前節の3本の式を定義に代入すると
$$ \gamma^2 = \frac{|H S_{xx}|^2}{(S_{xx}+S_{uu})(|H|^2 S_{xx}+S_{vv})} $$
分子は $|H|^2 S_{xx}^2$ です。分母の第1因子を $S_{xx}$ で、第2因子を $|H|^2 S_{xx}$ でくくると
$$ \begin{equation} \gamma^2 = \underbrace{\frac{1}{1 + S_{uu}/S_{xx}}}_{\displaystyle \gamma^2_{\text{in}}} \cdot \underbrace{\frac{1}{1 + S_{vv}/(|H|^2 S_{xx})}}_{\displaystyle \gamma^2_{\text{out}}} \end{equation} $$
という、驚くほどきれいな積の形になります。$\gamma^2_{\text{in}}$ は入力側の劣化、$\gamma^2_{\text{out}}$ は出力側の劣化を表します。そしてバイアスの式と見比べると
$$ H_1 = \gamma^2_{\text{in}} \, H, \qquad H_2 = \frac{H}{\gamma^2_{\text{out}}} $$
と対応していることが分かります。積をとれば $H_1/H_2 = \gamma^2_{\text{in}}\gamma^2_{\text{out}} = \gamma^2$ で、先ほどの恒等式と一致します。
この分解が教えてくれる、実務上いちばん大事な事実はこれです。$\gamma^2$ の値だけを見ても、劣化が入力側で起きているのか出力側で起きているのかは分かりません。 $\gamma^2 = 0.5$ という同じ数字が、「入力雑音だけで $\gamma^2_{\text{in}}=0.5$」の場合も「出力雑音だけで $\gamma^2_{\text{out}}=0.5$」の場合もありえます。前者なら $H_2$ が正解で、後者なら $H_1$ が正解です。答えは正反対なのに、コヒーレンスは同じ値を示す。だからコヒーレンスは「どれだけ信じられないか」の指標であって、「どちらを使うべきか」の指標ではありません。使い分けを決めるのは、あくまで測定系についての物理的な知識です。
横軸が $\gamma^2_{\text{in}}$、縦軸が $\gamma^2_{\text{out}}$ で、色が観測されるコヒーレンス $\gamma^2 = \gamma^2_{\text{in}}\gamma^2_{\text{out}}$ です。白い破線は $\gamma^2 = 0.5$ の等高線で、右下の丸（出力側だけが悪く $H_1$ が正解）と左上の四角（入力側だけが悪く $H_2$ が正解）はどちらもこの線の上に乗っています。つまり測定器が表示する $\gamma^2$ の値だけからは、正反対の2つの状況を区別できません。等高線に沿って動く自由度が残ることこそが、「コヒーレンスは量を教えるが向きは教えない」の図的な意味です。
コヒーレンスが下がる原因は雑音だけではない
なお、$\gamma^2 < 1$ になる原因は測定雑音以外にもあります。実験で低コヒーレンスに出会ったら、次の候補も疑ってください。
- 非線形性 — がたつき、摩擦、大振幅での剛性変化。入力レベルを変えて測り直すと、$\gamma^2$ の下がり方が変わることで見分けられます。
- 漏れ誤差 — 前節で触れた窓関数由来の系統誤差。共振の鋭いピークで $\gamma^2$ が局所的にへこむのが典型的な症状です。ブロック長を変えて再測定すると変化します。
- 多入力 — 測っていない別の加振源が同時に働いている。$y$ のうち $x$ で説明できない成分が増えます。
- 時間遅れ — 遅延がブロック長に対して無視できないと、ブロック端で相関が切れます。
平均回数が足りないとコヒーレンスは嘘をつく
最後に、絶対に外せない落とし穴を1つ。平均回数 $n_d = 1$ のとき、コヒーレンスは恒等的に1になります。 実際、1ブロックだけなら
$$ |\hat S_{xy}|^2 = |X^* Y|^2 = |X|^2 |Y|^2 = \hat S_{xx} \hat S_{yy} $$
なので、定義から $\hat\gamma^2 \equiv 1$ です。同時に $H_1 = X^*Y/|X|^2 = Y/X$、$H_2 = |Y|^2/(Y^*X) = Y/X$ となって、$H_1$ と $H_2$ も完全に一致してしまいます。これは「システムが完璧に線形で雑音がない」ことを意味しません。単に情報が足りていないだけです。
より一般に、$x$ と $y$ がまったく無相関でも、有限の $n_d$ で計算したコヒーレンスの期待値は
$$ E[\hat\gamma^2] \approx \frac{1}{n_d} $$
だけ持ち上がります。$n_d = 10$ なら、無相関な信号でも平均 $0.1$ のコヒーレンスが出るということです。この上げ底を「そこそこ相関がある」と誤読しないよう、平均回数は必ず記録し、できれば $n_d \ge 30$ を確保してください。この関係は後の数値実験で確かめます。
コヒーレンスによって「どれだけ信じられないか」は分かるようになりました。ただし、その不確かさにはバイアス（系統誤差）とばらつき（ランダム誤差）の2種類が混ざっています。次はその内訳を分けて見ます。
バイアスと分散のトレードオフ
$H_1$ は出力雑音に対して不偏でした。しかし「不偏」と「精度が高い」は別物です。不偏でもばらつきが大きければ、1回の測定値は当てになりません。
Bendat と Piersol による古典的な結果として、$|H_1|$ の正規化ランダム誤差（推定値の標準偏差を期待値で割ったもの）は
$$ \begin{equation} \varepsilon\!\left[|\hat H_1|\right] = \frac{\sqrt{1-\gamma^2}}{|\gamma|\,\sqrt{2 n_d}} \end{equation} $$
で与えられます。同じくコヒーレンス自身の正規化ランダム誤差は
$$ \varepsilon\!\left[\hat\gamma^2\right] = \frac{\sqrt{2}\,(1-\gamma^2)}{|\gamma|\,\sqrt{n_d}} $$
です。$H_1$ の式を読むと、精度を決める要素が2つあることが分かります。
1つ目は平均回数で、$1/\sqrt{n_d}$ で改善します。 誤差を半分にしたければ測定時間を4倍にする必要があり、これは避けられません。
2つ目はコヒーレンスで、$\sqrt{1-\gamma^2}/\gamma$ という形で効きます。 $\gamma^2$ が1に近ければこの係数は0に近づきますが、$\gamma^2$ が小さくなると急激に発散します。具体的に計算すると、$\gamma^2 = 0.9$、$n_d = 30$ で $\varepsilon = 4.3\%$。$\gamma^2 = 0.5$、$n_d = 30$ なら $12.9\%$。$\gamma^2 = 0.1$、$n_d = 30$ では $38.7\%$（約 $\pm 3$ dB）です。
ここが重要なポイントです。反共振や高周波域でコヒーレンスが $0.01$ 程度まで落ちると、$H_1$ は不偏であっても、1回の測定値は数dB単位で暴れます。 つまり $H_1$ が反共振でも真値の近くに見えるのは「平均的にはそうなる」という話であって、実測のカーブは低コヒーレンス帯でガタガタになります。後の数値実験でも、出力雑音のみのケース $\mathrm{A}$ において、反共振点で $H_1$ が $-3.0$ dB ずれる様子が観察されます。これはバイアスではなくランダム誤差で、$\gamma^2 = 0.011$、$n_d \approx 240$ を上の式に入れると $\varepsilon \approx 43\%$（約 $3$ dB）となり、ぴったり説明がつきます。同じ帯域での $H_2$ の $+36$ dB は、これとは桁が違う系統的なずれです。
4本の曲線は平均回数 $n_d$ の違いで、いずれも左（低コヒーレンス）へ行くほど急激に立ち上がります。$n_d = 30$ の青い線では $\gamma^2 = 0.9$ で4.3%だった誤差が、$\gamma^2 = 0.1$ では38.7%と9倍に膨れます。$n_d$ を10倍にしても曲線は $1/\sqrt{10} \approx 0.32$ 倍に下がるだけなので、コヒーレンスの低さは測定時間ではなかなか埋め合わせられません。オレンジの点は後の数値実験のケースAの反共振点（$\gamma^2 = 0.0109$、$n_d \approx 240$）で、43%＝振幅で約3 dBのばらつきが予測されます。
実務的には、次のように切り分けて考えると整理できます。
- バイアス（系統誤差）：平均回数を増やしても消えない。雑音の入り方に対する推定量の選択で対処する。
- ランダム誤差（分散）：平均回数を増やせば $1/\sqrt{n_d}$ で減る。測定時間で買える。
- 漏れ誤差：どちらでもない系統誤差。窓関数・分解能・加振波形の設計で対処する。
さて、ここまでで「$H_1$ か $H_2$ か、雑音の入り方で選ぶ」という指針は立ちました。しかし現実の測定では、入力にも出力にも雑音があります。両方あるとき、$H_1$ も $H_2$ も偏るなら、その中間を狙う推定量は作れないのでしょうか。
Hv推定 — 総和最小二乗による総合推定
$H_1$ は $y$ 方向の残差だけを、$H_2$ は $x$ 方向の残差だけを最小化していました。両方に誤差があるなら、両方向の残差を同時に最小化すればよいはずです。これが総和最小二乗（Total Least Squares, TLS）の発想で、得られる推定量を $H_v$（$v$ は vector の頭文字。固有ベクトルとして得られることに由来します）と呼びます。
直交距離としての定式化
複素平面上の点 $(X_k, Y_k)$ から直線 $Y = HX$ までの「垂直距離」を考えます。点と直線の距離の公式から、直線 $Y – HX = 0$ までの距離の二乗は
$$ d_k^2 = \frac{|Y_k – H X_k|^2}{1 + |H|^2} $$
です。分母の $1+|H|^2$ が、縦方向の残差を直交方向の残差に変換する係数になっています。これを全ブロックについて足したものが目的関数です。
$$ J_v(H) = \frac{1}{1+|H|^2}\sum_{k=1}^{n_d}\left|Y_k – H X_k\right|^2 $$
和の部分をスペクトルで書き直します。$\sum|Y_k – HX_k|^2$ を展開して $n_d$ で割れば
$$ S_{yy} – H S_{yx} – H^* S_{xy} + |H|^2 S_{xx} $$
なので、目的関数は
$$ J_v(H) = \frac{S_{yy} – H S_{yx} – H^* S_{xy} + |H|^2 S_{xx}}{1 + |H|^2} $$
となります。ここでも $H^*$ による微分をとります。分子を $N$、分母を $D = 1+HH^*$ と置くと、$\partial N/\partial H^* = -S_{xy} + H S_{xx}$、$\partial D/\partial H^* = H$ なので、商の微分から
$$ \frac{\partial J_v}{\partial H^*} = \frac{(H S_{xx} – S_{xy})D – N H}{D^2} = 0 $$
分子が0になる条件を書き下します。
$$ (H S_{xx} – S_{xy})(1 + |H|^2) = H\left(S_{yy} – H S_{yx} – H^* S_{xy} + |H|^2 S_{xx}\right) $$
左辺を展開すると $H S_{xx} + |H|^2 H S_{xx} – S_{xy} – |H|^2 S_{xy}$、右辺を展開すると $H S_{yy} – H^2 S_{yx} – |H|^2 S_{xy} + |H|^2 H S_{xx}$ です。両辺に共通する $|H|^2 H S_{xx}$ と $-|H|^2 S_{xy}$ が消えるので、残るのは
$$ H S_{xx} – S_{xy} = H S_{yy} – H^2 S_{yx} $$
移項して $H$ についての二次方程式に整理します。
$$ \begin{equation} S_{yx}\,H^2 + (S_{xx} – S_{yy})\,H – S_{xy} = 0 \end{equation} $$
見事に $|H|^2$ の項が消え、素直な複素二次方程式になりました。$S_{xy}S_{yx} = |S_{xy}|^2$ が実数であることを使えば、解の公式から
$$ H_v = \frac{(S_{yy} – S_{xx}) \pm \sqrt{(S_{yy}-S_{xx})^2 + 4|S_{xy}|^2}}{2 S_{yx}} $$
判別式の中身は非負なので平方根は実数で、最小値を与えるのは $+$ の側です。
スケーリングという落とし穴
ところが、この式をそのまま使うと期待外れの結果になります。理由は単位系です。入力が力 [N]、出力が変位 [m] だと、$S_{xx}$ と $S_{yy}$ は10桁近く違います。「$x$ 方向の残差と $y$ 方向の残差を同じ重みで足す」という発想は、両者が同じ単位・同じ雑音レベルであることを暗黙に仮定しているので、スケールが違うと片方が完全に無視されてしまうのです。実際、$S_{yy} \ll S_{xx}$ のときはこの $H_v$ が $H_1$ とほぼ同じ値になり、TLSの意味がなくなります。
正しくは、出力をスケール $c$ で割って $\tilde Y = Y/c$ とし、スケール後の空間でTLSを解いてから戻します。$\tilde S_{yy} = S_{yy}/c^2$、$\tilde S_{xy} = S_{xy}/c$、$\tilde H = H/c$ を先ほどの二次方程式に代入すると
$$ \frac{S_{yx}}{c}\cdot\frac{H^2}{c^2} + \left(S_{xx} – \frac{S_{yy}}{c^2}\right)\frac{H}{c} – \frac{S_{xy}}{c} = 0 $$
両辺に $c^3$ を掛けて整理すれば
$$ S_{yx}\,H^2 + (c^2 S_{xx} – S_{yy})\,H – c^2 S_{xy} = 0 $$
となり、解は
$$ \begin{equation} \boxed{\;H_v = \frac{(S_{yy} – c^2 S_{xx}) + \sqrt{(S_{yy} – c^2 S_{xx})^2 + 4 c^2 |S_{xy}|^2}}{2 S_{yx}}\;} \end{equation} $$
です。スケール $c$ の意味を確定させましょう。スケール後の出力雑音の分散は $S_{vv}/c^2$、入力雑音の分散は $S_{uu}$ です。TLSが「両方向の誤差を同じ重みで扱う」推定である以上、この2つが等しくなるように選ぶべきで
$$ c^2 = \frac{S_{vv}}{S_{uu}} $$
つまり $c$ は出力雑音と入力雑音のレベル比です。
H1・H2 は Hv の両極端
$c$ を極端な値にすると何が起きるか見てみましょう。$c \to \infty$（出力雑音が入力雑音に比べて圧倒的に大きい）とします。平方根の中身を $c^2 S_{xx}$ でくくると
$$ \sqrt{(c^2S_{xx})^2\left(1 + \frac{4|S_{xy}|^2}{c^2 S_{xx}^2}\right)} \approx c^2 S_{xx}\left(1 + \frac{2|S_{xy}|^2}{c^2 S_{xx}^2}\right) $$
と近似できます（$\sqrt{1+\epsilon}\approx 1+\epsilon/2$ を使いました）。分子は $-c^2 S_{xx} + c^2 S_{xx} + 2|S_{xy}|^2/S_{xx} = 2|S_{xy}|^2/S_{xx}$ となるので
$$ H_v \to \frac{2|S_{xy}|^2/S_{xx}}{2S_{yx}} = \frac{S_{xy}S_{yx}}{S_{xx}S_{yx}} = \frac{S_{xy}}{S_{xx}} = H_1 $$
逆に $c \to 0$（入力雑音が圧倒的）なら、分子は $S_{yy} + |S_{yy}| = 2S_{yy}$ になるので
$$ H_v \to \frac{2S_{yy}}{2S_{yx}} = H_2 $$
$H_1$ と $H_2$ は、$H_v$ という1つの族の両端だったわけです。3つの推定量がひとつの枠組みに収まりました。
横軸は正しい $c$ を1としたスケールの倍率で、緑の曲線が $H_v$ です。左端では赤い $H_2$ の水平線に、右端では青い $H_1$ の水平線にぴたりと重なり、$H_v$ が2つの推定量をつなぐ連続な族であることが確認できます。右パネル（反共振）では、正しい $c$（オレンジの点線）を選んだときだけ黒い真値の線の近くを通り、そこから1桁ずれるだけで誤差が数十dB変わってしまうことが読み取れます。紫の一点鎖線は次に述べる幾何平均型の $c$ で、正しい $c$ とは明らかに別の場所にあります。
正しい c を使えば Hv は不偏
$H_v$ の価値は、$c^2 = S_{vv}/S_{uu}$ を正しく与えれば、入力雑音と出力雑音が両方あっても不偏になる点にあります。証明しましょう。測定量のスペクトルを代入します。まず $S_{yy} – c^2 S_{xx}$ の部分を計算すると
$$ \left(|H|^2 S_{xx} + S_{vv}\right) – \frac{S_{vv}}{S_{uu}}\left(S_{xx} + S_{uu}\right) = |H|^2 S_{xx} + S_{vv} – c^2 S_{xx} – S_{vv} = S_{xx}\left(|H|^2 – c^2\right) $$
$S_{vv}$ が打ち消し合って消えるところが気持ちいいポイントです。次に平方根の中身を計算します。$|S_{xy}|^2 = |H|^2 S_{xx}^2$ を代入すると
$$ S_{xx}^2\left(|H|^2 – c^2\right)^2 + 4c^2 |H|^2 S_{xx}^2 = S_{xx}^2\left[\left(|H|^2-c^2\right)^2 + 4c^2|H|^2\right] $$
角括弧の中は $(|H|^2+c^2)^2$ の展開そのもの（$a^2-2ab+b^2+4ab = a^2+2ab+b^2$）ですから
$$ \sqrt{\;\cdot\;} = S_{xx}\left(|H|^2 + c^2\right) $$
したがって分子は $S_{xx}(|H|^2-c^2) + S_{xx}(|H|^2+c^2) = 2 S_{xx}|H|^2$ となり、分母は $2S_{yx} = 2H^* S_{xx}$ なので
$$ H_v = \frac{2S_{xx}|H|^2}{2H^*S_{xx}} = \frac{H H^*}{H^*} = H $$
雑音のレベル比さえ合っていれば、雑音の絶対量がどれだけ大きくても $H_v$ は不偏です。 これは強い結果です。$H_1$ が「出力雑音のみ（$c \to \infty$）」、$H_2$ が「入力雑音のみ（$c \to 0$）」という特別な比に対する不偏推定量だったのに対し、$H_v$ は任意の比に対応します。
幾何平均という便利な特別解
もうひとつ、実用上覚えておくと便利な特別なケースがあります。$c^2 = S_{yy}/S_{xx}$、つまり入力と出力をそれぞれのパワーで正規化する選び方です。このとき $S_{yy} – c^2S_{xx} = 0$ なので、平方根だけが残り
$$ H_v = \frac{2c|S_{xy}|}{2S_{yx}} = c\,\frac{|S_{xy}|}{S_{xy}^*} = \sqrt{\frac{S_{yy}}{S_{xx}}}\cdot\frac{S_{xy}}{|S_{xy}|} $$
振幅だけ取り出すと $|H_v| = \sqrt{S_{yy}/S_{xx}}$ です。ここで $|H_1||H_2| = \frac{|S_{xy}|}{S_{xx}}\cdot\frac{S_{yy}}{|S_{xy}|} = \frac{S_{yy}}{S_{xx}}$ であることに気づくと
$$ \begin{equation} H_v = \sqrt{H_1 H_2} \qquad \text{（位相は $H_1$、$H_2$ と共通）} \end{equation} $$
$H_1$ と $H_2$ の幾何平均です。$\gamma^2 = H_1/H_2$ を使えば $|H_v| = |H_1|/\gamma = |H_2|\,\gamma$ とも書けます。$H_1$ と $H_2$ をdBで表示したとき、ちょうど真ん中を通る曲線がこれです。
この選び方は「入力と出力の雑音対信号比が等しい」という仮定に対応します。両チャンネルに同じダイナミックレンジのADCを使い、どちらも同程度にレンジを使い切っているなら、悪くない近似です。ただし反共振のように出力だけがノイズフロアに沈む帯域では、この仮定が成り立たないので $H_v$ も外れます。後の数値実験でも、幾何平均型は反共振で $+16$ dB ずれるのに対し、正しい $c$ を使った $H_v$ は $-3.4$ dB で踏みとどまる、という差がはっきり出ます。
理論はここまでです。実際に2自由度の振動系を用意して、雑音の入り方を変えながら3つの推定量の振る舞いを見ていきましょう。
Pythonでの実装
真の周波数応答が解析的に分かっている系を使って、推定値と真値を直接比べます。題材は駆動点で加振・計測する2自由度振動系です。この系は共振を2つ、その間に反共振を1つ持つので、「共振で $H_1$ が崩れ、反共振で $H_2$ が崩れる」という現象を1枚の図に収められます。
対象システムの定義
質量 $m_1$ が地面にばね $k_1$ でつながれ、さらにばね $k_2$ を介して質量 $m_2$ がぶら下がる構成です。運動方程式を周波数領域に移すと、動剛性行列 $\bm{Z}(\omega) = -\omega^2 \bm{M} + j\omega\bm{C} + \bm{K}$ の逆行列がレセプタンス行列になります。今回は質量1に力を加えて質量1の変位を測る「駆動点」の応答 $H(\omega) = [\bm{Z}^{-1}]_{11}$ を対象にします。
import numpy as np
# 2自由度系: 地面 -k1- m1 -k2- m2
m1 = m2 = 1.0
k1 = k2 = 4.0e4
M = np.array([[m1, 0.0], [0.0, m2]])
K = np.array([[k1 + k2, -k2], [-k2, k2]])
# 固有角振動数
w_n = np.sqrt(np.sort(np.linalg.eigvals(np.linalg.inv(M) @ K).real))
# 比例減衰 C = alpha*M + beta*K を、両モードの減衰比が0.03になるよう決める
A = np.array([[1/w_n[0], w_n[0]], [1/w_n[1], w_n[1]]])
alpha, beta = np.linalg.solve(A, np.array([2*0.03, 2*0.03]))
C = alpha * M + beta * K
def frf(f):
"""駆動点レセプタンス H(f) = X1(f)/F1(f) [m/N]"""
w = 2 * np.pi * np.asarray(f, dtype=float)
H = np.empty(w.shape, dtype=complex)
for i, wi in enumerate(w):
Z = -wi**2 * M + 1j * wi * C + K # 動剛性行列
H[i] = np.linalg.inv(Z)[0, 0]
return H
f_res = w_n / (2*np.pi) # 共振周波数
f_anti = np.sqrt(k2/m2) / (2*np.pi) # 反共振周波数
print(f"共振: {f_res[0]:.2f} Hz, {f_res[1]:.2f} Hz / 反共振: {f_anti:.2f} Hz")
実行すると 共振: 19.67 Hz, 51.50 Hz / 反共振: 31.83 Hz と表示されます。反共振が $\sqrt{k_2/m_2}/2\pi = 31.83$ Hz にぴったり来ているのは偶然ではありません。駆動点レセプタンスの分子は $k_2 – \omega^2 m_2$ になり、これが0になる周波数、すなわち「質量2とばね2だけで作る動吸振器が共振して、質量1を完全に押さえ込む周波数」が反共振です。減衰比は両モードとも $0.03$ に設定したので、半値幅は共振1で約 $1.2$ Hz。後で使う周波数分解能 $0.125$ Hz に対して十分広く、漏れ誤差が無視できる条件になっています。
計測データの生成
真の入力と真の出力を作ります。入力は白色雑音、出力はそれを周波数応答で色付けしたものです。全記録長のFFTを使って一括で作ることで、ブロック境界での不連続を避けています。あわせて、力のドロップアウトを模擬するオプションも仕込んでおきます。
import numpy as np
fs = 512 # サンプリング周波数 [Hz]
nper = 4096 # 1ブロック長（周波数分解能 fs/nper = 0.125 Hz）
n_blk = 120 # ブロック数
N = nper * n_blk
def generate(seed=0, Zs=0.0):
"""真の入力 x（加振力）と真の出力 y（変位）を作る。
Zs>0 で加振器の内部インピーダンスによる『力のドロップアウト』を模擬。"""
rng = np.random.default_rng(seed)
f_all = np.fft.rfftfreq(N, 1/fs)
H_all = frf(np.maximum(f_all, 1e-9))
U = np.fft.rfft(rng.standard_normal(N)) # 白色の加振指令
if Zs > 0.0:
Y_mob = 1j * 2*np.pi*np.maximum(f_all, 1e-9) * H_all # モビリティ
U = U / (1.0 + Zs * Y_mob) # 共振で力が抜ける
x = np.fft.irfft(U, n=N)
y = np.fft.irfft(U * H_all, n=N)
return x, y, rng
x, y, rng = generate(seed=0)
print(f"入力RMS={np.std(x):.4f} N, 出力RMS={np.std(y):.3e} m")
出力は 入力RMS=1.0011 N, 出力RMS=3.018e-05 m となり、入力と出力でパワーが10桁近く違うことが確認できます。この巨大なスケール差こそが、前節で述べた「$H_v$ にスケーリングが必要な理由」の実体です。何も考えずに $S_{xx}$ と $S_{yy}$ を足し引きする定式化は、この時点で破綻することが数字から見て取れます。
3つの推定量を計算する関数
Welch法でPSDとCSDを求め、$H_1$、$H_2$、$H_v$、コヒーレンスを一括で返す関数を用意します。
import numpy as np
from scipy import signal
def estimate(xm, ym, c2=None):
"""H1 / H2 / Hv とコヒーレンスを返す。c2 は Svv/Suu（未指定なら幾何平均型）。"""
kw = dict(fs=fs, nperseg=nper, noverlap=nper//2, window='hann')
f, Sxx = signal.welch(xm, **kw)
_, Syy = signal.welch(ym, **kw)
_, Sxy = signal.csd(xm, ym, **kw) # scipy は E[X* Y] の規約
Syx = np.conj(Sxy)
H1 = Sxy / Sxx
H2 = Syy / Syx
coh = np.abs(Sxy)**2 / (Sxx * Syy)
if c2 is None:
Hv = np.sqrt(Syy/Sxx) * Sxy/np.abs(Sxy) # 幾何平均型
else:
d = Syy - c2*Sxx
Hv = (d + np.sqrt(d**2 + 4*c2*np.abs(Sxy)**2)) / (2*Syx)
return f, H1, H2, Hv, coh
def db_err(H_hat, f):
"""真値からのずれを dB で返す"""
return 20*np.log10(np.abs(H_hat) / np.abs(frf(np.maximum(f, 1e-9))))
f, H1, H2, Hv, coh = estimate(x, y)
band = (f > 5) & (f < 100)
print("雑音なし 5-100 Hz RMS誤差[dB]:",
f"H1={np.sqrt(np.mean(db_err(H1,f)[band]**2)):.3f}",
f"H2={np.sqrt(np.mean(db_err(H2,f)[band]**2)):.3f}",
f"平均コヒーレンス={coh[band].mean():.4f}")
雑音を一切加えない状態では H1=0.011 H2=0.016 平均コヒーレンス=0.9994 となります。両推定量とも真値から $0.02$ dB 以内、コヒーレンスもほぼ1です。この結果が意味するのは、この実験系では漏れ誤差が事実上ゼロに抑えられているということで、これ以降に観測されるずれはすべて雑音由来だと安心して解釈できます。ベースラインを先に取るのは、数値実験でも実測でも欠かせない手順です。
ケースA・B: 片側だけに雑音があるとき
いよいよ本題です。ケースAでは出力にだけ白色雑音（PSDが出力ピークの $-30$ dB）を、ケースBでは入力にだけ白色雑音（入力SNR $10$ dB）を加えます。
import numpy as np
import matplotlib, matplotlib.pyplot as plt
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
if any(cand == fo.name for fo in matplotlib.font_manager.fontManager.ttflist):
plt.rcParams["font.family"] = cand
break
plt.rcParams["axes.unicode_minus"] = False
rng = np.random.default_rng(100)
Ht = frf(np.maximum(f, 1e-9))
peak_psd = np.max((np.abs(Ht)**2 * np.mean(signal.welch(x, fs=fs, nperseg=nper)[1]))[band])
# ケースA: 出力にだけ白色雑音（PSDが出力ピークの -30 dB）
sig_v = np.sqrt(peak_psd * 10**(-30/10) * (fs/2))
ymA = y + sig_v * rng.standard_normal(N)
fA, H1A, H2A, HvA, cohA = estimate(x, ymA)
# ケースB: 入力にだけ白色雑音（入力SNR 10 dB）
sig_u = np.std(x) * 10**(-10/20)
xmB = x + sig_u * rng.standard_normal(N)
fB, H1B, H2B, HvB, cohB = estimate(xmB, y)
for tag, ff, Ha, Hb, cc in [("A", fA, H1A, H2A, cohA), ("B", fB, H1B, H2B, cohB)]:
ea, eb = db_err(Ha, ff), db_err(Hb, ff)
ir = np.argmin(np.abs(ff - f_res[0])); ia = np.argmin(np.abs(ff - f_anti))
print(f"[{tag}] 共振{ff[ir]:.2f}Hz: H1={ea[ir]:+.2f}dB H2={eb[ir]:+.2f}dB γ²={cc[ir]:.4f} | "
f"反共振{ff[ia]:.2f}Hz: H1={ea[ia]:+.2f}dB H2={eb[ia]:+.2f}dB γ²={cc[ia]:.4f} | "
f"RMS誤差 H1={np.sqrt(np.mean(ea[band]**2)):.2f}dB H2={np.sqrt(np.mean(eb[band]**2)):.2f}dB")
得られた出力は次のとおりです。
[A] 共振19.62Hz: H1=-0.09dB H2=+0.03dB γ²=0.9857 | 反共振31.88Hz: H1=-2.99dB H2=+36.29dB γ²=0.0109 | RMS誤差 H1=0.56dB H2=10.23dB
[B] 共振19.62Hz: H1=-0.72dB H2=+0.10dB γ²=0.9100 | 反共振31.88Hz: H1=-0.80dB H2=-0.04dB γ²=0.9162 | RMS誤差 H1=0.84dB H2=0.13dB
この2行に、本記事の理論のほぼすべてが凝縮されています。読み取れることを整理します。
ケースA（出力雑音のみ）では、$H_1$ が圧勝します。 5〜100 Hz でのRMS誤差は $H_1$ が $0.56$ dB に対し $H_2$ は $10.23$ dB。特に反共振点では $H_2$ が $+36.3$ dB、つまり真値の約65倍という致命的な過大評価をしています。この値も理論から予言できます。出力雑音のみなら $H_1 = H$ なので、恒等式 $\gamma^2 = H_1/H_2$ から $H_2 = H/\gamma^2$。振幅のdB表記では $-20\log_{10}\gamma^2 = -20\log_{10}(0.0109) = 39.3$ dB になります。実測の $+36.3$ dB との差 $3.0$ dB は、$H_1$ 側に現れている $-2.99$ dB のずれとちょうど打ち消し合う量です（$H_2/H_1$ の比は $\gamma^2$ で厳密に決まるので、$H_1$ が $3$ dB 下振れれば $H_2$ の見かけの誤差も $3$ dB 小さくなります）。理論と実測が完全に整合しています。
その $H_1$ の $-2.99$ dB は、バイアスではなくランダム誤差です。 $\gamma^2 = 0.0109$、$n_d \approx 240$（50%重複で120ブロック）をBendat-Piersolの式に入れると $\varepsilon = \sqrt{1-0.0109}/(0.104\sqrt{480}) \approx 0.44$、すなわち約44%（振幅で約 $3$ dB 相当）のばらつきが予測されます。平均回数を増やせば消える誤差である一方、$H_2$ の $+36$ dB はいくら平均しても消えません。
ケースB（入力雑音のみ）では立場が完全に逆転します。 $H_2$ のRMS誤差 $0.13$ dB に対し $H_1$ は $0.84$ dB。しかも $H_1$ のずれは共振で $-0.72$ dB、反共振で $-0.80$ dB と、周波数によらずほぼ一定です。これは理論式 $H_1 = H/(1+S_{uu}/S_{xx})$ から予測されるとおりで、入力SNRが $10$ dB（$S_{uu}/S_{xx} = 0.1$）なら $20\log_{10}(1/1.1) = -0.83$ dB。実測値とぴったり一致します。入力が白色で入力雑音も白色なら、$H_1$ のバイアスは全周波数で同じ、という予言が数値で確認できました。
そして両ケースを並べたときに最も重要なのは、コヒーレンスの見た目からは、どちらのケースか判別できないことです。ケースAの反共振で $\gamma^2 = 0.011$、ケースBの反共振で $\gamma^2 = 0.916$ と値は違いますが、仮にケースBの雑音レベルを上げて $\gamma^2 = 0.011$ にしても、コヒーレンスの曲線からは「入力側が悪い」とは読めません。前節で述べた「$\gamma^2$ は劣化の量を教えるが向きは教えない」が、実験でも成り立っています。
図にすると、この対比はさらに鮮明になります。
fig, ax = plt.subplots(2, 2, figsize=(13, 8))
for j, (tag, ff, Ha, Hb, cc) in enumerate(
[("ケースA: 出力にだけ雑音", fA, H1A, H2A, cohA),
("ケースB: 入力にだけ雑音", fB, H1B, H2B, cohB)]):
ax[0, j].semilogy(ff, np.abs(frf(np.maximum(ff, 1e-9))), 'k-', lw=2.5,
label='真の周波数応答')
ax[0, j].semilogy(ff, np.abs(Ha), 'C0-', lw=1.2, label='$H_1 = S_{xy}/S_{xx}$')
ax[0, j].semilogy(ff, np.abs(Hb), 'C3-', lw=1.2, alpha=0.8,
label='$H_2 = S_{yy}/S_{yx}$')
ax[0, j].set_xlim(5, 100); ax[0, j].set_ylim(1e-7, 1e-3)
ax[0, j].set_title(tag); ax[0, j].set_ylabel('振幅 |H| [m/N]')
ax[0, j].legend(fontsize=9); ax[0, j].grid(alpha=0.3)
ax[1, j].plot(ff, cc, 'C2-', lw=1.2)
ax[1, j].set_xlim(5, 100); ax[1, j].set_ylim(0, 1.05)
ax[1, j].set_xlabel('周波数 [Hz]')
ax[1, j].set_ylabel('コヒーレンス $\\gamma^2$')
ax[1, j].grid(alpha=0.3)
plt.tight_layout()
plt.show()
左列（ケースA）を見ると、$H_1$ の曲線は真値の黒線にぴたりと重なり、反共振の深い谷まで追従します。一方 $H_2$ の赤い曲線は、共振ピーク付近では黒線と重なるのに、反共振では谷がまったく再現されず、平らに埋まってしまいます。下段のコヒーレンスは反共振と高周波域で急降下しており、$H_2$ が破綻する帯域とぴったり一致します。右列（ケースB）では、$H_2$ が黒線に重なり、$H_1$ が全帯域でわずかに下側にずれた平行線を描きます。ずれ方が「一様に下」なのか「特定の周波数で暴れる」なのかを見るだけで、雑音がどちら側に乗っているかを推測する手がかりになります。
ケースC・D: 両方に雑音があるときと Hv の効果
現実の測定は両方に雑音が乗ります。ケースCでは入力SNR $10$ dB と出力ノイズフロア $-30$ dB を同時にかけ、正しい $c^2 = S_{vv}/S_{uu}$ を与えた $H_v$ の効果を見ます。ケースDではさらに力のドロップアウトを加え、実測に近い状況を作ります。
import numpy as np
rng = np.random.default_rng(200)
# ケースC: 入力SNR 10 dB + 出力雑音フロア -30 dB
xmC = x + sig_u * rng.standard_normal(N)
ymC = y + sig_v * rng.standard_normal(N)
c2 = sig_v**2 / sig_u**2 # Svv/Suu（白色同士なので分散比）
fC, H1C, H2C, HvC, cohC = estimate(xmC, ymC, c2=c2)
_, _, _, HgC, _ = estimate(xmC, ymC) # 幾何平均型（c を誤ったケース）
# ケースD: 力のドロップアウトを入れた実測に近い設定
xD, yD, _ = generate(seed=7, Zs=300.0)
rng = np.random.default_rng(300)
sig_uD = np.std(xD) * 10**(-20/20)
_, SyyD = signal.welch(yD, fs=fs, nperseg=nper, noverlap=nper//2, window='hann')
sig_vD = np.sqrt(np.max(SyyD[band]) * 10**(-30/10) * (fs/2))
xmD = xD + sig_uD * rng.standard_normal(N)
ymD = yD + sig_vD * rng.standard_normal(N)
fD, H1D, H2D, HvD, cohD = estimate(xmD, ymD, c2=sig_vD**2/sig_uD**2)
for tag, ff, Ha, Hb, Hc, cc in [("C", fC, H1C, H2C, HvC, cohC),
("D", fD, H1D, H2D, HvD, cohD)]:
ea, eb, ec = db_err(Ha, ff), db_err(Hb, ff), db_err(Hc, ff)
i1 = np.argmin(np.abs(ff-f_res[0])); ia = np.argmin(np.abs(ff-f_anti))
print(f"[{tag}] 共振1: H1={ea[i1]:+.2f} H2={eb[i1]:+.2f} Hv={ec[i1]:+.2f} γ²={cc[i1]:.3f} | "
f"反共振: H1={ea[ia]:+.2f} H2={eb[ia]:+.2f} Hv={ec[ia]:+.2f} γ²={cc[ia]:.3f}")
print(f" RMS誤差[dB] H1={np.sqrt(np.mean(ea[band]**2)):.2f} "
f"H2={np.sqrt(np.mean(eb[band]**2)):.2f} Hv={np.sqrt(np.mean(ec[band]**2)):.2f}")
出力は次のようになります。
[C] 共振1: H1=-0.89 H2=-0.02 Hv=-0.03 γ²=0.905 | 反共振: H1=-4.15 H2=+36.75 Hv=-3.40 γ²=0.009
RMS誤差[dB] H1=1.10 H2=10.22 Hv=0.75
[D] 共振1: H1=-6.57 H2=+0.18 Hv=+0.17 γ²=0.460 | 反共振: H1=-0.39 H2=+3.32 Hv=-0.32 γ²=0.653
RMS誤差[dB] H1=1.34 H2=0.59 Hv=0.13
ケースCでは、正しい $c$ を与えた $H_v$ が全帯域で最良（RMS $0.75$ dB）です。 共振では $H_2$ と同じく $-0.03$ dB の精度を保ちつつ、反共振では $H_2$ の $+36.8$ dB という破綻を避けて $-3.4$ dB に踏みとどまっています。理論で証明した「正しい雑音比を与えれば $H_v$ は不偏」が、そのまま数値に現れました。残る $-3.4$ dB はコヒーレンス $0.009$ に対応するランダム誤差です。
ケースDは、実測でいちばん厄介なパターンです。 力のドロップアウトによって加振力PSDが共振点で $-18.1$ dB 落ち込むため、共振点での入力SNRが極端に悪化し、$H_1$ は共振1で $-6.57$ dB、共振2で $-6.60$ dB もピークを低く見積もります。これは「$H_1$ で読んだ減衰比が実際の2倍以上になる」ことを意味し、モーダル試験としては致命的です。一方、同じ共振点で $H_2$ は $+0.18$ dB とほぼ正確。しかしその $H_2$ も、反共振では $+3.32$ dB 谷を埋めてしまいます。同じ1本の測定データの中で、周波数によって最良の推定量が入れ替わるという、教科書どおりの状況です。$H_v$ は共振で $+0.17$ dB、反共振で $-0.32$ dB、全帯域RMSで $0.13$ dB と、両方の良いところを取っています。
下段の誤差カーブを見ると、ケースC（左）では赤い $H_2$ が反共振で $+36.7$ dB、高周波でも右肩上がりに $+20$ dB近くまで外れるのに対し、緑の $H_v$ は0 dBの線に貼りついたままです。ケースD（右）では逆に青い $H_1$ が共振1で $-6.6$ dBの深い谷を作り、$H_2$ は反共振でわずかに持ち上がります。同じ1本の測定の中で「共振では $H_1$ が負け、反共振では $H_2$ が負ける」という役割分担が、そのまま誤差カーブの形に現れています。緑の $H_v$ だけが両方の穴を避けている点に注目してください。
比較のため、$c$ を誤って幾何平均型を使った場合も見ておきましょう。
egC = db_err(HgC, fC)
print(f"[C] 幾何平均型Hv（cを誤った場合）RMS誤差 = {np.sqrt(np.mean(egC[band]**2)):.2f} dB, "
f"反共振で {egC[np.argmin(np.abs(fC-f_anti))]:+.2f} dB")
_, SffD = signal.welch(xD, fs=fs, nperseg=nper, noverlap=nper//2, window='hann')
i_r, i_b = np.argmin(np.abs(f-f_res[0])), np.argmin(np.abs(f-25.0))
print(f"[D] 加振力PSD: 共振1で25Hz比 {10*np.log10(SffD[i_r]/SffD[i_b]):+.1f} dB")
結果は [C] 幾何平均型Hv RMS誤差 = 4.84 dB, 反共振で +16.30 dB と [D] 加振力PSD: 共振1で25Hz比 -18.1 dB です。幾何平均型は $H_1$ と $H_2$ のちょうど中間を通るので、$H_2$ の $+36.8$ dB のちょうど半分の $+16.3$ dB になっています。「よく分からないから真ん中を取る」は、真ん中が正しい保証がなければ、単に誤差を半分にするだけだということです。$H_v$ の威力は $c$ の正しさに全面的に依存する、という当たり前ですが見落としやすい事実がここに出ています。実務では、加振していない状態（暗騒音測定）で各チャンネルのノイズフロアを測り、その比から $c$ を決めるのが現実的な手順になります。
オレンジの幾何平均型は、反共振でも高周波でも赤い $H_2$ のちょうど半分の高さをなぞっています。dBで真ん中を取る操作は誤差も半分にするだけで、真値（黒の0 dB線）に戻す力はありません。一方、緑の正しい $c$ の $H_v$ は反共振の $-3.4$ dB以外はほぼ0 dB線上にあり、高周波でも外れません。「分からないから中間」ではなく「暗騒音を測って $c$ を決める」ことが、$H_v$ を使う条件だと分かります。
恒等式と平均回数の確認
最後に、理論で導いた関係が数値でも厳密に成り立つことを確認します。
import numpy as np
b2 = (f > 1) & (f < 200)
print("γ² と |H1/H2| の最大差 :", np.max(np.abs(cohC[b2] - np.abs(H1C[b2]/H2C[b2]))))
print("arg H1 と arg H2 の最大差:", np.max(np.abs(np.angle(H1C[b2]/H2C[b2]))))
print("|H1|<=|Hv|<=|H2| の成立率:",
np.mean((np.abs(H1C[b2]) <= np.abs(HvC[b2])*(1+1e-9)) &
(np.abs(HvC[b2]) <= np.abs(H2C[b2])*(1+1e-9))))
print("|Hgeo| と sqrt(|H1||H2|) の最大相対差:",
np.max(np.abs(np.abs(HgC[b2]) - np.sqrt(np.abs(H1C[b2]*H2C[b2])))/np.abs(HgC[b2])))
rng = np.random.default_rng(11)
for n_d in [1, 2, 5, 10, 50, 200]:
a = rng.standard_normal(nper*n_d); bb = rng.standard_normal(nper*n_d)
kw = dict(fs=fs, nperseg=nper, noverlap=0, window='hann')
_, Saa = signal.welch(a, **kw); _, Sbb = signal.welch(bb, **kw)
_, Sab = signal.csd(a, bb, **kw)
g2 = np.abs(Sab)**2/(Saa*Sbb)
print(f"n_d={n_d:3d}: 無相関な2信号の平均γ² = {g2.mean():.4f}（理論 1/n_d = {1/n_d:.4f}）")
出力は次のとおりです。
γ² と |H1/H2| の最大差 : 5.55e-16
arg H1 と arg H2 の最大差: 1.99e-16
|H1|<=|Hv|<=|H2| の成立率: 1.0
|Hgeo| と sqrt(|H1||H2|) の最大相対差: 4.36e-16
n_d= 1: 無相関な2信号の平均γ² = 1.0000（理論 1/n_d = 1.0000）
n_d= 2: 無相関な2信号の平均γ² = 0.4843（理論 1/n_d = 0.5000）
n_d= 5: 無相関な2信号の平均γ² = 0.2009（理論 1/n_d = 0.2000）
n_d= 10: 無相関な2信号の平均γ² = 0.0976（理論 1/n_d = 0.1000）
n_d= 50: 無相関な2信号の平均γ² = 0.0197（理論 1/n_d = 0.0200）
n_d=200: 無相関な2信号の平均γ² = 0.0051（理論 1/n_d = 0.0050）
上4行はすべて倍精度の丸め誤差レベル（$10^{-16}$）で一致しており、$\gamma^2 = H_1/H_2$、位相の一致、$|H_1| \le |H_v| \le |H_2|$、幾何平均の関係が恒等式であることが確認できます。 これらはデータの性質によらず、定義から代数的に従う関係なので、実測データでも必ず成り立ちます。逆に言えば、解析ソフトの出力で $\gamma^2 \ne |H_1/H_2|$ になっていたら、規約の取り違えかバグを疑うべきです。
下6行は、コヒーレンスの上げ底の実演です。 $x$ と $y$ をまったく無関係な白色雑音にしているので、真のコヒーレンスは0のはずです。それでも $n_d = 1$ では $\gamma^2$ がぴったり $1.0000$、$n_d = 10$ でも平均 $0.098$ が出ます。「コヒーレンスが1に近いから測定はうまくいっている」という判断が、平均回数を確認しないと完全に無意味になることが分かります。$n_d = 200$ でようやく $0.005$ まで落ち、理論値 $1/n_d$ とよく一致しています。コヒーレンス曲線を見るときは、必ず平均回数とセットで読む——これは実験室で最も守られていないルールのひとつです。
左の対数グラフでは、実測した平均コヒーレンス（紫の点）が理論値 $1/n_d$（灰色の破線）にきれいに重なっています。$n_d = 1$ の点がちょうど1.0000に乗っていることが、「1ブロックではコヒーレンスが恒等的に1になる」ことの実証です。右の棒グラフは3つの恒等式の実測誤差で、いずれも倍精度の計算機イプシロン $2.2\times10^{-16}$ と同じ桁に収まっています。これらが近似ではなく恒等式であることが、数値からも裏づけられました。
理論と数値実験がそろいました。最後に、現場での判断手順としてまとめておきます。
実務での使い分け
ここまでの結果を、測定現場での意思決定の順に並べ替えます。まず、4つのケースでの成績を1枚にまとめておきます。
ケースA（出力雑音のみ）では $H_1$ が0.56 dBで $H_2$ の10.23 dBに圧勝し、ケースB（入力雑音のみ）では0.13 dB対0.84 dBで立場が逆転します。どちらか一方が常に優れているという事実はありません。 一方で緑の $H_v$ は、正しい $c$ を与えたケースC・Dで0.75 dB・0.13 dBと最良になっています（ケースA・Bの $H_v$ は $c$ を測っていない幾何平均型なので中間の成績です）。この棒グラフが「雑音の入り方を知らないまま推定量を固定するのは賭けである」ことを端的に示しています。
ステップ1: 雑音がどちら側に乗りやすいかを、物理から先に決める。 コヒーレンス曲線を眺めてから決めるのではなく、測定系の構成から予測します。加振器＋力センサ＋加速度センサの一般的な構成では、加速度センサ側のSN比が良く、力センサ側が苦しいことが多いため、共振点では $H_2$ が有利です。逆にインパルスハンマ加振では、ハンマの力信号は大きく（入力SNRが良い）、応答が減衰しきる前に窓で切るため出力側が苦しいので、$H_1$ が標準です。
ステップ2: 暗騒音を測ってノイズフロアを把握する。 加振していない状態で両チャンネルのPSDを取ります。これが $S_{uu}$ と $S_{vv}$ の実測値になり、そのまま $c^2 = S_{vv}/S_{uu}$ が求まります。$H_v$ を使うなら必須の手順ですし、使わなくても「どの周波数から先が測定限界か」が分かります。
ステップ3: $H_1$ と $H_2$ を両方プロットし、開き具合を見る。 2本の間隔が真値の存在区間です。開きが $1$ dB 以内の帯域は、どちらを使っても実用上同じ。開きが $6$ dB を超える帯域は、推定量の選択が結果を左右する危険地帯だと認識します。$\gamma^2$ を見るのと本質的には同じ情報ですが、dBスケールで直接「不確かさ幅」として見えるぶん判断しやすくなります。
ステップ4: コヒーレンスが低い帯域では、まず原因を切り分ける。 平均回数は足りているか（$n_d \ge 30$）、共振ピークで局所的にへこんでいないか（漏れ誤差）、入力レベルを変えると変化しないか（非線形性）。原因が測定雑音だと分かってから、初めて推定量の選択の話になります。
ステップ5: モーダルパラメータを同定するなら、帯域ごとに使い分ける。 共振ピークの高さから減衰比を読むなら、その帯域で入力SNRが良い推定量を。反共振から零点を読むなら、その帯域で出力SNRが良い推定量を。1本の周波数応答曲線を作るなら $H_v$ が便利ですが、$c$ の設定を間違えれば単に誤差を半分にするだけだという点は忘れないでください。
なお、多入力多出力（MIMO）系への拡張では、$H_1$ は多変数の正規方程式
$$ \bm{H}\, E\!\left[\bm{X}\bm{X}^H\right] = E\!\left[\bm{Y}\bm{X}^H\right] $$
を解く形に自然に一般化されます。左辺の $E[\bm{X}\bm{X}^H]$ は入力どうしの相互スペクトル行列で、これが正則である（＝複数の加振点の信号が互いに独立である）ことが解の存在条件です。多点加振で無相関な加振信号を使うのは、この行列を良条件に保つためです。一方、$H_2$ の一般化は入力数と出力数が等しい場合に限られるなど制約が強く、実用上は $\bm{H}_1$ か、$\bm{S}$ の固有分解に基づく $\bm{H}_v$ が使われます。SISOで身につけた「雑音がどちら側か」という視点は、MIMOでもそのまま通用します。
まとめ
本記事では、実測データから周波数応答を推定する2つの標準的な推定量、$H_1$ と $H_2$ について解説しました。
- 周波数応答の推定は割り算ではなく回帰である — $Y/X$ という素朴な式は「入力も出力も完璧に測れている」という仮定を隠し持っており、実測では成り立たない
- $H_1 = S_{xy}/S_{xx}$ — 出力にだけ雑音があると仮定して $\sum|Y_k – HX_k|^2$ を最小化した正規方程式の解。出力雑音に対して不偏
- $H_2 = S_{yy}/S_{yx}$ — 入力にだけ雑音があると仮定して $x$ 方向の残差を最小化し、逆数をとったもの。入力雑音に対して不偏
- バイアスの向きは逆 — $H_1 = H/(1+S_{uu}/S_{xx})$ で必ず過小、$H_2 = H(1+S_{vv}/(|H|^2S_{xx}))$ で必ず過大。真値は常に $|H_1|$ と $|H_2|$ の間にあり、位相は完全に一致する
- 崩れる周波数が違う — 力のドロップアウトで入力パワーが落ちる共振点では $H_1$ が、出力パワーがノイズフロアに沈む反共振点・高周波域では $H_2$ が崩れる。数値実験では共振で $H_1$ が $-6.6$ dB、反共振で $H_2$ が $+36$ dB のずれを示した
- $\gamma^2 = H_1/H_2$ — コヒーレンスは2つの推定量の比そのもの。$\gamma^2 = \gamma^2_{\text{in}}\gamma^2_{\text{out}}$ と分解でき、劣化の「量」は教えるが「どちら側か」は教えない
- $H_v$ は $H_1$ と $H_2$ を結ぶ族 — スケール $c^2 = S_{vv}/S_{uu}$ 付きのTLS解で、$c\to\infty$ で $H_1$、$c\to 0$ で $H_2$、$c^2 = S_{yy}/S_{xx}$ で幾何平均 $\sqrt{H_1H_2}$。正しい $c$ を与えれば両側に雑音があっても不偏
- 平均回数を必ず確認する — $n_d=1$ ではコヒーレンスは恒等的に1になり、無相関な信号でも $E[\hat\gamma^2] \approx 1/n_d$ の上げ底が乗る
この枠組みは、振動試験に限らず、スピーカと室内音響の伝達特性測定、回路の周波数特性測定、制御対象のプラント同定など、「入出力を測って伝達関数を推定する」あらゆる場面に同じ形で現れます。次にコヒーレンスの低い測定データに出会ったら、$H_1$ と $H_2$ を両方描いて、その開きの中に真値を探してみてください。
次のステップとして、以下の記事も参考にしてください。