dsPIC33FJ64GP802 (17) --- 周波数シフター (1)

　2018 年 11 月に開催された「アナログシンセ・ビルダーズ・サミット 2018」で、「公生32+」さんと「yama」さんが、「Bode Frequency Shifter」(以下「周波数シフター」と略記) の構成に基づいた (アナログ) シンセ・モジュールを出品されていました。
　この周波数シフターを dsPIC33FJ64GP802 を使って「ディジタル的」に実現しようと試みていましたが、「性能」や「実用性」はともかく、「周波数シフト」という基本機能については動作することが確認できました。
　内蔵 12 ビット ADC / 内蔵 16 ビット DAC を使用し、入出力の「シグナル・コンディショニング」用の OP アンプを除けば、主要なチップとしては dsPIC33FJ64GP802 ワンチップだけで実現できています。
　外部 16 ビット ADC を接続すれば性能の向上が期待できます。
　外部 ADC 対応、使用メモリ量の低減など、プログラムの整理がついたら公開したいと思っています。
　周波数シフターの構成について、簡略化したブロック・ダイアグラムをディジタル実現の場合について下に示します。

　「実信号」である入力信号を「ヒルベルト変換器」に通して元の信号と 90° の位相差を持つ虚部の信号を得て、「複素信号」である「解析信号」(analytic signal) に変換します。
　90° 位相差の 2 系統サイン波出力を持つ VCO 出力と解析信号とを複素乗算して出力信号を得ます。
　この VCO の発振周波数分だけ出力信号の周波数がシフトします。
　アナログ回路での実現の簡略化したブロック・ダイアグラムを次に示します。

　アナログ的にヒルベルト変換器を実現するのは難しく、Bode Frequency Shifter では Dome filter と称する Phase Shift Network (PSN) によって 90° 位相差の信号を得ています。
　90° 位相差というと、直接的には微分回路や積分回路で実現できますが、振幅が 20 dB/dec のレートで変化し、実用的ではないので、1 次移相回路の縦続接続により実現しています。
　1 次移相回路では、0° ～ -180° までの位相変化が得られますが、位相が変化する周波数帯域が限られているため、広帯域の信号に対応するには、それぞれコーナー周波数を変化させた 1 次移相回路を多段に接続する必要があります。
　多段に接続した移相回路を 2 系統用意し、両者の出力の位相差が 90° になるように各段のコーナー周波数を設計しなければなりません。
　製作にあたっても、特性が相互に干渉する各段のコーナー周波数を調整していく必要があります。
　「MOOG Bode Frequency Shifter Model 1630」では、1 次移相回路 5 段で構成されているので、合計の調整箇所は 2 × 5 = 10 箇所ということになります。
　その点、ディジタル・ヒルベルト変換器については等リプル近似の設計プログラムが使え、設計段階で性能がほぼ決まり、手動で調整する必要はありません。
　等リプル近似の FIR フィルタ設計プログラムとして有名な Parks-McClellan 法のプログラムにヒルベルト変換器の設計機能がついていたこともあり、その流れを汲んだ多くの FIR 設計プログラムにもヒルベルト変換器の設計機能が付いています。
　そのひとつとして、GNU Octave の remez() 関数で設計した例を示します。

# load "signal" package
pkg load signal;
# 254 tap FIR filter for Hilbert transform
N_tap = 254;
# sampling frequency
f_s = 24000;
# passband edge frequency
f_0 = 90;
# gain
gain = 1.0;
# weight
weight = 1;
# plot points
N_plot = 2048;
# Parks-McClellan optimal FIR filter design
hb = remez(N_tap-1, [f_0/(f_s/2), 1.0], [gain, gain], [weight], "hilbert");
# frequency plot
freqz(hb, 1, N_plot, f_s);

　これは、必要とするサンプリング周波数 f_s = 48 kHz の半分の周波数 f_s/2 = 24 k Hz を設計上のサンプリング周波数として、 90 Hz ～ 12 kHz までを通過帯域とする偶数タップ (254) のヒルベルト変換器を設計し、その周波数特性のグラフを freqz() 関数で表示するものです。
　通過域リプルがよく見えるように拡大したグラフを下に示します。 (図をクリックすると拡大します)

　これを、フィルタ係数には手を加えずに、2 倍のアップサンプリング、つまりサンプル間に「ゼロ」を補ってサンプル数を倍にすると、フィルタ係数はそのままなので周波数特性の形状は変化しません。
　ベースバンド領域の幅が 2 倍になるだけで、もともとはエイリアス領域だった部分がベースバンド領域に組み込まれることになります。
　具体的には、

サンプリング周波数 f_s/2 = 24 kHz で、

通過域下端 = 90 Hz

通過域上端 = 12 kHz

の偶数タップの FIR フィルタが、

サンプリング周波数 f_s = 48 kHz で、

通過域下端 = 90 Hz

通過域上端 = 24 kHz -90 Hz = 23.91 kHz

の奇数タップの FIR フィルタとなります。
　プロトタイプのフィルタのタップ数は 254 であり、そのタップ間に全部で (254 - 1) = 253 個の「ゼロ」を挿入するので、48 kHz サンプリングの (254 + 253) = 507 タップのフィルタとなります。
　507 タップのフィルタ係数のうち約半数の 253 タップ分の係数はゼロなので、入力データの遅延要素は必要ですが、タップの「計算」は必要でなく、手間を減らせます。
　また、奇数タップとなったので、解析信号の実部として入力信号をそのまま使うためのフィルタリング (実態は遅延時間合わせ) は値が「1」の「センタ・タップ」ひとつで実現できます。
　ヒルベルト変換器では、その定義上からフィルタ係数はセンタ・タップ位置を中心に「奇対称」となるので、係数はセンタ・タップの前側か後ろ側のいずれか片方あれば計算できます。
　それらを踏まえて、計算時間およびフィルタ係数の数を削減した、解析信号への変換フィルタをアセンブラで記述して使っています。
　WaveSpectra のリサジュー波形表示機能を使って固定周波数の正弦波を入力したヒルベルト変換器の入出力の位相差を観測した図を下に示します。

　90° 位相差を示す「真円」となっています。
　ヒルベルト変換器の 90° 移相特性はフィルタ係数の対称性であらかじめ決まっており、フィルタ設計の「近似」の対称ではありません。
　「等リプル近似」するのは出力振幅で、周波数スイープする正弦波に対するリサジュー波形を観察しているとリプルに応じて「傾く」のではなく、軸方向に伸び縮みするのがわかります。
　WaveGene + WaveSpectra でリニア周波数スイープ + ピークホールドにより観測したヒルベルト変換器の通過域リプルのようすを下に示します。　(図をクリックすると拡大します)

　上の図で、高域でリプルの下側のレベルが上がっているように見えるのは、周波数をログスケールで表示する際のデータの扱いによるものと思われ、拡大画面の方ではあまり変化がなく見えています。
　 10 kHz 以上でレベルが下がっているのは dsPIC 内蔵の DAC の特性によるものと思われます。