初めてのGNU Radioフローグラフ

これまでのレッスンでは、numpy配列全体を一度に処理してきました。実際の無線機はそうは動きません。サンプルはADCから途切れなく届きます。無線機が動いている間ずっと、処理はブロックごとにそれに追いつかなければなりません。GNU Radioは、このストリーミングモデルを中心に作られたオープンソースのソフトウェア無線フレームワークです。このレッスンは、レッスン6から8のnumpyによるリンクを理解していて、同じ種類のリンクをGNU Radioで作りたい読者を対象とします。フローグラフをPythonで書き、グラフィカルインタフェースなしで実行してビット誤りを数えます。そのあと、同じグラフがGNU Radio Companionでどう見えるかを確認します。本稿の内容はすべてGNU Radio 3.10.12で実行しました。

ブロック・ストリーム・トップブロック

GNU Radioのプログラムはフローグラフ、つまりストリームで接続されたブロックの集まりです。各ブロックは、フィルタリング、雑音の付加、ビットからシンボルへのマッピングといった一つの仕事をします。ストリームは、出力ポートから入力ポートへ、一つの型のアイテムを運びます。よく使うアイテム型は、バイト、32 ビット浮動小数点数、そして32 ビット浮動小数点数2つからなる複素数です。この複素数をGNU Radioではgr_complex、numpyではcomplex64と呼びます。ブロック名の末尾には型を表す接尾辞が付きます。_bcはバイトを受け取って複素数を出力します。_ccfは複素数を受け取って複素数を出力し、フィルタタップはfloatです。_ccは入出力ともに複素数です。

コンテナはgr.top_blockです。ブロックを作り、connectで接続し、run()を呼びます。スケジューラはアイテムのバッファをブロックからブロックへ渡します。GNU Radio APIリファレンスによれば、run()はグラフを開始してその終了を待ちます。フローグラフは、すべてのブロックが完了を報告したとき、または何かがstop()を呼んだときに終了します。つまりシミュレーションでは、終わりのあるソースがグラフに必要です。repeat=Falseのblocks.vector_source_*はデータを出し終えると終了し、blocks.headは決まった数のアイテムを通したあと停止します。終了したソースの下流にあるものもすべて終了します。blocks.vector_sink_*は届いたものをすべて集めるので、あとでリストとして読み出せます。

ほとんどのブロックは、レートによって3つのグループに分かれます。同期ブロック（sync block）は、加算器のように入力1アイテムにつき出力1アイテムを生成します。デシメータは入力D個につき出力1個を生成します。インターポレータは入力1個につき出力I個を生成します。ストリーム自体にはサンプルレートの設定がありません。レートは自分で追跡する数値であり、純粋なシミュレーションではグラフはCPUが許す限りの速さで動きます。時間の進み方を決めるのは、実際のサンプルクロックに結び付いたハードウェアのソースかシンクだけです。

作るリンク

このフローグラフは、ランダムなビットをQPSK変調器、ルートレイズドコサイン（RRC）パルス整形フィルタ、加法性白色ガウス雑音、RRC整合フィルタ、判定器の順に通し、ビットに戻します。レッスン6のリンクのBPSKをQPSKに置き換え、GNU Radioのブロックで組んだものです。

Pythonコードの各ブロックには、グラフィカルエディタであるGNU Radio Companion（GRC）に対応するブロックがあり、その名前を図に示しています。レートは公称値です。シンボルレート125 kbaud、1シンボルあたり8サンプルでは、フィルタ後のストリームは1 MHzの複素帯域幅を表します。ただし、シミュレーション自体はCPUが処理できる速さで動きます。

Pythonで書くフローグラフ

下のコードブロックはそれだけで完結しています。先頭行の# gnuradioは、当サイトのコードチェッカーにフローグラフであることを示す印です。雑音レベルはレッスン6と同じ単位エネルギーの規約に従います。RRCタップは二乗和が1になるようにスケーリングしているので、各QPSKシンボルはEs = 1、各ビットはEb = 0.5を持ちます。このとき1サンプルあたりの複素雑音分散はN0となり、analog.noise_source_cはその標準偏差を振幅として受け取ります。

# gnuradio
import numpy as np
from scipy.special import erfc
from gnuradio import gr, blocks, analog, digital, filter
from gnuradio.filter import firdes

sps = 8                     # samples per symbol
sym_rate = 125e3            # symbol rate, so the sample rate is 1 MHz
samp_rate = sps * sym_rate
n_bits = 200_000
ebn0_db = 6.0

# Gray-mapped QPSK: symbol index -> point. Indices 0..3 carry two bits each.
points = [(1 + 1j) / np.sqrt(2), (-1 + 1j) / np.sqrt(2), (1 - 1j) / np.sqrt(2), (-1 - 1j) / np.sqrt(2)]
const = digital.constellation_calcdist(points, [], 4, 1, digital.constellation.NO_NORMALIZATION)

# Root-raised-cosine taps, roll-off 0.35, 11 symbols long, scaled to unit energy.
ntaps = 11 * sps + 1
rrc = np.array(firdes.root_raised_cosine(1.0, samp_rate, sym_rate, 0.35, ntaps))
rrc /= np.sqrt(np.sum(rrc**2))
delay_sym = (ntaps - 1) // sps          # two filters, each (ntaps - 1) / 2 samples, = 11 symbols

# With unit-energy pulses each symbol carries Es = 1, so Eb = 0.5. Complex noise variance per sample = N0.
n0 = 0.5 / 10 ** (ebn0_db / 10)

rng = np.random.default_rng(9)
bits = rng.integers(0, 2, n_bits).astype(np.uint8)

tb = gr.top_block()
src = blocks.vector_source_b(bits.tolist(), False)          # one bit per byte, finite, so the graph ends
pack = blocks.pack_k_bits_bb(2)                             # 2 bits -> one symbol index 0..3
mapper = digital.chunks_to_symbols_bc(points)               # index -> complex point
tx_filt = filter.interp_fir_filter_ccf(sps, rrc.tolist())   # upsample by 8 and pulse-shape
noise = analog.noise_source_c(analog.GR_GAUSSIAN, float(np.sqrt(n0)), 42)
add = blocks.add_cc()                                       # channel: signal + noise
rx_filt = filter.fir_filter_ccf(sps, rrc.tolist())          # matched filter, keep 1 sample per symbol
decoder = digital.constellation_decoder_cb(const.base())    # nearest point -> index 0..3
unpack = blocks.unpack_k_bits_bb(2)                         # index -> 2 bits
sym_sink = blocks.vector_sink_c()                           # matched-filter output, for the SNR estimate
bit_sink = blocks.vector_sink_b()

tb.connect(src, pack, mapper, tx_filt)
tb.connect(tx_filt, (add, 0))
tb.connect(noise, (add, 1))
tb.connect(add, rx_filt, decoder, unpack, bit_sink)
tb.connect(rx_filt, sym_sink)
tb.run()

rx_bits = np.array(bit_sink.data(), dtype=np.uint8)
rx_sym = np.array(sym_sink.data())
skip = 2 * delay_sym                                       # filter delay in bits (2 bits per symbol)
n = min(rx_bits.size - skip, bits.size)
errors = int(np.sum(rx_bits[skip : skip + n] != bits[:n]))
ber = errors / n
ber_theory = 0.5 * erfc(np.sqrt(10 ** (ebn0_db / 10)))

tx_sym = np.array(points)[bits[: 2 * (n // 2)].reshape(-1, 2) @ [2, 1]]
err_vec = rx_sym[delay_sym : delay_sym + tx_sym.size] - tx_sym
snr_db = 10 * np.log10(1 / np.mean(np.abs(err_vec) ** 2))

print(f"sample rate {samp_rate / 1e6:.1f} MHz, symbol rate {sym_rate / 1e3:.0f} kbaud, filter delay {delay_sym} symbols")
print(f"{n} bits compared, {errors} errors")
print(f"BER {ber:.2e} at Eb/N0 = {ebn0_db} dB (theory {ber_theory:.2e})")
print(f"measured Es/N0 after the matched filter: {snr_db:.2f} dB (expected {ebn0_db + 10 * np.log10(2):.2f} dB)")

sample rate 1.0 MHz, symbol rate 125 kbaud, filter delay 11 symbols
199978 bits compared, 474 errors
BER 2.37e-03 at Eb/N0 = 6.0 dB (theory 2.39e-03)
measured Es/N0 after the matched filter: 9.01 dB (expected 9.01 dB)

Eb/N0が6 dBでのビット誤り率は2.37 × 10⁻³で、理論値は2.39 × 10⁻³です。474個の誤りを数えた結果なので、この差は統計的なばらつきの範囲に十分収まります。整合フィルタ後に測ったシンボルSNRは9.01 dBです。これはEb/N0の6 dBに、1シンボル2ビット分の3.01 dBを加えた値にちょうど一致します。この確認はBERより重要です。雑音振幅、フィルタのスケーリング、サンプリング時点がすべて意図どおりであることを裏付けるからです。

コードのうち4つの点を詳しく見ておきます。

遅延。 各RRCフィルタは89タップで、信号を44サンプル遅らせます。2つ合わせると88サンプル、つまり1シンボル8サンプルでちょうど11シンボルの遅延です。合計がシンボルの整数倍なので、間引きフィルタの出力は追加の位置合わせなしでシンボルのピークに乗ります。受信ビット列は22ビット遅れて始まるため、コードは比較の前に22ビットを読み飛ばします。合計遅延が8の倍数にならないようにタップ数を変えると、整合フィルタはピークの間でサンプリングし、BERが上がります。

コンステレーション。 digital.chunks_to_symbols_bcは、シンボルインデックスiをリストのi番目の点に対応付けます。digital.constellation_decoder_cbはその逆で、最も近い点のインデックスを返します。デコーダのコンステレーションは、同じリストからdigital.constellation_calcdistで作ります。このとき正規化をオフにして、GNU Radioが信号点を再スケーリングしないようにします。リストの順序がそのままビットマッピングになります。ここでは、隣り合う点が1ビットだけ異なるようにグレイ符号化しています。

ビットとシンボル。 blocks.pack_k_bits_bb(2)は、1ビットずつ入った2バイトを、0から3の値を持つ1バイトにまとめます。最初のビットが最上位です。blocks.unpack_k_bits_bb(2)は、デコーダの後でこれを元に戻します。それに応じて、レートは半分になり、また2倍に戻ります。

終了。 雑音源は終わりませんが、フィルタ後の信号が届かなくなると加算器は出力を生成できません。ベクトルソースのデータが尽きると、その下流のブロックはすべて終了し、run()が戻ります。

クラスで描くBER曲線

GRCが生成するコードは、gr.top_blockから派生したクラスにフローグラフを入れ、ブロックを属性として持たせます。次のコードは同じことを手で書き、Eb/N0を0 dBから10 dBまで掃引します。点ごとに新しいグラフを作ります。7 dBまでは200 000 ビットを使い、誤りがまれになるそれより上ではさらに多くのビットを使います。

# gnuradio
import numpy as np
from scipy.special import erfc
from gnuradio import gr, blocks, analog, digital, filter
from gnuradio.filter import firdes

SPS = 8
POINTS = [(1 + 1j) / np.sqrt(2), (-1 + 1j) / np.sqrt(2), (1 - 1j) / np.sqrt(2), (-1 - 1j) / np.sqrt(2)]
NTAPS = 11 * SPS + 1
RRC = np.array(firdes.root_raised_cosine(1.0, SPS, 1.0, 0.35, NTAPS))
RRC /= np.sqrt(np.sum(RRC**2))


class QpskLink(gr.top_block):
    """Bits -> QPSK -> RRC -> AWGN -> matched filter -> decisions -> bits, all headless."""

    def __init__(self, bits: np.ndarray, ebn0_db: float, seed: int):
        gr.top_block.__init__(self, "QPSK link")
        n0 = 0.5 / 10 ** (ebn0_db / 10)                       # Es = 1, Eb = 0.5
        const = digital.constellation_calcdist(POINTS, [], 4, 1, digital.constellation.NO_NORMALIZATION)
        self.src = blocks.vector_source_b(bits.tolist(), False)
        self.pack = blocks.pack_k_bits_bb(2)
        self.mapper = digital.chunks_to_symbols_bc(POINTS)
        self.tx_filt = filter.interp_fir_filter_ccf(SPS, RRC.tolist())
        self.noise = analog.noise_source_c(analog.GR_GAUSSIAN, float(np.sqrt(n0)), seed)
        self.add = blocks.add_cc()
        self.rx_filt = filter.fir_filter_ccf(SPS, RRC.tolist())
        self.decoder = digital.constellation_decoder_cb(const.base())
        self.unpack = blocks.unpack_k_bits_bb(2)
        self.sink = blocks.vector_sink_b()
        self.connect(self.src, self.pack, self.mapper, self.tx_filt, (self.add, 0))
        self.connect(self.noise, (self.add, 1))
        self.connect(self.add, self.rx_filt, self.decoder, self.unpack, self.sink)


rng = np.random.default_rng(9)
skip = 2 * (NTAPS - 1) // SPS                                  # filter delay in bits
print("Eb/N0 dB, bits, errors, BER sim, BER theory")
for k, ebn0_db in enumerate(range(0, 11)):
    n_bits = 200_000 if ebn0_db < 8 else 2_000_000 if ebn0_db < 10 else 8_000_000   # more bits where errors are rare
    bits = rng.integers(0, 2, n_bits).astype(np.uint8)
    tb = QpskLink(bits, ebn0_db, seed=100 + k)
    tb.run()
    rx = np.array(tb.sink.data(), dtype=np.uint8)
    n = min(rx.size - skip, n_bits)
    errors = int(np.sum(rx[skip : skip + n] != bits[:n]))
    theory = 0.5 * erfc(np.sqrt(10 ** (ebn0_db / 10)))
    print(f"{ebn0_db}, {n}, {errors}, {errors / n:.2e}, {theory:.2e}")

Eb/N0 dB, bits, errors, BER sim, BER theory
0, 199978, 15852, 7.93e-02, 7.86e-02
1, 199978, 11152, 5.58e-02, 5.63e-02
2, 199978, 7554, 3.78e-02, 3.75e-02
3, 199978, 4517, 2.26e-02, 2.29e-02
4, 199978, 2495, 1.25e-02, 1.25e-02
5, 199978, 1177, 5.89e-03, 5.95e-03
6, 199978, 520, 2.60e-03, 2.39e-03
7, 199978, 162, 8.10e-04, 7.73e-04
8, 1999978, 376, 1.88e-04, 1.91e-04
9, 1999978, 75, 3.75e-05, 3.36e-05
10, 7999978, 29, 3.63e-06, 3.87e-06

シミュレーションは全範囲で理論に沿っています。10 dBでは800万ビットで29個の誤りが出て、BERは3.63 × 10⁻⁶、理論値は3.87 × 10⁻⁶でした。このようにまれな事象の点には多くのビットが必要です。誤りを100個数えると、推定値の相対標準偏差は約1/√100 = 10%になります。少なくともその数を目安にしてください。

ブロックは信頼する前に調べる

GNU Radioのchannels.channel_modelは、雑音、周波数オフセット、タイミングオフセット（epsilon）、マルチパスFIR（taps）を一つのブロックにまとめたものです。GNU Radio APIリファレンスでは、基本的なチャネルシミュレータと説明されています。便利ですが、雑音を加える以上のことをします。なじみのないブロックをBER試験に入れる前に、インパルスと無音を入力してください。

# gnuradio
import numpy as np
from gnuradio import gr, blocks, channels

# Probe channels.channel_model with an impulse and with silence before trusting it in a link.
x = np.zeros(64, dtype=np.complex64)
x[10] = 1.0
tb = gr.top_block()
src = blocks.vector_source_c(x.tolist(), False)
chan = channels.channel_model(noise_voltage=0.0, frequency_offset=0.0, epsilon=1.0, taps=[1.0], noise_seed=0)
sink = blocks.vector_sink_c()
tb.connect(src, chan, sink)
tb.run()
y = np.array(sink.data())
print(f"in: {x.size} samples, impulse at index 10; out: {y.size} samples, peak at index {int(np.argmax(np.abs(y)))}")

tb = gr.top_block()
src = blocks.vector_source_c([0j] * 100_000, False)
chan = channels.channel_model(noise_voltage=0.5, frequency_offset=0.0, epsilon=1.0, taps=[1.0], noise_seed=7)
sink = blocks.vector_sink_c()
tb.connect(src, chan, sink)
tb.run()
w = np.array(sink.data())
print(f"noise_voltage 0.5: complex noise variance {np.var(w):.3f}, per rail {np.var(w.real):.3f}")

in: 64 samples, impulse at index 10; out: 56 samples, peak at index 7
noise_voltage 0.5: complex noise variance 0.248, per rail 0.126

2つのことが分かります。まず、GNU Radio 3.10.12では、このブロックは64サンプルの入力に対して56サンプルを返し、インパルスはインデックス10から7に移りました。タイミングオフセットを実装する内部のリサンプラが、ストリームをずらし、いくつかのサンプルを保持するためです。フィルタを使った本レッスンのリンクでは、この3サンプルのずれによって整合フィルタがシンボルのピークから外れます。そのため、ストリームの位置を合わせ直す必要があります。たとえば3サンプルのblocks.delayを入れて、ずれを打ち消します。次に、noise_voltageは複素雑音の標準偏差です。0.5を指定すると全分散は0.248となり、IとQに均等に分かれました。最初のフローグラフで代わりにanalog.noise_source_cとblocks.add_ccを使ったのは、タイミングが正確に保たれ、雑音レベルが明示的になるからです。

GNU Radio Companionで見る同じグラフ

GRCはドラッグ＆ドロップ式のエディタで、フローグラフを.grcファイルとして保存します。このファイルはYAMLです。パレットの各ブロックは、GNU Radioに同梱される.block.ymlファイルで記述されています。このファイルには、ブロックのラベル、パラメータ、makeテンプレートが入っています。makeテンプレートは、GRCが生成スクリプトに書き込むPythonの呼び出しです。GNU Radio 3.10.12では、Channel Modelブロックのテンプレートはchannels.channel_model(...)を生成します。Addブロックはblocks.add_vcc(1)を生成します。これは加算器のベクトル版で、ベクトル長は1です。

このレッスンのPython呼び出し	GRCブロックのラベル
blocks.vector_source_b	Vector Source
blocks.pack_k_bits_bb	Pack K Bits
digital.chunks_to_symbols_bc	Chunks to Symbols
filter.interp_fir_filter_ccf	Interpolating FIR Filter
firdes.root_raised_cosine	RRC Filter Taps（変数ブロック）
analog.noise_source_c	Noise Source
blocks.add_cc	Add
filter.fir_filter_ccf	Decimating FIR Filter
digital.constellation_decoder_cb	Constellation Decoder
blocks.unpack_k_bits_bb	Unpack K Bits
blocks.vector_sink_b	Vector Sink
channels.channel_model	Channel Model

すべてのGRCフローグラフにはOptionsブロックがあります。そのgenerate_options設定で、GRCが書き出すスクリプトの種類を選びます。GNU Radio 3.10.12に同梱されるサンプルフローグラフでは、値はグラフィカルアプリケーション用のqt_gui、ヘッドレススクリプト用のno_gui、階層ブロック用のhbです。samp_rateなどの変数は、生成されるクラスの属性になります。実験台用にGRCで作ったフローグラフは、通常、スペクトラムやコンステレーションの表示といったGUIシンクで終わります。自動テストでは、それらを入れずに本レッスンのようにベクトルシンクで終えるか、no_guiで生成してください。

多くのGRCシミュレーションには、Throttleブロックも見られます。APIリファレンスによれば、これはほかにレートを制限するブロックがないGUIアプリケーション専用であり、レートを正確には制御しません。速く終わるべきヘッドレステストには入れません。ソフトウェア無線（SDR）とやり取りするフローグラフでは、無線機のサンプルクロックがペースを決めるので、まったく必要ありません。次のレッスンでは、そのような無線機からサンプルを記録し、標準的なファイル形式で保存します。

演習

1. ntapsを11 * sps（偶数）に変えてください。遅延とBERに何が起きるかを説明し、タップを変えずに位置合わせを直してください。
2. 最初のコードブロックの雑音源と加算器をchannels.channel_modelに置き換え、その後にblocks.delay(gr.sizeof_gr_complex, 3)を追加してください。BERが理論値に戻ることを確認してください。
3. channels.channel_modelに、正規化したfrequency_offsetとして0.001を与えてください。整合フィルタの出力をベクトルシンクで集め、コンステレーションが時間とともにどう変化するかを説明してください。

まとめ

• GNU Radioのフローグラフは、gr.top_blockの中で型付きのストリームによって接続されたブロックです。run()はすべてのブロックが完了すると戻るので、シミュレーションには有限のソースかblocks.headが必要です。
• サンプルレートは自分で管理する数値です。インターポレータとデシメータがレートを変えます。シミュレーションでは、グラフはCPUが許す限りの速さで動きます。
• 本レッスンのQPSKリンクでは、Eb/N0 6 dBで2.37 × 10⁻³を測定し、理論値は2.39 × 10⁻³でした。測定したシンボルSNRは、設定値と0.01 dBの精度で一致しました。
• なじみのないブロックは、インパルスと無音で調べてください。channels.channel_modelはストリームを3サンプルずらし、雑音の大きさを複素数の標準偏差で指定します。
• GRCはブロックごとのYAMLテンプレートからPythonを書き出します。そのため、どのGRCブロックにも、直接使えるPythonの呼び出しがあります。

さらに読む

• GNU Radio Manual and C++ API Reference: gr::top_blockやgr::channels::channel_modelを含む、全ブロックの公式リファレンスです。本稿執筆時点で、公開されているマニュアルはGNU Radio 3.10.9.1から生成されています。
• GNU Radio v3.10.12.0 release: このページのすべてのブロックで使用したリリースです。
• J. G. Proakis and M. Salehi, Digital Communications, 5th edition, McGraw-Hill, 2008: シミュレーションの照合に使ったQPSKの誤り率については、章「Optimum Receivers for AWGN Channels」を参照してください。