QAMのpython実装

AI_qammod_demod.py

import numpy as np
import matplotlib.pyplot as plt

# =====================================
# Multi-QAM (16, 32, 64) Modem Functions
# =====================================

def bits_to_gray(n):
    """整数をGrayコードに変換"""
    return n ^ (n >> 1)

def gray_levels(n_bits):
    """Gray符号化されたPAMレベルを生成"""
    M = 2 ** n_bits
    levels = np.arange(-(M - 1), M, 2)
    bin_to_level = np.zeros(M, dtype=float)
    for i in range(M):
        gray_i = bits_to_gray(i)
        bin_to_level[i] = levels[gray_i]
    return bin_to_level

def make_qam_constellation(M, normalize=True):
    """
    任意M-QAMコンステレーション（Gray符号化）を生成
    対応：M = 16, 32, 64
    """
    if M == 16:
        n_I, n_Q = 2, 2  # 4x4
    elif M == 32:
        n_I, n_Q = 3, 2  # 8x4 (近似的な正方格子)
    elif M == 64:
        n_I, n_Q = 3, 3  # 8x8
    else:
        raise ValueError("M must be one of 16, 32, 64")

    I_levels = gray_levels(n_I)
    Q_levels = gray_levels(n_Q)
    const_points = []
    bit_labels = []

    for i in range(2**n_I):
        for q in range(2**n_Q):
            sym = I_levels[i] + 1j * Q_levels[q]
            bits_i = [int(b) for b in format(i, f'0{n_I}b')]
            bits_q = [int(b) for b in format(q, f'0{n_Q}b')]
            bit_labels.append(bits_i + bits_q)
            const_points.append(sym)

    const_points = np.array(const_points)
    bit_labels = np.array(bit_labels)

    # 正規化（平均電力 = 1）
    if normalize:
        Es = np.mean(np.abs(const_points)**2)
        const_points /= np.sqrt(Es)

    # ビット列→シンボル辞書を作成
    mapping = {tuple(bits.tolist()): sym for bits, sym in zip(bit_labels, const_points)}

    return const_points, bit_labels, mapping, n_I + n_Q

# ----------------------------------
# 変調・復調・BER・チャネル関数群
# ----------------------------------

def qam_modulate(bits, M):
    """ビット列をM-QAMシンボルに変換"""
    const_points, bit_labels, mapping, k = make_qam_constellation(M)
    bits = np.asarray(bits).astype(int)
    assert bits.size % k == 0, f"Bits length must be multiple of {k} for {M}-QAM"
    n_sym = bits.size // k
    syms = np.zeros(n_sym, dtype=complex)
    for i in range(n_sym):
        btuple = tuple(bits[i*k:(i+1)*k].tolist())
        syms[i] = mapping[btuple]
    return syms, const_points, bit_labels, k

def qam_demodulate(symbols, const_points, bit_labels):
    """M-QAMをハード判定復調"""
    dists = np.abs(symbols.reshape(-1, 1) - const_points.reshape(1, -1))
    idx = np.argmin(dists, axis=1)
    bits_out = bit_labels[idx].reshape(-1)
    return bits_out.astype(int)

def awgn_channel(symbols, EbN0_dB, k):
    """AWGNチャネル"""
    EbN0 = 10**(EbN0_dB/10)
    N0 = 1.0 / (EbN0 * k)
    sigma2 = N0/2.0
    noise = np.sqrt(sigma2) * (np.random.randn(len(symbols)) + 1j*np.random.randn(len(symbols)))
    return symbols + noise

def ber(bits_tx, bits_rx):
    """BER計算"""
    bits_tx = np.asarray(bits_tx).astype(int)
    bits_rx = np.asarray(bits_rx).astype(int)
    assert bits_tx.shape == bits_rx.shape
    errors = np.sum(bits_tx != bits_rx)
    return errors / bits_tx.size, errors

def plot_constellation(const_points, rx_symbols=None, title="QAM Constellation"):
    """コンステレーション描画"""
    plt.figure(figsize=(6,6))
    plt.scatter(const_points.real, const_points.imag, marker='o', label='Ideal Points')
    if rx_symbols is not None:
        plt.scatter(rx_symbols.real, rx_symbols.imag, marker='.', s=10, alpha=0.3, label='Received')
    plt.axhline(0, color='gray', linewidth=0.5)
    plt.axvline(0, color='gray', linewidth=0.5)
    plt.xlabel("In-phase")
    plt.ylabel("Quadrature")
    plt.title(title)
    plt.legend()
    plt.grid(True)
    plt.gca().set_aspect('equal', adjustable='box')
    plt.show()

# -------------------------------
# 実行例（16, 32, 64-QAMテスト）
# -------------------------------
if __name__ == "__main__":
    for M in [16, 32, 64]:
        print(f"\n=== {M}-QAM Simulation ===")

        n_bits = 300000
        EbN0_dB = 10.0
        bits_tx = np.random.randint(0, 2, n_bits)

        # 変調
        symbols, const_points, bit_labels, k = qam_modulate(bits_tx, M)
        # チャネル
        rx = awgn_channel(symbols, EbN0_dB=EbN0_dB, k=k)
        # 復調
        bits_rx = qam_demodulate(rx, const_points, bit_labels)
        # BER
        ber_val, errors = ber(bits_tx, bits_rx)

        print(f"Eb/N0 = {EbN0_dB:.1f} dB, Bits = {n_bits}, Errors = {errors}, BER = {ber_val:.5e}")
        plot_constellation(const_points, rx_symbols=rx[:2000], title=f"{M}-QAM Constellation (Eb/N0={EbN0_dB} dB)")