衛星から送られてくるディジタル信号を地上で復調するとき、受信機はまず「いま届いている搬送波の位相がどこにあるのか」を知らなければなりません。ところが BPSK や QPSK といった位相変調では、送信側で電力を効率よく使うために搬送波そのものの成分が消し去られています。スペクトルを見ても、本来あるべき搬送波周波数のところにスパイク（純粋な正弦波の成分）が立っていないのです。これでは、搬送波周波数に同調しようとする普通の位相同期ループ（PLL）はロックする相手を失ってしまいます。それでもなお、受信機はその「見えない搬送波」の位相を寸分の狂いもなく復元しなければ、データを正しく読み取れません。
この一見すると不可能に思える問題を、1956 年に John P. Costas が驚くほど巧妙な方法で解決しました。それが本記事で扱う Costasループ（Costas loop） です。Costasループを理解すると、次のような場面で何が起きているのかが手に取るように見えてきます。
- 衛星通信のモデム: DVB-S2 や CCSDS の地上局受信機は、抑圧搬送波の BPSK/QPSK/APSK 信号からキャリアを再生するためにこの種のループを内蔵しています
- GNSS 受信機: GPS の航法データは BPSK で変調されており、受信機内のキャリアトラッキングループは Costasループ（あるいはその一般化）として実装されます
- ソフトウェア無線（SDR）: GNU Radio をはじめとする SDR フレームワークの位相同期ブロックは、まさに Costasループのアルゴリズムをそのままコードにしたものです
本記事の内容
- なぜ抑圧搬送波信号には普通の PLL がロックできないのか
- Costasループの I・Q 分岐と位相検出器が誤差信号 $\sin 2\Delta\phi$ を生む仕組みの導出
- ループフィルタと VCO（数値制御発振器）の役割と次数の選び方
- 位相の $\pi$ 曖昧性がなぜ生じ、差動符号化でどう解決するか
- QPSK 用の四相位相検出器への拡張
- Python による BPSK Costasループの実装、位相誤差の収束、周波数オフセット下の引き込み、コンスタレーションの収束の可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が格段に深まります。
抑圧搬送波とは — なぜ普通のPLLが効かないのか
搬送波が「抑圧される」とはどういうことか
AM ラジオを思い浮かべてください。音声信号 $m(t)$ を搬送波 $\cos\omega_c t$ に乗せるとき、古典的な AM では $s(t) = [A + m(t)]\cos\omega_c t$ のように、信号にバイアス $A$ を足してから掛けます。このバイアスがあるおかげで、$m(t)$ がゼロのときでも $A\cos\omega_c t$ という純粋な搬送波成分が残り続けます。スペクトル上では搬送波周波数 $f_c$ にくっきりとしたスパイクが立ち、受信機はこのスパイクを目印にして同調できます。
しかし、このバイアス $A$ を送り続けるのは電力の無駄です。情報を何も運んでいない純粋な搬送波に、送信電力のかなりの割合が費やされてしまうからです。そこでバイアスを取り去り、$s(t) = m(t)\cos\omega_c t$ として送るのが 両側波帯抑圧搬送波（DSB-SC: Double Sideband Suppressed Carrier） です。ここで「抑圧」とは、搬送波周波数のスパイクが消えていることを指します。BPSK はこの DSB-SC の特別な場合で、$m(t)$ が $+1$ か $-1$ のどちらかを取るディジタル信号になっています。
BPSK信号の数式とスペクトル
BPSK 信号を式で書くと、ビットに応じて符号 $d(t) \in \{+1, -1\}$ が変わる次の形になります。
$$ s(t) = d(t) A \cos(\omega_c t + \theta) $$
ここで $A$ は振幅、$\theta$ は搬送波の（受信機にとっては未知の）位相です。$d(t) = -1$ のとき $-\cos = \cos(\cdot + \pi)$ ですから、BPSK は「位相が $0$ か $\pi$ かでビットを表す変調」とも言えます。
このとき搬送波成分が消えていることを確認しましょう。受信機が搬送波を再生するには、$s(t)$ の中に $\cos\omega_c t$ という確定的な（データに依存しない）成分があってほしいのですが、$d(t)$ が $+1$ と $-1$ を等確率でランダムに取るなら、長時間平均すると、
$$ \langle s(t) \rangle = A\cos(\omega_c t + \theta)\,\langle d(t) \rangle = A\cos(\omega_c t + \theta)\cdot 0 = 0 $$
となり、確定的な搬送波成分の期待値はゼロです。データのランダムな符号反転によって搬送波が「平均的に打ち消されている」のです。これが抑圧搬送波の本質です。
普通のPLLが失敗する理由
通常の PLL は、入力信号の中にある正弦波成分に対して、内部の発振器（VCO）の位相を合わせ込みます。位相検出器が「入力と VCO の位相差」に比例した誤差信号を作り、それをフィードバックして VCO を引き込むわけです。ところが、ロックすべき相手である搬送波成分が平均ゼロで存在しないのですから、PLL は誤差信号を作れません。
もう少し具体的に見てみましょう。普通の PLL の位相検出器は、入力 $s(t)$ と VCO 出力 $\sin(\omega_c t + \hat\theta)$ を掛けてローパスフィルタを通します。掛け算すると、
$$ d(t) A \cos(\omega_c t + \theta)\sin(\omega_c t + \hat\theta) = \frac{d(t)A}{2}\big[\sin(\hat\theta – \theta) + \sin(2\omega_c t + \theta + \hat\theta)\big] $$
倍角項（$2\omega_c$ の成分）をローパスで落とすと、誤差信号は $\frac{d(t)A}{2}\sin(\hat\theta – \theta)$ になります。問題はここに $d(t)$ が残っていることです。データビットが反転するたびに誤差信号の符号も反転してしまい、フィードバックの向きが定まりません。VCO は引き込もうとするたびに逆向きに蹴飛ばされ、決してロックできないのです。
つまり、本質的な障害は「誤差信号がデータ $d(t)$ に汚染されている」ことにあります。Costasループの天才的なアイデアは、データの影響を打ち消す形で誤差信号を作ることにありました。次節でその仕組みを丁寧に組み立てていきましょう。
Costasループの基本構成 — I・Q二つの腕でデータを消す
全体像をイメージでつかむ
Costasループの核心は「二人がかりで観測する」ことだと考えると腑に落ちます。受信信号を、VCO が作る互いに $90^\circ$ ずれた二つの参照波（同相 $\cos$ と直交 $\sin$）にそれぞれ掛けて、二本の腕（I 腕と Q 腕）を作ります。片方の腕にはデータと位相誤差の両方が、もう片方の腕にも同じデータと位相誤差が現れます。この二本を掛け合わせると、データの符号 $d(t)$ は二乗されて $d(t)^2 = 1$ になり、きれいに消えてくれるのです。残るのは位相誤差だけ。これが Costas のマジックの種明かしです。
ブロック図としては次の要素から成ります。
- 入力信号を二分配する
- I 腕: 入力に VCO の $\cos$ を掛け、ローパスフィルタ（LPF）を通す
- Q 腕: 入力に VCO の $\sin$ を掛け、ローパスフィルタを通す
- 位相検出器: I 腕と Q 腕の出力を掛け算する
- ループフィルタ: 掛け算結果を平滑化して制御信号にする
- VCO: 制御信号に応じて自分の位相・周波数を調整する
I腕とQ腕の出力を導出する
受信信号を、これまで通り次のように置きます（時間引数は省略します）。
$$ s = d \cdot A\cos(\omega_c t + \theta) $$
VCO が生成する二つの参照波は、推定位相 $\hat\theta$ を使って、
$$ v_I = 2\cos(\omega_c t + \hat\theta), \qquad v_Q = -2\sin(\omega_c t + \hat\theta) $$
とします（振幅 $2$ は後で $1/2$ が出るのを相殺するための便宜的な係数で、本質ではありません）。まず I 腕で $s$ と $v_I$ を掛けます。積和の公式 $2\cos\alpha\cos\beta = \cos(\alpha-\beta) + \cos(\alpha+\beta)$ を使うと、
$$ s \cdot v_I = d A \cdot 2\cos(\omega_c t + \theta)\cos(\omega_c t + \hat\theta) = dA\big[\cos(\theta – \hat\theta) + \cos(2\omega_c t + \theta + \hat\theta)\big] $$
第 2 項は $2\omega_c$ の高周波なのでローパスフィルタで除去されます。$2\omega_c$ 成分が落ちることを使うと、I 腕の出力は次のようになります。
$$ I = dA\cos(\theta – \hat\theta) $$
ここで位相誤差を $\Delta\phi = \theta – \hat\theta$ と定義すると、$I = dA\cos\Delta\phi$ です。
同じことを Q 腕でも行います。$s$ と $v_Q = -2\sin(\omega_c t + \hat\theta)$ を掛け、積和の公式 $2\cos\alpha\sin\beta = \sin(\alpha+\beta) – \sin(\alpha – \beta)$ を $\alpha = \omega_c t + \theta,\ \beta = \omega_c t + \hat\theta$ に適用します。
$$ s \cdot v_Q = -dA\big[\sin(2\omega_c t + \theta + \hat\theta) – \sin(\theta – \hat\theta)\big] $$
ここでも $2\omega_c$ の項をローパスフィルタで落とすと、
$$ Q = dA\sin(\theta – \hat\theta) = dA\sin\Delta\phi $$
が残ります。I 腕には $\cos\Delta\phi$ が、Q 腕には $\sin\Delta\phi$ が現れる、というのが Costasループの第一の重要な結果です。直感的に言えば、I 腕は「VCO が正しく合っていればデータがそのまま出てくる腕」、Q 腕は「位相がずれているときだけ信号が漏れ出す腕」です。位相が完全に合えば $\Delta\phi = 0$ で $Q = 0$ となり、Q 腕の信号が消えます。
位相検出器 — I·Qが誤差信号を作る
さて、ここからが Costas の真骨頂です。I と Q を単純に掛け算します。
$$ e = I \cdot Q = (dA\cos\Delta\phi)(dA\sin\Delta\phi) = d^2 A^2 \cos\Delta\phi\sin\Delta\phi $$
BPSK では $d \in \{+1, -1\}$ なので、$d^2 = 1$ です。データの符号が二乗によって完全に消えました。これがデータ変調を打ち消すトリックの正体です。残った積に倍角公式 $2\sin\Delta\phi\cos\Delta\phi = \sin 2\Delta\phi$ を使うと、
$$ e = A^2\cos\Delta\phi\sin\Delta\phi = \frac{A^2}{2}\sin 2\Delta\phi $$
これが Costasループの 位相検出器（PD）の出力 であり、ループの誤差信号です。重要な性質を二つ確認しておきましょう。
第一に、誤差信号 $e$ は データ $d$ に依存しません。ビットがどう反転しようと $d^2 = 1$ なので、誤差信号の符号が暴れることはありません。普通の PLL が抱えていた「誤差信号がデータに汚染される」問題が、ここで根本的に解消されています。
第二に、$\Delta\phi$ が小さいとき $\sin 2\Delta\phi \approx 2\Delta\phi$ なので、誤差信号は位相誤差にほぼ比例します。これは普通の PLL の位相検出器特性 $\sin\Delta\phi$ とよく似た振る舞いで、フィードバック制御の対象として理想的です。$\Delta\phi = 0$ で $e = 0$（平衡点）、$\Delta\phi$ がわずかに正なら $e > 0$ で VCO を正しい向きに押し戻す、という負帰還が成立します。
ここまでで「データを消しつつ位相誤差を取り出す」という Costasループの心臓部が完成しました。次は、この誤差信号を使ってどのように VCO を制御し、ループとして閉じるのかを見ていきます。
ループフィルタとVCO — ループを閉じる
誤差信号からVCO制御へ
位相検出器が $e = \frac{A^2}{2}\sin 2\Delta\phi$ を吐き出しても、それだけでは VCO を動かせません。誤差信号には雑音や、ローパスで取り切れなかった成分が残っているため、これを平滑化し、かつループの動特性（応答の速さ・安定性）を整形する必要があります。その役割を担うのが ループフィルタ です。
ループフィルタの出力を $u(t)$ とすると、VCO はこの $u(t)$ に比例した瞬時周波数で位相を進めます。離散時間（サンプル周期 $T_s$、サンプル番号 $k$）で書くと、VCO の位相更新は次の積分になります。
$$ \hat\theta_{k+1} = \hat\theta_k + \omega_0 T_s + u_k $$
ここで $\omega_0 T_s$ は VCO の中心周波数に対応する位相進み、$u_k$ がループフィルタからの制御による補正です。VCO は本質的に「制御入力を積分して位相を作る」素子であり、これがループの中に必ず一つ積分器（極が原点にある要素）を持つことを意味します。
比例・積分ループフィルタ（2次ループ）
最もよく使われるのは、比例項と積分項を持つ PI 型ループフィルタ です。連続時間では伝達関数 $F(s) = K_p + K_i/s$、離散時間では次の更新式で表せます。
$$ u_k = K_p\, e_k + K_i \sum_{j=0}^{k} e_j $$
右辺第 1 項（比例項 $K_p e_k$）は、いまの瞬間の位相誤差に即座に反応してループを速く引き込みます。第 2 項（積分項）は、過去の誤差を累積することで、定常的な周波数オフセットを吸収する役割を果たします。
なぜ積分項が周波数オフセットの吸収に必要なのか、直感を補っておきましょう。送信側と受信側の発振器周波数がわずかにずれている（あるいはドップラーで $\omega_c$ がずれている）と、位相誤差 $\Delta\phi$ は時間とともに一定の傾きで増えていきます。これを止めるには、VCO の周波数自体を一定量だけ恒久的にシフトさせ続ける必要があります。比例項だけ（1 次ループ）では、誤差がゼロでなければ補正が出ないため、必ず一定の残留位相誤差（あるいは周波数追従できない状態）が生じます。積分項があると、誤差がゼロになるまで補正を溜め込み続けるので、最終的に周波数オフセットを完全に打ち消せるのです。
VCO の積分とループフィルタの積分を合わせると、ループ全体は 2 つの積分器を持つ 2 次ループ になります。2 次ループは、一定の周波数オフセットに対して定常位相誤差ゼロで追従できる、という極めて有用な性質を持ちます。
自然角周波数と減衰係数
2 次ループの設計では、制御工学でおなじみの 自然角周波数 $\omega_n$ と 減衰係数 $\zeta$ という二つのパラメータで挙動を特徴づけます。連続時間の閉ループ伝達関数を 2 次系の標準形に合わせると、ループゲイン $K$（位相検出器ゲイン $\times$ VCO ゲイン）とフィルタ係数の間に次の関係が成り立ちます。
$$ \omega_n = \sqrt{K K_i}, \qquad \zeta = \frac{1}{2}\sqrt{\frac{K K_p^2}{K_i}} = \frac{K K_p}{2\omega_n} $$
実装では、所望の $\omega_n$（引き込みの速さ）と $\zeta$（オーバーシュートの少なさ、通常 $\zeta = 0.707$ がよく選ばれる）を先に決め、そこから係数 $K_p, K_i$ を逆算します。離散時間のループ帯域 $B_L$ とサンプル周波数 $f_s$ を使った実用的な係数設計式は、規格化帯域 $\theta = B_L T_s$ を用いて、
$$ K_p = \frac{4\zeta\theta}{1 + 2\zeta\theta + \theta^2}, \qquad K_i = \frac{4\theta^2}{1 + 2\zeta\theta + \theta^2} $$
のように書けます（VCO とPDのゲインを 1 に正規化した場合）。ループ帯域 $B_L$ を広く取れば引き込みは速くなりますが、雑音を多く通すので位相ジッタが増えます。逆に狭くすれば雑音には強くなりますが、引き込みが遅く、速い周波数変動（高いドップラーレート）に追従できなくなります。このトレードオフが受信機設計の腕の見せどころです。
ここまでで Costasループの動的な振る舞いを決める部品が揃いました。しかし、ループがロックしても実は一つ厄介な問題が残ります。それが「位相の曖昧性」です。次節で詳しく見ましょう。
位相の π 曖昧性とその解決
なぜ π ずれてもロックしてしまうのか
誤差信号 $e = \frac{A^2}{2}\sin 2\Delta\phi$ をもう一度眺めてください。$\sin 2\Delta\phi$ がゼロになるのは $\Delta\phi = 0$ のときだけではありません。$2\Delta\phi = 0, \pi, 2\pi, \dots$、つまり $\Delta\phi = 0, \pi/2, \pi, 3\pi/2, \dots$ でもゼロになります。このうち安定平衡点（誤差信号が引き戻す向きを持つ点）は、$\sin 2\Delta\phi$ の傾きが正になる $\Delta\phi = 0$ と $\Delta\phi = \pi$ の二つです。
つまり、Costasループは $\Delta\phi = 0$ にロックすることもあれば、$\Delta\phi = \pi$ にロックすることもあります。$\Delta\phi = \pi$ にロックすると何が起こるでしょうか。I 腕の出力は $I = dA\cos\pi = -dA$ となり、全データビットの符号が反転して復調されてしまうのです。$+1$ が $-1$ に、$-1$ が $+1$ に化けます。これが $\pi$ 位相曖昧性（phase ambiguity） です。
これは Costasループの宿命とも言える性質です。BPSK の信号 $d\cos(\omega_c t + \theta)$ は、$d \to -d$ かつ $\theta \to \theta + \pi$ としても全く同じ波形になります。受信機には、いま見ている搬送波の位相が「本来の $\theta$」なのか「$\theta + \pi$」なのかを、波形だけからは区別する手段が原理的に存在しないのです。Costasループがデータを二乗して消してしまった以上、その代償として絶対位相の情報も失われます。
差動符号化による解決
この曖昧性を解決する標準的な方法が 差動符号化（differential encoding） です。アイデアは単純で、「絶対的なビットの値」ではなく「ビットが前回から変化したかどうか」に情報を載せるのです。
送信側では、送りたい情報ビット列 $\{a_k\}$ から、送信シンボル列 $\{b_k\}$ を次の漸化式で作ります。
$$ b_k = a_k \oplus b_{k-1} $$
ここで $\oplus$ は排他的論理和（XOR）です。たとえば $a_k = 0$ なら $b_k = b_{k-1}$（前と同じ）、$a_k = 1$ なら $b_k = \overline{b_{k-1}}$（前と反転）となります。受信側では逆に、復調したシンボル $\{\hat b_k\}$ から、
$$ \hat a_k = \hat b_k \oplus \hat b_{k-1} $$
として情報ビットを復元します。なぜこれで曖昧性が解決するのか、考えてみましょう。仮に Costasループが $\pi$ ずれてロックし、全シンボルが反転して $\hat b_k \to \overline{\hat b_k}$ になったとします。すると差動復号の出力は、
$$ \overline{\hat b_k} \oplus \overline{\hat b_{k-1}} = \hat b_k \oplus \hat b_{k-1} = \hat a_k $$
となります。連続する二つのシンボルが両方とも同じように反転するので、その XOR を取ると反転が打ち消し合い、元の情報ビットがそのまま復元されるのです。差動符号化の本質は「絶対位相ではなく位相の差分に情報を載せる」ことであり、これにより絶対位相の曖昧性は無害になります。
代償として、1 ビットの誤りが隣接する 2 つの復号ビットに波及するため、ビット誤り率はおおよそ 2 倍に悪化します。しかし、$\pi$ 曖昧性で全データが反転して使い物にならなくなることを思えば、これは十分許容できるコストです。実際の衛星通信や GNSS では、差動符号化に加えて、既知のパターン（ユニークワードやプリアンブル）を周期的に挿入して絶対位相を解決する方法も併用されます。
ここまでは BPSK（二相）に話を絞ってきました。実際の高効率な通信では、一度に 2 ビットを送る QPSK（四相）が広く使われます。次節では、Costasループを QPSK へ拡張する方法を見ていきます。
QPSK用Costasループ — 四相位相検出器
QPSKでは I も Q も両方データを運ぶ
QPSK では、I 腕と Q 腕の両方に独立したデータが乗ります。送信信号は二つのビット列 $d_I, d_Q \in \{+1, -1\}$ を使って、
$$ s(t) = A\big[d_I\cos(\omega_c t + \theta) – d_Q\sin(\omega_c t + \theta)\big] $$
と書けます。コンスタレーション上では、四つの点 $(\pm 1, \pm 1)$ が対角線上に配置された姿になります。受信機が VCO で復調すると、位相誤差 $\Delta\phi$ がある場合、I 腕・Q 腕の出力は回転行列を掛けた形になります。導出すると、
$$ I = A\big(d_I\cos\Delta\phi + d_Q\sin\Delta\phi\big), \qquad Q = A\big(d_Q\cos\Delta\phi – d_I\sin\Delta\phi\big) $$
となります。これはコンスタレーションが角度 $\Delta\phi$ だけ回転して見える、という幾何学的な事実をそのまま表しています。BPSK の Costasループのように単純に $I \cdot Q$ を掛けただけでは、$d_I d_Q$ が残ってしまい誤差信号が汚れます。QPSK では別の工夫が要ります。
四相位相検出器の誤差信号
QPSK 用 Costasループで広く使われる位相検出器は、次の形の誤差信号を作ります。
$$ e = Q \cdot \mathrm{sgn}(I) – I \cdot \mathrm{sgn}(Q) $$
ここで $\mathrm{sgn}(\cdot)$ は符号関数です。$\mathrm{sgn}(I)$ は I 腕で判定したビットの推定値（硬判定）に相当し、$\mathrm{sgn}(Q)$ は Q 腕の判定値です。なぜこれが位相誤差を表すのか、小さな $\Delta\phi$ で確かめましょう。$\Delta\phi$ が小さければ判定は正しく、$\mathrm{sgn}(I) \approx d_I$、$\mathrm{sgn}(Q) \approx d_Q$ です。これらを代入すると、
$$ e \approx Q d_I – I d_Q = A(d_Q\cos\Delta\phi – d_I\sin\Delta\phi)d_I – A(d_I\cos\Delta\phi + d_Q\sin\Delta\phi)d_Q $$
右辺を展開します。$d_I^2 = d_Q^2 = 1$ を使うと $d_I d_Q\cos\Delta\phi$ の項は相殺し、$-\sin\Delta\phi$ の項だけが二回足し合わされて残ります。
$$ e \approx A\big(d_I d_Q\cos\Delta\phi – \sin\Delta\phi – \cos\Delta\phi\, d_I d_Q – \sin\Delta\phi\big) = -2A\sin\Delta\phi $$
このように、データ $d_I, d_Q$ がきれいに消え、$-2A\sin\Delta\phi$ という位相誤差に比例した誤差信号が得られます。BPSK のときと同じく、データを打ち消して位相誤差だけを抽出する構造になっているわけです。
この「判定指向（decision-directed）」型の位相検出器は、判定が正しい限り良好に働きますが、SNR が低くて判定が頻繁に誤るとループが乱れます。そのため実装では、低 SNR でも安定する別の位相検出器（たとえば $I^3 Q – IQ^3$ の形など）が使われることもあります。いずれにせよ、QPSK 版では曖昧性が $\pi/2$ ごと（90 度の四通り）に広がるため、差動符号化も二相ビットそれぞれに対して行う必要があります。
ここまで理論を組み立ててきました。あとは実際に動かして、誤差信号がどう収束し、コンスタレーションがどう整列していくかを目で確かめましょう。次節は Python による実装です。
Pythonでの実装
BPSK信号の生成
まず、未知の位相と周波数オフセットを持つ BPSK 信号を生成します。受信機側はこの位相・周波数を知らないという設定で、Costasループに復元させます。
import numpy as np
import matplotlib.pyplot as plt
# サンプリングと信号のパラメータ
fs = 8000.0 # サンプリング周波数 [Hz]
N = 4000 # サンプル数
sps = 8 # 1シンボルあたりのサンプル数
n_sym = N // sps # シンボル数
# ランダムなBPSKシンボル列 (+1 / -1)
rng = np.random.default_rng(0)
bits = rng.integers(0, 2, n_sym)
symbols = 2 * bits - 1 # 0/1 -> -1/+1
d = np.repeat(symbols, sps) # 各シンボルをsps回繰り返してベースバンド波形に
# 受信機が知らない真の搬送波: 位相オフセットと周波数オフセット
true_phase = 1.2 # 初期位相オフセット [rad]
foff = 20.0 # 周波数オフセット [Hz]
t = np.arange(N) / fs
carrier_phase = 2 * np.pi * foff * t + true_phase
# ベースバンドBPSK信号(複素表現)。搬送波位相で回転させる
rx = d * np.exp(1j * carrier_phase)
ここでは複素ベースバンド表現を使っています。実際の RF 信号を直交ダウンコンバートした後の信号がこの複素信号 rx に相当し、d がデータ（$\pm 1$）、carrier_phase が受信機にとって未知の残留位相・周波数です。次にこの rx を Costasループに通して位相を復元します。
Costasループの実装
PI 型ループフィルタを持つ 2 次の Costasループを実装します。複素表現では、VCO の参照波を掛けて回転を打ち消すと、その実部が I 腕、虚部が Q 腕に対応します。
import numpy as np
def costas_loop_bpsk(rx, Kp, Ki):
"""BPSK用Costasループ。複素ベースバンド信号rxを位相補正する。"""
N = len(rx)
phase = 0.0 # VCOの推定位相
freq = 0.0 # VCOの推定周波数(積分器の状態)
integ = 0.0 # 積分項の累積
out = np.zeros(N, dtype=complex) # 位相補正後の信号
err_log = np.zeros(N) # 誤差信号の記録
phase_log = np.zeros(N) # 推定位相の記録
for k in range(N):
# VCO参照波で回転を打ち消す(復調)
y = rx[k] * np.exp(-1j * phase)
out[k] = y
I, Q = y.real, y.imag
# 位相検出器: e = I * Q (BPSK Costas)
err = I * Q
err_log[k] = err
# PIループフィルタ
integ += Ki * err
freq = integ
phase += freq + Kp * err
phase_log[k] = phase
return out, err_log, phase_log
このループの肝は err = I * Q の一行です。前節で導出した $e = \frac{A^2}{2}\sin 2\Delta\phi$ を、複素信号の実部と虚部の積として計算しています。Kp * err が比例項で即応的な引き込みを、integ（Ki で累積）が積分項で周波数オフセットの吸収を担います。freq が周波数の推定値（積分器の状態）として保持され続ける点に注目してください。
ループ係数の設計と実行
ループ帯域と減衰係数から係数 $K_p, K_i$ を計算して実行します。
import numpy as np
# ループ係数の設計(規格化帯域から)
zeta = 0.707 # 減衰係数(臨界減衰に近い)
BL = 0.01 # 規格化ループ帯域 (B_L * Ts)
theta = BL
denom = 1 + 2 * zeta * theta + theta**2
Kp = (4 * zeta * theta) / denom
Ki = (4 * theta**2) / denom
print(f"Kp = {Kp:.5f}, Ki = {Ki:.6f}")
out, err_log, phase_log = costas_loop_bpsk(rx, Kp, Ki)
このコードはループ帯域 $B_L T_s = 0.01$、減衰係数 $\zeta = 0.707$ から係数を計算します。出力される Kp は約 0.028、Ki は約 0.0004 程度の小さな値になります。これらが小さいほどループは狭帯域で雑音に強くなりますが、引き込みは遅くなります。次に、ループがどのように位相誤差を収束させていくかを可視化します。
位相誤差の収束を可視化する
真の位相 carrier_phase と、Costasループが推定した位相 phase_log を比較し、その誤差がゼロに収束する様子を描きます。
import numpy as np
import matplotlib.pyplot as plt
# 真の位相と推定位相の差(πの整数倍の曖昧性を考慮してラップ)
true_phase_arr = carrier_phase
phase_error = true_phase_arr - phase_log
# BPSKのπ曖昧性: 誤差をπで畳んで表示
phase_error_wrapped = np.mod(phase_error + np.pi/2, np.pi) - np.pi/2
fig, ax = plt.subplots(2, 1, figsize=(10, 7))
ax[0].plot(err_log, color='tab:blue', lw=0.8)
ax[0].set_title('Costas Loop Error Signal e = I*Q')
ax[0].set_xlabel('Sample')
ax[0].set_ylabel('Error')
ax[0].grid(alpha=0.3)
ax[1].plot(phase_error_wrapped, color='tab:red', lw=1.0)
ax[1].axhline(0, color='k', ls='--', lw=0.8)
ax[1].set_title('Residual Phase Error (wrapped to [-pi/2, pi/2])')
ax[1].set_xlabel('Sample')
ax[1].set_ylabel('Phase error [rad]')
ax[1].grid(alpha=0.3)
plt.tight_layout()
plt.savefig('costas_convergence.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフからは二つのことが読み取れます。第一に、上段の誤差信号 $e = I\cdot Q$ は、初めのうちは大きく振れていますが、ループが引き込むにつれて振幅が小さくなり、平均的にゼロ付近へ収束していきます。データの符号反転によるスパイク的な揺らぎは残りますが、その「平均」がゼロに向かう点が重要です。第二に、下段の残留位相誤差は、初期位相オフセット（1.2 rad）と周波数オフセット（20 Hz）があるにもかかわらず、数百サンプルで $0$（または $\pi$）付近に収束しています。これは 2 次ループの積分項が周波数オフセットを吸収し、定常位相誤差をゼロにできていることの証拠です。
コンスタレーションの収束
最後に、位相補正前と補正後のコンスタレーション（IQ 平面上の信号点配置）を比較します。ループがロックすると、回転していた信号点が実軸上の $\pm 1$ にきれいに整列するはずです。
import numpy as np
import matplotlib.pyplot as plt
# シンボル中央のサンプルを抽出してコンスタレーションを作る
mid = sps // 2
idx = np.arange(mid, N, sps)
const_before = rx[idx] # 補正前(回転している)
const_after = out[idx] # 補正後
# 後半(ロック後)のみを取り出す
lock = idx > N // 2
fig, ax = plt.subplots(1, 2, figsize=(11, 5))
ax[0].scatter(const_before.real, const_before.imag, s=8,
color='tab:gray', alpha=0.6)
ax[0].set_title('Before: rotating constellation')
ax[0].set_xlabel('I'); ax[0].set_ylabel('Q')
ax[0].set_aspect('equal'); ax[0].grid(alpha=0.3)
ax[0].set_xlim(-1.6, 1.6); ax[0].set_ylim(-1.6, 1.6)
ax[1].scatter(const_after[lock[:len(const_after)]].real,
const_after[lock[:len(const_after)]].imag,
s=8, color='tab:green', alpha=0.6)
ax[1].set_title('After Costas lock: aligned to +/-1')
ax[1].set_xlabel('I'); ax[1].set_ylabel('Q')
ax[1].set_aspect('equal'); ax[1].grid(alpha=0.3)
ax[1].set_xlim(-1.6, 1.6); ax[1].set_ylim(-1.6, 1.6)
plt.tight_layout()
plt.savefig('costas_constellation.png', dpi=150, bbox_inches='tight')
plt.show()
左図（補正前）では、周波数オフセットのために信号点が原点まわりにぐるぐると回り、円環状にばらけています。位相がまったく定まっていない状態です。一方、右図（Costasロック後）では、信号点が実軸上の $+1$ と $-1$ の二点に集中しています。Q 軸（虚部）方向の広がりがほぼ消えているのは、$\Delta\phi \to 0$ で $Q = dA\sin\Delta\phi \to 0$ となったことの直接の現れです。これこそ Costasループが「見えない搬送波」の位相を復元できた証拠であり、ここから先は I 腕の符号を読むだけでビットを判定できます。
周波数オフセット下の引き込みを確かめる
周波数オフセットを変えながら、ループが推定する周波数が真の値に追従するかを確認します。これは 2 次ループの真価を見るための実験です。
import numpy as np
import matplotlib.pyplot as plt
def run_with_foff(foff, fs, N, sps, Kp, Ki, seed=1):
rng = np.random.default_rng(seed)
n_sym = N // sps
symbols = 2 * rng.integers(0, 2, n_sym) - 1
d = np.repeat(symbols, sps)
t = np.arange(N) / fs
rx = d * np.exp(1j * (2*np.pi*foff*t + 0.5))
# ループ内で周波数推定の履歴を取得
phase = 0.0; integ = 0.0
freq_log = np.zeros(N)
for k in range(N):
y = rx[k] * np.exp(-1j * phase)
err = y.real * y.imag
integ += Ki * err
phase += integ + Kp * err
freq_log[k] = integ
# 推定周波数 [Hz] に換算 (1サンプルあたりの位相増分 -> Hz)
return freq_log * fs / (2*np.pi)
plt.figure(figsize=(10, 6))
for foff in [10.0, 25.0, -15.0]:
fhat = run_with_foff(foff, fs, N, sps, Kp, Ki)
plt.plot(fhat, lw=1.2, label=f'true foff = {foff:.0f} Hz')
plt.axhline(foff, color='k', ls=':', lw=0.7)
plt.title('Estimated Frequency vs Sample (2nd-order Costas loop)')
plt.xlabel('Sample'); plt.ylabel('Estimated frequency [Hz]')
plt.legend(); plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('costas_freq_pullin.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、3 つの異なる周波数オフセット（+10 Hz、+25 Hz、−15 Hz）に対して、ループの周波数推定値がそれぞれ真の値（点線）へ収束していく様子が読み取れます。積分器に蓄えられた状態 integ が、まさに周波数オフセットを表す量に育っていくのです。オフセットが大きいほど引き込みに時間がかかること、そして負のオフセットにもきちんと追従できることが確認できます。もしループ帯域 $B_L$ を狭くしすぎると、大きなオフセットに対しては引き込みに失敗する（ロックレンジを超える）ことも、BL の値を変えて試すと体感できます。
これらの実験により、Costasループが抑圧搬送波の BPSK 信号から、位相と周波数の両方を雑音なしの理想条件で正確に復元できることが確かめられました。
まとめ
本記事では、抑圧搬送波信号から搬送波位相を取り戻す Costasループについて、原理から実装まで解説しました。要点を整理します。
- 抑圧搬送波（BPSK/DSB-SC）には搬送波成分が平均ゼロで存在せず、普通の PLL は誤差信号がデータ $d(t)$ に汚染されてロックできない
- Costasループは I 腕（$\cos$）と Q 腕（$\sin$）の二本に分け、$I = dA\cos\Delta\phi$、$Q = dA\sin\Delta\phi$ を得る
- 位相検出器 $e = I\cdot Q = \frac{A^2}{2}\sin 2\Delta\phi$ はデータが二乗されて $d^2 = 1$ となり消えるため、位相誤差だけを取り出せる
- PI 型ループフィルタ + VCO で 2 次ループを構成し、比例項で速い引き込み、積分項で周波数オフセットの吸収を実現する
- $\sin 2\Delta\phi$ の平衡点が $\Delta\phi = 0$ と $\pi$ の二つあるため $\pi$ 位相曖昧性が生じ、差動符号化で位相の差分に情報を載せることで無害化する
- QPSK では判定指向型の四相位相検出器 $e = Q\,\mathrm{sgn}(I) – I\,\mathrm{sgn}(Q)$ を用い、データを消して $-2A\sin\Delta\phi$ を得る
Costasループは、ディジタル通信受信機におけるキャリア同期の中核技術です。ここで身につけた「データを二乗・判定で消し、位相誤差だけを残す」という発想は、シンボルタイミング同期（Gardner 検出器など）やより高次の APSK のキャリア再生にも通じる普遍的なものです。受信機の信号処理チェーン全体を理解するための足がかりになるでしょう。
次のステップとして、以下の記事も参考にしてください。