Teensy（13）Keiths' SDRのFilterConvオブジェクト

Keiths' SDRのRx信号処理連鎖ブロック図

Keiths' SDR（K7MDL局Mike OM版）のRx信号処理連鎖をコードから読み解いて抜き出したもの（α版）を下記に再録します。

今回は、朱色点線内★でマーキングしたFilterConvオブジェクトを調べました。

Block diagram (alpha version) of the Rx signal processing chain and the FilterConv object under consideration here.

FilterConvオブジェクト（AudioFilterConvolution_F32）

Keiths' SDRのRx信号処理連鎖では、water fall 用のFFT処理を除いて、おそらく唯一の周波数領域での信号処理ブロックです。Teensy 4.1を用いたもう一つのSDRプロジェクトであるT41-EPが周波数領域を中心とした処理系であるのとは対照的です。

経緯

ソースコードのコメント文に記された経緯を拾うと、IK8YFW局Giuseppe Callipo OMが2021/02/06にTeensy 3.6用のOpenAudio_ArduinoLibraryに実装し、W7PUA局 Bob Larkin OMが2022/03/06にTeensy 4.x用に改良した比較的新しいオブジェクトのようです。

オリジナルは、Brian Millier氏が2017/03に開発したギターアンプ用のTeensy 3.6向けのコードとのことです。Teensyがオーディオ信号処理の市場を主要なターゲットの１つにしていたことが分かります。下記に分かり易い解説が有ります。

疑問

先に調査したRX_Hilbertオブジェクト（AudioFilterFIR_F32）もFilterでした。

二重にFilterを適用しているのでしょうか？　使い分けはどうなっているのでしょうか？

GUIからの追跡

K7MDL局Mike OMが開発したタッチGUIは「Filter」タッチボタンを備えます。「Filter」タッチボタンが呼び出すコードを追跡し、フィルタ係数を切り換えている対象がRX_Hilbertオブジェクトなのか、あるいはこのFilterConvオブジェクトなのかを調査しました。結論は、このFilterConvオブジェクトでした。

Tracking results of the Filter object invoked by the "Filter" touch button in the GUI.

GUIの「Filter」ボタンに対するタッチイベントは、まずControls.cppファイルにあるFilter関数を呼んでいます。引数dirはフィルタの帯域を拡大する方向か、縮小する方向かを指示します。まだ確認できていませんが、タッチジェスチャを読み取っているものと思います。美しいGUIと合わせて、タッチイベント管理に関するMike OMの実装力には敬服するばかりです。

Filter関数は、Bandwidth2.cppファイルにあるselectBandwidth関数を呼んでいます。引数bndxは帯域設定番号であり、タッチ操作で増減させ、予め登録済みの帯域を逐次に切り換えているようです。帯域設定番号には、以下の帯域および中心周波数が割り振られています。CWモードの場合は帯域設定番号1から4の間で、SSBモードの場合は帯域設定番号5から9の間で切り換えているようです。これ以外に好みのフィルタ設定があれば、コンパイル前に増設や入れ換えも容易です。

1. 帯域：250 Hz、中心：CWピッチ（可変）
2. 帯域：500 Hz、中心：CWピッチ（可変）
3. 帯域：700 Hz、中心：CWピッチ（可変）
4. 帯域：1.0 kHz、中心：CWピッチ（可変）
5. 帯域：1.8 kHz、中心：1.850/2 kHz
6. 帯域：2.3 kHZ、中心：2.400/2 kHz
7. 帯域：2.8 kHZ、中心：2.900/2 kHz
8. 帯域：3.2 kHZ、中心：3.300/2 kHz
9. 帯域：4.0 kHZ、中心：4.100/2 kHz

selectBandwidth関数のコードの中でRX_Hilbertオブジェクトの初期化がコメントアウトされ、代わりに上記のフィルタ仕様の１つを設定した後にSDR_RA8875.inoメインファイルにあるSetFilter関数を呼んでいます。

Bandwidth2.cppファイル名にはバージョン番号と思われる「2」が付いています。おそらく当初は、RX_Hilbertオブジェクトの初期化関数によってHilbert FIRフィルタ係数を切り換え、帯域設定をしていたものと推定されます。FilterConvオブジェクトがライブラリに導入されたことによって、性能の良いFilterConvオブジェクトに乗り換えたものと思われます。その証拠に、次のようなコメント文を発見しました。

For convolutional filter method, just set a single fixed Hibert filter width. Rest is taken care of after the summer.

畳み込みフィルタ方式では、Hilbertフィルタの幅を1つ固定で設定するだけです。残りは夏以降に整理します。

RX_Hilbertオブジェクトは、最も帯域の広い4.0kHzで初期化されています。それ以外のHilbertフィルタ係数はもはや不要ということのようです。設計の一般論としては、１つの要素に２つの機能を割り付けるのは悪手とされています。二兎を追う者は一兎をも得ず、２つの機能の同時最適化が困難になることから導かれた教訓です。RX_Hilbertオブジェクトの最重要機能はイメージ抑圧のための位相制御です。位相制御と帯域制御がフィルタ設計で競合するかどうかの知見は無いのですが、少なくとも機能が分離していた方が信号処理連鎖も理解し易いと思います。

SDR_RA8875.inoメインファイルにあるSetFilter()関数は、そのままリダイレクトして、FilterConvオブジェクトのクラスを定義したAudioFilterConvolution_F32.cppライブラリファイルのinitFilterメソッドを呼んでいます。引数は下記の通りです。

• fc : filter center

• Astop : stop-band attenuation from which to calculate Kaiser-Bessel window shape factor beta（90dBに固定）

• type : filter type（2=BANDPASSに固定）

• dfc : filter band width

initFilterメソッドは、calc_FIR_coeffsメソッドとimpulseメソッドを呼んでいます。calc_FIR_coeffsメソッドは、フィルタ切換の都度、フィルタ係数をその場（on the fly）で再計算します。これにより、CWピッチをユーザの好みでコンパイル後も変更可能になります。腕に覚えがあれば、GUIによる帯域シフト機能や帯域幅調整機能を実装することも可能でしょう。

impulseメソッドは、FFTによりフィルタ係数を周波数領域のマスクに変換します。時間領域の畳み込み演算は周波数領域のマスク乗算になります。

信号処理連鎖の更新処理

update()関数は4段階の状態遷移マシンになっています。つまり、eDMA割込みによるデータブロック処理を4回実行して一周します。これは、割込み処理の4回につき1回だけ512（=4*128）データ長のFFT処理を実行するためです。各状態で以下の処理を行っています。

1. 状態=0：畳み込み計算と逆FFT、データ出力（1回目）
   ・周波数領域のデータとフィルタマスクとの畳み込み演算を実行。
   　（データ長は1024*2）
   ・畳み込み演算結果を逆FFTにより時間領域データに変換。
   　（データ長は1024*2）
   ・前半（0～511）のオーバラップ長512のデータを前回の結果と加算。
   　後半（1024～1535）のオーバラップ長512のデータを次回の準備として保存。
   ・最初の128（0～127）個の時間領域データのブロックを次の信号処理オブジェクトへ送信。
   ・メモリブロックの解放。
2. 状態=1：データ出力（2回目）
   ・2つ目の128（128～255）個の時間領域データのブロックを次の信号処理オブジェクトへ送信。
   ・メモリブロックの解放。
3. 状態=2：データ出力（3回目）、FFTバッファの準備
   ・3つ目の128（256～383）個の時間領域データのブロックを次の信号処理オブジェクトへ送信。
   ・メモリブロックの解放。
   ・FT用バッファの後半（1024～2047）をゼロで埋める。
4. 状態=3：データ出力（4回目）、FFT
   ・4つ目の128（256～383）個の時間領域データのブロックを次の信号処理オブジェクトへ送信。
   ・メモリブロックの解放。
   ・FFT用バッファの前半（0～1023）にeDMA割込み4回分の512（=4*128）個の時間領域データを格納。
   ・時間領域データをFFTにより周波数領域データに変換。

FFTは512（=4*128）データ長の周期信号を仮定します。実際の信号はそのような周期信号ではないため、オーバラップ処理がFFTのエイリアスを緩和しているようです。

カイザー窓 FIR フィルタのシミュレーション

カイザー窓 FIR フィルタ

Kaiser-Bessel windowを計算するための上記引数Astop（stop-band attenuation from which to calculate Kaiser-Bessel window shape factor beta）の存在から明らかなように、FilterConvオブジェクトが採用しているフィルタは「カイザー窓 FIR 帯域通過フィルタ」です。この方法の利点はFIRタップ数も同時に推定してくれる点ですが、Keiths' SDRではFIRタップ数は513に固定されています。FIRタップ数固定により、スロープの傾斜は帯域設定に係わらず同じになります。

Design method for Kaiser-Bessel window FIR bandpass filter on the fly.

Wheatley, M. (2011): CuteSDR Technical Manual. www.metronix.com, pages 118 - 120, FIR with Kaiser-Bessel Window.

initFilterメソッドのカイザー窓 FIR フィルタ係数初期化コードは、上記CuteSDRの文献の計算式を参照して、WMXZ（ハンドル名？）OMとDD4WH局OMが2016年にTeensy 3.6に実装したコードをオリジナルとして、W7PUA局 Bob Larkin OMがOpenAudio_ArduinoLibraryに移植し、K7MDL局Mike OMがKeiths' SDR用に実装したコードのようです。

まず、dfc（filter bandwidth）を遮断周波数にしたLPFを設計してから、fc（filter center）をシフト周波数にしてBPFに変換しています。これによって、LPFの遮断周波数がBPFの通過帯域に変換されるようです。

既に2011年に、Technical Manualが付いたSDRが存在していたんですね。今やSDRのTechnicalな側面（理論と工学の狭間）は一般常識化したのか、CuteSDRのようなTechnical Manualが付いたSDRは見かけません。WDSPライブラリにはProgrammer's Guideが付いていますが、APIの説明が主であり、Technicalな側面はほとんど分かりません。T41-EPのAmazon解説本は貴重な存在ですが、信号処理のTechnicalな側面は他の文献を参照する方式です。一般的には、ソースコードのコメント文が唯一貴重なTechnicalリソースです。

その場フィルタ設計のシミュレーション方法

AudioFilterConvolution_F32::calc_FIR_coeffsメソッドのC++コードをPython（Jupyter Notebook）に移植しました。カイザー窓 FIR フィルタの設計に特殊な関数は必要なく、mathライブラリの三角関数と平方根関数で計算できます。ただし、0次修正ベッセル関数の計算に収束演算が必要です。Pythonには変数の型宣言がないため、int変換関数やfloat変換関数を適宜適用する必要がありました。

各フィルタの帯域中心fc（filter center）と帯域幅dfc（filter band width）を与えてBPFを設計し、513タップFIRフィルタ係数を出力しました。FIRフィルタ係数をFFTによりフーリエ変換して周波数特性を求めました。下図の左からインパルス応答、周波数特性（横軸固定）、周波数特性（横軸拡大）を示します。帯域の広いフィルタから順にシミュレーション結果を示します。

その場フィルタ設計のシミュレーション結果

9.  帯域：4.0 kHZ、中心：4.100/2 kHz

8. 帯域：3.2 kHZ、中心：3.300/2 kHz

7. 帯域：2.8 kHZ、中心：2.900/2 kHz

6. 帯域：2.3 kHZ、中心：2.400/2 kHz

5. 帯域：1.8 kHz、中心：1.850/2 kHz

4. 帯域：1.0 kHz、中心：CWピッチ（750Hz）

3. 帯域：700 Hz、中心：CWピッチ（750Hz）

2. 帯域：500 Hz、中心：CWピッチ（750Hz）

1. 帯域：250 Hz、中心：CWピッチ（750Hz）

気付きを以下にまとめます。

1. Hilbertフィルタとは異なり、インパルス応答（フィルタ係数）は左右対称。
2. 周波数帯域を狭くするほど、インパルス応答（フィルタ係数）の幅は広くなり、振幅は小さくなる。
3. dfc=250HzのCW用最狭帯域フィルタを除いて、-6dBcの帯域幅は目分量で設計仕様を満たしている。
4. FIRタップ数固定により、スロープの傾斜は帯域設定に係わらず同じ。これにより、CW用狭帯域フィルタの通過域の平坦度が影響を受けているかもしれない。
   ※ゲインが減少した分は、別途補正されている。

悩み再び！？（2022/08/12）

AudioFilterConvolution_F32::calc_FIR_coeffsメソッドによるカイザー窓 FIR フィルタのその場設計では、サンプリング周波数がパラメータとして必要です。信号連鎖の全てのオブジェクトのコンストラクタは、サンプリング周波数=96kHzとデータブロック長=128の２つを引数として取り、これらをオブジェクトのprivate変数として保持し、必要に応じて計算に使用しています。つまり、カイザー窓 FIR フィルタはサンプリング周波数=96kHzで設計されているのです。このフィルタが設計仕様通りに動作しているなら、Keiths' SDR（K7MDL局Mike OM版）のサンプリング周波数は96kHzで正しいということになります。

逆に、Hilbert FIR FilterはPJRC社ディフォルトのオーディオ周波数44.1kHzで設計され、そのまま放置されているのではないか？との疑念が浮上しました。44.1kHzで設計したHilbert FIR Filterを96kHzで使用すると通過帯域が約2倍に拡大します。しかし、後段のFilterConvオブジェクトで帯域を絞るため、実害はノイズ低減機能に影響する可能性に止まり、特にノイズが大きくない環境では設計のミスマッチに気付かない可能性が高いのです。

Keiths' SDRのコミュニティでコーディックのサンプリング周波数をオーバークロックさせたりの検討をしていた模様ですが、その都度Hilbert FIR Filterを再設計していたかどうかは定かではありません。Hilbert FIR Filterのフィルタ係数はコンパイル時にコードに埋め込んでおり、その場設計ではないため、変更していない可能性が高いのではないかと推定しています。

最終的にはハードウェアを組んで、信号処理連鎖の起動周期を測定しないと白黒が付かない状況になりました。

Teensy（12）Keiths' SDRのLMS_Notchオブジェクト

Keiths' SDRのRx信号処理連鎖ブロック図

Keiths' SDR（K7MDL局Mike OM版）のRx信号処理連鎖をコードから読み解いて抜き出したもの（α版）を下記に再録します。

今回は、朱色点線内★でマーキングしたLMS_Notchオブジェクトを調べました。

Block diagram (alpha version) of the Rx signal processing chain and the LMS_Notch object under consideration here.

LMS_Notchオブジェクト（AudioLMSDenoiseNotch_F32）

LMS_Notchオブジェクトは、学習アルゴリズムを作り込むことができるSDRの優位点が出る機能の１つと思います。都市ノイズが増えた今日、自環境に合わせて最適化すれば有用なオブジェクトの１つになるかもしれません。

役割

自動ノッチ処理を行う信号処理オブジェクトです。前回のNoiseBlankerオブジェクトが間欠的なパルスノイズの除去を対象にしていたのに対して、LMS_Notchオブジェクトは定常的なノイズの除去を対象にしています。頭を切り換えると、同じ構成で白色ノイズ低減にも使えます。

トーンノイズのように、音声信号よりも時間的コヒーレンス（自己相関）が強いノイズの場合は、自動ノッチ機能として働きます。音声信号よりもコヒーレンスの強いトーンノイズに学習適応したFIRフィルタによってノイズを選択出力し、入力信号から差し引いてノイズを除去します。

反対に、白色ノイズのように、音声信号よりも時間的コヒーレンスが弱いノイズの場合は、ノイズよりもコヒーレンスが相対的に強い音声信号を通過させるようにFIRフィルタを学習適応させます。これにより、音声信号を選択出力してノイズを低減します。

Two modes of operation for the LMS_Notch object.

前者がトーンノイズ自動ノッチ、後者が白色ノイズ低減の動作になります。一見、異なる機能のようですが、コヒーレンス特性が強い信号にFIRフィルタが同調するというアルゴリズムの働きは同じです。ただし、FIRフィルタの出力が180度変わることに頭の切り換えが必要です。オブジェクトを2重直列にして、自動ノッチの後にノイズ低減を作用させたらどうなるのか、信号遅延に見合う効果がえられるのか、ノイズも信号も全て消えてしまうのか、興味のあるところです。

信号処理連鎖の更新処理

update()関数では以下の信号処理を行っています。

1. データブロックの128サンプルに以下の処理を逐次適用。
   ・フィルタ係数の1つにLeakage（忘却係数を乗じる）処理を逐次適用。
   ・遅延循環バッファの最後尾に、入力データブロックから新規データを挿入。
   ・新規データを反映した正規化パワーを計算。
   ・FIRフィルタ内部循環バッファの最後尾に、遅延循環バッファの最前部のデータを挿入。
   ・FIRフィルタ内部循環バッファのデータにFIRフィルタ係数をタップ毎に乗じて総和を出力。
   ・入力データとFIRフィルタ出力の差分エラーを計算。
   ・差分エラーから正規化学習係数を計算。
   ・FIRフィルタ内部循環バッファのデータに学習係数を乗じてFIRフィルタの係数に加算し、FIRフィルタ係数を適応学習。
   ・ノイズ低減モードの場合はFIRフィルタの出力をLMS_Notchオブジェクトの出力に設定し、自動ノッチモードの場合は差分エラーをLMS_Notchオブジェクトの出力に設定。
2. 次の信号処理オブジェクトへのデータ送信。
3. メモリ・ブロックの解放。
   注意：FIRフィルタの係数と内部循環バッファのデータは、次回の信号処理連鎖起動に向けて保持される。

信号処理連鎖のデータブロック単位の再考

IQ信号データは、128サンプルのブロック単位で信号処理連鎖の各オブジェクトに送信されて処理されます。遅ればせながら、このLMS_Notchオブジェクトの調査段階で認識しました。

遡ると、コーディックからeDMAで転送されるIQ信号データも128サンプルのブロック単位になります。よって、信号処理連鎖が起動される周波数は、98kHz / 2 / 128 = 375Hzとなり、周期は2.7msecになります。データブロックが1サンプルだけの極端な場合は周期が21μsecになることから、CPUの割込み負荷には大きな違いが生じます。データブロック単位で信号処理を行うのは、CPUの負荷の変動があっても、オーディオ処理のストリーミングが途切れないようにする処置と思います。

ブロック長が128サンプルであることは、16bit整数のTeensy Audio Libraryのディフォルトの仕様を踏襲しています。元祖Teensy Audio Libraryは既に任意のブロック長に対応しています。しかし、32bit浮動小数点のOpenAudio_ArduinoLibraryを使用するKeiths' SDRでは、ブロック長を128サンプルから変えないように注意コメントが付いています。

const int   audio_block_samples = 128;          // do not change this!

どこかのオブジェクトで「128」がハードコーディングされている（コンパイル後に変更できない数値としてコードに埋め込まれている）ことが理由と思われます。Teensy 4.1の計算能力は向上しているため、ブロック長を短くしても良いと思いますが、ハードコーディングでは致し方なしです。

LMS_Notchオブジェクトはハードコーディング・オブジェクトの１つでした。128を埋め込んだfor文でFIRフィルタ計算を行っています。LMS_Notchオブジェクトでは、CMSIS DSP LibraryのFIRフィルタ関数ではなく、オリジナル・コードで計算しています。適応学習機能を組み込むためと思われます。

LMSアルゴリズム

歴史

LMS_Notchオブジェクトは、LMS（Least mean squares：最小平均二乗）学習アルゴリズムに基づく適応フィルタが名称の由来になっています。ソースコードのコメントを参照すると、今年に入ってから（January 2022）、W7PUA局 Bob Larkin OMによって実装された新しいオブジェクトのようです。1999年に発表されたDSP-10（All-Mode 2m Transceiver）からの移植ではないかと思われます。

歴史を紐解くと、ベースとなっているLMSアルゴリズムは Bernard Widrow博士らによって1960年に提案されたものとのことです。博士は、IEEE Centennial Medal (1984)、IEEE Alexander Graham Bell Medal (1986)、IEEE Neural Networks Pioneer Medal (1991)、IEEE Millennium Medal (2000)、Benjamin Franklin Medal (2001)を授与されているレジェンドです。LMSアルゴリズム発見の経緯については、次の文献にて博士自ら回顧されています。

B. Widrow "Thinking about thinking: the discovery of the LMS algorithm," in IEEE Signal Processing Magazine, vol. 22, no. 1, pp. 100-106, Jan. 2005, doi: 10.1109/MSP.2005.1407720.

Widrow博士は学位論文を書き上げた後のMITの教壇に立つまでの余暇に、彼の有名なダートマス会議（1956/7-8）に参加されていたとのこと。彼の地で人工知能を夢見ていた人達の裾野の広さには驚きます。ダートマス会議にはポスドクの聴講者として参加されていたようですが、後にIEEE Neural Networks Pioneer Medalを授与されています。

Stanford大に移ってから書き上げた1960年の原著論文ではLMSアルゴリズムを「Neuron」という言葉で論じており（下記Fig. 9参照）、ダートマス会議で触発された可能性が推測されます。もっとも、上記回顧録では、2値化Neuron導入の契機はMIT時代の修士院生の提案ということになっています。最急降下法アルゴリズムが完成したのは、Stanford大に移ってからの博士院生の指導の過程だったようです。原著論文はその博士院生との共著になっています。

The LMS algorithm was discussed under the term "Neuron", for eample, in Figure 9 of the original paper: WIDROW B. "Adaptive Switching Circuits," 1960 IRE WESCON Convention Record, pp. 96-104  (1960)

1959年の週末にパーツショップで買い集めた部品で製作したアナログコンピュータ（ADALINE）でアイディアを検証したというエピソードにも驚きます。 Intel 4004がリリースされたのは1971年です。1959年当時は、黒板とアナログコンピュータが研究ツールだったんですね。

ADALINE（ADAptive LInear NEuron）

ノイズ除去用FIRフィルタのLMSアルゴリズム

ノイズ除去用FIRフィルタのLMSアルゴリズムについては、下記QEX誌に詳しい解説があります。

Steven E. Reyer, WA9VNJ, and David L. Hershberger, W9GR: “Using The LMS Algorithm For QRM and QRN Reduction,” QEX, September 1992, pp 3-8.

ノイズ除去用FIRフィルタのLMSアルゴリズムは、事前の設定が必要なパラメータを4つ持ちます。時間的コヒーレンス特性を持つ疑似参照信号を作るためのdelay（遅延バッファ長）、FIRフィルタのタップ数、フィルタ係数を更新する際の学習速度を定めるbeta（論文ではμ）、信号が無いときにフィルタ効果を漸近的にOFFにする忘却？係数decayです。LMS_Notchオブジェクトでは、これら4つのパラメータの初期値は以下のように設定されています。

• delay バッファ長さ = 4（Max. 16）
• FIRフィルタのタップ数 = 64（Max. 128）
• beta = 0.03
• decay = 0.995　（1-decay = 0.005）

速くノイズをキャンセルしたいがためにbetaを大きくしたいところですが、発散を防ぐためには匍匐前進が必要になります。betaの最適値は自局のノイズ環境に依存するため、アルゴリズムが自動化さても最適な学習速度を見出すという運用スキルの余地は残るようです。

自動ノッチの対象になるトーン・ノイズに対して、音声信号は相対的にコーヒレンス（相関）が弱いと見なします。逆に、ノイズ低減の対象になる白色ノイズに対しては、音声信号は相対的にコーヒレンスが強いと見なします。具体的にアルゴリズムの動作を切り換えるためには出力を取る位置を変え、時間的コヒーレンスを持つ疑似参照信号を作るためのdelayや学習速度を定めるbetaを調整する必要があります。

自動ノッチのシミュレーション

音声データ

PC上で、ノイズ除去用FIRフィルタのLMSアルゴリズムのシミュレーションを試みました。bit数等をTeensy 4.1に合わせ込んだ厳密なシミュレーションではありません。

試験用の音声データを下記に示します。振幅は16bit整数最大値（32,767）で正規化してあります。サンプリング周波数は、48kHzではなく、CDクオリティの44.1kHzです。約90万点（21秒）の音声データから、信号が連続する期間を抜き出して利用しました。

Voice data for testing.

シミュレーションに用いた音声データ

上記の音声データの100,000点目を始点として1,024点（23msec）を切り出し、自動ノッチのシミュレーション用データとして利用しました。切り出した音声データの時系列波形、周波数スペクトル、自己相関関数を下記に示します。

Voice data for simulation of automatic notch.

Spectrum of voice data for simulation of automatic notch.

Autocorrelation function of voice data for simulation of automatic notch.

参照信号を作る遅延バッファのサイズはディフォルトで4 点（91μsec）です。自己相関関数から、僅か4点の遅延では音声データの相関がまだまだ高いことが分かります。トーンノイズとコーヒレンスの強さにおいて、FIRフィルタの適応を競うことになるのではないかと心配されます。

シミュレーションに用いたトーンノイズ

音声信号の最大振幅に対して3倍の振幅を持つ1kHzのトーンノイズを加えて、シミュレーション用の試験データとしました。時系列波形データおよび周波数スペクトルを下記に示します。

Time-series waveform of the test data with tone noise added.

音声信号の時系列波形データはトーンノイズに埋もれてしまいました。

Frequency spectrum of the test data with tone noise added.

周波数スペクトルでは、1kHzの周波数のみ強度が伸長すると思っていましたが、1kHzの周囲のサポート部も強度が大きくなっています。また、5k～8kHzの高周波領域の強度も大きくなっています。

1kHzのトーンノイズの周期が1024点と一致していないため、境界で不連続点が発生しているのかもしれません。不連続点を解消する窓関数がディフォルトで適用されていると思ったのですが、今回利用したnumpyのfft関数のマニュアルには言及が無いようでした。

自動ノッチのシミュレーション結果

参照信号を作るdelayバッファおよびFIRフィルタの内部バッファの初期値はゼロとして、1024点のトーンノイズを重畳した音声信号データに対して自動ノッチのシミュレーションを行いました。アルゴリズムのパラメータは、LMS_Notchオブジェクトの初期値に合わせて下記としました。

• delay バッファ長さ = 4
• FIRフィルタのタップ数 = 64
• beta = 0.03
• decay = 0.995　（1-decay = 0.005）

FIRフィルタ係数の学習過程を下記に示します。左（１）は全1024ステップ、右（２）は最初の100ステップの拡大です。

Variation of FIR filter coefficients by adaptive learning.

1024ステップ目でも厳密には収束していませんが、約200ステップで大きな変動は無くなっています。途中で脈動が見られるのは、decayパラメータによる忘却作用と思われます。

右（２）の最初の100ステップの拡大を見ると、バッファ初期値ゼロによる遅れ時間の後、一斉にプラスの係数が増大していきます。20ステップ後半でマイナスの係数が出現し、その中から50ステップ過ぎにプラスに転じる係数も出てくるようです。

最後の1024ステップ目のFIRフィルタ係数とその周波数特性を下記に示します。

Characteristics of the FIR filter coefficients for the 1024th step.

右（２）の周波数応答から、1kHzのトーンノイズを通過させるようにFIRフィルタが適応していることが分かります。ただし、若干1kHzより高めにシフトしています。また、シャープなノッチ特性ではなく、期待していたよりブロードな特性を示し、多少の音声データも道ずれに抑制してしまいそうな特性となりました。

左（１）のFIRフィルタ係数（インパルス応答）を見ると、自己相関関数の様相を示していることに気付きます。1kHzのトーンノイズの半周期は22 samplesとなり、FIRフィルタ係数の極小値と近い位置になります。1周期は44 samplesになり、極大値と近い位置になります。

下記にトーンノイズの自己相関関数を示します。参考のため、音声信号の自己相関関数も横軸をFIRフィルタのタップ数に合わせて再掲します。

Comparison of the autocorrelation function of the tone noise with the autocorrelation function of the voice signal.

左（１）のトーンノイズの自己相関関数の極小値と極大値の位置が、FIRフィルタ係数の極小値と極大値の位置と近いことが分かります。右（２）の音声信号はその限りではありません。LMSアルゴリズムの数理最適解であるWiener解は自己相関関数と相互相関関数から計算されるとのことなので、詳細は不肖ながら何らかの関係があるのかもしれません。

1kHzトーンノイズを加えた音声信号に自動ノッチを適用した後の時系列波形データを下記に示します。FIRフィルタ係数の学習過程で見たように、約200ステップから後はトーンノイズが抑圧されています。ただし、音声信号の振幅も小さくなっているようです。

Time series waveform data after auto-notch is adapted to the voice signal with 1 kHz tone noise added.

トーンノイズが十分に抑圧された400 sampleより後ろのデータを切り出し、周波数スペクトルを調べた結果を下記に示します。

Frequency spectrum after auto-notch is adapted to the voice signal with 1 kHz tone noise added.

上記の自動ノッチ適用前の周波数スペクトルと比較すると、1kHzトーンノイズの強度は、約17dBから8dBまで11dBc低減しています。5k～8kHzの高周波領域への影響も消滅しています。ただし、1kHz以下の音声スペクトルの強度も小さくなっています。また、2kHz強にもゲインの落ち込みがあります。これらが、時系列波形データの振幅の低下に呼応しているものと思われます。

予めトーンノイズの周波数が分かっている場合は、参照信号の周波数を固定して、位相と振幅のみ適応させるように構成した方が良いかもしれません。SDRをソフトウェアも含めて自家で構築すれば、ハンダ鏝をキーボードに代えて、そのような実験試行も自由に検討できるようになります。

ノイズ低減のシミュレーション

シミュレーションに用いた音声データ

事前の試行錯誤により、トーンノイズに比べてコーヒレンスの小さい音声データにFIRフィルタを適応させるためには学習ステップが長くなると推定されたため、トーンノイズのシミュレーションよりも長い区間の音声データを切り出しました。具体的には、170,000点目を始点として4,096点（93msec）を切り出し、ノイズ低減のシミュレーション用データとして利用しました。切り出した音声データの時系列波形、周波数スペクトル、自己相関関数を下記に示します。

Voice data ffor simulation of automatic noise reduction.

Spectrum of voice data for simulation of automatic noise reduction.

Autocorrelation function of voice data for simulation of automatic noise reduction.

切り出し区間を長くした結果、音声信号には２つの音素が含まれているように見えます。周波数スペクトルは分解能が上がりましたが、低域の挙動は自動ノッチの試験データに似ていると思います。高域は3kHzと6.5kHz付近の信号強度が高くなっています。自己相関関数は音素の違いから周期が異なりますが、大勢は自動ノッチの試験データに似ています。

シミュレーションに用いた白色ガウスノイズ

音声信号の最大振幅に対して1/10の振幅を持つ白色ガウスノイズを加えて、シミュレーション用の試験データとしました。時系列波形データおよび周波数スペクトルを下記に示します。

Time-series waveform of the test data with white Gaussian noise added.

音声信号の時系列波形データに高周波信号が付加されたイメージです。

Frequency spectrum of the test data with white Gaussian noise added.

周波数スペクトルでは、4kHz以上の高域の強度がノイズで持ち上がっています。1.2kHz以下の低周波帯域ではノイズの影響を受けていません。

ノイズ低減のシミュレーション結果

参照信号を作るdelayバッファおよびFIRフィルタの内部バッファの初期値はゼロとして、4096点の白色ガウスノイズを重畳した音声信号データに対してノイズ低減のシミュレーションを行いました。アルゴリズムのパラメータは、試行錯誤で調整し下記としました。

• delay バッファ長さ = 4
• FIRフィルタのタップ数 = 128
• beta = 0.015
• decay = 0.995　（1-decay = 0.005）

delayバッファ長さとdecay忘却係数は自動ノッチと同じですが、FIRフィルタのタップ数を2倍の128（最大）とし、学習係数betaを1/2の0.015に設定しました。

FIRフィルタ係数の学習過程を下記に示します。左（１）は全4096ステップ、右（２）は最初の100ステップの拡大です。

Variation of FIR filter coefficients by adaptive learning.

4096ステップ目でも収束していませんが、最初の100ステップだけ拡大して見ると、100ステップで収束しているように見えます。

最後の4096ステップ目のFIRフィルタ係数とその周波数特性を下記に示します。

Characteristics of the FIR filter coefficients for the 4096th step.

右（２）の周波数応答から、1.2kHz以上の高域を徐々に絞るようにFIRフィルタが適応していることが分かります。ただし、1.2kHz以下の低域も平坦ではありません。音素の数が高々2個の音声データに適応した結果と思います。さらに長時間学習することで平坦な特性に近づく可能性はあります。

左（１）のFIRフィルタ係数（インパルス応答）を見ると、デルタ関数の様相を示していることに気付きます。白色ノイズの自己相関関数は、理論通りならデルタ関数になると思います。

下記に白色ガウスノイズの自己相関関数を示します。参考のため、音声信号の自己相関関数も再掲します。

Comparison of the autocorrelation function of the white Gaussian noise with the autocorrelation function of the voice signal.

左（１）の白色ガウスノイズの自己相関関数がデルタ関数になっていることが分かります。ただし、FIRフィルタ係数（インパルス応答）がデルタ関数になってしまうと、入力データがそのまま出力されてしまうため、フィルタ効果はなくなります。FIRフィルタ係数のスロープ部分の畳み込み演算で移動平均処理を行い、白色ガウスノイズを打ち消して音声信号のみ出力するように適応学習が進んでいると思います。

白色ガウスノイズを加えた音声信号に自動ノイズ低減を適用した後の時系列波形データを下記に示します。目視による定性的な印象になりますが、大部分のノイズは削減されているようです。ただし、音声信号の振幅も75%程度に低減されてしまったようです。

Time series waveform data after automatic noise reduction is adapted to the voice signal with white Gaussian noise added.

周波数スペクトルを調べた結果を下記に示します。

Frequency spectrum after automatic noise reduction is adapted to the voice signal with white Gaussian noise added.

上記の自動ノイズ低減適用前の周波数スペクトルと比較すると、4kHz以上の高域のノイズが低減されていることが分かります。1.2kHz以下の音声信号の強度は変化が無いようです。ただし、3kHzと6.5kHz付近の音声信号の強度が低減されています。これが、音声信号の時系列波形の振幅が75%程度に低減されてしまった理由かもしれません。今回、再生音には変換していませんが、音声の伸びやメリハリが低減されてしまっているかもしれません。

ノイズ低減の効果は期待できますが、完璧ではなく、工夫の余地が残されている・・・という感じでしょうか。

ノイズ低減機能の展開

NR0V局Warren C. Pratt OMが、WDSPと称するSDRのためのopen-source Digital Signal Processing libraryをgithubに公開されています。経緯を記した情報が少なく、不確かなのですが、OpenHPSDR（High Performance SDR）のために開発され、mcHF QRP transceiverのfirmwareであるUHSDR（UniversalHamSDR）等にも再利用されているようです。

WDSPにはLMS Noise Reductionの他に、Spectral Noise Reductionが用意されています。未読ですが、少し昔の1984年の下記参考文献が引用されています。半導体の進歩によって、ようやく手の届く所に来たアルゴリズムという感じでしょうか・・・。

Yariv Ephraim and David Malah, "Speech Enhancement Using a Minimum Mean‐Square Error Short‐Time Spectral Amplitude Estimator," 1984.

The WDSP Guide Using WDSP ‐ for Software Developersには、信号をFFTで周波数領域に変換した後、周波数bin毎に音声信号とノイズの強度を推定してgainを調整すると記されています。音声信号をガウス分布もしくはガンマ分布で推定するとのことですが、そのような単峰対称分布の仮定で良いのか少し懐疑的に思います。音声信号とノイズの強度の推定が課題ですが、それが上手く行くなら、今回のシミュレーションで課題になった3kHzと6.5kHz付近の音声信号の強度低減を解消できるかもしれません。

Keiths' SDRの今まで見てきた信号処理オブジェクトは、全て時間領域で処理されています。（実は最後のFilterConvだけは周波数領域の処理オブジェクトです。）一方、同じTeensy 4.1を用いたもう一つのSDRプロジェクトであるT41-EPは、対照的に中核の信号処理を周波数領域で実行しています。コードはまだ見ていませんが、下記解説本には「UHSDRのSpectral Noise Reductionを実装した」旨の記述があります。転じて、Keiths' SDRにも実装できそうですね。

Teensy（11）Keiths' SDRのHilbertオブジェクト

修正（2022/07/28）
信号処理連鎖を流れるデータの単位は128サンプルのデータブロックであることを認識し、修正しました。

修正（2022/07/23）
フィルタ特性をLPFとして見ていましたが、BPFとして見るのが適切でした。特にCW用の狭帯域フィルタはBPFとして見る必要があります。通過帯域幅以外のフィルタ設計諸元は不明ですが、周波数特性図に帯域中央周波数を推定して記入し、コメントも修正しました。

Keiths' SDRのRx信号処理連鎖ブロック図

Keiths' SDR（K7MDL局Mike OM版）のRx信号処理連鎖をコードから読み解いて抜き出したもの（α版）を下記に再録します。

今回は、朱色点線内★でマーキングした３つのオブジェクトを調べました。ノイズブランカ、ヒルベルトフィルタ、加算器です。

Block diagram (alpha version) of the Rx signal processing chain and the NoiseBlanker, Hilbert, and Summer objects under consideration here.

NoiseBlankerオブジェクト（radioNoiseBlanker_F32クラス）

役割

衝撃性パルス雑音を空白（blank）に置換して抑制します。

空白置換は急峻な信号変化となるため、周波数帯域が広がってしまいます。そのため、Noise Blankerは高域をカットするHilbert FIRフィルタの前に挿入されています。

信号処理連鎖の更新処理

update()関数では以下の信号処理を行っています。

1. データブロックの128サンプルに以下の処理を逐次適用
   ・長さ256（5.3msec）のリング・バッファの1stデータに新規データの絶対値を積み、その総和を算出
   ・総和にしきい値（1.0E3f）を乗じた値よりデータの絶対値が大きければ、衝撃性ノイズと見做し、データ値をゼロに置換（blanked out）
2. 次の信号処理オブジェクトにデータを送信
3. データブロックの解放

ディフォルトでは、平均信号強度の約4倍のパルス入力をノイズと見做す設定のようです。しきい値を変更する関数も準備されています。

RX_Hilbert_...オブジェクト（AudioFilterFIR_F32クラス）

役割

±45度のHilbert FIRフィルタのオブジェクトです。信号の低域通過によるエイリアス等の除去と、位相シフトによるイメージ除去の準備を行います。（後段の加算器Summerにてイメージ除去が行われます。）

PJRC社のオーディオ・ボードはアナログ・アンチエイリアス・フィルタを備えていません。RFフロントエンドにRS-HFIQ（HobbyPCB）を使用した場合は、RS-HFIQ側でQSDの出力をオペアンプのLPFに通してから出力しています。カットオフ周波数はfc=145kHz（C=220pF、R=4.99kΩ）です。販売元のHobbyPCBが推奨するUSBサウンドカード（StarTech - ICUSBAUDIO2D）は96KHz、24-bitの仕様です。これに対して、PJRC社のオーディオ・ボードは48KHz、16-bitの仕様です。サンプリング周波数fs=48kHz、ナイキスト周波数fn=24kHzに対してfc=145kHzのLPFは、有効なアンチエイリアス・フィルタにはならないようです。将来、エイリアス雑音が課題になる可能性があります。

Hilbert FIRフィルタの係数はオブジェクト宣言時ではなく、後にbegin()関数を呼んで設定しています。予め複数の遮断周波数（500Hz、700Hz、1.0kHz、1.8kHz、2.3kHz、2.8kHz、3.2kHz、4.0kHz）に対応したフィルタ係数が、ヘッダ・ファイルHilbert.hの中に定義されています。タップ数は151タップで共通です。GUI操作を通して適切なバンド幅のフィルタ係数に切り換えられるようになっているようです。SDRの利点として、好みのフィルタを追加することも容易と思います。

立ち上げ時には、最も帯域の広い4.0kHzのフィルタを設定しています。

    RX_Hilbert_Plus_45.begin(Hilbert_Plus45_40K,151);   // Left channel Rx
    RX_Hilbert_Minus_45.begin(Hilbert_Minus45_40K,151); // Right channel Rx

信号処理連鎖の更新処理

update()関数では以下の信号処理を行っています。

1. 128サンプルのデータブロックに以下の処理を一括適用
   ・入力データブロック（のポインタ）の取得
   ・FIRフィルタのタップ数に合わせてデータブロック長を調整
   ・arm_fir_f32()関数のコール
2. 次の信号処理オブジェクトへのデータ送信
3. メモリ・ブロックの解放

単純な積和演算ですが、メモリ・ブロックの生成、解放に注意を払った実装がされています。余計なデータの複写を避けて、可能な限り演算時間を最小にする工夫と思われます。

CMSIS DSP Library: arm_fir_f32()

RX_Hilbert_...オブジェクトのupdate()関数では大量の積和演算が発生するため、CMSIS（Cortex Microcontroller Software Interface Standard）DSP Libraryのarm_fir_f32()関数をコールしています。このライブラリ関数はCortex-MプロセッサのDSP演算に最適化されてコンパイル済みです。

入力データブロック長128サンプルに対して、フィルタのタップ数は151で一致しません。不一致の場合は、内部のデータブロック長を調整するCMSIS DSP Libraryのarm_fir_init_f32()関数を呼んで再初期化をしています。フィルタを切り換えて再びタップ数が変更されたら、arm_fir_init_f32()関数を呼んで再々初期化を行います。このように、タップ数の異なるフィルタ切換を任意に適用することが保証されています。

CMSIS DSP Libraryのarm_fir_f32()関数は、フィルタのタップ数に合わせた内部データバッファの状態をarm_fir_instance_f32インスタンスとして、関数から戻った後も保持します。次回の信号処理連鎖が起動した時に、次回のデータブロックに含まれる128サンプルの頭からフィルタ計算を継続して行うために、現在の内部データバッファの状態を引き継ぐためです。

CMSIS DSP Libraryの関数はコンパイル済みのため、中身を調査することはできないと思っていましたが、githubにコードが公開されていました。

https://github.com/ARM-software/CMSIS/blob/master/CMSIS/DSP_Lib/Source/FilteringFunctions/arm_fir_f32.c

フィルタ特性

Hilbert.hのフィルタ係数をPython（Jupyter Notebook）にExportし、各フィルタのカットオフ周波数やスロープ特性等を確認するために可視化しました。上段に位相シフト+45度のフィルタの特性を、下段に位相シフト-45度のフィルタの特性を示します。左からインパルス応答、周波数特性（ゲイン）、周波数特性（位相）を示します。

カットオフ周波数狙い値4.0kHz

Hilbert FIR filter with a cutoff frequency of 4.0 kHz.

カットオフ周波数狙い値4.0kHzにおいて、-3dBの減衰が達成されているように見えます。-3dBの通過帯域幅は設計仕様通りの4.0kHzです。帯域中央周波数は2.0kHzより若干大きい値です。スロープは急峻です。

カットオフ周波数狙い値3.2kHz

Hilbert FIR filter with a cutoff frequency of 3.2 kHz.

カットオフ周波数狙い値3.2kHzにおいて、-10dB以上の減衰過達に見えます。-3dBの通過帯域幅は3.2kHzより狭いようです。

カットオフ周波数狙い値2.8kHz

Hilbert FIR filter with a cutoff frequency of 2.8 kHz.

カットオフ周波数狙い値2.8kHzにおいて、-10dB減衰が達成されているように見えます。-3dBの通過帯域幅は2.8kHzより狭いようです。

カットオフ周波数狙い値2.3kHz

Hilbert FIR filter with a cutoff frequency of 2.3 kHz.

カットオフ周波数狙い値2.3kHzにおいて、-10dB弱の減衰が達成されているように見えます。-3dBの通過帯域幅は2.3kHzより狭いようです。

カットオフ周波数狙い値1.8kHz

Hilbert FIR filter with a cutoff frequency of 1.8 kHz.

カットオフ周波数狙い値1.8kHzにおいて、-3dBの減衰が達成されているように見えます。-3dBの通過帯域幅は設計仕様通りの1.8kHz。帯域中央周波数は約800Hz。スロープは緩慢です。

カットオフ周波数狙い値1.0kHz

Hilbert FIR filter with a cutoff frequency of 1 kHz.

カットオフ周波数狙い値1kHzにおいて減衰はなく、通過域の中央のように見えます。-3dBの通過帯域幅は設計仕様通りの1kHz。帯域中央周波数も約1kHz。スロープは緩慢です。

カットオフ周波数狙い値700Hz

Hilbert FIR filter with a cutoff frequency of 700 Hz.

カットオフ周波数狙い値700Hzにおいて減衰はなく、通過域の中央のように見えます。-3dBの通過帯域幅は設計仕様通りの700Hz。帯域中央周波数も約700Hz。スロープは緩慢です。

カットオフ周波数狙い値500Hz

Hilbert FIR filter with a cutoff frequency of 500 Hz.

カットオフ周波数狙い値500Hzにおいて減衰はなく、通過域の中央のように見えます。-3dBの通過帯域幅は設計仕様通りの500Hz。帯域中央周波数は約700Hz。スロープは緩慢です。

フィルタ特性のまとめ

1. +45度と-45度のフィルタのステップ応答は鏡像関係にある。
2. +45度と-45度のフィルタの位相は90度シフトしている（明瞭ではない）。
3. カットオフ周波数通過帯域幅狙い値4.0kHzのフィルタが最もスロープが急峻である。
4. カットオフ周波数が低く通過帯域幅狙い値が狭くなるほどスロープが緩慢になり、設計が難しくなるように見える。

位相シフトのシミュレーション

+45度と-45度のフィルタの位相は90度シフトしているようですが、明瞭ではありません。そこで、シミュレーションにより確認しました。

2kHzのSin信号を生成し、+45度と-45度のFIRフィルタ（カットオフ周波数4.0kHz）に入力しました。計算にはscipy.signal.lfilterを用いました。S（青）は入力した2kHzのSin信号、S_Plus（橙）は+45度のFIRフィルタの出力、S_Minus（緑）は-45度のFIRフィルタの出力です。

esults of checking phase shifts by inputting 2 kHz Sin signal to plus and minus 45 degree Hilbert FIR filters (cutoff frequency 4.0 kHz).

S_Plus（橙）は125サンプル点目（2.6msec）までに入力信号Sに目視上重なります。カットオフ周波数4.0kHzのフィルタのインパルス応答は125サンプル点目までに目視上収束していることと呼応しています。+45度のFIRフィルタの位相シフトは0度でした。これに対して、S_Minus（緑）の-45度のFIRフィルタの位相シフトは-90度でした。±45度のFIRフィルタと言うより、0度と-90度のFIRフィルタと言う方が適切と思います。

RX_Summerオブジェクト（AudioMixer4_F32クラス）

役割

Audioミキサは信号加算器です。RFミキサの乗算器とは異なります。0度と-90度のFIRフィルタを通したIQ信号を加算し、イメージ信号を抑制することが役割です。

オブジェクトはゲインを持ち、ゲインを乗じてから加算します。2チャンネルのゲインは通常（1,1）であり、LSBモードの時に（1,-1）に切り換えます。-90度のFIRフィルタの出力に対して、180度の位相反転を行うことに相当します。通信モード毎のゲイン切り換えは、Mode.cppファイルに定義されています。

信号処理連鎖の更新処理

update()関数では以下の信号処理を行っています。

1. 128サンプルのデータブロックに以下の処理を一括適用
   ・最初のチャネルの信号（I信号）データを取得し、共有出力変数に格納し、ゲインを乗算
   ・次のチャネルの信号（Q信号）データを取得し、ゲインを乗算し、一時的変数に格納
   ・共有出力変数と一時的変数を加算し、結果を共有出力変数に格納
2. 次の信号処理オブジェクトへのデータ送信
3. メモリ・ブロックの解放

単純な加算演算のため、RX_Summerオブジェクトのupdate()関数では大量の積和演算は発生しませんが、乗算はCMSIS DSP Libraryのarm_scale_f32()関数を、加算はarm_add_f32()関数をコールし、128サンプルのデータブロックに一括適用しています。

イメージ信号抑制のシミュレーション

イメージ信号の抑制を確認するシミュレーションのために、Sin/CosのIQ信号を生成し、それぞれ+45度と-45度（0度と-90度）のHilbert FIRフィルタ（カットオフ周波数4.0kHz）に入力し、その出力を加算しました。

チューニング信号

チューニング周波数のUSB側の対象信号と想定したIQ信号を生成しました。下記IQ試験信号は2kHzで半時計方向に回転します。

IQ signal assumed as the target signal on the USB side of the tuning frequency.

このIQ試験信号をHilbert FIRフィルタに入力し、出力を加算した結果を下記に示します。

Result of inputting the IQ test signal to the Hilbert FIR filter and adding its outputs.

2kHzのIQ試験信号はフィルタの通過帯域内にあるため、フィルタのタップ数に依存した過渡応答経過後は位相シフトのみが適用されます。位相シフトによってIQ信号の位相が揃った結果、定常状態の振幅は加算によって2倍になっています。フィルタのタップ数は151タップですが、目視上は125サンプル点目（2.6msec）で収束しています。先に見たように、カットオフ周波数4.0kHzのフィルタのインパルス応答は125サンプル点目までに目視上収束していることと呼応しています。

イメージ信号

チューニング周波数より高い周波数のLSB側のイメージ信号と想定したIQ信号を生成しました。下記IQ試験信号は時計方向に回転します。

IQ signal assumed to be the image signal on the LSB side at a frequency higher than the tuning frequency.

このIQ試験信号をHilbert FIRフィルタに入力し、出力を加算した結果を下記に示します。

The result of inputting the image-assumed IQ test signal to the Hilbert FIR filter and adding its outputs.

位相シフトによってIQ信号が逆位相の関係になった結果、定常状態の振幅は加算によってゼロに収束しています。フィルタのタップ数は151タップですが、目視上は125サンプル点目（2.6msec）までに収束しています。

Keiths' SDR（K7MDL局Mike OM版）に予め実装されたHilbert FIRフィルタによって、2.6msec程度の遅延を伴ってイメージの抑制が図られることがシミュレーションによって確認されました。

Teensy（10）Keiths' SDRのInputオブジェクト

注記（2022/10/30）
Keiths' SDRのサンプリング周波数は48KHzであることが実機による実験で確認できました。

修正（2022/07/28）
信号処理連鎖を流れるデータの単位は128サンプルのデータブロックであることを認識し、修正しました。

Keiths' SDRのRx信号処理連鎖ブロック図

Keiths' SDR（K7MDL局Mike OM版）のRx信号処理連鎖をコードから読み解いて抜き出したもの（α版）を下記に再録します。

今回は、Rx信号処理連鎖の先頭のInputオブジェクト（朱色点線内★）を調べました。

Block diagram (alpha version) of the Rx signal processing chain and the Input object under consideration here.

Inputオブジェクト（AudioInputI2S_F32クラス）

Teensy Audio BoardからF32（32bit浮動小数点）型のステレオ・オーディオ信号（IQ信号）を取り込むオブジェクトです。Rx信号連鎖の起点になります。 なお、SGTL5000を搭載するディフォルトのPJRC社のオーディオ・ボードには、IQ信号が順不同になるバグがあるらしいのですが、その解決策等を含めたコミュニティの成果の詳細は別の機会に報告します。

Inputオブジェクトが、IQ信号入力に係わるハードウェアの初期化も請け負っていると思われます。IQ信号入力に係わるハードウェア（青）とソフトウェア（赤）の関係図を下記に示します。なお、この図はコードとi.MX RT1060データシートから読み取ったものであり、将来修正の可能性があります。

Diagram of the relationship between hardware and software involved in IQ signal input.

コンストラクタ

IQ信号2channel分のサンプリング周波数は96/2=48kHzです。サンプリングはMCU（i.MX RT1060）のPLLがマスターとなってオーディオ・ボードを管理します。Inputオブジェクトのコンストラクタはbegin()関数をコールし、その関数の中で以下のMCU（i.MX RT1060）の初期化を実行しています。

1. DMAチャネル・オブジェクトへのチャネルの割り当て
2. MCUとオーディオ・ボード間のPLL Clockを含めたI2S通信の設定
3. DMAチャネル・オブジェクトへのISR（割込みサービス・ルーチン）の割り当て

MCUの機能ブロック図を下記に示します。

Functional block diagram of the MCU (i.MX RT1060) and modules involved in the Input object (red dotted circle).

Inputオブジェクトに係わる朱色点線のMCUモジュールについて、以下にメモします。

eDMA

MCU（i.MX RT1060）のDMA（Direct Memory Access）モジュールは、eDMA（enhanced DMA）と称されており、TCD（Transfer Control Descriptor）を用いて柔軟なデータ転送を可能にしています。Keiths' SDRは、DMAチャネル・クラスのbegin()関数内で、IQサンプリングデータ転送のためのDMAチャネルをTCDに設定しています。eDMAが転送するデータの単位は、128サンプルのデータブロックです。I信号128サンプルとQ信号128サンプルを交互に転送しています。

MCU内の eDMAがCPUコア（Arm Cortex-M7）とは独立にIQサンプリングデータの転送を取り仕切るため、低速なI/O処理に対するCPUコアの空転は発生しません。Teensyが搭載するMCUのCPUコアは600MHzで動作する高速1コアであるため、CPUコアを空転させない実装が求められます。

I2S/SAI

MCU（i.MX RT1060）のI2S/SAI（Inter-IC Sound / Synchronous Audio Interface）モジュールがI2Sバスを介して、オーディオ・ボードからIQ信号サンプリング・データを受け取り、eDMAでメモリに転送しているものと思われます。クラス名のAudioInputI2S_F32は、この「I2S」を冠しています。

begin()関数内で、スレーブ動作のオーディオ・ボードに対してI2S/SAIに必要なClock信号を提供するために、CCM内のAnalog System PLLを設定しています。

CCM

MCU（i.MX RT1060）のCCM（Clock Controller Module）がAnalog System PLLの中のAudio / Video PLLを用いて、オーディオ・ボードに供給する様々なClock信号を生成します。

I2S/SAIに必要なClock信号は、①audio master clock、②bit clock、③bus clockの三つです。それぞれTeensy 4.1のピンに、①23：MCLK、②21：BCLK、③20：LRCLKとして割り当てられ、オーディオ・ボードに供給されます。③bus clock（LRCLK）がサンプリング周波数96kHz48KHzに設定されます。

NVIC

IQ信号サンプリング・データのDMA転送完で割込を発生させ、ISR（割込みサービス・ルーチン）の中のupdate_all関数によって信号処理連鎖を起動しているものと思われます。ただし、この関数の中身は下記1行のみです。

static void update_all(void) { NVIC_SET_PENDING(IRQ_SOFTWARE); }

CPUコア（Arm Cortex-M7）の割込コントローラNVIC（Nested Vectored Interrupt Controller）に対して、IRQ_SOFTWARE番号の割込を保留状態にセットしているだけです。IRQ_SOFTWARE番号にはsoftware_isr()関数を割り当てています。この関数が信号処理連鎖起動の割込処理本体になります。

ここでの「保留」は実行中断という意味ではなく、実行予約という意味で使われていると思われます。割込優先順位に応じて順番が回って来れば実行されると思います。IRQ_SOFTWAREの割込優先順位は208/255番と低く設定されています。ハードウェア割込を優先するためと思われます。

規模の小さい組込ソフトウェアでは、割込でフラグセットのみ行い、ポーリングループでフラグを見て割込処理を起動するといった実装方法を取ることがあります。割込処理を長くすると多重割込処理が複雑になるからです。ここでは、そのような処理をシステマチックに行っている印象です。

IRQ_SOFTWAREはIRQ70番のエイリアスのようですが、ARM v7-M Architecture Reference Manualでは、IRQ70番は「Reserved」になっていてユーザには解放されていないようです。使用されていない割込番号をHackして独自に活用しているのでしょうか・・・？。OSが載っている訳ではないため、将来、「Reserved」が「Used」になっても問題はないと思いますが、Teensy 5？にバージョンアップした時には要注意です。別のハードウェア・モジュールがIRQ70番を使用している可能性もあり得ます。

保留設定後に実行されるsoftware_isr()関数の中身は、信号処理連鎖を構成するオブジェクトに対する更新処理update()関数の順次適用の繰り返しです。Inputオブジェクトに対する更新は、I信号とQ信号の個別チャネルの更新です。チャネル更新の中では以下の処理を行っています。

1. F32型オーディオ・データのスケール正規化
2. 次の信号処理オブジェクトへのデータ送信
3. メモリ・ブロックの解放

メモリ解放と言っても参照カウンタを減じるだけで、メモリ・ブロックを参照する他の信号処理オブジェクトが存在する限り解放はされません。こうして安全に、メモリ・ブロック上のIQ信号データは信号処理連鎖を構成するオブジェクトに転送されて行きます。

サンプリング周波数の再考

悩み

Rx信号処理連鎖の後段のHilbert FIR Filterを調べているときに、フィルタの設計仕様と通過特性の確認計算値が合わないことに悩みました。確認計算した遮断周波数が設計仕様に対して倍異なるのです。

信号処理の骨幹を成すサンプリング周波数の間違いが原因として濃厚と仮説を立てました。Decimation（Down Sampling）等の絡みも疑ったのですが、遮断周波数が倍異なることに対してロジックが合いません。（Keiths' SDRではDecimationを採用していないと暫定的に結論付けました。）

気付き

コードを精査している中で、サンプリング周波数を言及する時に2つの視点があることに気付きました。ADCの視点と信号処理の視点です。

通常、「サンプリング周波数」と言えば、後者の信号処理の視点で言及していると思い込んでいました。しかし、Inputオブジェクトのコードでは、暗黙の中に前者のADCの視点で言及していることに気付きました。

確信

InputオブジェクトのISR（割込みサービス・ルーチン）の中で、128サンプルのデータブロックを格納するi2s_rx_bufferに対するIOサンプリン・データの充填状況を調べています。バッファの前半部が充填中か、あるいは後半部が充填中かを調べ、前半部が充填中の時のみupdate_all関数を起動しています。つまり、ADCサンプリング周波数96kHz48KHz/(128x2)サンプルに同期した割込みの中の2回に1回しかRx信号処理連鎖の更新を起動していません。これにより、信号処理視点のサンプリング周波数は48KHzになります。（煩雑ですが、信号処理の起動周波数は96kHz48KHz/128サンプル/2回=375187.5Hz、信号処理の起動周期は2.75.4msecになります。）

具体的に見てみます。バッファの前半部にI信号を格納し、後半部にQ信号を格納すると仮定します。割込がかかった後にバッファ充填状況を確認に行き、バッファ前半部に時刻nのI信号を格納中であれば、時刻n-1のIQ信号は揃っていることを意味します。時刻n-1のIQ信号をRx信号処理連鎖の後段に送致することが可能です。

IQ data acquisition timing diagram and sampling frequency of the signal processing chain.

逆に、割込がかかった時にバッファ後半部にQ信号を格納中であれば、時刻n-1のIQ信号はまだ揃っていません。次の割込を待つ必要があります。

もっとスッキリと実装できないものかとも思ったのですが、eDMAがI信号あるいはQ信号のどちらのデータをメモリに転送したかを知りようがないため、このような実装になっているものと思います。MCU（600MHz）とI2S/SAI（96kHz48KHz）の速度が大きく乖離しているために成り立つ実装とも思います。

Teensy 4.1のIOピン数は多いため、PJRC社のオーディオボードを2枚スタックしてI信号およびQ信号専用のADCにすれば、サンプリング周波数96kHzを容易に実現できると思うのですが、そのような議論はないようです。サンプリング周波数と同時にbit数も拡張したいとのことから、別のADCへの換装が議論されています。

結論（暫定）（決定）

PJRC社のオーディオボードを採用したKeiths' SDRのRx信号処理連鎖のサンプリング周波数は48KHzである・・・と思います。

（Keiths' SDRのソフトウェア開発に係わるOM諸氏には周知の暗黙知と思いますが、コードも含めて明示されている情報はまだ見つかっていません。）

悩み再び！？（2022/08/12）

信号連鎖ブロック図の後段に位置するFilterConvオブジェクトではサンプリング周波数96kHz48KHzでBPFが設計されていることが分かりました。もしかしたら、Hilbert FIR FilterはTeensyディフォルトのオーディオ周波数44.1kHz48KHzで設計されているのではないか？との疑念が浮上しました。44.1kHz48KHzで設計したHilbert FIR Filterを96kHzで使用すると帯域が約2倍に拡大します。しかし、後段のFilterConvオブジェクトで帯域を絞るため、実害は無く仕様のミスマッチに気付かない可能性があるのではなかと推測しています。最終的にはハードウェアを組んで、信号処理連鎖の起動周期を測定しないと白黒が付かない状況です。

Teensy（９）Keiths' SDRのRx信号処理連鎖

注記（2022/10/30）
Keiths' SDRのサンプリング周波数は48KHzであることが実機による実験で確認できました。

背景

Keiths' SDR（K7MDL局Mike OM版）の信号処理連鎖を学び中です。Mike OMの光速開発には付いて行けていません。最近、Mike OMはUSBオーディオとCATをTeensy 4.1に実装し、PC上のWSJT-XとUSBケーブル1本でつながるようになった模様です。本命のSDR信号処理に加えて、使い勝手を左右するユーティリティ機能を光速で実装するスキルには敬服させられます。現在は、オーディオIQ信号のbit数を上げるためにコーディック換装の検討をされているようです。

Mike OMが夏のバカンスに出かけている間に、出来るだけキャッチアップしたいと思います。遅ればせながら、Rx信号処理連鎖の解読に着手しました。

Keiths' SDRの信号処理実装方法

オブジェクト指向による実装

Keiths' SDRの信号処理は「OpenAudio_ArduinoLibrary」と称するクラス・ライブラリを使用して実装されています。拡張子cppのC言語実装ではなく、C++のオブジェクト指向で実装されています。

オブジェクト指向はアルゴリズムを考える抽象化の階層を切り換えられるため、クラスを設計した当人には便利なのだろうと思います。しかし、個人的な感想としては、低レベルのハードウェアの操作がどこで行われているか分かり難いと感じます。

例えば、

• コーディックのサンプリング周波数の設定はどのオブジェクトのコンストラクターで暗黙のうちに実行されるのか？
  （コンストラクターは全てのクラスが持っている）、
• サンプリング割込みはどこでトリガーされるのか？
  （実はタイマー割込みではなく、PLLで自動起動されるA/D変換処理の終点のDMA転送で割込みがトリガーされている模様）、

・・・等々のハードウェア処理を追うのが大変です。

IQ信号サンプリング周波数の設定

なぜハードウェア設定を気にしているのかと言いますと、Teensyのオーディオ信号処理はCDクオリティのサンプリング周波数44.1kHzをディフォルトにしているからです。開発元のPJRC社が標準コーディックSGTL5000用に提供するクラス・ライブラリ「AudioControlSGTL5000」はそのenable関数で定数44.1kHzをディフォルト引数にしています。

無線IQ信号処理に44.1kHzは低いため、Keiths' SDRではSGTL5000の仕様上限の96kHzに設定しているはずです。最初はどこで設定しているのか分からず、もしかしたら44.1kHzのままなのかと思ったほどです。Rx信号処理連鎖の先頭のオブジェクトのAudioInputI2S_F32クラスが怪しいと踏んでコンストラクターを追い駆けると、下記の通り4階層目にやっとハードウェアのi.MX RT1060チップ（MCU）のPLLに対してサンプリング・クロック96kHz48KHzを設定するset_audioClock関数が出現しました。

1. AudioInputI2S_F32クラスのコンストラクター
2. AudioInputI2S_F32クラスのbegin関数
3. AudioOutputI2S_F32クラスのconfig_i2s関数
4. imxrt_hw.cppのset_audioClock関数

SGTL5000はMCUのスレーブモードで動作する設定です。MCUのPLLからSGTL5000に96kHz48KHzのサンプリングClockを供給しています。混乱させることに、SGTL5000もPLLを持っており、マスタモードでも動作できるようです。

コンストラクター内で、InputクラスがOutputクラスの関数を呼んでいるというのが分かり難いですね。ADCとDACを備えたコーディックなのでどちらかで設定すれば良いということですが、Rx信号処理連鎖で動いている時間の方が長いことから、Inputクラスに設定関数があった方が分かり易い気がします。SDRラジオの場合だったとしたら、Outputクラスが出てくるのは変ですし・・・。

Audio Adaptor Board for Teensy 4.x

注記（2022/07/09）
コードの精査から、96kHz/1channel、48kHz/2channelのようです。ADCの視点では96kHz/1channelであっても、2channel分のオーディオIQ信号を取り込む必要のある信号連鎖の視点からはサンプリングタイムfs=48kHzになります。次回に詳細を報告します。

OpenAudio_ArduinoLibrary

SDR信号処理連鎖の実装に用いられているオブジェクト指向コアライブラリの「OpenAudio_ArduinoLibrary」について簡単にメモします。

Teensy Development Boardの開発元であるPJRC社が提供する「Teensy Audio Library」が継承ベースになっています。PJRC社はオーディオ信号処理をTeensyの主要なアプリ分野の１つに据えている模様で、専用のコーディック・ドーター・ボードも提供しています。Keiths' SDRはこれを標準？コーディックとして使用しています。

オーディオ信号処理アプリ開発の敷居を下げるために、PJRC社はTeensy Audio Libraryを使用したアプリをビジュアル設計できるAudio System Design Tool for Teensy Audio Libraryも提供しています。ビジュアルプログラミング用ツールのNode-REDにインスパイアされて開発した模様です。

しかしながら大きな制約として、Teensy Audio Libraryは「I16」（16bit整数）型でしか信号を扱えません。通常のオーディオ信号処理アプリにはこれで十分と判断したのかもしれませんが、無線IQ信号処理にはダイナミックレンジが不足します。Teensy 4.xが搭載するi.MX RT1060チップのArm Cortex-M7コアは、FPUとDSP extensionsによって「F32」（32bit浮動小数点）型の積和演算を1クロックで実行する能力があるのに、宝の持ち腐れになってしまいます。なお、DSP extensionsはMCUと別のDSPコアではなく、MCUの内部でDSP積和演算を実行可能にする拡張機能を指します。

そこで、TeensyユーザのOM諸氏がF32型の信号を扱える上位互換の継承ライブラリを開発して公開しました。それが「OpenAudio_ArduinoLibrary」です。「OpenAudio_ArduinoLibrary」のクラス名には全て「_F32」の識別子が末尾に付き、継承している元のTeensy Audio Libraryと混在して使用することができます。下記に示すように、ビジュアルプログラミングも可能です。

Rx信号処理連鎖

ブロック図

Keiths' SDR（K7MDL局Mike OM版）のRxの信号処理連鎖をコードから読み解いて抜き出したものを下記に示します。α版です。

Rx signal processing chain implemented using the OpenAudio_ArduinoLibrary.

着色したブロックが信号処理に係わり、無着色のブロックはスイッチ等の信号連鎖の管理に係わっています。ブロックを表すオブジェクトのインスタンス同士はpatchCord...関数が結びつけます。ビジュアル設計ツール（Audio System Design Tool for Teensy Audio Library）を使用すれば、自動的に生成してくれるようです。

ビジュアル設計ツール上にブロック図の一部をトレースして、コードのExportを実験してみました。

Experiment with tracing a portion of a block diagram on the Visual Design Tool and exporting the code.

Exported code from the Visual Design Tool.

Exportされたコードを見ると、必要なヘッダーファイルの宣言、OpenAudio_ArduinoLibraryのオブジェクトの宣言、オブジェクト間を連結するpatchCord...関数の宣言が自動生成されていました。ビジュアル設計ツール上の座標はコメントとして付加され、Import時に配置が再現されるようです。

次回より、着色した信号処理ブロックの中身を調べて行きたいと思います。Teensy4.1に搭載されているArm v7-Mベースのi.MX RT1060は複雑なMCUです。機能設定のレジスタの数も膨大です。ハードウェアがオブジェクト指向によって抽象化されていることもあり、まだ全貌を把握するには至っていません。分かった範囲、あるいは予想がついた範囲の内容をメモし、後で正誤が判明した場合は修正したいと思います。

AFP-FSK Transceiver（13）FSK送信試験

試験方法

FT8等のFSKモードのあい路の１つは、送信電波の品質チェックが難しい点を挙げられると思います。SSB帯域外は従来のスプリアス確認で間に合いますが、オーディオ帯域内の品質、特にオーディオ高調波を見逃すと他局に混信を与えてしまう可能性があります。WSJT-Xは「疑似スプリット」まで用意し、送信電波の品質確保に注力していることが分かります。

AFP-FSK Transceiverは原理的にオーディオ高調波を生じないはずですが、ソフトウェアの実装に依存して意図しない動作をしないか予め確認しておく必要があります。試験システムの系統図を下記に示します。前記受信試験の13TR-FT8とAFP-FSK Transceiverを入れ替えるだけです。Dummy Loadにつなげて送信した際の近接場漏洩電磁波を評価します。漏洩電磁波の強度を試験用最小限にするため、AFP-FSK Transceiverの電源は12VDCとしました。

Systematic diagram of the test system for transmitting FSK radio waves.

WSJT-Xチューン・モード試験

500Hz チューン

まず、WSJT-Xチューン・モードを試しました。指定したTx周波数のAUDIO信号をWSJT-XがAFP-FSK Transceiverに入力してくれます。Tx周波数を500Hzに設定してチューンをクリックしました。

Results of ASP-FSK Transceiver transmission test using WSJT-X tune mode.

一般論では、混合器等の非線形素子によって低周波の500Hzから発生したオーディオ高調波（1k、1.5k、2k、2.5k、3kHz）は狭帯域フィルタ（＜3kHz）を掻い潜って送信されてしまうリスクが高くなります。

しかし、AFP-FSK Transceiverからはオーディオ高調波に由来するRF高調波が一切出ていないことが確認できました。AUDIO周波数と搬送波周波数の混合はAFP-FSK演算の中で単なるスカラ値の加算として行われているため、原理的にRF高調波は出ません。その原理的帰結が念のため確認されたことになります。

気付き

送信開始時に周波数が少し右にシフトしています。周波数軸を拡大して、チューンを複数回試行した結果を下記に示します。

Minute frequency shift at the start of the tune.

間違いなく、右に約10Hzのシフトが繰り返し再現されています。WSJT-Xのチューン・モードの処理に依存しているのか、あるいはFSKの周波数遷移時の高調波を防止するWSJT-Xのフィルタ処理の副作用か、あるいはAFP-FSK TransceiverのAUDIO周波数計測の移動平均処理の遅れか、・・・原因は今のところ不明です。

AFP-FSK Transceiver送信と13TR-FT8受信

400Hz、1kHz、2kHz、4kHz AUDIO信号

WaveGenから400Hz、1kHz、2kHz、4kHzのAUDIO信号をAFP-FSK Transceiverに入力しました。AFP-FSK Transceiverがダミーロードに送信した電力からの漏洩電磁波をSDRplayで受信した結果を下記に示します。WaveGenの指示通り、周波数がシフトしていることが分かります。

Receiving results by SDRplay of the leaked electromagnetic wave transmitted by the AFP-FSK Transceiver to the dummy load show that the frequency is shifted as instructed by WaveGen.

この時、13TR-FT8が受信した漏洩電磁波のWSJT-Xワイドグラフを下記に示します。目視ですが、周波数シフトはSDRplayの受信結果と一致します。若干左にずれているのは、13TR-FT8の前記校正結果と一致します。

The WSJT-X wide graph received by 13TR-FT8 shows that multiple shifts of the AFP-FSK Transceiver's transmit frequency match the results received by SDRplay.

なお、AFP-FSK Transceiverはオーディーオ出力の帯域が許す限り可聴域を超えた4kHzでも送信可能ですが、国内で許可されるHFバンドの狭帯域データ・モードとするためには3kHz以下（副搬送波は2.95kHz以下）に制限する必要があります。

気付き

Windowsのオーディオ関係のAPIは、以下の順に進化してきたようです。

1. MME（Multi-Media Extensions）
2. Direct Sound
3. ASIO（Audio Stream Input Output）（他社規格）
4. WASAPI（Windows Audio Session API）

WaveGenの再生デバイス設定では、MMEを使用し、オプションの「EXTENSIBLEを使う」にチェックを入れていました。ヘルプでは「多チャンネル多ビット長対応の拡張フォーマットである WAVE_FORMAT_EXTENSIBLE 形式でデバイスをオープン」となっています。

この設定で13TR-FT8は動作していましたが、AFP-FSK Transceiverは送信ON/OFFが頻繁に切り替わるなど動作が不安定になり送信周波数も安定しませんでした。「EXTENSIBLEを使う」のチェックを外すと安定になりました。USBサウンドカードとの不整合が発生し、おそらくAUDIO信号が間欠的になっていたと現象から推定されます。アナログの13TR-FT8の方がこの点ではロバストでした。ディジタルのAFP-FSK Transceiverは過敏だったようです。

Message送信試験

SDRplay受信

次に、User messageの送信試験として、WSJT-XからAFP-FSK Transceiverに、1kHz副搬送波に対するCQメッセージのModulationを指示しました。SDRplayによるModulation漏洩電磁波の受信結果を下記に示します。

Receiving results of electromagnetic leakage by SDRplay when the modulation of CQ message for 1 kHz subcarrier wave was instructed to AFP-FSK Transceiver from WSJT-X as a transmission test of User message.

SDRplayのRefresh rateを最大にしてコントラストを調整すると、ウオーター・フォールに8-GFSK変調波が浮かび上がってきました。7,075kHzからの占有周波数帯域幅が50Hzであることが確認できます。この中に「CQ」が隠されているはずです。

上側のFFTスペクトルは瞬時周波数スペクトルを示している訳ではありません。40mバンドに対するSDRunoのFFTの仕様は、データ長65,536点、サンプリング周波数333,333Hz、平均回数3回です。FFTが対象にする時系列データ区間は0.59秒になります。FT8の送信速度は6.25-baudであることから、0.59秒には3.7回の周波数遷移が含まれていると思われます。これにより、前記チューン・モードよりも幅の広いスペクトルになります。

50dBcより下には、お椀形状のスペクトルが占有周波数帯域幅50Hzの外側に広がっていることが確認できます。さらに、お椀の左右外側にも50dBc以下のスプリアスが見えます。このスペクトルがどこから来るのか、必然なのか、意図しない漏洩なのかはまだ不明です。

FT8プロトコルの深耕

AFP-FSK Transceiverのようなシンプルなトランシーバによって従来にも増して遠方と交信できる秘密は、FT8のプロトコルの中に隠されているはずです。FT8がこれだけ普及しても、プロトコルの分かり易い解説は流布が少ないようです。Taylor博士らの比較的新しい文献*1が見つかりましたので、Message送受信のベースにあるFT8プロトコルの概要を調べました。

Processing steps of the FT8 protocol on which the message transmission test is based.

FEC: forward error correction
LDPC: Low Density Parity Check
CRC: Cyclic Redundancy Check
GFSK: Gaussian Frequency-Shift Keying
BP: Belief Propagation
OSD: Ordered Statistics Decoding
AP: a priori

プロトコルに盛り込まれた多くの技術は、惑星探査機の通信技術で実績があるようです。FT8の生い立ちを考えると、納得の行く話です。

一番の特徴はMessageの内容が約束事で固められている点と思います。これによって誤り訂正を施し易くなるとともに、無味乾燥の印象を与えることがあるようです。Taylor博士らの工夫が随所に盛り込まれていることに反比例して、ユーザの工夫点や運用スキルを入れる余地が小さくなり、運用体験も一様化され、無味乾燥となってしまうのかもしれません。自らの工夫点を入れて、WSJT-Xに代わるソフトウェアを作るのは自由なようですが、敷居は遥かに高くなります。

送信174-bitの中の97-bit（56%）は誤り訂正用符号です。通信の信頼性を上げる基本はレポート交換時のように2回繰り返すことですから、bit数が2倍に膨れていても驚くことではないのかもしれません。デジタルの恩恵により、2回繰り返すよりも大幅に信頼性が上がっているはずです。

13TR-FT8受信

AFP-FSK Transceiverが正しくMessageを送信しているかどうかを試験するためには、送信Messageを受信してデコードする必要があります。AFP-FSK Transceiver送信用のWSJT-Xと、13TR-FT8受信用のWSJT-Xの2系統が必要になります。そこで、ジャンク箱に眠っていたRaspberry Pi 3Bを発掘して送信系に用いることにしました。（レガシーRaspberry Pi 3Bを立ち上げるに際に、何も調査せずに試行錯誤をしたところ時間を浪費してしまったため、立ち上げの記録を付録に残します。）

Messageの送受信を検証する試験システムの構成と写真を下記に示します。

System configuration for testing whether the AFP-FSK Transceiver is correctly sending Messages.

Photo of the system used to test whether the AFP-FSK Transceiver is correctly sending a Message.

CQメッセージを送受信した試験結果を下記に示します。Raspberry Pi 3B上のWSJT-XにてエンコードされたCQメッセージがAFP-FSK Transceiverから送信され、そのCQメッセージは13TR-FT8で受信され、Windows PC上のWSJT-Xにて正しくデコードされました。

A CQ message encoded in WSJT-X on the Raspberry Pi 3B was sent from the AFP-FSK Transceiver, which was received by the 13TR-FT8 and correctly decoded in WSJT-X on the Windows PC.

気付きとして、受信信号のワイドグラフを見ると、左下に舌状のスペクトルが出ています。前記チューン・モード試験で気付いた現象と同様に、送信開始時点の周波数が正しく送信されていないか、あるいは受信されていない可能性があります。

気付きとSDRplayによる目視デコード

CQメッセージ送信開始時点の頭出しに注意して、再びSDRplayで変調信号を受信しました。今回は、Bias-TをOFFにして余分な漏洩電磁波を拾わないように注意しました。13TR-FT8をケースに入れた効果と相まって、前回のようなお椀形状の外乱スペクトルは影を潜めました。

Result of receiving the modulated signal in SDRplay again, including the start of CQ message transmission.

先に調べたFT8のプロトコル仕様から、エンコード・トーン・シンボルの頭には同期用の7トーンのCostas array（3,1,4,0,6,5,2）が付加されているはずです。周波数シフトを目で追うと、「3,1,4,0,6,5」トーンまでは確認できますが、トーン「2」が約30Hz右にオフセットしたと仮定しないと辻褄が合いません。Tx周波数が2.000kHzであったことと、占有周波数帯域幅50Hzを考慮すると、逆に「3,1,4,0,6,5」トーンまでは本来より左にオフセットしていたと考えられ、これが舌状スペクトルの正体と思われます。

Costas arrayの次は、3-bitのMessage Type（Std Msg Type=001）か、28-bitの無指定CQメッセージ（0000000000000000000000000010）が来るはずです。どちらが先に来るかは、前述のプロトコル解説書には記載が無いようでした。3bitsを束ねたトーン・シンボルは、Message Type（1）もしくはCQメッセージ（000000001..）になります。Message Typeが先に来るとすると、目視デコードが間違っていたことになり、舌状スペクトルの解釈も怪しくなります。逆に、CQメッセージが先に来るとすると、目視デコードは正しく、舌状スペクトルの解釈も正しくなります。

FT8プロトコルの机上エンコード

WSJT-Xのマニュアルを見ても、Message Typeがメッセージの前に付くのか、あるいは後ろに付くのか分かりません。しかし、マニュアルに次の記述がありました。

Additional utility programs jt4code, jt9code, and jt65code let you explore the conversion of user-level messages into channel symbols or “tone numbers,” and back again. 

（また、ユーティリティプログラムjt4code、jt9code、jt65codeにより、ユーザーレベルのメッセージをチャンネル記号や「トーン番号」に変換したり、また元に戻したりすることができます。）

残念ながら、ft8codeユーティリティは無いようですが、マニュアルの改訂を失念しているのかもしれません。ダメ元で打ってみると幸いにもUsageが出てきました。「トーン番号」への変換機能はft8code.exeでもサポートされているようです。

CQメッセージのエンコード結果を下記に示します。緑がメッセージ、青がFEC符号、赤が「トーン番号」です。

Complete encoding result of the CQ message. Green bits are the message bits, blue bits is the FEC bits, and red simbols are the "tone symbols".

赤の「トーン番号」の"Sync"は7トーンのCostas array（3,1,4,0,6,5,2）です。その次は（000000001..）なので、CQメッセージです。緑のメッセージの最後は（001）bitsであり、Message Type（Std Msg Type=001）を表していると思われます。Source encodeの要の1つであるMessage Typeが語尾に付くことが確認されました。

これにより、SDRplayの目視エンコードが正しいことが確認され、舌状スペクトルの原因がAFP-FSK Transceiverの周波数計測遅れ等によるオフセットにあるのではないかという仮説が濃厚になりました。

しかし、30Hzぐらいの出だしのオフセットの影響を受けずにデコードできている結果から、FT8のロバスト性に改めて感服する次第です。同期用のCostas arrayを前、中、後の3ヵ所に配置した効果でしょう。

残された課題

1. 送信開始時点の周波数オフセットの原因究明（例えば、移動平均処理の影響？）
2. バンド・モジュール追加によるHFオールバンド化

追記（2022/07/02）

課題１の送信開始時点の周波数オフセットの原因究明のためにコードを熟読しましたが、初期の周波数が10～30Hz程度低くなる原因は見当たりませんでした。オーディオ周波数の計測において、重み付き移動平均の計算をしている箇所はありますが、計算値をリップルの判定に用いているだけであり、オーディオ周波数の計測に反映されることはないようです。FT8のAFP-FSK計算だけでなく、強制送信試験ピンによる単純搬送波出力の際にも初期周波数の低下が見られるため、別の所に原因があると考えるのが妥当です。

次に仮説を立てたのは、PLLシンセサイザーMS5351MのPLLの周波数ロック遅延です。チャージポンプの時定数等は分からないのですが、PLLリセットコマンドをMS5351Mに送信すると負帰還動作を始め、幾何かの遅延は生じる可能性があります。PLLリセットコマンドで受信音に雑音が入るという話もWeb上に見られます。

FT8のFSKチャネルのシンボル送信速度は6.25baudであるため、周波数がずれているCostas arrayの最初の6トーンの送信時間は約1秒（0.96 sec）です。PLLのロックに1秒かかるのは遅いような気がしますが、30Hzは7MHzに対して4ppm程度なのでチャージポンプのチャージに時間がかかるのでしょうか。

AFP-FSK TrasceiverのPLLシンセサイザーの制御コードは、PLLリセットを頻繁に発生させないように注意深く組まれているようです。PLLの後段の分周器設定に乖離が発生しない範囲の周波数変更については、PLLリセットを回避するコードになっています。試算してみると、40mバンドにおいて、0～1kHzの副搬送波変更に対してPLLリセットなし、1.5kHzに変更する時にPLLリセットが発生し、1.5～3kHzの変更に対してPLLリセットなしとなりました。FT8送信中の50Hz範囲の周波数シフトでPLLリセットが発生することはなさそうです。今回の同じTx周波数に対する送信試験で舌状の周波数ずれが繰り返し見られたのは、PLLリセットが原因ではなかったことになります。

振り出しに戻ってしまいました。Clock信号を止めないでQCXのようにゲート制御する方法を試してみたいところですが、改良が大掛かりになってしまいそうです。

付録：レガシーRaspberry Pi 3B 再立ち上げの覚え書き

今回の試験のために、2016年頃に触っていたRaspberry Pi 3Bをジャックボックスから引っ張り出して再立ち上げしました。何も調査せずに試行錯誤をしたところ時間を浪費してしまったため、再立ち上げの記録を付録に残します。

回顧のRaspbian

まず、当時のOS（Raspbian）のセキュリティを向上するためにupdateやupgradeを繰り返しましたが、404エラーが多発してパッチが当てられません。どうやら、パッチファイルが既にサーバから撤去されている模様。

調べると、Raspberry Pi OSは最新版と1つ前のLegacy版しかサポートされないことが分かりました。

最新版のBullseye

素直に最新版のRaspberry Pi OSを焼き直すことにしました。公式のSDカード作成ツールImagerのトップに表示される現時点の最新版は通称「Bullseye」でした。何も考えずに「Bullseye」を書き込み、時間を掛けて初期設定を行いました。これが無駄な徒労でした。

現時点のWSJT-XのRaspberry Pi用コンパイル済みパッケージは、一つ前のLegacy版の通称「Buster」用しかありませんでした。一応「Bullseye」へのインストールをトライしましたが、どうにもこうにも多数のlib***のバージョンが不適合でインストールできませんでした。ソースからコンパイルするのもRaspberry Pi 3Bには荷が重そうです。

Legacy版のBuster

ここも素直にImagerの下の階層の中にある「Buster」を焼き直し、初期設定をやり直しました。「Buster」を焼き直した甲斐あって、WSJT-Xはトラブル無く簡単にインストールできました。

WSJT-X installed on Raspberry Pi OS "Buster".

その他

次にTime.isを見ると「あなたの時計はちょうどぴったりです」と表示。Raspberry PiはRTCを持っていないはずですが、インターネットに接続して時間が経てば自動的に時計を合わせてくれるようです。

USBオーディオカードは、（１）廉価版と、（２）96kHzに対応した高級版？の2種類を保有しています。（２）高級版は、Raspberry Pi 用のドライバが無く、上手く立ち上がりませんでした。IQ信号用に96kHz対応カードを準備していたのですが、オーディオ用は44.1kHzで十分です。（１）廉価版は、プルダウンリストに表示される名前がWindowsより多く特定が厄介でしたが、何とかつながりました。

*1:Steve Franke, Bill Somerville and Joe Taylor, "The FT4 and FT8 Communication Protocols - Motivation and design of the digital modes FT4 and FT8, and some details of how they are implemented in WSJT-X", QEX July/August 2020（https://physics.princeton.edu/pulsar/k1jt/FT4_FT8_QEX.pdf）

AFP-FSK Transceiver（12）FSK受信試験

試験方法

AFP-FSK Transceiverの前に組み立てた13TR-FT8（ディスクリート・アナログ・FT8トランシーバ）を活用した試験方法を思いつきました。13TR-FT8からダミーロードに送信して、近接場から漏洩した電磁波をAFP-FSK Transceiverで受信します。都市ノイズやPCノイズの海の中になりますが、自由にAUDIO周波数を設定できるため、ノイズとの切り分けも容易と予想しました。

試験システムの系統図と試験機材の写真を下記に示します。近接場の漏洩電磁波を捉えるためにSDRplayのアクティブ・ミニホイップ・アンテナはダミーロードの近くに置き、必要に応じてBias-Tをアクティブにしました。AFP-FSK Transceiverのモービル・ホイップもエレメントが近くに来るように室内の床に置きました。アンテナから送信はしないため、SWRは気にしないことにします。13TR-FT8はLOのトリマ再調査が必要になり、ケースから出す必要がありました。漏洩電磁波がどこからでも回り込みそうな系ではあります。

Systematic diagram of the test system for receiving FSK radio waves.

Photograph of the equipment for the test to receive FSK radio waves.

送信機13TR-FT8の周波数校正

まず、SDRplayを原器として、13TR-FT8の周波数を再校正しました。

校正前

両トランシーバの電源を投入すると、各々のLO（Local Oscillator）の信号が見えました。

When both transceivers were powered up, the respective LO (Local Oscillator) signals were visible.

AFP-FSK TransceiverのLO信号は、外部発振器をTCXOとしたMS5351Mの受信Clock信号です。一方、13TR-FT8のLO信号は、ノーマルXOによるコルピッツ発振回路の出力です。ディスクリートSBMの入力とするためか、13TR-FT8のLO信号の強度の方が高いようでした。もっとも、単にSDRplayのMini-Whip Antennaとの距離が影響しているだけかもしれません。

校正後

前回、13TR-FT8のLOも7,074,004Hzまで追い込んだはずですが、経時変化かあるいは暖機運転が足りないのか、約-40Hz程ずれているようです。これを±10Hz以内に再調整しました。

Readjust the LO of 13TR-FT8 to within ±10 Hz.

13TR-FT8送信とAFP-FSK Transceiver受信

2.5kHz AUDIO信号

WaveGenから2.5kHzのAUDIO信号を13TR-FT8に入力しました。13TR-FT8がダミーロードに送信した電力の漏洩電磁波をSDRplayで受信した結果を下記に示します。7,076.5（=7,074+2.5）kHzになる予定でしたが、約-2.5Hz程低くなりました。

Spectrum of electromagnetic leakage received by SDRplay when 13TR-FT8, with 2.5 kHz AUDIO signal input from WaveGen, outputs FSK power to a dummy load.

この時、AFP-FSK Transceiverが受信したWSJT-Xのワイドグラフを下記に示します。WSJT-Xは近接場の漏洩電磁波も捉えてくれました。目視ですが、約2.47kHzを受信しており、SDRplayの受信結果と一致します。約-2.5Hzは13TR-FT8のLOの校正誤差と思われます。

Wide graph of WSJT-X showing the results of the AFP-FSK Transceiver receiving the electromagnetic leakage when the 13TR-FT8 outputs FSK power to a dummy load with a 2.5 kHz AUDIO signal input from WaveGen.

2.3kHz AUDIO信号（-200Hzシフト）

次に、AUDIO周波数に対して2.5kHzから2.3kHzへの-200Hzの周波数シフトを行いました。7,076.3（=7,074+2.3）kHzになる予定でしたが、同様に約-2.5Hz程低くなりました。LO校正誤差を維持したまま、正確に周波数シフトした模様です。

A frequency shift of -200 Hz from 2.5 kHz to 2.3 kHz for the AUDIO frequency resulted in an accurate amount of frequency shift at the RF transmit frequency while maintaining the LO calibration error.

AFP-FSK Transceiverが受信したWSJT-Xのワイドグラフでも同じ結果を確認できました。AFP-FSK Transceiverの受信部が正常に稼働していることを確認できました。

The same results can be seen in the WSJT-X wide graph received by the AFP-FSK Transceiver, confirming that the receiving part of the AFP-FSK Transceiver is working properly.

AFP-FSK Transceiverによる交信受信試験

そのままAFP-FSK Transceiverを放置し、WSJT-Xの継続受信を行いました。室内床置き仮設のモービル・ホイップでは出力の大きい局やローカル局の電波しか捉えられないのですが、待っていると数局の電波を拾えました。ただし、DX（7.074MHz）も国内（7.041MHz）も交信相手局の電波は捉えられませんでした。

Receiving test results of QSOs on DX frequency and domestic frequency in the 40m band with AFP-FSK Transceiver.

DXにも国内にも両方に対応できるのがVFO Moduleを搭載したAFP-FSK Transceiverの利点です。受信用のVFO機能も正常に動作していることが確認できました。