低軌道（LEO）の衛星から地上に降ってくる電波の周波数は、頭上を通過するわずか数分間に最大で $\pm 300\,\mathrm{kHz}$ も揺れ動きます。さらに厄介なことに、衛星が天頂付近を通過する瞬間にはドップラー周波数そのものが急峻に変化し、その変化率（Hz/s）は数キロHz/sにも達します。受信機は、こうした揺れ動く搬送波の位相を「ぴったり」追跡し続けなければ、復調器に渡すべき同相成分 $I$ と直交成分 $Q$ を正しく分離できません。この「揺れ動く搬送波に追従し続ける仕事」を担うのが、搬送波同期（carrier synchronization）であり、その心臓部に位置するのがPLL（Phase-Locked Loop）とCostasループです。
搬送波同期を理解すると、たとえば次のような実問題に踏み込めるようになります。LEOコンステレーション端末でなぜ三次（タイプIII）ループが必須なのか、なぜ受信機の起動時には FLL（Frequency-Locked Loop）から PLL への切替シーケンスが組まれるのか、そしてループ帯域 $B_L$ を広げると追跡力は強くなるのに位相ジッタ分散 $\sigma_\phi^2$ が悪化するという「ダイナミクス vs 雑音」のトレードオフはどう設計するか。これらは衛星通信端末、GNSS受信機、深宇宙探査機の地上系などで日々向き合うエンジニアリング問題そのものです。
本記事の内容
- 位相検出器（PD）→ ループフィルタ（LF）→ 数値制御発振器（NCO）の3要素構成と直感
- 線形化PLLモデルと閉ループ伝達関数 $H(s)$ の導出
- 二次PLLの減衰比 $\zeta$ と自然周波数 $\omega_n$、定常誤差表
- 三次PLLがドップラー変化率（ランプ位相入力）を定常誤差ゼロで追跡できる理由
- ループ帯域 $B_L$ と位相ジッタ分散 $\sigma_\phi^2 = B_L / (C/N_0)$ のトレードオフ
- Costasループ（BPSK/QPSK向け）の誤差信号と $\pi$ 位相曖昧
- Pythonによる離散時間PLLの実装と、ドップラー変化率追跡の比較実験
- FLL→PLL切替シーケンスの実装例
前提知識
この記事を読む前に、以下の記事を読んでおくと理解がスムーズです。
なぜ搬送波同期が必要か（直感）
衛星から送られてくる信号を素朴に書くと、次のような形をしています。
$$ s(t) = A \cos\bigl(2\pi f_c t + \phi(t)\bigr) + n(t) $$
ここで $f_c$ は搬送波周波数（たとえば Ku帯なら $12\,\mathrm{GHz}$）、$\phi(t)$ には変調による位相変化（BPSKなら $0$ または $\pi$）と、ドップラー効果による位相のドリフトの両方が含まれます。
受信機は、この信号にローカル発振器の正弦波 $\cos(2\pi \hat f_c t + \hat\phi)$ と $-\sin(2\pi \hat f_c t + \hat\phi)$ をそれぞれ乗算して低域通過することで、同相成分 $I$ と直交成分 $Q$ を取り出します。問題は、ローカル発振器の位相 $\hat\phi$ が真の搬送波位相 $\phi$ と $\Delta\phi$ だけずれていると、本来 $I$ チャネルに乗るべき信号エネルギーの一部が $Q$ チャネルに漏れ込むことです。具体的に、$I$ 出力には係数 $\cos(\Delta\phi)$ が、$Q$ 出力には $\sin(\Delta\phi)$ が掛かります。
$\Delta\phi = 30^\circ$ では $\cos(30^\circ) \approx 0.87$ ですから、信号エネルギーが約 $1.3\,\mathrm{dB}$ 失われ、雑音は逆に $Q$ 側から漏れ込みます。$\Delta\phi = 90^\circ$ になれば信号は完全に $Q$ 側に移ってしまい、$I$ 復調器は何も拾えません。これがまず「位相を合わせる」必要性です。
しかし話はここで終わりません。LEO衛星の信号は時間と共に周波数そのものがドップラー効果で動きます。さらに、ドップラー周波数を時間で微分した「変化率」も無視できません。LEO衛星が天頂を通過する瞬間、視線速度はサインカーブの最大傾き地点を通り、ドップラー変化率は毎秒数キロHzに達します。受信機はこの「位相だけでなく周波数も、しかも周波数の変化率まで」追跡できなければ、いずれ位相ロックを失います。
ここで自然に湧き上がる疑問があります。位相を追えるなら自動的に周波数も追えるのではないか？ 周波数を追えるならその変化率も追えるのではないか？ 答えは「ループの次数によって追跡能力が決まる」という非常にきれいな数学的事実にあります。これを以降で順に紐解いていきましょう。
PLLの基本構造
3つの構成要素
PLLは、フィードバック制御理論の最も美しい応用例の一つです。最も基本的な構成は3つの要素から成ります。
- 位相検出器（Phase Detector, PD）: 入力信号の位相 $\phi$ と推定位相 $\hat\phi$ の差 $e = \phi – \hat\phi$ を出力する
- ループフィルタ（Loop Filter, LF）: 検出された位相誤差を平滑化し、NCOへの制御信号を生成する
- 数値制御発振器（Numerically Controlled Oscillator, NCO）: 制御信号に応じて出力位相を更新する。デジタル実装では数値で位相を蓄積する
このループは、ちょうど自転車のハンドルで方向を補正する動作によく似ています。位相検出器は「目標との進路ずれ」を計測し、ループフィルタは「揺れすぎないように補正量を平滑化」し、NCOは「実際にハンドルを切る」役割を担います。重要なのは、PDで瞬時の誤差を見るだけでは振動してしまうため、LFで時間方向に平均化することと、過去の誤差の積分項を持つことで定常的なバイアスを除去することです。
ブロック図の流れ
入力信号と NCO 出力を位相検出器に入れると、その出力は位相誤差 $e(t)$ に比例します。これをループフィルタ $F(s)$ に通して制御信号 $v(t)$ を作り、NCO は $v(t)$ を周波数偏差として位相を積分します。NCO の位相積分が「もう1つの積分器」として働くため、ループには少なくとも1つの自然な積分が組み込まれています。
ループフィルタが「比例のみ（P制御）」であれば、NCOの積分と合わせてループは一次ループになります。ループフィルタが「比例+積分（PI制御）」であれば、NCO の積分と合わせて積分器が2つ並び、二次ループになります。さらに「比例+積分+積分（PII制御）」を加えれば三次ループです。後で見るように、この「ループに含まれる積分器の数」がそのまま「追跡できる位相入力の次数」を決めます。
ここまでで PLL の構造を直感的に把握できました。次に、これを数式モデルで定量的に解析していきましょう。
線形化PLLの伝達関数
位相検出器の線形化
位相検出器の物理的な実装は、入力 $\sin(\phi)$ と NCO 出力 $\cos(\hat\phi)$ を乗算する形になることが多く、その出力は
$$ \sin(\phi)\cos(\hat\phi) = \frac{1}{2}\bigl[\sin(\phi – \hat\phi) + \sin(\phi + \hat\phi)\bigr] $$
の高域成分を低域通過で除いた $\frac{1}{2}\sin(\phi – \hat\phi)$ になります。位相誤差 $\Delta\phi = \phi – \hat\phi$ が小さいとき、テイラー展開で $\sin(\Delta\phi) \approx \Delta\phi$ とできるため、PD の出力は位相誤差に比例する線形要素 $K_d \cdot \Delta\phi$ として扱えます。これが「線形化PLL」の出発点です。
実際の受信機で位相誤差が大きく揺れているとき、$\sin$ の非線形性が効いて引き込み（acquisition）特性を悪化させますが、いったんロックされた後の追跡特性（tracking）は線形モデルで十分良く近似できます。本記事では主に追跡特性に焦点を当てます。
閉ループ伝達関数
線形化されたPLLは、入力位相 $\Phi(s)$ から推定位相 $\hat\Phi(s)$ への線形システムとして書けます。PDのゲインを $K_d$、ループフィルタを $F(s)$、NCOのゲインを $K_o/s$（積分器）とすると、開ループ伝達関数は
$$ G(s) = \frac{K_d K_o F(s)}{s} $$
になります。閉ループの位相伝達関数 $H(s) = \hat\Phi(s)/\Phi(s)$ は、フィードバック系の標準的な公式から導出します。まず $\hat\Phi(s) = G(s)\bigl(\Phi(s) – \hat\Phi(s)\bigr)$ を整理すると、
$$ \hat\Phi(s)\bigl(1 + G(s)\bigr) = G(s)\Phi(s) $$
になり、両辺を $(1+G(s))$ で割って $\Phi(s)$ で除すると、
$$ H(s) = \frac{G(s)}{1 + G(s)} = \frac{K_d K_o F(s)}{s + K_d K_o F(s)} $$
が得られます。一方、位相誤差伝達関数 $H_e(s) = (\Phi(s) – \hat\Phi(s))/\Phi(s)$ は、$H_e = 1 – H$ より
$$ H_e(s) = \frac{s}{s + K_d K_o F(s)} $$
になります。$H_e(s)$ は「入力位相のどの成分が誤差として残るか」を表す重要な関数です。
一次・二次・三次ループ
ループフィルタ $F(s)$ の選択でループの次数が決まります。一次ループは $F(s) = 1$（純粋な比例ゲイン）で、閉ループ伝達関数は
$$ H(s) = \frac{K}{s + K}, \quad K = K_d K_o $$
の単純な一次低域通過になります。二次ループは $F(s) = K_p + K_i/s$（PI 制御）で、
$$ H(s) = \frac{2\zeta \omega_n s + \omega_n^2}{s^2 + 2\zeta \omega_n s + \omega_n^2} $$
の標準的な二次系になります。ここで $\omega_n = \sqrt{K_d K_o K_i}$ は自然周波数、$\zeta = \frac{K_p}{2}\sqrt{K_d K_o / K_i}$ は減衰比です。$\zeta = 1/\sqrt{2} \approx 0.707$ がベッセル特性（過渡応答とジッタのバランス最適）、$\zeta = 1$ が臨界制動です。三次ループは $F(s) = K_p + K_i/s + K_{ii}/s^2$（PII 制御）として実現されます。
これらの設計パラメータが、後で見る「定常誤差」「ループ帯域」「ジッタ分散」を決めます。ここまで伝達関数の構造が見えてきました。次に、PLLがどんな位相入力に対してどんな誤差を残すのか、定常誤差を整理しましょう。
定常誤差とループ次数
ステップ・ランプ・パラボラ入力
衛星通信受信機が直面する位相入力は、典型的に次の3パターンに分解できます。
- 位相ステップ $\phi(t) = \Delta\phi \cdot u(t)$（突然の位相飛び。ラプラス変換: $\Delta\phi/s$）
- 位相ランプ $\phi(t) = \Delta\omega \cdot t$（一定周波数オフセット、つまり一定のドップラー周波数。$\Delta\omega/s^2$）
- 位相パラボラ $\phi(t) = \frac{1}{2}\dot\omega \cdot t^2$（一定の周波数変化率、つまり一定のドップラー変化率。$\dot\omega/s^3$）
LEO衛星の通過では、ドップラー周波数自体が「ランプ」のように変化する区間と、その変化率自体（パラボラ位相入力）が無視できない区間が両方存在します。
終値定理による定常誤差
位相誤差の定常値は、ラプラス変換の終値定理 $e_{\infty} = \lim_{s \to 0} s \cdot E(s) = \lim_{s \to 0} s \cdot H_e(s) \Phi(s)$ で計算できます。一次ループに位相ランプ $\Phi(s) = \Delta\omega/s^2$ を入れて、$H_e(s) = s/(s+K)$ を代入すると、
$$ e_{\infty} = \lim_{s \to 0} s \cdot \frac{s}{s + K} \cdot \frac{\Delta\omega}{s^2} $$
になります。$s$ を約分して整理すると、
$$ e_{\infty} = \lim_{s \to 0} \frac{\Delta\omega}{s + K} = \frac{\Delta\omega}{K} $$
となり、ゼロにならないことが分かります。すなわち一次ループはドップラー周波数オフセットに対して定常位相誤差を残します。
ここで「ループの次数」と「追跡できる入力次数」の関係を一気に表にまとめると次のようになります。
| ループ次数 | 含まれる積分器 | 位相ステップ | 位相ランプ（周波数オフセット） | 位相パラボラ（周波数変化率） |
|---|---|---|---|---|
| 一次（タイプI） | 1 | 0 | 有限 | 発散 |
| 二次（タイプII） | 2 | 0 | 0 | 有限 |
| 三次（タイプIII） | 3 | 0 | 0 | 0 |
この表は PLL設計の核 です。LEO衛星のドップラー変化率が無視できる程度なら二次ループで十分ですが、変化率が大きく、それを定常誤差ゼロで追跡したいなら、三次ループが必要になります。これが「LEO衛星端末になぜ三次PLLが採用されるのか」の数学的な答えです。
二次ループでパラボラ入力を入れた場合
二次ループでパラボラ入力を試してみましょう。$H_e(s) = s^2/(s^2 + 2\zeta\omega_n s + \omega_n^2)$ と $\Phi(s) = \dot\omega/s^3$ を終値定理に入れると、
$$ e_{\infty} = \lim_{s \to 0} s \cdot \frac{s^2}{s^2 + 2\zeta\omega_n s + \omega_n^2} \cdot \frac{\dot\omega}{s^3} $$
になります。$s$ の三乗をすべて約分すると、
$$ e_{\infty} = \lim_{s \to 0} \frac{\dot\omega}{s^2 + 2\zeta\omega_n s + \omega_n^2} = \frac{\dot\omega}{\omega_n^2} $$
になります。すなわち二次ループでドップラー変化率を追跡すると、$\dot\omega/\omega_n^2$ の定常位相誤差が残り続けます。$\dot\omega$ が大きい LEO 通過時には、これが復調品質を大きく劣化させます。
たとえば LEO 衛星の Ku 帯通過で $\dot f = 1\,\mathrm{kHz/s}$ 程度のドップラー変化率を想定し、ループの $\omega_n = 2\pi \cdot 10\,\mathrm{Hz}$ で設計すると、定常位相誤差は $\dot\omega/\omega_n^2 = (2\pi \cdot 1000) / (2\pi \cdot 10)^2 \approx 1.59\,\mathrm{rad} \approx 91^\circ$ にもなり、完全にロックを外します。これは机上で予測できる「設計失敗」です。
定常誤差の話で「追跡能力」が掴めました。次は、雑音耐性の側面から、ループ帯域とジッタのトレードオフを見ていきます。
ループ帯域と位相ジッタのトレードオフ
ループ帯域 $B_L$ の定義
PLL の雑音特性を一つのスカラーで表す指標が等価ループ帯域 $B_L$ です。これは閉ループ伝達関数 $H(s)$ の振幅二乗を全周波数で積分したもので、定義式は
$$ B_L = \int_0^{\infty} |H(j 2\pi f)|^2 \, df \quad [\mathrm{Hz}] $$
になります。物理的には「PLLが入力の白色位相雑音に対してどれだけの帯域で雑音を通すか」を表します。二次ループの場合、留数計算と部分分数分解で積分を実行すると
$$ B_L = \frac{\omega_n}{2}\left(\zeta + \frac{1}{4\zeta}\right) $$
になります。$\zeta = 0.707$ では $B_L \approx 0.53\omega_n / (2\pi)$ になります。
位相ジッタ分散
PLLの出力位相は、ループ帯域内に落ちる雑音電力で揺らぎます。受信信号の搬送波対雑音電力密度比を $C/N_0$（単位は Hz、つまり $1\,\mathrm{Hz}$ 帯域あたりの SNR）とすると、PLLの位相ジッタ分散は
$$ \sigma_\phi^2 = \frac{B_L}{C/N_0} \quad [\mathrm{rad}^2] $$
になります。これはPLL設計で最も使われる公式の一つで、極めてシンプルな形をしています。たとえば $C/N_0 = 40\,\mathrm{dBHz}$（典型的な衛星受信）、$B_L = 10\,\mathrm{Hz}$ では $\sigma_\phi^2 = 10/10^4 = 10^{-3}\,\mathrm{rad}^2$、つまり RMS位相ジッタは $\sqrt{10^{-3}} \approx 32\,\mathrm{mrad} = 1.8^\circ$ になります。
ダイナミクス vs ジッタのトレードオフ
ここから直接的な設計トレードオフが現れます。$B_L$ を広く取るとループは速いダイナミクス（大きなドップラー変化率）を追えますが、ジッタが増えます。逆に $B_L$ を狭く取るとジッタは小さくなりますが、追跡が遅くなり、急峻なドップラー変化に追従できずロックを失います。
LEO衛星端末では、$B_L$ をたとえば $10\,\mathrm{Hz}$ から $50\,\mathrm{Hz}$ の範囲で衛星通過の状況に応じて動的に切り替える「適応的ループ帯域制御」がよく使われます。深宇宙探査機の受信では $C/N_0$ が非常に低いため、$B_L$ を $0.1\,\mathrm{Hz}$ 以下に絞り、その代わりドップラー予測モデルでループ外から大まかな周波数補正を与えて、PLL自体は微小な残留誤差だけを追跡する設計が採られます。
このトレードオフを実際に Python で可視化してみると、感覚がさらに掴めます。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# C/N0 を 30〜55 dBHz でスイープし、各BLでの位相ジッタを計算
CN0_dB = np.linspace(30, 55, 100)
CN0 = 10**(CN0_dB / 10)
BL_list = [1, 5, 20, 100]
fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))
for BL in BL_list:
sigma_phi = np.sqrt(BL / CN0) * 180 / np.pi # deg
axes[0].semilogy(CN0_dB, sigma_phi, linewidth=2, label=f'$B_L$ = {BL} Hz')
axes[0].axhline(15, color='gray', linestyle=':', label='Typical lock threshold ~15deg')
axes[0].set_xlabel('C/N0 [dBHz]')
axes[0].set_ylabel('RMS Phase Jitter [deg]')
axes[0].set_title('Phase Jitter vs C/N0 for various $B_L$')
axes[0].legend()
axes[0].grid(alpha=0.3, which='both')
# 二次ループの周波数応答
zeta = 0.707
for BL in [5, 20, 100]:
wn = BL * 8 * zeta / (4 * zeta**2 + 1)
num = [2*zeta*wn, wn**2]
den = [1, 2*zeta*wn, wn**2]
sys = signal.TransferFunction(num, den)
w, mag, _ = signal.bode(sys, np.logspace(-1, 3, 500))
axes[1].semilogx(w/(2*np.pi), mag, linewidth=2, label=f'$B_L$ = {BL} Hz')
axes[1].axhline(-3, color='gray', linestyle=':')
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('|H(j2$\\pi$f)| [dB]')
axes[1].set_title('2nd-order PLL Closed-loop Response')
axes[1].legend()
axes[1].grid(alpha=0.3, which='both')
plt.tight_layout()
plt.savefig('pll_jitter_tradeoff.png', dpi=130, bbox_inches='tight')
plt.show()
左図から、同じ $C/N_0$ でも $B_L$ を $1\,\mathrm{Hz}$ に絞ると $B_L = 100\,\mathrm{Hz}$ より位相ジッタが $10$ 倍（$\sqrt{100}$ 倍）小さくなることが読み取れます。一方、右図の閉ループ応答を見ると $B_L$ が広いほど高周波の位相成分も通すため、ドップラー変化率のような高周波成分も追従できます。これがまさに「狭い $B_L$ で雑音は減るが、速い変動には追従できない」というトレードオフの可視化です。実用では、たとえば $C/N_0 = 35\,\mathrm{dBHz}$ で安定ロック（ジッタ $15^\circ$ 以下）を取るためには $B_L \lesssim 20\,\mathrm{Hz}$ に抑える必要がある、といった設計指針が一目で得られます。
ジッタの話が出たところで、ディジタル通信受信機で頻繁に登場する Costas ループに話を進めましょう。なぜ通常のPLLではダメなのか、その理由から見ていきます。
Costasループ（BPSK/QPSK向け）
BPSK信号でPLLが使えない理由
BPSK 信号は、データビットに応じて搬送波位相を $0$ または $\pi$ で変調します。すなわち送信信号は
$$ s(t) = A \cdot d(t) \cdot \cos(2\pi f_c t + \phi) $$
の形で、ここで $d(t) \in \{-1, +1\}$ がデータです。この信号を普通の PLL（PDが $\sin(\phi – \hat\phi)$ を出すタイプ）に入れると、データの符号反転が位相 $\pi$ の飛びとして PD に見え、PLL は安定にロックできません。
ここで Costas は「データに依存しない位相誤差信号を作る」アイデアを提案しました。それがCostasループです。
Costasループの誤差信号
Costas ループは、入力信号を NCO 出力の同相成分と直交成分の両方で乗算し、それぞれを低域通過した後の I, Q を用いて位相誤差を作ります。具体的には、
$$ I(t) = A \cdot d(t) \cos(\Delta\phi) + n_I(t) $$
$$ Q(t) = A \cdot d(t) \sin(\Delta\phi) + n_Q(t) $$
になります。$d(t) \in \{-1, +1\}$ なので、$I \cdot Q$ という積を計算します。積を素直に展開すると
$$ I \cdot Q = A^2 d^2(t) \sin(\Delta\phi)\cos(\Delta\phi) + (\text{雑音項}) $$
になり、ここで倍角公式 $\sin(\Delta\phi)\cos(\Delta\phi) = \frac{1}{2}\sin(2\Delta\phi)$ と $d^2(t) = 1$ を使うと
$$ I \cdot Q = \frac{A^2}{2}\sin(2\Delta\phi) + (\text{雑音項}) $$
となり、データに依存しない位相誤差信号 $\sin(2\Delta\phi)$ が得られます。小さい位相誤差では $\sin(2\Delta\phi) \approx 2\Delta\phi$ なので、これはPDゲインが2倍の線形PDとして扱えます。
$\pi$ 位相曖昧（phase ambiguity）
しかし $\sin(2\Delta\phi)$ という形は、$\Delta\phi$ と $\Delta\phi + \pi$ の両方で同じ値（ゼロ）を取ります。したがって Costas ループは $\pi$ ラジアンの位相曖昧性を持ち、ロックされたときに $\Delta\phi = 0$ か $\Delta\phi = \pi$ かを区別できません。この曖昧性はビットの極性反転に等価であり、フレーム同期や差動符号化（DBPSK）で解決します。
QPSK では $\sin(2\Delta\phi)$ ではなく $\sin(4\Delta\phi)$ を作るため、$\pi/2$ ラジアンの4値曖昧性が生じます。これも同様にグレイ符号と差動符号化で吸収するのが標準です。
Costas ループの誤差信号設計まで来ました。理論パートはここまでにして、いよいよ Python で離散時間 PLL を実装し、ドップラー変化率追跡の実験をしていきましょう。
Pythonで離散PLLを実装する
離散時間化と係数設計
連続時間の二次PLLを離散化するとき、最もシンプルな実装はループフィルタを後退オイラー法で離散化し、NCOを位相累算器で実現する形です。サンプリング周波数を $f_s$、自然周波数を $\omega_n$、減衰比を $\zeta$ とすると、Gardnerの教科書で標準的に紹介されている係数は
$$ K_p = \frac{8 \zeta \omega_n / f_s}{1 + 2\zeta\omega_n/f_s + (\omega_n/f_s)^2} $$
$$ K_i = \frac{4 (\omega_n/f_s)^2}{1 + 2\zeta\omega_n/f_s + (\omega_n/f_s)^2} $$
になります。これらの係数を使うと、連続時間ループに良く対応した離散時間挙動が得られます。
それでは、ドップラー変化率を持つ複素ベースバンド信号を生成し、二次PLLで追跡する実験をやってみましょう。
import numpy as np
import matplotlib.pyplot as plt
# 信号パラメータ
fs = 10_000.0 # サンプリング周波数 [Hz]
T = 5.0 # 観測時間 [s]
N = int(fs * T)
t = np.arange(N) / fs
# 真のドップラー: 初期周波数 + 変化率
f0 = 200.0 # 初期周波数オフセット [Hz]
fdot = 50.0 # ドップラー変化率 [Hz/s]
phi_true = 2 * np.pi * (f0 * t + 0.5 * fdot * t**2)
# 受信複素ベースバンド信号 (雑音付き, 純CW)
CN0_dBHz = 45.0
A = 1.0
N0 = A**2 / (10**(CN0_dBHz / 10))
noise_std = np.sqrt(N0 * fs / 2)
np.random.seed(0)
noise = noise_std * (np.random.randn(N) + 1j * np.random.randn(N))
rx = A * np.exp(1j * phi_true) + noise
# 二次PLL
def pll_second_order(rx, fs, BL, zeta=0.707):
wn = BL * 8 * zeta / (4 * zeta**2 + 1)
Ts = 1 / fs
Kp = (8 * zeta * wn * Ts) / (1 + 2 * zeta * wn * Ts + (wn * Ts)**2)
Ki = (4 * (wn * Ts)**2) / (1 + 2 * zeta * wn * Ts + (wn * Ts)**2)
N = len(rx)
freq_hat = np.zeros(N)
err = np.zeros(N)
integ = 0.0
phase = 0.0
for n in range(N):
prod = rx[n] * np.exp(-1j * phase)
e = np.arctan2(prod.imag, prod.real)
err[n] = e
integ += Ki * e
freq_update = Kp * e + integ
phase += freq_update
freq_hat[n] = freq_update * fs / (2 * np.pi)
return freq_hat, err
freq_20, err_20 = pll_second_order(rx, fs, BL=20.0)
freq_100, err_100 = pll_second_order(rx, fs, BL=100.0)
f_true = f0 + fdot * t
fig, axes = plt.subplots(2, 2, figsize=(13, 8))
axes[0, 0].plot(t, f_true, 'k-', linewidth=2, label='True Doppler')
axes[0, 0].plot(t, freq_20, 'b-', alpha=0.7, label='$B_L$=20Hz')
axes[0, 0].plot(t, freq_100, 'r-', alpha=0.5, label='$B_L$=100Hz')
axes[0, 0].set_xlabel('Time [s]'); axes[0, 0].set_ylabel('Frequency [Hz]')
axes[0, 0].set_title('Frequency Tracking by 2nd-order PLL')
axes[0, 0].legend(); axes[0, 0].grid(alpha=0.3)
axes[0, 1].plot(t, f_true - freq_20, 'b-', alpha=0.7, label='$B_L$=20Hz')
axes[0, 1].plot(t, f_true - freq_100, 'r-', alpha=0.5, label='$B_L$=100Hz')
axes[0, 1].axhline(0, color='k', linestyle=':')
axes[0, 1].set_xlabel('Time [s]'); axes[0, 1].set_ylabel('Freq Error [Hz]')
axes[0, 1].set_title('Tracking Error (steady-state bias remains)')
axes[0, 1].legend(); axes[0, 1].grid(alpha=0.3)
axes[1, 0].plot(t, err_20 * 180 / np.pi, 'b-', alpha=0.7, label='$B_L$=20Hz')
axes[1, 0].plot(t, err_100 * 180 / np.pi, 'r-', alpha=0.5, label='$B_L$=100Hz')
axes[1, 0].set_xlabel('Time [s]'); axes[1, 0].set_ylabel('Phase Error [deg]')
axes[1, 0].set_title('Phase Detector Output')
axes[1, 0].legend(); axes[1, 0].grid(alpha=0.3)
mask = t > 2.0
axes[1, 1].hist(err_20[mask] * 180/np.pi, bins=50, alpha=0.6, label='$B_L$=20Hz', color='b')
axes[1, 1].hist(err_100[mask] * 180/np.pi, bins=50, alpha=0.6, label='$B_L$=100Hz', color='r')
axes[1, 1].set_xlabel('Phase Error [deg]'); axes[1, 1].set_ylabel('Count')
axes[1, 1].set_title('Phase Jitter Distribution')
axes[1, 1].legend(); axes[1, 1].grid(alpha=0.3)
plt.tight_layout()
plt.savefig('pll_second_order_tracking.png', dpi=130, bbox_inches='tight')
plt.show()
print(f"BL=20Hz freq bias: {np.mean(f_true[mask] - freq_20[mask]):.3f} Hz, "
f"phase jitter RMS: {np.std(err_20[mask])*180/np.pi:.3f} deg")
print(f"BL=100Hz freq bias: {np.mean(f_true[mask] - freq_100[mask]):.3f} Hz, "
f"phase jitter RMS: {np.std(err_100[mask])*180/np.pi:.3f} deg")
このグラフから二次PLLの限界が明確に見て取れます。左上のサブプロットでは、真のドップラー（直線的に増加）に対して $B_L=20\,\mathrm{Hz}$ と $B_L=100\,\mathrm{Hz}$ のいずれも一見追跡できているように見えますが、右上の周波数誤差を見ると、両者ともゼロに収束しない定常バイアスを残しています。これは前節で導いた「二次ループのパラボラ入力に対する定常誤差 $\dot\omega/\omega_n^2$」そのものです。
左下の位相誤差も同様にゼロ平均にはならず、定常的にオフセットを持っています。$B_L$ を広げる（赤）と追跡力は強くなり定常誤差は小さくなりますが、右下のヒストグラムを見れば分かるように位相ジッタの広がりは大幅に増加しています。これがダイナミクス vs ジッタのトレードオフの実物です。
三次PLLで定常誤差をゼロに
二次ループでは消せない定常バイアスを、三次ループで除去してみましょう。
import numpy as np
import matplotlib.pyplot as plt
def pll_third_order(rx, fs, BL):
"""三次PLL: 周波数変化率も追跡(経験的ゲイン)"""
Ts = 1 / fs
wn = BL * 1.2 # 経験的に良好な自然周波数
K1 = 2.4 * wn * Ts
K2 = 1.1 * (wn * Ts)**2
K3 = (wn * Ts)**3
N = len(rx)
freq_hat = np.zeros(N)
err = np.zeros(N)
v1 = 0.0 # 周波数積分器
v2 = 0.0 # 周波数変化率積分器
phase = 0.0
for n in range(N):
prod = rx[n] * np.exp(-1j * phase)
e = np.arctan2(prod.imag, prod.real)
err[n] = e
v2 += K3 * e
v1 += K2 * e + v2
update = K1 * e + v1
phase += update
freq_hat[n] = update * fs / (2 * np.pi)
return freq_hat, err
# 信号生成 (独立実行用)
fs = 10_000.0; T = 5.0; N = int(fs*T); t = np.arange(N)/fs
f0 = 200.0; fdot = 50.0
phi_true = 2*np.pi*(f0*t + 0.5*fdot*t**2)
np.random.seed(0)
A = 1.0; N0 = A**2 / (10**(45/10))
noise = np.sqrt(N0*fs/2)*(np.random.randn(N)+1j*np.random.randn(N))
rx = A*np.exp(1j*phi_true) + noise
f_true = f0 + fdot*t
def pll2(rx, fs, BL, zeta=0.707):
wn = BL * 8 * zeta / (4 * zeta**2 + 1); Ts = 1/fs
Kp = (8*zeta*wn*Ts)/(1+2*zeta*wn*Ts+(wn*Ts)**2)
Ki = (4*(wn*Ts)**2)/(1+2*zeta*wn*Ts+(wn*Ts)**2)
N = len(rx); freq = np.zeros(N); integ = 0.0; phase = 0.0
for n in range(N):
prod = rx[n]*np.exp(-1j*phase)
e = np.arctan2(prod.imag, prod.real)
integ += Ki*e
upd = Kp*e + integ
phase += upd
freq[n] = upd*fs/(2*np.pi)
return freq
freq_20 = pll2(rx, fs, 20.0)
freq_3rd, err_3rd = pll_third_order(rx, fs, BL=20.0)
fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))
axes[0].plot(t, f_true, 'k-', linewidth=2, label='True')
axes[0].plot(t, freq_20, 'b-', alpha=0.6, label='2nd-order')
axes[0].plot(t, freq_3rd, 'g-', alpha=0.7, label='3rd-order')
axes[0].set_xlabel('Time [s]'); axes[0].set_ylabel('Frequency [Hz]')
axes[0].set_title('2nd vs 3rd-order PLL (BL=20Hz)')
axes[0].legend(); axes[0].grid(alpha=0.3)
axes[1].plot(t, f_true - freq_20, 'b-', alpha=0.7, label='2nd-order')
axes[1].plot(t, f_true - freq_3rd, 'g-', alpha=0.7, label='3rd-order')
axes[1].axhline(0, color='k', linestyle=':')
axes[1].set_xlabel('Time [s]'); axes[1].set_ylabel('Frequency Error [Hz]')
axes[1].set_title('3rd-order eliminates steady-state bias')
axes[1].legend(); axes[1].grid(alpha=0.3)
axes[1].set_ylim([-20, 20])
plt.tight_layout()
plt.savefig('pll_third_order_comparison.png', dpi=130, bbox_inches='tight')
plt.show()
mask = t > 2.0
print(f"2nd-order freq bias: {np.mean(f_true[mask] - freq_20[mask]):.3f} Hz")
print(f"3rd-order freq bias: {np.mean(f_true[mask] - freq_3rd[mask]):.3f} Hz")
このグラフが本記事の核心結果です。左図では二次ループと三次ループの追跡が一見似たように見えますが、右図の周波数誤差を比較すれば違いは歴然です。二次ループ（青）は約 $-10\,\mathrm{Hz}$ 前後の一定のバイアスを残し続けるのに対し、三次ループ（緑）は過渡を抜けた後ほぼゼロ平均で揺らいでいます。これが「三次ループはパラボラ入力（一定のドップラー変化率）に対して定常誤差ゼロ」という解析結果の Python による実証です。
ただし三次ループには引換えがあります。三つ目の積分器がループの位相余裕を減らすため、設計が悪いと振動的な挙動を示します。実装時は経験的に検証された係数セット（Stephens & Thomas 1995、Kaplan & Hegarty の GNSS 教科書など）から出発するのが安全です。
Costasループの実装と $\pi$ 曖昧性
最後に、BPSK 信号に対する Costas ループの動作を確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import freqz
# BPSK 変調信号
fs = 10_000.0
np.random.seed(42)
N_bits = 500
Nsym = 100 # 1ビットあたりのサンプル
N_total = N_bits * Nsym
t2 = np.arange(N_total) / fs
data = 2*np.random.randint(0, 2, N_bits) - 1
bpsk = np.repeat(data, Nsym).astype(float)
# ドップラー
f0_b = 50.0; fdot_b = 10.0
phi_b = 2*np.pi*(f0_b*t2 + 0.5*fdot_b*t2**2)
A = 1.0; CN0 = 50.0
N0 = A**2 / (10**(CN0/10))
noise2 = np.sqrt(N0*fs/2)*(np.random.randn(N_total)+1j*np.random.randn(N_total))
rx_bpsk = A * bpsk * np.exp(1j * phi_b) + noise2
def costas_loop(rx, fs, BL, zeta=0.707):
wn = BL * 8 * zeta / (4 * zeta**2 + 1); Ts = 1/fs
Kp = (8*zeta*wn*Ts)/(1+2*zeta*wn*Ts+(wn*Ts)**2)
Ki = (4*(wn*Ts)**2)/(1+2*zeta*wn*Ts+(wn*Ts)**2)
N = len(rx); I_out = np.zeros(N); Q_out = np.zeros(N)
err = np.zeros(N); integ = 0.0; phase = 0.0
for n in range(N):
prod = rx[n] * np.exp(-1j * phase)
I_n = prod.real; Q_n = prod.imag
amp = I_n**2 + Q_n**2 + 1e-9
e = I_n * Q_n / amp # 正規化Costas検出
I_out[n] = I_n; Q_out[n] = Q_n; err[n] = e
integ += Ki * e
phase += Kp * e + integ
return I_out, Q_out, err
I_c, Q_c, err_c = costas_loop(rx_bpsk, fs, BL=30.0)
fig, axes = plt.subplots(2, 2, figsize=(13, 8))
mask2 = t2 > 1.0
axes[0, 0].scatter(I_c[mask2][::20], Q_c[mask2][::20], s=2, alpha=0.4)
axes[0, 0].axhline(0, color='k', linewidth=0.5)
axes[0, 0].axvline(0, color='k', linewidth=0.5)
axes[0, 0].set_xlabel('I'); axes[0, 0].set_ylabel('Q')
axes[0, 0].set_title('BPSK Constellation after Costas Lock')
axes[0, 0].set_aspect('equal'); axes[0, 0].grid(alpha=0.3)
axes[0, 1].plot(t2, I_c, 'b-', alpha=0.6, linewidth=0.5, label='I')
axes[0, 1].plot(t2, Q_c, 'r-', alpha=0.6, linewidth=0.5, label='Q')
axes[0, 1].plot(t2, bpsk, 'k--', alpha=0.5, label='True data')
axes[0, 1].set_xlabel('Time [s]'); axes[0, 1].set_ylabel('Amplitude')
axes[0, 1].set_title('I (data) and Q (~0)'); axes[0, 1].set_xlim([2.0, 2.3])
axes[0, 1].legend(); axes[0, 1].grid(alpha=0.3)
axes[1, 0].plot(t2, err_c, 'g-', alpha=0.6, linewidth=0.5)
axes[1, 0].set_xlabel('Time [s]'); axes[1, 0].set_ylabel('Costas Error')
axes[1, 0].set_title('Costas PD output (data-independent)')
axes[1, 0].grid(alpha=0.3)
# ループフィルタ周波数応答 (PI)
wn_ = 30.0*8*0.707/(4*0.707**2+1); Ts_ = 1/fs
Kp_ = (8*0.707*wn_*Ts_)/(1+2*0.707*wn_*Ts_+(wn_*Ts_)**2)
Ki_ = (4*(wn_*Ts_)**2)/(1+2*0.707*wn_*Ts_+(wn_*Ts_)**2)
b = [Kp_ + Ki_, -Kp_]; a = [1, -1]
w, h = freqz(b, a, worN=2048, fs=fs)
axes[1, 1].semilogx(w, 20*np.log10(np.abs(h)+1e-12))
axes[1, 1].set_xlabel('Frequency [Hz]'); axes[1, 1].set_ylabel('Mag [dB]')
axes[1, 1].set_title('Loop Filter (PI) Frequency Response')
axes[1, 1].grid(alpha=0.3, which='both'); axes[1, 1].set_xlim([0.1, fs/2])
plt.tight_layout()
plt.savefig('costas_loop_bpsk.png', dpi=130, bbox_inches='tight')
plt.show()
このグラフから Costas ループの動作が確認できます。左上のコンスタレーションは、ロック後の I/Q 平面上のシンボル分布で、$\pm 1$ の2点に明瞭に集中していることが分かります。これは位相同期が確立し、$Q$ 側の信号エネルギーがほぼゼロに追い込まれたことを意味します（$\pi$ 曖昧性で時々左右が反転する場合もありますが、本ランでは初期位相がゼロから出発したのでうまくロックしています）。
右上の時系列を見ると、$I$ が真のデータ波形（黒破線）と一致して符号反転を再現する一方、$Q$ は雑音レベルにとどまっています。左下のCostas誤差信号はデータの符号反転に影響されず、零点周辺で揺らぐだけです。これが $I \cdot Q$ という積によってデータ依存が消えた証拠です。右下のループフィルタの周波数応答は、低域で大きなゲインを持ち、高周波の雑音を抑制する PI フィルタの典型的な特性を示しています。
実装まで来ました。あとは衛星通信受信機の起動シーケンスで実際に使われる FLL→PLL 切替について見ておきましょう。
FLL→PLL切替シーケンス
引き込み（acquisition）と追跡（tracking）
PLL は線形化された範囲（位相誤差が $\pi/4$ 程度以下）で初めてきちんと動きます。受信機の起動直後は位相どころか周波数すら大きくずれているため、いきなり PLL を起動するとロックに数十秒かかったり、最悪ロックしません。
そこで実用の受信機は2段階の同期戦略を取ります。第一段階は粗い周波数同期（FLL）で、周波数を周波数誤差を直接検出するループで合わせます。位相は気にしません。広い周波数範囲で素早く引き込めるのが利点です。第二段階は精密な位相同期（PLL）で、FLL である程度周波数が合ったら PLL に切り替え、位相まで追い込みます。
FLL の位相検出器
FLLは、連続する2サンプル間の位相変化量を使って周波数誤差を出します。具体的には、
$$ e_{\mathrm{FLL}}[n] = \frac{1}{T_s}\arg\bigl( z[n] \cdot z^*[n-1] \bigr) $$
の形で、ここで $z[n]$ はミキシング後の複素サンプル、$T_s$ はサンプリング周期です。この検出器は位相のオフセットに対しては鈍感で、純粋に「周波数のずれ」だけを抽出します。連続するサンプルの位相差は、ちょうど周波数 $\Delta f$ だけ離れた成分なら $2\pi \Delta f T_s$ になるため、これを $T_s$ で割れば周波数誤差 $\Delta f$ そのものになります。
切替判定の戦略
FLL から PLL への切替判定は、典型的には以下のいずれか（または組み合わせ）で行います。第一に、周波数誤差の絶対値が閾値（例: $5\,\mathrm{Hz}$）以下になることです。第二に、$I/Q$ の比から推定した位相誤差の分散が小さくなることです。第三に、一定時間経過に加えてロック検出器（lock detector）の確信度が十分高いことです。
切替の瞬間にも工夫があり、FLL で確立した周波数推定値を PLL の積分器初期値として引き継ぐことで、PLL の引き込みを高速化します。これを「ウォーム・スタート」と呼びます。逆に、衛星の電波が一時的に途切れて（fade）位相ロックを失った場合は、PLL を一旦停止して FLL に戻り、再引き込みを試みる「リコンフィギュレーション」が走ります。
簡易シミュレーション
最後に、FLL→PLL 切替の概念を Python で確認しましょう。
import numpy as np
import matplotlib.pyplot as plt
fs = 10_000.0; T = 3.0; N = int(fs*T); t = np.arange(N)/fs
# 初期に大きな周波数誤差 +500Hz
f0 = 500.0
phi_true = 2 * np.pi * f0 * t
np.random.seed(7)
noise = 0.05 * (np.random.randn(N) + 1j*np.random.randn(N))
rx = np.exp(1j * phi_true) + noise
# ハイブリッド: 最初は FLL, その後 PLL
def hybrid(rx, fs, switch_t, BL_fll, BL_pll):
Ts = 1/fs
zeta = 0.707
wn = BL_pll * 8 * zeta / (4 * zeta**2 + 1)
Kp = (8*zeta*wn*Ts)/(1+2*zeta*wn*Ts+(wn*Ts)**2)
Ki = (4*(wn*Ts)**2)/(1+2*zeta*wn*Ts+(wn*Ts)**2)
K_fll = 2 * BL_fll
N = len(rx); freq = 0.0; phase = 0.0; integ = 0.0
freq_hist = np.zeros(N); prev = 1+0j
switch_n = int(switch_t * fs)
for n in range(N):
cur = rx[n] * np.exp(-1j * phase)
if n < switch_n:
# FLL: 連続サンプルの位相差を周波数誤差とする
if n > 0:
d = cur * np.conj(prev)
e = np.arctan2(d.imag, d.real) / Ts
freq += K_fll * Ts * e
phase += 2*np.pi*freq*Ts
else:
# PLL (ウォームスタート: FLL の周波数推定を積分器に引き継ぐ)
if n == switch_n:
integ = 2 * np.pi * freq * Ts
e = np.arctan2(cur.imag, cur.real)
integ += Ki * e
upd = Kp * e + integ
phase += upd
freq = upd * fs / (2*np.pi)
freq_hist[n] = freq
prev = cur
return freq_hist
# FLL のみ (比較用)
def fll_only(rx, fs, BL):
Ts = 1/fs
K = 2 * BL
N = len(rx); freq = 0.0; phase = 0.0
freq_hist = np.zeros(N); prev = 1+0j
for n in range(N):
cur = rx[n] * np.exp(-1j * phase)
if n > 0:
d = cur * np.conj(prev)
e = np.arctan2(d.imag, d.real) / Ts
freq += K * Ts * e
phase += 2*np.pi*freq*Ts
freq_hist[n] = freq
prev = cur
return freq_hist
freq_fll = fll_only(rx, fs, BL=50.0)
freq_hybrid = hybrid(rx, fs, switch_t=1.0, BL_fll=50.0, BL_pll=20.0)
fig, ax = plt.subplots(figsize=(11, 4.5))
ax.axhline(f0, color='k', linewidth=2, linestyle='--', label='True freq (500Hz)')
ax.plot(t, freq_fll, 'b-', alpha=0.7, label='FLL only ($B_L$=50Hz)')
ax.plot(t, freq_hybrid, 'g-', alpha=0.8, label='FLL$\\to$PLL switch at 1s')
ax.axvline(1.0, color='r', linestyle=':', label='Switch time')
ax.set_xlabel('Time [s]'); ax.set_ylabel('Frequency Estimate [Hz]')
ax.set_title('FLL$\\to$PLL Hybrid Acquisition')
ax.legend(); ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('fll_pll_handoff.png', dpi=130, bbox_inches='tight')
plt.show()
このグラフから FLL→PLL の切替戦略の効果が見て取れます。FLLのみ（青）は周波数誤差を素早く $500\,\mathrm{Hz}$ 近傍まで詰めますが、その後も雑音由来のジッタが大きく残ります。一方ハイブリッド方式（緑）は最初の $1$ 秒で FLL により大まかに周波数を引き込み、$1$ 秒の時点で PLL に切り替えると、$B_L = 20\,\mathrm{Hz}$ の狭帯域 PLL が引き継いで周波数推定がぐっと滑らかになります。切替時に FLL の推定値を PLL 積分器の初期値として渡しているため、切替時の不連続が小さく抑えられています。
実用受信機では、これに加えてさらに「ロック判定」「アンロック検出時のFLLへのフォールバック」「ループ帯域のスケジュール（最初は広く、その後狭く）」といったロジックが組み合わされ、ステートマシンとして実装されます。
ここまで衛星通信受信機で実際に使われる搬送波同期の全景を見てきました。最後にまとめましょう。
まとめ
本記事では、衛星通信受信機の心臓部である搬送波同期を、PLL と Costas ループの理論・解析・実装の三方向から詳しく解説しました。
- 3要素の構成: 位相検出器（PD）、ループフィルタ（LF）、数値制御発振器（NCO）の3要素で PLL は構築される。NCO は積分器として機能し、ループ全体の積分器の数がループ次数を決める
- ループ次数と定常誤差: 一次は周波数オフセットを追えず、二次はオフセットは追えるが変化率に有限誤差を残す。三次（タイプIII）ループは一定の周波数変化率まで定常誤差ゼロで追跡できる
- ループ帯域とジッタ: 位相ジッタ分散は $\sigma_\phi^2 = B_L / (C/N_0)$ で与えられ、$B_L$ を広げるとダイナミクス追跡力は上がるがジッタは増える。LEO通過のように $\dot f$ が大きい局面では適応的に $B_L$ を切り替える
- Costasループ: BPSK 信号のデータ反転に対し $I \cdot Q$ という形で位相誤差を作ることでデータ依存を消す。$\pi$ 位相曖昧性は差動符号化やフレーム同期で吸収する
- FLL→PLL切替: 受信機起動直後はFLLで広範囲の周波数引き込みを行い、誤差が小さくなったらPLLに切り替えて位相まで追い込む。FLLの推定値を PLL の積分器初期値として引き継ぐ「ウォームスタート」で切替を滑らかにする
これらの設計概念は衛星通信端末だけでなく、GNSS受信機、深宇宙テレメトリ受信、地上のクロック同期、レーダー信号処理など、無線受信機のほぼ全ての分野で共通の基礎を成します。本記事のPython実装をベースに、自分で雑音電力や $\dot f$ を変えて挙動を観察すると、設計トレードオフへの直感が一段深まります。
次のステップとして、以下の記事も参考にしてください。