パイロットによるチャネル推定を備えたOFDMリンク

作るもの

このプロジェクトでは、NumPyのファイル1つで直交周波数分割多重（OFDM）リンクを構築し、マルチパスにどこまで耐えられるかを測定します。送信機は48本のサブキャリアにQPSKデータを、さらに16本に既知のパイロットを配置します。続いて64点IFFTを実行し、16サンプルのサイクリックプレフィックスを付加します。チャネルは8タップのレイリーフェージングフィルタです。OFDMシンボルごとに新たに生成し、白色ガウス雑音を加えます。受信機はプレフィックスを除去してFFTを実行し、パイロット位置で最小二乗法によりチャネルを推定します。その推定値をデータサブキャリアに補間し、各サブキャリアをその推定値で除算します。

プログラムは4種類の受信機を比較します。1つはチャネルを正確に知っている受信機です。残る3つはパイロットからチャネルを推定し、それぞれ異なる補間方式を使います。各チャネル推定の平均二乗誤差（MSE）と各受信機のBERを、SNR 0 dBから30 dBの範囲で表示します。筆者らの実行では、チャネルを完全に知っている受信機は教科書どおりのレイリー曲線に一致しました。パイロット間の線形補間は、補間誤差のためにBER約1.4 × 10⁻²で頭打ちになりました。先頭8タップのみを残すDFTベースの補間は30 dBで7.42 × 10⁻⁴に達し、完全なチャネル情報の場合のBERの約1.5倍でした。

背景はゼロから作るOFDMとフェージングチャネルのレッスンで扱っています。QPSKモデムのプロジェクトは、同じ測定のシングルキャリア版です。

必要なもの

• Python 3.10以降とNumPy。筆者らはPython 3.12とNumPy 2.5.4で実行しました。
• ほかには何も要りません。処理全体がベクトル化されており、数秒で終わります。
• 設計を読み、プログラムを実行し、変更を試すのに約3時間。

設計

数値パラメータ

サンプルレートは20 MHz、FFTは64点なので、サブキャリア間隔は312.5 kHzです。プレフィックスを含む1シンボルは80サンプル、つまり4.0 µsです。サイクリックプレフィックスは16サンプル（0.8 µs）です。サブキャリア0から始めて4本に1本のサブキャリアが値1のパイロットを運び、パイロットは16本になります。残る48本のサブキャリアは、パイロットと同じく平均エネルギー1のQPSKを運びます。パイロットのオーバーヘッドは25%です。

チャネル

チャネルは1サンプル（50 ns）間隔の8タップで構成され、最大過剰遅延は350 nsです。各タップの平均電力は、1/eまでの減衰が2タップの指数関数に従って減少し、合計が1になるよう正規化しています。このときrms遅延広がりは82 nsです。各タップはこの平均電力をもつ複素ガウス変数なので、各タップの振幅はレイリー分布になります。プログラムはOFDMシンボルごとに新しいチャネルを生成します。したがって、各シンボルは時間方向の平均化なしに、そのシンボル自身のパイロットから推定しなければなりません。

チャネルがプレフィックスより短いため、各シンボルの末尾は次のシンボルのプレフィックス内に収まり、受信機はその部分を捨てます。FFTの後は、どのサブキャリアも単一の複素利得を受けるだけになり、Y[k] = H[k] X[k] + W[k]となります。これこそOFDMが存在する理由です。サブキャリアごとに1タップの等化器で、マルチパスチャネル全体を打ち消せます。

SNRとIFFTのスケール

コードでは√64でスケーリングしたユニタリなIFFTとFFTを使います。こうすると時間サンプルあたりの平均電力がサブキャリアあたりの平均電力に等しくなります。時間領域で分散N0の白色雑音は、どのサブキャリアでも分散N0になります。表中のSNRはサブキャリアあたりのEs/N0です。QPSKではEb/N0はその半分、つまり3 dB低い値です。

最小二乗推定

パイロットサブキャリアkでは受信機はX[k] = 1を知っているので、最小二乗（LS）推定値はY[k] / X[k]です。この推定値は不偏で、誤差はそのサブキャリア上の雑音そのもの、分散はN0です。LSだけでは、雑音を含むHのサンプルが16個得られるだけです。受信機は64本すべてのサブキャリアでHを必要とします。

補間

線形。 実部と虚部それぞれで、隣接するパイロット間を直線で結びます。サブキャリア61〜63は最後のパイロットより外側にあるので、コードは最後の区間を延長します。線形補間は計算量が小さい一方、Hの曲率を無視して直線で近似します。そのため雑音がなくても誤差が残ります。

DFT。 パイロットは4サブキャリアごとにHをサンプリングしています。チャネル長が16サンプル以下であれば、これらのサンプルの16点IDFTはインパルス応答の先頭16タップを正確に返します。そのタップを64までゼロパディングして64点FFTを実行すると、補間誤差のまったくないHが全サブキャリアで得られます。サイクリックプレフィックスがすでに有効なチャネルを16サンプルまでに制限しているので、この条件による追加の制約はありません。

打ち切り付きDFT。 チャネル長が最大8タップだと受信機が知っていれば、FFTの前にタップ8〜15をゼロにできます。これらのタップは雑音しか含まないので、MSEは半分になります。実際には、受信機はチャネル長を推定するか、安全な上限としてサイクリックプレフィックス長を使うことになります。

完全なチャネル情報での理論値

チャネルを完全に知っている場合、各データサブキャリアはフラットなレイリーフェージング利得を通したQPSKになります。1ビットあたりの平均BERは0.5 (1 − √(γ / (1 + γ)))で、γ = Eb/N0です。プログラムはチェーン全体の検証として、この値を測定値の横に表示します。

コード

これをofdm_link.pyとして保存してください。パラメータ、送信機、チャネル、受信機、3つの推定器、掃引を1つのブロックにまとめています。

# ofdm_link.py: OFDM over a multipath channel with comb pilots, LS channel estimation,
# interpolation and one-tap equalisation. Prints channel-estimate MSE and BER against SNR.
import numpy as np

N_FFT = 64            # subcarriers
N_CP = 16             # cyclic prefix in samples
FS = 20e6             # sample rate, so subcarrier spacing is 312.5 kHz
PILOT_STEP = 4        # a pilot on every 4th subcarrier
N_TAPS = 8            # channel length in samples (0.4 us at 20 MHz)
DECAY_TAPS = 2.0      # exponential power delay profile, 1/e decay in samples
N_SYM = 20_000        # OFDM symbols per SNR point
rng = np.random.default_rng(11)

pilot_idx = np.arange(0, N_FFT, PILOT_STEP)                 # 16 pilots: 0, 4, ..., 60
data_idx = np.setdiff1d(np.arange(N_FFT), pilot_idx)        # 48 data subcarriers
pdp = np.exp(-np.arange(N_TAPS) / DECAY_TAPS)
pdp /= pdp.sum()                                            # unit average channel power


def qpsk(bits: np.ndarray) -> np.ndarray:
    b = bits.reshape(-1, 2)
    return ((1 - 2 * b[:, 0]) + 1j * (1 - 2 * b[:, 1])) / np.sqrt(2)


def random_channels(n: int) -> np.ndarray:
    """n independent Rayleigh channels, one per OFDM symbol, shape (n, N_TAPS)."""
    g = (rng.standard_normal((n, N_TAPS)) + 1j * rng.standard_normal((n, N_TAPS))) / np.sqrt(2)
    return g * np.sqrt(pdp)


def ofdm_tx(grid: np.ndarray) -> np.ndarray:
    """grid (n, N_FFT) of subcarrier values -> time samples with CP, shape (n, N_FFT + N_CP)."""
    x = np.fft.ifft(grid, axis=1) * np.sqrt(N_FFT)          # unitary IFFT: power per sample = per subcarrier
    return np.concatenate([x[:, -N_CP:], x], axis=1)


def multipath(x: np.ndarray, h: np.ndarray) -> np.ndarray:
    """Convolve each symbol with its own channel. The last N_TAPS - 1 output samples would land in
    the next symbol's cyclic prefix, which the receiver discards, so they are dropped here."""
    y = np.zeros((x.shape[0], x.shape[1] + N_TAPS - 1), dtype=complex)
    for k in range(N_TAPS):
        y[:, k:k + x.shape[1]] += h[:, [k]] * x
    return y[:, :x.shape[1]]


def ofdm_rx(y: np.ndarray) -> np.ndarray:
    return np.fft.fft(y[:, N_CP:], axis=1) / np.sqrt(N_FFT)


def interp_linear(h_p: np.ndarray) -> np.ndarray:
    """Linear interpolation between pilots, linear extrapolation past the last pilot."""
    k = np.arange(N_FFT)
    re = np.array([np.interp(k, pilot_idx, v) for v in h_p.real])
    im = np.array([np.interp(k, pilot_idx, v) for v in h_p.imag])
    out = re + 1j * im
    tail = k > pilot_idx[-1]                                # subcarriers 61, 62, 63
    slope = (h_p[:, -1] - h_p[:, -2]) / PILOT_STEP
    out[:, tail] = h_p[:, [-1]] + slope[:, None] * (k[tail] - pilot_idx[-1])
    return out


def interp_dft(h_p: np.ndarray, n_keep: int) -> np.ndarray:
    """DFT interpolation: pilots -> impulse response, keep n_keep taps, zero-pad, back to 64 bins."""
    g = np.fft.ifft(h_p, axis=1)                            # len(pilot_idx) taps
    g[:, n_keep:] = 0
    full = np.zeros((h_p.shape[0], N_FFT), dtype=complex)
    full[:, :len(pilot_idx)] = g
    return np.fft.fft(full, axis=1)


def run(snr_db: float) -> dict[str, tuple[float, float]]:
    n_bits = N_SYM * len(data_idx) * 2
    bits = rng.integers(0, 2, n_bits, dtype=np.int8)
    grid = np.zeros((N_SYM, N_FFT), dtype=complex)
    grid[:, data_idx] = qpsk(bits).reshape(N_SYM, -1)
    grid[:, pilot_idx] = 1.0                                # known pilot value, same power as data
    h = random_channels(N_SYM)
    H = np.fft.fft(h, n=N_FFT, axis=1)                      # true frequency response
    n0 = 10 ** (-snr_db / 10)                               # SNR = Es/N0 per subcarrier
    y = multipath(ofdm_tx(grid), h)
    y += np.sqrt(n0 / 2) * (rng.standard_normal(y.shape) + 1j * rng.standard_normal(y.shape))
    Y = ofdm_rx(y)
    h_ls = Y[:, pilot_idx] / grid[:, pilot_idx]             # least squares at the pilots
    estimates = {"perfect CSI": H,
                 "LS + linear": interp_linear(h_ls),
                 "LS + DFT, 16 taps": interp_dft(h_ls, 16),
                 "LS + DFT, 8 taps": interp_dft(h_ls, N_TAPS)}
    out = {}
    for name, H_hat in estimates.items():
        Z = Y[:, data_idx] / H_hat[:, data_idx]             # one-tap zero-forcing equaliser
        rx_bits = np.column_stack([(Z.real < 0).ravel(), (Z.imag < 0).ravel()]).astype(np.int8).ravel()
        ber = np.count_nonzero(rx_bits != bits) / n_bits
        mse = np.mean(np.abs(H_hat - H) ** 2)
        out[name] = (mse, ber)
    return out


print(f"{N_FFT} subcarriers, CP {N_CP}, {len(pilot_idx)} pilots, {len(data_idx)} data, "
      f"spacing {FS / N_FFT / 1e3:.1f} kHz, symbol {(N_FFT + N_CP) / FS * 1e6:.1f} us")
tau = np.arange(N_TAPS) / FS
rms_ds = np.sqrt(np.sum(pdp * tau ** 2) - np.sum(pdp * tau) ** 2)
print(f"channel: {N_TAPS} taps, rms delay spread {rms_ds * 1e9:.0f} ns, max excess delay {(N_TAPS - 1) / FS * 1e9:.0f} ns")
# Interpolation error alone: feed the true channel at the pilots, no noise.
h = random_channels(N_SYM)
H = np.fft.fft(h, n=N_FFT, axis=1)
for name, H_hat in [("linear", interp_linear(H[:, pilot_idx])), ("DFT, 8 taps", interp_dft(H[:, pilot_idx], N_TAPS))]:
    print(f"noise-free interpolation MSE, {name}: {np.mean(np.abs(H_hat - H) ** 2):.2e}")

names = ["perfect CSI", "LS + linear", "LS + DFT, 16 taps", "LS + DFT, 8 taps"]
results = {snr_db: run(snr_db) for snr_db in range(0, 31, 3)}
print("channel estimate MSE")
print("SNR (dB)  " + "".join(f"{n:>19}" for n in names[1:]))
for snr_db, r in results.items():
    print(f"{snr_db:>8}  " + "".join(f"{r[n][0]:>19.2e}" for n in names[1:]))
print(f"BER ({N_SYM * len(data_idx) * 2} bits per point)")
print("SNR (dB)  " + "".join(f"{n:>19}" for n in names) + "   theory (Rayleigh)")
for snr_db, r in results.items():
    g = 10 ** (snr_db / 10) / 2                              # Eb/N0 = SNR / 2 for QPSK
    theory = 0.5 * (1 - np.sqrt(g / (1 + g)))
    print(f"{snr_db:>8}  " + "".join(f"{r[n][1]:>19.2e}" for n in names) + f"   {theory:>17.2e}")

すべての演算は2次元配列に対して行い、1行が1つのOFDMシンボルに対応します。Pythonのループは、8つのチャネルタップに対するものとSNR点に対するものだけです。

実行

python -m pip install numpy
python ofdm_link.py

掃引では、11点のSNRそれぞれで20,000個のOFDMシンボル、192万ビットをシミュレーションします。所要時間は約3秒です。

テスト

以下は、シードを11に固定した筆者らの実行結果です。

64 subcarriers, CP 16, 16 pilots, 48 data, spacing 312.5 kHz, symbol 4.0 us
channel: 8 taps, rms delay spread 82 ns, max excess delay 350 ns
noise-free interpolation MSE, linear: 2.71e-02
noise-free interpolation MSE, DFT, 8 taps: 2.50e-32
channel estimate MSE
SNR (dB)          LS + linear  LS + DFT, 16 taps   LS + DFT, 8 taps
       0             8.10e-01           1.00e+00           5.01e-01
       3             4.20e-01           5.03e-01           2.51e-01
       6             2.23e-01           2.52e-01           1.26e-01
       9             1.25e-01           1.25e-01           6.28e-02
      12             7.68e-02           6.31e-02           3.15e-02
      15             5.20e-02           3.15e-02           1.58e-02
      18             3.97e-02           1.58e-02           7.91e-03
      21             3.38e-02           7.93e-03           3.97e-03
      24             3.03e-02           3.99e-03           1.99e-03
      27             2.87e-02           2.00e-03           9.99e-04
      30             2.83e-02           9.99e-04           4.99e-04
BER (1920000 bits per point)
SNR (dB)          perfect CSI        LS + linear  LS + DFT, 16 taps   LS + DFT, 8 taps   theory (Rayleigh)
       0             2.12e-01           2.96e-01           3.11e-01           2.77e-01            2.11e-01
       3             1.47e-01           2.19e-01           2.33e-01           1.99e-01            1.47e-01
       6             9.20e-02           1.49e-01           1.57e-01           1.29e-01            9.21e-02
       9             5.29e-02           9.40e-02           9.63e-02           7.62e-02            5.31e-02
      12             2.87e-02           5.79e-02           5.43e-02           4.20e-02            2.88e-02
      15             1.51e-02           3.72e-02           2.93e-02           2.25e-02            1.51e-02
      18             7.74e-03           2.55e-02           1.53e-02           1.15e-02            7.74e-03
      21             3.97e-03           1.96e-02           7.76e-03           5.92e-03            3.92e-03
      24             1.92e-03           1.61e-02           3.88e-03           2.88e-03            1.98e-03
      27             9.75e-04           1.46e-02           1.99e-03           1.46e-03            9.95e-04
      30             4.94e-04           1.38e-02           1.02e-03           7.42e-04            4.99e-04

出力で次の点を確認してください。

補間のみ。 パイロットに真のチャネルを与え、雑音を加えない場合、線形補間ではMSE 2.71 × 10⁻²が残ります。DFT補間で残るのは2.50 × 10⁻³²で、これは丸め誤差です。この1行が、両方のグラフに現れるフロアを説明しています。

完全なチャネル情報。 測定したBERは、すべての点で統計的なばらつきの範囲内でレイリー理論の列に一致します。たとえば27 dBでは9.75 × 10⁻⁴に対して9.95 × 10⁻⁴です。一致しない場合は、雑音のスケーリングかIFFTの正規化が誤っています。曲線の下がり方が緩やかな点に注目してください。SNR 10 dBあたり、BERは1桁しか下がりません。フェージングチャネルでは、誤りは深いフェードに落ちたサブキャリアから生じます。

DFT推定器。 16タップDFT推定器のMSEはN0に等しく（0 dBで1.00、30 dBで9.99 × 10⁻⁴）、LSの雑音がそのまま通過しています。8タップを残すと、30 dBで4.99 × 10⁻⁴に半減します。完全なチャネル情報に対するBERの劣化は、高SNRでほぼ一定です。30 dBでは、完全な情報の4.94 × 10⁻⁴に対し、8タップで7.42 × 10⁻⁴、16タップで1.02 × 10⁻³です。

線形補間。 低SNRでは、線形補間が2つのパイロット間で雑音を平均化するため、16タップDFT推定器を上回ります（0 dBでMSE 0.810対1.00）。約12 dBを超えると補間誤差が支配的になります。電力をいくら増やしても、MSEは2.8 × 10⁻²付近、BERは1.4 × 10⁻²付近で横ばいになります。

実際の無線機では

無線で伝送する場合、チャネル推定の前に、このシミュレーションが省略している3つの処理が入ります。

• タイミング。 受信機は各OFDMシンボルの開始位置を見つける必要があります。プレフィックスのうち前のシンボルの末尾がかからない部分であれば、どこを開始点にしても動作します。これはシングルキャリアのリンクにはないOFDMの許容幅です。開始点が遅すぎると次のシンボルが漏れ込みます。一般的な手法は、繰り返されるプレフィックスとの相関、または前半と後半が同一のトレーニングシンボルとの相関です。
• 搬送波周波数オフセット。 312.5 kHz間隔の数分の一のオフセットでもサブキャリアの直交性が崩れ、1タップ等化器では除去できないキャリア間干渉が加わります。トレーニングシンボルの繰り返し部分間の位相回転からオフセットを推定し、FFTの前に時間領域で補正してください。用語集の搬送波周波数オフセットの項を参照してください。
• 変化するチャネル。 ここでは各シンボルが独立したチャネルをもちます。実際のチャネルはコヒーレンス時間にわたって変化します。そのため受信機は複数シンボルにわたってパイロットを平均し、パイロットは時間と周波数の両方向に分散して配置されます。当サイトのリファレンスライブラリにあるRelease 19版の3GPP TS 38.211（NR physical channels and modulation）は、5G NRがこの目的に使う復調用参照信号を7.4.1.1節で定義しています。

実際のシステムでは、DCサブキャリアと帯域端も空けておきます。その場合、DFT補間は不完全なパイロットの集合で動作し、打ち切りによって帯域端でエネルギーが漏れます。上記のきれいな結果は劣化すると考え、端のサブキャリアは別途テストしてください。

無線でリンクを試すには、時間サンプルをcomplex64で書き出し、ソフトウェア無線（SDR）からサンプルレート20 MHzで送信して（または無線機が対応するレートに合わせてすべての数値を縮小して）、結果を記録してください。その周波数の免許をもっていない場合は、減衰器付きのケーブルかダミーロードを使ってください。タイミングと周波数の推定用に、各フレームの先頭にトレーニングシンボルを追加します。ファイル側の扱いはIQの記録と読み込みのレッスンで扱っています。

試してみること

1. パイロット間隔を広げる。 PILOT_STEP = 8に設定します。Hをサンプリングするパイロットは8本になり、8タップのチャネルはDFT補間の限界ちょうどにあります。筆者らの実行では、DFT補間は依然として正確でした（30 dBでBER 1.05 × 10⁻³）。一方、線形補間の雑音なしMSEは0.244に、BERフロアは約9 × 10⁻²に上昇しました。続けてN_TAPS = 12も設定します。インパルス応答が8点IDFTでエイリアシングを起こし、DFT推定器の雑音なしMSEは3.19 × 10⁻²になり、BERは1.8 × 10⁻²付近で頭打ちになります。
2. 遅延広がりを長くする。 DECAY_TAPS = 4.0に設定します（rms遅延広がり104 ns）。Hは周波数方向により速く変化します。筆者らの実行では、線形補間の雑音なしMSEは2.71 × 10⁻²から6.80 × 10⁻²に上昇し、30 dBでのBERは2.93 × 10⁻²でした。チャネルはまだ8タップに収まるので、DFT推定器は影響を受けませんでした。
3. プレフィックスを短くしすぎる。 チャネルを8タップのままN_CP = 4に設定します。効果を確認するには、multipathを連続ストリームを処理するように変更し、あるシンボルの末尾が実際に次のシンボルにかかるようにします。完全なチャネル情報があってもBERは頭打ちになります。
4. MMSE等化。 ゼロフォーシングの除算をH* Y / (|H|² + N0)に置き換えます。サブキャリアごとに硬判定するQPSKでは、スケーリングが正の値なので判定結果は変わりません。違いが出てくるのは、振幅が情報を運ぶ16-QAMです。
5. 16-QAM。 サブキャリアあたり4ビットをマッピングし、4 × 4の格子で判定します。こうするとチャネル推定値の振幅が重要になり、推定器間の差が広がります。
6. 時間方向に平均する。 各チャネルを10シンボルの間保持し、その間のLS推定値を平均します。DFT推定器ではMSEが10分の1に下がるはずです。

まとめ

• チャネルより長いサイクリックプレフィックスにより、マルチパスはサブキャリアごとに1つの複素利得になります。サブキャリアごとに1回の除算で等化できます。
• パイロットでのLS推定は計算コストがかかりませんが、各推定値に雑音分散N0がそのまま残ります。
• 線形補間は補間誤差（ここではMSE 2.71 × 10⁻²）を残し、SNRによらずBERを1.4 × 10⁻²付近で頭打ちにします。
• チャネルがパイロット数より短ければ、DFT補間は正確です。今回の実行では、既知のチャネル長に打ち切ることで推定雑音が半減しました。
• まず完全なチャネル情報の曲線を測定してください。それがレイリーの式に一致すれば、送信機・チャネル・雑音のスケーリングは正しく、残る差は推定器によるものです。

背景：OFDMとマルチキャリアのトピックハブとOFDMとNRヌメロロジーの早見表。