GPSの信号を捜索する

はじめに

本記事はオープンソース半導体 Advent Calendar 2025の12/20の記事です。

オープンソースツール、オープンソースPDKを用いてLSI設計を行っているグループ「ISHI会」では、活動の一つとしてオープンハードウェアのGPS受信機をロケットに搭載することを目指して日々GPS受信機の設計を行っています。

この記事ではGPSの信号で自分の位置を知るまでの一番最初の部分である「GPSの信号を捜索する」までを書いてみたいと思います。せっかくなので、実際にGPSの信号を探すプロセスをPythonを用いて記述し、GPS衛星の信号を発見してみるまでステップ・バイ・ステップでやってみようと思います。

GPSの信号を使って自分の位置を知るまでの全てについてはISHI会で作成する同人誌に掲載される予定です。（自分へのプレッシャー）

What is GPS?

GPSはGloval Positioning Systemです。その名の通り、地球上で位置を知る（測位）ためのシステムです。人工衛星から受信した信号を使って位置を知ることができます。GPSはアメリカが運用しているシステムですGPS以外にも同様のシステムが運用されています。下記表は世界中で運用されている衛星測位システムです。

名前	運用国
GPS	アメリカ合衆国
GLONASS	ロシア連邦
BeiDou	中華人民共和国
Galileo	欧州連合
NavIC	インド共和国
QZSS	日本

NavICとQZSSは地域コンステレーションと呼ばれる、地域限定の測位システムです。このように世界中でGPSと同様のシステムが運用されています。このような人工衛星を利用した測位システムをGloval Navigation Satelite Systemの頭文字を取ってGNSSと言います。

ここからはなにか特定のシステムを話すときにはそれぞれのシステムの名称を、特定のシステムに限らず人工衛星を用いたシステムを話すときはGNSSと呼称します。

GPSの信号

GPSではいくつかの周波数帯で信号を送信しています。いくつかある中でスマートフォンやカーナビなど、一般の測位用途で使用する信号をL1C/Aといいます。C/AはCoarse Acquisitionのことで、（他の信号と比べて）精度の粗い測位に使用する信号です。このあとに詳細を紹介しますが、疑似ランダムコードのコード長が1023しかないので精度が少し劣ります。ただ、このコード長の短さのおかげてハードウェアが比較的かんたんになるので、安価なFPGAでも実装することができるようになります。

通信でも、放送でも、同時に情報を送るには異なる周波数を使うと思います。移動通信ではバンドと言う形でキャリアごとに異なる周波数が割り当てられ、その周波数を時分割で共有しています。テレビやラジオ放送ではチャンネルと言う形で異なる周波数が割り当てられています。

しかし、GPSでは全衛星が同じ周波数（L1C/Aでは1575.42MHz)で情報を放送しています。しかも、すべての衛星がバラバラの情報を送信しています。受信機ではどのようにこの信号を受信して、どの衛星からの信号なのか識別しているのでしょうか?

ここからは、GPSの信号がどのように送信されているか順を追ってみていきます。GPS衛星には次のような送信機が搭載されています。

一番左側の10.23MHzの基準周波数源となる発振器から内部の信号が全て作られています。この基準信号の発振器には原子時計が使われます。Rb、Cs、Hメーザーなどが使われます。GNSSの根幹は時間差を用いて距離を測定することなので、時間差を高精度に計算できるように高安定な発振器が必要になります。

それでは一番下の航法データから見ていきます。この航法データがGPS衛星が放送する情報になります。詳細はこの記事では立ち入りませんが、軌道要素と呼ばれる自身の軌道を計算するための情報や、時刻が含まれています。航法データは50bpsになります。

航法データは疑似ランダム符号を用いて変調されます。これをスペクトラム拡散と言います。航法データのみでは、ある周波数のスペクトルが生じます。航法データは50bpsなので、その信号周期に対応したところにスペクトルが生じてしまいます。全衛星が同じデータレートの航法メッセージを送信するので、この状態で複数の衛星の情報が重なると、混信してしまって信号が識別できなくなってしまうわけです。

そこで、ランダムな符号で変調を行います。これによって航法データの周期性が小さくなるため、様々な周波数成分を持つようになります。つまりランダムな符号で変調を行うことで情報を広範囲な周波数に薄める事ができます。広範囲の周波数に情報が薄められるので、先ほどとは異なりピークが生じなく（限りなく小さく）なります。

全くランダムな符号で変調してしまうと、復調するための手段がなくなってしまいます。そこで、周期性があり計算によって再現できる疑似ランダム符号で変調します。受信機では同じ疑似ランダム符号で復調することで、拡散されたデータを復元することができます。衛星ごとに異なる疑似ランダム符号を用いることで、どの衛星から送信された信号なのか識別することができます。完全なランダムではなく数式から復元できる擬似ランダム符号を使うことで、鍵のように衛星の識別に使うことができるわけです。

疑似ランダム符号で暗号化された信号は1575.42MHzのキャリア信号でBPSK(2位相シフトキーイング）で変調してアンテナから地球に向かって送信されます。

疑似ランダム符号とは一体何者か?

L1C/A信号で使う疑似ランダム符号はcoarse / acquisition コードといいます。coarseというのは雑なとか粗いという意味があります。

もう少し具体的に言うとゴールド系列と呼ばれる、疑似ランダム符号系列です。ゴールド系列は相互相関が非常に小さくなるランダム符号になります。

上記2行で全く情報が増えていませんが、大事なことは相互相関が非常に小さくなるランダム符号ということです。相互相関とは、2つの情報を比べたときにどれくらい似ているか、という指標のことです。したがって、相互相関が小さいということは似ていないということを意味します。L1C/Aで使うゴールド符号系列はそれぞれが全く似ていない疑似ランダム信号を作り出すことができるということです。

相互相関の計算方法は一旦おいておいて、ゴールド符号系列を作り出す回路を見ていきます。L1C/Aのゴールド符号系列は2つの10ビット線形帰還シフトレジスタで作ります。次のような回路になります。

2つのシフトレジスタは1.023MHzのクロックで更新されます。上部のシフトレジスタは2番目と10番目を取り出して、排他的論理和を取って1番目のレジスタにフィードバックします。

下部のシフトレジスタは2,3,6,8,9,10番目のレジスタを取り出して、排他的論理和を取って1番目にフィードバックします。

C/Aコードは下部のレジスタの取り出す位置（タップ）が衛星ごとにが異なります。衛星番号1の衛星は2番目と6番目をタップするというルールになっています。タップした2つのレジスタの排他的論理和を、上部のシフトレジスタの10番目のレジスタと排他的論理和を取ることでC/Aコード（疑似ランダムコード）が生成されます。ここからは衛星番号はPRN（Pseudo Random Number)とします。

C/Aコードの相互相関が小さくなることを確認してみましょう。ここではPythonを使って確かめてみます。

下のプログラムはC/Aコードを生成するプログラムです。

import numpy as np
import matplotlib.pyplot as plt

def shift(g1, g2):
  fb1 = g1[2]^g1[9]
  fb2 = g2[1]^g2[2]^g2[5]^g2[7]^g2[8]^g2[9]

  g1 = np.roll(g1, 1)
  g1[0] = fb1
  g2 = np.roll(g2, 1)
  g2[0] = fb2

  return (g1, g2)

def cacode(g1, g2, sat):
  return g1[9]^g2[sat[0]]^g2[sat[1]]

g1 = np.ones(10).astype('uint8')
g2 = np.ones(10).astype('uint8')

prn1 = np.zeros(1023)
prn2 = np.zeros(1023)

for n in range(1023):
  prn1[n] = cacode(g1, g2, [1,5])
  prn2[n] = cacode(g1, g2, [2,6])
  g1, g2 = shift(g1, g2)

x = np.arange(40)
fig = plt.figure(figsize=(6,4))
a1 = fig.add_subplot(211)
a1.step(x, prn1[:40])
a1.set_title("PRN1")
a2 = fig.add_subplot(212)
a2.step(x, prn2[:40])
a2.set_title("PRN2")
fig.tight_layout()
plt.show()

shift 関数が2つの線形帰還シフトレジスタを実現しています。関数 cacode にg1とg2、PRN番号を渡すことで、その時のC/Aコードの値を求めることができます。

2つの線形帰還シフトレジスタの初期値はすべて1です。そこで、numpyのones関数を用いて、初期化された配列を用意します。

C/Aコードの周期は1023なので、C/Aコードを求めてから、シフトする操作を1023回繰り返します。ここではPRN1とPRN2のC/Aコードを求めるために、PRN1のタップ位置2、6とPRN2のタップ位置3、7から1引いた値をタップ位置としてC/Aコードを求めています。1引く理由はPythonの配列がゼロオリジンのためです。Matlabなどの1オリジンの言語ではそのままの値を指定できます。shift関数内でも同様の理由で1引いたインデックスによりタップ位置を指定しています。

先頭40個のC/Aコードを表示すると下の図のようになります。

パット見で2つの符号は異なっていそうなことがわかります。

相互相関を計算する前にC/Aコードに少し細工をします。生成したC/Aコードは1と0の値を取っています。これが少し不便なので、0→1、1→-1に変換します。これでC/Aコードは(1,-1)の値を取ります。

変換のコードは下記のようになります。

prn1 = (1 - 2*prn1)
prn2 = (1 - 2*prn2)

同じようにグラフを書くと下の図のようになります。

相互相関は相関を計算したい2つの要素のそれぞれ掛け算して、掛け算の結果を足し合わせたものです。これを1個ずつずらしながら計算していきます。ずらしていって同じ状態のものが現れるとその時の値が極端に大きくなります。それではPRN1とPRN2の相互相関を計算してみます。計算のコードは下記のようになります。

coef = np.zeros(1023)

for n in range(1023):
  coef[n] = np.sum(prn1*np.roll(prn2, n))

x = np.arange(1023)
fig = plt.figure(figsize=(6,4))
a1 = fig.add_subplot(111)
a1.step(x, coef)
a1.set_title("Cross-correlation result")

fig.tight_layout()
plt.show()

コードのコアはnp.sum(prn1*np.roll(prn2, n))になります。np.roll(prn2, n)がずらす処理をしています。ここではprn2をn個ずらす処理をしています。その結果をprn1と掛け算し、np.sumで掛け算の結果をすべて足し合わせています。下の方はグラフを描く処理です。

相互相関の結果は次の図のようになります。うん、さっぱりわかりませんね。なんとなく全体に似ている箇所があるようなないような…。最大値と最小値は(60,-60)くらいです。先に述べておきますと、この値は全く似ていないときの値です。

それでは次に自己相関を見ていきたいと思います。自己相関は自分自身で行う相互相関のことです。自己相関は先程のnp.rollの中身を書き換えるだけで実現できます。下記に自己相関のコードを示します。

auto_coef1 = np.zeros(1023)
auto_coef2 = np.zeros(1023)

for n in range(1023):
  auto_coef1[n] = np.sum(prn1*np.roll(prn1, n))
  auto_coef2[n] = np.sum(prn2*np.roll(prn2, n))

x = np.arange(1023)
fig = plt.figure(figsize=(6,4))
a1 = fig.add_subplot(211)
a1.step(x, auto_coef1)
a1.set_title("Auto-correlation for PRN1")
a2 = fig.add_subplot(212)
a2.step(x, auto_coef2)
a2.set_title("Auto-correration for PRN2")
fig.tight_layout()
plt.show()

ここではPRN1とPRN2の自己相関を計算してみました。これをグラフにすると次のような図になります。

ひと目見てわかるように0、つまりズレが無いところで1000以上の値を叩き出しています。それ以外では先程と同じように60程度の値です。このように1000個以上のずらして比較しても、自身と全く同じものでないと相関の値は小さくなります。

C/Aコードの性質をまとめます。

• 異なるPRNの符号列は相互相関が小さくなる。
• 同じPRNの自己相関はズレが一致したときに最大になる。

この性質があることで、航法メッセージをC/Aコードで符号化すると、同じキャリア周波数で送信しても、受信機では混信せずにそれぞれの航法メッセージを復調することができます。

ちなみに衛星ごとのタップ位置はGPSの仕様書に記載されています。GPSの仕様書はGPS Interface Specificationといいます。そこからPRNごとのタップ位置を示す表を引用します。表の中のCode phase selectionがシフトレジスタのタップする位置を示しています。PRN1は2と6です。

GPS受信機の構成

だいぶ長くなってきて疲れてきました。しかしここからが本番です。

ここからはGPSの受信機について見ていきます。GPSの受信機は下の図のような構成をしています。

GPS受信機は大きく、

• RFフロントエンド
• ベースバンドプロセッサ
• PVTエンジン

という3つのブロックからできています。

それぞれを見ていきます。

RFフロントエンド

RFフロントエンドはアンテナで捉えられた信号が真っ先に入ってくる場所です。ここでは1547.42MHzの高周波を処理するブロックです。はじめにLNA（Low noise amplifier)で信号を増幅します。GPS衛星はMEO（Medium Earth Orbit）という高度20,200kmという果てしなく遠い場所を飛んでいます。そこから飛んでくる電波はかなり弱い電波です。そこでまずはアンプで増幅してあげるわけです。

増幅された信号はBPF（Band pass filter)を通して必要な信号だけを取り出します。今回はL1C/Aを取り出すので、1575.42MHzだけを通すようなフィルターを通過させます。

その後にあるマルにばってん2個は掛け算機です。回路的にはミキサーと呼ばれる回路になります。2個あるのは位相0度の信号と90度の信号2つをそれぞれ掛け算するためです。この0度と90度の信号は1575.42MHzの信号です。難しい話は省略しますが、これによってGPSの信号を低い周波数、今回ではほぼ直流の信号に変換できます。このように周波数を変換することをコンバージョンといい、今回は周波数を下げるのでダウンコンバージョンといいます。高い周波数の信号を扱うのはとても大変です。しかも、GPSの信号はほぼC/Aコードのチップレートが1.023MHzなので、帯域4MHz程度もあれば十分すぎるくらいです。そこで扱いやすい周波数までダウンコンバージョンします。ダウン子バージョン後の信号はIとQの2つになります。それぞれIn-phase、Quadrature phaseの略で2つ合わせてIQ信号と呼ばれます。

0度と90度2つで周波数変換するこの回路はIQ復調器やIQ変調器と呼ばれます。無線通信では必須の回路の1つです。90度位相差のある信号は直交しているので、とても良い性質があります。まぁ、そのへんの話は長くなるので、ここでは省略します。気になる方は信号処理、無線通信、線形代数の本などを読んでもらうといいのかなと思います。

ダウンコンバージョンした信号はまだ高い周波数成分も含まれているので、LPF(Low pass filter)で低い周波数成分のみを残すようにフィルタリングします。

最後にADC(Analog to digital converter)でデジタル信号に変換します。

このように必要な信号成分のみが含まれる状態になった信号をよくベースバンドといいます。今回はIQ復調器でキャリア信号の1575.42MHzを直接掛け算しましたが、もう少し高い周波数成分に変換することもあります。そのときのミキサー出力の信号はIF(Intermidiate Frequency)信号といいます。また、今回のように直接キャリア信号の1575.42MHzを掛け算することをIFがゼロになることから、Zero-IFと呼んだりします。

ベースバンドプロセッサ

ベースバンドプロセッサではその名の通りベースバンド信号を色々処理して、GPS衛星の送信する航法メッセージを復調するまでを行います。ベースバンドプロセッサはすべてをデジタル回路として構成されています。ほとんどの受信機ではASICとして専用のデジタル回路で処理されています。この図ではベースバンドプロセッサが1つしか書かれていませんが、受信する衛星の数だけこの回路が必要になります。測量に用いたりRTKを行うような受信機はこのベースバンドプロセッサが百個の単位で積まれていますが、安価な受信機は十個程度になります。同時受信できる衛星数がそのまま測位精度に直結するため、どれだけベースバンドプロセッサを積めるかがそのまま性能にもなります。また、最近ではADC以降は直接コンピュータに取り込んで処理するSDR(Software defined radio)という技術も使われることがあります。コンピュータに取り込むとGPUが使えるためかなり高速化したり、処理を柔軟に変えたり色々できるようになります。

ベースバンドプロセッサに入ってきた信号はもう一度IQミキサーによる周波数変換を行います。これはドップラーシフトのためです。下の図は地球上の受信機とGPS衛星の位置関係を示しています。

地球上から見て地平線から昇ってきたGPS衛星は真上で一番近くなり、また地平線に沈むときに遠くなっていきます。電波も波ですからこの変化でドップラーシフトが起こります。救急車が近づいて遠ざかっていくときのサイレンのように周波数が変化するわけです。このドップラーシフトの補正のためにベースバンドプロセッサではもう1回周波数変換を行います。ドップラーシフトによる周波数変化は±5kHz程度です。この図では受信機が固定されていますが、受信機も移動する場合それにより生じるドップラーシフトの補正が必要になります。

ドップラーシフトが補正された信号はC/Aコードと掛け算されて復調されます。この部分は掛け算機とシグマ(数学の足し合わせの記号です)で構成されています。つまり相関器となっているわけです。C/A code generator で受信したい衛星のC/Aコードを発生させて、相関器を通すわけです。相関器は1つIQそれぞれの信号に3つずつあります。それぞれ、E(Early)、P(Punctual)、L(Late)となっています。Punctualがぴったり一致したC/Aコード、EarlyとLateが少し進んだ、少し遅れたC/Aコードを生成します。EarlyとLateの2つの信号はPunctualのC/Aコードがどの遅れているのか、進んでいるのかを判断するために使用します。そして判断した結果としてC/Aコードを遅らせたり、進めたりすることでPunctualのC/Aコードがピッタリ一致するように調整しています。これをDelay lock loopといいます。

Punctual の信号はFLL/PLLに使用されます。これは先程述べたドップラーシフトを調整する処理のことです。それぞれFrequency lock loopとPhase lock loopです。航法メッセージはBPSKで変調されているので、位相をピッタリ合わせる必要があります。そのためにCostasループというものも使われます。CostasループはPLLの一種です。

PunctualのI信号はPLLが位相同期している状態で航法メッセージを出力します。Navigation dataブロックでは受信した信号から同期信号を見つけて、航法メッセージとして意味ある状態まで整えてから後段のPVTエンジンに情報を渡します。

Position, Velocity, Timing engine

このブロックはPVTエンジンと呼ばれます。その名の通り航法メッセージやDLLのデータを用いて測位を行い、速度を求め、時刻を作り出す処理を行います。ここは難しい話がたくさんあるので、ISHI会の同人誌にお譲りします。

GPS受信機のまとめ

GPSの受信機は大まかに3つのブロックがあることを説明しました。高周波信号を処理するためのRFフロントエンド、ベースバンド信号を処理して航法メッセージを求めるベースバンドプロセッサ、航法メッセージを用いて測位を行うPVTエンジンです。

ベースバンドプロセッサ以降は衛星が見つかってからの処理になります。まずは衛星を見つけないと行けません。GPS衛星の数だけ専用のベースバンドプロセッサを持てると良いですが、そういうわけにもいきません。GPSは衛星は基本的に30機以上が運用されています。流石にシリコンの面積としてもこれだけの数のベースバンドプロセッサを乗せることは難しいです。また、新しい衛星が増えたり、古い衛星が減ったりするとベースバンドブロックが使えなくなったり不足したりしてしまいます。

そこで現在自分の視野内にある衛星を捜索して、その情報からベースバンドブロックに追跡させるということを行います。

GPSの信号を捜索する

ここからはやっとタイトルにある「GPSの信号を捜索する」ということをやっていきます。

GPS受信機が起動した直後に自身がわからない情報をまとめてみます。

まず現在時刻がわかりません。つまり、どこにGPS衛星が飛んでいるのかわからないので、今見える衛星は不明です。

見えている衛星がどの程度のドップラーシフトがあるかわかりません。

現在時刻が不明なので、見えている衛星のC/Aコードの始点がわかりません。つまり、どの程度の遅延でC/Aコードを設定したらいいのか不明です。

つまり調べたいことは次の3つです。

• 見えいている衛星
• その衛星のドップラーシフト
• その衛星のC/Aコードの遅延

これを探します。

どのように捜索したらいいのか?

答えは力技です。すべての衛星のC/Aコードを作って、相互相関を考えられるドップラーシフトの範囲すべてでやります。

C/Aコードの遅延、ドップラーシフトが一致したときに相関の値が大きくなるので、それが今見えている衛星、C/Aコードの遅延、ドップラーシフトになります。

スマートにやる方法がないわけでもないのですが、C/Aコード、受信した信号をフーリエ変換してから相互相関をするという処理です。FFTを使うことで並列処理が可能になります。

今回は最も単純に時間領域での相互相関で衛星の捜索を行いたいと思います。処理時間はかかりますが、ハードウェアでFFTを行うより回路が単純になるので、面積が小さくなります。シリコンの面積＝お金です。安くなるようにいきましょう。

また今回は原理検証なので、Pythonとnumpyを最大限に利用します。SystemVerilog等で書くことはできますが（現に書いています）、RTLレベルのシミュレーションはとんでもなく時間がかかります。ええ、FPGAで実行するほうが早くね???と思うレベルです。

原理検証を行うにはデータが必要になります。今回は東京海洋大学の高須先生が開発されているPocketSDR及びそのソフトウェアに同梱されているデータを利用させていただきます。

データの読み込み

データはサンプリングレート4MHz、Zero-IFのIQデータです。(1,-1)のIQデータ交互に1バイトで格納されています。奇数番号のデータがQ、偶数番号のデータがIのように取り出していくと簡単に使うことができます。

データの読み込みは次のように行います。

import numpy as np
from google.colab import drive
drive.mount('/content/drive', force_remount=True)

fs = 4.0e6

raw = np.fromfile('/content/drive/MyDrive/L1_20211202_084700_4MHz_IQ(1).bin', dtype=np.int8)
i = raw[0::2]
q = raw[1::2]

samples = len(i)
print("Loaded {0:.1f} ms".format(samples/fs*1000))

今回はGoogle Colabを利用しているので、そのためのコードがいくつか入っています。np.fromfileでファイルを読み込みます。1つ目の引数がファイル名になるので、自分の環境に合わせて決めてください。

fsがデータのサンプリングレートで今回は4MHzのため、4e6としています。

IQの順番で交互にデータが入っているので、偶数番のインデックスをI、奇数番のインデックスをQのデータとします。

今回のデータでは500ms分のデータが入っています。

C/Aコードの生成

次にC/Aコードを生成するコードを書きます。C/Aコード生成はさっき書いたのでそのまま使います。今回はとりあえずPRN31の衛星に対して処理を行ってみます。ちなみにPRN31はこのベースバンド信号に含まれている衛星です。PRN31のタップは3と8です。そこで、2と7をcacode関数には指定します。

prn31 = np.zeros(1023)

for n in range(1023):
  prn31[n] = cacode(g1, g2, [2,7])
  g1, g2 = shift(g1, g2)

これでC/Aコードが生成できました。

C/Aコードのチップレートは1.023MHzで1023で1周するなので、1周期は1msです。したがって1ms分のデータで相互相関を取れば良いことになります。このデータのサンプリングレートは4MHzなので、4000点分のデータで相関を取るとぴったりです。C/Aコードもちょうど4000点になるように引き伸ばしておきましょう。ただ単純に引き伸ばしてしまうと毎回C/Aコードを生成し直さなくてはならないので、そのタイミングで参照するインデックスのリストを作成することで対応します。そのコードは次のようになります。

chip_rate = 1.023e6 # C/Aコードのチップレート(Hz)
code_len = 1023 # C/Aコード1周期の長さ
fcw = chip_rate / fs # ベースバンド1サンプルに対するC/Aコードの更新速度
coherent_time = 1e-3 # 相互相関する時間 (s)
num_coherent_data_sample = int(fs*coherent_time) # 相互相関するサンプル数→4000点
chip_index = (np.floor(np.arange(num_coherent_data_sample) * fcw) % code_len).astype(int)

とても複雑に見えますが、一番最後の行以外は定数を指定しているところです。それぞれの定数の意味は#以降にコメントとして入れてあります。

np.arangeは0から3999までの配列を作ります。これにfcwを掛けることで、あるタイミングのサンプルにおけるC/Aコードのインデックスになります。ただ、インデックスは整数でなければならないので、np.floorで切り捨てて整数にしてしまいます。最後のcode_lenで余りを取っている部分はもしもの時用です。インデックスに1022より大きな数が入るとアクセスエラーが生じます。そこで、コードの長さで余りを取ることで1022超えていたら、0からに戻す役割をしています。1022を超えないことが保証されていれば不要です。

ドップラーシフトの処理

次にドップラーシフトを補正する処理を書きます。ただ、ドップラーシフトがいくつくらいかわからないので、±5kHzの範囲を500Hzステップで補正するような処理とします。ドップラーシフト補正用のキャリア信号は次のように生成します。

fd = 500 # ドップラー周波数
n = np.arange(samples)
carrier_i = (np.sign(np.cos(2*np.pi*fd/fs*n))).astype(np.int8)
carrier_q = (np.sign(np.sin(2*np.pi*fd/fs*n))).astype(np.int8)

fdがドップラー周波数です。捜索するときはここを書き換えていきます。

carrier_iとcarrier_qがそれぞれのドップラーシフト補正用の信号です。全サンプルのデータを生成しています。np.signは符号を求める関数です。ベースバンド信号は1ビットや2ビットです。そうなると、キャリア信号も同じくらい少ないビット数でも航法メッセージの復調には十分です。そこで符号だけにしてしまいます。

ドップラーシフトの補正はIQミキサーで行うので、2つの信号を掛け算します。コードは次のようになります。

i_mixed = carrier_i*i
q_mixed = carrier_q*q

とても単純ですね。ただの掛け算です。

相互相関処理

相互相関はドップラーシフトを補正したIQそれぞれの信号に行います。相互相関は次のように記述します。

non_cohnum = 8 # ノンコヒーレント積分時間(ms)

power_sum = 0
for code_delay in range(code_len):
  local_code = np.roll(prn31, code_delay)[chip_index]

  power_sum = 0.0
  for blk in range(non_cohnum):
    start = blk * int(fs * coherent_time)
    stop = start + int(fs * coherent_time)
    i_corr = np.sum(i_mixed[start:stop] * local_code)
    q_corr = np.sum(q_mixed[start:stop] * local_code)

    power_sum += i_corr*i_corr + q_corr*q_corr

単純な相互相関でないことがわかるかと思います。外側のforはC/Aコードに遅延を与える処理です。これは前に説明した相互相関の処理と変わりません。local_codeはこのときのずらしたC/Aコードです。np.rollでずらしてから、前に作成したchip_indexでアクセスすることで、1023点のC/Aコードを4000点まで引き伸ばしています。

内側のforが少し複雑になっています。これは1msの相互相関をデータを1msずつずらしながら8回行う処理です。8回はコードの先頭にあるnon_cohnum変数で決めます。内側のforのstartとstopが相互相関する範囲を決めています。4000点ずつforの各回でずれるようになっています。

相互相関の結果はIとQそれぞれあるので、それをベクトルのように捉えて2乗してから足し算することでベクトルの大きさを出しています。各回の処理結果はpower_sumに加算されます。このような処理をノンコヒーレント積分といいます。

なぜ単純に8ms分のデータを相互相関しては駄目なのでしょうか?それはRFフロントエンドに用いる発振器の安定度のためです。ここで使っている発振器は水晶発振器です。ジッターは少ないですが、安定度はそこまで高くありません。また周波数がズレている可能性もあります。すると正確に4MHzでサンプリングできていない可能性があるのです。というより、できていません。したがって8msもの長期間で相互相関を行うと、計算上のC/Aコードとベースバンド信号のC/Aコードがズレてきてしまうのです。そうなると正確な相互相関ができず、相関が小さいままになってしまい衛星を発見できません。そこでコード1周期分の相互相関を8回繰り返すことでズレをリセットしながら相互相関をするのです。このようにズレがリセットされる、つまり位相が保存されないのでノンコヒーレントといいます。コヒーレントとは位相が揃っている状態を指します。つまり通常の相互相関は位相が揃っている前提なので、コヒーレント積分になります。ここでは8回のコヒーレント積分をノンコヒーレント積分しているわけです。

全部の処理をガッチャンコしてPRN31の衛星があるか確かめる

ここまでで必要な要素が全て揃いました。それではすべてをガッチャンコして衛星を捜索してみましょう。ガッチャンコしたコードはこのようになります。定数を上にまとめたり、ドップラーシフトの総当り、グラフの描画を追加しています。

import numpy as np
import matplotlib.pyplot as plt
from google.colab import drive
drive.mount('/content/drive', force_remount=True)

# 定数とか
fs = 4.0e6 # ベースバンド信号のサンプリングレート(Hz)
chip_rate = 1.023e6 # C/Aコードのチップレート(Hz)
code_len = 1023 # C/Aコード1周期の長さ
fcw = chip_rate / fs # ベースバンド1サンプルに対するC/Aコードの更新速度
coherent_time = 1e-3 # 相互相関する時間 (s)
num_coherent_data_sample = int(fs*coherent_time) # 相互相関するサンプル数→4000点
non_cohnum = 8 # ノンコヒーレント積分時間(ms)
f_doppler_candidate = np.arange(-5000, 5000, 500) # ドップラー周波数の探索範囲

# C/Aコードの生成
def shift(g1, g2):
  fb1 = g1[2]^g1[9]
  fb2 = g2[1]^g2[2]^g2[5]^g2[7]^g2[8]^g2[9]

  g1 = np.roll(g1, 1)
  g1[0] = fb1
  g2 = np.roll(g2, 1)
  g2[0] = fb2

  return (g1, g2)

def cacode(g1, g2, sat):
  return g1[9]^g2[sat[0]]^g2[sat[1]]

g1 = np.ones(10).astype('int8')
g2 = np.ones(10).astype('int8')

prn31 = np.zeros(1023)

for n in range(1023):
  prn31[n] = cacode(g1, g2, [2,7])
  g1, g2 = shift(g1, g2)

# C/Aコードを4000点まで引き伸ばす
chip_index = (np.floor(np.arange(num_coherent_data_sample) * fcw) % code_len).astype(int)

# ベースバンドデータの読み込み
raw = np.fromfile('/content/drive/MyDrive/L1_20211202_084700_4MHz_IQ(1).bin', dtype=np.int8)
i = raw[0::2]
q = raw[1::2]

samples = len(i)
print("Loaded {0:.1f} ms".format(samples/fs*1000))

# 結果の格納用変数の作成。ドップラー周波数とコード遅延の2軸あるので、2次元の配列にする。
corr_map = np.zeros((len(f_doppler_candidate), code_len))

# 力技処理の部分。ドップラーシフトとコード遅延を総当りする。
for fi, f_doppler in enumerate(f_doppler_candidate): # f_dopplerがドップラー周波数
  print(f"Search Doppler frequency {f_doppler} Hz")
  n = np.arange(samples)
  # ドップラーシフト補正用キャリア信号生成
  carrier_i = (np.sign(np.cos(2*np.pi*f_doppler/fs*n))).astype(np.int8)
  carrier_q = (np.sign(np.sin(2*np.pi*f_doppler/fs*n))).astype(np.int8)

  # ドップラーシフトの補正
  i_mixed = carrier_i*i
  q_mixed = carrier_q*q

  power_sum = 0
  # コヒーレント積分
  for code_delay in range(code_len):
    local_code = np.roll(prn31, code_delay)[chip_index]

    power_sum = 0.0
    # ノンコヒーレント積分
    for blk in range(non_cohnum):
      start = blk * int(fs * coherent_time)
      stop = start + int(fs * coherent_time)
      # 相関処理
      i_corr = np.sum(i_mixed[start:stop] * local_code)
      q_corr = np.sum(q_mixed[start:stop] * local_code)
      # ノンコヒーレント積分の足し合わせ処理
      power_sum += i_corr*i_corr + q_corr*q_corr
    corr_map[fi, code_delay] = power_sum

# 結果表示部分
max_fd_index, max_code_index = np.unravel_index(np.argmax(corr_map), corr_map.shape)
print(f"Detected peak: Doppler={f_doppler_candidate[max_fd_index]} Hz, Code phase = {max_code_index}, Corr = {corr_map[max_fd_index, max_code_index]}")

fig = plt.figure()
ax = fig.add_subplot(211)
im1 = ax.imshow(corr_map.T, aspect='auto',
                extent = [f_doppler_candidate[0], f_doppler_candidate[-1], 0, 1023],
                origin='lower', cmap = 'inferno')
ax.set_xlabel("Doppler (Hz)")
ax.set_ylabel("Code phase (chips)")
fig.colorbar(im1, label='Correlation')

ay = fig.add_subplot(212)
ay.plot(np.arange(1023), corr_map[max_fd_index, :])
ay.set_xlabel("Code phase (chips)")
ay.set_ylabel("Correlation")
fig.tight_layout()
plt.show()

このコードを実行すると…

Mounted at /content/drive
Loaded 500.0 ms
Search Doppler frequency -5000 Hz
Search Doppler frequency -4500 Hz
Search Doppler frequency -4000 Hz
Search Doppler frequency -3500 Hz
Search Doppler frequency -3000 Hz
Search Doppler frequency -2500 Hz
Search Doppler frequency -2000 Hz
Search Doppler frequency -1500 Hz
Search Doppler frequency -1000 Hz
Search Doppler frequency -500 Hz
Search Doppler frequency 0 Hz
Search Doppler frequency 500 Hz
Search Doppler frequency 1000 Hz
Search Doppler frequency 1500 Hz
Search Doppler frequency 2000 Hz
Search Doppler frequency 2500 Hz
Search Doppler frequency 3000 Hz
Search Doppler frequency 3500 Hz
Search Doppler frequency 4000 Hz
Search Doppler frequency 4500 Hz
Detected peak: Doppler=-500 Hz, Code phase = 296, Corr = 1492903.0

このように表示されます。そしてグラフは…

このようになります!!

上のマッピングはわかりにくいですが、ドップラー周波数0Hz、コード300付近に明るい細い線があるのが目を凝らすとわかると思います。

この部分をグラフを垂直に切り抜いたのが下のグラフになります。

なんとコード位相296のところにピークが現れています。これが相関が大きくなる。つまり、このベースバンド信号にPRN31の衛星が含まれていることを示しているのです。

まとめ

GPS衛星の信号を捜索するというテーマでGPSの信号ってなに?ってところから受信機の構成を説明して、acquisition(アクイジションと読みます)の処理をPythonで書いてみました。この後には発見した衛星をトラッキング、航法メッセージの復調と擬似距離の計算、PVTエンジンでの測位とまだまだやらないといけないことがたくさんあります。普通にはGPSモジュールを買ってきてUARTで出力されるNMEAをマイコンで読み取れば事足りると思います。でもそれじゃ面白くない。車輪の再発明上等じゃいってことでオレオレGNSS受信機を作りたいのです。

来年にQZSSの7号機が打ち上がります。これにより日本上空に常時4機以上のQZSS衛星が滞在することになります。これが何を意味するかというと、QZSS単体で測位ができるようになるのです。これまでQZSSはGPSやそれ以外のGNSSの補完衛星でしかありませんでした。つまり物理的にQZSSしか受信しない受信機を作ることも可能になるのです。これってすごく面白いと思いませんか?

協定世界時というものがあります。UTCと呼ばれるものです。これは世界中の原子時計を集合して作られる国際原子時TAIに天文時UT1から決まる閏秒を加算したものになります。この世界中の原子時計を集合する、つまり比較するためにGNSSが使われています。GNSSコモンビューといいます。これを行うための専用の受信機も売られたりしています。GNSS衛星は原子時計を搭載しているため、その周波数の確度は非常に高いです。つまり、GNSSからの信号は基準周波数として使うことができます。また、航法メッセージに含まれる時刻はUTCとほとんど同期しています。つまりGNSSを受信するだけで高精度な時刻・周波数源ともなるわけです。しかし、市販の受信機は100MHz程度のクロックで動いているので1ppsの信号には10~50ns程度のジッターがあります。もしオレオレ受信機ICを作れるなら1pps専用の遅延回路(digital to time converter)を搭載して、psオーダの精度で1pps信号を出力したりもできるかもしれません。IC自作にはそんな夢を見ることができる魅力的な世界なのです。

ISHI会ではIC自作に夢を見るそんな人達がたくさん活動しています。全くの初心者がチュートリアルで自分のICを作る機会も定期的に開催されています。集積回路が専門の大学の先生による大学院レベルの講義をタダで聴講して、それで勉強した内容でICを作ったりもできます。IC自作に興味があるそんな人は是非ISHI会Discordに参加してみてください。

参考文献

• GNSS 受信機の概要 I アンテナ信号入力から、擬似距離出力まで
• 新しいフェーズに入った衛星測位技術を加速させる人材育成　プログラム及び教材
• PocketSDR
• Homemade GPS Receiver
• KiwiSDR design review
• Analog-GPS-data-receiver
• A homemade receiver for GPS & GLONASS satellites
• Understanding GPS Principles and Application, Second Edtion
• A Software-Defined GPS and Galileo Receiver - A Single Frequency Approach
• FUNDAMENTALS OF GLOBAL POSITIONING SYSTEM RECEIVERS
• GLOBAL POSITIONING SYSTEM STANDARD POSITIONING SERVICE SIGNAL SPECIFICATION
• Quasi-Zenith Satellite System Interface Specification Satellite Positioning, Navigation and Timing Service (IS-QZSS-PNT-006)
• MAX2769B