0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@digital_amplifier_1

FPGA風 DDS / NCO サイン波生成シミュレータ

0
Posted at 
import numpy as np
import matplotlib.pyplot as plt

# ============================================================
# FPGA風 DDS / NCO サイン波生成シミュレータ
# - 位相加算器
# - FCW
# - 位相打切り
# - 1/4波LUT
# - 振幅量子化
# - DAC量子化
# - パイプライン遅延
# - ZOH
# - FFT
# ============================================================

# -----------------------------
# 設定値
# -----------------------------
FS_CLK = 100e6
F_OUT_REQ = 1.234e6
PHASE_BITS = 32
ADDR_BITS = 10
LUT_AMP_BITS = 14
DAC_BITS = 12
N_SAMPLES = 65536
PIPELINE_DELAY = 2
SHOW_TIME_SAMPLES = 400
ZOH_OVERSAMPLE = 16

ENABLE_PHASE_DITHER = False
DITHER_BITS = 4


# -----------------------------
# 1/4波LUTを作る
# -----------------------------
def build_quarter_wave_lut(addr_bits: int, amp_bits: int) -> np.ndarray:
    lut_size = 2 ** addr_bits
    amp_max = (2 ** (amp_bits - 1)) - 1

    theta = np.linspace(0, np.pi / 2, lut_size, endpoint=False)
    lut = np.round(amp_max * np.sin(theta)).astype(np.int32)
    return lut


# -----------------------------
# 位相から1/4波LUT参照値を取り出す
# -----------------------------
def quarter_wave_lookup(
    phase_word: int,
    phase_bits: int,
    addr_bits: int,
    lut: np.ndarray
) -> int:
    total_lookup_bits = 2 + addr_bits
    shift = phase_bits - total_lookup_bits
    lookup_word = (phase_word >> shift) & ((1 << total_lookup_bits) - 1)

    quadrant = (lookup_word >> addr_bits) & 0b11
    idx = lookup_word & ((1 << addr_bits) - 1)

    last = (1 << addr_bits) - 1

    if quadrant == 0:
        val = lut[idx]
    elif quadrant == 1:
        val = lut[last - idx]
    elif quadrant == 2:
        val = -lut[idx]
    else:
        val = -lut[last - idx]

    return int(val)


# -----------------------------
# signed量子化
# -----------------------------
def quantize_signed(x: np.ndarray, bits: int) -> np.ndarray:
    maxv = (2 ** (bits - 1)) - 1
    minv = -(2 ** (bits - 1))
    y = np.round(x).astype(np.int64)
    y = np.clip(y, minv, maxv)
    return y.astype(np.int32)


# -----------------------------
# DDS本体
# -----------------------------
def simulate_dds(
    fs_clk: float,
    f_out_req: float,
    phase_bits: int,
    addr_bits: int,
    lut_amp_bits: int,
    dac_bits: int,
    n_samples: int,
    pipeline_delay: int = 0,
    enable_phase_dither: bool = False,
    dither_bits: int = 0
):
    phase_mod = 1 << phase_bits
    fcw = int(round((f_out_req / fs_clk) * phase_mod))
    f_out_actual = fcw * fs_clk / phase_mod

    lut = build_quarter_wave_lut(addr_bits, lut_amp_bits)

    phase_acc = 0
    phase_hist = np.zeros(n_samples, dtype=np.uint64)
    lut_out_hist = np.zeros(n_samples, dtype=np.int32)
    dac_out_hist = np.zeros(n_samples, dtype=np.int32)

    pipe = [0 for _ in range(max(0, pipeline_delay))]
    rng = np.random.default_rng(0)

    for n in range(n_samples):
        # 位相加算器
        phase_acc = (phase_acc + fcw) & (phase_mod - 1)

        # 微小ディザ
        phase_for_lookup = phase_acc
        if enable_phase_dither and dither_bits > 0:
            dither = rng.integers(0, 1 << dither_bits)
            phase_for_lookup = (phase_for_lookup + dither) & (phase_mod - 1)

        # LUT参照
        lut_val = quarter_wave_lookup(
            phase_word=phase_for_lookup,
            phase_bits=phase_bits,
            addr_bits=addr_bits,
            lut=lut
        )

        # LUT出力をDACビット数へ丸める
        shift_down = max(lut_amp_bits - dac_bits, 0)
        dac_val = int(np.round(lut_val / (2 ** shift_down)))

        # パイプライン遅延
        if pipeline_delay > 0:
            pipe.append(dac_val)
            dac_val_delayed = pipe.pop(0)
        else:
            dac_val_delayed = dac_val

        phase_hist[n] = phase_acc
        lut_out_hist[n] = lut_val
        dac_out_hist[n] = dac_val_delayed

    t = np.arange(n_samples) / fs_clk

    return {
        "t": t,
        "fcw": fcw,
        "f_out_actual": f_out_actual,
        "phase_hist": phase_hist,
        "lut_out": lut_out_hist,
        "dac_out": dac_out_hist,
        "lut": lut,
    }


# -----------------------------
# ZOH波形を表示用に作る
# -----------------------------
def make_zoh_wave(t: np.ndarray, x: np.ndarray, oversample: int = 8):
    dt = t[1] - t[0]
    tz = np.repeat(t, oversample) + np.tile(
        np.arange(oversample) * dt / oversample,
        len(t)
    )
    xz = np.repeat(x, oversample)
    return tz, xz


# -----------------------------
# FFT
# -----------------------------
def calc_spectrum_dbfs(x: np.ndarray, fs: float):
    x = x.astype(np.float64)
    x = x - np.mean(x)

    window = np.hanning(len(x))
    xw = x * window
    X = np.fft.rfft(xw)
    f = np.fft.rfftfreq(len(x), d=1 / fs)

    mag = np.abs(X)
    mag = mag / (np.max(mag) + 1e-30)
    mag_db = 20 * np.log10(mag + 1e-15)
    return f, mag_db


# -----------------------------
# 実行
# -----------------------------
res = simulate_dds(
    fs_clk=FS_CLK,
    f_out_req=F_OUT_REQ,
    phase_bits=PHASE_BITS,
    addr_bits=ADDR_BITS,
    lut_amp_bits=LUT_AMP_BITS,
    dac_bits=DAC_BITS,
    n_samples=N_SAMPLES,
    pipeline_delay=PIPELINE_DELAY,
    enable_phase_dither=ENABLE_PHASE_DITHER,
    dither_bits=DITHER_BITS,
)

t = res["t"]
dac_out = res["dac_out"]
lut_out = res["lut_out"]
phase_hist = res["phase_hist"]

dac_norm = dac_out / np.max(np.abs(dac_out))
lut_norm = lut_out / np.max(np.abs(lut_out))

tz, xz = make_zoh_wave(
    t[:SHOW_TIME_SAMPLES],
    dac_norm[:SHOW_TIME_SAMPLES],
    oversample=ZOH_OVERSAMPLE
)

freq, mag_db = calc_spectrum_dbfs(dac_out, FS_CLK)

# 時間波形
plt.figure(figsize=(10, 4))
plt.plot(
    t[:SHOW_TIME_SAMPLES] * 1e6,
    lut_norm[:SHOW_TIME_SAMPLES],
    label="LUT output (normalized)",
    linewidth=1.0
)
plt.step(
    tz * 1e6,
    xz,
    where="post",
    label="DAC output with ZOH",
    linewidth=1.0
)
plt.xlabel("Time [us]")
plt.ylabel("Amplitude")
plt.title("FPGA-like DDS Sine Wave (time domain)")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

# スペクトル
plt.figure(figsize=(10, 4))
plt.plot(freq / 1e6, mag_db, linewidth=1.0)
plt.xlim(0, 10)
plt.ylim(-140, 5)
plt.xlabel("Frequency [MHz]")
plt.ylabel("Magnitude [dBFS]")
plt.title("DDS Spectrum")
plt.grid(True)
plt.tight_layout()
plt.show()

# 位相加算器
plt.figure(figsize=(10, 3.5))
plt.plot(np.arange(120), phase_hist[:120], linewidth=1.0)
plt.xlabel("Sample index")
plt.ylabel("Phase accumulator")
plt.title("Phase accumulator behavior")
plt.grid(True)
plt.tight_layout()
plt.show()

print("========== DDS Parameters ==========")
print(f"Requested output frequency  : {F_OUT_REQ:.6f} Hz")
print(f"Actual output frequency     : {res['f_out_actual']:.6f} Hz")
print(f"Clock frequency             : {FS_CLK:.6f} Hz")
print(f"FCW                         : {res['fcw']}")
print(f"Phase accumulator bits      : {PHASE_BITS}")
print(f"LUT address bits            : {ADDR_BITS}")
print(f"LUT amplitude bits          : {LUT_AMP_BITS}")
print(f"DAC bits                    : {DAC_BITS}")
print(f"Pipeline delay              : {PIPELINE_DELAY} clocks")
print(f"Samples                     : {N_SAMPLES}")
print("====================================")
0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?