8ビットADCのFFT解析と性能評価

1. 概要

本記事では、8ビットADCのFFT解析を行い、SNDR（信号対ノイズ歪み比）、ENOB（有効ビット数）、SFDR（スプリアスフリー・ダイナミックレンジ）を計算します。Pythonを用いたシミュレーションを通じて、ADCの性能評価を実施します。

2. シミュレーション環境

サンプリング周波数: 1 MSPS（1MHz）
入力信号周波数: 124.511 kHz
振幅: 8ビットADCのフルスケール範囲
FFTポイント数: 4096
時間領域のシミュレーション: 4.096 ms
補足として、入力信号はコヒーレントサンプリングに基づいて計算します

3. 実装コード

以下のPythonコードでは、信号の生成、量子化、FFT解析を行い、SNDR・ENOB・SFDRの計算を実施します。

import numpy as np
import matplotlib.pyplot as plt

# -------------------- 関数定義 / Function Definitions --------------------

def freq_search(N: int, fs: float, finRate=8):
    """Calculate input frequency for coherent sampling"""
    freq = []
    FFT_points = 2**N
    for n in range(2, N):
        prime = np.exp(n * np.log(2)) - 1
        fin_per_fs = prime / FFT_points
        freq.append(fin_per_fs)
    array = np.array(freq)
    index = (np.abs(array - (1 / finRate))).argmin()
    fin = freq[index] * fs
    return fin

def wave_sort(data, fs: float, f_signal: float):
    """Sort data for one cycle"""
    data_num = len(data)
    t_signal = 1 / f_signal
    Ts = 1 / fs
    dx_list = [(Ts * n) % t_signal for n in range(data_num)]
    sorted_indices = sorted(range(len(dx_list)), key=lambda k: dx_list[k])
    sort = [data[i] for i in sorted_indices]
    return sort

def ideal_adc(input_voltage, n_bits: int, v_min: float, v_max: float):
    """Ideal ADC conversion"""
    LSB = (v_max - v_min) / (2**n_bits)
    digital_level = ((input_voltage - v_min) / LSB).astype(int)
    digital_level = np.clip(digital_level, 0, 2**n_bits - 1)
    return digital_level

def sine_wave_generator(N: int, f_signal: float, fs: float, amplitude: float, offset: float):
    """Generate sine wave"""
    t_sampling = 1 / fs
    duration = N * t_sampling
    t = np.arange(0, duration, t_sampling)
    signal = amplitude * np.sin(2 * np.pi * f_signal * t) + offset
    return signal

def AC_evaluation(bit: int, data, N: int, fs: float, f_in: float):
    """Evaluate ADC's AC performance"""
    f = np.divide(data, 2 ** (bit - 1))
    dt = 1 / fs
    t = np.arange(0, N * dt, dt)
    freq = np.fft.fftfreq(N, dt)
    freq_norm = freq / fs
    F = np.fft.fft(f) * 2 / N
    Power = np.abs(F) ** 2
    Pow_dB = 10 * np.log10(Power / 1.0)

    signal_index = np.argmax(Power[1 : int(N / 2)]) + 1
    Signal = Power[signal_index]
    Noise = np.sum(Power[1 : int(N / 2)]) - Signal
    SNDR = 10 * np.log10(Signal / Noise)
    ENOB = (SNDR - 1.76) / 6.02
    spurious = [np.max(Pow_dB[(signal_index * i) - 2 : signal_index * i + 2]) for i in range(2, 5)]
    SFDR = Pow_dB[signal_index] - np.max(spurious)

    # --- 結果表示 / Print results ---
    print("===== AC Evaluation Results =====")
    print("FFT Peak Frequency =", round(freq[signal_index], 2), "Hz")
    print("FFT Peak Power =", round(Pow_dB[signal_index], 2), "dB")
    print("SNDR =", round(SNDR, 2), "dB")
    print("ENOB =", round(ENOB, 2), "bit")
    print("SFDR =", round(SFDR, 2), "dB")
    print("=================================")

    # --- グラフ表示 / Plotting ---
    sorted_data = wave_sort(data, fs, f_in)
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))

    ax1.set_title("Sorted Wave")
    ax1.set_xlabel("Time [sec]")
    ax1.set_ylabel("Code")
    ax1.grid()
    ax1.set_ylim(0, 2**bit)
    ax1.plot(t, sorted_data, marker=".", linestyle="", markersize=4)

    ax2.set_title("Power Spectrum")
    ax2.set_xlabel("f/fs")
    ax2.set_ylabel("Power [dB]")
    ax2.grid()
    ax2.plot(freq_norm[1:int(N/2)], Pow_dB[1:int(N/2)])
    ax2.text(0.3, -5, f"SNDR = {round(SNDR, 2)} dB")
    ax2.text(0.3, -10, f"ENOB = {round(ENOB, 2)} bit")
    ax2.text(0.3, -15, f"SFDR = {round(SFDR, 2)} dB")

    plt.tight_layout()
    plt.show()

# -------------------- メイン処理 / Main Execution --------------------

N = 12                # FFTポイントの指数（2^Nポイント）
fs = 1000000          # サンプリング周波数
bit = 8               # ADCビット数
fin_rate = 8          # 入力周波数に対するfsの比（整合サンプリングのため）
v_amplitude = 2.5     # 振幅
v_offset = 2.5        # オフセット
v_min = 0             # ADCの最小電圧
v_max = 5             # ADCの最大電圧
noise = 0         # ノイズ振幅

# 入力周波数と信号生成
f_in = freq_search(N, fs, fin_rate)
print("Input signal frequency:", round(f_in, 2), "Hz")

# 入力信号の周期を計算 / Calculate signal period
T_signal = 1 / f_in
print("Input signal period T =", T_signal, "s")


# 入力信号の周期を計算 / Calculate signal period
T_signal = 1 / f_in

# 入力信号の周期を表示 / Print the input signal period
print("Input signal period T =", T_signal, "s")


# FFTポイントとシミュレーション時間
FFT_points = 2 ** N
Ts = 1 / fs
simulation_time = FFT_points * Ts

print("FFT Points =", FFT_points)
print("Sampling Period Ts =", Ts, "s")
print("Simulation Time =", simulation_time, "s ( = FFT Points × Ts )")

# パラメータ確認
print("===== Experiment Parameters =====")
print("N =", N, "->", 2**N, "points")
print("fs =", fs, "Hz")
print("bit =", bit)
print("fin_rate =", fin_rate)
print("v_amplitude =", v_amplitude, "V")
print("v_offset =", v_offset, "V")
print("v_min =", v_min, "V")
print("v_max =", v_max, "V")
print("noise amplitude =", noise, "V")
print("=================================")

# 入力周波数と信号生成
f_in = freq_search(N, fs, fin_rate)
print("Input signal frequency:", round(f_in, 2), "Hz")

v_signal = sine_wave_generator(2**N, f_in, fs, v_amplitude, v_offset)
v_signal += noise * np.random.rand(2**N)  # ノイズ加算



# ADC変換と評価
data = ideal_adc(v_signal, bit, v_min, v_max)
AC_evaluation(bit, data, 2**N, fs, f_in)

# 時間軸生成
t = np.arange(0, 2**N) / fs

# 入力信号全体プロット（アナログ信号）
plt.figure(figsize=(10, 4))
plt.plot(t, v_signal, label="Input Signal")
plt.title("Full Input Signal")
plt.xlabel("Time [s]")
plt.ylabel("Voltage [V]")
plt.grid()
plt.tight_layout()
plt.show()

# 入力信号全体プロットでADCコードをプロット
# ADCコードによる全体波形をプロット / Plot full ADC output (Code)
plt.figure(figsize=(10, 4))
plt.plot(t, data, label="ADC Output Code")  # ← data に変更
plt.title("Full ADC Output")
plt.xlabel("Time [s]")
plt.ylabel("Code")  # ← Voltage → Code に変更
plt.grid()
plt.tight_layout()
plt.show()


# 1周期のサンプル数でADCコードをプロット
cycle_length = int(fs / f_in)  # 1周期のサンプル数
data_cycle = data[:cycle_length]
t_cycle = t[:cycle_length]

plt.figure(figsize=(10, 4))
plt.plot(t_cycle, data_cycle, marker='o', linestyle='-', label="1 Cycle Code")
plt.title("One Cycle of ADC Output")
plt.xlabel("Time [s]")
plt.ylabel("Code")  # Codeとして縦軸表示
plt.grid()
plt.tight_layout()
plt.show()

# 5周期のサンプル数でADCコードをプロット

# -------------------- ADCコードの5周期分をプロット / Plot 5 cycles of ADC code --------------------

# 1周期あたりのサンプル数を計算 / Calculate samples per cycle
samples_per_cycle = int(fs / f_in)

# 5周期分のインデックス範囲 / Total number of samples for 5 cycles
num_samples_5cycles = samples_per_cycle * 5

# 時間軸とデータを抽出 / Extract time and ADC data for 5 cycles
t_5cycle = t[:num_samples_5cycles]
data_5cycle = data[:num_samples_5cycles]

# プロット / Plot
plt.figure(figsize=(10, 4))
plt.plot(t_5cycle, data_5cycle, marker='o', linestyle='-', label="5 Cycles of ADC Code")
plt.title("Five Cycles of ADC Output")
plt.xlabel("Time [s]")
plt.ylabel("Code")
plt.grid()
plt.tight_layout()
plt.show()

# -------------------- 入力信号とADCコードを先頭と末尾10点表示 / Display first and last 10 points of signal and ADC code --------------------

print("===== First 10 Points =====")
print(" Index |     Time [s]     |  Voltage [V]  |  ADC Code")
print("-----------------------------------------------")
for i in range(10):
    print(f"{i:6d} | {t[i]:.6f} s | {v_signal[i]:.6f} V | {data[i]:9d}")

print("\n===== Last 10 Points =====")
print(" Index |     Time [s]     |  Voltage [V]  |  ADC Code")
print("-----------------------------------------------")
for i in range(10):
    idx = -10 + i
    print(f"{len(t)+idx:6d} | {t[idx]:.6f} s | {v_signal[idx]:.6f} V | {data[idx]:9d}")




# --- 並び替え前後の波形を個別にプロット / Plot before and after sorting (separately) ---
sorted_data = wave_sort(data, fs, f_in)

# 並び替え前 / Original Wave
plt.figure(figsize=(10, 4))
plt.title("Original Wave")
plt.xlabel("Time [sec]")
plt.ylabel("Code")
plt.grid()
plt.ylim(0, 2**bit)
plt.plot(t, data, marker=".", linestyle="", markersize=4)
plt.tight_layout()
plt.show()

# 並び替え後 / Sorted Wave
plt.figure(figsize=(10, 4))
plt.title("Sorted Wave")
plt.xlabel("Time [sec]")
plt.ylabel("Code")
plt.grid()
plt.ylim(0, 2**bit)
plt.plot(t, sorted_data, marker=".", linestyle="", markersize=4)
plt.tight_layout()
plt.show()

結果

追加CSVファイルで読み込んだ波形をFFT

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from google.colab import files
from scipy.signal import find_peaks
from scipy.fft import fft

# ========================================
# 📥 ① CSVファイルのアップロード / Upload CSV
# ========================================
uploaded = files.upload()
filename = next(iter(uploaded))
df = pd.read_csv(filename)

# ========================================
# ⏱ ② 信号データの抽出と前処理 / Extract & preprocess signal
# ========================================
time = df["time_data"]
signal = df["signal_data"]
adjusted_signal = signal - 127  # DCオフセット除去 / Remove DC offset

# ========================================
# 📊 ③ 波形の表示 / Plot original & adjusted signals
# ========================================
plt.figure(figsize=(12, 4))
plt.plot(time, signal, label="Original Signal", color='blue', alpha=0.6)
plt.plot(time, adjusted_signal, label="Adjusted Signal", color='green', alpha=0.6)
plt.title("Original vs Adjusted Signal")
plt.xlabel("Time [s]")
plt.ylabel("Amplitude")
plt.grid(True)
plt.legend()
plt.show()

# ========================================
# 🧾 ④ 先頭と末尾の値を確認 / Print head & tail of signals
# ========================================
print("=== Original Signal: First 10 values ===")
print(signal.head(10))
print("=== Original Signal: Last 10 values ===")
print(signal.tail(10))

print("=== Adjusted Signal: First 10 values ===")
print(adjusted_signal.head(10))
print("=== Adjusted Signal: Last 10 values ===")
print(adjusted_signal.tail(10))

# ========================================
# 🎯 ⑤ サンプリング処理 / Sampling
# ========================================
sampling_freq = 10_000_000  # 10 MHz
sampling_period = 1 / sampling_freq
sampled_times = np.arange(time.iloc[0], time.iloc[-1], sampling_period)

# インデックスを検索してサンプリング / Find nearest indices
time_np = time.to_numpy()
idx = np.searchsorted(time_np, sampled_times, side="left")
idx = np.clip(idx, 0, len(time_np) - 1)

sampled_original = signal.to_numpy()[idx]
sampled_adjusted = adjusted_signal.to_numpy()[idx]

# プロット / Plot sampled signals
plt.figure(figsize=(12, 4))
plt.plot(time, signal, alpha=0.3, label="Original")
plt.scatter(sampled_times, sampled_original, s=10, c='blue', label="Sampled Original")
plt.plot(time, adjusted_signal, alpha=0.3, label="Adjusted")
plt.scatter(sampled_times, sampled_adjusted, s=10, c='green', label="Sampled Adjusted")
plt.title("Sampled Signals")
plt.xlabel("Time [s]")
plt.ylabel("Amplitude")
plt.legend()
plt.grid(True)
plt.show()

# ========================================
# ⏱ ⑥ 周波数推定 / Frequency estimation
# ========================================
def estimate_frequency(signal, label):
    peaks, _ = find_peaks(signal, height=10)
    if len(peaks) < 2:
        print(f"{label} → Not enough peaks.")
        return None
    peak_times = time.iloc[peaks].to_numpy()
    periods = np.diff(peak_times)
    avg_period = np.mean(periods)
    freq = 1 / avg_period
    print(f"{label} Frequency: {freq:.2f} Hz, Period: {avg_period:.2e} s")
    return avg_period

period_orig = estimate_frequency(signal, "Original")
period_adj = estimate_frequency(adjusted_signal, "Adjusted")

# ========================================
# 🔄 ⑦ 1周期で並び替え / Sort samples by phase
# ========================================
def sort_by_phase(t, y, period):
    phase = t % period
    sorted_idx = np.argsort(phase)
    return phase[sorted_idx], y[sorted_idx]

if period_orig and period_adj:
    phase_orig, sorted_orig = sort_by_phase(sampled_times, sampled_original, period_orig)
    phase_adj, sorted_adj = sort_by_phase(sampled_times, sampled_adjusted, period_adj)

    plt.figure(figsize=(12, 4))
    plt.plot(phase_orig, sorted_orig, '.', markersize=3, label="Original")
    plt.plot(phase_adj, sorted_adj, '.', markersize=3, label="Adjusted")
    plt.title("Sorted by Phase")
    plt.xlabel("Phase [s]")
    plt.ylabel("Amplitude")
    plt.legend()
    plt.grid(True)
    plt.show()

# ========================================
# 🔍 ⑧ FFTによるAC性能評価 / AC performance via FFT
# ========================================
def fft_ac_analysis(signal, fs, label):
    N = len(signal)
    bit_resolution = 8
    v_min = np.min(signal)
    v_max = np.max(signal)

    # 正規化（8ビットADC相当）/ Normalize to 8-bit range
    normalized = (signal - v_min) / (v_max - v_min) * (2**bit_resolution - 1)

    # FFT & パワースペクトル / FFT and power spectrum
    F = fft(normalized) * 2 / N
    Power = np.abs(F)**2
    freq = np.fft.fftfreq(N, d=1/fs)        # 実周波数 / Actual frequency
    freq_norm = freq / fs                   # 正規化周波数 / Normalized frequency
    Pow_dB = 10 * np.log10(Power + 1e-12)   # Avoid log(0)

    # 基本波成分 / Fundamental component
    signal_idx = np.argmax(Power[1:N//2]) + 1
    Signal = Power[signal_idx]
    Noise = np.sum(Power[1:N//2]) - Signal

    # SNDR, ENOB計算 / Compute SNDR and ENOB
    SNDR = 10 * np.log10(Signal / Noise)
    ENOB = (SNDR - 1.76) / 6.02

    # SFDR（スプリアス除去比）/ Spurious-Free Dynamic Range
    spurious = [np.max(Pow_dB[signal_idx * i - 2: signal_idx * i + 2]) for i in range(2, 5)]
    SFDR = Pow_dB[signal_idx] - np.max(spurious)

    # 正規化：最大成分を0dBに調整 / Normalize peak to 0 dB
    Pow_dB -= Pow_dB[signal_idx]

    # プロット / Plot
    plt.figure(figsize=(10, 4))
    plt.plot(freq_norm[1:N//2], Pow_dB[1:N//2])
    plt.title(f"{label} - Power Spectrum (Normalized to 0 dB)")
    plt.xlabel("f/fs")
    plt.ylabel("Power [dB]")
    plt.grid(True)
    plt.text(0.3, -5, f"SNDR = {SNDR:.2f} dB")
    plt.text(0.3, -10, f"ENOB = {ENOB:.2f} bit")
    plt.text(0.3, -15, f"SFDR = {SFDR:.2f} dB")
    plt.tight_layout()
    plt.show()

    # 結果の表示 / Print metrics
    print(f"===== {label} FFT Results =====")
    print(f"Peak Frequency Index: {signal_idx}")
    print(f"→ Normalized Freq: {freq_norm[signal_idx]:.6f} f/fs")
    print(f"→ Actual Freq:     {freq[signal_idx]:.2f} Hz")
    print(f"SNDR: {SNDR:.2f} dB")
    print(f"ENOB: {ENOB:.2f} bit")
    print(f"SFDR: {SFDR:.2f} dB")
    print("================================")

# FFT解析の実行 / Run FFT analysis
fft_ac_analysis(sampled_original, sampling_freq, "Original")
fft_ac_analysis(sampled_adjusted, sampling_freq, "Adjusted")


import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from google.colab import files
from scipy.signal import find_peaks
from scipy.fft import fft

# ========================================
# 📥 ① CSVファイルのアップロード / Upload CSV
# ========================================
# Google Colab で CSV ファイルをアップロードする
# Upload a CSV file from your local machine into Google Colab
uploaded = files.upload()
filename = next(iter(uploaded))
df = pd.read_csv(filename)

# ========================================
# ⏱ ② 信号データの抽出と前処理 / Extract & preprocess signal
# ========================================
# time_data列を取得し、全体を0.0001秒だけ前にシフト
# Get time_data and subtract 0.0001s from all values
time = df["time_data"] - 0.0001

# 信号データを取得
# Get the raw signal values
signal = df["signal_data"]

# DCオフセット（127）を除去して中央を0に揃える
# Remove DC offset assuming center at 127 (e.g., 8-bit ADC)
adjusted_signal = signal - 127

# 時間が0以上のデータだけを残す（負の時間は無効とする）
# Filter out negative time values (keep only time >= 0)
mask = time >= 0
time = time[mask].reset_index(drop=True)
signal = signal[mask].reset_index(drop=True)
adjusted_signal = adjusted_signal[mask].reset_index(drop=True)

# ========================================
# 📊 ③ 波形の表示 / Plot original & adjusted signals
# ========================================
# 元の信号とDCオフセット除去後の信号を同時にプロット
# Plot both original and adjusted signal on the same graph
plt.figure(figsize=(12, 4))
plt.plot(time, signal, label="Original Signal", color='blue', alpha=0.6)
plt.plot(time, adjusted_signal, label="Adjusted Signal", color='green', alpha=0.6)
plt.title("Original vs Adjusted Signal")
plt.xlabel("Time [s]")
plt.ylabel("Amplitude")
plt.grid(True)
plt.legend()
plt.show()

# ========================================
# 🧾 ④ 先頭と末尾の値を確認 / Print head & tail of signals
# ========================================
# 信号と調整後信号の最初と最後の10点を表示
# Print first and last 10 values of both original and adjusted signals
print("=== Original Signal: First 10 values ===")
print(signal.head(10))
print("=== Original Signal: Last 10 values ===")
print(signal.tail(10))

print("=== Adjusted Signal: First 10 values ===")
print(adjusted_signal.head(10))
print("=== Adjusted Signal: Last 10 values ===")
print(adjusted_signal.tail(10))

# ========================================
# 🎯 ⑤ サンプリング処理 / Sampling
# ========================================
# サンプリング周波数を10MHz（=0.1μs間隔）に設定
# Set sampling frequency to 10 MHz (0.1 µs interval)
sampling_freq = 10_000_000
sampling_period = 1 / sampling_freq

# サンプリング時間軸を作成（等間隔）
# Create evenly spaced sampling timestamps
sampled_times = np.arange(time.iloc[0], time.iloc[-1], sampling_period)

# 最も近い実データ点を検索し、インデックスとして取得
# Find indices of actual data points nearest to the sampling times
time_np = time.to_numpy()
idx = np.searchsorted(time_np, sampled_times, side="left")
idx = np.clip(idx, 0, len(time_np) - 1)  # 範囲外を防ぐ / Prevent out-of-bounds

# サンプリングした信号の配列を作成
# Extract sampled signals using the indices
sampled_original = signal.to_numpy()[idx]
sampled_adjusted = adjusted_signal.to_numpy()[idx]

# サンプリング結果のプロット
# Plot sampled points over the original and adjusted signals
plt.figure(figsize=(12, 4))
plt.plot(time, signal, alpha=0.3, label="Original")
plt.scatter(sampled_times, sampled_original, s=10, c='blue', label="Sampled Original")
plt.plot(time, adjusted_signal, alpha=0.3, label="Adjusted")
plt.scatter(sampled_times, sampled_adjusted, s=10, c='green', label="Sampled Adjusted")
plt.title("Sampled Signals")
plt.xlabel("Time [s]")
plt.ylabel("Amplitude")
plt.legend()
plt.grid(True)
plt.show()

# ========================================
# ⏱ ⑥ 周波数推定 / Frequency estimation
# ========================================
# ピーク検出から信号周期・周波数を推定する関数
# Function to estimate frequency from signal peaks
def estimate_frequency(signal, label):
    peaks, _ = find_peaks(signal, height=10)
    if len(peaks) < 2:
        print(f"{label} → Not enough peaks.")
        return None
    peak_times = time.iloc[peaks].to_numpy()
    periods = np.diff(peak_times)
    avg_period = np.mean(periods)
    freq = 1 / avg_period
    print(f"{label} Frequency: {freq:.2f} Hz, Period: {avg_period:.2e} s")
    return avg_period

# オリジナルと調整後の周期を取得
# Estimate the period for original and adjusted signals
period_orig = estimate_frequency(signal, "Original")
period_adj = estimate_frequency(adjusted_signal, "Adjusted")

# ========================================
# 🔄 ⑦ 1周期で並び替え / Sort samples by phase
# ========================================
# 位相順にデータを並べ替える関数
# Sort signal data by phase modulo period
def sort_by_phase(t, y, period):
    phase = t % period
    sorted_idx = np.argsort(phase)
    return phase[sorted_idx], y[sorted_idx]

# 周期が得られていれば、位相順で並べてプロット
# If period is valid, sort and plot by phase
if period_orig and period_adj:
    phase_orig, sorted_orig = sort_by_phase(sampled_times, sampled_original, period_orig)
    phase_adj, sorted_adj = sort_by_phase(sampled_times, sampled_adjusted, period_adj)

    plt.figure(figsize=(12, 4))
    plt.plot(phase_orig, sorted_orig, '.', markersize=3, label="Original")
    plt.plot(phase_adj, sorted_adj, '.', markersize=3, label="Adjusted")
    plt.title("Sorted by Phase")
    plt.xlabel("Phase [s]")
    plt.ylabel("Amplitude")
    plt.legend()
    plt.grid(True)
    plt.show()

# ========================================
# 🔍 ⑧ FFTによるAC性能評価 / AC performance via FFT
# ========================================
# FFTを使ってSNDR, ENOB, SFDRを計算しプロットする関数
# Analyze AC performance (SNDR, ENOB, SFDR) using FFT
def fft_ac_analysis(signal, fs, label):
    N = len(signal)
    bit_resolution = 8
    v_min = np.min(signal)
    v_max = np.max(signal)

    # 0〜255に正規化（8ビットADC相当）
    # Normalize signal to 8-bit digital levels
    normalized = (signal - v_min) / (v_max - v_min) * (2**bit_resolution - 1)

    # FFTの実行とパワースペクトルの計算
    # Perform FFT and compute power spectrum
    F = fft(normalized) * 2 / N
    Power = np.abs(F)**2
    freq = np.fft.fftfreq(N, d=1/fs)
    freq_norm = freq / fs
    Pow_dB = 10 * np.log10(Power + 1e-12)  # log(0)回避 / Avoid log(0)

    # 基本波のインデックスを探す
    # Identify fundamental frequency component
    signal_idx = np.argmax(Power[1:N//2]) + 1
    Signal = Power[signal_idx]
    Noise = np.sum(Power[1:N//2]) - Signal

    # SNDRとENOBの計算
    SNDR = 10 * np.log10(Signal / Noise)
    ENOB = (SNDR - 1.76) / 6.02

    # SFDR（スプリアスの最大値との差）
    spurious = [np.max(Pow_dB[signal_idx * i - 2: signal_idx * i + 2]) for i in range(2, 5)]
    SFDR = Pow_dB[signal_idx] - np.max(spurious)

    # 最大成分を0dBに正規化
    Pow_dB -= Pow_dB[signal_idx]

    # パワースペクトルのプロット
    plt.figure(figsize=(10, 4))
    plt.plot(freq_norm[1:N//2], Pow_dB[1:N//2])
    plt.title(f"{label} - Power Spectrum (Normalized to 0 dB)")
    plt.xlabel("f/fs")
    plt.ylabel("Power [dB]")
    plt.grid(True)
    plt.text(0.3, -5, f"SNDR = {SNDR:.2f} dB")
    plt.text(0.3, -10, f"ENOB = {ENOB:.2f} bit")
    plt.text(0.3, -15, f"SFDR = {SFDR:.2f} dB")
    plt.tight_layout()
    plt.show()

    # 結果出力
    print(f"===== {label} FFT Results =====")
    print(f"Peak Frequency Index: {signal_idx}")
    print(f"→ Normalized Freq: {freq_norm[signal_idx]:.6f} f/fs")
    print(f"→ Actual Freq:     {freq[signal_idx]:.2f} Hz")
    print(f"SNDR: {SNDR:.2f} dB")
    print(f"ENOB: {ENOB:.2f} bit")
    print(f"SFDR: {SFDR:.2f} dB")
    print("================================")

# FFTによる解析を実行 / Run FFT analysis
fft_ac_analysis(sampled_original, sampling_freq, "Original")
fft_ac_analysis(sampled_adjusted, sampling_freq, "Adjusted")