長波JJY受信専用のダイレクトサンプリング方式のソフトウェアラジオ(Software Defined Radio)を作りました． 入力した電波をマイコンのADコンバータにより直接サンプリングして，その後の信号処理はすべてマイコン上で行います． STM32F103C8T6を搭載した2ドル弱程度(最近半導体不足のため値上がりしていますが)のBlue Pillと呼ばれるマイコンボードを使い，Arduino IDEでソフトの開発を行いました． しばらく前にSTM32F4マイコンを使ったソフトウェアラジオを作りましたが，長波帯のCW受信に適した仕様ではなかったので，JJYを受信することに特化したラジオを作りました．
仕様は次のようになります:
今回使用したSTM32F103はSTM32マイコンの中では非力な部類に入りますが，前回STM32F4でSDRを作成した際に計算速度には余裕があったので，安価に出回っているBlue PillボードでSDRを作ってみたいと思いこれを選択しました． 計算速度よりもRAM容量の方がむしろ問題になり，現状のコードではグローバル変数によるRAMの使用率は91%です． このマイコンにはFPUがないため固定小数点演算で信号処理を行います． STM32F4には固定小数点演算を効率的に行うのに便利なDSP命令が搭載されていますが，STM32F103にはそれもありません． ただし固定小数点演算に便利なSSAT命令という桁あふれを扱う命令はあるので多用しています． コード中のコメントでq15(q0.15)やq1.14という表記がありますが，前者は実際の値を32768倍したものを，後者は16384倍したものをint16_t型で保持していることを意味します．
STM32F103C8T6のクロック周波数は最大で72MHzですが，112MHzにオーバークロックして使いました． これはADコンバータのサンプリング周波数を最大にする(ADCCLKを12MHzにする)には，分周比の都合でクロック周波数を56MHzか112MHzにする必要があるからです． Web上では128MHzでの動作例もいくつかあり，これまでのところオーバークロックによる問題は起きていません． サンプリング周波数をキリの良い1Mspsにしたかったのでオーバークロックしましたが，それにこだわらなければ処理時間にはまだ余裕があるので定格内のクロックで動くと思います．
マイコンのポートの使用状況は次のようになっています:
| PA0 | ADC入力 |
| PA1 | PWM出力 |
| PA8 | クロック出力(クロック周波数測定用) |
| PA9 | UARTデバッグ出力 |
| PB2 | 周波数切り替えスイッチ |
| PB10 | I/Q信号出力(モニター用) |
受信周波数(40kHz/60kHz)を切り替えるスイッチはPB2を使っていて，これはマイコンボード上でBOOT 1のジャンパースイッチに使われていますが，このジャンパースイッチはフラッシュ書き込み時以外(BOOT 0がLの時)は使われないため，ピン端子にスライドスイッチを取り付けて使いました． マイコンのクロック周波数測定用に，PA8にクロック源の8MHzを出力しています．今回使ったボードでは8,000,634Hzだったことから，60kHzの信号は実際には59,995Hzとして読み込まれていることになりますが，無視できる誤差だと思います．
Blue Pillボードには互換チップや偽チップを使ったものが出回っていることが知られていますが，今回もその問題にはまって時間を費やしてしまいました． 正確な50kHzの矩形波を入力しているのに，ADCの結果計算される振幅が計算よりも短いという問題が起きました． 色々試したところ，STM32F103のADCは(サンプリングタイム + 12.5)のADCCLKで変換するはずが，(サンプリングタイム + 13.5)で計算すると実測結果と一致することが分かりました． 試しに別の店で購入したBlue Pillボードに交換すると問題が解決しました．このような検証が面倒な微妙な非互換性は勘弁して欲しいと思います． また，問題があったボードでは112MHzにオーバークロックすると動きませんでした． 以前問題が起きたBlue Pillボードはチップの表記自体がSTM32とは異なっている互換チップでしたが，今回の偽チップはマーキングを見ても本物と区別がつきませんでした(下の写真参照)． もっとも，今回正常に動いたチップもよくできた偽物という可能性もありますが． ちなみにこれらの3種類のボードはAliExpressで異なる業者から購入したものですが，商品説明にはどれも"STM32F103C8T6"と明記してありました． 問題のあるBlue Pillボードでも，サンプリング周波数を書き換えるだけで本プログラムを動作させることは可能ですが，ADCのアナログ面の特性等にどこまで互換性があるのか不明です．
CWの受信だけを考えるのでダイレクトコンバージョン方式としました． つまり信号周波数を音声(CWピッチ)周波数に直接変換します． 単純なダイレクトコンバージョン受信機では両側のサイドバンドが混じりノイズが増加してしまいますが，本機では直交ミキサーを使い片側の成分だけを取り出します． ADCでは1Mspsで信号を取り込みますが，その後サンプリング周波数を256分の1に落として3906.25Hzとします． CWピッチ信号を保持できる範囲内でサンプリング周波数をできるだけ低くした方がBPFを設計しやすく，またデシメーション比が高くなりSNRを改善できます．
今回作成するSDRの構成図は次のとおりです:
アンテナから入力された信号は，高周波アンプで増幅されてマイコンに入力されます． マイコンでは，ADコンバータ(ADC)でデジタル信号に変換した後に直流成分(DC)の除去を行います． そして数値制御発信器(NCO)と直交ミキサー(I/Qミキサー)によりCW信号の周波数をピッチ周波数に変換します． ピッチ周波数より高い信号は不要となるので，CICフィルターとFIRフィルターによりサンプリング周波数を256分の1にします． 次にヒルベルト変換によりイメージ信号を除去します． そしてバンドパスフィルター(BPF)でピッチ信号以外の不要な信号を除去して，自動利得制御(AGC)で信号レベルを調節します． アナログ信号出力時のスプリアスを減らすためサンプリング周波数を8倍に増やして(インターポレーション)，PWMでマイコンから出力します． 最後に高調波除去のためのフィルターを通してイヤホンから出力します．
ADCからバッファに取り込んで一度に処理するデータ量が多い方が，プログラム実行時のオーバーヘッドを削減できますが，入力から出力までの遅延が増加して使用するRAMも増大します． 今回はRAMのサイズに合わせて2,048サンプルをADCから一度に取得します．
アンテナから入力した信号を増幅しますが，高周波といっても40〜60kHzしかありません． そのためオペアンプのNJM2122Dを使いました． これは入力換算雑音が1.5nV√Hzという超低ノイズのアンプでGB積も12MHzありますが，ゲインが30dB以上で安定するというのがやや面倒です． ここでは100倍(40dB)のゲインとしています． 単電源で動かすため回路がシンプルになる反転増幅回路を使い，120kHzのLPF特性も持たせてあります． 最低電源電圧は4Vなのでオペアンプの電源には5Vを使用し，出力は3.3VフルスケールのADC入力のために1.65Vのオフセットを与えます． デュアルオペアンプなので，余った方は簡易的な(あまり推奨されない)方法で処理しています．
回路図と完成した基板は次のとおりです(この回路には高周波アンプ部だけではなく出力フィルター部も含まれています)．
マイコン内蔵ADCにより1Mspsで入力信号をサンプリングします． STM32F103C8T6は，2つのADCを交互に動作させるインターリーブモードを使うことで最大2Mspsで信号を取り込めますが，ここではシンプルに1つのADCだけを使っています． ところで，インターリーブモードを使った場合は入力データ(ハーフワード)の順序がビッグエンディアンのように逆になるので少し面倒になります． バッファを2つ用意することで(ダブルバッファリング)，片方のバッファの処理中に別のバッファにDMA転送でデータを読み込みます． 120kHz以上のサンプリング周波数であれば60kHzの信号を処理できますが，デシメーション処理によるSNRの改善も考慮して1Mspsという高めのサンプリング周波数を選んでいます．
ADCでデータを取り込んだ後は，直流成分を除去します．
直交ミキサーにより目的のCW信号(40kHz/60kHz)をピッチ周波数に直接変換します． 今回CWのピッチ周波数は40kHz受信時は937.5Hz，60kHz受信時は918.0Hzとしたため，40kHzおよび60kHzからこれらのピッチ周波数を引いた差分の周波数の信号をNCOで生成します． このように中途半端なピッチ周波数に決めた理由は，NCO結果を保存するテーブルのサイズに関係します． 今回実装した直交ミキサーでは，計算を効率的に行うためにNCOテーブルのサイズをADCバッファのサイズと同じにしています． そのため，NCOで生成する信号の周期の整数倍がNCOテーブルのサイズと等しくなる必要があります． 具体的には，fを生成したい目標の周波数，f'を実際に生成される周波数とすると，
f' = ⌊ f * 2,048 / 1,000,000 + 0.5 ⌋ * 1,000,000 / 2,048
となりますが，ここで⌊⌋は床関数，2,048はNCOテーブルのサイズ，1,000,000はサンプリング周波数です． つまり，NCOの周波数分解能を上げるには，NCOテーブルのサイズを大きくするかサンプリング周波数を小さくする必要があります． 前述のピッチ周波数は，目標となるピッチ周波数を800Hzとした場合にこの式で得られる値です．
デシメーション処理によりサンプリング周波数を低い周波数に変換します． 直交ミキサーにより得られたCWのピッチ周波数が1kHz弱で，これより高い周波数はもう不要です． CICフィルターを使って32分の1に落としてから，FIRフィルターを使ってさらに8分の1に落とします． デシメーション処理を行う理由の一つはデータ量を削減して計算効率を上げることですが，SNRを上げることもできます． 具体的には今回の場合，NをADCの有効ビット(12)，Rをデシメーション比(256)として，SNR=6.02 N + 1.76 + 10 log10 R=98.1dBとなります．
ヒルベルト変換により負の周波数成分を除去して，イメージ成分を取り除いた信号を取り出します．
なお，デバッグ用にヒルベルト変換前のI/Q信号をシリアルインターフェースで外部に出力しています． Spectrum LabはIQ12というフォーマットでシリアルポートからI/Q信号を入力する機能を持っているので，マイコンに取り込んだ信号のスペクトラムやウォーターフォール出力を簡単に見ることができます． IQ12は1サンプルが4バイトのフォーマットであり，startビットとstopビットが1ビットで周波数が3906.25Hzの場合は通信速度は156,250 baud以上必要なので，500,000 baudで通信しています． メモリ上に用意したデータをDMAでペリフェラルに送ります．
ピッチ周波数を中心として幅が50Hzのバンドパスフィルターで不要な信号を除去します． フィルター幅が狭いので，40kHと60kHzのそれぞれのピッチ周波数でフィルターを分けています． 急峻な特性を得るためにbiquad (双2次)フィルターを3段重ねています． biquadフィルターの実装方法には直接型Iと直接型IIがあり，最初はシンプルな直接型IIを試しましたがシミュレーション結果がおかしくなりました． 固定小数点演算で直接型IIを使うと誤差が大きくなるようなので，直接型Iに変えたところ期待通りの動作をしました．
このバンドパスフィルター，FIRデシメーションフィルター，ヒルベルト変換フィルター，及び後述のFIRインターポレーションフィルターはPythonのscipyパッケージを使って設計しました(使用したスクリプト)． このバンドパスフィルターの特性は下のようになりました．
設定されたアタック時間，ホールド時間，リリース時間に基づいて，出力する信号レベルを制御します．
この時点で信号のサンプリング周波数は約4kHzであり，これは可聴周波数帯域です． このままアナログ信号に変換して出力するとノイズの原因になります． そのため，出力前にサンプリング周波数を8倍にして可聴周波数帯域よりも高くします．
なおFIRデシメーションとこのFIRインターポレーションでは比率が同じ(それぞれ8分の1と8倍)であり，FIRフィルターとしては同じなので係数を共有できます． しかしながらシミュレーション結果では，32タップFIRフィルターだとインターポレーションではスプリアスが残るので64タップとしました． 8倍のインターポレーション処理では，1つの信号を入力するたびに7個のゼロを挿入して出力して，それをフィルタリングします．ゼロの値に対する係数の乗算は不要なので計算を効率化できます．
今回使用するSTM32F103C8T6はDACを搭載していないため，音声信号の出力にはPWMを使いました．分解能は8ビットとしています． タイマーのDMA転送でダブルバッファリングされたメモリからデータを転送します． DACと比べて，PWMではアンプを使わずに直接イヤホンを駆動できたので外付け回路がシンプルになりました．
PWM信号からアナログ信号を得るために簡単なRCローパスフィルターを通した後で，直流カット用コンデンサーを介してイヤホンへ信号が送られます． 回路図は，前述の高周波増幅のものに含まれているのでそちらを参照して下さい． 回路図ではフィルターのカットオフ周波数が2kHzではなく4kHzとなっているのは，最初の設計ではデシメーション後の周波数を約8kHzとしていたからです． インターポレーション処理によりカットオフ周波数より十分高いサンプリング周波数に変換されているので，フィルターは簡単なもので済みます．
作成したプログラムはページ末尾にリンクのあるGitHub上に置いてあります． 前回SDRを作成した時と同様に，DSP処理部分はシンプルでバグの入りにくい実装とアルゴリズムを最適化した実装を用意して，PC上でシミュレーションしながらデバッグしました． 各処理内容の実行時間を測定した結果は下記の表のようになりました(単位はマイクロ秒)． 最適化してもあまり速くならなかった処理や，元々速くて最適化の必要がない処理については最適化版は用意していません(表の中では括弧で示してあります)．
| 処理 | 最適化前 | 最適化後 |
| 直流成分除去 | 220 | (220) |
| 直交ミキサーとCICデシメーション | 1072 | 608 |
| FIRデシメーション | 987 | 42 |
| I/Q信号外部出力 | 4 | (4) |
| ヒルベルト変換 | 97 | 17 |
| バンドパスフィルター | 14 | (14) |
| 自動利得制御 | 4 | (4) |
| インターポレーション | 1493 | 25 |
| 合計 | 3891 | 936 |
ADCの結果を1Mspsで2,048サンプル取り込むので2,048usの処理時間がありますが，最適化後の処理は十分にそれに収まっています． 今回の構成には無駄があって，直交ミキサーの前に周波数を8分の1に減らすデシメーション処理を行えば計算量は減らせます．
同調型のループアンテナをつけて実際にJJYを受信してみました． 40kHzは室内でも強力に受信できましたが，60kHzは以前受信したときと同様に近所の公園でないと受信が困難でした． アップコンバーターとワンセグ受信機を組み合わせてPC上で受信した場合と比べて，今回作成した受信機の方が信号が聞きやすいと感じました． サイドバンド除去と狭帯域BPFが効いたのかもしれませんが，定量的な評価はしていません． きれいな信号が受信できているので，JJYのタイムコードをデコードして時刻情報を取り出すことも可能だと思いますが，今回はJJYの受信だけが目的だったのでそこまではやっていません．
今回ソフトウェアラジオの他に，JJY受信用のアンテナも新たに一つ作成しました． 前回作成したループアンテナは巻線の数があまり多くなかったので，もっと巻線数を増やしたいと思っていました． サイズは420x320mmの長方形で，0.2mmの100mのワイヤーを巻いたので巻数は67回ほどあると思います． 前回のループアンテナもそうでしたが，外に持ち出せるように大きめのリュックサックに入るサイズにしてあります． フレームをすべて3Dプリンターで作ろうとしたところ剛性がなかったので，9mm角の木材を3Dプリンターで作ったジョイント兼ワイヤーのガイドで固定して作りました． 前回同様にFRMSで同調用コンデンサーの値を求めました． 今回作成したSDRで使用したところ，40kHzも60kHzも問題なく受信できました． しかしPC用のSDRでAGCを切って使用してみたところ，以前作ったループアンテナに比べて信号が弱くやや静かなアンテナでした．
2021年10月にJJYのはがね山送信所(60kHz)が開局20周年を迎えて記念のベリカードを発行しているため，受信して受信報告書を送りました． 40kHzのおおたかどや山送信所の受信報告書もついでに出しましたが，下記のベリカードが届きました． 今回JJY受信専用ラジオを作った大きな目的は，この受信報告書を書くことでした． 以前書きましたが，10年前にはがね山送信所が開局10周年を迎えたときにはどうしても受信できず悔しい思いをしました． 今回は受信に成功し，また(ソフトウェアも)完全自作のラジオで受信するという目標も達成できました． 10年前はアンテナや受信機をいくつか試しても受信できませんでしたが，結局マンションの自室というロケーションが大きな問題だったと思います． 今回は自作SDRの他にも秘密兵器がありました． それは最近買い換えたPCに付属のサウンドカードで，24bit 192kHzのADCを持つRealtek ALC887チップを使っているため，60kHzの電波も直接サウンドカードに入力してPCで受信できます． しかし，これでもやはり室内では受信できませんでした．
2つのBlue Pillボード(左: 問題のあるボード, 右: 正常品)
回路図
完成した基板(表)
完成した基板(裏)
バンドパスフィルターの特性
配線後の基板
ケースに収めた状態
イヤホン出力をSpectrum Labで表示したもの(40kHz受信時)
新しいループアンテナ
40kHzに同調させた時
60kHzに同調させた時
JJYのベリカード