OFDM実装の詳細 — IFFT/FFT・サイクリックプレフィクス・チャネル推定を実装する

5G NR のダウンリンクで動画をストリーミングしているとき、電波は無数のビルや壁に反射し、何十本ものパスを通って端末に到達します。受信機はこの「めちゃくちゃに混ざった信号」から、元のデータを正しく復元しなければなりません。では、受信機はチャネルの状態をどうやって知るのでしょうか？ そして、知った情報をどう使って信号を修復するのでしょうか？

その答えが、パイロットシンボルによるチャネル推定と周波数領域等化です。OFDMの基本原理（サブキャリアの直交性、IFFT/FFTによる変復調）はOFDMの原理と実装で解説しました。本記事はその続編として、実際のOFDMシステムを「動かす」ために不可欠な実装の詳細に踏み込みます。

OFDMの実装詳細を理解すると、以下のような分野が見えてきます。

• 無線通信規格の設計: 5G NR や Wi-Fi 6 のパラメータ設計（CP長、パイロット配置、サブキャリア間隔）の根拠がわかる
• ソフトウェア無線（SDR）: GNURadio 等でのOFDM送受信機の自作に必要な知識が揃う
• 通信システムシミュレーション: リンクレベルシミュレータを自分の手で構築できるようになる

本記事の内容

• OFDMの基本構造の復習と記号の整理
• IFFT/FFTによる変復調の実装上のポイント
• サイクリックプレフィクス（CP）がISIを除去し、チャネルを巡回畳み込みに変換する数学的証明
• CP長の設計指針
• パイロットシンボルの配置戦略（コムタイプ・ブロックタイプ）
• チャネル推定（LS推定・MMSE推定）の理論と導出
• 1タップ等化（ZF等化・MMSE等化）
• PAPR問題とその低減手法（クリッピング、SLM法）
• Python: OFDM送受信のフルシミュレーション（パイロット・チャネル推定込み）

前提知識

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

• OFDM（直交周波数分割多重）の原理と実装を解説
• 離散フーリエ変換(DFT)の理論と実装
• フェージングとダイバーシチ

OFDMの基本構造の復習

送受信の全体像

まず記号を整理しましょう。OFDMシステムの主要パラメータは以下のとおりです。

送信側の処理フローを改めて確認します。

1. 送信ビット列を変調（QPSK、16-QAMなど）し、周波数領域のデータシンボル $X[k]$ を生成
2. パイロットシンボルを所定のサブキャリアに挿入
3. $N$ 点IFFTで時間領域信号 $x[n]$ を生成
4. サイクリックプレフィクスを付加
5. D/A変換・アップコンバートして送信

受信側はこの逆順で処理を行います。CP除去、FFT、チャネル推定、等化、デマッピングです。

周波数領域の入出力関係

チャネルのインパルス応答を $h[n]$（$n = 0, 1, \dots, L-1$）とします。CP付加と除去が正しく行われた場合、FFT後の受信シンボル $Y[k]$ は

$$ Y[k] = H[k] \cdot X[k] + W[k], \quad k = 0, 1, \dots, N-1 $$

という非常にシンプルな形をとります。ここで $H[k] = \sum_{n=0}^{L-1} h[n] e^{-j2\pi kn/N}$ はチャネルの周波数応答、$W[k]$ はFFT後の雑音です。

この式はOFDMの最も重要な性質です。マルチパスチャネルを通過した信号が、周波数領域では各サブキャリアごとに独立なスカラー乗算になる、ということを意味しています。この性質のおかげで、複雑な時間領域等化器が不要になり、単純な割り算（1タップ等化）でデータを復元できます。

では、この性質がなぜ成り立つのか、その鍵を握るサイクリックプレフィクスの役割を数学的に厳密に証明しましょう。

サイクリックプレフィクスの数学的役割

直感的な理解

サイクリックプレフィクス（CP）とは、OFDMシンボルの末尾をコピーして先頭に付けたものです。なぜわざわざ冗長なデータを付加するのでしょうか？

日常的なアナロジーで考えましょう。トランプのカードを円形に並べて輪を作ったとします。どこから数え始めても、同じカードの並びが繰り返されます。これが「巡回」の意味です。OFDMシンボルの末尾を先頭にコピーすることで、チャネルの遅延がシンボルの境界をまたいでも、あたかもシンボルが巡回的に繰り返されているかのように見えます。

CPの役割は2つあります。

1. ISI（シンボル間干渉）の除去: 前のシンボルからの遅延成分がCP区間に収まる
2. 線形畳み込みを巡回畳み込みに変換: FFTによる1タップ等化を可能にする

それぞれを数学的に証明していきます。

CP付きOFDMシンボルの構成

$N$ 点のOFDMシンボル $x[n]$（$n = 0, 1, \dots, N-1$）に対して、CP長を $N_{\text{CP}}$ とすると、送信される信号は

$$ \tilde{x}[n] = x[(n)_N], \quad n = -N_{\text{CP}}, -N_{\text{CP}}+1, \dots, N-1 $$

です。ここで $(n)_N = n \bmod N$ は $N$ を法とした剰余を表します。具体的に書き下すと

$$ \underbrace{x[N-N_{\text{CP}}], \dots, x[N-1]}_{\text{CP（末尾のコピー）}}, \underbrace{x[0], x[1], \dots, x[N-1]}_{\text{本体}} $$

という $N + N_{\text{CP}}$ サンプルの信号列になります。

定理: CPによるISI除去

第 $m$ 番目のOFDMシンボルの送信信号を $\tilde{x}_m[n]$、チャネルを $h[n]$（$n = 0, \dots, L-1$）とします。

定理: $N_{\text{CP}} \geq L – 1$ のとき、CP除去後の第 $m$ 番目シンボルの受信信号は、第 $m$ 番目の送信シンボルのみに依存し、他のシンボルからの干渉（ISI）は含まれない。

証明:

チャネル出力 $r[n]$ は線形畳み込みで表されます。

$$ r[n] = \sum_{\ell=0}^{L-1} h[\ell] \, \tilde{x}[n – \ell] + w[n] $$

受信側ではCP区間を捨て、$n = 0, 1, \dots, N-1$ の部分のみを取り出します。このとき、$n – \ell$ の最小値は $n = 0, \ell = L – 1$ のとき $-(L-1)$ です。

もし $N_{\text{CP}} \geq L – 1$ ならば、$n – \ell \geq -(L-1) \geq -N_{\text{CP}}$ が成り立ちます。つまり、受信信号の本体部分（$n = 0, \dots, N-1$）に影響を与えるのは、CP区間と本体部分だけであり、前のOFDMシンボルの信号は一切含まれません。

これがISI除去の条件です。$\square$

定理: 巡回畳み込みへの変換

定理: CP除去後の受信信号 $y[n]$（$n = 0, 1, \dots, N-1$）は、送信シンボル $x[n]$ とチャネル $h[n]$ の $N$ 点巡回畳み込みに等しい。

$$ y[n] = (h \circledast_N x)[n] + w[n] $$

証明:

CP除去後（$n = 0, 1, \dots, N-1$）の受信信号は

$$ y[n] = \sum_{\ell=0}^{L-1} h[\ell] \, \tilde{x}[n – \ell] $$

です（雑音項は省略）。ここで $n \geq 0$ かつ $\ell \leq L – 1 \leq N_{\text{CP}}$ なので、$n – \ell \geq -N_{\text{CP}}$ です。

CPの定義より $\tilde{x}[n – \ell] = x[(n – \ell)_N]$ が成り立ちます。なぜなら、$n – \ell < 0$ の場合でも、CPが存在するため $\tilde{x}[n - \ell] = x[N + (n - \ell)] = x[(n - \ell) \bmod N]$ となるからです。

したがって

$$ y[n] = \sum_{\ell=0}^{L-1} h[\ell] \, x[(n – \ell) \bmod N] $$

$h[\ell] = 0$（$\ell \geq L$）と定義すれば、$\ell$ の範囲を $0$ から $N-1$ まで拡張しても和は変わりません。

$$ y[n] = \sum_{\ell=0}^{N-1} h[\ell] \, x[(n – \ell) \bmod N] = (h \circledast_N x)[n] $$

これはまさに $N$ 点巡回畳み込みの定義式です。$\square$

巡回畳み込み定理による周波数領域表現

巡回畳み込みの最も重要な性質は、DFT領域では要素ごとの積に変換されることです。

$$ \text{DFT}\{h \circledast_N x\} = H[k] \cdot X[k] $$

ここで $H[k] = \sum_{n=0}^{N-1} h[n] e^{-j2\pi kn/N}$ です。

$y[n]$ の $N$ 点DFTを $Y[k]$ とすると

$$ Y[k] = H[k] \cdot X[k] + W[k] $$

が得られます。これが冒頭で述べた周波数領域の入出力関係です。CPの付加・除去という一見単純な操作が、この強力な性質を生み出しているのです。

CPの数学的意味が明らかになったところで、実際の設計においてCP長をどう決めるかを考えましょう。

CP長の設計

基本条件

前節で示したとおり、ISI除去と巡回畳み込み化の条件は

$$ N_{\text{CP}} \geq L – 1 $$

です。ここで $L$ はチャネルのタップ数（インパルス応答の長さ）で、最大遅延拡がり $\tau_{\max}$ とサンプリング周期 $T_s$ から

$$ L = \left\lceil \frac{\tau_{\max}}{T_s} \right\rceil + 1 $$

と求まります。

周波数効率とのトレードオフ

CPは冗長データなので、帯域効率を低下させます。CP付きOFDMシンボルの帯域効率は

$$ \eta = \frac{N}{N + N_{\text{CP}}} $$

です。$N_{\text{CP}} = N/4$（LTEのNormal CPに相当）なら $\eta = 0.8$、つまり20%のオーバーヘッドです。

設計のポイント

CP長を決める際には、以下のポイントを考慮します。

1. チャネルの最大遅延拡がり: 屋内環境は $\tau_{\max} \approx 50 \sim 100\,\text{ns}$、都市部は $\tau_{\max} \approx 1 \sim 5\,\mu\text{s}$、山岳部は $\tau_{\max} \approx 10 \sim 20\,\mu\text{s}$
2. 帯域効率: CPが長すぎると帯域効率が大きく低下する
3. サブキャリア間隔との関係: $T = 1/\Delta f$ なので、$\Delta f$ を大きくすれば $T$ が短くなり、同じ $\tau_{\max}$ に対して相対的にCPのオーバーヘッドが大きくなる
4. ドップラー拡がりとの兼ね合い: $\Delta f$ を小さくするとドップラー拡がりに対するロバスト性が低下する（ICI増大）

5G NR では、用途（屋内/屋外、高移動体/低移動体）に応じて $\mu = 0, 1, 2, 3$ の4つのニューメロロジーが定義されており、$\Delta f = 15 \times 2^\mu\,\text{kHz}$ です。$\mu$ が大きいほどシンボル長が短くなり、高速移動体に適しますが、CPのオーバーヘッドは相対的に増加します。

ここまでで、チャネルの影響が周波数領域で $Y[k] = H[k] X[k] + W[k]$ というシンプルな形になることを示しました。しかし、受信機は $H[k]$ を事前に知りません。チャネルの周波数応答を推定するために使われるのがパイロットシンボルです。

パイロットシンボルの配置

パイロットシンボルとは

実際の通信システムでは、受信機はチャネルの状態を知りません。チャネル推定のために、送信側が受信側に既知のシンボル（パイロットシンボル）を埋め込みます。受信機は「送ったはずの値」と「実際に受け取った値」を比較することで、チャネルの影響を逆算します。

これは、カメラの色補正でグレーカードを使うのに似ています。既知の色（18%グレー）を撮影しておき、撮影結果と理想の色の差分から照明の影響を推定するわけです。

パイロットサブキャリアの集合を $\mathcal{P}$、パイロットシンボルの既知の値を $X_p[k]$（$k \in \mathcal{P}$）とします。パイロットサブキャリアでの受信値は

$$ Y_p[k] = H[k] \cdot X_p[k] + W[k], \quad k \in \mathcal{P} $$

です。$X_p[k]$ が既知なので、ここから $H[k]$ を推定できます。

配置パターン

パイロットの配置パターンは大きく3つに分類されます。

コムタイプ（Comb-type）配置

すべてのOFDMシンボル（時間方向）にパイロットを挿入し、周波数方向に $D_f$ サブキャリアごとに配置します。

• 利点: 時間変動するチャネルへの追従性が高い
• 欠点: 周波数方向のチャネル変動が大きい場合、推定精度が低下
• 用途: 高速移動体通信

ナイキストの標本化定理から、チャネルのインパルス応答の長さが $L$ タップの場合、周波数方向のパイロット間隔は

$$ D_f \leq \frac{N}{L} $$

を満たす必要があります。これを超えると、チャネルの周波数応答をエイリアシングなしに復元できません。

ブロックタイプ（Block-type）配置

$D_t$ シンボルごとに、すべてのサブキャリアにパイロットを挿入します。

• 利点: 周波数選択性チャネルでの推定精度が高い
• 欠点: シンボル間のチャネル変動に追従しにくい
• 用途: 低速移動・固定通信

時間方向のパイロット間隔は、チャネルの最大ドップラー周波数 $f_D$ とOFDMシンボル周期 $T_{\text{sym}} = (N + N_{\text{CP}})T_s$ に対して

$$ D_t \leq \frac{1}{2 f_D T_{\text{sym}}} $$

を満たす必要があります（ナイキスト条件）。

格子タイプ（Lattice-type / Scattered）配置

時間方向と周波数方向の両方に一定間隔でパイロットを散布します。LTEや5G NRはこの方式を採用しています。

• 利点: 時間変動と周波数選択性の両方に対応
• 欠点: 2次元補間が必要で計算量が増える

本記事では、実装のわかりやすさを優先してコムタイプ配置を用います。

パイロットの配置が決まれば、次はパイロット位置でのチャネル応答をどう推定するかです。LS推定とMMSE推定の2つの代表的手法を導出しましょう。

チャネル推定: LS推定

最小二乗（LS）推定の導出

最も基本的なチャネル推定がLS（Least Squares）推定です。パイロットサブキャリア $k \in \mathcal{P}$ での受信値は

$$ Y_p[k] = H[k] \cdot X_p[k] + W[k] $$

です。LS推定では、雑音の統計的性質を一切使わず、単に「受信値 ÷ 送信値」でチャネルを推定します。

LS推定量を求めるために、以下の二乗誤差を最小化します。

$$ J(\hat{H}) = \sum_{k \in \mathcal{P}} |Y_p[k] – \hat{H}[k] X_p[k]|^2 $$

$\hat{H}[k]$ で微分してゼロとおくと（$\hat{H}[k]$ は複素数なので共役微分を使います）

$$ \frac{\partial J}{\partial \hat{H}^*[k]} = -X_p^*[k](Y_p[k] – \hat{H}[k] X_p[k]) = 0 $$

整理すると

$$ \hat{H}[k] X_p[k] X_p^*[k] = Y_p[k] X_p^*[k] $$

$|X_p[k]|^2 = X_p[k] X_p^*[k]$ なので

$$ \hat{H}_{\text{LS}}[k] = \frac{Y_p[k]}{X_p[k]}, \quad k \in \mathcal{P} $$

非常にシンプルな結果です。行列形式で書くと

$$ \hat{\bm{H}}_{\text{LS}} = \bm{X}_p^{-1} \bm{Y}_p $$

ここで $\bm{X}_p = \text{diag}(X_p[k_1], X_p[k_2], \dots, X_p[k_{N_p}])$ はパイロットシンボルの対角行列、$\bm{Y}_p = (Y_p[k_1], \dots, Y_p[k_{N_p}])^T$ はパイロット位置の受信ベクトルです。

LS推定の推定誤差

LS推定量の推定誤差を評価しましょう。

$$ \hat{H}_{\text{LS}}[k] – H[k] = \frac{Y_p[k]}{X_p[k]} – H[k] = \frac{H[k] X_p[k] + W[k]}{X_p[k]} – H[k] = \frac{W[k]}{X_p[k]} $$

推定誤差の分散は

$$ \text{Var}(\hat{H}_{\text{LS}}[k]) = \frac{\sigma_w^2}{|X_p[k]|^2} $$

ここで $\sigma_w^2 = E[|W[k]|^2]$ は雑音の分散です。パイロットシンボルの電力が $|X_p[k]|^2 = 1$ のとき

$$ \text{MSE}_{\text{LS}} = \frac{\sigma_w^2}{1} = \frac{1}{\text{SNR}} $$

となります。LS推定は実装が極めて簡単ですが、雑音に対する耐性が低いことがわかります。SNRが低い環境では推定精度が大きく劣化します。

パイロット間の補間

コムタイプ配置の場合、パイロットは一部のサブキャリアにしか存在しません。データサブキャリアのチャネル応答は、パイロット位置の推定値から補間して求めます。代表的な補間方法は以下のとおりです。

線形補間: パイロット位置 $k_i$ と $k_{i+1}$ の間のサブキャリア $k$ に対して

$$ \hat{H}[k] = \hat{H}_{\text{LS}}[k_i] + \frac{k – k_i}{k_{i+1} – k_i} \left( \hat{H}_{\text{LS}}[k_{i+1}] – \hat{H}_{\text{LS}}[k_i] \right) $$

DFT補間: パイロット位置のチャネル推定値をIDFTで時間領域に変換し、チャネル長 $L$ を超える成分をゼロにして（ノイズ除去）、DFTで周波数領域に戻します。この方法は線形補間より高精度です。

LS推定はシンプルですが、低SNR環境では精度が不十分です。チャネルの統計的性質（相関構造）を活用すれば、さらに高精度な推定が可能です。それがMMSE推定です。

チャネル推定: MMSE推定

MMSE推定の導出

MMSE（Minimum Mean Square Error）推定は、LS推定にチャネルの事前情報（周波数相関行列と雑音分散）を組み合わせることで、推定精度を向上させます。

LS推定量 $\hat{\bm{H}}_{\text{LS}}$ をさらに線形変換して、MSEを最小化する推定量を求めます。

$$ \hat{\bm{H}}_{\text{MMSE}} = \bm{W} \hat{\bm{H}}_{\text{LS}} $$

ここで $\bm{W}$ は $N \times N_p$ の重み行列です。MSEは

$$ \text{MSE} = E\left[\|\bm{H} – \bm{W} \hat{\bm{H}}_{\text{LS}}\|^2\right] $$

です。ここで $\bm{H} = (H[0], H[1], \dots, H[N-1])^T$ は全サブキャリアのチャネル応答ベクトルです。

この最適化問題を解くために、ウィーナー・ホップフ方程式を使います。$\bm{W}$ で微分してゼロとおくと

$$ \bm{W}_{\text{opt}} = \bm{R}_{H\hat{H}} \bm{R}_{\hat{H}\hat{H}}^{-1} $$

ここで

• $\bm{R}_{H\hat{H}} = E[\bm{H} \hat{\bm{H}}_{\text{LS}}^H]$: 真のチャネルとLS推定値の相互相関行列
• $\bm{R}_{\hat{H}\hat{H}} = E[\hat{\bm{H}}_{\text{LS}} \hat{\bm{H}}_{\text{LS}}^H]$: LS推定値の自己相関行列

これらを具体的に計算しましょう。パイロット位置のみに注目した場合

$$ \hat{\bm{H}}_{\text{LS}} = \bm{H}_p + \bm{X}_p^{-1} \bm{W}_p $$

ここで $\bm{H}_p$ はパイロット位置の真のチャネル応答ベクトル、$\bm{W}_p$ はパイロット位置の雑音ベクトルです。

チャネルと雑音が無相関であることを使うと

$$ \bm{R}_{H\hat{H}} = E[\bm{H} (\bm{H}_p + \bm{X}_p^{-1}\bm{W}_p)^H] = \bm{R}_{HH_p} $$

次に、LS 推定値の自己相関行列を同様に展開します:

$$ \bm{R}_{\hat{H}\hat{H}} = E[(\bm{H}_p + \bm{X}_p^{-1}\bm{W}_p)(\bm{H}_p + \bm{X}_p^{-1}\bm{W}_p)^H] $$

雑音とチャネルの独立性より

$$ \bm{R}_{\hat{H}\hat{H}} = \bm{R}_{H_pH_p} + \sigma_w^2 (\bm{X}_p \bm{X}_p^H)^{-1} $$

パイロットシンボルの電力が等しく $|X_p[k]|^2 = \sigma_p^2$ の場合

$$ \bm{R}_{\hat{H}\hat{H}} = \bm{R}_{H_pH_p} + \frac{\sigma_w^2}{\sigma_p^2} \bm{I} $$

したがって、MMSE推定量は

$$ \hat{\bm{H}}_{\text{MMSE}} = \bm{R}_{HH_p} \left( \bm{R}_{H_pH_p} + \frac{1}{\text{SNR}} \bm{I} \right)^{-1} \hat{\bm{H}}_{\text{LS}} $$

となります。

MMSE推定の直感的解釈

MMSE推定の式を直感的に理解しましょう。

• 高SNR極限（$\text{SNR} \to \infty$）: $\frac{1}{\text{SNR}} \bm{I} \to \bm{0}$ なので $\hat{\bm{H}}_{\text{MMSE}} \to \bm{R}_{HH_p} \bm{R}_{H_pH_p}^{-1} \hat{\bm{H}}_{\text{LS}}$、つまりチャネルの相関構造に基づく補間になります
• 低SNR極限（$\text{SNR} \to 0$）: 雑音項が支配的になり、推定値はチャネルの事前平均に引き寄せられます

MMSE推定はチャネルの相関行列 $\bm{R}_{HH}$ を必要とします。これは、チャネルの遅延プロファイルから計算できます。指数減衰プロファイル $E[|h[\ell]|^2] = \sigma_h^2 e^{-\ell/\tau_{\text{rms}}}$ の場合、周波数領域の相関関数は

$$ R_H[\Delta k] = E[H[k] H^*[k + \Delta k]] = \sum_{\ell=0}^{L-1} E[|h[\ell]|^2] e^{-j2\pi \Delta k \ell / N} $$

です。これはフーリエ変換の関係にあるため、遅延プロファイルが既知であれば相関行列を構成できます。

LS推定とMMSE推定の比較

実際のシステムでは、LS推定 + DFT補間や、簡略化MMSE推定（相関行列の近似）が広く用いられます。

チャネル推定でチャネルの周波数応答 $\hat{H}[k]$ が得られたら、次はこの情報を使ってデータシンボルを復元する等化を行います。

1タップ等化

ゼロフォーシング（ZF）等化

最もシンプルな等化がZF等化です。各サブキャリアで推定チャネル応答の逆数を掛けます。

$$ \hat{X}_{\text{ZF}}[k] = \frac{Y[k]}{\hat{H}[k]} $$

ZF等化は $H[k] X[k]$ の影響を完全に除去しますが、$|H[k]|$ が小さい（深いフェージングの）サブキャリアでは、雑音が大きく増幅されるという問題があります。

等化後のSNRは

$$ \text{SNR}_{\text{ZF}}[k] = |H[k]|^2 \cdot \text{SNR} $$

となり、$|H[k]| \to 0$ のサブキャリアではSNRがゼロに近づきます。

MMSE等化

ZF等化の雑音増幅問題を改善するのがMMSE等化です。チャネル反転と雑音抑圧のバランスをとります。

等化後のシンボルは

$$ \hat{X}_{\text{MMSE}}[k] = \frac{H^*[k]}{|H[k]|^2 + \sigma_w^2 / \sigma_x^2} Y[k] $$

です。ここで $\sigma_x^2 = E[|X[k]|^2]$ はデータシンボルの平均電力です。

この式を導出しましょう。線形MMSE推定の一般形から、スカラーの場合に $\hat{X}[k] = c[k] Y[k]$ とおいて、$E[|X[k] – c[k] Y[k]|^2]$ を最小化する $c[k]$ を求めます。

$Y[k] = H[k] X[k] + W[k]$ を代入すると

$$ E[|X[k] – c[k] (H[k] X[k] + W[k])|^2] $$

$X[k]$ と $W[k]$ が無相関であることを使って展開します。

$$ = E[|X[k]|^2] |1 – c[k] H[k]|^2 + |c[k]|^2 E[|W[k]|^2] $$

$c[k]$ の共役で微分してゼロとおくと

$$ \sigma_x^2 H^*[k] (c[k] H[k] – 1)(-1) + c[k] \sigma_w^2 = 0 $$

$c[k]$ を含む項を左辺にまとめて整理すると、次のようになります:

$$ c[k] (\sigma_x^2 |H[k]|^2 + \sigma_w^2) = \sigma_x^2 H^*[k] $$

$c[k]$ について解くと

$$ c_{\text{MMSE}}[k] = \frac{\sigma_x^2 H^*[k]}{\sigma_x^2 |H[k]|^2 + \sigma_w^2} = \frac{H^*[k]}{|H[k]|^2 + 1/\text{SNR}} $$

最後の等号では $\text{SNR} = \sigma_x^2 / \sigma_w^2$ を使いました。

MMSE等化後のSNRは

$$ \text{SNR}_{\text{MMSE}}[k] = \frac{|H[k]|^2}{1/\text{SNR}} = |H[k]|^2 \cdot \text{SNR} $$

ではなく、正確には

$$ \text{SINR}_{\text{MMSE}}[k] = \frac{|H[k]|^2 \cdot \text{SNR}}{1 + |H[k]|^2 \cdot \text{SNR}} \cdot (|H[k]|^2 \cdot \text{SNR}) $$

ではなく、等化後の残留干渉と雑音を考慮した次の式になります。

$$ \text{SINR}_{\text{MMSE}}[k] = |H[k]|^2 \cdot \text{SNR} $$

実は出力SINRの値自体はZFと形式的に同じですが、MMSE等化の利点はバイアスと雑音増幅のトレードオフを最適に制御することにあります。深いフェード点では、ZFが大きな雑音を出力するのに対し、MMSEは雑音を抑える代わりに少量の残留ISIを許容します。全サブキャリアの平均BERで見ると、MMSEはZFよりも優れた性能を示します。

ここまでで、パイロットからチャネルを推定し、等化でデータを復元する一連の処理が揃いました。次に、OFDMのもう一つの課題であるPAPR問題と、その低減手法を見ていきましょう。

PAPR問題とその低減手法

PAPRの問題点を再確認する

OFDM信号の瞬時電力は、$N$ 個のサブキャリアの信号が重ね合わされた結果です。すべてのサブキャリアが同じ位相で加算されると、ピーク電力は平均電力の $N$ 倍に達します。PAPRは

$$ \text{PAPR} = \frac{\max_{0 \leq n \leq N-1} |x[n]|^2}{E[|x[n]|^2]} $$

で定義され、最悪ケースでは $\text{PAPR}_{\max} = N$ です。

なぜこれが問題なのでしょうか。送信機の電力増幅器（PA）には線形動作範囲の上限があります。PAPRが高いと、ピーク電力がPAの飽和領域に入り、信号が非線形歪みを受けます。非線形歪みは帯域外放射（隣接チャネルへの干渉）とシンボル間の干渉を引き起こします。

これを避けるために、PAの動作点を下げる（入力バックオフを大きくする）必要がありますが、それはPAの電力効率を低下させます。バッテリー駆動の端末ではこれが致命的です。

クリッピング

最もシンプルなPAPR低減手法がクリッピングです。時間領域信号のピークを閾値 $A$ で制限します。

$$ x_{\text{clip}}[n] = \begin{cases} x[n], & |x[n]| \leq A \\ A \cdot \frac{x[n]}{|x[n]|}, & |x[n]| > A \end{cases} $$

クリッピング比（Clipping Ratio）$\gamma$ をデシベルで定義すると

$$ \gamma = 10 \log_{10} \frac{A^2}{E[|x[n]|^2]} $$

クリッピングは非線形操作なので、2つの副作用があります。

1. 帯域内歪み（EVM劣化）: 信号点がずれてBERが悪化
2. 帯域外放射: クリッピングによる高調波成分がスペクトルの外に漏れる

帯域外放射はフィルタリングで抑制できますが、フィルタリング後にPAPRが再び上がる（ピーク再成長）ため、「クリッピング → フィルタリング」を繰り返す反復法が実用的です。

SLM（Selected Mapping）法

SLMは、同じデータを複数の異なる位相回転パターンで変調し、最もPAPRが低いものを選んで送信する手法です。

具体的な手順は以下のとおりです。

1. データシンボルベクトル $\bm{X} = (X[0], X[1], \dots, X[N-1])^T$ を用意
2. $U$ 個の位相回転ベクトル $\bm{B}^{(u)} = (B^{(u)}[0], B^{(u)}[1], \dots, B^{(u)}[N-1])^T$（$u = 1, \dots, U$）を生成。各 $B^{(u)}[k] = e^{j\phi_k^{(u)}}$ で、$\phi_k^{(u)}$ は $\{0, \pi/2, \pi, 3\pi/2\}$ などから選ぶ
3. 各パターンに対して $\bm{X}^{(u)} = \bm{X} \odot \bm{B}^{(u)}$（要素ごとの積）を計算
4. $\bm{x}^{(u)} = \text{IFFT}(\bm{X}^{(u)})$ を計算
5. $\text{PAPR}(\bm{x}^{(u)})$ が最小の $u^*$ を選択
6. $\bm{x}^{(u^*)}$ を送信

受信側では、どの位相回転パターンが使われたかを知る必要があります。$u^*$ のインデックス（$\lceil \log_2 U \rceil$ ビット）をサイド情報として伝送します。

SLM法の効果を定量的に評価しましょう。$U$ 個の候補のPAPRが独立で同一の分布に従うとすると、PAPRのCCDF（相補累積分布関数）は

$$ \Pr(\text{PAPR}_{\text{SLM}} > \gamma) = \left[\Pr(\text{PAPR} > \gamma)\right]^U $$

つまり、元のCCDFの $U$ 乗になります。これは大幅なPAPR低減を意味します。例えば、CCDF $= 10^{-3}$ のPAPRが、$U = 4$ で約 $1.5\,\text{dB}$、$U = 16$ で約 $3\,\text{dB}$ 低減されます。

SLMは信号の歪みを一切発生させない（歪みフリー）という利点がありますが、$U$ 回のIFFT演算が必要で計算量が増加します。

ここまでの理論をすべて統合して、OFDMシステムのフルシミュレーションをPythonで実装しましょう。

Pythonによるフルシミュレーション

以下のシミュレーションでは、パイロットシンボルの挿入、チャネル推定（LS/MMSE）、1タップ等化、BER評価、PAPR低減までを一貫して実装します。

まず、OFDMシステムの基本パラメータとユーティリティ関数を定義します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import toeplitz

# ===== システムパラメータ =====
N = 64                    # FFTサイズ（サブキャリア数）
N_cp = 16                 # サイクリックプレフィクス長
N_pilot = 8               # パイロットサブキャリア数
pilot_interval = N // N_pilot  # パイロット間隔（コムタイプ）
pilot_indices = np.arange(0, N, pilot_interval)  # パイロット位置
data_indices = np.array([i for i in range(N) if i not in pilot_indices])
N_data = len(data_indices)  # データサブキャリア数

# チャネルパラメータ
L = 10                    # チャネルタップ数
tau_rms = 3.0             # RMS遅延拡がり（サンプル単位）

# パイロットシンボル（既知、BPSK: +1）
pilot_value = np.ones(N_pilot, dtype=complex)

このコードでは、64サブキャリアのOFDMシステムを構成し、8サブキャリアごとにコムタイプのパイロットを配置しています。パイロット間隔 $64/8 = 8$ はチャネルタップ数 $L = 10$ に対して $N/L = 6.4$ を若干上回っていますが、指数減衰チャネルでは主要エネルギーが前方に集中するため、実用上は十分です。

次に、変復調関数とOFDM送受信関数を定義します。

# ===== QPSK変復調 =====
def qpsk_mod(bits):
    """ビット列をQPSKシンボルに変換（Gray符号化）"""
    symbols = (1 - 2*bits[0::2] + 1j*(1 - 2*bits[1::2])) / np.sqrt(2)
    return symbols

def qpsk_demod(symbols):
    """QPSKシンボルをビット列に復元"""
    bits = np.zeros(2 * len(symbols), dtype=int)
    bits[0::2] = (np.real(symbols) < 0).astype(int)
    bits[1::2] = (np.imag(symbols) < 0).astype(int)
    return bits

# ===== チャネル生成 =====
def generate_channel(L, tau_rms):
    """指数減衰プロファイルのレイリーフェージングチャネルを生成"""
    power_profile = np.exp(-np.arange(L) / tau_rms)
    power_profile /= np.sum(power_profile)  # 電力正規化
    h = (np.random.randn(L) + 1j * np.random.randn(L)) * np.sqrt(power_profile / 2)
    return h

# ===== OFDM送信 =====
def ofdm_tx(data_symbols, pilot_value, pilot_indices, data_indices, N, N_cp):
    """周波数領域シンボル構成 → IFFT → CP付加"""
    X = np.zeros(N, dtype=complex)
    X[pilot_indices] = pilot_value
    X[data_indices] = data_symbols
    x = np.fft.ifft(X, N)
    x_cp = np.concatenate([x[-N_cp:], x])
    return x_cp, X

# ===== OFDM受信 =====
def ofdm_rx(rx_signal, N, N_cp):
    """CP除去 → FFT"""
    y = rx_signal[N_cp:N_cp + N]
    Y = np.fft.fft(y, N)
    return Y

ofdm_tx 関数では、データシンボルとパイロットシンボルを周波数領域で所定の位置に配置し、IFFTで時間領域信号に変換してからCPを付加しています。ofdm_rx はその逆で、CP除去後にFFTを適用します。

続いて、チャネル推定関数を実装します。

# ===== チャネル推定 =====
def estimate_channel_ls(Y, pilot_indices, pilot_value, N):
    """LS推定 + 線形補間"""
    H_pilot = Y[pilot_indices] / pilot_value
    # 線形補間で全サブキャリアの推定値を得る
    H_est = np.interp(
        np.arange(N),
        pilot_indices,
        np.real(H_pilot)
    ) + 1j * np.interp(
        np.arange(N),
        pilot_indices,
        np.imag(H_pilot)
    )
    return H_est

def estimate_channel_mmse(Y, pilot_indices, pilot_value, N, snr_linear, L, tau_rms):
    """MMSE推定（パイロット位置 → 全サブキャリア）"""
    # パイロット位置でのLS推定
    H_ls_pilot = Y[pilot_indices] / pilot_value
    N_p = len(pilot_indices)

    # チャネル周波数相関行列の構成
    # 遅延プロファイルから周波数相関を計算
    power_profile = np.exp(-np.arange(L) / tau_rms)
    power_profile /= np.sum(power_profile)

    # パイロット位置間の相関行列 R_HpHp
    R_HpHp = np.zeros((N_p, N_p), dtype=complex)
    for i in range(N_p):
        for j in range(N_p):
            dk = pilot_indices[i] - pilot_indices[j]
            R_HpHp[i, j] = np.sum(
                power_profile * np.exp(-1j * 2 * np.pi * dk * np.arange(L) / N)
            )

    # 全サブキャリアとパイロット位置間の相互相関行列 R_HHp
    R_HHp = np.zeros((N, N_p), dtype=complex)
    for i in range(N):
        for j in range(N_p):
            dk = i - pilot_indices[j]
            R_HHp[i, j] = np.sum(
                power_profile * np.exp(-1j * 2 * np.pi * dk * np.arange(L) / N)
            )

    # MMSE推定
    H_est = R_HHp @ np.linalg.inv(R_HpHp + (1/snr_linear) * np.eye(N_p)) @ H_ls_pilot
    return H_est

estimate_channel_ls はLS推定値をパイロット位置で計算し、実部と虚部を別々に線形補間して全サブキャリアの推定値を得ています。estimate_channel_mmse は前節で導出したMMSE推定の式をそのまま実装しています。遅延プロファイルから周波数相関行列を構成し、ウィーナーフィルタを適用しています。

次に、PAPR低減関数を実装します。

# ===== PAPR計算 =====
def calc_papr(x):
    """PAPR [dB] を計算"""
    peak = np.max(np.abs(x)**2)
    avg = np.mean(np.abs(x)**2)
    return 10 * np.log10(peak / avg)

# ===== クリッピング =====
def clipping(x, clip_ratio_db):
    """クリッピングによるPAPR低減"""
    avg_power = np.mean(np.abs(x)**2)
    threshold = np.sqrt(avg_power * 10**(clip_ratio_db / 10))
    x_clip = x.copy()
    mask = np.abs(x_clip) > threshold
    x_clip[mask] = threshold * x_clip[mask] / np.abs(x_clip[mask])
    return x_clip

# ===== SLM法 =====
def slm(X, N, N_cp, U=8):
    """SLM法によるPAPR低減"""
    best_papr = np.inf
    best_x = None
    best_u = 0

    # 位相回転候補 {1, j, -1, -j}
    phase_set = np.array([1, 1j, -1, -1j])

    for u in range(U):
        np.random.seed(u + 1000)  # 再現性のための固定シード
        B = phase_set[np.random.randint(0, 4, N)]
        X_rotated = X * B
        x_rotated = np.fft.ifft(X_rotated, N)
        x_cp = np.concatenate([x_rotated[-N_cp:], x_rotated])
        papr = calc_papr(x_cp)
        if papr < best_papr:
            best_papr = papr
            best_x = x_cp
            best_u = u

    return best_x, best_u, best_papr

clipping 関数は、閾値を超えるサンプルの振幅を閾値に制限しつつ位相を保存します。slm 関数は $U$ 個の位相回転パターンを試し、最もPAPRが低いものを選択します。

ここまでの関数を組み合わせて、BERシミュレーションを実行します。LS推定とMMSE推定の性能を比較します。

# ===== BERシミュレーション =====
np.random.seed(42)
SNR_dB_list = np.arange(0, 31, 3)
N_ofdm_symbols = 500  # OFDMシンボル数

ber_ls = []
ber_mmse = []
ber_perfect = []

for snr_db in SNR_dB_list:
    snr_linear = 10**(snr_db / 10)
    errors_ls = 0
    errors_mmse = 0
    errors_perfect = 0
    total_bits = 0

    for _ in range(N_ofdm_symbols):
        # チャネル生成（シンボルごとに変化）
        h = generate_channel(L, tau_rms)
        H_true = np.fft.fft(h, N)

        # 送信データ生成
        tx_bits = np.random.randint(0, 2, N_data * 2)
        tx_symbols = qpsk_mod(tx_bits)

        # OFDM送信
        tx_signal, X_freq = ofdm_tx(
            tx_symbols, pilot_value, pilot_indices, data_indices, N, N_cp
        )

        # チャネル通過
        rx_channel = np.convolve(tx_signal, h)[:len(tx_signal)]

        # AWGN付加
        signal_power = np.mean(np.abs(rx_channel)**2)
        noise_power = signal_power / snr_linear
        noise = np.sqrt(noise_power / 2) * (
            np.random.randn(len(rx_channel)) + 1j * np.random.randn(len(rx_channel))
        )
        rx_signal = rx_channel + noise

        # OFDM受信
        Y = ofdm_rx(rx_signal, N, N_cp)

        # チャネル推定
        H_est_ls = estimate_channel_ls(Y, pilot_indices, pilot_value, N)
        H_est_mmse = estimate_channel_mmse(
            Y, pilot_indices, pilot_value, N, snr_linear, L, tau_rms
        )

        # ZF等化 + 復調
        for H_est, errors_ref in [
            (H_est_ls, 'ls'), (H_est_mmse, 'mmse'), (H_true, 'perfect')
        ]:
            X_eq = Y[data_indices] / H_est[data_indices]
            rx_bits = qpsk_demod(X_eq)
            err = np.sum(tx_bits != rx_bits)
            if errors_ref == 'ls':
                errors_ls += err
            elif errors_ref == 'mmse':
                errors_mmse += err
            else:
                errors_perfect += err

        total_bits += len(tx_bits)

    ber_ls.append(errors_ls / total_bits)
    ber_mmse.append(errors_mmse / total_bits)
    ber_perfect.append(errors_perfect / total_bits)
    print(f"SNR={snr_db:2d}dB  BER(LS)={ber_ls[-1]:.5f}  "
          f"BER(MMSE)={ber_mmse[-1]:.5f}  BER(Perfect)={ber_perfect[-1]:.5f}")

このシミュレーションでは、OFDMシンボルごとに新しいレイリーフェージングチャネルを生成し、LS推定・MMSE推定・完全チャネル情報の3条件でBERを比較しています。完全チャネル情報（Perfect CSI）は推定誤差がゼロの理想条件であり、上界として機能します。

結果をプロットしましょう。

fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# (1) BER比較
ax1 = axes[0]
ax1.semilogy(SNR_dB_list, ber_ls, 'bs-', label='LS estimation', markersize=6)
ax1.semilogy(SNR_dB_list, ber_mmse, 'ro-', label='MMSE estimation', markersize=6)
ax1.semilogy(SNR_dB_list, ber_perfect, 'g^--', label='Perfect CSI', markersize=6)
ax1.set_xlabel('SNR [dB]')
ax1.set_ylabel('BER')
ax1.set_title('BER: LS vs MMSE vs Perfect CSI')
ax1.legend()
ax1.grid(True, which='both')
ax1.set_ylim([1e-5, 1])

# (2) チャネル推定MSE
ax2 = axes[1]
mse_ls_list = []
mse_mmse_list = []
for snr_db in SNR_dB_list:
    snr_linear = 10**(snr_db / 10)
    mse_ls_sum = 0
    mse_mmse_sum = 0
    N_trial = 200
    for _ in range(N_trial):
        h = generate_channel(L, tau_rms)
        H_true = np.fft.fft(h, N)
        tx_bits = np.random.randint(0, 2, N_data * 2)
        tx_symbols = qpsk_mod(tx_bits)
        tx_signal, X_freq = ofdm_tx(
            tx_symbols, pilot_value, pilot_indices, data_indices, N, N_cp
        )
        rx_channel = np.convolve(tx_signal, h)[:len(tx_signal)]
        signal_power = np.mean(np.abs(rx_channel)**2)
        noise_power = signal_power / snr_linear
        noise = np.sqrt(noise_power / 2) * (
            np.random.randn(len(rx_channel)) + 1j * np.random.randn(len(rx_channel))
        )
        rx_signal = rx_channel + noise
        Y = ofdm_rx(rx_signal, N, N_cp)
        H_ls = estimate_channel_ls(Y, pilot_indices, pilot_value, N)
        H_mmse = estimate_channel_mmse(
            Y, pilot_indices, pilot_value, N, snr_linear, L, tau_rms
        )
        mse_ls_sum += np.mean(np.abs(H_true - H_ls)**2)
        mse_mmse_sum += np.mean(np.abs(H_true - H_mmse)**2)
    mse_ls_list.append(mse_ls_sum / N_trial)
    mse_mmse_list.append(mse_mmse_sum / N_trial)

ax2.semilogy(SNR_dB_list, mse_ls_list, 'bs-', label='LS + interpolation')
ax2.semilogy(SNR_dB_list, mse_mmse_list, 'ro-', label='MMSE')
ax2.set_xlabel('SNR [dB]')
ax2.set_ylabel('MSE of channel estimation')
ax2.set_title('Channel Estimation MSE')
ax2.legend()
ax2.grid(True, which='both')

# (3) コンスタレーション比較（SNR=20dB）
ax3 = axes[2]
h_demo = generate_channel(L, tau_rms)
H_demo = np.fft.fft(h_demo, N)
tx_bits_demo = np.random.randint(0, 2, N_data * 2)
tx_sym_demo = qpsk_mod(tx_bits_demo)
tx_sig_demo, X_demo = ofdm_tx(
    tx_sym_demo, pilot_value, pilot_indices, data_indices, N, N_cp
)
rx_ch_demo = np.convolve(tx_sig_demo, h_demo)[:len(tx_sig_demo)]
snr_demo = 10**(20/10)
sig_pow = np.mean(np.abs(rx_ch_demo)**2)
n_pow = sig_pow / snr_demo
noise_demo = np.sqrt(n_pow/2) * (
    np.random.randn(len(rx_ch_demo)) + 1j * np.random.randn(len(rx_ch_demo))
)
Y_demo = ofdm_rx(rx_ch_demo + noise_demo, N, N_cp)
H_ls_demo = estimate_channel_ls(Y_demo, pilot_indices, pilot_value, N)
X_eq_demo = Y_demo[data_indices] / H_ls_demo[data_indices]
ax3.scatter(np.real(tx_sym_demo), np.imag(tx_sym_demo),
            s=30, marker='x', c='green', label='TX (ideal)', zorder=3)
ax3.scatter(np.real(X_eq_demo), np.imag(X_eq_demo),
            s=10, alpha=0.6, c='blue', label='RX (after EQ)')
ax3.set_xlabel('In-phase')
ax3.set_ylabel('Quadrature')
ax3.set_title('Constellation at SNR=20dB (LS+ZF)')
ax3.legend()
ax3.grid(True)
ax3.set_aspect('equal')

plt.tight_layout()
plt.savefig('ofdm_detail_ber.png', dpi=150, bbox_inches='tight')
plt.show()

上のグラフから、3つの重要な結果が読み取れます。

1. BER曲線（左図）: MMSE推定はLS推定に比べて約2〜4 dBのSNRゲインを達成しています。特にSNR 10〜20 dBの実用的な範囲でその差が顕著です。これは、MMSE推定がチャネルの周波数相関を利用して雑音を効果的に抑制しているためです。Perfect CSIとMMSE推定の差は1〜2 dB程度で、MMSE推定が理想に近い性能を持つことがわかります。
2. チャネル推定MSE（中央図）: LS推定のMSEは $1/\text{SNR}$ に比例して減少しますが、MMSE推定はチャネルの相関構造を利用することでMSEがより速く減少します。高SNR領域では両者の差が縮まりますが、低SNR領域ではMMSEの優位性が明確です。
3. コンスタレーション（右図）: SNR 20 dBでのLS推定 + ZF等化後のコンスタレーションは、理想的なQPSK信号点（緑色の×）の周りに青い点が集まっています。点の広がりが推定誤差と雑音の影響を反映しています。

PAPR低減の効果を可視化する

クリッピングとSLM法の効果を、PAPRのCCDFで比較します。

# ===== PAPR低減シミュレーション =====
N_trial_papr = 5000
papr_original = []
papr_clipped_6dB = []
papr_clipped_4dB = []
papr_slm_4 = []
papr_slm_16 = []

for _ in range(N_trial_papr):
    bits = np.random.randint(0, 2, N_data * 2)
    symbols = qpsk_mod(bits)

    # 通常のOFDM信号
    tx_signal, X_freq = ofdm_tx(
        symbols, pilot_value, pilot_indices, data_indices, N, N_cp
    )
    papr_original.append(calc_papr(tx_signal))

    # クリッピング（CR=6dB）
    x_clip_6 = clipping(tx_signal, 6.0)
    papr_clipped_6dB.append(calc_papr(x_clip_6))

    # クリッピング（CR=4dB）
    x_clip_4 = clipping(tx_signal, 4.0)
    papr_clipped_4dB.append(calc_papr(x_clip_4))

    # SLM（U=4）
    _, _, papr_slm = slm(X_freq, N, N_cp, U=4)
    papr_slm_4.append(papr_slm)

    # SLM（U=16）
    _, _, papr_slm = slm(X_freq, N, N_cp, U=16)
    papr_slm_16.append(papr_slm)

# CCDFプロット
def plot_ccdf(values, label, ax, linestyle='-'):
    sorted_vals = np.sort(values)
    ccdf = 1 - np.arange(1, len(sorted_vals) + 1) / len(sorted_vals)
    ax.semilogy(sorted_vals, ccdf, linestyle, label=label, linewidth=1.5)

fig, ax = plt.subplots(figsize=(8, 5))
plot_ccdf(papr_original, 'Original OFDM', ax)
plot_ccdf(papr_clipped_6dB, 'Clipping (CR=6dB)', ax, '--')
plot_ccdf(papr_clipped_4dB, 'Clipping (CR=4dB)', ax, '-.')
plot_ccdf(papr_slm_4, 'SLM (U=4)', ax, '--')
plot_ccdf(papr_slm_16, 'SLM (U=16)', ax, '-.')
ax.set_xlabel('PAPR [dB]')
ax.set_ylabel('CCDF = Pr(PAPR > x)')
ax.set_title('CCDF of PAPR: Clipping vs SLM')
ax.legend()
ax.grid(True, which='both')
ax.set_xlim([2, 12])
ax.set_ylim([1e-4, 1])
plt.tight_layout()
plt.savefig('ofdm_papr_reduction.png', dpi=150, bbox_inches='tight')
plt.show()

このCCDFプロットから、以下の点が読み取れます。

1. Original OFDM: $N = 64$ のQPSK-OFDMでは、CCDF $= 10^{-3}$ でのPAPRが約9〜10 dBに達します。つまり、1000シンボルに1回はこの値を超えるピークが発生します。
2. クリッピング: CR = 6 dBでは約1.5 dBのPAPR低減、CR = 4 dBでは約3 dBの低減が得られます。ただし、クリッピング比を下げすぎるとBERが悪化するため、通常は5〜6 dB程度に設定します。
3. SLM法: $U = 4$ で約1.5 dBの低減、$U = 16$ で約2.5〜3 dBの低減が得られます。SLMは歪みフリーの手法であり、BER劣化を起こさない点でクリッピングより優れています。ただし、計算量が $U$ 倍になる点と、サイド情報の伝送が必要な点がトレードオフです。

チャネル推定精度の可視化

最後に、1つのOFDMシンボルについて、真のチャネル応答とLS/MMSE推定値を周波数領域で重ねてプロットし、推定の様子を視覚的に確認します。

# ===== チャネル推定の可視化 =====
np.random.seed(123)
h_vis = generate_channel(L, tau_rms)
H_vis = np.fft.fft(h_vis, N)

# 送信・受信
tx_bits_vis = np.random.randint(0, 2, N_data * 2)
tx_sym_vis = qpsk_mod(tx_bits_vis)
tx_sig_vis, X_vis = ofdm_tx(
    tx_sym_vis, pilot_value, pilot_indices, data_indices, N, N_cp
)
rx_ch_vis = np.convolve(tx_sig_vis, h_vis)[:len(tx_sig_vis)]

# SNR = 15 dBでの推定
snr_vis = 10**(15/10)
sig_pow_vis = np.mean(np.abs(rx_ch_vis)**2)
n_pow_vis = sig_pow_vis / snr_vis
noise_vis = np.sqrt(n_pow_vis / 2) * (
    np.random.randn(len(rx_ch_vis)) + 1j * np.random.randn(len(rx_ch_vis))
)
Y_vis = ofdm_rx(rx_ch_vis + noise_vis, N, N_cp)

H_ls_vis = estimate_channel_ls(Y_vis, pilot_indices, pilot_value, N)
H_mmse_vis = estimate_channel_mmse(
    Y_vis, pilot_indices, pilot_value, N, snr_vis, L, tau_rms
)

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

# 振幅
ax1 = axes[0]
ax1.plot(np.arange(N), 20*np.log10(np.abs(H_vis)), 'k-', linewidth=2, label='True H[k]')
ax1.plot(np.arange(N), 20*np.log10(np.abs(H_ls_vis) + 1e-10),
         'b--', alpha=0.7, label='LS estimate')
ax1.plot(np.arange(N), 20*np.log10(np.abs(H_mmse_vis) + 1e-10),
         'r-.', alpha=0.7, label='MMSE estimate')
ax1.scatter(pilot_indices, 20*np.log10(np.abs(H_vis[pilot_indices])),
            c='green', s=50, zorder=5, label='Pilot positions')
ax1.set_xlabel('Subcarrier index k')
ax1.set_ylabel('|H[k]| [dB]')
ax1.set_title('Channel Frequency Response (Magnitude) at SNR=15dB')
ax1.legend()
ax1.grid(True)

# 位相
ax2 = axes[1]
ax2.plot(np.arange(N), np.angle(H_vis), 'k-', linewidth=2, label='True H[k]')
ax2.plot(np.arange(N), np.angle(H_ls_vis), 'b--', alpha=0.7, label='LS estimate')
ax2.plot(np.arange(N), np.angle(H_mmse_vis), 'r-.', alpha=0.7, label='MMSE estimate')
ax2.set_xlabel('Subcarrier index k')
ax2.set_ylabel('Phase [rad]')
ax2.set_title('Channel Frequency Response (Phase) at SNR=15dB')
ax2.legend()
ax2.grid(True)

plt.tight_layout()
plt.savefig('ofdm_channel_estimation.png', dpi=150, bbox_inches='tight')
plt.show()

このグラフから、LS推定とMMSE推定の違いが視覚的にわかります。

1. 振幅（左図）: 黒い実線が真のチャネル応答 $|H[k]|$ です。LS推定（青い破線）は特にパイロット間のサブキャリアで揺らぎ（ギザギザ）が見られます。これはLS推定が雑音の影響をそのまま受け、線形補間では雑音を平滑化しきれないためです。一方、MMSE推定（赤い一点鎖線）は真の値に非常に近い滑らかな曲線を描いています。これはチャネルの周波数相関を利用した雑音抑制の効果です。
2. 位相（右図）: 位相の推定でも同様の傾向が見られます。LS推定は雑音による位相のばらつきが大きく、MMSE推定はより安定した推定を行っています。位相誤差はコンスタレーションの回転として現れるため、位相推定の精度はBERに直結します。
3. パイロット位置: 緑色の点で示したパイロット位置では、LS推定値とMMSE推定値がともに真の値に近い値を示しています。パイロット間のサブキャリアでの補間精度の差が、2つの手法のBER差を生んでいることが確認できます。

まとめ

本記事では、OFDMの実装に不可欠な詳細技術について解説しました。

• サイクリックプレフィクスの数学的役割: CPがISIを除去し、かつ線形畳み込みを巡回畳み込みに変換することを厳密に証明しました。これにより $Y[k] = H[k] X[k] + W[k]$ という周波数領域の入出力関係が成立し、1タップ等化が可能になります
• CP長の設計: $N_{\text{CP}} \geq L – 1$ が基本条件で、帯域効率 $\eta = N/(N + N_{\text{CP}})$ とのトレードオフがあります。5G NR ではニューメロロジーによって用途に応じたCP設計を行います
• パイロット配置: コムタイプ・ブロックタイプ・格子タイプの3種類があり、チャネルの時間変動と周波数選択性に応じて選択します。パイロット間隔はナイキスト条件で決まります
• チャネル推定: LS推定は $\hat{H}[k] = Y_p[k] / X_p[k]$ で計算量は最小ですが雑音に弱い。MMSE推定はチャネルの相関構造を利用して2〜4 dBのSNRゲインを達成します
• 1タップ等化: ZF等化 $\hat{X}[k] = Y[k] / H[k]$ は単純ですが深いフェード点で雑音増幅。MMSE等化は $\hat{X}[k] = H^*[k] Y[k] / (|H[k]|^2 + 1/\text{SNR})$ でバランスをとります
• PAPR低減: クリッピングは実装が簡単ですが信号歪みが発生。SLM法は歪みフリーでCCDF $10^{-3}$ のPAPRを2〜3 dB低減できます

OFDMの基本原理から実装詳細まで、2記事にわたって体系的に学びました。これらの知識は、5G NR や Wi-Fi 6 の物理層設計を理解するうえで不可欠です。

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

• 多重化技術（FDM・TDM・CDM）を体系的に理解する
• 等化器（イコライザ）の設計