QAM(直交振幅変調)の理論と性能を解説

QAM(Quadrature Amplitude Modulation: 直交振幅変調)は、搬送波の振幅と位相の両方を同時に変化させるディジタル変調方式です。PSKが位相のみを用いるのに対し、QAMは2次元の信号空間(振幅と位相)をフルに活用するため、同じ帯域幅でより多くの情報を伝送できます。

Wi-Fi(802.11ac/ax)、LTE/5G、ケーブルテレビ、ADSLなど、現代の高速ディジタル通信システムの多くがQAMを基盤としています。特に256QAMや1024QAMといった高次QAMは、高い帯域効率を実現する鍵となる技術です。

本記事の内容

  • QAMの原理(振幅+位相の2次元変調)
  • M-QAMの信号点配置(矩形QAM)
  • 平均シンボルエネルギーの計算
  • BER近似式の導出(最近接距離とQ関数)
  • グレイ符号化によるBER改善
  • 帯域効率とシャノン限界との比較
  • Pythonでの星座図・BER曲線・雑音シミュレーション

前提知識

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

QAMの原理

PSKの限界とQAMの発想

PSKでは、すべての信号点が同じ振幅(一定半径の円周上)に配置されます。$M$ が大きくなると信号点間距離が急激に小さくなり、雑音耐性が悪化します。

QAMの発想は、位相だけでなく振幅も変えて、IQ平面を2次元的に広く使うことです。これにより、信号点を格子状に配置でき、PSKに比べて信号点間距離を大きく保てます。

QAM信号の数学的表現

QAM信号は、I成分(同相成分)とQ成分(直交成分)の独立な振幅変調の和として表されます。

$$ \boxed{ s(t) = I(t) \cos(2\pi f_c t) – Q(t) \sin(2\pi f_c t) } $$

ここで $I(t)$ と $Q(t)$ はそれぞれ離散的な振幅値を取ります。複素表現では、

$$ s(t) = \text{Re}\!\left[ \left(I(t) + jQ(t)\right) e^{j2\pi f_c t} \right] $$

$I(t) + jQ(t)$ がIQ平面上の信号点を表し、これを複素シンボルと呼びます。

M-QAMの信号点配置

矩形QAM(Square QAM)

$M = 2^{2k}$($M = 4, 16, 64, 256, \ldots$)のとき、$\sqrt{M} \times \sqrt{M}$ の正方格子状に信号点を配置できます。これを矩形QAM(Square QAM)と呼びます。

$L = \sqrt{M}$ として、I成分とQ成分がそれぞれ以下の $L$ 個の値を取ります。

$$ I_i, Q_j \in \left\{ -(L-1)d, \, -(L-3)d, \, \ldots, \, (L-3)d, \, (L-1)d \right\} $$

ここで $d$ は信号点の間隔の半分であり、隣接信号点間の距離は $2d$ です。

$L$ 個の値を一般的に書くと、

$$ I_i = (2i – L + 1)d, \quad i = 0, 1, \ldots, L-1 $$

$$ Q_j = (2j – L + 1)d, \quad j = 0, 1, \ldots, L-1 $$

具体的な信号点配置

4QAM(= QPSK): $L = 2$、信号点は $(\pm d, \pm d)$ の4点

16QAM: $L = 4$、信号点は $\{-3d, -d, d, 3d\} \times \{-3d, -d, d, 3d\}$ の16点

64QAM: $L = 8$、信号点は $\{-7d, -5d, -3d, -d, d, 3d, 5d, 7d\} \times \{-7d, \ldots, 7d\}$ の64点

平均シンボルエネルギーの計算

一般式の導出

矩形M-QAMの平均シンボルエネルギーを計算します。各信号点のエネルギーは $|I_i + jQ_j|^2 = I_i^2 + Q_j^2$ に比例します。

$L = \sqrt{M}$ 個の等確率のI成分値について、I成分の平均二乗値は

$$ \overline{I^2} = \frac{1}{L} \sum_{i=0}^{L-1} \left[(2i – L + 1)d\right]^2 $$

$m = 2i – L + 1$ とおくと、$m$ は $-(L-1), -(L-3), \ldots, (L-3), (L-1)$ を取ります。これらは奇数で、$|m|$ は $1, 3, 5, \ldots, L-1$ です。

$$ \overline{I^2} = \frac{d^2}{L} \sum_{i=0}^{L-1} (2i – L + 1)^2 $$

この和を計算します。$k = i – (L-1)/2$ とおくと $2k$ の範囲は $-(L-1), -(L-3), \ldots, (L-1)$ です。

$$ \sum_{i=0}^{L-1}(2i – L + 1)^2 = \sum_{\ell=0}^{L-1}(2\ell – L + 1)^2 $$

この和は、奇数の二乗和として計算できます。$L$ 個の値 $m = -(L-1), -(L-3), \ldots, (L-1)$ について

$$ \sum m^2 = 2\sum_{k=0}^{L/2-1}(2k+1)^2 = 2\sum_{k=0}^{L/2-1}(4k^2 + 4k + 1) $$

$$ = 2\left[ 4 \cdot \frac{(L/2-1)(L/2)(L-1)}{6} + 4 \cdot \frac{(L/2-1)(L/2)}{2} + \frac{L}{2} \right] $$

これを整理すると($L$ が偶数の場合)、

$$ \sum m^2 = \frac{L(L^2 – 1)}{3} $$

したがって、

$$ \overline{I^2} = \frac{d^2}{L} \cdot \frac{L(L^2 – 1)}{3} = \frac{d^2(L^2 – 1)}{3} $$

対称性より $\overline{Q^2} = \overline{I^2}$ です。平均シンボルエネルギーは

$$ E_s = \overline{I^2} + \overline{Q^2} = \frac{2d^2(L^2 – 1)}{3} = \frac{2d^2(M – 1)}{3} $$

$$ \boxed{ E_s = \frac{2d^2(M – 1)}{3} } $$

ビットエネルギーとの関係

1シンボルには $\log_2 M$ ビットが含まれるので、

$$ E_b = \frac{E_s}{\log_2 M} = \frac{2d^2(M – 1)}{3\log_2 M} $$

逆に、$d$ を $E_b$ で表すと、

$$ d^2 = \frac{3\log_2 M \cdot E_b}{2(M – 1)} $$

最近接信号点間距離

隣接する信号点間の距離(最近接距離)は $d_{\min} = 2d$ です。

$$ d_{\min} = 2d = 2\sqrt{\frac{3\log_2 M \cdot E_b}{2(M – 1)}} = \sqrt{\frac{6\log_2 M \cdot E_b}{M – 1}} $$

BER近似式の導出

シンボル誤り率(SER)の導出

矩形M-QAMのSERを、最近接シンボルへの誤り確率から導出します。

各軸(I, Q)を独立に考えます。$L = \sqrt{M}$ 値のPAM(パルス振幅変調)の誤り確率をまず求めます。

内部の信号点(端でない信号点)は、両側に隣接点があるため、一方の隣接点への誤り確率が $Q(d_{\min}/\sqrt{N_0})$ であり、両側合わせると

$$ P_{\text{内部}} = 2Q\!\left(\frac{d}{\sqrt{N_0/2}}\right) = 2Q\!\left(\sqrt{\frac{2d^2}{N_0}}\right) $$

端の信号点は、片側のみに隣接点があるため、

$$ P_{\text{端}} = Q\!\left(\sqrt{\frac{2d^2}{N_0}}\right) $$

$L$ 個の信号点のうち、端は2個、内部は $L – 2$ 個です。平均すると、

$$ P_{\text{PAM}} = \frac{2(L-2)Q(\cdot) + 2 \cdot Q(\cdot)}{L} \cdot \frac{1}{1} = \frac{2(L-1)}{L} Q\!\left(\sqrt{\frac{2d^2}{N_0}}\right) $$

すなわち、各軸の誤り確率は

$$ P_{\text{axis}} = \frac{2(\sqrt{M} – 1)}{\sqrt{M}} Q\!\left(\sqrt{\frac{2d^2}{N_0}}\right) $$

I軸とQ軸は独立なので、シンボルが正しい確率は

$$ P_{\text{correct}} = (1 – P_{\text{axis}})^2 $$

SERは

$$ P_s = 1 – (1 – P_{\text{axis}})^2 = 2P_{\text{axis}} – P_{\text{axis}}^2 \approx 2P_{\text{axis}} \quad (\text{高SNRでは } P_{\text{axis}} \ll 1) $$

$d^2$ を $E_b$ で表す式を代入します。

$$ \frac{2d^2}{N_0} = \frac{3\log_2 M}{M – 1} \cdot \frac{E_b}{N_0} $$

したがって、

$$ P_s \approx 4 \cdot \frac{\sqrt{M} – 1}{\sqrt{M}} \cdot Q\!\left(\sqrt{\frac{3\log_2 M}{M – 1} \cdot \frac{E_b}{N_0}}\right) $$

ビット誤り率(BER)

グレイ符号化の下では、隣接シンボルへの誤りは1ビット誤りに対応するため、

$$ P_b \approx \frac{P_s}{\log_2 M} $$

$$ \boxed{ P_b^{\text{M-QAM}} \approx \frac{4(\sqrt{M} – 1)}{\sqrt{M} \log_2 M} \cdot Q\!\left(\sqrt{\frac{3\log_2 M}{M – 1} \cdot \frac{E_b}{N_0}}\right) } $$

各QAMについて具体的に書くと:

4QAM(= QPSK) ($M = 4$):

$$ P_b = Q\!\left(\sqrt{\frac{2E_b}{N_0}}\right) $$

16QAM ($M = 16$):

$$ P_b \approx \frac{3}{4} Q\!\left(\sqrt{\frac{4E_b}{5N_0}}\right) $$

64QAM ($M = 64$):

$$ P_b \approx \frac{7}{12} Q\!\left(\sqrt{\frac{2E_b}{7N_0}}\right) \cdot \frac{1}{1} \approx \frac{7}{12} Q\!\left(\sqrt{\frac{E_b}{7N_0} \cdot 2}\right) $$

ただし正確には、

$$ P_b^{\text{64QAM}} \approx \frac{4 \cdot 7}{8 \cdot 6} Q\!\left(\sqrt{\frac{3 \cdot 6}{63} \cdot \frac{E_b}{N_0}}\right) = \frac{7}{12} Q\!\left(\sqrt{\frac{2E_b}{7N_0}}\right) $$

グレイ符号化によるBER改善

グレイ符号化の原理

矩形QAMでは、I軸とQ軸のそれぞれに独立にグレイ符号を適用します。$L = \sqrt{M}$ 値の各軸について、隣接するPAMレベルが1ビットだけ異なるようにビットを割り当てます。

例えば16QAMの場合、I軸の4レベル $\{-3d, -d, d, 3d\}$ に対して、

$$ -3d \to 00, \quad -d \to 01, \quad d \to 11, \quad 3d \to 10 $$

Q軸も同様にグレイ符号化します。各シンボルのビット列は I軸のビット(MSB側)とQ軸のビット(LSB側)を連結したものになります。

グレイ符号化の効果

グレイ符号化なしの場合、隣接シンボルへの誤りで複数ビットが誤る可能性があり、BERが悪化します。グレイ符号化により、隣接誤りが常に1ビット誤りに対応するため、BERが $P_s / \log_2 M$ で近似できるようになります。

自然2進符号と比較すると、グレイ符号化は高SNR領域で約 $1/\log_2 M$ 倍のBER改善をもたらします。

帯域効率

M-QAMの帯域効率

シンボルレート $R_s$ [symbol/s] で帯域幅 $B \approx R_s$ [Hz](ナイキスト帯域幅)の場合、ビットレートは

$$ R_b = R_s \log_2 M $$

帯域効率は

$$ \eta = \frac{R_b}{B} = \log_2 M \quad [\text{bit/s/Hz}] $$

方式 $M$ 帯域効率 [bit/s/Hz]
4QAM (QPSK) 4 2
16QAM 16 4
64QAM 64 6
256QAM 256 8
1024QAM 1024 10

$M$ を増やすと帯域効率は対数的に向上しますが、同じBERを達成するために必要な $E_b/N_0$ が増加します。

シャノン限界との比較

シャノンの通信容量定理によれば、帯域幅 $B$ のAWGNチャネルで達成可能な最大ビットレートは

$$ C = B \log_2\!\left(1 + \frac{S}{N}\right) = B \log_2\!\left(1 + \frac{E_b R_b}{N_0 B}\right) $$

帯域効率 $\eta = C/B$ を最大とする条件は $R_b = C$ のとき、

$$ \eta = \log_2\!\left(1 + \eta \cdot \frac{E_b}{N_0}\right) $$

これを $E_b/N_0$ について解くと、

$$ \frac{E_b}{N_0} = \frac{2^\eta – 1}{\eta} $$

この式はシャノン限界曲線を与えます。$\eta \to 0$ の極限で $E_b/N_0 \to \ln 2 \approx -1.59$ dB となり、これが理論的な最低限の $E_b/N_0$ です。

実際のM-QAMの性能は、BER $= 10^{-5}$ の条件でシャノン限界から約 1.5〜3 dB の差があります。この差をシャノンギャップと呼びます。誤り訂正符号(ターボ符号、LDPC符号)を併用することで、この差を大幅に縮小できます。

M-PSKとM-QAMの性能比較

同じ $M$ に対して、M-QAMはM-PSKよりも高い雑音耐性を持ちます。

$M = 16$ の場合を比較します。

16PSK: 信号点は半径 $r$ の円周上に等間隔に16点配置されます。最近接距離は

$$ d_{\min}^{\text{16PSK}} = 2r\sin\!\left(\frac{\pi}{16}\right) \approx 0.390 r $$

16QAM: 同じ平均電力($r^2 = E_s$)のとき、最近接距離は

$$ d_{\min}^{\text{16QAM}} = 2d = 2\sqrt{\frac{E_s}{10}} \approx 0.632 \sqrt{E_s} $$

16QAMの方が最近接距離が大きいため、同じ $E_b/N_0$ に対してBERが小さくなります。これが高次変調でQAMが好まれる理由です。

ただし、QAMは振幅が一定でないため、非線形増幅器(飽和動作の増幅器)では歪みが生じます。この点ではPSKが有利です。通信システムの設計では、帯域効率、雑音耐性、増幅器の線形性を総合的に考慮する必要があります。

Pythonでの実装

各種QAMの星座図

import numpy as np
import matplotlib.pyplot as plt

def generate_qam_constellation(M):
    """
    矩形M-QAMの星座図を生成する
    """
    L = int(np.sqrt(M))
    assert L * L == M, "M must be a perfect square"

    # PAMレベル
    levels = 2 * np.arange(L) - (L - 1)  # -(L-1), ..., (L-1)

    # 全信号点
    I_grid, Q_grid = np.meshgrid(levels, levels)
    symbols = (I_grid + 1j * Q_grid).flatten()

    return symbols, levels

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

M_list = [4, 16, 64, 256]
titles = ['4QAM (QPSK)', '16QAM', '64QAM', '256QAM']

for idx, (M, title) in enumerate(zip(M_list, titles)):
    symbols, levels = generate_qam_constellation(M)
    ax = axes[idx]

    # 信号点のプロット
    ax.scatter(symbols.real, symbols.imag, s=50 if M <= 64 else 10,
               c='red', zorder=5, edgecolors='darkred',
               linewidths=0.5 if M <= 64 else 0.2)

    # グレイ符号ラベル(16QAM以下)
    if M <= 16:
        L = int(np.sqrt(M))
        bits_per_axis = int(np.log2(L))
        for sym in symbols:
            i_idx = int((sym.real + L - 1) / 2)
            q_idx = int((sym.imag + L - 1) / 2)
            gray_i = i_idx ^ (i_idx >> 1)
            gray_q = q_idx ^ (q_idx >> 1)
            label_i = format(gray_i, f'0{bits_per_axis}b')
            label_q = format(gray_q, f'0{bits_per_axis}b')
            label = label_i + label_q
            ax.annotate(label, (sym.real, sym.imag),
                        textcoords="offset points", xytext=(5, 5),
                        fontsize=7, color='blue')

    max_val = int(np.sqrt(M)) - 1
    margin = 1.5
    ax.set_xlim(-max_val - margin, max_val + margin)
    ax.set_ylim(-max_val - margin, max_val + margin)
    ax.set_aspect('equal')
    ax.grid(True, alpha=0.3)
    ax.axhline(y=0, color='k', linewidth=0.5)
    ax.axvline(x=0, color='k', linewidth=0.5)
    ax.set_xlabel('In-Phase (I)')
    ax.set_ylabel('Quadrature (Q)')
    ax.set_title(f'{title} ({int(np.log2(M))} bits/symbol)')

plt.suptitle('QAM Constellation Diagrams', fontsize=14)
plt.tight_layout()
plt.show()

$M$ が大きくなるにつれて信号点が密に配置され、隣接点間の距離が相対的に小さくなることが確認できます。

BER理論曲線の比較

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

def Q_func(x):
    """Q関数"""
    return 0.5 * erfc(x / np.sqrt(2))

# Eb/N0の範囲
EbN0_dB = np.linspace(0, 25, 300)
EbN0 = 10**(EbN0_dB / 10)

# M-QAM の BER理論式
def ber_mqam(M, EbN0):
    """矩形M-QAMのBER近似式"""
    L = np.sqrt(M)
    k = np.log2(M)
    coeff = 4 * (L - 1) / (L * k)
    arg = np.sqrt(3 * k / (M - 1) * EbN0)
    return coeff * Q_func(arg)

# M-PSK の BER近似式
def ber_mpsk(M, EbN0):
    """M-PSKのBER近似式"""
    k = np.log2(M)
    if M == 2:
        return Q_func(np.sqrt(2 * EbN0))
    else:
        return (2 / k) * Q_func(np.sqrt(2 * k * EbN0) * np.sin(np.pi / M))

fig, axes = plt.subplots(1, 2, figsize=(16, 7))

# QAM の BER比較
M_qam = [4, 16, 64, 256]
colors = ['blue', 'red', 'green', 'purple']

for M, c in zip(M_qam, colors):
    ber = ber_mqam(M, EbN0)
    axes[0].semilogy(EbN0_dB, ber, color=c, linewidth=2,
                     label=f'{M}-QAM ({int(np.log2(M))} bit/s/Hz)')

axes[0].set_xlabel('$E_b/N_0$ [dB]', fontsize=12)
axes[0].set_ylabel('Bit Error Rate (BER)', fontsize=12)
axes[0].set_title('BER vs $E_b/N_0$ for M-QAM', fontsize=13)
axes[0].legend(fontsize=10)
axes[0].grid(True, which='both', alpha=0.3)
axes[0].set_xlim(0, 25)
axes[0].set_ylim(1e-6, 1)
axes[0].axhline(y=1e-3, color='gray', linestyle=':', alpha=0.5)

# QAM vs PSK の比較(M=16)
M = 16
ber_qam16 = ber_mqam(16, EbN0)
ber_psk16 = ber_mpsk(16, EbN0)

axes[1].semilogy(EbN0_dB, ber_qam16, 'b-', linewidth=2, label='16-QAM')
axes[1].semilogy(EbN0_dB, ber_psk16, 'r--', linewidth=2, label='16-PSK')

# M=64
ber_qam64 = ber_mqam(64, EbN0)
ber_psk64 = ber_mpsk(64, EbN0)
axes[1].semilogy(EbN0_dB, ber_qam64, 'g-', linewidth=2, label='64-QAM')
axes[1].semilogy(EbN0_dB, ber_psk64, 'm--', linewidth=2, label='64-PSK')

axes[1].set_xlabel('$E_b/N_0$ [dB]', fontsize=12)
axes[1].set_ylabel('Bit Error Rate (BER)', fontsize=12)
axes[1].set_title('QAM vs PSK Comparison', fontsize=13)
axes[1].legend(fontsize=10)
axes[1].grid(True, which='both', alpha=0.3)
axes[1].set_xlim(0, 25)
axes[1].set_ylim(1e-6, 1)

plt.tight_layout()
plt.show()

# BER = 10^-5 に必要な Eb/N0
print("=== BER = 10^-5 に必要な Eb/N0 ===")
target = 1e-5
for M in [4, 16, 64, 256]:
    ber_vals = ber_mqam(M, EbN0)
    idx = np.argmin(np.abs(ber_vals - target))
    print(f"  {M:>4d}-QAM ({int(np.log2(M)):>2d} bit/s/Hz): "
          f"Eb/N0 = {EbN0_dB[idx]:.1f} dB")

QAMの方がPSKよりもBER性能が良いこと、$M$ が大きくなるにつれて必要な $E_b/N_0$ が増加することが確認できます。

シャノン限界との比較

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

def Q_func(x):
    """Q関数"""
    return 0.5 * erfc(x / np.sqrt(2))

def ber_mqam(M, EbN0):
    """矩形M-QAMのBER近似式"""
    L = np.sqrt(M)
    k = np.log2(M)
    coeff = 4 * (L - 1) / (L * k)
    arg = np.sqrt(3 * k / (M - 1) * EbN0)
    return coeff * Q_func(arg)

# シャノン限界曲線
eta_range = np.linspace(0.1, 12, 500)
EbN0_shannon = (2**eta_range - 1) / eta_range
EbN0_shannon_dB = 10 * np.log10(EbN0_shannon)

# 各M-QAMがBER=10^-5を達成するEb/N0と帯域効率
M_list = [4, 16, 64, 256, 1024]
target_ber = 1e-5
EbN0_test = np.linspace(0, 35, 5000)
EbN0_test_lin = 10**(EbN0_test / 10)

qam_points = []
for M in M_list:
    ber = ber_mqam(M, EbN0_test_lin)
    idx = np.argmin(np.abs(ber - target_ber))
    eta = np.log2(M)
    qam_points.append((EbN0_test[idx], eta, M))

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

# シャノン限界
ax.plot(EbN0_shannon_dB, eta_range, 'k-', linewidth=2,
        label='Shannon Limit')
ax.fill_betweenx(eta_range, -2, EbN0_shannon_dB, alpha=0.1, color='gray')
ax.text(2, 6, 'Achievable Region\n(Shannon Limit)',
        fontsize=10, color='gray', ha='center')

# QAM点のプロット
for ebn0, eta, M in qam_points:
    ax.plot(ebn0, eta, 'ro', markersize=10)
    ax.annotate(f'{M}-QAM', (ebn0, eta),
                textcoords="offset points", xytext=(10, 5),
                fontsize=10, fontweight='bold', color='red')

# シャノン限界との差(ギャップ)
print("=== シャノンギャップ (BER = 10^-5) ===")
for ebn0, eta, M in qam_points:
    shannon_ebn0_dB = 10 * np.log10((2**eta - 1) / eta)
    gap = ebn0 - shannon_ebn0_dB
    print(f"  {M:>4d}-QAM: Eb/N0 = {ebn0:.1f} dB, "
          f"Shannon = {shannon_ebn0_dB:.1f} dB, Gap = {gap:.1f} dB")
    # ギャップを矢印で表示
    ax.annotate('', xy=(shannon_ebn0_dB, eta),
                xytext=(ebn0, eta),
                arrowprops=dict(arrowstyle='<->', color='blue', lw=1.5))

ax.axvline(x=10*np.log10(np.log(2)), color='k', linestyle=':',
           alpha=0.5, label=f'$E_b/N_0$ min = {10*np.log10(np.log(2)):.2f} dB')

ax.set_xlabel('$E_b/N_0$ [dB]', fontsize=12)
ax.set_ylabel('Spectral Efficiency $\\eta$ [bit/s/Hz]', fontsize=12)
ax.set_title('M-QAM Performance vs Shannon Limit (BER = $10^{-5}$)',
             fontsize=13)
ax.legend(fontsize=10, loc='upper left')
ax.grid(True, alpha=0.3)
ax.set_xlim(-2, 30)
ax.set_ylim(0, 12)

plt.tight_layout()
plt.show()

各M-QAMの動作点がシャノン限界曲線の右側(達成不可能領域の外側)に位置していること、そしてシャノンギャップが存在することが確認できます。

雑音耐性シミュレーション(モンテカルロ法)

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

def Q_func(x):
    """Q関数"""
    return 0.5 * erfc(x / np.sqrt(2))

def ber_mqam(M, EbN0):
    """矩形M-QAMのBER近似式"""
    L = np.sqrt(M)
    k = np.log2(M)
    coeff = 4 * (L - 1) / (L * k)
    arg = np.sqrt(3 * k / (M - 1) * EbN0)
    return coeff * Q_func(arg)

def qam_monte_carlo(M, EbN0_dB_range, num_symbols=200000):
    """
    M-QAMのモンテカルロBERシミュレーション
    """
    L = int(np.sqrt(M))
    k = int(np.log2(M))
    bits_per_axis = int(np.log2(L))

    # PAMレベル
    levels = 2 * np.arange(L) - (L - 1)

    # 全星座点
    I_grid, Q_grid = np.meshgrid(levels, levels)
    constellation = (I_grid + 1j * Q_grid).flatten()

    # グレイ符号マッピング(各軸独立)
    def gray_encode(n, bits):
        return n ^ (n >> 1)

    # 平均シンボルエネルギー
    Es_avg = np.mean(np.abs(constellation)**2)

    ber_list = []

    for EbN0_dB in EbN0_dB_range:
        EbN0 = 10**(EbN0_dB / 10)
        Es = EbN0 * k  # Eb/N0 * bits_per_symbol
        # 正規化のためのスケール
        scale = np.sqrt(Es / Es_avg)

        # ランダムシンボル生成
        sym_idx = np.random.randint(0, M, num_symbols)
        tx = constellation[sym_idx] * scale

        # AWGN付加(N0 = 1に正規化)
        noise = (np.random.randn(num_symbols)
                 + 1j * np.random.randn(num_symbols)) / np.sqrt(2)
        rx = tx + noise

        # 最尤判定
        rx_rescaled = rx / scale
        # I軸とQ軸を独立に判定
        rx_I = np.round((rx_rescaled.real + L - 1) / 2).astype(int)
        rx_Q = np.round((rx_rescaled.imag + L - 1) / 2).astype(int)
        rx_I = np.clip(rx_I, 0, L - 1)
        rx_Q = np.clip(rx_Q, 0, L - 1)
        detected_idx = rx_Q * L + rx_I

        # ビット誤りのカウント(グレイ符号化)
        bit_errors = 0
        for s_tx, i_rx, q_rx in zip(sym_idx, rx_I, rx_Q):
            i_tx = s_tx % L
            q_tx = s_tx // L
            # グレイ符号でのビット比較
            gray_i_tx = gray_encode(i_tx, bits_per_axis)
            gray_q_tx = gray_encode(q_tx, bits_per_axis)
            gray_i_rx = gray_encode(i_rx, bits_per_axis)
            gray_q_rx = gray_encode(q_rx, bits_per_axis)
            # ビットごとの誤りカウント
            bit_errors += bin(gray_i_tx ^ gray_i_rx).count('1')
            bit_errors += bin(gray_q_tx ^ gray_q_rx).count('1')

        total_bits = num_symbols * k
        ber_list.append(bit_errors / total_bits)

    return np.array(ber_list)

# シミュレーション実行
EbN0_sim = np.arange(0, 22, 2)
EbN0_theory = np.linspace(0, 24, 300)
EbN0_theory_lin = 10**(EbN0_theory / 10)

results = {}
for M in [4, 16, 64]:
    print(f"{M}-QAM シミュレーション実行中...")
    results[M] = qam_monte_carlo(M, EbN0_sim, num_symbols=100000)

# プロット
fig, ax = plt.subplots(figsize=(10, 8))

colors = {4: 'blue', 16: 'red', 64: 'green'}

for M in [4, 16, 64]:
    c = colors[M]
    # 理論曲線
    ber_th = ber_mqam(M, EbN0_theory_lin)
    ax.semilogy(EbN0_theory, ber_th, f'{c[0]}-', linewidth=2,
                color=c, label=f'{M}-QAM (Theory)')

    # シミュレーション結果
    mask = results[M] > 0
    ax.semilogy(EbN0_sim[mask], results[M][mask], 'o',
                color=c, markersize=8,
                label=f'{M}-QAM (Simulation)')

ax.set_xlabel('$E_b/N_0$ [dB]', fontsize=12)
ax.set_ylabel('Bit Error Rate (BER)', fontsize=12)
ax.set_title('M-QAM BER: Theory vs Monte Carlo Simulation', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, which='both', alpha=0.3)
ax.set_xlim(0, 22)
ax.set_ylim(1e-5, 1)

plt.tight_layout()
plt.show()

シミュレーション結果が理論曲線と良好に一致し、$M$ の増加に伴うBER劣化が定量的に確認できます。

雑音環境下の星座図

import numpy as np
import matplotlib.pyplot as plt

def plot_noisy_qam(M, EbN0_dB, num_symbols=5000, ax=None):
    """雑音環境下のQAM星座図をプロットする"""
    L = int(np.sqrt(M))
    k = int(np.log2(M))

    # 星座図
    levels = 2 * np.arange(L) - (L - 1)
    I_grid, Q_grid = np.meshgrid(levels, levels)
    constellation = (I_grid + 1j * Q_grid).flatten()

    # 平均エネルギー
    Es_avg = np.mean(np.abs(constellation)**2)
    EbN0 = 10**(EbN0_dB / 10)
    Es = EbN0 * k
    scale = np.sqrt(Es / Es_avg)

    # ランダムシンボル生成
    sym_idx = np.random.randint(0, M, num_symbols)
    tx = constellation[sym_idx] * scale

    # AWGN付加
    noise = (np.random.randn(num_symbols)
             + 1j * np.random.randn(num_symbols)) / np.sqrt(2)
    rx = tx + noise

    if ax is None:
        fig_local, ax = plt.subplots(figsize=(6, 6))

    ax.scatter(rx.real, rx.imag, s=1, alpha=0.3, c='blue')

    # 理想星座点
    ideal = constellation * scale
    ax.scatter(ideal.real, ideal.imag, s=30, c='red',
               zorder=5, edgecolors='darkred', linewidths=0.5)

    max_range = np.max(np.abs(ideal)) * 1.5
    ax.set_xlim(-max_range, max_range)
    ax.set_ylim(-max_range, max_range)
    ax.set_aspect('equal')
    ax.grid(True, alpha=0.3)
    ax.set_xlabel('I')
    ax.set_ylabel('Q')

fig, axes = plt.subplots(3, 3, figsize=(15, 15))

configs = [
    (4, 5), (4, 10), (4, 15),
    (16, 5), (16, 10), (16, 15),
    (64, 10), (64, 15), (64, 20),
]

for idx, (M, snr) in enumerate(configs):
    row, col = divmod(idx, 3)
    plot_noisy_qam(M, snr, ax=axes[row, col])
    axes[row, col].set_title(f'{M}-QAM, $E_b/N_0$ = {snr} dB')

plt.suptitle('Noisy QAM Constellation Diagrams', fontsize=14, y=1.01)
plt.tight_layout()
plt.show()

$E_b/N_0$ が低い環境では受信信号点が広く散らばり、隣接シンボルとの判定が困難になることが確認できます。64QAMでは信号点が密に配置されるため、十分な $E_b/N_0$ がないと通信品質が著しく劣化します。

まとめ

本記事では、QAM(直交振幅変調)の理論と性能について解説しました。

  • QAM信号は $s(t) = I(t)\cos(2\pi f_c t) – Q(t)\sin(2\pi f_c t)$ で表される2次元変調方式である
  • 矩形M-QAMの平均シンボルエネルギーは $E_s = 2d^2(M-1)/3$ で計算される
  • BER近似式は $P_b \approx \frac{4(\sqrt{M}-1)}{\sqrt{M}\log_2 M} Q\!\left(\sqrt{\frac{3\log_2 M}{M-1} \cdot \frac{E_b}{N_0}}\right)$ で与えられる
  • グレイ符号化により、隣接シンボル誤りが1ビット誤りに限定され、BERが最小化される
  • 帯域効率は $\log_2 M$ [bit/s/Hz] であり、$M$ の増加で向上するが雑音耐性は低下する
  • 同じ $M$ に対して、QAMはPSKよりも高い雑音耐性を持つ(最近接距離が大きいため)
  • シャノン限界との差(シャノンギャップ)は約1.5〜3 dBで、誤り訂正符号との併用で縮小可能

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