- Keiths' SDRのRx信号処理連鎖ブロック図
- FilterConvオブジェクト（AudioFilterConvolution_F32）
- カイザー窓 FIR フィルタのシミュレーション
- 悩み再び！？（2022/08/12）
Keiths' SDRのRx信号処理連鎖ブロック図
Keiths' SDR（K7MDL局Mike OM版）のRx信号処理連鎖をコードから読み解いて抜き出したもの（α版）を下記に再録します。
今回は、朱色点線内★でマーキングしたFilterConvオブジェクトを調べました。
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オブジェクトでした。
GUIの「Filter」ボタンに対するタッチイベントは、まずControls.cppファイルにあるFilter関数を呼んでいます。引数dirはフィルタの帯域を拡大する方向か、縮小する方向かを指示します。まだ確認できていませんが、タッチジェスチャを読み取っているものと思います。美しいGUIと合わせて、タッチイベント管理に関するMike OMの実装力には敬服するばかりです。
Filter関数は、Bandwidth2.cppファイルにあるselectBandwidth関数を呼んでいます。引数bndxは帯域設定番号であり、タッチ操作で増減させ、予め登録済みの帯域を逐次に切り換えているようです。帯域設定番号には、以下の帯域および中心周波数が割り振られています。CWモードの場合は帯域設定番号1から4の間で、SSBモードの場合は帯域設定番号5から9の間で切り換えているようです。これ以外に好みのフィルタ設定があれば、コンパイル前に増設や入れ換えも容易です。
- 帯域：250 Hz、中心：CWピッチ（可変）
- 帯域：500 Hz、中心：CWピッチ（可変）
- 帯域：700 Hz、中心：CWピッチ（可変）
- 帯域：1.0 kHz、中心：CWピッチ（可変）
- 帯域：1.8 kHz、中心：1.850/2 kHz
- 帯域：2.3 kHZ、中心：2.400/2 kHz
- 帯域：2.8 kHZ、中心：2.900/2 kHz
- 帯域：3.2 kHZ、中心：3.300/2 kHz
- 帯域：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処理を実行するためです。各状態で以下の処理を行っています。
- 状態=0：畳み込み計算と逆FFT、データ出力（1回目）
・周波数領域のデータとフィルタマスクとの畳み込み演算を実行。
（データ長は1024*2）
・畳み込み演算結果を逆FFTにより時間領域データに変換。
（データ長は1024*2）
・前半（0～511）のオーバラップ長512のデータを前回の結果と加算。
後半（1024～1535）のオーバラップ長512のデータを次回の準備として保存。
・最初の128（0～127）個の時間領域データのブロックを次の信号処理オブジェクトへ送信。
・メモリブロックの解放。 - 状態=1：データ出力（2回目）
・2つ目の128（128～255）個の時間領域データのブロックを次の信号処理オブジェクトへ送信。
・メモリブロックの解放。 - 状態=2：データ出力（3回目）、FFTバッファの準備
・3つ目の128（256～383）個の時間領域データのブロックを次の信号処理オブジェクトへ送信。
・メモリブロックの解放。
・FT用バッファの後半（1024～2047）をゼロで埋める。 - 状態=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タップ数固定により、スロープの傾斜は帯域設定に係わらず同じになります。
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によりフーリエ変換して周波数特性を求めました。下図の左からインパルス応答、周波数特性（横軸固定）、周波数特性（横軸拡大）を示します。帯域の広いフィルタから順にシミュレーション結果を示します。
その場フィルタ設計のシミュレーション結果
気付きを以下にまとめます。
- Hilbertフィルタとは異なり、インパルス応答（フィルタ係数）は左右対称。
- 周波数帯域を狭くするほど、インパルス応答（フィルタ係数）の幅は広くなり、振幅は小さくなる。
- dfc=250HzのCW用最狭帯域フィルタを除いて、-6dBcの帯域幅は目分量で設計仕様を満たしている。
- 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のフィルタ係数はコンパイル時にコードに埋め込んでおり、その場設計ではないため、変更していない可能性が高いのではないかと推定しています。
最終的にはハードウェアを組んで、信号処理連鎖の起動周期を測定しないと白黒が付かない状況になりました。
