ゼロから作るOFDM

前回のレッスンは、ある問題を残して終わりました。マルチパスチャネル上の広帯域信号は、帯域内で20 dB以上変動する周波数応答を受けます。これは時間領域では符号間干渉です。シングルキャリアの受信機でこれを取り除くには、多数のタップを持つ等化器が必要です。直交周波数分割多重（OFDM）は別の方法を取ります。帯域を多数の狭いサブキャリアに分割し、各サブキャリアがフラットフェージングとみなせるほど狭くします。そしてこの分割を1回のFFTで行います。LTEの下りリンク（3GPP TS 36.211）と5G NR（TS 38.211）は、どちらもOFDMを使っています。このレッスンでは、numpy約60行で完全なOFDMリンクを構築します。各ステップを確認し、意図的に壊して試せるようにするためです。

1つの式で見る考え方

N個の複素シンボルX0〜X(N−1)を、サブキャリアごとに1つずつ用意します。送信機はこれらの逆DFTを計算し、N個の時間サンプルを順に送ります。これで各シンボルは、1サンプルあたりk/Nサイクルの周波数を持つ複素指数関数を変調します。N サンプルにわたって、これらの複素指数関数は互いに直交します。受信機はN サンプルを取り込み、DFTを計算します。チャネルがなければ、シンボルは正確に復元されます。

次に、長さLのインパルス応答hを持つマルチパスチャネルを加えます。受信サンプルは、送信ブロックとhの線形畳み込みです。DFTは巡回畳み込みを乗算に変えますが、線形畳み込みは巡回畳み込みではありません。各ブロックの末尾は次のブロックにはみ出します。そこで使うのがサイクリックプレフィックス（CP）です。各ブロックを送る前に、その末尾Ncp サンプルを先頭にコピーします。NcpがL − 1以上であれば、受信機は受信した各ブロックの先頭Ncp サンプルを捨てられます。残りは、元のブロックとhの巡回畳み込みに等しくなります。FFTの後、サブキャリアkは次の値を運びます。

Y_k = H_k X_k + W_k

ここでH_kはサブキャリアkにおけるチャネルの周波数応答、W_kは雑音です。シンボルを多数のサンプルにわたって広げていたチャネルが、N個の独立したフラットなチャネルになりました。受信機はそれぞれをH_kで割ります。これが1タップ等化器です。OFDMがブロードバンド無線の主流になった理由は、まさにこの点にあります。

DFTを使って周波数分割多重を構成する考え方は、少なくとも1971年のWeinsteinとEbertの論文にまでさかのぼります。この論文は「さらに読む」に挙げています。

numpyで作るリンク

パラメータは意図的に小さくしています。FFTは64点で、すべてのサブキャリアがQPSKシンボルを1つ運びます。サンプルレートは1 MHzなので、サブキャリア間隔は1 MHz / 64 = 15.625 kHzです。各OFDMシンボルの長さは、プレフィックスを除いて64 µsです。サイクリックプレフィックスは16 サンプル、つまり16 µsです。チャネルは直接波と3 µs・8 µs・12 µsの遅延波を持ちます。したがってチャネルは13 サンプルにわたり、プレフィックスに収まります。

各フレームは、既知のトレーニングシンボル4個で始まります。受信機は各サブキャリアで、受信した値を送信されたとわかっている値で割ってH_kを推定し、4つの推定値を平均します。これは最小二乗推定であり、チャネル推定の最も単純な形です。その後にデータシンボル200個が続きます。このコードブロックはEb/N0の各点につき40フレーム、つまり1 024 000 ビットを処理し、4通りの結果を求めます。CPあり・真のチャネル、CPあり・推定チャネル、CPなし・真のチャネル、そして理論BERです。

Ebにはプレフィックスに費やすエネルギーも含めます。送信機はそのエネルギーを負担するからです。理論曲線は、64本のサブキャリアそれぞれの利得|H_k|²でQPSKの誤り率を求め、それを平均したものです。受信機が捨てるプレフィックスのエネルギーの分だけ、係数64/80で割り引いています。

import numpy as np
from scipy.special import erfc

rng = np.random.default_rng(8)
n_sc = 64                       # subcarriers (FFT size)
n_cp = 16                       # cyclic prefix, samples
n_sym = 200                     # data OFDM symbols per frame
n_frames = 40                   # frames per Eb/N0 point (1 024 000 bits)
n_train = 4                     # training symbols for the channel estimate
fs = 1.0e6                      # sample rate 1 MHz, so subcarrier spacing = fs / n_sc = 15.625 kHz

# A fixed multipath channel (complex128): echoes at 3, 8 and 12 samples (3, 8 and 12 us), unit total power.
h = np.zeros(13, dtype=complex)
h[[0, 3, 8, 12]] = [1.0, 0.45 * np.exp(0.8j), 0.25 * np.exp(-1.2j), 0.15j]
h /= np.linalg.norm(h)
H = np.fft.fft(h, n_sc)         # true frequency response on each subcarrier

def qpsk(bits: np.ndarray) -> np.ndarray:
    """Gray-mapped QPSK with unit average energy."""
    b = bits.reshape(-1, 2)
    return ((1 - 2 * b[:, 0]) + 1j * (1 - 2 * b[:, 1])) / np.sqrt(2)

def ofdm_frame(ebn0_db: float, cp: int, estimate: bool) -> float:
    """One frame: n_train training symbols then n_sym data symbols. Returns the bit error rate."""
    bits = rng.integers(0, 2, 2 * n_sc * n_sym)
    X = qpsk(bits).reshape(n_sym, n_sc)                       # one row per OFDM symbol
    x = np.fft.ifft(X, axis=1) * np.sqrt(n_sc)                # IFFT: unit power per sample
    x = np.hstack([x[:, n_sc - cp:], x]).ravel()              # prepend the cyclic prefix, serialise
    pilots = qpsk(rng.integers(0, 2, 2 * n_sc))               # known training symbol, sent n_train times
    p = np.fft.ifft(pilots) * np.sqrt(n_sc)
    tx = np.concatenate([np.tile(np.concatenate([p[n_sc - cp:], p]), n_train), x])
    y = np.convolve(tx, h)[: tx.size]                         # multipath channel
    eb = (n_sc + cp) / (2 * n_sc)                             # energy per bit, CP included
    n0 = eb / 10 ** (ebn0_db / 10)
    y = y + np.sqrt(n0 / 2) * (rng.normal(size=y.size) + 1j * rng.normal(size=y.size))
    blocks = y.reshape(n_train + n_sym, n_sc + cp)[:, cp:]          # drop the CP of every symbol
    Y = np.fft.fft(blocks, axis=1) / np.sqrt(n_sc)            # FFT back to subcarriers
    H_hat = (Y[:n_train] / pilots).mean(axis=0) if estimate else H   # least-squares estimate, averaged
    X_hat = Y[n_train:] / H_hat                                     # one-tap equaliser per subcarrier
    bits_hat = np.empty(bits.size, dtype=int)
    bits_hat[0::2] = (X_hat.real < 0).ravel()
    bits_hat[1::2] = (X_hat.imag < 0).ravel()
    return float(np.mean(bits_hat != bits))

def ofdm_link(ebn0_db: float, cp: int, estimate: bool) -> float:
    return float(np.mean([ofdm_frame(ebn0_db, cp, estimate) for _ in range(n_frames)]))

def ber_theory(ebn0_db: float, cp: int) -> float:
    g = 10 ** (ebn0_db / 10) * n_sc / (n_sc + cp)             # CP energy is not used for detection
    return float(np.mean(0.5 * erfc(np.sqrt(g * np.abs(H) ** 2))))

g_db = 20 * np.log10(np.abs(H))
print(f"subcarrier spacing {fs / n_sc / 1e3:.3f} kHz, CP {n_cp / fs * 1e6:.0f} us, channel spans {h.size} samples ({h.size / fs * 1e6:.0f} us)")
print(f"channel gain across subcarriers: {g_db.min():.1f} dB to {g_db.max():.1f} dB")
print("Eb/N0 dB, CP known H, CP pilot LS, no CP known H, theory CP, AWGN")
for e in range(0, 21, 2):
    awgn = 0.5 * erfc(np.sqrt(10 ** (e / 10)))
    print(f"{e:2d}, {ofdm_link(e, n_cp, False):.2e}, {ofdm_link(e, n_cp, True):.2e}, "
          f"{ofdm_link(e, 0, False):.2e}, {ber_theory(e, n_cp):.2e}, {awgn:.2e}")

subcarrier spacing 15.625 kHz, CP 16 us, channel spans 13 samples (13 us)
channel gain across subcarriers: -12.6 dB to 4.0 dB
Eb/N0 dB, CP known H, CP pilot LS, no CP known H, theory CP, AWGN
 0, 1.34e-01, 1.63e-01, 1.18e-01, 1.34e-01, 7.86e-02
 2, 9.08e-02, 1.12e-01, 7.94e-02, 9.07e-02, 3.75e-02
 4, 5.65e-02, 7.08e-02, 5.05e-02, 5.64e-02, 1.25e-02
 6, 3.25e-02, 4.16e-02, 3.04e-02, 3.25e-02, 2.39e-03
 8, 1.74e-02, 2.29e-02, 1.94e-02, 1.75e-02, 1.91e-04
10, 8.97e-03, 1.17e-02, 1.28e-02, 8.86e-03, 3.87e-06
12, 4.25e-03, 5.52e-03, 9.32e-03, 4.19e-03, 9.01e-09
14, 1.80e-03, 2.35e-03, 7.34e-03, 1.78e-03, 6.81e-13
16, 6.42e-04, 8.89e-04, 6.10e-03, 6.26e-04, 2.27e-19
18, 1.39e-04, 2.43e-04, 5.30e-03, 1.61e-04, 1.40e-29
20, 2.64e-05, 4.69e-05, 4.89e-03, 2.47e-05, 1.04e-45

結果の読み方

真のチャネルを使ったシミュレーションのBERは、理論曲線によく一致します。10 dBでは8.86 × 10⁻³に対して8.97 × 10⁻³、20 dBでは2.47 × 10⁻⁵に対して2.64 × 10⁻⁵です。1タップ等化器は、式が示したとおりに動作しています。各サブキャリアは、それぞれ固有のSNRを持つフラットなチャネルとして振る舞います。

曲線から3つの効果が読み取れます。

周波数選択性はSNRを犠牲にします。 チャネル利得は、64本のサブキャリアにわたって−12.6 dBから+4.0 dBまで変動します。誤りの数は最も悪いサブキャリアで決まるため、同じBERを得るには、OFDMリンクはフラットなAWGNチャネルよりはるかに大きなEb/N0を必要とします。10 dBではBERが8.97 × 10⁻³であり、フラットなチャネルのQPSKなら3.87 × 10⁻⁶です。実際のシステムは、サブキャリアをまたいだ通信路符号化でこの損失の大部分を取り戻します。フェージングを受けたサブキャリアのビットを、強いサブキャリアのビットから復元するのです。符号化がなければ、周波数選択性チャネル上のOFDMは、前回のレッスンのフェージングチャネルとほぼ同じように振る舞います。

プレフィックスがないと誤りフロアが生じます。 低SNRでは、CPなしの受信機のほうがわずかに良好です。エネルギーの20%をプレフィックスに使わないからです。これは10·log10(80/64) = 0.97 dBの節約になります。約8 dBを超えると、関係は逆転します。各シンボルの末尾が次のシンボルに漏れ込みます（符号間干渉）。さらに畳み込みが巡回ではなくなるため、サブキャリア同士も漏れ込み合います（キャリア間干渉）。どちらも雑音を減らしても小さくなりません。BERは5 × 10⁻³付近で平坦になり、そこにとどまります。20 dBでのBERは4.89 × 10⁻³で、プレフィックスありでは2.64 × 10⁻⁵です。

チャネル推定のコストはわずかです。 最小二乗推定値は、それ自体が雑音を含みます。トレーニングシンボル4個を平均すると、推定値の雑音分散はデータの雑音の4分の1になり、損失は約1 dBです。10 dBでは、推定チャネルを使う受信機のBERが1.17 × 10⁻²、真のチャネルを使う場合は8.97 × 10⁻³です。実際のシステムはパイロットを時間と周波数にわたって分散配置し、その間の推定値を平滑化します。これはOFDMチャネル推定プロジェクトで扱います。

パラメータの選び方

設計規則はレッスン7から導かれます。プレフィックスは、チャネルのRMS遅延広がりだけでなく、遅延の全幅より長くなければなりません。RMS遅延広がりが300 nsの指数プロファイルでも、1 µsを超える位置に無視できないエネルギーが残ります。各サブキャリアがフラットフェージングとなるよう、サブキャリア間隔はコヒーレンス帯域幅より十分小さくする必要があります。また、OFDMシンボルはコヒーレンス時間に比べて短くなければなりません。そうでないとシンボルの途中でチャネルが変化し、サブキャリアの直交性が失われます。サブキャリア間隔に対して無視できない割合のドップラーシフトは、キャリア間干渉を引き起こします。このため、広いサブキャリア間隔は高速で移動するユーザーや高い搬送波周波数に適しています。

プレフィックスは純粋なオーバーヘッドです。長くすればより大きな遅延広がりに耐えられますが、無駄になる時間とエネルギーも増えます。64 サンプルに対して16 サンプルを使うこのトイリンクでは、20%を費やしています。実際のシステムは、より小さい割合とより多くのサブキャリアを使います。

5G NRのヌメロロジー

5G NRは、このトレードオフを調整可能にしています。3GPP TS 38.211（NR physical channels and modulation）は、μで番号付けされた一連のヌメロロジーを定義しています。サブキャリア間隔は15 × 2^μ kHzです。仕様のRelease 19版（V19.5.0）では、Table 4.2-1がμ = 0〜6、すなわち15、30、60、120、240、480、960 kHzを挙げています。すべてノーマルサイクリックプレフィックスを使い、μ = 2（60 kHz）では拡張プレフィックスも使えます。ノーマルプレフィックスでは、1スロットに14 OFDMシンボルが入ります。1 msのサブフレームには2^μ個のスロットが入ります。

仕様は、シンボル長とプレフィックス長を基本時間単位Tc = 1/(480 000 × 4096) sを単位として与えています（κ = 64）。有効シンボル長は2048κ × 2^−μ Tc、ノーマルプレフィックスは144κ × 2^−μ Tcです。半サブフレーム（0.5 ms）ごとの先頭シンボルには、16κ Tcが加算されます。次のコードブロックは、これらの式を評価します。

# NR normal-CP numerology from the formulas in 3GPP TS 38.211 V19.5.0 (Release 19), clauses 4.1 to 4.3.
tc = 1 / (480e3 * 4096)          # basic time unit Tc, seconds
kappa = 64
print("mu, SCS kHz, useful symbol us, normal CP us, CP of first symbol in each half-subframe us, slot ms")
for mu in range(7):
    scs = 15 * 2**mu
    n_u = 2048 * kappa * 2.0**-mu
    n_cp = 144 * kappa * 2.0**-mu
    n_cp_long = n_cp + 16 * kappa
    print(f"{mu}, {scs}, {n_u * tc * 1e6:.2f}, {n_cp * tc * 1e6:.3f}, {n_cp_long * tc * 1e6:.3f}, {1 / 2**mu:.4g}")

mu, SCS kHz, useful symbol us, normal CP us, CP of first symbol in each half-subframe us, slot ms
0, 15, 66.67, 4.688, 5.208, 1
1, 30, 33.33, 2.344, 2.865, 0.5
2, 60, 16.67, 1.172, 1.693, 0.25
3, 120, 8.33, 0.586, 1.107, 0.125
4, 240, 4.17, 0.293, 0.814, 0.0625
5, 480, 2.08, 0.146, 0.667, 0.03125
6, 960, 1.04, 0.073, 0.594, 0.01562

15 kHzでは、ノーマルプレフィックスは4.688 µs、有効シンボルは66.67 µsです。つまりプレフィックスはシンボルの約7%です。30 kHzではプレフィックスが半分の2.344 µsになりますが、それでもTR 38.901の「very long」遅延広がりの例である1000 nsを大きく上回ります。TS 38.101-2でFR2（24.25 GHz〜71 GHz）向けに規定されている120 kHzでは、プレフィックスは0.586 µsです。ネットワークは、キャリアと帯域幅パートごとにヌメロロジーを選びます。どのバンドがどのサブキャリア間隔をサポートするかは、TS 38.211には記載されていません。FR1（410 MHz〜7.125 GHz）については、TS 38.101-1（Release 19、V19.7.0）のTable 5.3.5-1が、各バンドのチャネル帯域幅とサブキャリア間隔を挙げています。OFDMとNRヌメロロジーのチートシートに、これらの値を1つの表にまとめています。

このトイリンクとNRにはいくつかの違いがあり、今後のレッスンやプロジェクトで扱います。NRは帯域端のサブキャリアをガードバンドとして空けます。チャネル推定用の参照信号は疎なグリッドに配置します。データチャネルはLDPC符号（TS 38.212）で符号化します。そして直交周波数分割多元接続（OFDMA）でサブキャリアをユーザー間で共有します。また、OFDMのピーク対平均電力比（PAPR）は高く、電力増幅器にとって問題になります。

演習

1. n_cp = 8に設定し、13 サンプルのチャネルがプレフィックスに収まらないようにします。コードブロックを実行し、CPなしの場合と比べてBER曲線がどこで平坦になるかを説明してください。
2. 固定チャネルを、指数型電力遅延プロファイルから生成したランダムなチャネルに置き換えてください（レッスン7の関数を再利用します）。100回のチャネル生成についてBERを平均し、結果をレッスン7のレイリー曲線と比較してください。
3. n_trainを4から1および16に変更してください。それぞれについて14 dBでのBERを測定し、その変化を推定値の雑音と関連付けてください。
4. 2番目のコードブロックを使い、各μについてプレフィックスのオーバーヘッドを全シンボル長に対する割合として計算してください。長い先頭シンボルを除いて、すべてのμで同じ値になる理由を説明してください。

まとめ

• OFDMはIFFTで各サブキャリアに1つのシンボルを載せ、FFTでそれらを復元します。チャネルより長いサイクリックプレフィックスを使うと、チャネルはサブキャリアごとに1つの複素利得として作用します。
• 1タップ等化器は、各サブキャリアをそのチャネル利得で割ります。このシミュレーションは、計数誤差の範囲内で理論と一致しました。
• プレフィックスがないと、符号間干渉とキャリア間干渉により、リンクは5 × 10⁻³付近の誤りフロアに達します。
• プレフィックスは時間とエネルギーを消費します。64 + 16 サンプルのこの設計では、20%と0.97 dBです。
• 5G NRはサブキャリア間隔を15 × 2^μ kHzでスケーリングし、Release 19ではμ = 0〜6です。ノーマルプレフィックスはそれに比例して短くなり、15 kHzで4.688 µs、30 kHzで2.344 µsです。

さらに読む

• S. Weinstein and P. Ebert, “Data transmission by frequency-division multiplexing using the discrete Fourier transform”, IEEE Transactions on Communication Technology, vol. 19, no. 5, 1971：DFTで周波数分割多重を生成する方法を扱った初期の論文です。
• 3GPP TS 38.211, NR; Physical channels and modulation：5G NRのヌメロロジー、フレーム構造、OFDM信号生成を規定しています。
• J. G. Proakis and M. Salehi, Digital Communications, 5th edition, McGraw-Hill, 2008：「Multichannel and Multicarrier Systems」の章を参照してください。
• A. Goldsmith, Wireless Communications, Cambridge University Press, 2005：「Multicarrier Modulation」の章を参照してください。