非職業的技師の覚え書き

JK1EJPの技術的検討事項を中心に記録を残します。

AFP-FSK Transceiver(16)20mバンドモジュールの測定

QRPGuys AFP-FSK (Audio Frequency Processed - Frequency Shift Keying) Digital Transceiver III に標準で付いてくる40m/30m/20mバンドモジュールの中で最後に残った20mバンド(14MHz帯)モジュールの測定を行いました。

送信電力の測定

バンドモジュールを20m用に換装すると、バンドID抵抗の分圧値をADCで読み取り、周波数を自動的に14.074MHzに設定しました。バンドID抵抗の公差に問題は無いようです。

測定方法

送信電力はUSBオシロスコープ(AnalogDiscovery2)で測定しました。DC電源は13.8Vです。

Measurement system of transmit power by USB oscilloscope (AnalogDiscovery2).

測定結果

Measurement result of dummy load voltage by USB oscilloscope (AnalogDiscovery2).

ダミーロード電圧のAC RMS(実効値)から送信電力を計算すると6.3W、ピーク電圧から計算すると6.1Wになりました。40m/30mとほぼ同じ値です。E級増幅による送信電力にはハイバンド化の影響は出ていないようです。

スプリアス領域における不要発射の強度の測定

測定方法

tinySAを用いた「スプリアス領域における不要発射の強度」の測定系を下記に示します。

Measurement system for the strength of unwanted emissions in the spurious region.

送信出力6.3(W)は38.0(dBm)となるため、手持ちのMax. 41dBのステップアッチネータで39dB減衰して-1.0(dBm)をtinySAに入力しました。

tinySAはPC上のtinySA-Appから制御しました。周波数掃引範囲は10~75MHzとしました。分解能帯域幅RBWは自動設定とし、周波数掃引点数を変えて測定を繰り返しました。

測定結果

Measurement result of unwanted emission intensity in the spurious region by 500 frequency sweeps.

Measurement result of unwanted emission intensity in the spurious region by 1,000 frequency sweeps.

Measurement result of unwanted emission intensity in the spurious region by 3,000 frequency sweeps.

Measurement result of unwanted emission intensity in the spurious region by 10,000 frequency sweeps.

Summary of unwanted emission intensity measurements in the spurious region.

最も強度が高い不要発射は掃引3,000点時の第三高調波-51.3(dBc)でしたが、スプリアス規格を満たします。

第二高調波は平均-57.0(dBc)となり、無調整のウェーブトラップによって良く抑えられています。

スプリアス規格を満たしましたが、バンド周波数が上がると余裕が小さくなることが分かりました。次の17mは、第三高調波が境界条件上に比掛かるか、規格をオーバしそうな予感がします。経験上、第三高調波が規格を満たさない場合は、簡易なウェーブトラップの調整ではなく、LPFの段数を増やす改造が必要になります。

帯域外領域におけるスプリアス発射の強度の測定

測定方法

USBドングルSDRのSDRplay(RSP1A)を用いて「帯域外領域におけるスプリアス発射の強度」を測定しました。測定系統図と測定の様子を下記に示します。

Measurement system diagram of spurious emission intensity in the out-of-band region.

送信出力6.3(W)の38.0(dBm)に対して、20dBカプラで分岐し、手持ちのMax. 41dBのステップアッチネータと20dB固定アッチネータを通して減衰し、-43.0(dBm)の信号をSDRplayに入力しました。

測定結果

Measurement of spurious emission intensity in the out-of-band region.

SDRplayへの入力信号の測定強度は-44.4(dBm)でした。公称値による机上の概算予想-43.0(dBm)との乖離は1.4(dBm)でした。終段の温度上昇によって出力が減少する影響もあると思います。受信強度としては、S9+30dB弱の表示です。

帯域外±10kHzの領域のスプリアス発射の強度は-92.6(dBm)以下でした。帯域外領域におけるスプリアス発射の強度は基本波より48.2(dB)低い値となり、許容値に収まります。絶対値は1.2(mW)となり、こちらも許容値に収まります。

SDRplayとPCソフトウェアSDRunoの組み合わせによる20mバンドのRBWは最小5.93Hzになります。周波数軸を最大限拡大した結果を下記に示します。スペクトル中心のオフセットはRBW以下でした。

Expansion of the frequency axis around the fundamental frequency in the spectrum measurement result.

以上により、QRPGuys AFP-FSK Digital Transceiver III の3バンド(40m/30m/20m)化の準備が整いました。20mバンドモジュールも、LPFのトロイダルコイルについて巻き方等の調整は必要ありませんでした。ただし、40m/30mバンドよりも余裕が少なくなっているため、ハイバンドでどうなるか予断を許しません。

AFP-FSK Transceiver(15)30mバンドモジュールの測定

QRPGuys AFP-FSK (Audio Frequency Processed - Frequency Shift Keying) Digital Transceiver III に標準で付いてくる30mバンド(10MHz帯)モジュールの測定を行いました。

送信電力の測定

バンドモジュールを30m用に換装すると、バンドID抵抗の分圧値をADCで読み取り、周波数を自動的に10.136MHzに設定しました。

送信電力はUSBオシロスコープ(AnalogDiscovery2)で測定しました。DC電源は13.8Vです。

Measurement system of transmit power by USB oscilloscope (AnalogDiscovery2).

Measurement result of dummy load voltage by USB oscilloscope (AnalogDiscovery2).

ダミーロード電圧のAC RMS(実効値)から送信電力を計算すると6.0W、ピーク電圧から計算すると6.2Wになりました。40mとほぼ同じ値です。E級増幅QRPトランシーバに採用例の多い3並列FET(BS170)の終段回路として優秀な部類に入ると思います。

実効値は多点の数値積分に基づいて計測していると思いますが、ピーク電圧は最大値と最小値の2点から計測していると思います。大きな差はありませんが、実効値の方が信頼性が高いと思われます。

E級増幅スイッチングタイミングの測定

高出力の秘密はE級増幅のスイッチングタイミングにあると思われるため、その測定を試みました。下図の左にゲートソース間電圧Vgsと出力電圧V(ant)を示し、右にVgsとドレインソース間電圧Vdsを示します。

Measurement results of switching timing of FETs for class-E amplification.

Vgsは同じテストピンから同じ配線で測定しているにもかかわらず、リップルが大きく異なります。狭隘部から測定のための配線を引き回しているため、正確な測定ができていません。

ゲートのON/OFFのタイミングは想像で記入しました。測定配線に誘導されているリップルが大きいため心眼で判定するしかありませんが、Vdsがゼロのタイミングでスイッチングできているように思います。

スプリアス領域における不要発射の強度の測定

測定方法

tinySAを用いた「スプリアス領域における不要発射の強度」の測定系を下記に示します。

Measurement system for the strength of unwanted emissions in the spurious region.

送信出力6(W)は37.8(dBm)となるため、手持ちのMax. 41dBのステップアッチネータで減衰すると-3.2(dBm)となり、tinySAの安全入力範囲に収まります。

tinySAはPC上のtinySA-Appから制御しました。周波数掃引範囲は5~55MHzとしました。分解能帯域幅RBWは自動設定とし、周波数掃引点数を変えて測定を繰り返しました。

周波数掃引500点

Measurement result of unwanted emission intensity in the spurious region by 500 frequency sweeps.

周波数掃引1,000点

Measurement result of unwanted emission intensity in the spurious region by 1,000 frequency sweeps.

周波数掃引3,000点

Measurement result of unwanted emission intensity in the spurious region by 3,000 frequency sweeps.

周波数掃引10,000点

Measurement result of unwanted emission intensity in the spurious region by 10,000 frequency sweeps.

まとめ

基本波の強度は平均2.05(dBm)となり、事前想定からの差異は1.15(dBm)に収まりました。ステップアッチネータキットの抵抗値の公差から考えて、誤差の範囲と思います。

tinySAの周波数掃引点数を増やすと、自動設定のRBW(kHz)が小さくなることによって周波数解像度が上がり、ノイズフロアが下がります。強度の測定結果に系統的な変動(掃引点数に比例して強度が下がるなど)は見られないため、平均を取ってまとめても良いかと思われます。

「スプリアス領域における不要発射の強度」の測定結果を下表にまとめます。

Summary of unwanted emission intensity measurements in the spurious region.

最も強度が高い不要発射は掃引10,000点時の第三高調波-54.3(dBc)でしたが、スプリアス規格を満たします。経験上、第三高調波が規格を満たさない場合は、簡易なウェーブトラップの調整ではなく、LPFの段数を増やす改造が必要になります。

第二高調波は平均-56.2(dBc)となり、無調整のウェーブトラップによって良く抑えられています。

帯域外領域におけるスプリアス発射の強度の測定

測定方法

USBドングルSDRのSDRplay(RSP1A)を用いて「帯域外領域におけるスプリアス発射の強度」を測定しました。測定系統図と測定の様子を下記に示します。

Measurement system diagram of spurious emission intensity in the out-of-band region.

送信出力6(W)の37.8(dBm)に対して、20dBカプラで分岐し、手持ちのMax. 41dBのステップアッチネータと20dB固定アッチネータを通して減衰し、-43.2(dBm)の信号をSDRplayに入力しました。

測定結果

Measurement of spurious emission intensity in the out-of-band region.

SDRplayへの入力信号の測定強度は-44(dBm)でした。公称値による机上の概算予想-43.2(dBm)との乖離は僅かに0.8(dBm)でした。受信強度としては、S9+30dB弱の表示です。もっと減衰させた方が良いかもしれませんが、手持ちのアッチネータは全て投入してしまいました。

帯域外±10kHzの領域のスプリアス発射の強度は-92(dBm)以下でした。帯域外領域におけるスプリアス発射の強度は基本波より48(dB)低い値となり、許容値に収まります。絶対値は1.2(mW)となり、こちらも許容値に収まります。

SDRplayとPCソフトウェアSDRunoの組み合わせによる30mバンドのRBWは最小0.95Hzになります。周波数軸を最大限拡大した結果を下記に示します。

Expansion of the frequency axis around the fundamental frequency in the spectrum measurement result.

綺麗なスペクトルになっていると思います。スペクトルの中心は約6Hzオフセットしていました。0.59ppm程度の誤差です。TCXOでも±2ppm程度の変動はあるため、許容範囲です。

以上により、QRPGuys AFP-FSK Digital Transceiver III の2バンド(40m/30m)化の準備が整いました。30mバンドモジュールも、LPFのトロイダルコイルについて巻き方等の調整は必要ありませんでした。ただし、40mバンドよりも余裕が少なくなっているため、ハイバンドでどうなるか予断を許しません。

AFP-FSK Transceiver(14)マルチバンド化の準備

日乗

KCJコンテストの結果速報が届きました。KCJコンテストは、提出ログを照合し一致したQSOのみ得点を認める厳格な規約になっています。

当局は、コンテストの終盤にCQを出されている局を競合局のいないタイミングで呼ぶスタイルのため、QSO数は僅かです。その中の1局には再送を繰り返して辛うじて取って頂いたため、正確にログが一致するかどうか心配していました。当局のコールサインでは「E」が鬼門と思います。ノイズに埋もれてしまう場合は、スペースから逆算して「E」を感じ取って頂く必要があるかもしれません。

心配は杞憂に終わり、コールサイン、RST、マルチの互いのミスコピーは皆無でした。厳格な規約も振り返りには好適です。「相手のログが提出されていない」による1ポイント減点で済みました。ログ未提出だとコールサインが相手方にその旨伝達されることが分かりました。気をつけねば・・・。

プロローグ

しばらく工作から遠ざかってSDRソフトウェアの勉強をしてきましたが、猛暑も収まり換気のために窓を解放して工作を行うのに適した季節になって来ました。そこで、QRPGuys AFP-FSK (Audio Frequency Processed - Frequency Shift Keying) Digital Transceiver III のマルチバンド化に取り組むことにしました。

QRPGuys AFP-FSK Digital Transceiver III 

AFP-FSK Transceiver はバンドID抵抗を備えたバンドモジュール(送信LPFと受信HPF)を差し換えることで、HF帯(160~10m)のバンドQSYが可能です。標準で40m/30m/20mのバンドモジュールのPCBと部品が付属します。他のバンドモジュールの設計データも公開されているため、部品を集めれば製作可能です。

Band Module Specifications

当局の環境ではアンテナの実現が困難な160m/80mバンドは諦めます。60mは国内では許可されていません。標準添付の残りの30m/20mのバンドモジュールに加えて、17m/15m/12m/10mのバンドモジュールを製作し、合わせて7バンド化を試みることにしました。

部品集め

適切なユニバーサル基板が有れば、ロジスティックスが心配なUSPS経由でQRPGuysから部品を取り寄せなくても済みます。aitendoで塩梅の良い80×20穴のユニバーサル基板(UPCB80X20D)を見つけました。秋月でも同様の仕様の安価な基板がありましたが、コスト削減のためか4隅のRが付いていない仕様になっていました。バンドQSY時に素手で換装する基板になるため、今回はaitendoの基板にしました。

Universal PCB (UPCB80X20D) with 80 x 20 holes.

コンデンサは、LPFのウェーブトラップ調整用に、250V C0G の(2012サイズ)チップ積層セラミックコンデンサの系列をストックしているため、そこから適宜使用することにします。

最も入手が困難だった部品は、受信HPF用の固定インダクタ(1.2uHと680nH)です。近接放送局が無ければ、受信HPFは省略しても良いかとも思いましたが、他の部品と抱き合わせでDigi-keyに発注しました。

Fixed inductors for HPF of received signals that may not be distributed domestically.

配置設計

フリーソフトウェアとして公開されているPasS(Parts Arrange Support System)という名称のユニバーサルプリント基板エディタを使って、追加バンドモジュールの配置設計を検討しました。

PasSは、部品のビットマップイメージの端子ドットの色で配線の結合を判定していると思われるユニークな基板エディタです。トロイダルコイルのイメージが無かったため、似た寸法の電解コンデンサのイメージを流用しました。

Design of band module placement on a universal board.

配線の自由度はコネクタ端子に拘束されます。どうしてもGND配線と交差する信号配線(赤色)が1ヵ所生じるため、裏面配線(青色)だけでは収まりませんでした。

配置設計の過程で、QRPGuysオリジナルのバンドモジュールの配置設計に対して気付きがありました。同じように巻いた2つのトロイダルコイルのINとOUTの割り振りを変えていることです。つまり、IN1⇒OUT1⇒OUT2⇒IN2というパターン配線になっています。これによって磁束の方向が反転し、互いに漏れ磁束を打ち消し合うようになります。追加バンドモジュールの配置設計でも踏襲することにしました。回路図には表れないノウハウです。本来は設計思想は漏れなく図面指示すべきなのですが・・・。

LPFの事前シミュレーション

第二高調波がスプリアス規格を満たすようにするためには、ウェーブトラップをLPFに付加する方法が有効であることを経験しています。QRPGuys設計のLPFにはウェーブトラップが初めから付いています。

共振を利用するウェーブトラップでは、手持ち系列のコンデンサ単品で組んだだけでは、スプリアス規格を達成できる保証がありません。そこで、事前にシミュレーションによる容量感度の検討を行いました。ウェーブトラップのC21はチップコンデンサ2個並列まで許容する方針としました。

17m(18MHz)

LTspiceによる17mバンドLPFのAC解析結果を示します。ウェーブトラップのコンデンサC21は、22pF/33pF/39pFとして感度を調べまし。設計指示値は22pFです。

トロイダルコイルのインピーダンスは実測すべきですが、手持ちのLCRメータでは低インピーダンス領域になり、精度の高い測定は出来ませんでした。toroids.infoの計算値を用いています。

AC analysis result of 7m-band LPF by LTspice.

18MHzの第二高調波36MHzで共振するC21は、設計指示値22pFより大きな33pFであることが分かりました。

初段のコンデンサC24は「No C24」(実装不要)との指示がありました。これではLPFがチェビシェフ型にならない上に、E級増幅にも影響を及ぼしそうです。

LPFを搭載したバンドモジュールはコネクタ経由でメイン基板と結合するため、浮遊容量が隠れている可能性があります。それでも、1pF以下のオーダではないかと思われます。逆に、E級増幅の条件を最適化するために、LPFの仕様を二の次にしているのかもしれません。

15m(21MHz)

LTspiceによる15mバンドLPFのAC解析結果を示します。ウェーブトラップのコンデンサC21は、22pF/24pF/33pFとして感度を調べまし。設計指示値は22pFです。

AC analysis result of 15m-band LPF by LTspice.

21MHzの第二高調波42MHz付近で共振するC21は、設計指示値22pFより大きな24pF(12pF+12pF)であることが分かりました。

手持ちのコンデンサでは、22pFの次の系列は33pFになるため、12pFを2個組み合わせる必要があります。

12m(24MHz)

LTspiceによる12mバンドLPFのAC解析結果を示します。ウェーブトラップのコンデンサC21は、15pF/18pF/22pFとして感度を調べまし。設計指示値は15pFです。

AC analysis result of 12m-band LPF by LTspice.

24MHzの第二高調波48MHz付近で共振するC21は、設計指示値15pFより大きな18pFであることが分かりました。

10m(28MHz)

LTspiceによる10mバンドLPFのAC解析結果を示します。ウェーブトラップのコンデンサC21は、15pF/18pF/22pFとして感度を調べまし。設計指示値は15pFです。

AC analysis result of 10m-band LPF by LTspice.

28MHzの第二高調波56MHz付近で共振するC21は、設計指示値15pFより大きな18pFであることが分かりました。

事前検討のまとめ

総じて、ウェーブトラップのコンデンサC21を設計指示値よりも大きな容量にした方が良いという結果になりました。

わざわざ適正ポイントを外す設計をするとは思えないため、予想外に大容量の浮遊容量への対応等、設計の秘密が何か隠されているような気がします。

設計指示値通りに製作して、性能が未達となった時に、シミュレーションが示す適正値に変更する方針の方が良いかもしれません。

欠品 #26AWG UEW

コイルを巻き始めたのですが、QRPGuys支給のポリウレタン銅線(UEW線)はさすがに途中で費えました。手持ちのUEW線はΦ0.29mmとΦ0.5mmですが、指定は#26AWGでした。American Wire GaugeをUEW線の仕様で使う例を国内ではあまり見たことがありません。小分け売りのUEW線のパッケージにもAWGの表記はありません。

調べてみると、#26AWGはΦ0.4mm、安全電流0.32Aでした。5Wだと実効値で0.32A程度は流れそうなのでΦ0.29mmは危険と判断。Φ0.5mmだと巻いた時の曲げRが大きくなり、トロイダルコアとの密着が悪くなりそうです。設計指定通り#26AWG=Φ0.4mmを調達することにしました。しばし中断です。

Teensy(15)Keiths' SDRの"Twin Peaks"問題解決方法の調査

"Twin Peaks"問題とは

概要

SGTL5000コーディック(ADコンバータ、DAコンバータ)を搭載した、PJRC社のTeensy用オーディオアダプタには、起動時にIQ信号の片方の先頭データが欠落することがあるという問題が指摘されています。音声や音楽のステレオオーディオ信号を扱う場合は問題にならないようですが、IQ信号では大問題になります。これによって、それ以後のIQ信号相互の位相関係がシフトするため、IQ信号演算によるイメージ抑圧が正常に働かなくなります。その結果、搬送波の両側に信号ピークが立つ状況になるため、"Twin Peaks"と称されているようです。

この問題が存在することを知らないと、ある日突然、昨日まで順調に動いていたトランシーバの復調が出来なくなり、無暗にPCの再起動や再コンパイルを繰り返す等の都市伝説的対処に走ることになります。

幸いなことに、コミュニティにて原因の特定とSDRソフトによる解決策の造り込みが図られているため、調査した結果をここにまとめます。

経緯

  1. 2017/02/24
    Convolution SDRで著名なDD4WH局Frank OMが、PJRC社のフォーラムに「Reset audio board codec SGTL5000 in realtime processing」という名称のスレッドを立て、コーディックの問題提起をされたのが最初かと思います。
    Frank OMはこの時点で既に、IQ信号が欠落することが"Twin Peaks"の原因だと気付いておられます。mcHFトランシーバのコーディックWM8731でも同じ問題があり、コーディックのリセットで解決したことも記されています。コーディックSGTL5000でも同じ手段が取れるかどうかに問題が絞られています。
    人が気付いて対処するか(手動)、Teensyが気付いて対処するか(自動)が課題になります。
  2. 2018/04/16
    W7PUA局 Bob Larkin OMが議論に参加。様々なサンプリング周波数の設定で"Twin Peaks"問題が発生することを報告。
  3. 2018/05/22
    TeensyからSGTL5000に出力するClockのジッターが根本原因ではないかとの仮説を立て、これを減らすことを検討してスレッドの更新は止まっています。
    Clockのジッターが"Twin Peaks"の根本原因ではなかったようです。根本原因はTeensy側ではなくSGTL5000側に内包されているとして、対処療法で解決することが必要になりました。
  4. 2022/03/09
    W7PUA局 Bob Larkin OMからPJRC社のフォーラムの「 Floating-Point Audio Library Extension」スレッドに投稿があり、"Twin Peaks"問題を解決するAudioAlignLR_F32クラスをOpenAudio_ArduinoLibraryに追加したとのアナウンスがありました。提案された解決方法は後述します。
  5. 2022/03/10
    W7PUA局 Bob Larkin OMからkeithsdr@groups.ioに「Codec mis-aligned L-R timing」と題する投稿があり、"Twin Peaks"問題を解決するAudioAlignLR_F32クラスをOpenAudio_ArduinoLibraryに追加したとの同じアナウンスがこちらにもありました。
  6. 2022/04/09
    N0CTL局 John Scherer OMからkeithsdr@groups.ioの「PCB for Teensy + Teensy audio shield + RA8875」スレッドに投稿があり、Scherer OM設計のPCBに"Twin Peaks"の能動的自動検出のための回路パターンを搭載した旨の報告がありました。
  7. 2022/05/19
    K7MDL局Mike OMからkeithsdr@groups.ioの「Codec mis-aligned L-R timing」スレッドに投稿があり、AudioAlignLR_F32クラスを用いた"Twin Peaks"問題解決機能をKeiths' SDR(K7MDL局Mike OM版)ソフトウェアに標準搭載した旨の報告がありました。

以上で、ハードウェアとソフトウェア(ライブラリとアプリ)の両面で"Twin Peaks"問題は解決されました。こうして経緯をまとめると、集合知による課題解決絵巻のようです。

"Twin Peaks"問題の解決方法

問題

コーディックSGTL5000からのDMA転送データには、起動時にIQ信号の片方の先頭データが欠落することがあるという問題が存在します。Hilbertオブジェクトの調査の際にイメージ抑圧のシミュレーションに用いたPythonプログラムを流用して、具体的に何が起こるか見てみましょう。

I/Q信号の欠落シフトが無い健全な場合のイメージ抑圧結果を再録します。目的のIn-phase信号はHilbertフィルタ後段のSummerを通過し、Image信号はSummerで相殺されることが分かります。

Image suppression results for a healthy case with no shift due to missing I or Q signals.

I信号の欠落シフトが有る不健全な場合のイメージ抑圧結果を示します。目的のIn-phase信号はHilbertフィルタ後段のSummerを歪を伴って通過しています。ただし、その歪は目視で確認できるほど大きくはありません。

一方、Image信号はSummerで完全には相殺され得ないことが分かります。これにより、In-phase側とImage側の両方にピークが立つことになります。これがTwin Peaksです。

Image suppression results for the unhealthy case where there is a shift due to the missing I signal.

Q信号の欠落シフトが有る不健全な場合のイメージ抑圧結果を示します。Imageの位相が異なるだけで、同様にTwin Peaksが発生することが分かります。

Image suppression results for the unhealthy case where there is a shift due to the missing Q signal.

僅かにデータ1個分シフトしただけで、このような大きな影響が出ることが確認できました。

解決方法の原理

"Twin Peaks"問題が顕在化するタイミングは不定です。よって、SDR起動時に能動的に検出する方法が探られました。

IQサンプリングデータの欠落の有無を検出するために、起動時にコーディックの両L/Rチャンネルに同じテストデータを能動的に注入して比較する方法が採用されました。Teensy 4.x自体はDACを搭載していないため、サイン波を注入することは断念し、有り余るディジタル出力ピンの1つ(22番)から矩形波を注入しています。テスト時にRF信号を切断できると良いのですが、リレー回路等が必要になり複雑化するため、RF信号にテスト信号を重畳する方法が取られています。したがって、起動時に静かなバンドで"Twin Peaks"検出を行う必要があります。

サンプリング周波数の1/4の周波数の矩形波をコーディックに注入しています。下記に示すように、1周期に4点のサンプリングを行うことになります。コーディックからのデータ転送ブロック長は128ワードのため、32個の矩形波をサンプリングして判定します。

Normal conversion data of SGTL5000 for test square wave from Teensy's digital output pin 22.

Lチャンネルに対してRチャンネルのサンプリングシフト3までの相関関数を計算して、データ欠損によるシフトの有無を判定します。サンプリングシフトによってRチャンネルの最後の矩形波は欠損するため、相関関数の計算に使用するデータブロック長は124ワードとしています。31個の矩形波を使用することでノイズへの耐性を高めていることになります。

能動測定によって、LチャンネルもしくはRチャンネルのどちらかに欠落シフトがあることが検出された場合は、Rx信号処理連鎖に送出するデータブロックに整列補正シフトを施す対策を実施しています。コーディックのリセットは行っていないようです。SGTL5000コーディックがリセット命令を備えていないのかもしれません。

以下、サンプリングシフトの能動的測定方法の詳細を場合分けして記します。

Case 1:欠落シフトが無い場合

Case 1:  Normal L-channel data without missing and normal R-channel data without missing.

Lチャンネル(I信号)とRチャンネル(Q信号)の両方に欠落シフトが無い場合は、Rチャンネルのシフトが0の時の相互相関(正規化していないので相互相関係数ではなくベクトル乗算)xc[0]が最大の62になり、シフトが2の時のxc[2]が最小の-62になります。その他のxc[1]とxc[3]は0になります。ノイズ耐性を高めるために差分を取って正規化すると、xc[0]-xc[2]が1となり、他は0です。しきい値0.7と比較して、欠落シフト無し(In phase)と判定しています。

Case 2:Lチャンネルが欠落シフトしている場合

Case 2: Shifted L-channel data with one missing and normal R-channel data without missing.

Lチャンネル(I信号)が欠落シフトしている場合は、Rチャンネルのシフトが1の時の相互相関xc[1]が最大の62になり、シフトが3の時のxc[3]が最小の-62になります。差分xc[1]-xc[3]の正規化値が最大の1となり、Lチャンネル(I信号)の欠落シフト有り(Shift I)と判定しています。

Case 3:2つ分欠落シフトしている場合

Case 3: Shifted L-channel data with two missing and normal R-channel data without missing.

Lチャンネル(I信号)が2つ分欠落シフトしている場合は、Rチャンネルのシフトが2の時の相互相関xc[2]が最大の62になり、シフトが0の時のxc[0]が最小の-62になります。差分xc[2]-xc[0]の正規化値が最大の1となりますが、DNApplyと判定してエラーを返しています。DNApplyは「Do not apply」ということでしょうか。この場合は、サンプリングシフトの再測定を行うようです。

Case 4:Rチャンネルが欠落シフトしている場合

Case 4: Normal L-channel data without missing and shifted R-channel data with one missing.

Rチャンネル(Q信号)が欠落シフトしている場合は、Rチャンネルのシフトが3の時の相互相関xc[3]が最大の62になり、シフトが1の時のxc[1]が最小の-62になります。差分xc[3]-xc[1]の正規化値が最大の1となり、Rチャンネル(Q信号)の欠落シフト有り(Shift Q)と判定しています。

ハードウェアの要件

Teensyのdigital output pin 22から矩形波を出力して、コーディックSGTL5000のLINE INの2つのチャンネル(L/R)に同時に入力する配線が必要になります。インピーダンス整合のために、当初は100kΩ抵抗を挟んでいましたが、10kΩ抵抗が最適であることがMike OMによって見出されています。

Hardware wiring for active detection of "Twin Peaks".

N0CTL局 John Scherer OMの設計による「Teensy 4.1 Ra8875 Backpack with Audio Shield」(REV.3)ボードには、"Twin Peaks"の能動的測定のための配線パターンが準備されています。

信号の名称は「I2S_COR」となっていますが、I2S信号ではないことに注意する必要があります。I2S通信プロトコルに乗るIQ信号(L/R信号)の相関を能動的に調べるためのディジタル出力信号という意味からの名称と思います。

公開されているBackpack PCBのガーバーファイルから、裏面のパターンを可視化して下記に示します。黄色の文字は当局の注記です。"Twin Peaks"の能動的測定のための配線パターンが確認できます。

Backside wiring of the Backpack PCB.
(You can see the wiring from the Tennsy 4.1 digital output to the audio adapter's LINE L/R.)

組立中のBackpack PCBの裏面写真を下記に示します。透視で表面から見たガーバーファイル可視化イメージとは鏡像関係にあることにご注意ください。

Photo of the backside of the Backpack PCB during assembly.

信号配線は表面にあり、本来、裏面にはバイパスコンデンサとプルアップ/ダウン抵抗が配置されているだけだったようですが、そこに重要な"Twin Peaks"能動的検出パターンが追加配置されたようです。

ソフトウェアへの実装

"Twin Peaks"の能動的測定のためのハードウェアの仕掛けに合わせて、ソフトウェアの準備も必要になります。Keiths' SDR(K7MDL局Mike OM版)における実装を見ていきたいと思います。

W7PUA局 Bob Larkin OMが開発したAudioAlignLR_F32クラスのTwinPeakオブジェクトを、K7MDL局Mike OMがRx信号処理連鎖の先頭に挿入しました。このオブジェクトによって、先頭データの欠落シフトの能動的測定と、必要に応じて後段のRx信号処理連鎖に送出するデータブロックの整列補正シフトを実施しています。

Insertion of a TwinPeak object at the beginning of the Rx signal processing chain.

必要な信号処理オブジェクトをコード上で結線して行くPJRC社の作法を取り入れたKeiths' SDRのコーディングスタイルによって、TwinPeakオブジェクトのようなオプション機能を提供するオブジェクトの後付け挿入が容易になっています。コンパイラディレクティブ("W7PUA_I2S_CORRECTION")によって、上記ハードウェアの要件を満たしている場合はTwinPeakオブジェクトへの挿入結線をコンパイラ前処理に指示し、要件を満たしていない場合は点線の迂回結線を指示します。

mainプログラムに対するTwinPeakオブジェクトの動作を下記に示します。左側がsetup()関数とloop()関数から成るArduinoのmainプログラムを、右側がTwinPeakオブジェクトを示します。TwinPeakオブジェクトは、「計測モード」と「実行モード」の2つのモードを持ちます。

TwinPeak object behavior.

setup()関数からTwinPeakオブジェクトを「計測モード」に設定します。この時、上述の矩形波をTeensyのdigital output pin 22から出力します。タイマー等を用いることなく、mainプログラム内で経過アイドル時間を監視して、所定の周期になるように出力を反転させています。Inputオブジェクト経由でSGTL5000コーディックからI2S/DMA転送されてきたデータブロックに対して、欠落シフトの有無/種類を検出し、補正方法を選択します。「計測モード」が完了したら矩形波出力を終了して、TwinPeakオブジェクトを「実行モード」に切換えます。

mainプログラムがloop()関数に制御を移すと、SGTL5000コーディックからのI2S/DMA割込みによってRx信号連鎖を構成する各信号処理オブジェクトのupdate()関数が順番に起動することを、所定のサンプリング周期で繰り返します。「実行モード」にあるTwinPeakオブジェクトのupdate()関数は、「計測モード」で選択された欠落シフト補正をデータブロックに施し、後段の信号処理オブジェクト(I_Switch/Q_Switch)に送出します。

Teensy(14)Keiths' SDRのOutputオブジェクト

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

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

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

Rx信号処理連鎖の調査も佳境です。今回は、朱色点線内★でマーキングしたOutputオブジェクトを調べました。

前段のAm1_LおよびAm1_Rの調査も行いましが、モノラル信号を2分岐してスカラー倍しているだけです。左右で音量バランスを調整したい場合には有用です。出力段で16bit整数に丸めることを考えれば、スカラー倍のゲインはできるだけ大きくして、アナログAFアンプで音量を調整した方が良いように思います。有限bit長の制限から、SDRでは飽和に対してゲインを上手く調整する必要が生じるかもしれません。

Outputオブジェクト(AudioOutputI2S_F32)

端的に言えば、Inputオブジェクトと逆の処理を行います。ディジタル信号処理したオーディオ信号をコーディックのDACによってアナログ信号に変換するオブジェクトです。

MPUの負荷の時間的偏在によってオーディオ信号のストリームが途切れることを防ぐためか、数々の工夫が見られます。しかし、その設計思想が明示されていないため、解読が難しいオブジェクトでした。

ひとまず、解釈した結果を暫定版として記すことにします。暫定版のデータフローを下記に示します。

Data flow to output audio data (tentative version).

コンストラク

Outputオブジェクト宣言時に、暗黙のうちに実行されるコンストラクタの中でbegin()関数をコールして、ハードウェアの初期化を行っています。MCU(i.MX RT1060)のハードウェア・モジュールを下記に示します。

Hardware module of MCU (i.MX RT1060).

この中で特に重要なモジュールはeDMA(enhanced Direct Memory Access)と思います。ソフトウェアの世界とハードウェアの世界の境界を跨ぐ橋渡しの役割をeDMAが担っています。Rx信号処理アプリが出力したオーディオデータのコーディックへの転送を、CPUの負荷に関係なく独立にeDMAが実行します。

コンストラクタの中のbegin()関数では、以下の順番で初期化を実行しています。

eDMA(1)

最初に、DMAChannelオブジェクトのbegin()関数をコールして、eDMAにチャネルの割付を行っています。

ハードウェアであるeDMAモジュールのTCD(Transfer Control Descriptor)レジスタと対になるソフトウェア側のTCD構造体を用いてレジスタの設定を行っていますが、その設定は複雑なため追跡は断念しました。元々はTeensy販売元のPJRC社が開発したコードを踏襲していると思います。

コーディックに転送する出力データバッファi2s_tx_buffer(16bit整数x2、128ワード)のRAM上のアドレスをTCDに設定しています。出力データバッファi2s_tx_bufferは、Rx信号処理アプリが宣言したRAM上の変数です。アプリの変数の任意のアドレスをTCDに設定することで、アプリのデータをコーディックに転送することが可能になります。この設定が、ソフトウェアとハードウェアの境界の橋渡しになります。

I2S/SAI

Outputオブジェクトのconfig_i2s()関数をコールして、MCUのI2S/SAI(Inter-IC Sound / Synchronous Audio Interface)モジュールの設定を行っています。この関数はサンプリング周波数96kHzを引数に取り、コーディックに供給するClockを生成するMCUのPLLの設定を行っています。I2S/SAIのレジスタ設定も複雑なため追跡を断念しました。非同期通信を行っているらしいことはコメントから分かりました。

確か、Inputオブジェクトのコンストラクタでも同じconfig_i2s()関数をコールしていたはずです・・・? メモリ割り付けを行っていないなら、同じ関数を2回コールしても副作用はないと思いますが、分かり難い印象は残ります。

eDMA(2)

eDMAを起動するハードウェア・イベントとして「DMAMUX_SOURCE_SAI1_TX」(=20)を設定しています。MCU(i.MX RT1060)のデータシートと照合すると、Channel 20番はSAI1(Synchronous Audio Interface 1)モジュールの「SAI TX FIFO DMA Request」に割り付けられています。

DMA転送の要求はコーディックから出されると思っていたのですが、特に該当する信号線は見当たりませんでした。そこで、SAI1(SAI TX FIFO > コーディック)が送信の対象とするSAI TX FIFOバッファが空になると、自動的にeDMA(RAM上のi2s_tx_buffer  > SAI TX FIFO)を起動しているのではないかと推定しています。

NVIC

元々PJRC社が開発したAudioStreamクラスのupdate_setup()関数をコールして、割込コントローラNVIC(Nested Vectored Interrupt Controller)に、Rx信号処理連鎖の更新処理関数software_isr()をソフトウェア割込ベクタIRQ_SOFTWAREに設定しています。

eDMA(3)

DMAChannelオブジェクトのattachInterrupt()関数を用いて、DMA割込に割込サービスルーチンAudioOutputI2S_F32::isr()を割り当てています。eDMAの転送(RAM上のi2s_tx_buffer > SAI TX FIFO)が完了すると、割込が自動的に発生し、Outputオブジェクトのisr()を実行します。アプリが管理する出力データバッファi2s_tx_bufferを前半と後半に分け、それぞれの転送完了で割込を発生させているようです。次のDMA転送データをデータバッファi2s_tx_bufferに準備するのが、この割込サービスルーチンisr()の役割です。isr()の実行内容の詳細は、節を変えて記します。

コーディックはマスタのMCUに対してスレーブの位置付けになっています。MCUのSAI送信(SAI TX FIFO > コーディック)完了(①)からプルされてeDMA転送(RAM上のi2s_tx_buffer > SAI TX FIFO)が起動(②)し、eDMA転送完了(③)からプルされて割込サービスルーチン(RAM i2s_tx_bufferへの次回出力データの充填)が起動(④~⑧)する複雑な流れが見えてきました。

Diagram of the relationship between hardware and software involved in LR audio signal output.

ソフトウェアによる明示的な制御とハードウェアによる暗黙的な制御が絡み合っているため、全体の流れが分かり難くなっています。MCUのデータシートの理解が必要になりますが、データシートは可能なオプションの列挙に等しく、特定の機能の実現方法を理解するようには書かれていません。MCUが複雑で実施できるオプションが多過ぎることが理由の1つかと思います。SDRをコーディングしたOM諸氏も、Teensy開発元PJRC社のコードを可能な限り踏襲して拡張しているようです。

DMA割込サービスルーチンAudioOutputI2S_F32::isr()

  1. eDAMのTCDが指し示すメモリ上の現時点のDMA転送データのアドレス(SADDR:Source Address)を取得します。
  2. 取得したSADDRから、DMA転送用の出力データバッファi2s_tx_bufferの前半を転送中かどうかを確認します。
  3. 現時点で後半を転送中であれば、前回のデータブロックを転送途中と判定し、前半の先頭アドレスを次回充填先アドレスに設定する処理のみを行います。
  4. 現時点で前半を転送中であれば、前回の前半と後半のデータブロックを全て転送完了したと判定し、後半の先頭アドレスを次回充填先アドレスに設定します。
    続けて、update_all()関数をコールして割込コントローラNVICに対してソフトウェア割込ベクタIRQ_SOFTWAREを保留(PENDING)に設定します。この処理はInputオブジェクトでも実施していました。追跡できていませんが、InputとOutputの両方が完了するまで、更新処理を保留(PENDING)にしているのかもしれません。
  5. DMA転送用出力データバッファi2s_tx_bufferの次回充填先アドレスに出力データを充填します。この時に、32bit浮動小数点データを16bit整数データに型変換し、LチャンネルとRチャンネルのデータを合わせて32bitバッファに充填します。
  6. CPUのデータキャッシュをフラッシュして、キャッシュの内容をメモリに書き込みます。DMA転送の対象はメモリだからです。
  7. 出力データの前半をDMA転送用出力データバッファi2s_tx_bufferに充填したかどうかを調べます。
  8. 出力データの前半の充填が完了しているなら、出力データの後半の充填を次回に備えて設定します。
  9. 出力データの後半の充填が完了しているなら、次回に備えて、予備の2ndデータブロックblock_left_2ndとblock_right_2nd(のアドレス)を出力データblock_left_1stとblock_right_1st(のアドレス)に複写し、2ndデータブロック(のアドレス)をNULLに設定します。
    このような複雑な構成にした設計思想は不明ですが、データブロックを1stと2ndの2段構成にするのは、MPUの負荷の時間的偏在によってオーディオ信号のストリームが途切れることを防ぐためと推測しています。

update()関数による更新処理

上記DMA割込サービスルーチンのステップ4において、前半と後半のデータブロックを全て転送完了した場合は、更新処理を実行するsoftware_isr()関数を割り付けたソフトウェア割込ベクタIRQ_SOFTWAREを保留設定していました。保留が解ければ、全てのRx信号処理連鎖が更新され、Outputオブジェクトも更新されます。その更新内容を定義しているのが、このupdate()関数になります。

update()関数では、Rx信号処理連鎖の前段(オーディオ・ゲイン処理)の処理結果のデータブロックをQueueから取得します。LチャンネルとRチャンネルのデータを別々に処理します。16bit整数のスケール変換を行い、Double bufferの2ndデータブロックに格納します。スケール変換を行っただけで、型は32bit浮動小数点です。

update()関数の役割はここまでです。Double bufferから先のデータ転送は、前述したDMA割込サービスルーチンAudioOutputI2S_F32::isr()の役割になります。

感想

以上、Outputオブジェクトの全貌を理解するには、MCUのデータシートのレジスタ設定に精通する必要があります。SDRの改良や実験を行う上で、当面、それは必須ではないと思います。

RX信号処理連鎖のInputオブジェクトとOutputオブジェクトはファームウェアと見做し、その他をアプリケーション・ソフトウェアと見做して、信号処理アプリの改良実験に専念するのが良いかと思います。

 

RX信号処理連鎖の調査は今回で終了です。次回は、"Twin Peaks"問題の調査結果をまとめたいと思います。

 

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がオーディオ信号処理の市場を主要なターゲットの1つにしていたことが分かります。下記に分かり易い解説が有ります。

疑問

先に調査した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オブジェクトの初期化がコメントアウトされ、代わりに上記のフィルタ仕様の1つを設定した後に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フィルタ係数はもはや不要ということのようです。設計の一般論としては、1つの要素に2つの機能を割り付けるのは悪手とされています。二兎を追う者は一兎をも得ず、2つの機能の同時最適化が困難になることから導かれた教訓です。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の2つを引数として取り、これらをオブジェクトの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の優位点が出る機能の1つと思います。都市ノイズが増えた今日、自環境に合わせて最適化すれば有用なオブジェクトの1つになるかもしれません。

役割

自動ノッチ処理を行う信号処理オブジェクトです。前回の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オブジェクトはハードコーディング・オブジェクトの1つでした。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フィルタ係数の学習過程を下記に示します。左(1)は全1024ステップ、右(2)は最初の100ステップの拡大です。

Variation of FIR filter coefficients by adaptive learning.

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

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

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

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

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

左(1)の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.

左(1)のトーンノイズの自己相関関数の極小値と極大値の位置が、FIRフィルタ係数の極小値と極大値の位置と近いことが分かります。右(2)の音声信号はその限りではありません。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.

切り出し区間を長くした結果、音声信号には2つの音素が含まれているように見えます。周波数スペクトルは分解能が上がりましたが、低域の挙動は自動ノッチの試験データに似ていると思います。高域は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フィルタ係数の学習過程を下記に示します。左(1)は全4096ステップ、右(2)は最初の100ステップの拡大です。

Variation of FIR filter coefficients by adaptive learning.

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

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

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

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

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

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

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

左(1)の白色ガウスノイズの自己相関関数がデルタ関数になっていることが分かります。ただし、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にも実装できそうですね。