21世紀仕様の正弦波発生装置を作る・デジタルオシレータの製作（１）

やりたいこと

20世紀の回路設計を現代流に再構築し、個人レベルで実現可能な最上級スペックを狙った正弦波発生装置（オシレータ）制作にチャレンジしてみようというプロジェクトの一環です。

ベンチマークとするのはオーディオなど低周波・音声帯域系のアナログ回路設計の現場で定番だった VP-7201A という古の伝説的オシレータ。

レトロ計測器ではあるものの、5Hzから500kHzの正弦波を、0.005% 以下という低歪み率（ただし、この歪率を満足するのは 20Hz～20kHz の可聴帯域のみですが）で発生させることのできる、かなり魅力的な性能を持つ信号発生器です。

VP-7201A は状態変数型RC発振器（State Variable Oscillator: SVO）の構成を採用しており、これを今回、現代流のリファインを加えつつ、同じ方式で正弦波発生器を作ってみよう、というのが本プロジェクトの趣旨です。

一方、現代流のアプローチは DSP や FPGA を使ったデジタル信号処理技術を使う方法で、純アナログ方式より高精度・低歪み・広帯域、かつ安定した信号を、職人芸的な調整無しで得ることができます。
アナログ方式では（無線帯域の特殊分野の発振器は別として）どうやってもこの方式には太刀打ちできないのですが、それでもどこまで迫れるか見定めたいという、完全に自己満足を満たすためだけの、趣味の電子工作となります。

とは言え、大幅なモディファイを加えるにせよ昭和の回路をただ組むだけでは面白くありません。比較用のベンチマークとして、このデジタルオシレータもこの際作ってしまおう……というのが、今回のお題となります。

目標スペックと基本構成

VP-7201Aと競わせるのですから、当然のことながらスペックはこれと同じ

• 周波数: 5Hzから500kHz
• 生成波形: 正弦波
• 歪み率（THD）: 0.005% 以下

に設定します。
アッテネータ回路は状態変数型仕様のものと同じ回路を使用するので、出力仕様もこれに合わせます（出力電圧約 5Vrms ）。

実現方法（構成）ですが、生成周波数は500kHzオーバーを目指しますので、全くもってマイコンの領域ではない。かといって、個人の趣味のレベルでの製作ですので、あまり高価なデバイス、例えば高性能な DAC は使いたくないです。
そこで、基本構成はこうすることにしました。

• FPGAによるSoC構成
  • RISC-Vプロセッサ上のソフトウェアでマネジメント
    周波数設定用UIの提供と状態方程式パラメータ算出。
  • FPGA上のロジック回路で積和演算器を構成
    2次IIRフィルタを構成。$f_S$は 500kHz×40 で 20MHz あたりを目安とする。
  • 32ビット符号付固定小数点
    整数1ビット＋少数31ビットのQ1.31フォーマットとする。
• ΔΣ変調方式（DSM/PDM）の1ビットDACで波形出力
  • 3次 1-bit デルタ・シグマ変調器を採用
    ノイズシェーピングを活用しノイズを低減。
  • LVDS出力による差動出力
    後段のアナログLPF回路に差動信号を供給しSN比を稼ぐ。
    高速信号とすることで数百MHzオーダーのビットストリーム出力に対応。
  • PDMのクロック（キャリア周波数）は 200MHz 程度を狙う
    正弦波500KHzの400倍とし、ノイズ周波数を高周波側に追いやる。
• 外部アナログLPF
  • 初段パッシブRCフィルタ
    搬送波のエネルギーを落とし、次段のアクティブ回路の発熱を抑える。
  • 次段オペアンプによる2次～4次アクティブフィルタ
    この程度の傾きでも十分量子化ノイズを消せる。
• 外部アナログバッファアンプ

FPGA（FPGAボード）の選定については、SoC構成で、しかも20MHzサンプリングでの32ビット積和演算と200MHzオーダーのクロックを必要としますので、そこそこのスペックのものが必要になります。
もうしばらく検討を進めてからの選定となります。

使えそうな手持ちのFPGAボードはEfinix の Trion（速度的に不安）、Microsemi の PolarFire SoC（少しオーバースペック）、Altera の MAX10（ちょっと爆熱）といった感じで、これといった決定打が無いので悩ましいところ。かといって新しいボードを調達する予算は……。

ここからは、具体的なデジタル発振器のアルゴリズム検討に入ります。

ダイレクト型デジタル発振器（Direct Form Digital Oscillator）

z平面の単位円上に極（Pole）を配置した2次IIRフィルタ（伝達関数）を構成した、最も基本的なデジタル発振器です。

発振周波数を $f_0$、サンプリング周波数を $f_s$、サンプリング時間を $T = \frac{1}{f_s}$、1サンプルあたりの回転角を $\theta = \frac{2\pi f_0}{f_s}$ とします。

z平面の単位円周上に共役複素極 $z = e^{\pm j\theta} = \cos\theta \pm j\sin\theta$ を配置した2次の特性方程式は以下のようになります。

$$D(z) = (1 - e^{j\theta}z^{-1})(1 - e^{-j\theta}z^{-1}) = 1 - 2\cos\theta z^{-1} + z^{-2}$$

これを伝達関数 $H(z) = \frac{1}{D(z)}$ とし、この2次遅れ系を時間ドメインの状態空間表現に落とし込むと、状態方程式と出力方程式は以下のように書けます。

$ \begin{bmatrix} x_1[nT+1] \\ x_2[nT+1] \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ -1 & a \end{bmatrix} \begin{bmatrix} x_1[nT] \\ x_2[nT] \end{bmatrix} + \begin{bmatrix} 0 \\ 1\end{bmatrix} u[nT]$

$y[nT] = \begin{bmatrix} b & 0 \end{bmatrix} \begin{bmatrix} x_1[nT] \\ x_2[nT] \end{bmatrix} + \begin{bmatrix} 0 \end{bmatrix} u[nT]$

ここで、パラメータ $a$ $b$ は以下で与えられる定数。

$a = 2\cos\theta = 2\cos\left(\frac{2\pi f_0}{f_s}\right)$
$b = \frac{1}{16} \sin\theta$

つまり、周波数を変える時にのみパラメータを書き換えるだけで済みます。サンプリング周期ごとに三角関数を計算する必要はありません。
言い方を変えると、FPGAでのDSPブロック（IIRフィルタ）ではサンプリング周期毎、1回の積和演算と1回の加算（減算）処理をすれば良いことになります。

出力方程式のパラメータ $b$ は出力の振幅を一定に保つもので、上の式は周波数による変動を相殺する働きをします。実際の実装ではこの形にはせず、ソフト的に出力の振幅を安定化させるAGCを構成するので、計算ブロックで使用するのはパラメータ1個だけとなります。

なお、入力 $u[nT]$ は発振開始時のユニットステップ入力です。

発振開始時のステップ入力をパラメータ $b$ の値にして発振周波数に応じて調整しても良いです。この場合は出力にゲインを与える必要は無いです。

固定小数点で実装する際の課題と対策

このダイレクト型は構造が極めて単純な反面、固定小数点で実現する場合に2つの実用上の問題が発生します。

① 低周波域での感度問題（係数あふれ）

発振周波数 $f_0$ がサンプリング周波数 $f_s$ に対して非常に低い場合（$f_0 \ll f_s$）、$\theta \approx 0$ となり、$\cos\theta \approx 1$ です。
つまり、パラメータ $a = 2\cos\theta$ は 「2に近い値」 に張り付きます。
固定小数点表現において、$a$ のわずかな量子化誤差（丸め誤差）が、発振周波数の大きなズレを引き起こします。

② 振幅の拡大・縮小（丸め誤差）

理論上は $z^{-2}$ の係数が正確に -1（減衰なし）のため単位円上に極が留まりますが、固定小数点の乗算丸め処理（切り捨て等）によって極がわずかに単位円の外側（$r > 1$）に行くと振幅が発散（クリップ）し、内側（$r < 1$）に行くと減衰（消滅）します。

---

ダイレクト型はディジタル信号処理の教科書的な内容で、このまま実装できると思っていたのですが、AI（Gemini）に具体的な実装を相談していく過程で、この問題に言及してきました。そこで、彼のアドバイスに沿って実用的な構成を探っていきます。

Chamberlin / 結合型デジタル状態変数発振器

上記のダイレクト型の欠点（低周波での精度悪化）を克服し、固定小数点で安定して動かすために考案されたのがChamberlin（チャンバリン）型デジタル状態変数発振器です。

1サンプルごとの更新式は以下の通りです。

$$v_1[n] = v_1[n-1] + f \cdot v_2[n-1]$$

$$v_2[n] = v_2[n-1] - f \cdot v_1[n]$$

ここで、周波数パラメータ $f$ は近似的に $f \approx 2\sin\left(\frac{\pi f_0}{f_s}\right) \approx \frac{2\pi f_0}{f_s}$ です。

正直、こう説明されてもちょっと分かり難いのですが、要するに2つのフィルタを合成した（？）ようなイメージです。「状態変数型発振器のデジタル版」といった感じでしょうか？

周波数が低いとき、$f$ の値は 0 に近くなります。固定小数点フォーマットにおいて、0 に近い値は精度をフルに活用できるため、ダイレクト型のような「2の付近での感度問題」が発生しません。

$v_1$ と $v_2$ から、互いに90度位相がずれた正弦波と余弦波が同時に取り出せます（とのことです）。

固定小数点フォーマットの選定

今回、周波数上限をできるだけ上の方に取りたいという要請から、高速性を優先し、浮動小数点演算ではなく固定小数点で演算回路を構成することにしました。
16ビットでは分解能に不安があるので、ビット幅は32ビットとします。32ビットであれば、シミュレーションのかけやすさや FPGA ロジックの組みやすさにも支障は無いはずです。
そこで 「符号付き 32-bit int」 での設計とします。

Qフォーマットの選定については、アルゴリズムに応じて以下のフォーマットを選択するのが最も効率的らしい。

A. Q2.30（整数2ビット、小数30ビット）

• 範囲: $-2.0 \le x < 2.0$（1LSB $\approx 9.31 \times 10^{-10}$）
• 適用対象: ダイレクト型
• 選定理由:
  パラメータ $a = 2\cos\theta$ は最大で 2.0（$f_0 \to 0\text{Hz}$ の時）に達するため、1.0を超える数値を表現できる Q2.30 が必須となる。

B. Q1.31（整数1ビット、小数31ビット）

• 範囲: $-1.0 \le x < 1.0$（1LSB $\approx 4.65 \times 10^{-10}$）
• 適用対象: Chamberlin型 や 状態変数（正弦波・余弦波の出力値）
• 選定理由:
  正規化された振幅（$-1.0 \sim +1.0$）のダイナミックレンジを最大限に活用できる。

32ビット幅における2つのアルゴリズム比較

項目	ダイレクト型 ($2\cos\theta$)	Chamberlin（状態変数）型
毎サンプルの乗算回数	1回	2回
Qフォーマット	Q2.30	Q1.31（または Q2.30）
低周波分解能	Chamberlin型に劣るが、 32ビットあれば実用レベル	極めて高精度
周波数変更時の挙動	振幅が変動	振幅が維持される
2相出力 (Sine/Cosine)	取り出せない	取り出せる

Chamberlin型の方が優れているようですが、先ずはダイレクト型で実装してみて、Chamberlin型に置き換えた時の改善効果を見る方向で進めたいと思います。

これまで馴染みの無かったChamberlin型の動作理論がまだどうも頭の中で整理できていない、というのもあります。

64ビット中間精度と丸め処理の方針

32ビット固定小数点の乗算を行う際、内部積は64ビットになります。これを適切に32ビットに戻す処理が、波形の品質（歪み・持続性）を左右します。

乗算とビットシフトの手順（Q1.31の場合）

1. 32bit $\times$ 32bit の乗算を行い、64bit整数を得る
   Q1.31 同士の乗算結果は Q2.62 になります。
2. 31ビット右シフトして Q1.31 に戻す
   この時に単純な切り捨て（Floor）を行うと、毎サンプル微小な負のバイアスがかかり、振幅が徐々に衰退します。
3. 四捨五入（丸め：Rounding） を行う
   シフト前に $+2^{30}$ を加算する。

振幅の長期安定化（デジタルAGC / ノルム補正）

32ビット精度であっても、丸め誤差の集積によって数百万サンプル動かし続けると、振幅がわずかにドリフトします。これを防止するための軽量な補正方法です。

• 手法：数千サンプルに1回だけ行うスカラ補正
  毎サンプル重い計算をする必要はありません。例えば 4096サンプルに1回 だけ以下の処理を挿入します。

1. 現在の振幅自乗和 $R^2 = x_1^2 + x_2^2$ を計算する。
2. $R^2$ が目標値（例: $1.0$）から外れている場合、非常に緩やかな係数（例: $\alpha = 0.0001$）で状態変数を微修正する。

$x_1 \leftarrow x_1 \times (1 + \alpha (1 - R^2))$
$x_2 \leftarrow x_2 \times (1 + \alpha (1 - R^2))$

この程度であれば、波形歪みはほぼゼロのまま、何日間連続稼働させても完全に振幅が維持されます。
CPU向きの処理なので、FPGA演算器の状態変数を周期イベントでモニタし、算出した補正データ（重みデータ）をレジスタ経由で加算する回路を割り込ませる感じでしょうか。

ここから、具体的な検討に入っていきます。

シミュレーションの実行（32bit固定小数点エミュレーション）

実装へと入る前に、Pythonエミュレータで確認しておきます。
32ビット整数の乗算・ビットシフト・丸め処理をビット再現し、波形・直交性・FFTスペクトルを出力します。

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ設定 ---
fs = 48000.0          # サンプリング周波数 (Hz)
f0 = 1000.0           # 発振周波数 (Hz)
duration = 0.05       # シミュレーション時間 (秒)
N = int(fs * duration)
t = np.arange(N) / fs

# --- Q1.31 演算補助関数 ---
SCALE_31 = 2**31
MAX_INT32 = 2**31 - 1
MIN_INT32 = -2**31

def float_to_q31(val):
    return int(np.clip(np.round(val * SCALE_31), MIN_INT32, MAX_INT32))

def q31_to_float(val):
    return val / SCALE_31

def mult_q31(a, b):
    """32bit x 32bit -> 64bit, 四捨五入丸め付き 31bit右シフト"""
    prod = int(a) * int(b)  # 64bit積 (Q2.62)
    prod += (1 << 30)       # 四捨五入用オフセット
    res = prod >> 31        # Q1.31へ復元
    return int(np.clip(res, MIN_INT32, MAX_INT32))

# --- 初期化 ---
# 周波数係数 f = 2 * sin(pi * f0 / fs)
f_float = 2.0 * np.sin(np.pi * f0 / fs)
f_q31 = float_to_q31(f_float)

# 状態変数 (初期値: Cosine=0.999, Sine=0.0)
v1 = 0                  # Sine
v2 = float_to_q31(0.999) # Cosine

sine_out = np.zeros(N)
cosine_out = np.zeros(N)

# --- 発振ループ (Chamberlin差分方程式) ---
for n in range(N):
    # v1[n] = v1[n-1] + f * v2[n-1]
    v1 = v1 + mult_q31(f_q31, v2)
    v1 = int(np.clip(v1, MIN_INT32, MAX_INT32))
    
    # v2[n] = v2[n-1] - f * v1[n]
    v2 = v2 - mult_q31(f_q31, v1)
    v2 = int(np.clip(v2, MIN_INT32, MAX_INT32))
    
    sine_out[n] = q31_to_float(v1)
    cosine_out[n] = q31_to_float(v2)

# --- 描画 ---
plt.figure(figsize=(12, 6))

# 時間波形
plt.subplot(2, 1, 1)
plt.plot(t[:200] * 1000, sine_out[:200], label="Sine (v1)", color="blue")
plt.plot(t[:200] * 1000, cosine_out[:200], label="Cosine (v2)", color="orange", linestyle="--")
plt.title(f"Chamberlin Oscillator (32-bit Fixed Point) - {f0:.1f}Hz @ {fs:.0f}Hz")
plt.xlabel("Time [ms]")
plt.ylabel("Amplitude")
plt.grid(True)
plt.legend()

# FFTスペクトル
plt.subplot(2, 1, 2)
fft_vals = np.fft.rfft(sine_out)
fft_freqs = np.fft.rfftfreq(N, 1/fs)
fft_mag = 20 * np.log10(np.abs(fft_vals) / (N / 2) + 1e-12)

plt.plot(fft_freqs, fft_mag, color="blue")
plt.xlim(0, fs / 2)
plt.ylim(-140, 10)
plt.title("Frequency Spectrum (Sine Output)")
plt.xlabel("Frequency [Hz]")
plt.ylabel("Magnitude [dBFS]")
plt.grid(True)

plt.tight_layout()
plt.show()

サンプリング周波数などは今回の構成に合わせていない。

実行結果がこちら。

きれいに出ています。周波数スペクトルを見ても十分ノイズは少ない。これなら固定小数点のビット数をもう少し減らしても大丈夫そう。
サンプリング周波数を実機に合わせ込むなどして、これをベースに基本設計を詰めていきます。

回路検討（C言語コード例）

Gemini が提示してきたマイコン実装時のお手本コードです。
64ビット中間積を活用した高精度演算と、長期稼働時の振幅ドリフトを防止する軽量AGCを実装しています。
実際のソフトウェア（信号処理部分）はFPGAにオフロードするため、このままとはなりません。あくまで参考です。

#include <stdint.h>
#include <math.h>

#define PI 3.14159265358979323846

typedef struct {
    int32_t v1;      // Sine状態 (Q1.31)
    int32_t v2;      // Cosine状態 (Q1.31)
    int32_t f;       // 周波数係数 (Q1.31)
    uint32_t count;  // AGCトリガ用カウンタ
} ChamberlinOsc;

// Q1.31 乗算器 (四捨五入丸め付き)
static inline int32_t mult_q31(int32_t a, int32_t b) {
    int64_t temp = (int64_t)a * (int64_t)b;
    temp += (1LL << 30); // 丸め加算
    return (int32_t)(temp >> 31);
}

// 初期化
void osc_init(ChamberlinOsc *osc, float f0, float fs) {
    float f_val = 2.0f * sinf((float)PI * f0 / fs);
    osc->f = (int32_t)(f_val * 2147483648.0f);
    
    // 初期振幅設定 (約0.9)
    osc->v1 = 0;
    osc->v2 = 1932735283; // 0.9 in Q1.31
    osc->count = 0;
}

// 周波数動的変更
void osc_set_frequency(ChamberlinOsc *osc, float f0, float fs) {
    float f_val = 2.0f * sinf((float)PI * f0 / fs);
    osc->f = (int32_t)(f_val * 2147483648.0f);
}

// 1サンプル更新関数 (オーディオ割り込み内で呼び出し)
void osc_process(ChamberlinOsc *osc, int32_t *sine_out, int32_t *cosine_out) {
    // 1. Chamberlin差分方程式の計算
    // v1[n] = v1[n-1] + f * v2[n-1]
    osc->v1 += mult_q31(osc->f, osc->v2);
    
    // v2[n] = v2[n-1] - f * v1[n]
    osc->v2 -= mult_q31(osc->f, osc->v1);

    // 2. 4096サンプルに1回だけ実行する超軽量振幅補正 (AGC)
    if ((++osc->count & 0x0FFF) == 0) {
        // ノルム自乗 R^2 = v1^2 + v2^2 (Q2.62 -> Q1.31)
        int32_t v1_sq = mult_q31(osc->v1, osc->v1);
        int32_t v2_sq = mult_q31(osc->v2, osc->v2);
        int32_t r2 = v1_sq + v2_sq;
        
        // ターゲット振幅 0.9^2 = 0.81 (1739461754 in Q1.31)
        int32_t target_r2 = 1739461754;
        int32_t err = target_r2 - r2;
        
        // 微小ゲイン調整: alpha = 1/4096
        int32_t delta = err >> 12;
        osc->v1 += mult_q31(osc->v1, delta);
        osc->v2 += mult_q31(osc->v2, delta);
    }

    *sine_out = osc->v1;
    *cosine_out = osc->v2;
}

回路検討（VHDL記述例）

こちらもAI提示のRTLです。実機では、これを参考に実現可能な構成で再構築します。
サンプリングイネーブルパルス（sample_ce）に同期して毎サンプル2回の乗算と更新を行う VHDL モジュールです。

library ieee;
use ieee.std_logic_1164.all;
use ieee.numeric_std.all;

entity chamberlin_osc_q31 is
    generic (
        INIT_COS : integer := 1932735283 -- 振幅0.9 (Q1.31)
    );
    port (
        clk        : in  std_logic;
        rst_n      : in  std_logic;
        sample_ce  : in  std_logic;                      -- サンプリング周波数同期パルス
        coeff_f    : in  signed(31 downto 0);             -- 周波数係数 f (Q1.31)
        sine_out   : out signed(31 downto 0);             -- 正弦波出力 (v1)
        cosine_out : out signed(31 downto 0)              -- 余弦波出力 (v2)
    );
end entity;

architecture rtl of chamberlin_osc_q31 is

    signal v1_reg : signed(31 downto 0) := (others => '0');
    signal v2_reg : signed(31 downto 0) := to_signed(INIT_COS, 32);

begin

    process(clk, rst_n)
        variable mult1_full : signed(63 downto 0);
        variable mult2_full : signed(63 downto 0);
        variable v1_next    : signed(31 downto 0);
        variable v2_next    : signed(31 downto 0);
    begin
        if rst_n = '0' then
            v1_reg <= (others => '0');
            v2_reg <= to_signed(INIT_COS, 32);
        elsif rising_edge(clk) then
            if sample_ce = '1' then
                -- 1. mult1 = f * v2[n-1] (丸め加算含む)
                mult1_full := (coeff_f * v2_reg) + to_signed(1073741824, 64); -- + 2^30
                v1_next    := v1_reg + mult1_full(62 downto 31);

                -- 2. mult2 = f * v1[n] (丸め加算含む)
                mult2_full := (coeff_f * v1_next) + to_signed(1073741824, 64); -- + 2^30
                v2_next    := v2_reg - mult2_full(62 downto 31);

                -- レジスタ更新
                v1_reg <= v1_next;
                v2_reg <= v2_next;
            end if;
        end if;
    end process;

    sine_out   <= v1_reg;
    cosine_out <= v2_reg;

end architecture;

キーポイントとなるのは "32bit×32bit" の乗算です。高速化のためにはFPGAのDSPモジュールを活用したいところですが、これは一般的に19ビット回路ですので、各FPGAメーカーの論理合成ツールが上手いことプリミティブに割り当てて論理合成してくれるか、見極めながらの作業となります。

テストベンチ

FPGAのロジック検証用テストベンチ（テストシナリオ）の雛形も作っておいてしまいます。

PythonベースでRTL検証が行える cocotb を用いたテストベンチ環境です。

cocotbを使うと、HDLシミュレータ上で動作するVHDLモジュールに対してPythonからダイレクトにクロック供給、入力データの注入、出力波形の取得、およびMatplotlibによる波形画像（PNG）の保存まで一括処理できます。

1. テストベンチ-1 (test_chamberlin_osc.py)

サンプリング周波数 $f_s = 192\text{kHz}$、発振周波数 $f_0 = 1\text{kHz}$ を想定したテストベンチです。これも実運用時は、実際の設計に合わせて調整していきます。

• パラメータ値の計算:

$$f = 2 \cdot \sin\left(\frac{\pi \times 1000}{192000}\right) \approx 0.03272185$$

Q1.31 変換: $0.03272185 \times 2^{31} \approx 70,269,382$ （16進数: x"042FFDC6"）

1kHzから5kHzへのステップ変化および動的スイープを行い、結果をグラフ画像 cocotb_sweep_result.png に自動出力するします。

import cocotb
from cocotb.clock import Clock
from cocotb.triggers import RisingEdge, Timer
import matplotlib.pyplot as plt


def q31_to_float(val_int: int) -> float:
    """VHDLから取得した32bit符号付き整数(Q1.31)を float (-1.0 ~ +1.0) に変換"""
    if val_int >= (1 << 31):
        val_int -= 1 << 32
    return val_int / 2147483648.0


async def generate_strobe(dut, interval_ns=5208):
    """サンプリングストローブ (192kHz ≒ 5.208us) 自動発生タスク"""
    while True:
        await Timer(interval_ns, units="ns")
        await RisingEdge(dut.clk)
        dut.sample_strobe.value = 1
        await RisingEdge(dut.clk)
        dut.sample_strobe.value = 0


@cocotb.test()
async def test_chamberlin_sweep(dut):
    """Chamberlin型DSVOのステップ変化 & リアルタイムスイープ検証"""

    # 1. システムクロック起動 (100MHz = 10ns 周期)
    cocotb.start_soon(Clock(dut.clk, 10, units="ns").start())

    # 2. リセット & 初期入力設定
    dut.reset.value = 1
    dut.sample_strobe.value = 0
    dut.f_param.value = 0x042FFDC6  # 初期周波数 1kHz @ 192kHz fs
    await Timer(100, units="ns")

    dut.reset.value = 0
    await RisingEdge(dut.clk)

    # 3. 192kHz サンプリングストローブタスクの起動
    cocotb.start_soon(generate_strobe(dut, 5208))

    sine_history = []
    cosine_history = []

    # -------------------------------------------------------------
    # Phase 1: 1kHz で定常動作 (100 サンプル取得)
    # -------------------------------------------------------------
    dut._log.info("Phase 1: Running at 1kHz...")
    for _ in range(100):
        await RisingEdge(dut.sample_strobe)
        sine_history.append(q31_to_float(int(dut.sine_out.value)))
        cosine_history.append(q31_to_float(int(dut.cosine_out.value)))

    # -------------------------------------------------------------
    # Phase 2: 5kHz へステップ変更 (100 サンプル取得)
    # -------------------------------------------------------------
    dut._log.info("Phase 2: Step change to 5kHz...")
    dut.f_param.value = 0x14EB1B38  # 5kHz
    for _ in range(100):
        await RisingEdge(dut.sample_strobe)
        sine_history.append(q31_to_float(int(dut.sine_out.value)))
        cosine_history.append(q31_to_float(int(dut.cosine_out.value)))

    # -------------------------------------------------------------
    # Phase 3: 連続リアルタイムスイープ (500 サンプル取得)
    # -------------------------------------------------------------
    dut._log.info("Phase 3: Starting real-time sweep...")
    f_current = 0x14EB1B38
    for _ in range(500):
        f_current += 0x000A0000  # 毎サンプル f_param を増分
        dut.f_param.value = f_current
        await RisingEdge(dut.sample_strobe)
        sine_history.append(q31_to_float(int(dut.sine_out.value)))
        cosine_history.append(q31_to_float(int(dut.cosine_out.value)))

    # -------------------------------------------------------------
    # 波形の描画とグラフ画像保存
    # -------------------------------------------------------------
    dut._log.info("Generating sweep plot...")
    plt.figure(figsize=(12, 6))
    plt.plot(sine_history, label="Sine (v1)", color="blue")
    plt.plot(cosine_history, label="Cosine (v2)", color="orange", linestyle="--")

    # フェーズ境界線の描画
    plt.axvline(x=100, color="gray", linestyle=":", label="1kHz -> 5kHz Step")
    plt.axvline(x=200, color="red", linestyle=":", label="Sweep Start")

    plt.title("cocotb Verification: Chamberlin DSVO Frequency Sweep")
    plt.xlabel("Sample Count (@192kHz)")
    plt.ylabel("Normalized Amplitude (Q1.31)")
    plt.grid(True)
    plt.legend(loc="upper right")
    plt.savefig("cocotb_sweep_result.png")

    dut._log.info(
        "Test completed successfully! Output saved to 'cocotb_sweep_result.png'"
    )

2. 実行用 Makefile (Makefile)

シミュレータとして GHDL を使用する場合の一般的な Makefile です。（cocotb と ghdl がインストールされた環境で使用します）

# シミュレータの指定 (GHDL)
SIM ?= ghdl
TOPLEVEL_LANG ?= vhdl

# VHDLソースファイルの指定
VHDL_SOURCES += $(PWD)/chamberlin_osc.vhd

# トップレベル・モジュール名およびテストベンチファイル名
TOPLEVEL = chamberlin_osc
MODULE   = test_chamberlin_osc

# GHDL用コンパイルオプション (VHDL 2008 規格対応)
EXTRA_ARGS += --std=08

# cocotb のビルド・ルールをロード
include $(shell cocotb-config --makefiles)/Makefile.sim

3. 実行手順

1. 必要なライブラリの準備:

pip install cocotb matplotlib numpy

1. 同じディレクトリに3ファイルを配置:

• chamberlin_osc.vhd（本体）
• test_chamberlin_osc.py（Pythonテストベンチ）
• Makefile（実行ルール）

1. テストの実行:

make

実行が完了するとターミナルにテスト合格ログが表示され、同ディレクトリに cocotb_sweep_result.png が生成されます。

4. テストベンチ-2 test_chamberlin_osc.py

生成波形のTHD（歪み率）自動測定をやらせます。

FFT（高速フーリエ変換）解析による自動THD（Total Harmonic Distortion: 全高調波歪率）計算とアサーション機能を追加した test_chamberlin_osc.py です。

定常状態の正弦波データにハニング窓（Hann Window）を適用してスペクトルリークを抑えた上でFFTを行い、基本波電力に対する高調波成分（第2〜第10高調波）の割合を算出して自動判定します。

import cocotb
from cocotb.clock import Clock
from cocotb.triggers import RisingEdge, Timer
import numpy as np
import matplotlib.pyplot as plt


def q31_to_float(val_int: int) -> float:
    """VHDLから取得した32bit符号付き整数(Q1.31)を float (-1.0 ~ +1.0) に変換"""
    if val_int >= (1 << 31):
        val_int -= 1 << 32
    return val_int / 2147483648.0


def calculate_thd_db(signal_data: np.ndarray, fs: float, f0: float) -> float:
    """
    ハニング窓＋FFTによる全高調波歪率(THD)の計算 (dB表記)
    """
    N = len(signal_data)
    # 窓関数（ハニング窓）を適用してスペクトルリークを防止
    window = np.hanning(N)
    windowed_signal = signal_data * window

    # FFT実行と振幅スペクトルの取得
    fft_vals = np.abs(np.fft.rfft(windowed_signal))
    fft_freqs = np.fft.rfftfreq(N, 1.0 / fs)

    # 基本波（f0）のインデックスを取得
    fundamental_bin = np.argmin(np.abs(fft_freqs - f0))
    v1 = fft_vals[fundamental_bin]

    # 第2〜第10高調波の振幅を合算 (ピーク近傍のビン漏れを拾うため 3ビン幅で補正)
    harmonic_powers = 0.0
    for h in range(2, 11):
        target_freq = h * f0
        if target_freq >= fs / 2.0:  # ナイキスト周波数超えはスキップ
            break
        h_bin = np.argmin(np.abs(fft_freqs - target_freq))
        # ピークのエネルギーを加算
        h_val = np.max(fft_vals[max(0, h_bin - 1) : min(len(fft_vals), h_bin + 2)])
        harmonic_powers += h_val**2

    # THD計算 (比率 -> dB)
    thd_ratio = np.sqrt(harmonic_powers) / v1
    thd_db = 20.0 * np.log10(thd_ratio)
    return thd_db


async def generate_strobe(dut, interval_ns=5208):
    """サンプリングストローブ (192kHz ≒ 5.208us) 自動発生タスク"""
    while True:
        await Timer(interval_ns, units="ns")
        await RisingEdge(dut.clk)
        dut.sample_strobe.value = 1
        await RisingEdge(dut.clk)
        dut.sample_strobe.value = 0


@cocotb.test()
async def test_chamberlin_thd(dut):
    """32ビット Chamberlin型DSVO の THD 自動判定付きテスト"""

    fs = 192000.0  # サンプリング周波数
    f0 = 1000.0  # テスト周波数 (1kHz)
    thd_threshold_db = -90.0  # THD パス判断閾値 (32bit Q1.31目標値: -90dB以下)

    # 1. システムクロック起動 (100MHz)
    cocotb.start_soon(Clock(dut.clk, 10, units="ns").start())

    # 2. リセット & 初期入力設定
    dut.reset.value = 1
    dut.sample_strobe.value = 0
    # 1kHz @ 192kHz fs (Q1.31)
    dut.f_param.value = 0x042FFDC6
    await Timer(100, units="ns")

    dut.reset.value = 0
    await RisingEdge(dut.clk)

    # 3. サンプリングストローブの開始
    cocotb.start_soon(generate_strobe(dut, 5208))

    # 4. 過渡状態（ウォームアップ）の待機 (1000サンプルスキップ)
    dut._log.info("Warming up (Waiting for steady state)...")
    for _ in range(1000):
        await RisingEdge(dut.sample_strobe)

    # 5. FFT解析用のサンプル収集 (4096サンプル)
    num_samples = 4096
    dut._log.info(
        f"Collecting {num_samples} steady-state samples for THD analysis..."
    )

    sine_data = []
    for _ in range(num_samples):
        await RisingEdge(dut.sample_strobe)
        sine_data.append(q31_to_float(int(dut.sine_out.value)))

    sine_arr = np.array(sine_data)

    # 6. THD の計算
    thd_db = calculate_thd_db(sine_arr, fs, f0)
    dut._log.info(
        f"=== THD Test Result: {thd_db:.2f} dB (Target < {thd_threshold_db:.1f} dB) ==="
    )

    # 7. 波形とFFTスペクトルの保存 (エビデンス画像生成)
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 7))

    # 時間波形 (最初の500サンプル)
    ax1.plot(sine_arr[:500])
    ax1.set_title("Steady State Sine Wave (Q1.31)")
    ax1.set_xlabel("Samples")
    ax1.set_ylabel("Amplitude")
    ax1.grid(True)

    # スペクトル
    fft_vals = np.abs(np.fft.rfft(sine_arr * np.hanning(num_samples)))
    fft_freqs = np.fft.rfftfreq(num_samples, 1.0 / fs)
    mag_db = 20 * np.log10(fft_vals / np.max(fft_vals) + 1e-12)

    ax2.plot(fft_freqs, mag_db)
    ax2.axhline(
        y=thd_threshold_db, color="r", linestyle="--", label="THD Pass Threshold"
    )
    ax2.set_xlim(0, 20000)
    ax2.set_ylim(-140, 10)
    ax2.set_title(f"FFT Spectrum (Calculated THD: {thd_db:.2f} dB)")
    ax2.set_xlabel("Frequency [Hz]")
    ax2.set_ylabel("Magnitude [dBFS]")
    ax2.grid(True)
    ax2.legend()

    plt.tight_layout()
    plt.savefig("cocotb_thd_analysis.png")

    # 8. 自動判定 (アサーション)
    assert (
        thd_db < thd_threshold_db
    ), f"THD Test FAILED! THD was {thd_db:.2f} dB, exceeding threshold of {thd_threshold_db} dB"

    dut._log.info("THD Assertion PASSED! Signal quality is excellent.")

追加機能のポイント

1. ハニング窓の適用

• 有限サンプル長（4096点）の不連続性によるスペクトルリーク（裾野の広がり）を抑え、高調波レベルを正確に算出します。

1. 自動評価（アサーション）

• 測定されたTHD値が目標値（例: -90dB）を超えた場合、assert が発動してcocotbのテスト結果が FAIL と判定され、CI/CD環境での回帰テストに活用できます。

1. エビデンス画像 (cocotb_thd_analysis.png) の自動生成

• テスト実行時、スペクトルチャート上に「パス赤破線（閾値）」と実際のノイズフロアを描画したPNG画像を同時に自動書き出しします。

GTKWave で正弦波（アナログ）表示にする設定

今回は簡易的に Python のグラフ描画機能を使い直接表示する方法で構築しましたが、RTLのシミュレーションでは波形ビューワ（ GTKWave ）で信号波形を表示させるのが基本です。補足として軽く記載します。

GTKWaveは信号レベルもしくはディジタル値の表示がデフォルトですが、下記方法でアナログ表示させることができます。

1. 左上の SST ウィンドウから dut を選択し、下の Signal パネルから sine_out[31:0] と cosine_out[31:0] を Signals パネル（中央）へドラッグ＆ドロップ。
2. 追加した信号を右クリックし、以下を設定：
   • Data Format $\to$ Signed Decimal （符号付き10進表示に変更）
   • Wave Display Style $\to$ Analog Extension $\to$ Interpolated （アナログ補間表示に変更）
3. 波形の高さを広げると、sine_out（正弦波）と cosine_out（余弦波）が位相90度ズレて綺麗なサイン波を描いていることが確認できます。

シミュレーション時の注目・検証ポイント

波形ビューア（Analog表示設定）で以下のポイントを拡大して確認します。

1. フェーズ 2, 3 の「切り替え瞬間」の波形（ステップ変化時）

• 一般的なDDS（ダイレクト・デジタル・シンプセサイザー）やダイレクト型発振器では、パラメータを唐突に変更すると出力電圧が飛んだり（スパイクノイズ）、振幅が急激に増大/減衰したりします。
• Chamberlin型では、切り替えた瞬間の出力電圧レベルとその傾き（1階微分）が完全に保持されたまま、周期だけが自然に収縮する様子が確認できます。

1. フェーズ 4 の「連続スイープ」波形

• 時間経過とともにサイン波の山の間隔（周期）がスムーズに狭まっていきます。
• 周波数が5kHzから高くなる過程でも、正弦波のピーク値（包絡線）が全く上下に揺れず、一定のエンベロープ（Q1.31 の 1.0 相当）を維持していることが確認できます。

1ビットDACの実装

デジタル発振器に続き、生成した正弦波をアナログ信号として出力するDAC部の実装検討に移ります。

今回使用するのは、FPGAの端子から高速な1ビットのデルタ・シグマ変調（DSM / PDM）ビットストリームを直接吐き出し、外部のアナログLPFだけで正弦波を復元する手法です。

FPGAの高速ロジックとデジタル信号処理の自由度を最大限に活かしつつ、アナログ回路との境界をシンプルにまとめたアプローチです。

1. 全体アーキテクチャのイメージ

1. 発振器段:
   上記 Chamberlin型もしくはダイレクト型で高精度な 32bit 正弦波データ生成。
2. 変調段（DSM）:
   発振器の出力をより高クロックで動作する 高次デルタ・シグマ変調器 に通し、1ビットのパルス列（パルス密度変調: PDM）に変換します。これにより、量子化ノイズは可聴域から高周波帯域へと強力に追い出されます（ノイズシェーピング）。
3. 出力段:
   FPGAのIOピンから高スルーレートで出力。
4. フィルタ段:
   アナログのローパスフィルタで高周波ノイズをそぎ落とし、目的の正弦波を取り出します。

2. LVDS出力の活用

多彩な出力信号レベルに対応したFPGAならではのアプローチです。

• 一般的なアプローチ（CMOS出力）:
  単一のFPGAピン（例: 3.3V CMOS）からパルスを出力し、GND基準のRCフィルタに通します。
  • 弱点: FPGA内部のスイッチングノイズ（グラウンドバウンスや電源の揺れ）がそのまま出力波形の振幅変動やジッタになり、歪み（THD）を悪化させます。
• LVDS活用のアプローチ（差動1-bit出力）:
  FPGAの LVDSペア（真と偽の出力ピン） から、1ビットの差動ビットストリームを相補的（コンプリメンタリー）に出力します。
  • メリット:
    1. コモンモードノイズのキャンセル:
       差動信号であるため、電源やGNDに載る共通のノイズが相殺されます。
    2. 振幅の安定性:
       CMOSのVCC変動による振幅スケーリングエラーに強く、差動終端抵抗できれいにアイソレーションできます。
    3. 後段を差動入力のオペアンプ回路で受けることで、理想的なバランス（平衡）伝送・変換回路が構成でき、SN比が飛躍的に向上します。

3. 実装上の重要なポイント・注意点

① デルタ・シグマ変調器（DSM）の次数

• 1次 / 2次 DSM: 回路はシンプルですが、ノイズシェーピングの傾きが緩いため、可聴域付近にノイズが残りやすく、LPFの次数を高くしないとSN比が稼げません。
• 推奨: 最低でも 3次〜5次のノイズシェーピング（CRFBトポロジーなど） を用いると、可聴域のノイズフロアを -100dB 以下に押さえ込むことができます。MATLABの Delta Sigma Toolbox などの係数設計手法が参考になります。

② 外部アナログ LPF の設計

1ビットのデジタルパルス（0V / 3.3V または LVDSの振幅）を元のきれいなアナログに戻すため、LPFは必須です。

• 受動RCフィルタだけでは不十分: 1次や2次のRCでは高周波ノイズ（数MHzのキャリア成分）が漏れ出し、後段のアンプをクリップさせたり発熱させたりします。
• 推奨:
• 初段にパッシブの LC フィルタ（または簡単な RC）を置いてキャリアのエネルギーを落とし、
• 後段にオペアンプを用いた 多重帰還型（Multiple Feedback: MFB）アクティブフィルタ や Sallen-Key型 3次/4次フィルタ を配置します。

③ 電源のクリーンさ（ジッタと電源変動）

1ビットDACの最大の弱点は、「出力パルスのハイレベルの電圧（VCCIO）が揺れると、それがそのままアナログ波形の振幅誤差になる」 点です。

• FPGAの I/O バンク電源（VCCO）には、スイッチングレギュレータのノイズが乗らないよう、LDOレギュレータ を専用に用意するか、十分にデカップリングコンデンサ（セラミック＋OS-CONなど）を配置してインピーダンスを極限まで下げる配慮が必要です。

続く

長くなってしまったので、実検討～製作は別記事とします。ここまでは解説編ということで。