畳み込み符号とビタビ復号 — トレリス線図で理解する最尤復号

import numpy as np

def convolutional_encode(bits, generators):
    """
    畳み込み符号化器
    bits: 入力ビット列 (0/1の配列)
    generators: 生成多項式のリスト（各多項式は係数のタプル）
      例: [(1,1,1), (1,0,1)] は (7,5)_8
    戻り値: 符号化されたビット列
    """
    n_outputs = len(generators)
    K = len(generators[0])  # 拘束長
    m = K - 1                # メモリ長

    # シフトレジスタ（初期状態は全ゼロ）
    register = [0] * m

    encoded = []
    for bit in bits:
        # 現在のレジスタ内容と入力ビットを結合
        state = [bit] + register

        # 各生成多項式に対して出力ビットを計算
        for g in generators:
            output_bit = 0
            for j in range(K):
                output_bit ^= (state[j] * g[j])
            encoded.append(output_bit)

        # シフトレジスタを更新（右シフト、新入力を先頭に）
        register = [bit] + register[:-1]

    return np.array(encoded)

# 動作確認: 具体例と一致するか検証
generators_75 = [(1, 1, 1), (1, 0, 1)]  # (7,5)_8
test_bits = np.array([1, 0, 1, 1, 0])
encoded = convolutional_encode(test_bits, generators_75)

print("入力:", test_bits)
print("符号化出力:", encoded)
print("出力ペア:", [(encoded[2*i], encoded[2*i+1]) for i in range(len(test_bits))])

import numpy as np
import matplotlib.pyplot as plt

def build_trellis(generators):
    """トレリスの遷移テーブルを構築"""
    n_outputs = len(generators)
    K = len(generators[0])
    m = K - 1
    n_states = 2 ** m

    # 遷移テーブル: next_state[state][input] = (next_state, output_bits)
    trellis = {}
    for state in range(n_states):
        trellis[state] = {}
        # 状態をビット列に変換
        state_bits = [(state >> (m - 1 - i)) & 1 for i in range(m)]

        for input_bit in [0, 1]:
            # レジスタの内容
            full_state = [input_bit] + state_bits

            # 出力計算
            outputs = []
            for g in generators:
                out = 0
                for j in range(K):
                    out ^= (full_state[j] * g[j])
                outputs.append(out)

            # 次の状態
            next_state_bits = [input_bit] + state_bits[:-1]
            next_state = 0
            for i, b in enumerate(next_state_bits):
                next_state |= (b << (m - 1 - i))

            trellis[state][input_bit] = (next_state, tuple(outputs))

    return trellis, n_states, m

# (7,5)_8 のトレリス構築
generators_75 = [(1, 1, 1), (1, 0, 1)]
trellis, n_states, m = build_trellis(generators_75)

# トレリス線図の可視化
fig, ax = plt.subplots(figsize=(14, 6))

T = 6  # 表示する時刻数
state_labels = ['00', '01', '10', '11']
y_positions = {0: 3, 1: 2, 2: 1, 3: 0}  # 状態の縦位置

# ノード（状態）の描画
for t in range(T + 1):
    for s in range(n_states):
        if t == 0 and s != 0:
            # 初期状態は00のみ
            ax.plot(t, y_positions[s], 'o', color='lightgray', markersize=15, zorder=3)
        else:
            ax.plot(t, y_positions[s], 'o', color='white', markersize=15, zorder=3,
                    markeredgecolor='black', markeredgewidth=1.5)
        if t == 0:
            ax.text(-0.5, y_positions[s], state_labels[s], ha='right', va='center',
                    fontsize=10, fontweight='bold')

# ブランチ（遷移）の描画
for t in range(T):
    for s in range(n_states):
        if t == 0 and s != 0:
            continue  # 初期は00のみ
        for inp in [0, 1]:
            ns, outputs = trellis[s][inp]
            output_str = ''.join(map(str, outputs))
            color = 'tab:blue' if inp == 0 else 'tab:red'
            linestyle = '-' if inp == 0 else '--'
            ax.plot([t, t+1], [y_positions[s], y_positions[ns]],
                    color=color, linestyle=linestyle, linewidth=1.5, alpha=0.6)
            # 出力ラベル
            mid_x = t + 0.5
            mid_y = (y_positions[s] + y_positions[ns]) / 2
            ax.text(mid_x, mid_y + 0.15, output_str, ha='center', va='bottom',
                    fontsize=7, color=color, alpha=0.8)

ax.set_xlabel('Time step t')
ax.set_ylabel('State')
ax.set_title('Trellis Diagram for (7,5)₈ Convolutional Code (K=3, R=1/2)')
ax.set_yticks([y_positions[s] for s in range(n_states)])
ax.set_yticklabels(state_labels)
ax.set_xlim(-1, T + 0.5)
ax.grid(True, alpha=0.2)

# 凡例
from matplotlib.lines import Line2D
legend_elements = [Line2D([0], [0], color='tab:blue', linestyle='-', label='Input = 0'),
                   Line2D([0], [0], color='tab:red', linestyle='--', label='Input = 1')]
ax.legend(handles=legend_elements, loc='upper right')

plt.tight_layout()
plt.savefig('trellis_diagram.png', dpi=150, bbox_inches='tight')
plt.show()

import numpy as np

def viterbi_decode(received, generators, traceback_depth=None):
    """
    ビタビ復号器（硬判定）
    received: 受信ビット列 (0/1の配列)
    generators: 生成多項式のリスト
    traceback_depth: トレースバック深さ（Noneなら全長）
    戻り値: 復号されたビット列
    """
    n_outputs = len(generators)
    K = len(generators[0])
    m = K - 1
    n_states = 2 ** m

    # トレリス構築
    trellis = {}
    for state in range(n_states):
        trellis[state] = {}
        state_bits = [(state >> (m - 1 - i)) & 1 for i in range(m)]
        for input_bit in [0, 1]:
            full_state = [input_bit] + state_bits
            outputs = []
            for g in generators:
                out = 0
                for j in range(K):
                    out ^= (full_state[j] * g[j])
                outputs.append(out)
            next_state_bits = [input_bit] + state_bits[:-1]
            next_state = 0
            for i, b in enumerate(next_state_bits):
                next_state |= (b << (m - 1 - i))
            trellis[state][input_bit] = (next_state, tuple(outputs))

    # 受信信号をn_outputsビットずつに分割
    n_steps = len(received) // n_outputs
    received_groups = received.reshape(n_steps, n_outputs)

    # パスメトリック（累積ハミング距離）
    INF = float('inf')
    path_metrics = np.full(n_states, INF)
    path_metrics[0] = 0  # 初期状態は00

    # 生き残りパスの履歴
    survivor_paths = [np.zeros((n_states, 0), dtype=int)]

    for t in range(n_steps):
        new_metrics = np.full(n_states, INF)
        new_predecessors = np.full(n_states, -1, dtype=int)
        new_inputs = np.full(n_states, -1, dtype=int)

        for state in range(n_states):
            if path_metrics[state] == INF:
                continue
            for input_bit in [0, 1]:
                next_state, outputs = trellis[state][input_bit]

                # ブランチメトリック（ハミング距離）
                branch_metric = np.sum(np.array(outputs) != received_groups[t])

                # 累積メトリック
                total_metric = path_metrics[state] + branch_metric

                # ACS: より良いパスを選択
                if total_metric < new_metrics[next_state]:
                    new_metrics[next_state] = total_metric
                    new_predecessors[next_state] = state
                    new_inputs[next_state] = input_bit

        path_metrics = new_metrics

        # 生き残りパスの更新
        old_paths = survivor_paths[-1]
        new_paths = np.zeros((n_states, old_paths.shape[1] + 1), dtype=int)
        for state in range(n_states):
            if new_predecessors[state] >= 0:
                pred = new_predecessors[state]
                if old_paths.shape[1] > 0:
                    new_paths[state, :-1] = old_paths[pred]
                new_paths[state, -1] = new_inputs[state]
        survivor_paths.append(new_paths)

    # トレースバック: 最小メトリックの状態から
    best_state = np.argmin(path_metrics)
    decoded = survivor_paths[-1][best_state]

    return decoded

# 動作確認: エラー訂正のテスト
generators_75 = [(1, 1, 1), (1, 0, 1)]

# 元の入力
original_bits = np.array([1, 0, 1, 1, 0])
# 符号化
encoded = convolutional_encode(original_bits, generators_75)
print("符号化出力:", encoded)

# 1ビットエラーを挿入
received = encoded.copy()
received[4] = 1 - received[4]  # 位置4のビットを反転
print("受信（エラーあり）:", received)

# ビタビ復号
decoded = viterbi_decode(received, generators_75)
print("復号結果:", decoded)
print("元の入力:", original_bits)
print("一致:", np.array_equal(decoded, original_bits))

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erfc

def simulate_conv_code_ber(generators, snr_db_range, n_info_bits=10000):
    """畳み込み符号+ビタビ復号のBERシミュレーション"""
    n_outputs = len(generators)
    R = 1.0 / n_outputs  # 符号化率

    ber_list = []
    for snr_db in snr_db_range:
        # Eb/N0からチャネルのSNRに変換
        snr_lin = 10 ** (snr_db / 10)
        # BPSK変調: ビットエネルギー Eb, コード化率Rを考慮
        # 符号化ビットあたりのエネルギー Ec = R * Eb
        # sigma^2 = N0/2 = Ec/(2*snr_lin*R)
        sigma = 1.0 / np.sqrt(2 * snr_lin * R)

        total_errors = 0
        total_bits = 0

        while total_bits < n_info_bits:
            # 情報ビット生成
            batch_size = min(100, n_info_bits - total_bits)
            info_bits = np.random.randint(0, 2, batch_size)

            # テイルビット追加（シフトレジスタをゼロに戻す）
            K = len(generators[0])
            tail_bits = np.zeros(K - 1, dtype=int)
            bits_with_tail = np.concatenate([info_bits, tail_bits])

            # 符号化
            encoded = convolutional_encode(bits_with_tail, generators)

            # BPSK変調: 0 -> -1, 1 -> +1
            modulated = 2.0 * encoded - 1.0

            # AWGN通信路
            noise = sigma * np.random.randn(len(modulated))
            received_soft = modulated + noise

            # 硬判定
            received_hard = (received_soft > 0).astype(int)

            # ビタビ復号
            decoded = viterbi_decode(received_hard, generators)

            # 情報ビット部分のみ比較
            decoded_info = decoded[:batch_size]
            bit_errors = np.sum(decoded_info != info_bits)
            total_errors += bit_errors
            total_bits += batch_size

        ber = total_errors / total_bits
        ber_list.append(max(ber, 1e-7))

    return np.array(ber_list)

# シミュレーション実行
snr_db = np.arange(0, 11, 1)

# (7,5) K=3
ber_75 = simulate_conv_code_ber(
    [(1, 1, 1), (1, 0, 1)], snr_db, n_info_bits=50000)

# (17,15) K=4
ber_1715 = simulate_conv_code_ber(
    [(1, 1, 1, 1), (1, 1, 0, 1)], snr_db, n_info_bits=50000)

# (35,23) K=5
ber_3523 = simulate_conv_code_ber(
    [(1, 1, 1, 0, 1), (1, 0, 0, 1, 1)], snr_db, n_info_bits=50000)

# 理論曲線
snr_fine = np.linspace(0, 10, 200)
snr_lin_fine = 10 ** (snr_fine / 10)

# 符号化なしBPSK
ber_bpsk = 0.5 * erfc(np.sqrt(snr_lin_fine))

fig, ax = plt.subplots(figsize=(10, 7))

ax.semilogy(snr_db, ber_75, 'bo-', markersize=6, linewidth=1.5,
            label='(7,5)₈ K=3, dfree=5 (sim)')
ax.semilogy(snr_db, ber_1715, 'rs-', markersize=6, linewidth=1.5,
            label='(17,15)₈ K=4, dfree=6 (sim)')
ax.semilogy(snr_db, ber_3523, 'g^-', markersize=6, linewidth=1.5,
            label='(35,23)₈ K=5, dfree=7 (sim)')
ax.semilogy(snr_fine, ber_bpsk, 'k--', linewidth=2, label='Uncoded BPSK (theory)')

ax.set_xlabel('Eb/N0 [dB]')
ax.set_ylabel('BER')
ax.set_title('BER of Convolutional Codes with Viterbi Decoding (Hard Decision)')
ax.legend()
ax.grid(True, which='both', alpha=0.3)
ax.set_ylim(1e-6, 1)
ax.set_xlim(0, 10)

plt.tight_layout()
plt.savefig('conv_code_ber.png', dpi=150, bbox_inches='tight')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 符号パラメータ
codes = [
    {'name': '(7,5)₈ K=3', 'dfree': 5, 'R': 0.5, 'color': 'tab:blue'},
    {'name': '(17,15)₈ K=4', 'dfree': 6, 'R': 0.5, 'color': 'tab:red'},
    {'name': '(35,23)₈ K=5', 'dfree': 7, 'R': 0.5, 'color': 'tab:green'},
    {'name': '(75,53)₈ K=6', 'dfree': 8, 'R': 0.5, 'color': 'tab:orange'},
    {'name': '(171,133)₈ K=7', 'dfree': 10, 'R': 0.5, 'color': 'tab:purple'},
]

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 左: 符号化利得 vs 拘束長
ax = axes[0]
K_values = [3, 4, 5, 6, 7]
dfree_values = [5, 6, 7, 8, 10]
# 硬判定の符号化利得（近似）: 10*log10(R*dfree)
coding_gain_hard = [10 * np.log10(0.5 * d) for d in dfree_values]
# 軟判定の符号化利得: 硬判定 + 約2dB
coding_gain_soft = [g + 2.0 for g in coding_gain_hard]

ax.plot(K_values, coding_gain_hard, 'bo-', markersize=8, linewidth=2,
        label='Hard decision')
ax.plot(K_values, coding_gain_soft, 'rs-', markersize=8, linewidth=2,
        label='Soft decision')
ax.set_xlabel('Constraint Length K')
ax.set_ylabel('Asymptotic Coding Gain [dB]')
ax.set_title('Coding Gain vs Constraint Length (R=1/2)')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xticks(K_values)

# 右: 復号の計算量 vs 拘束長
ax = axes[1]
m_values = [K - 1 for K in K_values]
states = [2 ** m for m in m_values]
# 1ビットあたりのACS演算回数
acs_per_bit = [2 * s for s in states]  # 各状態で2つの入力

ax.semilogy(K_values, acs_per_bit, 'ko-', markersize=8, linewidth=2)
ax.set_xlabel('Constraint Length K')
ax.set_ylabel('ACS Operations per Information Bit')
ax.set_title('Decoding Complexity vs Constraint Length')
ax.grid(True, which='both', alpha=0.3)
ax.set_xticks(K_values)

# 各点にラベル
for K, acs in zip(K_values, acs_per_bit):
    ax.annotate(f'{acs}', (K, acs), textcoords='offset points',
                xytext=(10, 5), fontsize=10)

plt.tight_layout()
plt.savefig('coding_gain_complexity.png', dpi=150, bbox_inches='tight')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erfc

def viterbi_decode_soft(received_soft, generators):
    """
    ビタビ復号器（軟判定）
    received_soft: 軟判定受信値（実数値の配列, BPSK: -1/+1 + noise）
    """
    n_outputs = len(generators)
    K = len(generators[0])
    m = K - 1
    n_states = 2 ** m

    # トレリス構築
    trellis = {}
    for state in range(n_states):
        trellis[state] = {}
        state_bits = [(state >> (m - 1 - i)) & 1 for i in range(m)]
        for input_bit in [0, 1]:
            full_state = [input_bit] + state_bits
            outputs = []
            for g in generators:
                out = 0
                for j in range(K):
                    out ^= (full_state[j] * g[j])
                outputs.append(out)
            next_state_bits = [input_bit] + state_bits[:-1]
            next_state = 0
            for i, b in enumerate(next_state_bits):
                next_state |= (b << (m - 1 - i))
            trellis[state][input_bit] = (next_state, tuple(outputs))

    n_steps = len(received_soft) // n_outputs
    received_groups = received_soft.reshape(n_steps, n_outputs)

    INF = float('inf')
    path_metrics = np.full(n_states, INF)
    path_metrics[0] = 0.0

    # 生き残りパスの履歴
    predecessors = []

    for t in range(n_steps):
        new_metrics = np.full(n_states, INF)
        new_pred = np.full(n_states, -1, dtype=int)
        new_inp = np.full(n_states, -1, dtype=int)

        for state in range(n_states):
            if path_metrics[state] == INF:
                continue
            for input_bit in [0, 1]:
                next_state, outputs = trellis[state][input_bit]

                # 軟判定ブランチメトリック（ユークリッド距離の二乗）
                expected = np.array([2.0 * o - 1.0 for o in outputs])  # BPSK: 0->-1, 1->+1
                branch_metric = np.sum((received_groups[t] - expected) ** 2)

                total = path_metrics[state] + branch_metric
                if total < new_metrics[next_state]:
                    new_metrics[next_state] = total
                    new_pred[next_state] = state
                    new_inp[next_state] = input_bit

        path_metrics = new_metrics
        predecessors.append((new_pred.copy(), new_inp.copy()))

    # トレースバック
    best_state = np.argmin(path_metrics)
    decoded = np.zeros(n_steps, dtype=int)
    state = best_state
    for t in range(n_steps - 1, -1, -1):
        pred, inp = predecessors[t]
        decoded[t] = inp[state]
        state = pred[state]

    return decoded

def simulate_conv_code_ber_both(generators, snr_db_range, n_info_bits=10000):
    """硬判定と軟判定の両方のBERをシミュレーション"""
    n_outputs = len(generators)
    R = 1.0 / n_outputs
    K = len(generators[0])

    ber_hard_list = []
    ber_soft_list = []

    for snr_db in snr_db_range:
        snr_lin = 10 ** (snr_db / 10)
        sigma = 1.0 / np.sqrt(2 * snr_lin * R)

        errors_hard = 0
        errors_soft = 0
        total_bits = 0

        while total_bits < n_info_bits:
            batch_size = min(100, n_info_bits - total_bits)
            info_bits = np.random.randint(0, 2, batch_size)
            tail_bits = np.zeros(K - 1, dtype=int)
            bits_with_tail = np.concatenate([info_bits, tail_bits])

            encoded = convolutional_encode(bits_with_tail, generators)
            modulated = 2.0 * encoded - 1.0
            noise = sigma * np.random.randn(len(modulated))
            received_soft = modulated + noise

            # 硬判定復号
            received_hard = (received_soft > 0).astype(int)
            decoded_hard = viterbi_decode(received_hard, generators)
            errors_hard += np.sum(decoded_hard[:batch_size] != info_bits)

            # 軟判定復号
            decoded_soft = viterbi_decode_soft(received_soft, generators)
            errors_soft += np.sum(decoded_soft[:batch_size] != info_bits)

            total_bits += batch_size

        ber_hard_list.append(max(errors_hard / total_bits, 1e-7))
        ber_soft_list.append(max(errors_soft / total_bits, 1e-7))

    return np.array(ber_hard_list), np.array(ber_soft_list)

# シミュレーション実行
snr_db = np.arange(0, 10, 1)
generators_75 = [(1, 1, 1), (1, 0, 1)]

ber_hard, ber_soft = simulate_conv_code_ber_both(generators_75, snr_db, n_info_bits=50000)

# 理論曲線
snr_fine = np.linspace(0, 10, 200)
snr_lin_fine = 10 ** (snr_fine / 10)
ber_bpsk = 0.5 * erfc(np.sqrt(snr_lin_fine))

fig, ax = plt.subplots(figsize=(10, 7))

ax.semilogy(snr_db, ber_hard, 'bo-', markersize=6, linewidth=1.5,
            label='(7,5)₈ Hard Decision')
ax.semilogy(snr_db, ber_soft, 'rs-', markersize=6, linewidth=1.5,
            label='(7,5)₈ Soft Decision')
ax.semilogy(snr_fine, ber_bpsk, 'k--', linewidth=2, label='Uncoded BPSK')

ax.set_xlabel('Eb/N0 [dB]')
ax.set_ylabel('BER')
ax.set_title('Hard vs Soft Decision Viterbi Decoding: (7,5)₈ Code')
ax.legend()
ax.grid(True, which='both', alpha=0.3)
ax.set_ylim(1e-6, 1)
ax.set_xlim(0, 10)

# 符号化利得を矢印で表示
target_ber = 1e-4
ax.axhline(y=target_ber, color='gray', linestyle=':', alpha=0.5)
ax.text(0.5, target_ber * 1.5, f'BER = {target_ber}', fontsize=9, color='gray')

plt.tight_layout()
plt.savefig('hard_vs_soft_viterbi.png', dpi=150, bbox_inches='tight')
plt.show()