衛星から地上へ送られてくる1枚の画像、深宇宙探査機が何億キロも先から細々と送ってくるテレメトリ、あるいはスマートフォンが基地局とやり取りするパケット——これらはすべてノイズだらけの通信路を通ってきます。受信機の手元に届くのは、送信されたビット列そのものではなく、雑音でぼやけたアナログ的な観測値です。ここで素朴な疑問が湧きます。「あるビットが0だったのか1だったのか、どうすれば最も確からしく当てられるのか？」
この問いに対して、2つの異なる「最適性」の答えがあります。1つは「受信系列全体として最も尤もらしい符号語はどれか」を探す系列MAP（最尤系列推定）で、これがViterbiアルゴリズムです。もう1つが「個々のビットについて、それが1である事後確率を最大化する」というビットMAPで、これを実現するのが本記事で扱うBCJRアルゴリズムです。BCJRは提案者のBahl, Cocke, Jelinek, Raviv の頭文字をとった名前で、1974年に発表されました。
長らくViterbiの陰に隠れていたBCJRが脚光を浴びたのは、1993年のターボ符号の登場がきっかけです。ターボ符号は2つの軟入力軟出力（SISO: Soft-In Soft-Out）復号器が「外部情報」を交換し合いながら反復的に推定を改善することでシャノン限界に肉薄しますが、この各復号器の中核を担うのがまさにBCJRなのです。Viterbiは硬い（0/1の）判定しか出せませんが、BCJRは「このビットが1である確からしさ」という軟らかい確率情報を出力できる——この違いが反復復号を成立させる決定的な鍵になります。
BCJRを理解すると、次のような分野が見えてきます。
- ターボ符号の反復復号: 3GPP LTE（4G）のデータチャネルや深宇宙通信規格CCSDSのターボ符号で、SISO復号の心臓部として使われています
- ターボ等化・ターボ受信機: 符号化と通信路（マルチパスやISI）を1つのトレリスとみなし、等化器と復号器の間で軟情報を交換する受信方式の基礎になります
- 連接符号・BICM-ID: 検出器と復号器の間で外部情報を反復交換するあらゆるシステムが、BCJRが出力する事後確率を土台にしています
本記事の内容
- ビットMAP復号の考え方と、系列MAP（Viterbi）との本質的な違い
- トレリス上の前向き再帰 $\alpha$、後ろ向き再帰 $\beta$、分岐メトリクス $\gamma$ の定義と導出（省略なし）
- ビット事後確率と対数尤度比（LLR）$L(u_k)$ の導出、外部情報（extrinsic）の分離
- 数値安定化のための log-MAP とヤコビ対数、max-log-MAP 近似
- ターボ符号での軟入力軟出力（SISO）反復復号におけるBCJRの役割
- 畳み込み符号のトレリス上で $\alpha/\beta/\gamma$ を計算しLLRを出力するPython実装、近似比較、BER評価
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 畳み込み符号とビタビアルゴリズムをわかりやすく解説 — トレリス図、状態遷移、Viterbi復号の基礎。本記事はこの上に立ちます
- ターボ符号の理論 — 反復復号と性能限界 — BCJRが部品として使われる全体像
- LDPC符号の理論 — スパースグラフ符号の原理と復号アルゴリズム — 同じく軟判定反復復号を行うもう一つのシャノン限界級符号
ビットMAP復号とは — 「系列」ではなく「ビット」を当てる
まず、どんな問題を解こうとしているのかを直感的につかみましょう。テスト問題に喩えてみます。100問の選択式テストがあり、採点には2つの流儀があるとします。1つは「全問の答案セットとして最も整合的なものを1つ選ぶ」やり方、もう1つは「各問について、正解である確率が最大の選択肢を独立に選ぶ」やり方です。前者は答案全体の一貫性を重視し、後者は1問1問の正答率を最大化します。
通信の文脈では、前者が系列MAP（MAP sequence estimation）で、送信された符号語系列 $\bm{u}$ 全体の事後確率 $P(\bm{u} \mid \bm{y})$ を最大にする系列を選びます。これは受信系列 $\bm{y}$ が与えられたときに最も尤もらしい「経路」をトレリス上で探す問題に帰着し、Viterbiアルゴリズムが効率的に解きます。系列誤り率（フレーム誤り率）を最小化するのが目的です。
後者がビットMAP（symbol-by-symbol MAP）で、各ビット $u_k$ について個別に事後確率 $P(u_k \mid \bm{y})$ を最大化します。
$$ \hat{u}_k = \arg\max_{u_k \in \{0,1\}} P(u_k \mid \bm{y}) $$
これはビット誤り率（BER）を最小化する判定です。BCJRアルゴリズムはこのビット事後確率 $P(u_k \mid \bm{y})$ を、全系列の情報を漏れなく使って効率的に計算します。
ここで重要なのは、ビットMAPの推定値を全ビット並べた系列 $(\hat{u}_1, \hat{u}_2, \dots)$ は、必ずしも正規の符号語になるとは限らないという点です。各ビットを独立に最尤判定するので、出来上がった系列がトレリス上の有効な経路に対応する保証はありません。一方Viterbiの出力は必ず有効な符号語です。「BERを最小化したいのか、フレーム誤り率を最小化したいのか」で使い分けるべき、というのが両者の関係です。
そして、ターボ符号やLDPC符号のように軟情報を反復交換するシステムでは、硬い判定ではなく「そのビットが1らしい度合い」という連続値の確信度が必要になります。BCJRはこの確信度を対数尤度比（LLR）という形で自然に出力できるため、SISO復号器として理想的なのです。
ではこの事後確率を、組合せ爆発を起こさずにどう計算するのか。鍵はトレリス構造にあります。次にその舞台となるトレリスとモデルを整理しましょう。
舞台設定 — トレリス、状態、観測モデル
BCJRが対象とするのは、送信側がマルコフ的な状態機械、すなわち畳み込み符号器のような有限状態機械であるシステムです。畳み込み符号とビタビアルゴリズムで見たように、符号器は時刻 $k$ で状態 $s_{k-1}$ から入力ビット $u_k$ に応じて次の状態 $s_k$ へ遷移し、同時に符号ビット $\bm{x}_k$ を出力します。この遷移を時間方向に並べたものがトレリスです。
記号を定義します。
- $u_k \in \{0, 1\}$: 時刻 $k$ の情報ビット（推定したい対象）
- $s_k \in \{0, 1, \dots, M-1\}$: 時刻 $k$ における符号器の状態（$M$ 個の状態）
- $\bm{x}_k$: 時刻 $k$ で符号器が出力する符号ビット（レート $1/n$ なら $n$ ビット）
- $\bm{y}_k$: 時刻 $k$ に対応する受信値（雑音が乗ったもの）
- $\bm{y} = (\bm{y}_1, \bm{y}_2, \dots, \bm{y}_N)$: 全受信系列（長さ $N$）
通信路はメモリレスのAWGN通信路を仮定します。BPSK変調 $x \in \{0,1\} \to (1-2x) \in \{+1, -1\}$ を施すと、受信値は
$$ y_{k,\ell} = (1 – 2 x_{k,\ell}) + n_{k,\ell}, \quad n_{k,\ell} \sim \mathcal{N}(0, \sigma^2) $$
となります（$\ell$ は時刻 $k$ 内の符号ビットのインデックス）。メモリレスなので、状態系列が与えられれば各時刻の観測は独立です。
ここでマルコフ性が効いてきます。状態 $s_k$ が与えられれば、それより前の系列と後の系列は条件付き独立になります。すなわち「現在の状態さえ分かれば、過去は未来に影響しない」という性質です。この性質こそが、後で全系列の確率を「前半パート」と「後半パート」に分解する根拠になります。
トレリス上の各「枝（branch）」は、状態対 $(s_{k-1}, s_k) = (s’, s)$ で指定される1つの遷移です。各枝には、その遷移を引き起こす入力ビット $u_k(s’, s)$ と、出力される符号ビット $\bm{x}_k(s’, s)$ が決定論的に対応づきます。BCJRはこの「枝」を確率の基本単位として扱います。
舞台が整いました。次は、求めたいビット事後確率を、この枝の確率を使ってどう書き下すかを見ていきます。
ビット事後確率の枝確率への分解
最終目標は $P(u_k = 0 \mid \bm{y})$ と $P(u_k = 1 \mid \bm{y})$ を求めることです。直感としては、「入力ビット $u_k$ が1になるような枝（状態遷移）すべての確率を足し合わせればよい」というものです。トレリス上で $u_k = 1$ に対応する枝の集合を $\mathcal{B}_1$、$u_k = 0$ に対応する枝の集合を $\mathcal{B}_0$ と書きます。
ベイズの定理から、事後確率は同時確率に比例します（$P(\bm{y})$ は $u_k$ によらない共通因子なので、最後にLLRを取るとき消えます）。
$$ P(u_k = b \mid \bm{y}) = \frac{P(u_k = b, \bm{y})}{P(\bm{y})} $$
そこで同時確率 $P(u_k = b, \bm{y})$ を求めます。ビット $u_k$ は枝 $(s’, s)$ によって決まるので、$u_k = b$ となる枝すべてについて同時確率 $P(s_{k-1} = s’, s_k = s, \bm{y})$ を足し合わせれば求まります。
$$ P(u_k = b, \bm{y}) = \sum_{(s’, s) \in \mathcal{B}_b} P(s_{k-1} = s’, s_k = s, \bm{y}) $$
ここで核心となる量を導入します。1つの枝の同時確率を
$$ \sigma_k(s’, s) \equiv P(s_{k-1} = s’, s_k = s, \bm{y}) $$
と定義します。これは「時刻 $k-1$ で状態 $s’$ にいて、時刻 $k$ で状態 $s$ に遷移し、かつ全受信系列 $\bm{y}$ が観測される」同時確率です。これさえ全枝について計算できれば、ビット事後確率は単なる足し算で得られます。
問題は $\sigma_k(s’, s)$ をどう計算するかです。ここでマルコフ性を使います。全受信系列を3つの時間区間に分けます。
$$
\bm{y} = (\underbrace{\bm{y}_1, \dots, \bm{y}_{k-1}}_{\bm{y}_{
同時確率を、過去・現在・未来に分解します。確率の乗法定理を繰り返し使うと、
$$
\sigma_k(s’, s) = P(s’, s, \bm{y}_{
ここでマルコフ性が威力を発揮します。状態 $s_k = s$ が分かっていれば、未来の観測 $\bm{y}_{>k}$ は過去 $(\bm{y}_{ $$
\sigma_k(s’, s) = \underbrace{P(s’, \bm{y}_{ と3つの因子の積に分解できます。この分解こそがBCJRの心臓部です。それぞれに名前を付けましょう。 $$
\alpha_{k-1}(s’) \equiv P(s_{k-1} = s’, \bm{y}_{ まとめると、 $$
\boxed{\;\sigma_k(s’, s) = \alpha_{k-1}(s’) \cdot \gamma_k(s’, s) \cdot \beta_k(s)\;}
$$ という美しい積の形になります。「過去 × 現在 × 未来」と覚えると忘れません。残る仕事は、$\alpha$ と $\beta$ をトレリスを1回ずつ掃くだけで計算する再帰式を導くことです。次節でこれを丁寧に導出します。 $\alpha_k(s)$ は「時刻 $k$ までの観測を使って、状態 $s$ にいる確からしさ」です。これを1つ前の時刻 $k-1$ の $\alpha_{k-1}$ から計算する再帰式を作ります。直感的には、「現在の状態 $s$ には、複数の前の状態 $s’$ から枝が入ってくる。それぞれの『これまでの確からしさ $\alpha_{k-1}(s’)$』に『その枝が起きる確率 $\gamma_k(s’, s)$』を掛けて足し合わせれば、状態 $s$ の新しい確からしさになる」ということです。 定義から出発します。 $$
\alpha_k(s) = P(s_k = s, \bm{y}_{\le k}) = P(s_k = s, \bm{y}_{ ここで、時刻 $k-1$ の状態 $s’$ について周辺化（あらゆる前状態の可能性を足し上げる）します。これは全確率の法則です。 $$
\alpha_k(s) = \sum_{s’} P(s_{k-1} = s’, s_k = s, \bm{y}_{ 次に、同時確率を乗法定理で分けます。$\bm{y}_{ $$
\alpha_k(s) = \sum_{s’} P(s_{k-1} = s’, \bm{y}_{ ここで再びマルコフ性を使います。状態 $s’$ が与えられれば、時刻 $k$ の遷移と観測 $\bm{y}_k$ は過去の観測 $\bm{y}_{ $$
P(s_k = s, \bm{y}_k \mid s_{k-1} = s’, \bm{y}_{ そして第1因子はまさに $\alpha_{k-1}(s’)$ の定義そのものです。したがって、 $$
\boxed{\;\alpha_k(s) = \sum_{s’} \alpha_{k-1}(s’) \, \gamma_k(s’, s)\;}
$$ という前向き再帰が得られました。トレリスを左から右へ1回掃きながら、各状態の $\alpha$ を更新していくだけです。 初期条件は、符号器が既知の状態（通常は全ゼロ状態 $s_0 = 0$）からスタートすることを反映して、 $$
\alpha_0(s) = \begin{cases} 1 & s = 0 \\ 0 & s \neq 0 \end{cases}
$$ とします。「最初は確実に状態0にいる」という事前知識です。 実装上の重要な注意点として、$\alpha_k(s)$ は確率の積を繰り返すため、時刻が進むにつれて指数的に小さくなり、すぐにアンダーフローします。これを防ぐため、各時刻で $\sum_s \alpha_k(s) = 1$ となるように正規化するのが定石です。正規化してもLLRを取るときに比で消えるので、結果は変わりません。 前向きに「これまで」を積み上げる仕組みが分かりました。次は時間を逆向きに掃いて「これから」を表す $\beta$ を導きます。構造は $\alpha$ と鏡写しです。 $\beta_k(s)$ は「時刻 $k$ で状態 $s$ にいるとき、それ以降に受信系列 $\bm{y}_{>k}$ が生じる尤度」です。今度はトレリスを右から左へ掃きます。直感は $\alpha$ の逆で、「状態 $s$ からは複数の次状態 $s’$ へ枝が出ていく。それぞれの枝の確率 $\gamma_{k+1}(s, s’)$ に、その先の未来の尤度 $\beta_{k+1}(s’)$ を掛けて足し合わせれば、状態 $s$ の未来尤度になる」というものです。 定義から始めます。 $$
\beta_k(s) = P(\bm{y}_{>k} \mid s_k = s) = P(\bm{y}_{k+1}, \bm{y}_{>k+1} \mid s_k = s)
$$ 次の時刻 $k+1$ の状態 $s’$ について周辺化します。 $$
\beta_k(s) = \sum_{s’} P(s_{k+1} = s’, \bm{y}_{k+1}, \bm{y}_{>k+1} \mid s_k = s)
$$ 同時確率を分解します。$\bm{y}_{>k+1}$ を後ろに取り出すと、 $$
\beta_k(s) = \sum_{s’} P(s_{k+1} = s’, \bm{y}_{k+1} \mid s_k = s) \cdot P(\bm{y}_{>k+1} \mid s_{k+1} = s’, s_k = s, \bm{y}_{k+1})
$$ マルコフ性により、次状態 $s’$ が与えられれば、その先の観測 $\bm{y}_{>k+1}$ は現在の状態 $s$ や観測 $\bm{y}_{k+1}$ に依存しません。条件部を簡約すると、 $$
P(\bm{y}_{>k+1} \mid s_{k+1} = s’, s_k = s, \bm{y}_{k+1}) = P(\bm{y}_{>k+1} \mid s_{k+1} = s’) = \beta_{k+1}(s’)
$$ そして第1因子は枝メトリクス $\gamma_{k+1}(s, s’)$ そのものです。よって、 $$
\boxed{\;\beta_k(s) = \sum_{s’} \beta_{k+1}(s’) \, \gamma_{k+1}(s, s’)\;}
$$ という後ろ向き再帰が得られました。$\alpha$ と完全に対称な形をしています。 初期条件（実際には終端条件）は、符号器を全ゼロ状態でターミネーションする（末尾にテールビットを入れて状態0に戻す）場合、 $$
\beta_N(s) = \begin{cases} 1 & s = 0 \\ 0 & s \neq 0 \end{cases}
$$ とします。ターミネーションしない場合は、すべての状態が等確率という事前知識を表す $\beta_N(s) = 1/M$ とします。$\beta$ も $\alpha$ 同様に各時刻で正規化してアンダーフローを防ぎます。 $\alpha$ は「過去から流れてくる確信」、$\beta$ は「未来から逆流してくる確信」で、トレリスの各点で両者が出会います。残るは両者をつなぐ枝メトリクス $\gamma$ の中身です。これがAWGN通信路や事前情報とどう結びつくかを次に見ます。 $\gamma_k(s’, s)$ は「状態 $s’$ から $s$ への1枝が、観測 $\bm{y}_k$ とともに起きる確率」です。この枝には決まった入力ビット $u_k$ と出力符号ビット $\bm{x}_k$ が対応するので、$\gamma$ は「入力ビットの事前確率」と「観測が出力ビットにどれだけ合致するかの尤度」の積に分解できます。これを丁寧に展開します。 定義を乗法定理で分けます。 $$
\gamma_k(s’, s) = P(s_k = s, \bm{y}_k \mid s_{k-1} = s’) = P(s_k = s \mid s_{k-1} = s’) \cdot P(\bm{y}_k \mid s_{k-1} = s’, s_k = s)
$$ 第1因子 $P(s_k = s \mid s_{k-1} = s’)$ は遷移確率です。状態対 $(s’, s)$ が有効な枝なら、それは入力ビット $u_k = u(s’, s)$ で一意に決まる遷移なので、この確率は入力ビットの事前確率 $P(u_k)$ に等しくなります（無効な遷移なら0）。 第2因子 $P(\bm{y}_k \mid s’, s)$ は、その枝が出力する符号ビット $\bm{x}_k(s’, s)$ が与えられたときの観測尤度です。AWGN通信路では各符号ビットが独立に雑音を受けるので、 $$
P(\bm{y}_k \mid \bm{x}_k) = \prod_{\ell=1}^{n} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\!\left(-\frac{(y_{k,\ell} – (1 – 2x_{k,\ell}))^2}{2\sigma^2}\right)
$$ となります。まとめると枝メトリクスは $$
\gamma_k(s’, s) = P(u_k) \cdot P(\bm{y}_k \mid \bm{x}_k(s’, s))
$$ です。ここで事前確率 $P(u_k)$ が独立した因子として現れる点が極めて重要です。これがターボ復号で「もう一方の復号器からの外部情報」を取り込む入口になります。最初の復号では情報ビットが等確率 $P(u_k=0) = P(u_k=1) = 1/2$ と仮定しますが、反復が進むと相手から渡された事前LLRに基づいてここが更新されます。 事前確率をLLR表現と結びつけておきましょう。事前LLRを $L_a(u_k) = \ln \frac{P(u_k=0)}{P(u_k=1)}$ とすると、 $$
P(u_k = b) = \frac{\exp(-b \, L_a(u_k))}{1 + \exp(-L_a(u_k))} \propto \exp\!\left(-b \, L_a(u_k)\right)
$$ と書け（$b \in \{0,1\}$）、比例係数は $u_k$ に依存しないので $\gamma$ の比を取るときに消えます。これで $\alpha, \beta, \gamma$ がすべて揃いました。次はこれらを組み合わせて、いよいよ目的のLLRを組み立てます。 すべての部品が揃ったので、ビット事後確率を計算します。前々節で導いた分解 $\sigma_k(s’, s) = \alpha_{k-1}(s’) \gamma_k(s’, s) \beta_k(s)$ を使うと、ビット $u_k = b$ の同時確率は $$
P(u_k = b, \bm{y}) = \sum_{(s’, s) \in \mathcal{B}_b} \alpha_{k-1}(s’) \, \gamma_k(s’, s) \, \beta_k(s)
$$ です。BCJRが出力するのは、これらの比の対数、すなわち事後対数尤度比です。 $$
L(u_k) = \ln \frac{P(u_k = 0 \mid \bm{y})}{P(u_k = 1 \mid \bm{y})} = \ln \frac{P(u_k = 0, \bm{y})}{P(u_k = 1, \bm{y})} = \ln \frac{\displaystyle\sum_{(s’, s) \in \mathcal{B}_0} \alpha_{k-1}(s’) \gamma_k(s’, s) \beta_k(s)}{\displaystyle\sum_{(s’, s) \in \mathcal{B}_1} \alpha_{k-1}(s’) \gamma_k(s’, s) \beta_k(s)}
$$ 共通因子 $P(\bm{y})$ が分母分子で約分されて消えていることに注意してください。$L(u_k) > 0$ なら $u_k = 0$ が確からしく、$< 0$ なら $u_k = 1$、絶対値が確信度です。硬判定は $\hat{u}_k = \text{sign}$ で取れますが、SISO復号ではこの軟らかい値そのものを次段に渡します。 ここからが反復復号の本質です。$\gamma_k$ には先ほど見たように事前確率の因子 $P(u_k)$（事前LLR $L_a(u_k)$）と通信路尤度が含まれていました。情報ビットそのものが受信値 $y_k^{(s)}$（系統ビット）として直接送られる系統符号の場合、$\gamma$ の中身はさらに分離できます。レート $1/2$ の系統的再帰畳み込み符号（ターボ符号の構成要素）を例にとると、時刻 $k$ の受信値は系統ビット $y_k^s$ とパリティビット $y_k^p$ から成り、 $$
\gamma_k(s’, s) \propto \exp\!\left(-b \, L_a(u_k)\right) \cdot \exp\!\left(\frac{1}{\sigma^2} y_k^s (1 – 2b)\right) \cdot \exp\!\left(\frac{1}{\sigma^2} y_k^p (1 – 2 x_k^p)\right)
$$ と3つの因子に分かれます（$b = u_k$）。AWGN尤度 $\exp\!\big(-(y-(1-2x))^2/(2\sigma^2)\big)$ を展開すると、$y^2$ と $(1-2x)^2=1$ の項は $b$ に依存しない共通因子として比例係数に吸収され、残る交差項が $\exp\!\big(y(1-2x)/\sigma^2\big)$ になります（係数は $1/\sigma^2$ であって $2/\sigma^2$ ではない点に注意）。ここで通信路信頼度 $L_c = 2/\sigma^2$ を導入すると $1/\sigma^2 = L_c/2$ なので、 $$
\gamma_k(s’, s) \propto \exp\!\left(-b \, L_a(u_k)\right) \cdot \exp\!\left(\frac{L_c}{2} y_k^s (1 – 2b)\right) \cdot \exp\!\left(\frac{L_c}{2} y_k^p (1 – 2 x_k^p)\right)
$$ と書けます。第1・第2因子は枝の集合 $\mathcal{B}_0, \mathcal{B}_1$ で共通（$b$ だけで決まり状態対によらない）なので、LLRの和から括り出せます。系統ビットの項は $\exp\!\big((L_c/2) y_k^s (1-2b)\big)$ で、$b=0$ なら $+(L_c/2)y_k^s$、$b=1$ なら $-(L_c/2)y_k^s$ の指数となるため、$\mathcal{B}_0$ と $\mathcal{B}_1$ の差を取ると $L_c\, y_k^s$ がそのまま現れます。結果として、事後LLRは3つの項に分解されます。 $$
\boxed{\;L(u_k) = \underbrace{L_c \, y_k^s}_{\text{通信路値}} + \underbrace{L_a(u_k)}_{\text{事前情報}} + \underbrace{L_e(u_k)}_{\text{外部情報}}\;}
$$ ここで第3項の外部情報（extrinsic information） $L_e(u_k)$ が、ターボ復号の通貨です。これは $$
L_e(u_k) = \ln \frac{\displaystyle\sum_{(s’, s) \in \mathcal{B}_0} \alpha_{k-1}(s’) \, \gamma_k^{e}(s’, s) \, \beta_k(s)}{\displaystyle\sum_{(s’, s) \in \mathcal{B}_1} \alpha_{k-1}(s’) \, \gamma_k^{e}(s’, s) \, \beta_k(s)}
$$ で、$\gamma_k^e$ は系統ビットと事前情報の寄与を除いた、パリティビットの尤度だけを含む枝メトリクスです。 なぜわざわざ「外部情報」だけを取り出すのでしょうか。直感は「自分が今出した結論のうち、自分が最初に受け取った入力（系統値と事前情報）に由来する部分を、もう一方の復号器に返してはいけない」というものです。もし全事後LLR $L(u_k)$ をそのまま相手に渡すと、相手はそれを新たな事前情報として使い、自分が以前渡した情報が増幅されて返ってくる——情報の二重計上（正のフィードバックループ）が起き、復号が不安定になります。外部情報 $L_e$ は「このビットについて、自分のパリティ検査だけが新たに知り得た、相手がまだ知らない情報」だけを含むので、これを交換することで健全な反復が成立します。 外部情報の概念が反復復号の鍵だと分かりました。しかしここまでの計算は確率の積と和で、実装するとアンダーフローや乗算コストが問題になります。次は対数領域での安定な計算法を導きます。 これまでの式は確率（小さな正の数）の積を多数含むため、有限精度の計算機ではすぐにアンダーフローします。また乗算は加算より高コストです。そこで全量を対数領域に移します。対数化した量を $$
\tilde{\alpha}_k(s) = \ln \alpha_k(s), \quad \tilde{\beta}_k(s) = \ln \beta_k(s), \quad \tilde{\gamma}_k(s’, s) = \ln \gamma_k(s’, s)
$$ と定義します。すると枝メトリクスの積 $\alpha \gamma \beta$ は和 $\tilde{\alpha} + \tilde{\gamma} + \tilde{\beta}$ になります。問題は前向き・後ろ向き再帰に現れる「確率の和」です。$\alpha_k(s) = \sum_{s’} \alpha_{k-1}(s’) \gamma_k(s’, s)$ を対数で書くと、 $$
\tilde{\alpha}_k(s) = \ln \sum_{s’} \exp\!\left(\tilde{\alpha}_{k-1}(s’) + \tilde{\gamma}_k(s’, s)\right)
$$ という「指数の和の対数（log-sum-exp）」が現れます。これを安定に計算する道具がヤコビ対数（Jacobian logarithm）です。2変数版は次の恒等式です。 $$
\ln(e^a + e^b) = \max(a, b) + \ln\!\left(1 + e^{-|a – b|}\right)
$$ この式は厳密です。導出は簡単で、$a \ge b$ とすると $\ln(e^a + e^b) = \ln\!\big(e^a(1 + e^{b-a})\big) = a + \ln(1 + e^{b-a})$ となり、$a < b$ の場合とまとめると上式になります。右辺第1項 $\max(a,b)$ がアンダーフローを防ぎ、第2項 $\ln(1 + e^{-|a-b|})$ は $0 \le \ln 2$ の範囲に収まる小さな補正項です。この補正項を補正関数 $f_c(|a-b|) = \ln(1 + e^{-|a-b|})$ と呼びます。 この演算を $\max^*$（max-star）と書きます。 $$
\max^*(a, b) \equiv \ln(e^a + e^b) = \max(a, b) + \ln\!\left(1 + e^{-|a-b|}\right)
$$ 3項以上の log-sum-exp は $\max^*$ を再帰的に適用すれば計算できます。$\max^*(a, b, c) = \max^*(\max^*(a, b), c)$ です。補正関数 $f_c$ は小さなルックアップテーブル（典型的に8段程度）で十分な精度が得られるため、ハードウェア実装でも乗算なしで処理できます。 この $\max^*$ を厳密に使うのがlog-MAPです。前向き再帰は $$
\tilde{\alpha}_k(s) = \max^*_{s’} \left(\tilde{\alpha}_{k-1}(s’) + \tilde{\gamma}_k(s’, s)\right)
$$ となり、事後LLRは $$
L(u_k) = \max^*_{(s’,s) \in \mathcal{B}_0}\!\big(\tilde{\alpha}_{k-1}(s’) + \tilde{\gamma}_k(s’,s) + \tilde{\beta}_k(s)\big) – \max^*_{(s’,s) \in \mathcal{B}_1}\!\big(\cdots\big)
$$ と、2つの $\max^*$ の差で書けます。log-MAPは数値的に安定で、かつ元のMAP（BCJR）と数学的に完全に等価です。 さらに計算を簡略化したのがmax-log-MAPです。これは補正項 $\ln(1 + e^{-|a-b|})$ を無視し、$\max^*(a, b) \approx \max(a, b)$ と近似します。すると $$
\tilde{\alpha}_k(s) \approx \max_{s’} \left(\tilde{\alpha}_{k-1}(s’) + \tilde{\gamma}_k(s’, s)\right)
$$ となり、補正関数のテーブルすら不要になります。この近似のもとでBCJRの $\alpha$ 再帰はViterbiのパスメトリクス更新と同じ「加算と最大値選択（add-compare-select）」になります。実は、max-log-MAPの前向きパスはViterbiの前向きパスメトリクスそのものであり、後ろ向きにも同じことを行ってビットLLRを近似する、という見方ができます。 max-log-MAPの代償は性能のわずかな劣化です。補正項を捨てることでLLRの大きさが過大評価される傾向があり、ターボ復号では典型的に0.3〜0.5 dB程度の損失が生じます。実用上はこの損失をスケーリング係数（外部情報に $0.7$ 程度を乗じる）で補償するscaled max-log-MAPが広く使われます。計算量と性能のトレードオフで、log-MAPは性能最良・計算重め、max-log-MAPは計算軽量・わずかに性能劣化、という関係です。 理論が一通り揃いました。次はBCJRが実際にどう反復復号の部品として働くか、ターボ符号の文脈で整理します。 ターボ符号がなぜシャノン限界に迫れるのか、その仕組みをBCJRの視点から見ましょう。ターボ符号は、同じ情報ビット列を2つの再帰的組織的畳み込み符号器（RSC）で符号化します。ただし2つ目の符号器にはインタリーバで順序を入れ替えた情報ビットを入れます。受信側には系統ビット $\bm{y}^s$、1つ目のパリティ $\bm{y}^{p1}$、2つ目のパリティ $\bm{y}^{p2}$ が届きます。 復号は2つのBCJR（SISO復号器）が交互に働く反復です。1サイクルは次のように進みます。 このループを数回（典型的に4〜8回）繰り返すと、2つの復号器が互いに「自分のパリティから新たに分かったこと」だけを教え合い、推定が単調に研ぎ澄まされていきます。最終的に、いずれかの復号器の事後LLR $L(u_k) = L_c y_k^s + L_{a} + L_e$ の符号で硬判定して復号ビットを得ます。 ここで、なぜViterbiではなくBCJRでなければならないのかが明確になります。Viterbiは硬判定（0/1）しか出せないので、ステップ1で次段に渡せる「確信度つきの軟情報」を生成できません。BCJRは各ビットの事後確率を連続値LLRで出力でき、しかもその中から外部情報だけをきれいに分離できる——この2点が反復復号の生命線です。 なぜインタリーバが必要かも直感的に押さえておきましょう。2つの復号器が「相関のない別の視点」から同じビットを見るためです。インタリーバで順序を入れ替えると、復号器1で隣接していたビット同士が復号器2では遠く離れます。一方の復号器でたまたま信頼度が低くなったビットも、もう一方では別の文脈に置かれるので、互いの弱点を補い合えます。これは2人の専門家がバラバラの順序で同じ答案を採点し、意見交換するようなものです。 反復が進むにつれて外部情報の信頼度がどう上がるかは、EXITチャート（相互情報量の伝達特性）で解析されますが、これはLDPC符号の記事でも触れた手法と同じ枠組みです。実際にこの反復がBERをどう改善するか、Pythonで確かめましょう。 ここからは、簡単な畳み込み符号のトレリス上でBCJRを実装し、$\alpha/\beta/\gamma$ を計算してLLRを出力します。題材として、メモリ2・レート $1/2$ の系統的再帰畳み込み符号（生成多項式 $(1, 5/7)_8$、ターボ符号の標準構成要素）を使います。4状態のトレリスです。まず符号器のトレリス構造を定義します。 このトレリス表は、各状態と入力ビットの組に対して、次状態・系統ビット・パリティビットを与えます。出力を見ると、各状態から入力0と1で2本の枝が出ており、4状態すべてで遷移が定義されていることが確認できます。系統符号なので系統ビットは入力 $u$ そのものになっています。この決定論的な遷移表が、$\gamma$ を計算するときの「どの枝がどの符号ビットを出すか」の参照表になります。 次に、このトレリスとAWGN受信値から $\gamma$（対数枝メトリクス）を計算する関数を実装します。 このコードは、各枝に対応する系統・パリティビットをBPSK信号に変換し、AWGN尤度の指数部（対数領域での値）と事前LLR項を足して対数枝メトリクスを返します。事前LLR 続いて前向き $\tilde\alpha$ と後ろ向き $\tilde\beta$ の再帰を、log-MAP（ このパーツを使って、BCJR本体（前向き・後ろ向き・LLR）を組み立てます。 ここまでで前向き再帰を計算しました。各時刻で全枝を走査し、入ってくる枝の 後ろ向き再帰とLLR計算を続けます。 後ろ向き再帰は前向きと鏡写しで、終端状態から逆方向に各状態の LLR計算では、$u_k=0$ の枝集合 $\mathcal{B}_0$ と $u_k=1$ の枝集合 $\mathcal{B}_1$ をそれぞれ 実際に動かして、軟出力LLRが出ることを確かめます。 この実行結果から、BCJRが各ビットについて符号付きの実数LLRを出力していることが分かります。LLRの絶対値が大きいビットは確信度が高く、0に近いビットは判定が際どいことを表しています。硬判定（LLRの符号）の系列が真のビット列とほぼ一致していれば、たった1回のBCJRパスでも低SNRでなければ正しく復号できることが読み取れます。 次に、log-MAPとmax-log-MAPでLLRがどう異なるかを比較します。 このグラフから、2つの方式の違いが読み取れます。両者のLLRは符号（=硬判定）はほぼ常に一致しますが、max-log-MAPの方が絶対値がやや大きめに出る傾向があります。これは補正項 $\ln(1+e^{-|a-b|})$ を捨てたために確信度が過大評価されるためで、1回のパスでは硬判定にほとんど影響しませんが、反復復号で外部情報を交換すると誤差が蓄積し、約0.3〜0.5 dBの性能差となって現れます。 最後に、ターボ風の反復をしない単一RSCのBCJR（log-MAP／max-log-MAP）と、軟判定をしないハード判定のBERをSNRに対して比較します。 このBER曲線から、いくつかの点が読み取れます。第一に、SNRが上がるにつれてBERが単調に下がり、軟判定MAP復号が確かに誤りを訂正していることが分かります。第二に、log-MAPとmax-log-MAPのBERはほぼ重なり、単一符号・1パスではmax-log-MAPの近似損失はごくわずかです。この差はターボ符号として複数回反復したときに顕在化します。単一RSCではターボ符号本来の急峻なウォーターフォールは現れませんが、これに第2の符号器とインタリーバを加え、外部情報を反復交換することで、シャノン限界に迫る性能が得られるのです。 本記事では、ターボ符号のSISO復号を支えるBCJR（MAP）アルゴリズムを、確率分解の原理から実装まで解説しました。 BCJRは、確率の分解という美しい原理が、反復復号という現代通信の基盤技術へ直結する好例です。トレリス上の前向き・後ろ向き再帰という発想は、隠れマルコフモデルのForward-Backwardアルゴリズムやカルマンスムーザとも数学的に同根であり、信号処理・機械学習の広い領域に通じています。 次のステップとして、以下の記事も参考にしてください。
前向き再帰 α の導出
後ろ向き再帰 β の導出
分岐メトリクス γ の中身
事後LLRの組み立てと外部情報の分離
log-MAP と max-log-MAP — 対数領域での安定化
ターボ符号でのSISO反復復号におけるBCJRの役割
Pythonでの実装
import numpy as np
import matplotlib.pyplot as plt
# 系統的再帰畳み込み符号 (1, 5/7)_8, メモリ2, 4状態, レート1/2
# 状態 = (レジスタ s1, s2) を整数 0..3 で表現
# 帰還多項式 7=(111), 順方向 5=(101)
def build_trellis():
"""トレリスの遷移表を構築する。
返り値 trans[state][u] = (next_state, systematic_bit, parity_bit)"""
trans = {}
for state in range(4):
s1 = (state >> 1) & 1 # 上位ビット
s2 = state & 1 # 下位ビット
for u in range(2):
# 帰還: fb = u XOR s1 XOR s2 (帰還多項式7=1+D+D^2)
fb = u ^ s1 ^ s2
# パリティ出力: p = fb XOR s2 (順方向5=1+D^2 を帰還構造で)
p = fb ^ s2
# 状態更新（シフトレジスタに fb を入力）
ns1 = fb
ns2 = s1
next_state = (ns1 << 1) | ns2
trans[(state, u)] = (next_state, u, p) # 系統ビットは u そのもの
return trans
TRELLIS = build_trellis()
for (st, u), (ns, sb, pb) in sorted(TRELLIS.items()):
print(f"状態{st} 入力{u} -> 次状態{ns}, 系統{sb}, パリティ{pb}")
def compute_log_gamma(ys, yp, La, Lc):
"""各時刻・各枝の対数枝メトリクス log_gamma[(state,u)] を計算。
ys, yp: 系統・パリティ受信値（1時刻分のスカラー）
La: 事前LLR, Lc: 通信路信頼度 = 2/sigma^2"""
log_gamma = {}
for (state, u), (ns, sb, pb) in TRELLIS.items():
# BPSK: bit b -> 信号 (1-2b)
x_s = 1 - 2 * sb
x_p = 1 - 2 * pb
# 事前項: -u*La を対数で（定数項は比で消える）
prior = -u * La
# 通信路項: (Lc/2)*(ys*x_s + yp*x_p)
channel = 0.5 * Lc * (ys * x_s + yp * x_p)
log_gamma[(state, u)] = prior + channel
return log_gamma
La が0なら通信路項のみ、反復が進めば La が効いてくる構造になっており、ターボ復号の事前情報入力口がここにあることが読み取れます。max*）と max-log-MAP（max）の両方で実装します。まず max* 演算を定義します。def maxstar(a, b, approx=False):
"""ヤコビ対数 max*(a,b) = ln(e^a + e^b)。
approx=True なら補正項を捨てる max-log-MAP 近似。"""
m = np.maximum(a, b)
if approx:
return m
# 補正項 ln(1 + e^{-|a-b|}) を加える（log-MAP）
return m + np.log1p(np.exp(-np.abs(a - b)))
def maxstar_list(vals, approx=False):
"""複数要素の max* を逐次適用で計算"""
acc = vals[0]
for v in vals[1:]:
acc = maxstar(acc, v, approx)
return acc
maxstar は2引数版で、approx=False なら補正項 log1p(exp(-|a-b|)) を加えて厳密なヤコビ対数を、True なら単なる max を返します。np.log1p を使うことで小さな引数での数値精度を保っています。この1つの関数の引数を切り替えるだけで log-MAP と max-log-MAP を統一的に扱える点が実装上のポイントです。def bcjr_decode(ys, yp, La, Lc, approx=False, terminated=True):
"""1つのRSC符号に対するBCJR(SISO)復号。
ys, yp: 系統・パリティ受信値の配列（長さ N）
La: 事前LLR配列（長さ N）, Lc: 通信路信頼度
返り値: 事後LLR L, 外部情報 Le（長さ N）"""
N = len(ys)
NEG = -1e300 # 対数領域の「ゼロ確率」
# 各時刻の対数gammaを事前計算
log_gamma = [compute_log_gamma(ys[k], yp[k], La[k], Lc) for k in range(N)]
# 前向き alpha: alpha[k][state]
log_alpha = np.full((N + 1, 4), NEG)
log_alpha[0, 0] = 0.0 # 初期状態0
for k in range(1, N + 1):
for (state, u), (ns, sb, pb) in TRELLIS.items():
cand = log_alpha[k - 1, state] + log_gamma[k - 1][(state, u)]
log_alpha[k, ns] = maxstar(log_alpha[k, ns], cand, approx)
log_alpha[k] -= np.max(log_alpha[k]) # 正規化（最大値を引く）
return log_alpha, log_gamma # 次のセルで続ける
alpha + gamma を max* で集約しています。各時刻の最後に最大値を引いて正規化し、対数領域でのオーバーフロー／アンダーフローを防いでいます。初期状態を log_alpha[0,0]=0（確率1）、他を NEG（確率0）としているのが「最初は確実に状態0」という初期条件の実装です。def bcjr_full(ys, yp, La, Lc, approx=False, terminated=True):
"""BCJR完全版: 前向き・後ろ向き・LLR・外部情報を返す"""
N = len(ys)
NEG = -1e300
log_gamma = [compute_log_gamma(ys[k], yp[k], La[k], Lc) for k in range(N)]
# 前向き
log_alpha = np.full((N + 1, 4), NEG)
log_alpha[0, 0] = 0.0
for k in range(1, N + 1):
new = np.full(4, NEG)
for (state, u), (ns, sb, pb) in TRELLIS.items():
cand = log_alpha[k - 1, state] + log_gamma[k - 1][(state, u)]
new[ns] = maxstar(new[ns], cand, approx)
log_alpha[k] = new - np.max(new)
# 後ろ向き
log_beta = np.full((N + 1, 4), NEG)
if terminated:
log_beta[N, 0] = 0.0 # 終端状態0
else:
log_beta[N, :] = 0.0 # 全状態等確率
for k in range(N - 1, -1, -1):
new = np.full(4, NEG)
for (state, u), (ns, sb, pb) in TRELLIS.items():
cand = log_beta[k + 1, ns] + log_gamma[k][(state, u)]
new[state] = maxstar(new[state], cand, approx)
log_beta[k] = new - np.max(new)
return log_alpha, log_beta, log_gamma
beta を集約します。terminated=True なら終端状態0に確率を集中させ、そうでなければ全状態を等確率にする初期条件の違いも実装に反映されています。これで $\tilde\alpha, \tilde\beta, \tilde\gamma$ が揃ったので、LLRを組み立てられます。def compute_llr(log_alpha, log_beta, log_gamma, ys, La, Lc, approx=False):
"""事後LLR L(u_k) と外部情報 Le(u_k) を計算"""
N = len(ys)
L = np.zeros(N)
for k in range(N):
terms0, terms1 = [], []
for (state, u), (ns, sb, pb) in TRELLIS.items():
val = log_alpha[k, state] + log_gamma[k][(state, u)] + log_beta[k + 1, ns]
if u == 0:
terms0.append(val)
else:
terms1.append(val)
# L = ln( sum_{u=0} ... / sum_{u=1} ... )
L[k] = maxstar_list(terms0, approx) - maxstar_list(terms1, approx)
# 外部情報 = 事後 - 通信路系統値項 - 事前情報
Le = L - Lc * ys - La
return L, Le
max* で集約し、その差を取って事後LLRを得ています。最後の行で事後LLRから通信路系統値項 Lc*ys と事前情報 La を引き、外部情報 Le を分離しています。理論で導いた $L = L_c y^s + L_a + L_e$ の関係がそのままコードに現れています。def encode_rsc(u_bits):
"""情報ビット列をRSC符号化（系統 + パリティ）"""
state = 0
sys_out, par_out = [], []
for u in u_bits:
ns, sb, pb = TRELLIS[(state, u)]
sys_out.append(sb)
par_out.append(pb)
state = ns
return np.array(sys_out), np.array(par_out)
np.random.seed(0)
N = 12
u_true = np.random.randint(0, 2, N)
sys_b, par_b = encode_rsc(u_true)
# AWGN通信路
EbN0_dB = 2.0
R = 0.5
sigma = np.sqrt(1.0 / (2 * R * 10 ** (EbN0_dB / 10)))
Lc = 2.0 / sigma ** 2
ys = (1 - 2 * sys_b) + sigma * np.random.randn(N)
yp = (1 - 2 * par_b) + sigma * np.random.randn(N)
La = np.zeros(N) # 初回は事前情報なし
la_a, lb_a, lg_a = bcjr_full(ys, yp, La, Lc, approx=False)
L, Le = compute_llr(la_a, lb_a, lg_a, ys, La, Lc, approx=False)
u_hat = (L < 0).astype(int)
print("真のビット :", u_true)
print("復号ビット :", u_hat)
print("事後LLR :", np.round(L, 2))
print("誤りビット数:", np.sum(u_hat != u_true))
# 同じ受信値で2つの近似を比較
L_logmap, _ = compute_llr(*[*bcjr_full(ys, yp, La, Lc, approx=False)], ys, La, Lc, approx=False)
L_maxlog, _ = compute_llr(*[*bcjr_full(ys, yp, La, Lc, approx=True)], ys, La, Lc, approx=True)
plt.figure(figsize=(9, 5))
idx = np.arange(N)
plt.bar(idx - 0.2, L_logmap, width=0.4, label='log-MAP (exact)')
plt.bar(idx + 0.2, L_maxlog, width=0.4, label='max-log-MAP (approx)')
plt.axhline(0, color='k', lw=0.8)
plt.xlabel('bit index $k$')
plt.ylabel('posterior LLR $L(u_k)$')
plt.title('log-MAP vs max-log-MAP LLR')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
def simulate_ber(EbN0_dB, num_frames=400, N=200, approx=False):
"""単一RSC + BCJRのBERを測定"""
R = 0.5
sigma = np.sqrt(1.0 / (2 * R * 10 ** (EbN0_dB / 10)))
Lc = 2.0 / sigma ** 2
errors, total = 0, 0
for _ in range(num_frames):
u = np.random.randint(0, 2, N)
sb, pb = encode_rsc(u)
ys = (1 - 2 * sb) + sigma * np.random.randn(N)
yp = (1 - 2 * pb) + sigma * np.random.randn(N)
La = np.zeros(N)
la, lb, lg = bcjr_full(ys, yp, La, Lc, approx=approx)
L, _ = compute_llr(la, lb, lg, ys, La, Lc, approx=approx)
u_hat = (L < 0).astype(int)
errors += np.sum(u_hat != u)
total += N
return errors / total
np.random.seed(1)
snr_range = np.arange(0, 5.5, 1.0)
ber_logmap = [simulate_ber(s, approx=False) for s in snr_range]
ber_maxlog = [simulate_ber(s, approx=True) for s in snr_range]
plt.figure(figsize=(8, 5))
plt.semilogy(snr_range, [max(b, 1e-5) for b in ber_logmap], 'o-', label='BCJR log-MAP')
plt.semilogy(snr_range, [max(b, 1e-5) for b in ber_maxlog], 's--', label='BCJR max-log-MAP')
plt.xlabel('$E_b/N_0$ (dB)')
plt.ylabel('BER')
plt.title('BER of single-RSC BCJR decoding')
plt.legend()
plt.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.show()
まとめ