合成IQで試験するヘッドレスGNU Radio FM受信機

作るもの

このプロジェクトでは、GNU Radioで広帯域FM受信機を作り、無線機を使わずに試験します。NumPyスクリプトで1 kHzのトーンを周波数変調し、5 kHzの搬送波オフセットと雑音を加え、結果をSigMF記録として書き出します。SigMF記録は、生のcomplex64データファイルとJSONメタデータファイルの組です。ヘッドレスのGNU Radioフローグラフがこのファイルを1回だけ読み、フィルタをかけ、直交弁別器で復調し、48 kHzのオーディオにデシメーションします。続いてスクリプトがオーディオをコードで測定します。測定項目はトーン周波数、ピーク偏移、搬送波オフセット、オーディオSNRです。最後にPASSまたはFAILを出力します。

自作の信号で試験するので、正解は事前にわかっています。筆者らの実行では、受信機はトーンを1000.000 Hz、偏移を60.14 kHz（送信は60 kHz）、搬送波オフセットを5.00 kHzと復元しました。続いて搬送波対雑音比（CNR）を掃引すると、FMスレッショルドが現れます。CNRが約10 dBを超える範囲では、CNRが1 dB上がるごとにオーディオSNRも1 dB上がります。それを下回ると、オーディオSNRは急速に低下します。

「実機での受信」の節では、ファイルの代わりにRTL-SDRと放送局を使います。この部分には無線機の接続が必要です。

背景は、レッスン「初めてのGNU Radioフローグラフ」「IQの記録と読み取り」「複素ベースバンドとIQ」で扱っています。インストール方法はGNU Radioのプラットフォームページにあります。

必要なもの

• PythonバインディングつきのGNU Radio 3.10とNumPy。筆者らは、Apple SiliconのMac mini上で、conda-forgeのPython 3.12、GNU Radio 3.10.12、NumPy 2.5.4を使って実行しました。フローグラフはblocks・filter・analogモジュールだけを使い、GUIブロックは使いません。
• 短い記録2つを保存するための、作業ディレクトリの約4 MBの空き容量。
• 任意のハードウェア節には、RTL-SDR受信機、gr-soapyブロック（GNU Radio 3.10に含まれる）、SoapyRTLSDRドライバモジュール。

設計

試験信号

メッセージは、フルスケールの0.8倍の1 kHz正弦波です。フルスケールは75 kHz偏移とします。これは勧告ITU-R BS.450-4が西欧と米国のFM音声放送に対して定める最大値なので、ピーク偏移は60 kHzになります。瞬時周波数は、60 kHz × sin(2π × 1000 t)に5 kHzのオフセットを加えたものです。このオフセットは、局からわずかにずれて同調した受信機を模擬します。位相は瞬時周波数の累積和で、IQサンプルはexp(jφ)、サンプルレートは240 kHzです。

カーソンの法則によれば、FM電力の大部分は2 × (60 + 1) = 122 kHzの帯域内に収まります。したがって信号は、オフセット分の余裕を残して240 kHzのサンプルレートに収まります。CNRは複素白色雑音で設定し、240 kHz全体で測ります。受信機のチャネルフィルタが100 kHzより上の雑音を除去するので、フィルタ内のCNRは表示値よりいくらか高くなります。

SigMF

記録は、同じベース名を持つ2つのファイルです。tone_fm.sigmf-dataには、生のサンプルがインターリーブされたリトルエンディアンの32 ビット浮動小数点数として入っています。これはnumpy.complex64.tofileが書き出す形式であり、GNU Radioのgr.sizeof_gr_complexアイテムの形式でもあります。tone_fm.sigmf-metaはJSONで、トップレベルのキーはglobal・captures・annotationsの3つです。SigMFスキーマは、global内のcore:datatypeとcore:versionを必須としています。データ型にはcf32_le、バージョンには本稿執筆時点で最新のSigMFリリースである1.2.6を指定します。さらに、サンプルレート、説明、サンプル0から始まるキャプチャセグメント1つ、トーンにラベルを付けるアノテーション1つを加えます。受信機は定数を信用せず、サンプルレートをメタデータから読みます。全フィールドはSigMF仕様に載っています。

フローグラフ

blocks.file_source(gr.sizeof_gr_complex, path, False)はデータファイルを1回だけ読みます。repeat=Falseの場合、ファイルの終わりでストリームの終了を通知するので、tb.run()が戻ります。29タップのローパスフィルタ（カットオフ100 kHz、遷移幅20 kHz）が帯域端付近の雑音を削ります。

analog.quadrature_demod_cf(gain)は、x[n] × conj(x[n − 1])の偏角、つまりサンプル間の位相ステップを計算し、利得を掛けます。位相ステップ2π f / fsラジアンは瞬時周波数fに対応するので、利得をfs / (2π × 75 kHz)にすると、75 kHz偏移で出力が1.0になります。GNU Radio自身のanalog.wfm_rcvブロックも同じスケーリングを使います。このとき出力のDCレベルは、搬送波オフセットを75 kHzで割った値になります。

15 kHzのローパスフィルタを持ち、係数5でデシメーションするfilter.fir_filter_fffが48 kHzのオーディオを生成します。それをblocks.vector_sink_fが収集します。送信側でプリエンファシスをかけていないので、受信側でもデエンファシスはかけません。

トーンの測定

フィルタの立ち上がり分として最初の50 msを捨てた後、スクリプトは平均（オフセット）を除去し、Hann窓をかけ、最大のFFTビンを探します。ビン間隔は1.05 Hzです。ピークビンとその両隣の2ビンの対数に放物線を当てはめ、ピーク位置をビン幅のごく一部の精度で求めます。偏移は、オーディオのRMSの√2倍に75 kHzを掛けた値です。オーディオSNRは、トーンから20 Hz以内の電力を、20 Hz〜15 kHz帯域の残りの電力で割った値です。

コード

これをfm_receiver.pyとして保存してください。1つのブロックで完結しています。先頭行の# gnuradioは当サイトのコードチェッカー用のマーカーで、ファイルを実行するときには何の効果もありません。

# gnuradio
# fm_receiver.py: make an FM test signal, save it as SigMF, demodulate it with GNU Radio,
# and check the recovered audio tone. Complex samples are numpy.complex64 at 240 kHz.
import json
import math
from pathlib import Path

import numpy as np
from gnuradio import analog, blocks, filter, gr

FS = 240_000            # complex baseband sample rate (Hz)
DURATION = 1.0          # seconds of IQ to generate
TONE_HZ = 1_000.0       # audio test tone
MAX_DEV = 75_000.0      # peak frequency deviation for a full-scale tone (Hz)
TONE_LEVEL = 0.8        # tone amplitude relative to full scale, so the deviation is 60 kHz
CARRIER_OFFSET = 5_000  # transmitter tuned 5 kHz away from the receiver (Hz)
CNR_DB = 20.0           # carrier-to-noise ratio over the full 240 kHz
AUDIO_DECIM = 5         # 240 kHz / 5 = 48 kHz audio
BASE = Path("tone_fm")


def make_recording(base: Path, cnr_db: float, seed: int = 7) -> int:
    """Frequency-modulate a tone, add a carrier offset and noise, write SigMF data and metadata."""
    rng = np.random.default_rng(seed)
    n = int(FS * DURATION)
    t = np.arange(n) / FS
    message = TONE_LEVEL * np.sin(2 * np.pi * TONE_HZ * t)
    inst_freq = MAX_DEV * message + CARRIER_OFFSET                   # Hz
    phase = 2 * np.pi * np.cumsum(inst_freq) / FS
    noise_std = math.sqrt(10 ** (-cnr_db / 10) / 2)
    iq = np.exp(1j * phase) + noise_std * (rng.standard_normal(n) + 1j * rng.standard_normal(n))
    iq.astype(np.complex64).tofile(base.with_suffix(".sigmf-data"))
    meta = {
        "global": {
            "core:datatype": "cf32_le",
            "core:sample_rate": FS,
            "core:version": "1.2.6",
            "core:description": f"Synthetic FM: 1 kHz tone, 60 kHz deviation, 5 kHz offset, {cnr_db:g} dB CNR",
            "core:recorder": "fm_receiver.py (numpy)",
        },
        "captures": [{"core:sample_start": 0}],
        "annotations": [{"core:sample_start": 0, "core:sample_count": n, "core:label": "1 kHz test tone"}],
    }
    base.with_suffix(".sigmf-meta").write_text(json.dumps(meta, indent=2))
    return n


def demodulate(base: Path) -> tuple[np.ndarray, float]:
    """Read the SigMF metadata, then run a finite GNU Radio flowgraph over the data file."""
    meta = json.loads(base.with_suffix(".sigmf-meta").read_text())
    assert meta["global"]["core:datatype"] == "cf32_le"
    fs_in = float(meta["global"]["core:sample_rate"])
    tb = gr.top_block()
    src = blocks.file_source(gr.sizeof_gr_complex, str(base.with_suffix(".sigmf-data")), False)
    chan = filter.fir_filter_ccf(1, filter.firdes.low_pass(1.0, fs_in, 100e3, 20e3))   # trim noise near fs/2
    demod = analog.quadrature_demod_cf(fs_in / (2 * math.pi * MAX_DEV))              # 1.0 = 75 kHz deviation
    audio = filter.fir_filter_fff(AUDIO_DECIM, filter.firdes.low_pass(1.0, fs_in, 15e3, 3e3))
    sink = blocks.vector_sink_f()
    tb.connect(src, chan, demod, audio, sink)
    tb.run()                                    # returns when file_source reaches the end of the file
    return np.array(sink.data(), dtype=np.float64), fs_in / AUDIO_DECIM


def analyse(y: np.ndarray, audio_rate: float) -> dict[str, float]:
    """Tone frequency (Hann window, parabolic peak), peak deviation, carrier offset and audio SNR."""
    y = y[int(0.05 * audio_rate):]                                   # drop 50 ms of filter start-up
    dc = float(np.mean(y))
    ac = y - dc
    spec = np.abs(np.fft.rfft(ac * np.hanning(len(ac)))) ** 2
    freqs = np.fft.rfftfreq(len(ac), 1 / audio_rate)
    k = int(np.argmax(spec[1:])) + 1
    a, b, c = np.log(spec[k - 1:k + 2])
    f_est = (k + 0.5 * (a - c) / (a - 2 * b + c)) * audio_rate / len(ac)
    tone = np.abs(freqs - f_est) < 20
    rest = ~tone & (freqs > 20) & (freqs < 15e3)
    return {"freq": f_est, "bin": audio_rate / len(ac), "dev": math.sqrt(2) * np.std(ac) * MAX_DEV,
            "offset": dc * MAX_DEV, "snr": 10 * np.log10(spec[tone].sum() / spec[rest].sum())}


n = make_recording(BASE, CNR_DB)
print(f"wrote {n} samples ({BASE.with_suffix('.sigmf-data').stat().st_size} bytes) and {BASE}.sigmf-meta")
y, audio_rate = demodulate(BASE)
print(f"flowgraph done: {len(y)} audio samples at {audio_rate:.0f} Hz")
r = analyse(y, audio_rate)
print(f"tone frequency: {r['freq']:.3f} Hz (FFT bin spacing {r['bin']:.2f} Hz)")
print(f"peak deviation: {r['dev'] / 1e3:.2f} kHz, carrier offset from DC: {r['offset'] / 1e3:.2f} kHz")
print(f"audio SNR in 20 Hz to 15 kHz: {r['snr']:.1f} dB")
ok = abs(r["freq"] - TONE_HZ) < 1.0 and abs(r["dev"] - TONE_LEVEL * MAX_DEV) < 1e3
print("PASS" if ok else "FAIL")

print("CNR (dB)  audio SNR (dB)  tone (Hz)  deviation (kHz)")
for cnr in range(0, 31, 2):
    make_recording(Path("sweep"), cnr, seed=100 + cnr)
    r = analyse(*demodulate(Path("sweep")))
    print(f"{cnr:>8}  {r['snr']:>14.1f}  {r['freq']:>9.2f}  {r['dev'] / 1e3:>15.2f}")

最後の掃引では、同じ3つの関数を再利用し、点ごとに異なるCNRと異なる雑音シードを与えます。

実行

conda create -n fmrx -c conda-forge gnuradio-core numpy
conda activate fmrx
python fm_receiver.py

16点の掃引を含む全体の実行時間は、筆者らのマシンで1 秒未満です。作業ディレクトリには、tone_fm.sigmf-data、tone_fm.sigmf-meta、最後の掃引の記録が残ります。各データファイルは1,920,000 バイトで、8 バイトのサンプル240,000個分です。

テスト

筆者らの実行での出力は次のとおりです。

wrote 240000 samples (1920000 bytes) and tone_fm.sigmf-meta
flowgraph done: 48000 audio samples at 48000 Hz
tone frequency: 1000.000 Hz (FFT bin spacing 1.05 Hz)
peak deviation: 60.14 kHz, carrier offset from DC: 5.00 kHz
audio SNR in 20 Hz to 15 kHz: 46.3 dB
PASS
CNR (dB)  audio SNR (dB)  tone (Hz)  deviation (kHz)
       0             5.2    1000.00            48.34
       2             9.6    1000.00            54.04
       4            15.2    1000.00            57.99
       6            22.1    1000.00            59.71
       8            31.4    1000.00            60.12
      10            35.6    1000.00            60.14
      12            38.3    1000.00            60.14
      14            40.4    1000.00            60.14
      16            42.3    1000.00            60.14
      18            44.3    1000.00            60.14
      20            46.4    1000.00            60.14
      22            48.3    1000.00            60.14
      24            50.4    1000.00            60.14
      26            52.3    1000.00            60.14
      28            54.3    1000.00            60.14
      30            56.4    1000.00            60.14

各行からわかることは次のとおりです。

• トーン周波数。 1000.000 Hz。1 kHzの誤差許容ではほとんどのバグが隠れてしまうので、スクリプトは1 Hzより良い精度を要求します。トーンの周波数がずれている場合、たいていはメタデータのサンプルレートと実際のレートが食い違っています。
• ピーク偏移。 送信した60 kHzに対して60.14 kHz。偏移が2π倍ずれている場合は、復調器の利得に問題があります。トーンがきれいなのに偏移が低すぎる場合は、チャネルフィルタがFM側波帯を削っています。
• 搬送波オフセット。 オーディオの平均から読み取った値は5.00 kHz。実際の受信機は、この値を使って同調を補正します。
• オーディオSNR。 CNR 20 dBで46.3 dB。FMは帯域幅とSNRを交換します。信号は約122 kHzを占有しますが、オーディオ帯域は15 kHzしかありません。
• 掃引。 CNR 10 dB〜30 dBでは、オーディオSNRは35.6 dBから56.4 dBへ、ほぼ1 dBあたり1 dBの割合で上昇します。8 dB以下ではより速く低下し、8 dBで31.4 dB、6 dBで22.1 dB、0 dBで5.2 dBになります。0 dBでは、測定偏移も48.34 kHzに下がります。スレッショルドを下回ると、雑音がときどき位相を原点の周りに一周させます。そのたびにオーディオにスパイク、つまりクリックが生じます。

実機での受信

この節には無線機の接続が必要です。以下のコードは、RTL-SDRを接続し、gr-soapyブロックとSoapyRTLSDRドライバモジュールをインストールしないと実行できません。

RTL-SDR Blogによれば、RTL-SDR Blog V4はR828Dチューナを搭載し、500 kHz〜1.766 GHzに同調できます。RTL-SDR BlogのRTL-SDR概要ページには、ネイティブのADC分解能が8 ビット、最大サンプルレートが3.2 MS/s、サンプルを欠落させない最高レートが2.56 MS/sと記載されており、RTL-SDRは送信できないとしています。これらのページにはチャネル数の記載がありません。詳しくはハードウェア表を参照してください。

テストベンチからの変更点は次のとおりです。

• ソース。 soapy.source("driver=rtlsdr", "fc32", 1, "", "bufflen=16384", [""], [""])がファイルソースに置き換わります。引数は、GNU Radio 3.10のソースツリーにあるgr-soapyのRTL-SDRブロックテンプレートに従っています。
• サンプルレート。 無線機は、サンプルを欠落させないレートの範囲内である1.92 MS/sで動作します。チャネルフィルタが8分の1にデシメーションしてテストベンチと同じ240 kHzにするので、以降のチェーンは変わりません。
• 有限の実行。 blocks.headが、5 秒分のサンプルを処理した後でフローグラフを止めます。
• デエンファシス。 放送局はオーディオにプリエンファシスをかけています。ITU-R BS.450-4では、欧州で50 µs、米国で75 µsです。analog.fm_deemph(fs, tau)でこれを元に戻します。デフォルトは75 µsなので、地域に合わせてtauを設定してください。
• ステレオ。 BS.450-4によれば、ステレオ局は19 kHzのパイロットトーンと、その上の帯域にステレオ差信号を加えます。15 kHzのオーディオフィルタが両方を除去し、モノラルのオーディオが残ります。

# needs-hardware
# fm_rtlsdr.py: receive 5 s of a broadcast FM station with an RTL-SDR through gr-soapy,
# demodulate it with the same chain as fm_receiver.py and save the audio as a WAV file.
import math
import wave

import numpy as np
from gnuradio import analog, blocks, filter, gr, soapy

STATION_HZ = 98.5e6     # replace with a station you can receive
SDR_RATE = 1.92e6       # RTL-SDR sample rate, decimated by 8 to 240 kHz
DECIM_RF = 8
FS = SDR_RATE / DECIM_RF
AUDIO_DECIM = 5         # 48 kHz audio
MAX_DEV = 75e3
TAU = 50e-6             # de-emphasis: 50 us in Europe, 75 us in the United States
SECONDS = 5

tb = gr.top_block()
src = soapy.source("driver=rtlsdr", "fc32", 1, "", "bufflen=16384", [""], [""])
src.set_sample_rate(0, SDR_RATE)
src.set_frequency(0, STATION_HZ)
src.set_gain_mode(0, False)
src.set_gain(0, 30.0)
head = blocks.head(gr.sizeof_gr_complex, int(SDR_RATE * SECONDS))
chan = filter.fir_filter_ccf(DECIM_RF, filter.firdes.low_pass(1.0, SDR_RATE, 100e3, 20e3))
demod = analog.quadrature_demod_cf(FS / (2 * math.pi * MAX_DEV))
audio = filter.fir_filter_fff(AUDIO_DECIM, filter.firdes.low_pass(1.0, FS, 15e3, 3e3))
deemph = analog.fm_deemph(FS / AUDIO_DECIM, TAU)
sink = blocks.vector_sink_f()
tb.connect(src, head, chan, demod, audio, deemph, sink)
tb.run()

y = np.array(sink.data())
pcm = np.int16(np.clip(y / (np.max(np.abs(y)) + 1e-12), -1, 1) * 32000)
with wave.open("station.wav", "wb") as w:
    w.setnchannels(1)
    w.setsampwidth(2)
    w.setframerate(int(FS / AUDIO_DECIM))
    w.writeframes(pcm.tobytes())
print(f"saved {len(pcm) / (FS / AUDIO_DECIM):.1f} s of audio to station.wav")

まず利得30 dBで始め、そこから調整してください。強い近隣局では、8 ビットADCは容易にクリップします。受信したオーディオにテストベンチのanalyseを実行すると、局の周波数におけるドングルの水晶による搬送波オフセットがわかります。

放送の受信や、受信内容の録音・共有に関する規則は国ごとに定められており、国によって異なります。録音する前に、規制当局の規則を確認してください（英国ではOfcom、米国ではFCC）。放送番組は著作物なので、録音は個人の範囲にとどめてください。この受信機は何も送信しません。

試してみること

1. wfm_rcvを使う。 復調器、オーディオフィルタ、欠けているデエンファシスをanalog.wfm_rcv(fs_in, AUDIO_DECIM)で置き換えてください。このブロックは75 µsのデエンファシスをかけ、1 kHzトーンをわずかに減衰させるので、測定偏移が変わるはずです。どれだけ変わるか測ってください。
2. オフセットを動かす。 CARRIER_OFFSET = 40_000に設定します。筆者らの実行では、トーンは依然として正確でしたが、CNR 20 dBでのオーディオSNRは46.3 dBから33.3 dBに下がりました。上側のFM側波帯が100 kHzのチャネルフィルタにかかるためです。測定したオフセット分だけ周波数をシフトする周波数変換フィルタ（filter.freq_xlating_fir_filter_ccf）を追加し、比較してください。
3. プリエンファシス。 送信側で50 µsの1次フィルタによりメッセージにプリエンファシスをかけ、受信側にanalog.fm_deemph(audio_rate, 50e-6)を追加します。偏移の測定値が60 kHzに戻ることを確認してください。
4. 2トーン。 低いレベルの3 kHzトーンを追加し、analyseを拡張して両方を報告させてください。チャネルフィルタを狭くしすぎたときに、2 kHzと4 kHzに相互変調積が現れるか調べてください。
5. 狭帯域FM。 MAX_DEV = 5_000、FS = 48_000とし、チャネルフィルタを8 kHz、オーディオフィルタを3 kHzに設定します。これは狭帯域FM信号です。CNR掃引でスレッショルドがどこに移るか確認してください。
6. ライブラリでSigMFを書く。 sigmf Pythonパッケージ（本稿執筆時点でバージョン1.14.0）をインストールし、それを使って記録を書き出します。メタデータファイルに対して、パッケージのsigmf_validateコマンドを実行してください。

まとめ

• 受信機は、自分で生成したIQで試験します。正解がわかっているので、試験を数値的かつ自動にできます。
• SigMFは、生のcf32_leデータファイルとJSONメタデータの組です。サンプルレートは受信機内の定数ではなく、メタデータから読みます。
• repeat=Falseのblocks.file_sourceはファイルの終わりでストリームを終えるので、ヘッドレスのフローグラフは最後まで実行されて戻ります。
• 利得をfs / (2π × 75 kHz)としたquadrature_demod_cfは、瞬時周波数をフルスケール偏移の単位で出力します。その平均が搬送波オフセットです。
• 筆者らの実行では、CNRが約10 dBを超える範囲でオーディオSNRは1 dBあたり1 dB上昇しました。それを下回るとFMスレッショルドに入り、オーディオSNRは急落しました。