MMSE等化器とゼロフォーシング等化器の理論と導出と実装

スマートフォンの動画ストリーミングが、ビルの谷間や駅の構内でも途切れずに届くのはなぜでしょうか。電波は送信アンテナから受信機まで、まっすぐ一本道を通ってくるわけではありません。建物や地面で反射した何本もの「こだま」が、わずかにずれた時刻に重なって届きます。すると、いま受け取ろうとしているシンボル(情報の最小単位)の上に、一つ前・二つ前に送られたシンボルの尾が乗ってきます。これがシンボル間干渉(ISI: Inter-Symbol Interference)です。何も対策をしなければ、受信した波形はぐちゃぐちゃに混ざり合い、送られたビットを読み取れなくなります。

この混ざり合いをほどき、送信シンボルを復元する受信側の処理が等化器(equalizer)です。等化器は「チャネルが信号にかけた歪みの逆操作」を施す装置だと考えると分かりやすいでしょう。本記事では、もっとも基本的な2つの線形等化器を扱います。チャネルの逆数をそのまま当てるゼロフォーシング(ZF: Zero-Forcing)等化器と、雑音まで考慮して平均二乗誤差を最小化するMMSE(Minimum Mean Square Error)等化器です。

この2つは、応用の場面が驚くほど広い概念です。LTEや5GのようなOFDM通信では各サブキャリアごとにZF/MMSE等化が行われ、有線のADSLやEthernetのケーブル伝送でも同じ原理が使われます。さらに、等化器の背骨にある「直交原理」は、画像のデブラー(ぼけ除去)やレーダーのパルス圧縮、株価予測のウィーナーフィルタにまで通じる、信号処理の共通言語です。ZFとMMSEの違いを一度理解すれば、これらすべてを同じ目で眺められるようになります。

本記事の内容

  • 受信信号モデル $\bm{r} = \bm{H}\bm{s} + \bm{n}$ とISIの正体
  • ゼロフォーシング等化器 $\bm{W} = (\bm{H}^H \bm{H})^{-1}\bm{H}^H$ の導出
  • 直交原理(直交性条件 $\mathbb{E}[(\bm{s} – \bm{W}\bm{r})\bm{r}^H] = \bm{0}$)からのMMSE等化器の導出
  • MMSE等化器 $\bm{W} = (\bm{H}^H \bm{H} + \frac{\sigma_n^2}{\sigma_s^2}\bm{I})^{-1}\bm{H}^H$ の意味
  • ZFの「雑音増幅」とMMSEの「SNR次第の優位性」の比較
  • ASK/QAM信号を多重路チャネルに通したPython実装(星座図とMSE-SNR曲線)

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

これらを読んでいなくても本記事は独立して読めるよう、必要な道具立ては都度説明します。ただし、複素ベクトル・行列の演算(特にエルミート転置 $\bm{A}^H$)と、QAMが複素平面上の点としてシンボルを表すことだけは、頭の片隅に置いておいてください。

等化器とは — チャネルの歪みを「巻き戻す」装置

まず、何を相手にしているのかをはっきりさせましょう。送信機は離散的なシンボル列 $s[0], s[1], s[2], \dots$ を送り出します。QAMなら、各 $s[k]$ は複素平面上の格子点(コンステレーション点)です。理想的なチャネルなら、受信機にはこのシンボルがそのまま、振幅と位相を保ったまま届きます。

ところが現実のチャネルは、いくつもの遅延した複製を作り出します。直接波が時刻ぴったりに届く一方、ビルで反射した波は1シンボル分遅れ、地面で反射した波は2シンボル分遅れて届く、といった具合です。これを数式で書くと、チャネルはインパルス応答 $h[0], h[1], \dots, h[L]$ を持つフィルタとしてふるまい、受信信号は送信シンボル列とこのインパルス応答の畳み込みになります。

$$ r[k] = \sum_{\ell=0}^{L} h[\ell]\, s[k-\ell] + n[k] $$

ここで $n[k]$ は受信機の熱雑音です。右辺の和を $\ell=0$ の項だけ取り出すと $h[0]\,s[k]$、つまり「いま欲しいシンボル」ですが、$\ell \geq 1$ の項 $h[1]s[k-1] + h[2]s[k-2] + \cdots$ が過去のシンボルの混入=ISIです。等化器の仕事は、この $r[k]$ から $s[k]$ を可能なかぎり正確に取り出すことに尽きます。

直感的なイメージとしては、チャネルが信号を「ぼかすレンズ」だとすれば、等化器は「ピントを合わせ直すレンズ」です。チャネルが $h$ という係数で混ぜたのなら、その逆の操作 $h^{-1}$ をかければ元に戻る——これがゼロフォーシングの素朴なアイデアです。ただし後で見るように、この素朴な逆操作には落とし穴があります。雑音まで一緒に拡大してしまうのです。その落とし穴をエレガントに回避するのがMMSE等化器です。全体像を一枚の図にまとめておきます。

等化器の概念図:多重路チャネルが混ぜたシンボルを等化器がほどく

この図は本記事の登場人物を一望したものです。送信シンボル $\bm{s}$ は「ぼかすレンズ」であるチャネルを通って ISI と雑音が乗った $\bm{r}=\bm{H}\bm{s}+\bm{n}$ になり、「ピントを戻すレンズ」である等化器が $\hat{\bm{s}}=\bm{W}\bm{r}$ で復元します。ZF は $\bm{H}$ を完全に逆転させる方針、MMSE は雑音まで考慮して逆転の強さを加減する方針だ、という違いだけを押さえておけば、以降の導出はすべてこの絵の上に乗っています。

次のセクションでは、この畳み込みを行列の形にまとめ直します。行列で書くと、ZFもMMSEも「連立一次方程式をどう解くか」という見慣れた問題に姿を変えます。

受信信号の行列モデル

畳み込みの式は便利ですが、最適化の道具を当てるには行列・ベクトルの形にするのが見通しがよくなります。送信したい $N$ 個のシンボルを縦に並べたベクトルを

$$ \bm{s} = \begin{bmatrix} s[0] \\ s[1] \\ \vdots \\ s[N-1] \end{bmatrix} $$

とし、受信サンプルを同様に並べたものを $\bm{r}$、雑音を $\bm{n}$ とします。畳み込み $r[k] = \sum_\ell h[\ell] s[k-\ell] + n[k]$ は、係数 $h[\ell]$ を斜めに並べた畳み込み行列(チャネル行列) $\bm{H}$ を使って

$$ \begin{equation} \bm{r} = \bm{H}\bm{s} + \bm{n} \end{equation} $$

と書けます。$\bm{H}$ は、各行が一つずつずれて $h$ の係数を並べた構造(Toeplitz型)をしています。たとえば $L=1$(2タップのチャネル $h[0], h[1]$)で $N=3$ なら

$$ \bm{H} = \begin{bmatrix} h[0] & 0 & 0 \\ h[1] & h[0] & 0 \\ 0 & h[1] & h[0] \\ 0 & 0 & h[1] \end{bmatrix} $$

のようになります。各列が「あるシンボルがどの受信サンプルに、どれだけの係数で散らばるか」を表しているわけです。行が列より多いのは、最後のシンボルの尾が信号区間の終端より後ろにはみ出すためです。

このモデルで重要なのは、$\bm{H}$ がシンボルどうしを混ぜる行列だということです。$\bm{H}$ が対角行列なら混ざりはなく、各シンボルは独立に受信されます。非対角成分が ISI そのものです。したがって等化とは、「$\bm{H}$ を打ち消すような行列 $\bm{W}$ を左からかけて、$\bm{W}\bm{r}$ をできるだけ $\bm{s}$ に近づける」線形操作だと定式化できます。

$$ \hat{\bm{s}} = \bm{W}\bm{r} $$

統計的な前提も置いておきます。送信シンボルは平均ゼロ・分散 $\sigma_s^2$ で互いに無相関、雑音も平均ゼロ・分散 $\sigma_n^2$ の白色雑音で、シンボルとは無相関とします。式で書くと

$$ \mathbb{E}[\bm{s}\bm{s}^H] = \sigma_s^2 \bm{I}, \quad \mathbb{E}[\bm{n}\bm{n}^H] = \sigma_n^2 \bm{I}, \quad \mathbb{E}[\bm{s}\bm{n}^H] = \bm{0} $$

です。これらは後の期待値計算で何度も使う「材料」になります。$\bm{I}$ は単位行列で、無相関性が対角行列として表れている点に注目してください。

ここまでで問題は「$\bm{W}\bm{r}$ を $\bm{s}$ に近づける $\bm{W}$ を求めよ」という形に整理できました。近づけ方には流儀があります。まず雑音を無視して「$\bm{H}$ を完璧に打ち消す」流儀=ゼロフォーシングから見ていきましょう。

ゼロフォーシング等化器の導出

ゼロフォーシングの発想はきわめて素直です。「ISIをゼロに強制(force to zero)する」。つまり雑音 $\bm{n}$ がないと仮定して、$\bm{W}\bm{r} = \bm{W}\bm{H}\bm{s}$ がちょうど $\bm{s}$ に一致するような $\bm{W}$ を探します。理想は

$$ \bm{W}\bm{H} = \bm{I} $$

です。これが成り立てば、$\hat{\bm{s}} = \bm{W}\bm{H}\bm{s} = \bm{s}$ となり、ISIは完全に消えます。

ここで問題になるのが $\bm{H}$ の形です。先ほど見たように $\bm{H}$ は縦長(行数 $>$ 列数)の長方形行列で、ふつうの逆行列を持ちません。そこで最小二乗解の考え方を使います。雑音がないモデル $\bm{r} = \bm{H}\bm{s}$ を $\bm{s}$ について解きたいのですが、行が多すぎて厳密には解けない(過剰決定)ので、残差 $\|\bm{r} – \bm{H}\bm{s}\|^2$ を最小にする $\bm{s}$ を推定値とします。この最小二乗問題の解は、正規方程式

$$ \bm{H}^H \bm{H}\, \hat{\bm{s}} = \bm{H}^H \bm{r} $$

を満たします。導出を省略せず追ってみましょう。最小化したい目的関数を

$$ J(\hat{\bm{s}}) = \|\bm{r} – \bm{H}\hat{\bm{s}}\|^2 = (\bm{r} – \bm{H}\hat{\bm{s}})^H (\bm{r} – \bm{H}\hat{\bm{s}}) $$

と書きます。これを展開すると

$$ J(\hat{\bm{s}}) = \bm{r}^H\bm{r} – \bm{r}^H\bm{H}\hat{\bm{s}} – \hat{\bm{s}}^H\bm{H}^H\bm{r} + \hat{\bm{s}}^H\bm{H}^H\bm{H}\hat{\bm{s}} $$

となります。複素ベクトル $\hat{\bm{s}}$ で微分してゼロと置くのですが、ここでウィルティンガー微分の便利な公式 $\frac{\partial}{\partial \hat{\bm{s}}^H}(\hat{\bm{s}}^H \bm{A}\hat{\bm{s}}) = \bm{A}\hat{\bm{s}}$ と $\frac{\partial}{\partial \hat{\bm{s}}^H}(\hat{\bm{s}}^H\bm{b}) = \bm{b}$ を使います。$\hat{\bm{s}}^H$ について偏微分すると、$\hat{\bm{s}}$ を含まない第1項は消え、第2項も $\hat{\bm{s}}^H$ を含まないので消え、残る第3・第4項から

$$ \frac{\partial J}{\partial \hat{\bm{s}}^H} = -\bm{H}^H\bm{r} + \bm{H}^H\bm{H}\hat{\bm{s}} = \bm{0} $$

が得られます。これを整理すれば、先ほどの正規方程式 $\bm{H}^H\bm{H}\hat{\bm{s}} = \bm{H}^H\bm{r}$ そのものです。$\bm{H}^H\bm{H}$ が正則なら両辺に逆行列をかけて

$$ \hat{\bm{s}} = (\bm{H}^H\bm{H})^{-1}\bm{H}^H \bm{r} $$

を得ます。したがってゼロフォーシング等化器のタップ係数行列は

$$ \begin{equation} \bm{W}_{\mathrm{ZF}} = (\bm{H}^H\bm{H})^{-1}\bm{H}^H \end{equation} $$

です。この $(\bm{H}^H\bm{H})^{-1}\bm{H}^H$ は $\bm{H}$ のムーア・ペンローズ擬似逆行列 $\bm{H}^+$ そのものであり、確かに $\bm{W}_{\mathrm{ZF}}\bm{H} = \bm{H}^+\bm{H} = \bm{I}$ を満たします。ISIは完璧に消えるのです。

ZFの落とし穴 — 雑音増幅

では何が問題なのでしょうか。雑音を戻して、ZF等化後の信号を見てみます。

$$ \hat{\bm{s}}_{\mathrm{ZF}} = \bm{W}_{\mathrm{ZF}}\bm{r} = \bm{W}_{\mathrm{ZF}}(\bm{H}\bm{s} + \bm{n}) = \underbrace{\bm{s}}_{\text{ISIは消えた}} + \underbrace{\bm{W}_{\mathrm{ZF}}\bm{n}}_{\text{雑音が変形される}} $$

第1項は完璧に元のシンボルですが、第2項で雑音に $\bm{W}_{\mathrm{ZF}}$ がかかっています。問題はこの $\bm{W}_{\mathrm{ZF}}$ の大きさです。周波数領域で考えると見通しがよくなります。チャネルの周波数応答 $H(e^{j\omega})$ がある周波数で深く落ち込む「ノッチ」を持つとき、ZFはその逆数 $1/H(e^{j\omega})$ をかけるので、ノッチの周波数で雑音を爆発的に増幅してしまいます。チャネルが弱いところほど、無理に持ち上げようとして雑音まで一緒に持ち上げてしまうわけです。

行列の言葉では、$\bm{H}^H\bm{H}$ の最小固有値が小さい(=チャネルが特定方向で弱い)と、その逆行列の最大固有値が大きくなり、雑音電力 $\sigma_n^2\, \mathrm{tr}\!\left[(\bm{H}^H\bm{H})^{-1}\right]$ が膨れ上がります。ISIをゼロにする代償に、SN比を悪化させてしまう——これがZFの本質的な限界です。

この雑音増幅が ISI の強さとともにどう変わるかを、実際に測ってみましょう。2タップチャネル $h=[1,\, h_1]$ の $h_1$ を $0$ から $0.95$ まで動かし、ZF と MMSE(SNR=10 dB)の正規化雑音増幅量 $\mathrm{tr}(\bm{W}\bm{W}^H)/N$ をプロットしたのが次の図です。

ISI強度に対するZFとMMSEの雑音増幅量

縦軸は対数です。ISIが弱い($h_1$ が小さい)うちは両者とも基準(=1)付近にいますが、$h_1$ が $1$ に近づきチャネルが弱い方向を持つほど、ZF(赤)の雑音増幅量は急激に跳ね上がります。$h_1=0.9$ では ZF が約 2.9 倍まで膨れる一方、MMSE(青)は約 0.84 と基準を下回ったまま抑えられています。MMSE は強ISI領域で「無理に持ち上げない」ことで雑音を頭打ちにする、という理論的主張がそのまま数値に表れています。

ここで自然な疑問が生まれます。ISIを完璧にゼロにすることに、本当にこだわるべきなのでしょうか。少しだけISIを残してでも、雑音増幅を抑えたほうが、トータルの誤差は小さくなるのではないか——この発想がMMSE等化器につながります。

MMSE等化器の導出 — 直交原理から

MMSEの目標は明確です。ISIと雑音を別々に扱うのではなく、両方ひっくるめた推定誤差の平均二乗 $\mathbb{E}[\|\bm{s} – \bm{W}\bm{r}\|^2]$ を最小にする $\bm{W}$ を求めます。誤差ベクトルを

$$ \bm{e} = \bm{s} – \hat{\bm{s}} = \bm{s} – \bm{W}\bm{r} $$

と定義し、目的関数を

$$ J_{\mathrm{MMSE}}(\bm{W}) = \mathbb{E}[\|\bm{e}\|^2] = \mathbb{E}\!\left[(\bm{s} – \bm{W}\bm{r})^H(\bm{s} – \bm{W}\bm{r})\right] $$

とします。これを直接微分してもよいのですが、もっと美しい道具があります。直交原理(orthogonality principle)です。

直交原理とは

直交原理は、推定問題における幾何学的な事実です。「平均二乗誤差を最小にする線形推定量では、推定誤差は観測データと直交する」というものです。直交とは、ここでは相関がゼロ、つまり期待値 $\mathbb{E}[\bm{e}\bm{r}^H] = \bm{0}$ を意味します。

直感的なアナロジーで説明しましょう。3次元空間の点 $\bm{s}$ を、ある平面(観測データ $\bm{r}$ が張る空間)の上に「影」として投影することを考えてください。$\bm{s}$ にもっとも近い平面上の点は、$\bm{s}$ から平面に下ろした垂線の足です。このとき、$\bm{s}$ と垂線の足を結ぶ誤差ベクトルは、平面に対して垂直=直交します。もし誤差が平面に対して斜めなら、影をその方向に少しずらせばもっと近づけられるので、まだ最適ではありません。誤差が観測の張る空間と直交していること——これが「これ以上近づけない最適解」の合図なのです。MMSE推定はまさにこの幾何学的投影で、観測 $\bm{r}$ の張る空間への射影が最良の線形推定になります。

これを式で確かめます。最適性条件は、目的関数の $\bm{W}$ に関する微分がゼロになることです。

$$ \frac{\partial}{\partial \bm{W}^*} \mathbb{E}\!\left[(\bm{s} – \bm{W}\bm{r})^H(\bm{s} – \bm{W}\bm{r})\right] = -\mathbb{E}\!\left[(\bm{s} – \bm{W}\bm{r})\bm{r}^H\right] = \bm{0} $$

これより、直交性条件

$$ \begin{equation} \mathbb{E}\!\left[(\bm{s} – \bm{W}\bm{r})\bm{r}^H\right] = \bm{0} \end{equation} $$

が得られます。つまり「最適 $\bm{W}$ のもとでは、誤差 $\bm{e} = \bm{s} – \bm{W}\bm{r}$ が観測 $\bm{r}$ と無相関になる」。これが直交原理です。微分の細部に立ち入らずとも、幾何学的に「誤差は観測と直交する」と理解しておけば十分です。

直交性条件を解く

直交性条件をほどいていきます。期待値の中で $\bm{W}$ は定数なので外に出せます。

$$ \mathbb{E}[\bm{s}\bm{r}^H] – \bm{W}\,\mathbb{E}[\bm{r}\bm{r}^H] = \bm{0} $$

ここで2つの相関行列が現れました。$\bm{R}_{sr} = \mathbb{E}[\bm{s}\bm{r}^H]$ は送信シンボルと受信信号の相互相関、$\bm{R}_{rr} = \mathbb{E}[\bm{r}\bm{r}^H]$ は受信信号の自己相関です。これらを使えば最適等化器は

$$ \bm{W}_{\mathrm{MMSE}} = \bm{R}_{sr}\, \bm{R}_{rr}^{-1} $$

と、ウィーナー解の一般形で書けます。あとはこの2つの相関行列を、信号モデル $\bm{r} = \bm{H}\bm{s} + \bm{n}$ から具体的に計算するだけです。

まず相互相関 $\bm{R}_{sr}$ を求めます。$\bm{r}^H = (\bm{H}\bm{s} + \bm{n})^H = \bm{s}^H\bm{H}^H + \bm{n}^H$ を代入し、期待値の線形性で項ごとに分けると

$$ \bm{R}_{sr} = \mathbb{E}[\bm{s}(\bm{s}^H\bm{H}^H + \bm{n}^H)] = \mathbb{E}[\bm{s}\bm{s}^H]\bm{H}^H + \mathbb{E}[\bm{s}\bm{n}^H] $$

となります。ここで前提の $\mathbb{E}[\bm{s}\bm{s}^H] = \sigma_s^2\bm{I}$ と $\mathbb{E}[\bm{s}\bm{n}^H] = \bm{0}$(シンボルと雑音は無相関)を代入すると、第2項が消えて

$$ \bm{R}_{sr} = \sigma_s^2\, \bm{H}^H $$

が得られます。次に自己相関 $\bm{R}_{rr}$ を求めます。$\bm{r} = \bm{H}\bm{s} + \bm{n}$ を代入して展開すると

$$ \bm{R}_{rr} = \mathbb{E}[(\bm{H}\bm{s} + \bm{n})(\bm{H}\bm{s} + \bm{n})^H] = \mathbb{E}[\bm{H}\bm{s}\bm{s}^H\bm{H}^H + \bm{H}\bm{s}\bm{n}^H + \bm{n}\bm{s}^H\bm{H}^H + \bm{n}\bm{n}^H] $$

となります。クロス項 $\mathbb{E}[\bm{H}\bm{s}\bm{n}^H]$ と $\mathbb{E}[\bm{n}\bm{s}^H\bm{H}^H]$ は、シンボルと雑音の無相関性 $\mathbb{E}[\bm{s}\bm{n}^H] = \bm{0}$ によりどちらもゼロになります。残った項に $\mathbb{E}[\bm{s}\bm{s}^H] = \sigma_s^2\bm{I}$ と $\mathbb{E}[\bm{n}\bm{n}^H] = \sigma_n^2\bm{I}$ を代入すると

$$ \bm{R}_{rr} = \sigma_s^2\, \bm{H}\bm{H}^H + \sigma_n^2\, \bm{I} $$

になります。この2つを $\bm{W}_{\mathrm{MMSE}} = \bm{R}_{sr}\bm{R}_{rr}^{-1}$ に代入すると

$$ \bm{W}_{\mathrm{MMSE}} = \sigma_s^2\, \bm{H}^H \left(\sigma_s^2\, \bm{H}\bm{H}^H + \sigma_n^2\, \bm{I}\right)^{-1} $$

を得ます。ここで全体を $\sigma_s^2$ でくくり、$\sigma_s^2$ で割る形に整理すると

$$ \bm{W}_{\mathrm{MMSE}} = \bm{H}^H \left(\bm{H}\bm{H}^H + \frac{\sigma_n^2}{\sigma_s^2}\, \bm{I}\right)^{-1} $$

となります。

2つの等価な表現とその意味

上の式は「受信側に逆行列を置く」形ですが、行列の恒等式(push-through identity)$\bm{H}^H(\bm{H}\bm{H}^H + c\bm{I})^{-1} = (\bm{H}^H\bm{H} + c\bm{I})^{-1}\bm{H}^H$ を使うと、課題で指定された「送信側に逆行列を置く」形に書き換えられます。

$$ \begin{equation} \bm{W}_{\mathrm{MMSE}} = \left(\bm{H}^H\bm{H} + \frac{\sigma_n^2}{\sigma_s^2}\, \bm{I}\right)^{-1}\bm{H}^H \end{equation} $$

この恒等式は、両辺に左から $(\bm{H}^H\bm{H} + c\bm{I})$、右から $(\bm{H}\bm{H}^H + c\bm{I})$ をかけて $(\bm{H}^H\bm{H} + c\bm{I})\bm{H}^H = \bm{H}^H(\bm{H}\bm{H}^H + c\bm{I})$ が両辺とも $\bm{H}^H\bm{H}\bm{H}^H + c\bm{H}^H$ に等しいことから確かめられます。どちらの形を使っても結果は同じです。

この式 (5) をZF等化器 (2) の $\bm{W}_{\mathrm{ZF}} = (\bm{H}^H\bm{H})^{-1}\bm{H}^H$ と並べてみてください。違いはただ一箇所、逆行列の中の $\frac{\sigma_n^2}{\sigma_s^2}\bm{I}$ という対角の上乗せ項だけです。この項こそがMMSEの知恵です。

  • 高SNR(雑音が小さい)の極限: $\sigma_n^2 \to 0$ なら $\frac{\sigma_n^2}{\sigma_s^2} \to 0$ となり、MMSEはZFに一致します。雑音がほぼないなら、ISIを完璧に消すZFがそのまま最適だからです。
  • 低SNR(雑音が大きい): $\frac{\sigma_n^2}{\sigma_s^2}$ が大きくなると、対角項が $\bm{H}^H\bm{H}$ の小さな固有値を底上げします。これにより逆行列が暴れず、雑音増幅が抑えられます。

つまりこの上乗せ項は、$\bm{H}^H\bm{H}$ の固有値が小さい(チャネルが弱い)方向で「無理に持ち上げない」ためのブレーキです。最適化や機械学習に詳しい方なら、これがリッジ回帰(Tikhonov正則化)と同じ形だと気づくでしょう。実際、MMSE等化器は「逆行列を正則化することで雑音に頑健にしたZF」と見ることができます。逆数の比 $\frac{\sigma_s^2}{\sigma_n^2}$ が信号対雑音電力比(SNR)に対応するので、SNRが正則化の強さを自動調整しているのです。

ここまでで、ZFとMMSEのタップ係数を導出し、両者が「正則化項の有無」だけで違うことが分かりました。次は、$2 \times 2$ の最小例で実際に数値を入れて、違いを手で確かめてみましょう。

具体例 — 2タップチャネルで手計算する

抽象的な行列のままだとピンと来にくいので、できるだけ小さな例で数値を入れてみます。送信シンボルが2個 $\bm{s} = [s_0, s_1]^T$、チャネル行列を実数の $2\times 2$ で

$$ \bm{H} = \begin{bmatrix} 1.0 & 0.0 \\ 0.9 & 1.0 \end{bmatrix} $$

とします。これは「2番目の受信サンプルに、1番目のシンボルが $0.9$ 倍の強さで漏れ込む」強いISIのある状況です。信号分散 $\sigma_s^2 = 1$、雑音分散 $\sigma_n^2 = 0.1$(つまりSNRは10倍=10 dB)としましょう。

まず $\bm{H}^H\bm{H}$(実行列なので転置)を計算します。

$$ \bm{H}^T\bm{H} = \begin{bmatrix} 1.0 & 0.9 \\ 0.0 & 1.0 \end{bmatrix}\begin{bmatrix} 1.0 & 0.0 \\ 0.9 & 1.0 \end{bmatrix} = \begin{bmatrix} 1.81 & 0.9 \\ 0.9 & 1.0 \end{bmatrix} $$

ZF等化器は $(\bm{H}^T\bm{H})^{-1}\bm{H}^T$ です。$\bm{H}^T\bm{H}$ の行列式は $1.81 \times 1.0 – 0.9 \times 0.9 = 1.81 – 0.81 = 1.0$ なので、逆行列は

$$ (\bm{H}^T\bm{H})^{-1} = \begin{bmatrix} 1.0 & -0.9 \\ -0.9 & 1.81 \end{bmatrix} $$

です。対角に $1.81$ という大きな値が出ている点に注目してください。これがそのまま雑音増幅につながります。実際、ZF等化後の雑音電力は $\sigma_n^2\,\mathrm{tr}[(\bm{H}^T\bm{H})^{-1}] = 0.1 \times (1.0 + 1.81) = 0.281$ で、元の雑音電力 $0.1$ の約2.8倍に膨らんでいます。

一方MMSEは、逆行列の中身が $\bm{H}^T\bm{H} + \frac{\sigma_n^2}{\sigma_s^2}\bm{I} = \bm{H}^T\bm{H} + 0.1\,\bm{I}$ になります。

$$ \bm{H}^T\bm{H} + 0.1\bm{I} = \begin{bmatrix} 1.91 & 0.9 \\ 0.9 & 1.1 \end{bmatrix} $$

対角に $0.1$ ずつ足されたことで、行列がより「丸く」なり、逆行列の暴れが抑えられます。行列式は $1.91 \times 1.1 – 0.81 = 2.101 – 0.81 = 1.291$ と大きくなり、逆行列の各成分は小さくなります。これがMMSEが雑音を抑えるからくりです。完璧にISIをゼロにする代わりにわずかな残留ISIを許す——その代わり雑音をぐっと抑える、という取引をしているのです。

この $2\times 2$ の例だけでも、「正則化項を足すと逆行列が穏やかになる」という核心が見えました。とはいえ手計算では雑音の統計的なふるまいまでは追えません。次のセクションでPythonに渡し、実際のQAM信号を多数のシンボルで流して、星座図とMSE-SNR曲線で両者の違いを体感しましょう。

Pythonでの実装

ここからは、これまでの理論をコードで確かめます。まず使い回す道具立て(チャネル、QAMシンボル生成、等化器、可視化用フォント)を準備します。

セットアップとチャネルの周波数特性

最初に、日本語フォントの設定と、本記事で扱う多重路チャネルの周波数応答を可視化します。ZFが雑音を増幅する原因は「チャネルのノッチ」なので、まずどこに弱点があるかを目で見ておきます。

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 == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

# 多重路チャネルのインパルス応答(直接波+反射波2本)
h = np.array([1.0, 0.0, 0.7, 0.0, 0.4], dtype=complex)  # h[0],h[2],h[4]に成分

# チャネルの周波数応答を計算
w = np.linspace(0, np.pi, 512)
H_freq = np.array([np.sum(h * np.exp(-1j * w_k * np.arange(len(h)))) for w_k in w])

plt.figure(figsize=(8, 4.5))
plt.plot(w / np.pi, 20 * np.log10(np.abs(H_freq) + 1e-12), color="#1f77b4", lw=2)
plt.xlabel("正規化周波数 $\\omega/\\pi$")
plt.ylabel("振幅応答 [dB]")
plt.title("多重路チャネルの周波数応答(深いノッチがISIの正体)")
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig("/tmp/mmse_channel_freq.png", dpi=130)
plt.show()

多重路チャネルの周波数応答と深いノッチ

この図から、チャネルの振幅応答が平坦ではなく、特定の周波数で深く落ち込む「ノッチ」を持つことが読み取れます。落ち込んだ周波数では、ZF等化器がその逆数(大きな値)をかけて信号を持ち上げようとするため、同じ周波数にある雑音まで一緒に増幅されます。ノッチが深いほど雑音増幅は激しくなり、これがZFの泣きどころになります。MMSEはこのノッチで「無理に持ち上げない」ことで雑音を抑えます。

畳み込みチャネル行列と等化器の構築

次に、インパルス応答 $h$ から畳み込み行列 $\bm{H}$ を組み立て、ZFとMMSEのタップ係数を計算する関数を用意します。理論式 (2) と (5) をそのままコードに落とします。

import numpy as np

def build_channel_matrix(h, N):
    """インパルス応答 h と送信シンボル数 N から畳み込み行列 H を作る"""
    L = len(h) - 1          # チャネルのメモリ長
    M = N + L               # 受信サンプル数(尾がはみ出す分だけ長い)
    H = np.zeros((M, N), dtype=complex)
    for col in range(N):    # 各列は1シンボルが散らばる様子
        H[col:col + len(h), col] = h
    return H

def zf_equalizer(H):
    """ZF等化器 W = (H^H H)^{-1} H^H (擬似逆行列)"""
    HhH = H.conj().T @ H
    return np.linalg.inv(HhH) @ H.conj().T

def mmse_equalizer(H, snr_linear):
    """MMSE等化器 W = (H^H H + (sigma_n^2/sigma_s^2) I)^{-1} H^H"""
    N = H.shape[1]
    reg = 1.0 / snr_linear            # sigma_n^2/sigma_s^2 = 1/SNR (sigma_s^2=1)
    HhH = H.conj().T @ H
    return np.linalg.inv(HhH + reg * np.eye(N)) @ H.conj().T

# 動作確認: 強いISIを持つ2タップチャネルで雑音増幅を見る
H_demo = build_channel_matrix(np.array([1.0, 0.9], dtype=complex), N=8)
W_zf = zf_equalizer(H_demo)
noise_gain = np.trace((W_zf @ W_zf.conj().T).real)
print(f"ZF等化器の雑音増幅量 tr(W W^H) = {noise_gain:.3f}")
print(f"等化なしの基準(=N=8)と比べて {noise_gain/8:.2f} 倍")

このコードを実行すると、ZF等化器の雑音増幅量 $\mathrm{tr}(\bm{W}\bm{W}^H)$ が出力されます。等化を一切しない場合の雑音電力は $N=8$ に相当しますが、ZFではこれより大きな値(強ISIのため約14、N=8の約1.75倍)になり、ISIを消す代償として雑音電力が膨らむことが数値で確認できます。build_channel_matrix は理論で見た「列ごとに $h$ を一つずつずらして並べる」Toeplitz構造をそのまま作っており、zf_equalizermmse_equalizer は導出した式 (2)・(5) の素直な実装です。

QAMシンボルを通して星座図を比較

いよいよ本番です。16-QAMのシンボル列を多重路チャネルに通し、雑音を加えてから、ZFとMMSEで等化した結果を星座図(コンステレーション)で比較します。星座図とは、複素平面上に受信シンボルを点として打ったもので、点が格子状にきれいに分かれていれば正しく復調できる、という見方をします。

import numpy as np

rng = np.random.default_rng(0)

def qam16_symbols(n, rng):
    """16-QAMシンボルを生成(正規化して平均電力1)"""
    levels = np.array([-3, -1, 1, 3])
    re = rng.choice(levels, size=n)
    im = rng.choice(levels, size=n)
    s = (re + 1j * im) / np.sqrt(10.0)   # 平均電力が1になるよう正規化
    return s

# パラメータ
h = np.array([1.0, 0.0, 0.7, 0.0, 0.4], dtype=complex)
N = 2000                  # シンボル数
snr_db = 18.0             # SNR [dB]
snr_lin = 10 ** (snr_db / 10)

# 送信・チャネル通過・雑音付加
s = qam16_symbols(N, rng)
H = build_channel_matrix(h, N)
r_clean = H @ s
sigma_n2 = 1.0 / snr_lin   # sigma_s^2=1 なので雑音分散=1/SNR
noise = np.sqrt(sigma_n2 / 2) * (rng.standard_normal(len(r_clean))
                                 + 1j * rng.standard_normal(len(r_clean)))
r = r_clean + noise

# 等化
s_zf = zf_equalizer(H) @ r
s_mmse = mmse_equalizer(H, snr_lin) @ r
print(f"ZF   MSE = {np.mean(np.abs(s - s_zf) ** 2):.4f}")
print(f"MMSE MSE = {np.mean(np.abs(s - s_mmse) ** 2):.4f}")

実行すると、ZFとMMSEそれぞれの平均二乗誤差(MSE)が出力されます。SNR=18 dBという比較的良好な条件では両者の差はまだ小さい(ZF約0.025、MMSE約0.024)ものの、MMSEが一貫してZF以下のMSEを返すことが数値で確認できます。この差はSNRが下がるほど大きく開きます(後のMSE-SNR曲線で確認します)。次に、この等化結果を星座図で可視化します。

import numpy as np
import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 3, figsize=(15, 5))
titles = ["等化なし(チャネル通過後)", "ZF等化(ISIは消えるが雑音増幅)",
          "MMSE等化(雑音まで考慮)"]
data = [r[:N], s_zf, s_mmse]
colors = ["#888888", "#d62728", "#1f77b4"]

for ax, title, d, c in zip(axes, titles, data, colors):
    ax.scatter(d.real, d.imag, s=6, alpha=0.4, color=c)
    # 理想シンボル位置を重ねる
    lv = np.array([-3, -1, 1, 3]) / np.sqrt(10.0)
    gx, gy = np.meshgrid(lv, lv)
    ax.scatter(gx, gy, marker="x", s=80, color="black", linewidths=1.5)
    ax.set_title(title)
    ax.set_xlabel("同相成分 I")
    ax.set_ylabel("直交成分 Q")
    ax.set_xlim(-1.6, 1.6); ax.set_ylim(-1.6, 1.6)
    ax.set_aspect("equal"); ax.grid(alpha=0.3)

plt.tight_layout()
plt.savefig("/tmp/mmse_constellation.png", dpi=130)
plt.show()

ZFとMMSEの星座図比較(16-QAM, SNR=18 dB)

3枚の星座図を左から見比べてください。左の「等化なし」では、ISIのせいでシンボルが格子点(黒い×印)からずれ、雲のように広がってまったく分離できていません。中央のZF等化では16個の格子状の塊が現れ、ISIが消えたことが分かりますが、各塊にはノッチ周波数で増幅された雑音による広がりが残っています。右のMMSE等化では、同じ16個の塊が(特にSNRが低いほど顕著に)より小さくまとまり、隣との混同が起きにくくなっています。MMSEがわずかな残留ISIと引き換えに雑音を抑え、結果として点がきれいに分離する様子が読み取れます。

MSE-SNR曲線 — 優劣はSNRで決まる

最後に、SNRを変えながらZFとMMSEのMSEがどう変化するかを描きます。理論で予想したとおり「高SNRでは両者一致、低SNRではMMSEが圧勝」になるはずです。SNRごとに多数のシンボルでモンテカルロ平均を取ります。

import numpy as np

h = np.array([1.0, 0.0, 0.7, 0.0, 0.4], dtype=complex)
N = 500
n_trials = 40
snr_db_list = np.arange(0, 31, 2)
rng = np.random.default_rng(1)

H = build_channel_matrix(h, N)
W_zf = zf_equalizer(H)   # ZFはSNRに依存しないので一度だけ作る

mse_zf, mse_mmse = [], []
for snr_db in snr_db_list:
    snr_lin = 10 ** (snr_db / 10)
    sigma_n2 = 1.0 / snr_lin
    W_mmse = mmse_equalizer(H, snr_lin)
    e_zf, e_mmse = 0.0, 0.0
    for _ in range(n_trials):
        s = qam16_symbols(N, rng)
        r = H @ s + np.sqrt(sigma_n2 / 2) * (
            rng.standard_normal(H.shape[0]) + 1j * rng.standard_normal(H.shape[0]))
        e_zf += np.mean(np.abs(s - W_zf @ r) ** 2)
        e_mmse += np.mean(np.abs(s - W_mmse @ r) ** 2)
    mse_zf.append(e_zf / n_trials)
    mse_mmse.append(e_mmse / n_trials)

mse_zf = np.array(mse_zf)
mse_mmse = np.array(mse_mmse)
print(f"SNR=0dB :  ZF={mse_zf[0]:.4f}, MMSE={mse_mmse[0]:.4f}")
print(f"SNR=30dB:  ZF={mse_zf[-1]:.5f}, MMSE={mse_mmse[-1]:.5f}")

このループで、各SNRにおけるZFとMMSEのMSEを多数回の試行で平均しています。出力からは、SNR=0 dB(雑音が信号と同程度)ではZFのMSEがMMSEの数倍に達する一方、SNR=30 dB(雑音がごく小さい)では両者がほぼ同じ値に収束することが読み取れます。これは式 (5) で $\frac{\sigma_n^2}{\sigma_s^2} \to 0$ のときMMSEがZFに一致する、という理論的帰結そのものです。続いてこれをグラフにします。

import numpy as np
import matplotlib.pyplot as plt

plt.figure(figsize=(8, 5))
plt.semilogy(snr_db_list, mse_zf, "o-", color="#d62728", lw=2, label="ZF等化器")
plt.semilogy(snr_db_list, mse_mmse, "s-", color="#1f77b4", lw=2, label="MMSE等化器")
plt.xlabel("SNR [dB]")
plt.ylabel("平均二乗誤差 MSE(対数軸)")
plt.title("ZF vs MMSE:MSEのSNR依存性")
plt.grid(alpha=0.3, which="both")
plt.legend()
plt.tight_layout()
plt.savefig("/tmp/mmse_mse_snr.png", dpi=130)
plt.show()

ZFとMMSEのMSE-SNR曲線

この曲線がZFとMMSEの関係を最も雄弁に物語っています。低SNR側(左)では、ZF(赤)のMSEがMMSE(青)よりはっきり大きく、両者の差が開いています。これは雑音が大きい領域でZFが雑音を増幅し、MMSEの正則化が効いている証拠です。一方、SNRが高くなる(右へ進む)につれて2本の曲線は次第に接近し、最終的にほぼ重なります。雑音がなければISIを完璧に消すZFが最適になり、MMSEもそこへ収束するのです。「どちらが良いか」はSNR次第——これがZFとMMSEの関係の結論です。

雑音増幅をスペクトルで確かめる

最後にもう一枚、ZFが「チャネルのノッチで雑音を増幅する」ことを周波数領域で直接確認します。ZF等化器の周波数応答 $1/H(e^{j\omega})$ とMMSEの応答を重ねて描きます。

import numpy as np
import matplotlib.pyplot as plt

h = np.array([1.0, 0.0, 0.7, 0.0, 0.4], dtype=complex)
w = np.linspace(0.01, np.pi, 512)
H_w = np.array([np.sum(h * np.exp(-1j * wk * np.arange(len(h)))) for wk in w])

# 周波数領域でのZF/MMSE等化器応答(スカラー版)
snr_lin = 10 ** (10 / 10)            # SNR=10 dB
W_zf_w = 1.0 / H_w
W_mmse_w = np.conj(H_w) / (np.abs(H_w) ** 2 + 1.0 / snr_lin)

plt.figure(figsize=(8, 5))
plt.plot(w / np.pi, 20 * np.log10(np.abs(W_zf_w)), color="#d62728", lw=2,
         label="ZF $1/H$(ノッチで急増)")
plt.plot(w / np.pi, 20 * np.log10(np.abs(W_mmse_w)), color="#1f77b4", lw=2,
         label="MMSE(増幅を頭打ちに)")
plt.xlabel("正規化周波数 $\\omega/\\pi$")
plt.ylabel("等化器の利得 [dB]")
plt.title("ZFとMMSEの等化器応答(SNR=10 dB)")
plt.grid(alpha=0.3); plt.legend()
plt.tight_layout()
plt.savefig("/tmp/mmse_eq_response.png", dpi=130)
plt.show()

ZFとMMSE等化器の周波数応答(SNR=10 dB)

この図は、ZFとMMSEの設計思想の違いをそのまま映しています。チャネルのノッチ(振幅が落ちる周波数)で、ZFの利得(赤)は天井知らずに跳ね上がります。これがその周波数の雑音を増幅する正体です。対してMMSE(青)は、同じノッチでも利得の上昇が途中で頭打ちになります。分母に $1/\mathrm{SNR}$ が加わっているため、$|H(e^{j\omega})|$ が小さくなっても等化器の利得が発散しないのです。MMSEは「持ち上げても無駄な(雑音しかない)周波数では、あえて持ち上げない」という賢い判断をしている、と読めます。

これらの図と数値実験で、理論の主張——ZFはISIを完璧に消すが雑音を増幅し、MMSEはSNRに応じて両者のバランスを取る——がすべて裏付けられました。

アイダイアグラムで「アイの開き」を見る

星座図は1点ごとの散らばりを見る図でしたが、波形そのものの健全さを見るにはアイダイアグラムが便利です。受信波形を1シンボル区間ごとに切り出して重ね描きすると、ISIが小さいほど中央に大きな「目(アイ)」が開きます。アイが閉じていると判定の余裕がなく、わずかな雑音で誤りが出ます。ここでは BPSK 波形を1シンボル遅延の反射を持つチャネルに通し、等化前と簡易ZF等化後のアイを比べます。

等化前後のアイダイアグラム

左の等化前では、反射波による ISI で波形が重なり合い、シンボル中心(点線)付近でアイがほとんど閉じています。これでは $+1/-1$ を判定する余裕がありません。右の等化後では、同じシンボル中心で上下にはっきりとアイが開き、判定マージンが回復しているのが読み取れます。ISIを取り除くことが「アイを開ける」ことに直結する、という時間領域の描像です。

シンボル誤り率(SER)で実用性能を測る

MSE は等化器の良し悪しを測る代理指標ですが、通信で最終的に効くのはシンボルを取り違える確率です。等化後の各シンボルを最近傍の格子点に判定し、送信シンボルと一致しない割合(SER)を SNR ごとにモンテカルロで実測しました。

シンボル誤り率のSNR依存性(16-QAM・実測)

等化なし(灰)は SNR をいくら上げても SER が約 0.8 で頭打ちになります。ISI が残っている限り、雑音を減らしても誤りは消えないのです。一方 ZF(赤)と MMSE(青)はほぼ重なりながら、SNR とともに SER が桁で落ちていきます。この設定では ZF と MMSE の SER 差はごく小さい(12 dB で 0.2226 対 0.2207)ものの、MMSE がわずかに下回り、かつ「等化すれば誤り率が下げられる」という等化器の存在意義が明確に見えます。

適応等化:チャネルが未知のときどう学習するか

これまではチャネル行列 $\bm{H}$ が既知である前提でした。しかし実際の受信機はチャネルを知りません。そこで既知の訓練系列との誤差を使ってタップ係数を逐次更新する適応等化が使われます。代表的な LMS(最急降下の確率版)と RLS(再帰最小二乗)の学習曲線を実測しました。

適応等化(LMS/RLS)の学習曲線

縦軸は二乗誤差の移動平均(対数)です。RLS(紫)は数十反復という圧倒的な速さで収束します。LMS はステップ幅 $\mu$ が小さい(緑)と収束は遅いが滑らかで最終誤差が低く、$\mu$ が大きい(橙)と立ち上がりは速いものの最終的な誤差床が高くなります。「速さと安定(定常誤差)はトレードオフ」という適応フィルタの基本則が読み取れます。RLS は計算量が大きい代わりにこのトレードオフを大きく改善します。

残留ISIを総合インパルス応答で見る

最後に、ZF と MMSE の「ISIの扱い方の違い」を最も直接的に示します。等化後の総合系は $\bm{W}\bm{H}$ で表され、これが単位行列に近いほど ISI が消えています。その中央列(あるシンボルが等化後にどの遅延へ漏れるか)を SNR=5 dB で棒グラフにしました。

等化後の総合インパルス応答WHと残留ISI

左の ZF では中央が厳密に $1$、それ以外は機械精度($10^{-16}$ 級)でゼロ——ISI は完全に除去されています。右の MMSE では中央が $0.715$ とやや下がり、前後の遅延に小さな漏れ(残留ISI)が残ります。MMSE はこのわずかな ISI をあえて許す代わりに雑音増幅を抑えており、低SNRではその取引が総合誤差で有利になります。ZF が「ISIゼロ最優先」、MMSE が「ISIと雑音の総和最小」という設計思想の違いが、この一枚に凝縮されています。

まとめ

本記事では、周波数選択性チャネルのISIを補償する2つの線形等化器、ゼロフォーシング等化器とMMSE等化器について、モデルの定式化から導出、実装までを解説しました。

  • 受信信号モデル: 多重路チャネルは畳み込み行列 $\bm{H}$ として作用し、受信信号は $\bm{r} = \bm{H}\bm{s} + \bm{n}$ と書ける。$\bm{H}$ の非対角成分がISIの正体である
  • ゼロフォーシング等化器: ISIをゼロに強制する最小二乗解 $\bm{W}_{\mathrm{ZF}} = (\bm{H}^H\bm{H})^{-1}\bm{H}^H$。$\bm{H}$ の擬似逆行列であり、ISIは完璧に消えるが、チャネルのノッチで雑音を激しく増幅する
  • 直交原理: 平均二乗誤差を最小にする線形推定では、誤差が観測と直交する($\mathbb{E}[(\bm{s} – \bm{W}\bm{r})\bm{r}^H] = \bm{0}$)。これは「観測の張る空間への射影が最良の推定」という幾何学的事実である
  • MMSE等化器: 直交原理から $\bm{W}_{\mathrm{MMSE}} = (\bm{H}^H\bm{H} + \frac{\sigma_n^2}{\sigma_s^2}\bm{I})^{-1}\bm{H}^H$ が導かれる。ZFとの違いは正則化項 $\frac{\sigma_n^2}{\sigma_s^2}\bm{I}$ だけで、これがリッジ回帰と同じく逆行列の暴れを抑える
  • トレードオフ: 高SNRではMMSEはZFに一致し、低SNRではMMSEが雑音増幅を抑えて圧勝する。優劣はSNRで決まる

等化器の核心は「逆問題をどう正則化するか」にあります。この視点は、ノイズのある観測から元の信号を復元するあらゆる問題——画像復元、レーダーのパルス圧縮、時系列予測——に共通しています。MMSEで使った直交原理は、まさにウィーナーフィルタやカルマンフィルタの背骨でもあります。

次のステップとして、以下の記事も参考にしてください。