ビット誤り率曲線を実測するQPSKモデム

作るもの

本プロジェクトでは、1つのNumPyファイルで完全なQPSKモデムを作り、そのビット誤り率（BER）を測定します。送信機はランダムなビットをグレイ符号化したQPSKシンボルに割り当て、ルートレイズドコサイン（RRC）フィルタでパルス整形します。チャネルは信号を数サンプル遅延させ、白色ガウス雑音を加えます。受信機は整合フィルタを通し、正しいサンプリング時点で1シンボルにつき1サンプルを取り出し、実部または虚部の符号から各ビットを判定します。

続いてプログラムはEb/N0を0 dBから10 dBまで掃引して誤りを数えます。測定したBERは、加法性白色ガウス雑音（AWGN）中のグレイ符号化QPSKの教科書的な理論値と並べて表示します。筆者らの実行では、すべての点で両者は誤り数の統計的なばらつきの範囲内で一致しました。6 dBでは、200,000ビット中482個の誤りが出ました（BERは2.41 × 10⁻³、理論値は2.39 × 10⁻³）。

最後のチェックでは、わざと誤ったサンプリング時点でサンプルを取ります。半シンボルのタイミング誤差があると、正常に動いていたリンクがコイン投げと同じになります。これがタイミング再生が必要な理由であり、実際の無線機で最初に必要になる機能です。

デジタル変調、パルス整形と整合フィルタ、雑音、SNR、通信路容量の各レッスンで、各要素を1つずつ解説しています。数式は変調とBERの早見表にまとめてあります。

必要なもの

• Python 3.10以降とNumPy。筆者らはApple Silicon搭載のMac miniで、Python 3.12とNumPy 2.5.4を使って実行しました。
• それ以外は不要です。SciPyもプロットライブラリも無線機も使いません。乱数はすべて、シードを固定した1つの生成器から得ます。
• 設計を読み、プログラムを実行し、各種の変更を試すのに約2時間。

設計

シンボルとエネルギー

2ビットの組ごとに、4点 (±1 ± j)/√2 のうち1点を選びます。1ビット目が実部の符号を、2ビット目が虚部の符号を決めます。これがグレイ符号化です。隣接する点どうしは1ビットしか違わないため、最も起こりやすいシンボル誤りで失うのは2ビットではなく1ビットです。どのシンボルもエネルギーはEs = 1なので、1ビットあたりEb = 0.5です。

シンボルレートは1 Msym/sで、1シンボルあたり8サンプルなので、サンプルレートは8 MHzです。これらの数値はラベルにしか影響しません。BERはEb/N0だけで決まります。

パルス整形と整合フィルタ

送信機は8サンプルごとのインパルスに各シンボルを載せます。このインパルス列を、ロールオフ率0.35、長さ10シンボル（81タップ）のRRCフィルタと畳み込みます。フィルタはエネルギーが1になるよう正規化します。占有帯域幅は (1 + 0.35) × 1 MHz = 1.35 MHz です。

受信機は同じフィルタを時間反転したものと畳み込みます。これが整合フィルタで、白色雑音中でサンプリング時点のSNRを最大にする線形フィルタです。2つのRRCフィルタを縦続接続するとレイズドコサインパルスになり、これは他のすべてのシンボル時点でゼロになります。そのため正しい時点でサンプリングすれば、見えるのは1シンボルだけで、符号間干渉はありません。プログラムでは雑音なしの実行でこれを確認します。

雑音のスケーリング

雑音レベルの設定は、初めて作るモデムで最も間違えやすい手順です。エネルギー1のシンボルとエネルギー1のフィルタを使うと、1シンボルあたりの送信エネルギーは、そのサンプル全体で合計してEs = 1になります。このとき、1サンプルあたり分散N0の複素白色雑音を加えれば正しいEs/N0が得られ、整合フィルタはこの比を変えません。コードでは、dB単位のEb/N0を、1シンボル2ビットとしてEs/N0 = 2 Eb/N0に換算します。そしてN0 = 1 / (Es/N0) と設定し、N0を実部と虚部に等分します。

既知のオフセットでのタイミング

このモデムはタイミングを再生しません。遅延を知っています。RRCフィルタは81タップなので、各フィルタは信号を40サンプル遅らせ、2つで80サンプル遅らせます。チャネルがさらに3サンプル加えます。受信機はインデックス83で最初のサンプルを取り、以降8サンプルごとに取ります。実際の受信機はこのオフセットを推定しなければなりません。最後のチェックでは、オフセットがずれるとBERがどれだけ急速に劣化するかを示します。

理論曲線

AWGN中でグレイ符号化QPSKを同期検波する場合、I軸とQ軸はそれぞれ独立したBPSK判定になります。ビット誤り確率はPb = Q(√(2 Eb/N0)) で、Qは標準正規分布の裾確率です。コードではmath.erfcからQ(x) = 0.5 erfc(x/√2) としてQを計算します。

コード

これをqpsk_modem.pyとして保存してください。フィルタ設計、変調器、チャネル、受信機、3つのチェックを1つのブロックにまとめています。

# qpsk_modem.py: QPSK modem with RRC pulse shaping, AWGN and a measured BER curve.
import math

import numpy as np

SPS = 8            # samples per symbol
ALPHA = 0.35       # RRC roll-off
SPAN = 10          # RRC length in symbols
CHANNEL_DELAY = 3  # extra delay in samples, known to the receiver
SYMBOL_RATE = 1e6  # 1 Msym/s, so the sample rate is 8 MHz
rng = np.random.default_rng(2026)


def rrc_taps(sps: int, alpha: float, span: int) -> np.ndarray:
    """Root-raised-cosine impulse response with unit energy."""
    t = np.arange(-span * sps // 2, span * sps // 2 + 1) / sps  # time in symbols
    h = np.zeros_like(t)
    for i, ti in enumerate(t):
        if abs(ti) < 1e-12:
            h[i] = 1.0 - alpha + 4 * alpha / math.pi
        elif abs(abs(4 * alpha * ti) - 1.0) < 1e-12:
            h[i] = (alpha / math.sqrt(2)) * ((1 + 2 / math.pi) * math.sin(math.pi / (4 * alpha))
                                            + (1 - 2 / math.pi) * math.cos(math.pi / (4 * alpha)))
        else:
            num = math.sin(math.pi * ti * (1 - alpha)) + 4 * alpha * ti * math.cos(math.pi * ti * (1 + alpha))
            den = math.pi * ti * (1 - (4 * alpha * ti) ** 2)
            h[i] = num / den
    return h / np.sqrt(np.sum(h ** 2))


def qpsk_map(bits: np.ndarray) -> np.ndarray:
    """Gray-coded QPSK, unit symbol energy: bit 0 -> +1, bit 1 -> -1 on each rail."""
    b = bits.reshape(-1, 2)
    return ((1 - 2 * b[:, 0]) + 1j * (1 - 2 * b[:, 1])) / np.sqrt(2)


def qpsk_demap(z: np.ndarray) -> np.ndarray:
    """Hard decisions: the sign of each rail."""
    return np.column_stack([(z.real < 0), (z.imag < 0)]).astype(np.int8).ravel()


def transmit(symbols: np.ndarray, h: np.ndarray) -> np.ndarray:
    up = np.zeros(len(symbols) * SPS, dtype=np.complex128)
    up[::SPS] = symbols                      # impulse train, one impulse per symbol
    return np.convolve(up, h)                # pulse shaping


def receive(rx: np.ndarray, h: np.ndarray, n_sym: int, offset: int) -> np.ndarray:
    mf = np.convolve(rx, h[::-1].conj())     # matched filter
    return mf[offset: offset + n_sym * SPS: SPS]  # one sample per symbol


def run_link(ebn0_db: float, n_bits: int, h: np.ndarray, timing_error: int = 0) -> tuple[int, int]:
    bits = rng.integers(0, 2, n_bits, dtype=np.int8)
    sym = qpsk_map(bits)
    tx = transmit(sym, h)
    tx = np.concatenate([np.zeros(CHANNEL_DELAY, dtype=np.complex128), tx])  # channel delay
    es_n0 = 10 ** (ebn0_db / 10) * 2         # 2 bits per symbol, Es = 1
    n0 = 1.0 / es_n0                         # noise power per complex sample
    noise = np.sqrt(n0 / 2) * (rng.standard_normal(len(tx)) + 1j * rng.standard_normal(len(tx)))
    rx = tx + noise
    offset = (len(h) - 1) + CHANNEL_DELAY + timing_error  # TX filter delay + MF delay + channel
    z = receive(rx, h, len(sym), offset)
    errors = int(np.count_nonzero(qpsk_demap(z) != bits))
    return errors, n_bits


def q_func(x: np.ndarray) -> np.ndarray:
    return 0.5 * np.array([math.erfc(v / math.sqrt(2)) for v in np.atleast_1d(x)])


h = rrc_taps(SPS, ALPHA, SPAN)
print(f"RRC: {len(h)} taps, alpha {ALPHA}, {SPS} samples/symbol, sample rate {SYMBOL_RATE * SPS / 1e6:.0f} MHz")

# Check 1: with no noise the RRC pair has no intersymbol interference.
bits = rng.integers(0, 2, 2000, dtype=np.int8)
sym = qpsk_map(bits)
z = receive(np.concatenate([np.zeros(CHANNEL_DELAY), transmit(sym, h)]), h, len(sym),
            len(h) - 1 + CHANNEL_DELAY)
print(f"noise-free: max |z - s| = {np.max(np.abs(z - sym)):.2e}, bit errors = {np.count_nonzero(qpsk_demap(z) != bits)}")

# Check 2: BER against Eb/N0, stop at 100 errors or 20 million bits.
print("Eb/N0 (dB)  bits       errors  BER sim    BER theory")
results = []
for ebn0_db in range(0, 11):
    errs, total = 0, 0
    while errs < 100 and total < 20_000_000:
        e, n = run_link(ebn0_db, 200_000, h)
        errs, total = errs + e, total + n
    ber = errs / total
    theory = q_func(np.sqrt(2 * 10 ** (ebn0_db / 10)))[0]
    results.append((ebn0_db, ber, theory))
    print(f"{ebn0_db:>6}      {total:>9}  {errs:>6}  {ber:.3e}  {theory:.3e}")

# Check 3: sample at the wrong instant (Eb/N0 = 6 dB).
print("timing error (samples)  BER at 6 dB")
for err in [0, 1, 2, 3, 4]:
    e, n = run_link(6.0, 400_000, h, timing_error=err)
    print(f"{err:>10}              {e / n:.3e}")

worst = max(abs(b / t - 1) for _, b, t in results)
print(f"largest relative gap between simulated and theoretical BER: {100 * worst:.1f} %")

掃引は、各Eb/N0で誤りを100個数えるか、2,000万ビットをシミュレーションした時点で停止します。100個の誤りで、BER推定値の相対標準誤差は約10%になります（誤り数の平方根の逆数）。2,000万ビットの上限により、10 dBの点も数秒以内に収まります。

実行

python -m pip install numpy
python qpsk_modem.py

筆者らのマシンでは、全体の実行に約8秒かかります。

テスト

筆者らの実行での出力を示します。シードを2026に固定し、同じNumPyバージョンを使えば、同じ数値が得られるはずです。

RRC: 81 taps, alpha 0.35, 8 samples/symbol, sample rate 8 MHz
noise-free: max |z - s| = 1.80e-02, bit errors = 0
Eb/N0 (dB)  bits       errors  BER sim    BER theory
     0         200000   15857  7.928e-02  7.865e-02
     1         200000   11082  5.541e-02  5.628e-02
     2         200000    7277  3.639e-02  3.751e-02
     3         200000    4472  2.236e-02  2.288e-02
     4         200000    2520  1.260e-02  1.250e-02
     5         200000    1185  5.925e-03  5.954e-03
     6         200000     482  2.410e-03  2.388e-03
     7         200000     174  8.700e-04  7.727e-04
     8         600000     104  1.733e-04  1.909e-04
     9        3000000     100  3.333e-05  3.363e-05
    10       20000000      85  4.250e-06  3.872e-06
timing error (samples)  BER at 6 dB
         0              2.452e-03
         1              6.370e-03
         2              3.231e-02
         3              1.130e-01
         4              2.498e-01
largest relative gap between simulated and theoretical BER: 12.6 %

出力を上から順に見ていきます。

雑音なしの行。 ビット誤りはゼロで、受信サンプルと送信シンボルの最大距離は0.018です。このわずかな残差は、RRCフィルタを10シンボルで打ち切ったことによる符号間干渉です。フィルタを長くすると小さくなります。この値が0.5に近い場合やビット誤りがある場合は、サンプリングオフセットが間違っているか、フィルタが正規化されていません。

BERの表。 相対差が最大になるのは7 dBの12.6%で、この推定は174個の誤りに基づいています。この誤り数では標準誤差は約7.6%なので、差は標準誤差の2倍以内です。10 dBでは85個の誤りで標準誤差は約11%となり、測定BERは理論値より10%高くなっています。曲線が理論値から一定の距離だけ右にずれている場合は、雑音のスケーリングを確認してください。N0が2倍違うと、曲線全体が3 dBずれます。

タイミングの掃引。 6 dBでは、BERは正しいタイミングでの2.45 × 10⁻³から、1サンプル（8分の1シンボル）の誤差で6.37 × 10⁻³に上がります。2サンプルの誤差では3.23 × 10⁻²に、4サンプル（半シンボル）では0.25になります。半シンボルずれると、サンプルは2つのシンボルの間に落ち、2つのパルスがほぼ同じだけ寄与します。

実際の無線機では

実際のリンクでは、このプログラムの3つの前提が崩れます。それぞれに受信機の新しいブロックが必要です。

• 未知のタイミング。 受信機は遅延を知りません。自己相関特性のよいプリアンブルでフレームの先頭を見つけます。次に、Gardner方式などのタイミング誤差検出器が、2つのクロックがずれていく中で最適なサンプリング時点を追跡します。これがないとどうなるかは、上のタイミング掃引が示しています。
• 搬送波オフセットと位相。 送信機と受信機の発振器は数ppm異なります。2.4 GHzでは10 ppmは24 kHzの搬送波周波数オフセットになり、1 Msym/sではコンステレーションが1シンボルあたり0.15 rad回転します。プリアンブルからの周波数推定と位相同期ループ（QPSKではCostasループ）でこれを除去します。QPSKには4重の位相不確定性もあり、既知のプリアンブルか差動符号化で解消します。
• 利得とアナログフロントエンド。 受信レベルは未知で、変化もします。QPSKの硬判定には振幅は不要ですが、高次のQAMでは必要なので、実際の受信機は自動利得制御を備えています。ミキサによるDCオフセットやIQインバランスもコンステレーションを歪ませます。

このモデムを無線で動かすには、送信サンプルをcomplex64としてファイルに書き出し、ソフトウェア無線（SDR）の送信機から8 MHzで再生します。受信サンプルを記録し、上の3つのブロックを追加したプログラムの受信側に入力します。使用する周波数の免許を持っていない限り、ダミーロードに向けるか、減衰器を入れたケーブル経由で行ってください。ファイル側はIQの記録と読み込みのレッスン、無線機側は初めてのGNU Radioフローグラフのレッスンで扱っています。GNU Radioのdigitalモジュールには、最初の2項目を実装した既製のシンボル同期ブロックとCostasループブロックがあります。

試してみること

1. ロールオフ率を変える。 ALPHA = 0.2、次にALPHA = 0.8に設定してください。正しいタイミングでのBERは変わりません。整合フィルタ出力のSNRはパルス形状に依存しないからです。一方、タイミング掃引の結果は変わります。筆者らの実行では、1サンプルのタイミング誤差でのBERは、ロールオフ率0.2で7.70 × 10⁻³、0.8で4.40 × 10⁻³でした（0.35では6.37 × 10⁻³）。ロールオフ率が小さいほどパルスの裾が大きくなるため、同じタイミング誤差でも損失が大きくなります。
2. フィルタを短くする。 SPAN = 4（33タップ）に設定してください。筆者らの実行では、雑音なしの残差は0.018から0.120に増えました。打ち切ったパルスがもはやナイキストパルスに近くないためです。曲線の高Eb/N0側が理論値から浮き上がります。10 dBでの測定BERは9.18 × 10⁻⁶で、理論値は3.87 × 10⁻⁶でした。
3. 整合フィルタを壊す。 受信フィルタを8サンプルの単純な移動平均に置き換えてください。曲線は右にずれます。何dBずれるか測定してください。
4. グレイ符号化をやめる。 ビットの組00、01、10、11を、円周上の順に4点へ割り当ててください。低Eb/N0ではBERが上がります。隣接シンボルへの誤りで、約半分の確率で2ビットを失うようになるためです。
5. 搬送波オフセットを加える。 整合フィルタの前で、受信サンプルにexp(j 2π f t) を掛けてください。f = 1 kHz、サンプルレートは8 MHzです。コンステレーションは1シンボルあたり0.36°回転し、45°回転したところでBERが破綻します。その後、簡単な判定帰還型の位相追従器を追加してください。
6. アイパターンを測る。 整合フィルタの出力を16サンプル（2シンボル）ずつの行に並べ替え、各列の値の広がりを出力してください。広がりが最も狭い列が最適なサンプリングオフセットです。目視でオフセットを見つけるときに使うのがアイパターンです。

まとめ

• QPSKモデムは、グレイマッピング、RRCパルス整形、整合フィルタ、符号判定器の4つの要素でできています。いずれもNumPyで数行に収まります。
• Eb/N0からの雑音のスケーリングは慎重に行ってください。エネルギー1のシンボルとフィルタを使う場合、QPSKでの複素サンプルあたりのN0は1 / (2 Eb/N0) です。
• 測定BERは、0 dBでの7.9 × 10⁻²から10 dBでの約4 × 10⁻⁶まで、誤り数から予想されるばらつきの範囲内でPb = Q(√(2 Eb/N0)) と一致しました。
• 誤りは十分な数を数えてください。100個の誤りで、BER推定値の相対不確かさは約10%です。
• タイミングは重要です。8分の1シンボルの誤差で6 dBでのBERは約2.6倍になり、半シンボルではリンクが使い物にならなくなりました。

関連ページ：変調と同期と推定のトピックハブ、および同じ雑音スケーリングをマルチパスチャネルで再利用するOFDMリンクのプロジェクト。