スマートフォンでWi-Fiにつながった瞬間や、人工衛星からの電波を地上局が受け取った瞬間、受信機はある厄介な現象と必ず戦っています。送信機と受信機はそれぞれ別の水晶発振器を持っており、その周波数はわずかにずれています。さらに送信機と受信機が相対的に動いていれば、ドップラー効果で受信周波数がずれます。たった数十ppm（百万分率）のずれでも、現代の高速通信、特にOFDM（直交周波数分割多重）にとっては致命的になりかねません。この「搬送波周波数のずれ」を搬送波周波数オフセット（CFO: Carrier Frequency Offset）と呼びます。
なぜCFOがそれほど深刻なのでしょうか。OFDMは多数の搬送波（サブキャリア）を直交させて密に並べることで高い周波数利用効率を得ています。この「直交」が成り立つのは、各サブキャリアの周波数がちょうど整数間隔で揃っているときだけです。受信機側で周波数が少しずれると、サブキャリア同士が互いに干渉し合い（ICI: Inter-Carrier Interference）、せっかく分離していた信号がぐちゃぐちゃに混ざってしまいます。だからこそ、復調の前に受信機はCFOを正確に推定し、補正しなければなりません。
CFO推定を理解すると、次のような分野への応用が見えてきます。
- 無線LAN（Wi-Fi）/ LTE / 5G: パケットの先頭に置かれた既知の繰り返し信号（プリアンブル）からCFOを推定し、ペイロード復調の前に補正します。Schmidl&Cox法やMoose法はこの目的で広く使われています。
- 衛星通信・深宇宙通信: 高速移動する衛星では大きなドップラーシフトが生じ、受信周波数が大きく動きます。受信機はまず周波数を引き込み、その後トラッキングします。
- GNSS（測位衛星）受信機: 衛星と受信機の相対速度によるドップラーを推定することは、信号捕捉と速度測定の両方の基礎になります。
上の図は、CFOがどこから来てどう信号に現れるかを一枚にまとめたものです。送信機の搬送波 $f_c$ と受信機発振器 $f_c’$ の差、そして相対運動によるドップラー $\Delta f_d$ が足し合わさって残留周波数差 $\Delta f$ になり、最終的に受信信号へ $e^{j2\pi\varepsilon n/N}$ という回転因子として掛かります。本記事はこの $\varepsilon$ をプリアンブルから測り取るまでが主題です。
本記事の内容
- CFOがOFDMにもたらす害（位相回転とICI）の直感的理解
- 受信信号モデル $r[n] = s[n] e^{j 2\pi \varepsilon n / N}$ から、繰り返し区間の相関位相を使った推定量を省略なく導出する（Moose法 / Schmidl&Cox法）
- 推定可能範囲（曖昧さ $\pm 1/2D$）の意味
- 推定分散のCramér-Rao下界（CRLB）の導出と、推定量がどこまで近づくか
- numpyでOFDM信号に既知CFOとAWGNを与え、RMSEをSNRに対してプロットして理論を実証する
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- OFDM（直交周波数分割多重）の原理 — サブキャリアの直交性とIFFT/FFTによる変復調。CFOがなぜ直交性を壊すかの土台になります
- IQ直交変調の仕組み — 複素ベースバンド信号の表現。CFOを $e^{j2\pi\varepsilon n/N}$ という複素指数で書く根拠になります
- AWGNチャネルモデル — 加法的白色ガウス雑音の扱い。推定分散やCRLBの計算で使います
搬送波周波数オフセットとは
ラジオのチューニングを思い浮かべてください。目的の放送局にダイヤルがぴったり合っていないと、音がこもったり混信したりします。受信機の「内部の物差し」である局部発振器の周波数が、送信側とぴったり合っていないと、信号をきれいに取り出せないのです。CFOとは、まさにこの「物差しのずれ」です。
もう少し物理的に言うと、送信機は搬送波周波数 $f_c$ で電波を出します。受信機はその電波をベースバンド（低い周波数帯）に落とすために、自分の発振器が作る $f_c’$ という周波数を掛け合わせます。理想的には $f_c = f_c’$ ですが、現実の水晶発振器には製造ばらつきや温度ドリフトがあり、数ppmのずれが残ります。さらに送受信機が相対速度 $v$ で動いていれば、ドップラー効果で $\Delta f_d = f_c \cdot v/c$ だけ周波数がずれます。これらを合わせた残留周波数差 $\Delta f = f_c – f_c’ + \Delta f_d$ がCFOの正体です。
このずれが受信ベースバンド信号にどう現れるかを見てみましょう。送信したベースバンド信号を $s(t)$ とすると、CFOがある場合、受信ベースバンド信号は
$$ r(t) = s(t) \, e^{j 2\pi \Delta f \, t} $$
となります。つまり、CFOは時間とともに回転する複素指数として信号に掛かります。時間 $t$ が進むほど位相 $2\pi \Delta f \, t$ がどんどん回っていくのです。サンプリング周期を $T_s$、サンプル番号を $n$（時刻 $t = n T_s$）とすると、離散信号では
$$ r[n] = s[n] \, e^{j 2\pi \Delta f \, n T_s} $$
となります。ここでOFDMの文脈では、サブキャリア間隔 $\Delta f_{\text{sc}} = 1/(N T_s)$（$N$ はFFTサイズ）で規格化した正規化周波数オフセット
$$ \varepsilon = \frac{\Delta f}{\Delta f_{\text{sc}}} = \Delta f \cdot N T_s $$
を使うのが便利です。$\varepsilon$ は「サブキャリア間隔何個分のずれか」を表す無次元量です。これを使うと $\Delta f \, n T_s = \varepsilon n / N$ となり、受信信号モデルは本記事の主役である次の形に整理されます。
$$ \begin{equation} r[n] = s[n] \, e^{j 2\pi \varepsilon n / N} \end{equation} $$
この $\varepsilon$ を、整数部分 $\varepsilon_I$ と小数部分 $\varepsilon_F$（$|\varepsilon_F| \le 1/2$）に分けて $\varepsilon = \varepsilon_I + \varepsilon_F$ と書くことがあります。整数部のずれはサブキャリアの「番号がまるごとずれる」効果を持ち、小数部のずれは直交性を壊してICIを生みます。本記事で扱うプリアンブル相関法は主にこの小数部 $\varepsilon_F$ を推定する手法です。
ここまでで、CFOが「回転する複素指数」として信号に掛かることがわかりました。次に、この回転がOFDMの直交性をどう破壊するのかを、もう少し定量的に見ていきましょう。
CFOがOFDMに与える害 — 位相回転とICI
CFOがない理想的なOFDMでは、$N$ 個のサブキャリアが直交しているため、FFTを取れば各サブキャリアのシンボルがきれいに分離されます。サブキャリア $k$ のシンボルを $X[k]$ とすると、送信時間信号は $s[n] = \frac{1}{N}\sum_{k} X[k] e^{j2\pi kn/N}$ です。受信側でCFO $\varepsilon$ が掛かると、受信信号は $r[n] = s[n] e^{j2\pi\varepsilon n/N}$ になります。
このときFFT出力（サブキャリア $m$ の成分）$Y[m] = \sum_{n=0}^{N-1} r[n] e^{-j2\pi mn/N}$ を計算すると、
$$ Y[m] = \frac{1}{N}\sum_{n=0}^{N-1} \sum_{k} X[k] \, e^{j2\pi kn/N} \, e^{j2\pi\varepsilon n/N} \, e^{-j2\pi mn/N} $$
となります。$n$ についての和を先に取るために、指数をまとめます。指数の肩は $\frac{2\pi n}{N}(k + \varepsilon – m)$ なので、
$$ Y[m] = \frac{1}{N}\sum_{k} X[k] \sum_{n=0}^{N-1} e^{j 2\pi n (k – m + \varepsilon)/N} $$
と書けます。ここで内側の幾何級数の和は、公比 $q = e^{j2\pi(k-m+\varepsilon)/N}$ の等比数列の和です。$\sum_{n=0}^{N-1} q^n = \frac{1-q^N}{1-q}$ を使うと、これは $k=m$ のとき（$\varepsilon=0$ なら）$N$ になり、それ以外では0になる——これが直交性です。ところが $\varepsilon \neq 0$ だと、$k=m$ の項ですら肩に $\varepsilon$ が残り、和は $N$ になりません。整理すると
$$ Y[m] = \underbrace{X[m] \, e^{j\pi\varepsilon(N-1)/N} \frac{\sin(\pi\varepsilon)}{N\sin(\pi\varepsilon/N)}}_{\text{自分の信号(振幅減衰+位相回転)}} + \underbrace{\sum_{k \neq m} X[k] (\cdots)}_{\text{ICI(他サブキャリアの漏れ込み)}} $$
という形になります。第1項は、本来の信号 $X[m]$ が $\frac{\sin(\pi\varepsilon)}{N\sin(\pi\varepsilon/N)}$ という係数で振幅が減衰し、さらに位相が回転していることを示します。第2項は、本来分離されていたはずの他サブキャリア $X[k]$ が漏れ込むICIです。$\varepsilon$ が小さいほど害は小さく、$\varepsilon=0$ なら完全に消えます。
要するに、CFOは（1）狙ったサブキャリアの振幅と位相を狂わせ、（2）隣のサブキャリアを混ぜ込んでSNRを劣化させます。たとえば $\varepsilon = 0.1$ 程度でもICIによるSNR劣化は無視できず、復調誤りが急増します。だからこそ、復調の前にこの $\varepsilon$ を推定して取り除く必要があるのです。
これを実際にOFDM信号で測ってみたのが次の図です。QPSKを載せたOFDMシンボルにCFOを与え、受信側でFFT復調した星座点を $\varepsilon=0$ と $\varepsilon=0.1$ で比較しました。
左の $\varepsilon=0$ では4つの判定点（赤い星）の上に受信点がきれいに集まっていますが、右の $\varepsilon=0.1$ では受信点が回転・収縮し、ICIで雲のように散らばっています。誤差ベクトルの大きさ（EVM）は0.1程度の小さなCFOだけで大きく悪化しており、わずかな $\varepsilon$ でも復調が崩れることが実測で確認できます。
さらに $\varepsilon$ を変えながら、狙ったサブキャリアの振幅減衰係数と信号対ICI比（SIR）を測ったのが次の図です。
青線は導出した振幅係数 $|\sin(\pi\varepsilon)/(N\sin(\pi\varepsilon/N))|$ で、$\varepsilon$ が大きいほど自分の信号が痩せていきます。オレンジの実測SIRは $\varepsilon=0.1$ で約9dB、$\varepsilon=0.3$ では約0dB（信号とICIがほぼ同電力）まで落ち込みます。$\varepsilon$ がわずか0.3でも信号がICIに埋もれることがわかり、CFO補正が必須である理由が定量的に見て取れます。
では、どうやって $\varepsilon$ を測ればよいのでしょうか。鍵になるのは「同じ信号を2回送る」という単純なアイデアです。次節でその原理を見ていきます。
繰り返しプリアンブルによる推定の原理
CFO推定の最も基本的で強力なアイデアは、既知の信号を一定間隔で2回（あるいは複数回）繰り返して送ることです。Wi-FiのショートプリアンブルやSchmidl&Cox法がこの構造を使っています。
直感を説明しましょう。いま、ある信号サンプル $s[n]$ を送り、その $D$ サンプル後にまったく同じ信号 $s[n+D] = s[n]$ を送ったとします（繰り返し）。CFOがなければ、受信側でもこの2つは完全に同一です。ところがCFOがあると、受信信号は $r[n] = s[n]e^{j2\pi\varepsilon n/N}$、$r[n+D] = s[n+D]e^{j2\pi\varepsilon (n+D)/N} = s[n]e^{j2\pi\varepsilon (n+D)/N}$ となります。両者の比を取ると、$s[n]$ がキャンセルして
$$ \frac{r[n+D]}{r[n]} = \frac{s[n]e^{j2\pi\varepsilon(n+D)/N}}{s[n]e^{j2\pi\varepsilon n/N}} = e^{j2\pi\varepsilon D/N} $$
となります。つまり、繰り返した2つのサンプルの間には、$D$ サンプル分だけ進んだ一定の位相差 $2\pi\varepsilon D/N$ が現れます。この位相差を測れば、そこから $\varepsilon$ を逆算できる——これがCFO推定の核心です。雑音のない理想状態なら、たった1組のサンプルの位相差から $\varepsilon$ がわかってしまいます。
ただし実際には雑音があるので、1点だけでなく繰り返し区間全体にわたって相関を取り、平均化して雑音を抑えます。具体的には、繰り返し区間の各サンプルについて $r[n+D] r^*[n]$ を計算し、これを足し合わせた相関量
$$ P = \sum_{n} r[n+D] \, r^*[n] $$
の位相を見ます（$r^*$ は複素共役）。各項 $r[n+D]r^*[n] = |s[n]|^2 e^{j2\pi\varepsilon D/N}$ は同じ位相 $2\pi\varepsilon D/N$ を持つので、足し合わせても位相は変わらず、振幅 $|s[n]|^2$ が同位相で積み上がります。雑音は位相がばらばらなので平均化で打ち消され、信号成分だけが強め合います。これがコヒーレント積分の威力です。
このアイデアを数式として厳密に定式化し、推定量を導出するのが次節です。
CFO推定量の導出（Moose法 / Schmidl&Cox法）
推定量の導出
ゴールを明確にしましょう。繰り返し区間の相関 $P$ の位相から $\varepsilon$ を推定する式を、雑音モデルのもとで導くことが目標です。
受信信号モデルは、雑音を加えて
$$ r[n] = s[n] \, e^{j 2\pi \varepsilon n / N} + w[n] $$
とします。ここで $w[n] \sim \mathcal{CN}(0, \sigma_w^2)$ は循環対称複素ガウス雑音です。プリアンブルは周期 $D$ で繰り返されているとし、$s[n+D] = s[n]$ が成り立つサンプル範囲 $n \in \{0, 1, \dots, L-1\}$（$L$ は相関に使うサンプル数）を考えます。
相関量を定義します。
$$ P = \sum_{n=0}^{L-1} r[n+D] \, r^*[n] $$
まず雑音のない理想ケースで中身を確認します。$r[n] = s[n]e^{j2\pi\varepsilon n/N}$、$r[n+D] = s[n]e^{j2\pi\varepsilon(n+D)/N}$ なので、共役を取って積を作ると
$$ r[n+D]\,r^*[n] = s[n]\,e^{j2\pi\varepsilon(n+D)/N} \cdot s^*[n]\,e^{-j2\pi\varepsilon n/N} $$
となります。$s[n]s^*[n] = |s[n]|^2$ にまとめ、指数の肩を計算します。肩は $\frac{2\pi\varepsilon}{N}\big((n+D) – n\big) = \frac{2\pi\varepsilon D}{N}$ で $n$ に依存しません。したがって
$$ r[n+D]\,r^*[n] = |s[n]|^2 \, e^{j 2\pi\varepsilon D / N} $$
となります。これを $n=0$ から $L-1$ まで足すと、共通の位相因子 $e^{j2\pi\varepsilon D/N}$ がくくり出せて
$$ P = \left(\sum_{n=0}^{L-1} |s[n]|^2\right) e^{j 2\pi \varepsilon D / N} $$
を得ます。$\sum |s[n]|^2$ は正の実数なので、$P$ の偏角（位相）はちょうど $2\pi\varepsilon D/N$ になります。式で書くと
$$ \arg\{P\} = \frac{2\pi \varepsilon D}{N} $$
です。両辺を $\frac{2\pi D}{N}$ で割って $\varepsilon$ について解くと、推定量
$$ \begin{equation} \hat{\varepsilon} = \frac{N}{2\pi D} \arg\left\{ \sum_{n=0}^{L-1} r[n+D] \, r^*[n] \right\} \end{equation} $$
が得られます。これがMoose法（A. P. Mooseが1994年に提案）の推定量です。OFDMの文脈ではFFTサイズ $N$ で正規化していますが、時間領域だけで考えるなら $\hat{\varepsilon}_{\text{Hz単位}} = \frac{1}{2\pi D T_s}\arg\{P\}$ と同じことです。
雑音 $w[n]$ がある場合は、$P$ に雑音項が加わって偏角がわずかに揺らぎます。$L$ 個のサンプルにわたる和（コヒーレント積分）が、この揺らぎを統計的に平均化して小さく抑えます。雑音がどの程度の推定誤差を生むかは、後の節でCRLBとして定量化します。
Schmidl&Cox法との関係
Schmidl&Cox法（1997年）は、Moose法の考え方をOFDMのプリアンブル設計と組み合わせて実用化した古典的アルゴリズムです。先頭OFDMシンボルを、前半 $N/2$ サンプルと後半 $N/2$ サンプルが同一になるように設計します（周波数領域で偶数サブキャリアだけに変調を載せ、奇数サブキャリアを0にすると、時間領域で前半と後半が一致するという性質を使います）。すると $D = N/2$ の繰り返し構造ができ、Moose法の推定量がそのまま使えます。
Schmidl&Coxではタイミング同期も同時に行います。スライディング窓で
$$ M(d) = \frac{|P(d)|^2}{\big(R(d)\big)^2}, \quad P(d) = \sum_{n=0}^{N/2-1} r^*[d+n]\,r[d+n+N/2], \quad R(d) = \sum_{n=0}^{N/2-1} |r[d+n+N/2]|^2 $$
というメトリックを計算し、$M(d)$ がプラトー（平坦域）になる位置でシンボル先頭を見つけます。タイミングが定まったら、その位置の $P(d)$ の偏角からCFOの小数部を
$$ \hat{\varepsilon}_F = \frac{1}{\pi}\arg\{P(d)\} $$
で推定します（$D=N/2$ を式(2)に代入すると $\frac{N}{2\pi \cdot N/2} = \frac{1}{\pi}$ になることに対応）。整数部 $\varepsilon_I$ は2番目のプリアンブル（既知パターン）を使って別途推定します。本記事では小数部推定の本質であるMoose法に焦点を当てます。
ここまでで推定式が手に入りました。しかし、この式には重要な限界があります。$\arg$（偏角）は $2\pi$ の不確定性を持つので、推定できる $\varepsilon$ の範囲には限りがあるのです。次節でこの「曖昧さ」を詳しく見ます。
推定可能範囲と曖昧さ（$\pm 1/2D$）
推定式 $\hat{\varepsilon} = \frac{N}{2\pi D}\arg\{P\}$ の弱点は、$\arg$ 関数の値域が $(-\pi, \pi]$ に限られることです。偏角は $2\pi$ だけずれても区別がつきません。$\arg\{P\}$ が $(-\pi, \pi]$ の範囲しか返せないということは、真の位相 $\frac{2\pi\varepsilon D}{N}$ がこの範囲に収まっていなければ正しく逆算できない、ということを意味します。
条件を式にすると、
$$ -\pi < \frac{2\pi\varepsilon D}{N} \le \pi $$
です。各辺を $\frac{2\pi D}{N}$ で割ると、推定可能な $\varepsilon$ の範囲は
$$ -\frac{N}{2D} < \varepsilon \le \frac{N}{2D} $$
となります。OFDMでサブキャリア間隔単位の $\varepsilon$ で考えると、$D=N/2$（Schmidl&Cox）のとき範囲は $-1 < \varepsilon \le 1$、すなわち $\pm 1$ サブキャリア間隔まで推定できます。一方、繰り返し間隔 $D$ を時間領域サンプル単位で見て、$\varepsilon$ を「正規化周波数 $f_{\text{norm}} = \Delta f \cdot T_s$」で測ると、推定範囲は $\pm \frac{1}{2D}$（cycle/sample）になります。これが見出しの「曖昧さ $\pm 1/2D$」です。
ここに本質的なトレードオフがあります。
- $D$ を大きくすると: 位相差 $2\pi\varepsilon D/N$ が大きくなるので、同じ $\varepsilon$ でも位相が大きく回り、雑音に対して測りやすくなります（推定の感度が上がり分散が下がる）。しかし推定範囲 $\pm N/2D$ は狭くなります。
- $D$ を小さくすると: 推定範囲は広がりますが、位相差が小さくなるので雑音に弱くなります（分散が増える）。
つまり「広い範囲を粗く測る（小さい $D$）」か「狭い範囲を精密に測る（大きい $D$）」かの選択です。実用システムでは、複数の異なる $D$ を組み合わせて、まず粗く引き込んでから精密化する多段推定を行うこともあります。Wi-Fiのショートプリアンブル（短い周期、広範囲・粗い）とロングプリアンブル（長い周期、狭範囲・精密）の組み合わせはまさにこの思想です。
範囲の限界を理解したところで、次は「雑音がある中で、どこまで正確に推定できるのか」という性能限界を考えます。推定誤差の理論的下限であるCramér-Rao下界を導きます。
推定分散の理論限界 — Cramér-Rao下界
Cramér-Rao下界とは
どんな推定アルゴリズムを工夫しても、雑音があれば推定値は真値の周りでばらつきます。このばらつき（分散）には、物理的に超えられない下限が存在します。それがCramér-Rao下界（CRLB: Cramér-Rao Lower Bound）です。CRLBは「与えられた観測モデルのもとで、不偏推定量が達成しうる分散の理論的な最小値」を与えます。推定量の良し悪しを測る物差しであり、CRLBに張り付いていれば「もうこれ以上良くできない最良の推定」だと言えます。
CRLBはフィッシャー情報量 $I(\varepsilon)$ の逆数として与えられます。
$$ \mathrm{Var}(\hat{\varepsilon}) \ge \frac{1}{I(\varepsilon)}, \quad I(\varepsilon) = -\mathbb{E}\left[\frac{\partial^2 \ln p(\bm{r} \mid \varepsilon)}{\partial \varepsilon^2}\right] $$
ここで $p(\bm{r}\mid\varepsilon)$ は観測 $\bm{r}$ の尤度（パラメータ $\varepsilon$ が与えられたときの確率密度）です。フィッシャー情報量が大きいほど、観測がパラメータについて多くの情報を含み、より正確に推定できることを意味します。
CFO推定のCRLBの導出
既知のプリアンブル $s[n]$（$n=0,\dots,L-1$）を、CFO $\varepsilon$ を介して観測する問題を考えます。観測モデルは
$$ r[n] = s[n]\,e^{j2\pi\varepsilon n/N} + w[n], \quad w[n]\sim\mathcal{CN}(0,\sigma_w^2) $$
です。雑音が独立な循環対称複素ガウス分布なので、対数尤度は（定数を除いて）
$$ \ln p(\bm{r}\mid\varepsilon) = -\frac{1}{\sigma_w^2}\sum_{n=0}^{L-1}\left| r[n] – s[n]e^{j2\pi\varepsilon n/N}\right|^2 $$
です。これを $\varepsilon$ で2回微分してフィッシャー情報量を求めます。途中の見通しを良くするため、信号成分を $\mu[n] = s[n]e^{j2\pi\varepsilon n/N}$ と置きます。$\mu[n]$ の $\varepsilon$ 微分は、肩に $\frac{j2\pi n}{N}$ が降りてきて
$$ \frac{\partial \mu[n]}{\partial \varepsilon} = \frac{j2\pi n}{N}\,s[n]\,e^{j2\pi\varepsilon n/N} = \frac{j2\pi n}{N}\,\mu[n] $$
となります。複素ガウス観測に対するフィッシャー情報量の一般公式
$$ I(\varepsilon) = \frac{2}{\sigma_w^2}\sum_{n=0}^{L-1}\left|\frac{\partial \mu[n]}{\partial\varepsilon}\right|^2 $$
に代入します。$\left|\frac{\partial\mu[n]}{\partial\varepsilon}\right|^2 = \left(\frac{2\pi n}{N}\right)^2 |s[n]|^2$ なので、
$$ I(\varepsilon) = \frac{2}{\sigma_w^2}\sum_{n=0}^{L-1}\left(\frac{2\pi n}{N}\right)^2 |s[n]|^2 = \frac{8\pi^2}{N^2 \sigma_w^2}\sum_{n=0}^{L-1} n^2 |s[n]|^2 $$
を得ます。プリアンブルの各サンプル電力が一定 $|s[n]|^2 = E_s$ だと仮定すると、$\sum_{n=0}^{L-1} n^2 = \frac{(L-1)L(2L-1)}{6}$ を使って
$$ I(\varepsilon) = \frac{8\pi^2 E_s}{N^2\sigma_w^2}\cdot\frac{(L-1)L(2L-1)}{6} $$
となります。したがってCRLBは
$$ \begin{equation} \mathrm{Var}(\hat{\varepsilon}) \ge \frac{N^2 \sigma_w^2}{8\pi^2 E_s}\cdot\frac{6}{(L-1)L(2L-1)} \end{equation} $$
です。ここから重要な3つの性質が読み取れます。第一に、CRLBは $\sigma_w^2/E_s$、すなわちSNRの逆数に比例します。SNRが10倍（10dB）上がれば分散は1/10になります。第二に、サンプル数 $L$ が大きいほど分母の $L^3$ オーダーの項が効いて、分散は急速に小さくなります（およそ $L^3$ に反比例）。長いプリアンブルほど精密に推定できるのです。第三に、$N^2$ に比例しますが、これは $\varepsilon$ をサブキャリア間隔で正規化していることによる見かけ上のスケールで、実周波数 $\Delta f$ の分散に直せばこの $N^2$ は相殺されます。
Moose推定量とCRLBの関係
実は、相関ベースのMoose推定量は、高SNR領域でCRLBにほぼ漸近的に到達することが知られています。直感的には、相関 $P = \sum r[n+D]r^*[n]$ を取って偏角を見る操作が、最尤推定（尤度を最大にする $\varepsilon$ を選ぶ操作）と高SNRで等価になるためです。後のPython実験で、推定量のRMSEがSNRに対してどのように下がり、CRLBの平方根（標準偏差の下限）にどこまで近づくかを実際に確認します。
理論の準備が整いました。導いた推定量と理論限界を、実際のOFDM信号でシミュレーションして確かめましょう。
Pythonでの実装
概念の可視化 — CFOによる星座点の回転
まず、CFOが信号にどう作用するかを直感的に掴むため、QPSK星座点がCFOによって時間とともに回転していく様子を可視化します。日本語ラベルのためにフォント設定を入れます。
import numpy as np
import matplotlib, matplotlib.pyplot as plt
# 日本語フォント設定
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
plt.rcParams["font.family"] = cand
break
plt.rcParams["axes.unicode_minus"] = False
rng = np.random.default_rng(0)
N = 64 # FFTサイズ
eps = 0.05 # 正規化CFO(サブキャリア間隔単位)
n = np.arange(2 * N) # 2シンボル分の時刻
# 固定のQPSKシンボル(+1+1j を例に)を連続送信したと仮定
s = (1 + 1j) / np.sqrt(2) * np.ones(len(n))
r = s * np.exp(1j * 2 * np.pi * eps * n / N) # CFOを掛ける
plt.figure(figsize=(6, 6))
sc = plt.scatter(r.real, r.imag, c=n, cmap="viridis", s=18)
plt.colorbar(sc, label="サンプル番号 $n$")
plt.scatter([s[0].real], [s[0].imag], c="red", s=120,
marker="*", label="本来の星座点")
plt.xlabel("同相成分 I"); plt.ylabel("直交成分 Q")
plt.title(f"CFO($\\varepsilon$={eps})による星座点の回転")
plt.axhline(0, color="gray", lw=0.5); plt.axvline(0, color="gray", lw=0.5)
plt.legend(); plt.axis("equal"); plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()
このグラフから、CFOの本質がはっきり見えます。本来は1点（赤い星）に留まるはずの星座点が、サンプル番号が進む（色が紫→黄に変わる）につれて反時計回りにぐるぐる回っていきます。回転の角速度はちょうど $2\pi\varepsilon/N$ で、$\varepsilon$ が大きいほど速く回ります。OFDM受信機はこの回転を打ち消さないと、復調時に各シンボルの位相がずれて誤判定が起きてしまいます。
同じデータの受信位相を時間に対してプロットすると、回転が「一定の速さ」であることがより直接的に見えます。
受信位相（青）はサンプル番号に対してまっすぐ増えていき、その傾きは理論値 $2\pi\varepsilon/N$ を引いた黒破線とぴったり重なります。つまりCFO推定とは、この「位相の傾き」を測る作業に他なりません。傾きが一定なので、離れた2点の位相差を取れば傾きが分かる——これが次の繰り返し相関のアイデアにつながります。
回転の速さを測る鍵が繰り返し相関でした。次に、繰り返しプリアンブルを作って実際に相関の位相からCFOを推定してみます。
繰り返しプリアンブルとMoose推定量の実装
Schmidl&Cox型の繰り返し構造（前半 $N/2$ と後半 $N/2$ が同一）を持つOFDMプリアンブルを作り、既知CFOとAWGNを与えてMoose法で推定します。
import numpy as np
def make_repeated_preamble(N, rng):
"""前半N/2と後半N/2が同一になるOFDMプリアンブルを生成"""
# 偶数サブキャリアのみにPSKを載せると時間領域で前半=後半になる
X = np.zeros(N, dtype=complex)
even = np.arange(0, N, 2)
# ±1 のBPSKを偶数サブキャリアに配置
X[even] = rng.choice([-1, 1], size=len(even)) * np.sqrt(2)
s = np.fft.ifft(X) * np.sqrt(N) # 時間領域信号(電力正規化)
return s
def add_cfo_awgn(s, eps, N, snr_db, rng):
"""信号にCFOとAWGNを付与"""
n = np.arange(len(s))
r = s * np.exp(1j * 2 * np.pi * eps * n / N) # CFO付与
sig_pow = np.mean(np.abs(s) ** 2)
noise_pow = sig_pow / (10 ** (snr_db / 10))
w = np.sqrt(noise_pow / 2) * (rng.standard_normal(len(s))
+ 1j * rng.standard_normal(len(s)))
return r + w
def moose_estimate(r, D, N, L=None):
"""繰り返し間隔Dの相関位相からCFOを推定(Moose法)"""
if L is None:
L = len(r) - D
P = np.sum(r[D:D + L] * np.conj(r[0:L])) # 相関量
return N / (2 * np.pi * D) * np.angle(P)
上の図は生成したプリアンブルの時間波形（実部・虚部）です。赤い破線で区切った前半 $N/2$ サンプルと後半 $N/2$ サンプルが完全に重なっており、最大差は $10^{-15}$ 程度（数値誤差レベル）です。偶数サブキャリアだけに変調を載せると時間領域で前半＝後半になる、という性質が実測で確認できました。この周期 $D=N/2$ の繰り返し構造がMoose法の土台になります。
上のコードでは3つの関数を用意しました。make_repeated_preamble は偶数サブキャリアだけにBPSKを載せることで、IFFT後の時間信号の前半と後半が一致する（周期 $D=N/2$ の繰り返し）プリアンブルを作ります。add_cfo_awgn は信号に $e^{j2\pi\varepsilon n/N}$ を掛けてCFOを与え、指定SNRのAWGNを加えます。moose_estimate は導出した式(2)そのもので、相関 $\sum r[n+D]r^*[n]$ の偏角に $N/(2\pi D)$ を掛けて $\hat\varepsilon$ を返します。
次に、これらを使って単発の推定がどれくらい合うかを確認します。
rng = np.random.default_rng(1)
N = 64
D = N // 2 # 繰り返し間隔
true_eps = 0.12 # 真のCFO
snr_db = 15
s = make_repeated_preamble(N, rng)
r = add_cfo_awgn(s, true_eps, N, snr_db, rng)
est = moose_estimate(r, D, N, L=N // 2)
print(f"真のCFO : {true_eps:.4f}")
print(f"推定CFO : {est:.4f}")
print(f"推定誤差 : {est - true_eps:+.4f}")
print(f"推定可能範囲: ±{N/(2*D):.2f} (サブキャリア間隔単位)")
実行すると、真のCFO 0.12 に対して推定値はおよそ 0.135（このseedでの実測値）になり、誤差は0.015程度に収まります。SNR 15dBという現実的な条件で、たった32サンプルの相関から小数CFOがかなり正確に測れていることがわかります。また推定可能範囲は $\pm N/2D = \pm 1$ サブキャリア間隔で、真値0.12はこの中に十分収まっています。
なぜ雑音があっても測れるのか、相関の中身を複素平面で見てみましょう。各相関項 $r[n+D]r^*[n]$ をベクトルとして原点から順に足し上げたのが次の図です。
個々の相関項（青い丸）は雑音で大きさも向きもばらつきますが、信号成分はすべて同じ位相 $2\pi\varepsilon D/N$（黒破線方向）を向いているため、足し合わせると赤い合計ベクトル $P$ が長く伸びます。一方、雑音成分は向きがバラバラなので積み上げで打ち消されます。これがコヒーレント積分の威力で、$P$ の偏角（ここでは約0.42 rad）が真の位相0.377 radに近く、そこからCFOが逆算できます。
単発では誤差が偶然に左右されるので、次は多数回試行してRMSEを統計的に評価し、SNR依存性を見ます。
RMSEのSNR依存性とCRLBへの接近
SNRを変えながらモンテカルロ試行でRMSE（二乗平均平方根誤差）を求め、理論限界であるCRLBの平方根と比較します。
import numpy as np
def crlb_std(N, D, L, snr_db):
"""CFO推定のCRLBの平方根(標準偏差の下限)。Es=1, sigma_w^2=1/SNR"""
snr = 10 ** (snr_db / 10)
sigma_w2 = 1.0 / snr # Es=1に正規化
Es = 1.0
I = (8 * np.pi**2 * Es) / (N**2 * sigma_w2) \
* ((L - 1) * L * (2 * L - 1) / 6)
return np.sqrt(1.0 / I)
N = 64
D = N // 2
L = N // 2
true_eps = 0.12
snr_list = np.arange(0, 31, 3)
num_trials = 4000
rmse = []
for snr_db in snr_list:
rng = np.random.default_rng(100 + snr_db)
errs = []
for _ in range(num_trials):
s = make_repeated_preamble(N, rng)
r = add_cfo_awgn(s, true_eps, N, snr_db, rng)
est = moose_estimate(r, D, N, L=L)
errs.append(est - true_eps)
errs = np.array(errs)
rmse.append(np.sqrt(np.mean(errs**2)))
rmse = np.array(rmse)
crlb = np.array([crlb_std(N, D, L, s) for s in snr_list])
plt.figure(figsize=(8, 5))
plt.semilogy(snr_list, rmse, "o-", label="Moose推定量のRMSE(実測)")
plt.semilogy(snr_list, crlb, "k--", label="CRLB(理論限界)")
plt.xlabel("SNR (dB)"); plt.ylabel("CFO推定RMSE(サブキャリア間隔単位)")
plt.title("CFO推定誤差のSNR依存性とCramér-Rao下界")
plt.legend(); plt.grid(True, which="both", alpha=0.3)
plt.tight_layout()
plt.show()
このグラフから、CFO推定の性能の本質が読み取れます。第一に、RMSEはSNRが上がるにつれて直線的（片対数軸で）に下がっていきます。SNRが3dB上がるごとにRMSEがおよそ半分になり（実測でSNR 0dBの約0.070から30dBの約0.0018まで低下）、これはCRLBが $\sigma_w^2$（SNRの逆数）に比例することと整合します。第二に、Moose推定量のRMSE曲線はCRLBの破線とほぼ平行で、ぴったり寄り添って下がっています。相関ベースの推定が理論限界に迫る最良の推定であることがわかります（ここでのRMSEがCRLB線をわずかに下回るのは、Moose法が区間 $n$ と $n+D$ の両方のサンプルを使うため、$L$ サンプル分で計算した破線よりも実効的に多くの観測を活用しているためです）。第三に、低SNR領域ではCRLBが前提とする高SNR近似が崩れ、雑音が強いと相関の偏角が時々大きく外れる（位相のアウトライア）ため、本来は曲線が破線から上に離れていきます。
推定値そのものの分布も見ておきましょう。同じ設定で多数回推定し、ヒストグラムにしたのが次の図です。
左のSNR 5dBでは推定値が真値0.12の周りに幅広く（標準偏差0.04前後）散らばりますが、右のSNR 20dBでは真値の周りに鋭く集中します（標準偏差は1桁小さい）。どちらも平均は真値にほぼ一致しており、Moose推定量がほぼ不偏で、SNRが上がるほど分散だけが小さくなる——CRLBが描く振る舞いそのものが分布の形として現れています。
次に、推定範囲の限界（曖昧さ）が実際にどう現れるかを見ます。真のCFOを範囲外まで振ってみましょう。
推定可能範囲の限界（曖昧さ）の可視化
真のCFOを $-1.5$ から $+1.5$ まで変化させ、推定値がどうなるかをプロットします。$D=N/2$ なので推定範囲は $\pm 1$ のはずです。
import numpy as np
N = 64
D = N // 2
snr_db = 25 # 曖昧さの効果を見るため高SNR
true_eps_list = np.linspace(-1.5, 1.5, 121)
est_list = []
for te in true_eps_list:
rng = np.random.default_rng(7)
s = make_repeated_preamble(N, rng)
r = add_cfo_awgn(s, te, N, snr_db, rng)
est_list.append(moose_estimate(r, D, N, L=N // 2))
est_list = np.array(est_list)
plt.figure(figsize=(8, 5))
plt.plot(true_eps_list, est_list, "o-", ms=3, label="推定CFO")
plt.plot(true_eps_list, true_eps_list, "k--", label="理想(推定=真値)")
plt.axvspan(-1, 1, color="green", alpha=0.1, label="推定可能範囲 ±1")
plt.axvline(1, color="red", ls=":"); plt.axvline(-1, color="red", ls=":")
plt.xlabel("真のCFO $\\varepsilon$"); plt.ylabel("推定CFO $\\hat{\\varepsilon}$")
plt.title("推定可能範囲の限界(曖昧さ $\\pm N/2D$)")
plt.legend(); plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()
このグラフは推定範囲の限界を鮮やかに示します。真のCFOが $-1 < \varepsilon \le 1$ の範囲（緑の帯）にある間は、推定値が理想直線（黒破線）にぴったり重なり、正確に推定できています。ところが $\varepsilon$ が $+1$ を超えた瞬間、推定値は $-1$ 側に折り返し（ラップアラウンド）、まったく違う値を返します。これは $\arg$ 関数が $\pm\pi$ で折り返すために起きる「曖昧さ」そのものです。実システムでは、まずこの範囲を超えないように粗い周波数引き込みを行ってからMoose法で精密化するか、繰り返し間隔 $D$ の異なる複数のプリアンブルで範囲と精度を両立させます。
最後に、繰り返し間隔 $D$ を変えると精度と範囲がどうトレードオフするかを確認します。
繰り返し間隔Dによる精度と範囲のトレードオフ
$D$ を $N/8, N/4, N/2$ と変えて、それぞれのRMSEと推定範囲を比較します。
import numpy as np
N = 128
true_eps = 0.1
snr_list = np.arange(0, 25, 3)
num_trials = 3000
D_list = [N // 8, N // 4, N // 2]
plt.figure(figsize=(8, 5))
for D in D_list:
L = D # 相関に使う長さをDに合わせる
rmse = []
for snr_db in snr_list:
rng = np.random.default_rng(200 + snr_db + D)
errs = []
for _ in range(num_trials):
# Dサンプル周期の繰り返し信号(BPSK系列をD周期で繰り返す)
base = rng.choice([-1, 1], size=D) + 0j
s = np.tile(base, 2) # 2周期分=長さ2D
n = np.arange(len(s))
r = s * np.exp(1j * 2 * np.pi * true_eps * n / N)
sp = np.mean(np.abs(s)**2)
npow = sp / (10 ** (snr_db / 10))
r = r + np.sqrt(npow/2)*(rng.standard_normal(len(s))
+ 1j*rng.standard_normal(len(s)))
est = N/(2*np.pi*D) * np.angle(np.sum(r[D:D+L]*np.conj(r[0:L])))
errs.append(est - true_eps)
rmse.append(np.sqrt(np.mean(np.array(errs)**2)))
plt.semilogy(snr_list, rmse, "o-",
label=f"D={D} (範囲±{N/(2*D):.1f})")
plt.xlabel("SNR (dB)"); plt.ylabel("CFO推定RMSE")
plt.title("繰り返し間隔 $D$ による精度と範囲のトレードオフ")
plt.legend(); plt.grid(True, which="both", alpha=0.3)
plt.tight_layout()
plt.show()
このグラフから、$D$ を巡るトレードオフがはっきり見えます。$D$ が大きいほど（$D=N/2$）、同じSNRでもRMSEが小さく、精密に推定できています。これは位相差 $2\pi\varepsilon D/N$ が大きく、雑音に対して測りやすいからです。一方、凡例に示したように $D=N/2$ では推定範囲が $\pm 1$ と狭く、$D=N/8$ では範囲が $\pm 4$ と広い代わりに同じSNRでのRMSEは大きくなっています。つまり「広く粗く測る」か「狭く精密に測る」かの選択であり、先に理論で述べたトレードオフが数値実験でそのまま確認できました。実システムが複数の $D$ を多段で使う理由がここにあります。
最後に、推定したCFOで実際に信号を補正すると復調がどれだけ回復するかを確かめます。
推定値による補正と復調の回復
プリアンブルから推定した $\hat\varepsilon$ を使って受信信号に逆回転 $e^{-j2\pi\hat\varepsilon n/N}$ を掛け、CFOを打ち消してからFFT復調します。
import numpy as np
N = 64
true_eps = 0.18
snr_db = 25
rng = np.random.default_rng(21)
# データOFDMシンボル(QPSK)にCFOとAWGNを付与
qpsk = (rng.choice([-1, 1], N) + 1j * rng.choice([-1, 1], N)) / np.sqrt(2)
s = np.fft.ifft(qpsk) * np.sqrt(N)
n = np.arange(N)
rx = s * np.exp(1j * 2 * np.pi * true_eps * n / N)
sp = np.mean(np.abs(s) ** 2); npow = sp / (10 ** (snr_db / 10))
rx = rx + np.sqrt(npow / 2) * (rng.standard_normal(N) + 1j * rng.standard_normal(N))
# 別の繰り返しプリアンブルで同じCFOを推定
rng2 = np.random.default_rng(22)
s_pre = make_repeated_preamble(N, rng2)
r_pre = add_cfo_awgn(s_pre, true_eps, N, snr_db, rng2)
eps_hat = moose_estimate(r_pre, N // 2, N, L=N // 2)
# 推定値で逆回転して補正
rx_corr = rx * np.exp(-1j * 2 * np.pi * eps_hat * n / N)
Y_bad = np.fft.fft(rx) / np.sqrt(N)
Y_good = np.fft.fft(rx_corr) / np.sqrt(N)
print(f"推定CFO : {eps_hat:.4f} (真値 {true_eps})")
print(f"補正前EVM : {np.sqrt(np.mean(np.abs(Y_bad - qpsk)**2)):.3f}")
print(f"補正後EVM : {np.sqrt(np.mean(np.abs(Y_good - qpsk)**2)):.3f}")
左の補正前は、CFO $\varepsilon=0.18$ によって星座点が回転・拡散し、4つの判定点（赤い星）がまったく見分けられません（EVM 0.67）。右の補正後は、プリアンブルから推定した $\hat\varepsilon \approx 0.178$ で逆回転を掛けただけで、受信点が4つの判定点にぴたりと戻り、EVMは0.05まで1桁以上改善します。推定がほぼ正確（誤差0.002程度）なので、補正後の劣化はもとのAWGNだけになっています。CFO推定→補正という一連の流れが、復調を救う決定的な処理であることがこの一枚でわかります。
まとめ
本記事では、搬送波周波数オフセット（CFO）の推定について、理論の導出から実装・検証まで通して解説しました。
- CFOの正体: 送受信機の発振器ずれやドップラーによる残留周波数差で、受信信号に $r[n] = s[n]e^{j2\pi\varepsilon n/N}$ という回転する複素指数として現れます。OFDMではサブキャリアの直交性を壊し、振幅減衰とICIを引き起こします
- 繰り返し相関の原理: 同一信号を間隔 $D$ で繰り返すと、CFOが一定位相差 $2\pi\varepsilon D/N$ として現れます。相関 $\sum r[n+D]r^*[n]$ を取ると信号が同位相で積み上がり、雑音が平均化されます
- Moose推定量: 相関の偏角から $\hat\varepsilon = \frac{N}{2\pi D}\arg\{\sum r[n+D]r^*[n]\}$ で推定できます。Schmidl&Cox法はこれをOFDMプリアンブル設計と組み合わせた実用形です
- 曖昧さ $\pm N/2D$: $\arg$ の値域が $(-\pi,\pi]$ のため、推定範囲は $\pm N/2D$ に限られます。$D$ が大きいほど精密ですが範囲は狭く、トレードオフがあります
- CRLBへの接近: 推定分散の理論限界はSNRの逆数とプリアンブル長 $L$ の3乗に反比例し、Moose推定量は高SNRでこの限界にほぼ到達することを数値実験で確認しました
CFO推定は、OFDM受信機の同期処理の中核です。タイミング同期（シンボルの先頭を見つける）、CFOの粗推定と精推定、そしてチャネル推定という一連の流れの中で、本記事のプリアンブル相関法は精推定の要として機能します。これらを学んだ上で、整数CFOの推定、位相雑音やサンプリング周波数オフセット（SFO）への拡張へと進むと、実際の受信機がどう動いているのかがさらに鮮明に見えてきます。
次のステップとして、以下の記事も参考にしてください。
- OFDM（直交周波数分割多重）の原理 — CFOが壊す直交性そのものの理解を深める
- IQ直交変調の仕組み — 複素ベースバンド信号と回転因子の関係
- AWGNチャネルモデル — 推定分散とCRLBの背景にある雑音モデル