Advent Calendar 2020: ソフトウェア無線
注意：LimeSDR などのソフトウェア無線は技術基準適合証明（通称，技適）を受けていないため，これらの機器を日本国内で無線機として利用することは電波法により禁じられています。無線機として使用するためには，実験試験局免許を取得するか電波暗室・シールドボックスなどの設備を使用する必要があります。または，アンテナの代わりにケーブルとアッテネータを用い有線接続をし電波を発しないようにすることで，無線通信ではなく有線通信とはなりますが実験することができます（この場合も電波が漏れないように注意してください）。このページを参照される方は，実験される国や地域の法令などを遵守するようにご注意ください。また，実験等はご自身の責任でお願いします。
Day 7: BPSKによる通信（基準信号の位相推定と復調）
6日目 ではBPSKで変調された信号を受信し，以下のような受信信号サンプルを描画しました。
今日はこの受信信号サンプルのプリアンブル，ユニークワードを分析することで基準信号の位相を推定し，それらに続くデータを復調します。今回はデータ内にデータ長などのヘッダが含まれていないので，とりあえずプリアンブル，ユニークワードに続くデータは（非常に効率が悪くなりますが）16ビットとします。
基準信号の位相推定
6日目 の受信信号サンプルを処理するループ内のうち，プリアンブルが含まれている可能性がある場合の処理で demodulate() 関数を呼ぶようにします。具体的には以下のように変更します。
prevSamples = np.array([], dtype=np.complex64)
while True:
# Receive samples
status = sdr.readStream(rxStream, [rxBuffer], rxBuffer.size)
if status.ret != rxBuffer.size:
sys.stderr.write("Failed to receive samples in readStream(): {}\n".format(status.ret))
return False
# Concatenate the previous samples and the current samples
samples = np.concatenate([prevSamples, rxBuffer])
# Calculate angles from previous symbols
cur = samples[args.samples_per_symbol:] # current symbols
prev = samples[0:-args.samples_per_symbol] # previous symbols
prev = np.where(prev==0, prev + 1e-9, prev) # to avoid zero division
x = cur / prev
diffAngles = np.arctan2(x.imag, x.real)
# Check if the preamble is included
if np.sum(np.absolute(diffAngles) > math.pi / 2) >= 8 * args.samples_per_symbol:
# Concatenate two sample buffers
samples = np.copy(rxBuffer)
status = sdr.readStream(rxStream, [rxBuffer], rxBuffer.size)
if status.ret != rxBuffer.size:
sys.stderr.write("Failed to receive samples in readStream(): {}\n".format(status.ret))
return False
samples = np.concatenate([samples, rxBuffer])
if demodulate(samples, args.samples_per_symbol):
break
prevSamples = np.copy(rxBuffer)
demodulate() 関数はプリアンブル，ユニークワードおよびデータが復調できた場合には True，それ以外の場合は False を返すようにします。今回はサンプルなので True が返された場合はそこで受信処理を終了します。
demodulate() 関数内ではまずは基準信号の位相推定と復調を行います。まずは，基準信号の位相推定について説明します。なお，6日目 では描画するだけでしたので，位相が反転する点からプリアンブルとして扱いましたが，今回はもう少し厳密に扱います。
基準信号の位相を推定するために，まずシンボルが切り替わるエッジを判定する機能を実装します。detectEdge() 関数は信号サンプルからシンボルのエッジのオフセットを検出します。具体的には，以下のコードのように，プリアンブルを含む信号サンプルのうち，位相反転が起こっているインデックスを最も起こるインデックスのオフセットを検出します。プリアンブルの先頭に 0 と 1 が交互に現れるパターンを入れたのはこのエッジの検出を行うためでもあります。
"""
Detect edge
"""
def detectEdge(samples, samples_per_symbol):
shifted = samples[1:]
orig = samples[0:-1]
orig = np.where(orig==0, orig + 1e-9, orig) # to avoid zero division
x = shifted / orig
diffAngles = np.arctan2(x.imag, x.real)
# Find the change points of a phase
changePoints = np.where(np.absolute(diffAngles) > math.pi / 2)
# Find the edge index
bc = np.bincount(changePoints % samples_per_symbol)
return np.argmax(bc)
NumPy の機能を活用しているため少し読みにくいものになっていますが，1サンプルずらした信号サンプル列 shifted と元のサンプル列 orig の位相差のリスト diffAngles を計算し，この位相差が ( 90^\circ ) よりも大きくなる点を changePoints として取り出します。この位相差が変化する点をシンボルあたりのサンプル数で剰余を取り，その最頻値を numpy.bincount と numpy.argmax により求めることで， samples に対するシンボル境界のオフセットを検出することができます。
この detectEdge() 関数を demodulate() 関数の最初に呼び，信号サンプル列からシンボルをデコードできるようにします。デフォルトの値の場合，1シンボルあたり10個の信号サンプルがあるので，複数のサンプルを使うことでエラー訂正に使えますが，今回は簡単のために中央のサンプル（5つ目のサンプル）をシンボルの値として使うこととします。コードとしては以下のようになります。
"""
Demodulate
"""
def demodulate(samples, samples_per_symbol):
# Detect the edge offset
edgeIndex = detectEdge(samples, samples_per_symbol)
# Decode symbols from the center signal
symbols = samples[range(edgeIndex + 4, samples.size, samples_per_symbol)]
次に，基準信号の位相を推定します。BPSK なので，受信信号の位相は基準信号から \( 0 \) または \( \pi \) ずれたものになります。基準信号の位相を推定するために，受信した信号サンプルの複素平面に対し主成分分析を行い，下図のように第1主成分（つまり分散が最も大きくなる軸）を取り出すことで，基準信号から \( 0 \) または \( \pi \) の位相（I-QダイアグラムにおけるI軸）を推定します。
上図からもわかるとおり時間変化によって位相が回転するため，主成分分析が最も適した手法とは言えませんが，今回は主成分分析を使いました。このPCAのコードは以下の通りです。第1主成分は base として得られます。
# PCA
c = np.cov(np.array([symbols.real, symbols.imag]))
e = np.linalg.eig(c)
v = e[1][:,np.argmax(e[0])]
base = v[0] + 1j * v[1]
if base == 0:
base = 1e-9
この base を基準信号の位相，または基準信号の位相が反転（ \( \pi \) ずれた）ものとして扱います。
復調
この base からの位相差により復調します。BPSK はI-Qダイアグラム上に2つの信号点のみを持つため，実部（I軸）を比較して，正であれば1，負であれば0として扱います。なお，位相が反転している場合は，正が0，負が1となります。まず，以下のコードでこの実部の正負を真偽値として得ます。
# Demodulate the symbols
demod = symbols / base
binary = np.where(demod.real >= 0, True, False)
次に以下のコードにより，シンボルからプリアンブルおよびユニークコードと一致する場所を探します。なお，位相が反転している場合を考慮して，プリアンブルおよびユニークコードを反転させたビット列とも比較しています。
# Detect the preamble
psync = np.array([False, True, False, True, False, True, True, True, False, False, True, False, True, False, True, True, True, True, False, True, False, True, False, False])
bstr = binary.tostring()
try:
idx0 = bstr.index(psync.tostring())
except:
idx0 = -1
try:
idx1 = bstr.index(np.logical_not(psync).tostring())
except:
idx1 = -1
プリアンブルおよびユニークコードと一致する場所が見つかった場合，位相が反転しているかどうかにより処理を変えます。位相が反転している場合は，シンボルの真偽値を反転させます。なお，位相が反転していないものと位相が反転したものの両方でプリアンブルおよびユニークコードと一致した場合は，シンボル列のより先頭で一致した方を優先しています。
if idx0 < 0 and idx1 < 0:
return False
elif idx0 < 0:
idx = idx1
binary = np.logical_not(binary)
elif idx1 < 0 or idx0 < idx1:
idx = idx0
else:
idx = idx1
binary = np.logical_not(binary)
最後に，以下のコードでプリアンブルおよびユニークコードよりも後ろの16シンボル（16ビット分）を出力して処理を終わります。
sys.stdout.write("Received data:")
for i in range(16):
if binary[idx + psync.size + i]:
sys.stdout.write(" 1")
else:
sys.stdout.write(" 0")
print("")
実験
以上の実装は bpsk_rx.py に置きました。送信機で bpsk_tx.py を実行後すぐに受信機で実行することで以下の出力を得ることができます。
$ python3 bpsk_rx.py
linux; GNU C++ version 7.3.0; Boost_106501; UHD_003.010.003.000-0-unknown
[INFO] Make connection: 'LimeSDR-USB [USB 3.0] XXXXXXXXXXXX'
[INFO] Reference clock 30.72 MHz
[INFO] Device name: LimeSDR-USB
[INFO] Reference: 30.72 MHz
[INFO] LMS7002M register cache: Disabled
[INFO] RX LPF configured
[INFO] RX LPF configured
[INFO] Rx calibration finished
Received data: 0 1 1 1 0 1 0 1 1 1 1 1 0 0 0 1
送信プログラムから送信したデータが正しく受信できていることがわかります。
今日のまとめと明日の予定
今日はBPSKの受信信号サンプル列から基準信号の位相推定と復調を行い，送信されたデータを復元しました。明日はこれまでのBPSKの送信および受信を組み合わせて送受信のプログラムを書いてみようと思います。