注意：LimeSDR などのソフトウェア無線は技術基準適合証明（通称，技適）を受けていないため，これらの機器を日本国内で無線機として利用することは電波法により禁じられています。無線機として使用するためには，実験試験局免許を取得するか電波暗室・シールドボックスなどの設備を使用する必要があります。または，アンテナの代わりにケーブルとアッテネータを用い有線接続をし電波を発しないようにすることで，無線通信ではなく有線通信とはなりますが実験することができます（この場合も電波が漏れないように注意してください）。このページを参照される方は，実験される国や地域の法令などを遵守するようにご注意ください。また，実験等はご自身の責任でお願いします。

5日目 で信号サンプルを受信するループを書きました。今日からこの受信した信号サンプルのプリアンブルとユニークワードから基準信号の位相を推定し復調をするプログラムを実装していきます。今日は位相を推定するため，プリアンブルを検出して受信信号の位相を可視化します。その可視化のため matplotlib を使います。プログラムの最初に以下を追加します。

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

3日目 で取り扱ったように，送信プログラムでは以下の信号を送っています。

• 基準信号との位相差0の正弦波を5000サンプル
• プリアンブル 0101010101010111 16ビット（160サンプル）
• ユニークワード 0010101111010100 16ビット（160サンプル）
• データ 0111010111110001 16ビット（160サンプル）
• 基準信号との位相差0の正弦波を5000サンプル

このうち，プリアンブルは 0 と 1 が14ビット分交互に続きます。プリアンブルの前の正弦波から位相を推定することもできますが，位相の変化がない信号は無信号（ノイズのみの信号）と区別が付きにくいため，プリアンブルを基準にします。

なお，readStream() で受信ストリームから読み出す信号サンプルは，読み出しタイミングやバッファサイズの関係でプリアンブルの途中となる可能性があります。そのため，2回分のバッファを連結して処理します。

5日目 で説明した通り，readStream() で読み出した信号サンプルは受信ストリームのバッファが埋まらないように高速で処理する必要があります。Pythonで1サンプルずつビットを推定すると最低サンプルレートまで落としても間に合わないため，以下の方針でプリアンブルの開始点および基準信号の位相を推定します。

1. 2回分のバッファを連結した信号サンプル列を取り出す
2. 1つ前のシンボル（デフォルトでは10サンプル前）との位相差を計算し，プリアンブルが含まれる可能性のある信号サンプル列か判定する
3. 位相が反転するサンプルからプリアンブルを検出し，プリアンブルの最後の 111 およびユニークワードのパターンから位相を推定する

まずは，ループ内で2回分のバッファを連結した信号サンプルを取り出します。コードとしては以下のようになります。

    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])

        ## (X) ここでプリアンブルの検出と位相推定を行う

        prevSamples = np.copy(rxBuffer)

(X) の場所にプリアンブルの検出と基準信号の位相推定のコードを書いていきます。

上記の samples は readStream() 2回分の信号サンプル列となります。以下のコードで，1つ前のシンボルである args.samples_per_symbol サンプル分前（デフォルトでは10サンプル前） prev との cur の位相差を計算します。位相差は複素平面における偏角なので， cur / prev の虚部と実部に対する numpy.arctan2 関数で計算できます。この位相差を diffAngles とします。

        # 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)

ここで，プリアンブルが信号サンプル列に含まれる場合は少なくとも 14 * args.samples_per_symbol だけの位相反転が含まれますが，プリアンブルの先頭を落とすことを想定して，8 * args.samples_per_symbol 以上の位相反転が含まれる場合にプリアンブルが含まれる可能性のある信号とします。これ以外の場合については，プリアンブルが含まれる可能性がないとして無視します。以下のコードのように，位相が (90^\circ ) より大きく変わるものを位相反転として扱い，プリアンブルが含まれる可能性がある場合の処理を実装します。

        # Check if the preamble is included
        if np.sum(np.absolute(diffAngles) > math.pi / 2) >= 8 * args.samples_per_symbol:

位相の反転が一定数以上含まれている信号サンプル列はプリアンブルが含まれている可能性があるため，今回はこの最初の位相反転から600サンプルを描画します。プリアンブルとユニークワードの検出，復調はは明日以降実装します。

まず，最初に位相差列から最初に位相が反転する（正確には (90^\circ\) 以上の位相変化がある）点を以下のコードにより検出します。

            firstIndex = np.where(np.absolute(diffAngles) > math.pi / 2)[0][0]

プリアンブルの先頭（の可能性がある位相反転したサンプル）の場所によっては，プリアンブル，ユニークワード，データが次のバッファに含まれていることもあるため，

            # Concatenate one more sample buffer
            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])

により，もう1回分 readStream() で読み出したバッファを追加します。

この信号サンプル列 sample から振幅と位相を以下のコードで求めます。

            amps = np.hypot(samples.real, samples.imag)
            angles = np.arctan2(samples.imag, samples.real)

以下のコードで位相反転した地点から600サンプル分の振幅と位相をそれぞれ緑と青の線により描画し，figure.png に保存し，ループを終了します。

            # Plot
            NR = 600
            ts = np.arange(NR)
            fig = plt.figure()
            plt.plot(ts, amps[firstIndex:firstIndex+NR], 'g', label='r')
            plt.plot(ts, angles[firstIndex:firstIndex+NR], 'bo-', label='theta')
            plt.xlabel('Samples')
            plt.ylabel('Amplitude / Phase')
            plt.savefig('figure.png')
            break

この受信コードは bpsk_rx_draw.py に置きました。送信プログラムである bpsk_tx.py を実行した直後に bpsk_rx_draw.py を実行すると，以下のようなグラフを得ることができます。

青色の線が位相でプリアンブル，ユニークワード，データの信号が現れています。基準信号の位相を0としていないため，+側の位相が0，-側の位相が1となっています。また，搬送波周波数とクロックの精度が原因と思われますが，位相が徐々に減少していることがこの図からわかります。このように位相がドリフトしていくため，信号処理により補正する必要があります。

今日は受信した信号サンプル列からプリアンブルの開始点を推定し，それ以降の信号を描画しました。明日はこれを用い，基準信号の位相推定と復調を実装します。