8ビットADCのFFT解析と性能評価のダミーシミュレーション(SNDR・ENOB・SFDRの計算)

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の計算を実施します。

# Program Name: ADC_Performance_Interactive_Architect
# Creation Date: 2026-01-22
# Purpose: Interactive simulation with PDF export capability for ADC AC Evaluation

import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
from IPython.display import display, clear_output
from google.colab import files
from datetime import datetime
import time

# =============================================================================
# 1. README（日本語・プレーンテキスト）
# =============================================================================
def print_documentation():
    print("""
-------------------------------------------------------------------------------
[シミュレータ概要 / Simulator Overview]
このシミュレータは、理想的なA/D変換器（ADC）のAC特性を評価するためのツールです。
コヒーレント・サンプリング理論に基づき、入力周波数を自動計算し、
SNDR（信号対雑音歪み比）、ENOB（有効ビット数）、SFDR（自由ダイナミックレンジ）を算出します。

[使い方 / How to use]
1. [Model Parameters] にADCの分解能や電圧範囲を入力します。
2. [Observed / Input Data] にサンプリング周波数やFFTポイント数を入力します。
3. [Environment / Options] でプロットの表示範囲を設定します。
4. "Run Simulation / 実行" をクリックして計算を開始します。
5. "Save PDF / PDF保存" をクリックすると、最新のグラフをダウンロードできます。
6. "Reset / リセット" で全ての状態を初期化します。

[支配方程式 / Governing Equations]
- コヒーレントサンプリング条件: fin = (M / N_fft) * fs (Mは素数)
- LSBサイズ: LSB = (Vmax - Vmin) / (2^bit)
- 理想SNR: SNR = 6.02 * bit + 1.76 [dB]
- ENOB: ENOB = (SNDR - 1.76) / 6.02
- 時間領域の並び替え（Sorted Wave）: サンプル時刻を信号周期で剰余計算(mod)してソート
-------------------------------------------------------------------------------
""")

# =============================================================================
# 2. Setup Code
# =============================================================================
# Seed 設定 / Set random seed
np.random.seed(42)
plt.style.use('seaborn-v0_8-muted')

# ログ保存用リスト / List for logging history
simulation_log = []
# 前回実行データ（ゴーストライン用） / Previous run data for ghost lines
ghost_data = None

# =============================================================================
# 3. TheoreticalModel
# =============================================================================
class ADCTheoreticalModel:
    @staticmethod
    def validate(params):
        """パラメータの妥当性検証 / Validate parameters"""
        if params['bit'] <= 0: raise ValueError("Bit depth must be positive. / ビット数は正である必要があります。")
        if params['fs'] <= 0: raise ValueError("Sampling frequency must be positive. / サンプリング周波数は正である必要があります。")
        if params['v_max'] <= params['v_min']: raise ValueError("Vmax must be greater than Vmin. / VmaxはVminより大きい必要があります。")
        if params['N_fft'] & (params['N_fft'] - 1) != 0: raise ValueError("FFT points must be power of 2. / FFTポイント数は2の累乗である必要があります。")

    @staticmethod
    def freq_search(N_fft, fs, fin_rate_target):
        """コヒーレントサンプリング用の周波数計算 / Calculate frequency for coherent sampling"""
        # N_fft / fin_rate_target に近い素数Mを探す
        M_target = N_fft / fin_rate_target
        # 簡易的に素数に近い奇数を選択
        M = int(M_target)
        if M % 2 == 0: M -= 1
        # MとN_fftが互いに素であることを確認（N_fftが2の累乗なのでMが奇数ならOK）
        fin = (M / N_fft) * fs
        return fin, M

    @staticmethod
    def print_calculation_steps(p, r):
        """計算過程の表示 / Print calculation steps"""
        print(f"\n[Calculation Steps / 計算過程]")
        print(f"1. LSB Size: (Vmax - Vmin) / 2^bit = ({p['v_max']} - {p['v_min']}) / 2^{p['bit']} = {r['lsb']:.6f} V")
        print(f"2. Input Freq (Coherent): (M / N_fft) * fs = ({r['M']} / {p['N_fft']}) * {p['fs']} = {r['f_in']:.2f} Hz")
        print(f"3. Ideal SNR: 6.02 * bit + 1.76 = 6.02 * {p['bit']} + 1.76 = {6.02*p['bit']+1.76:.2f} dB")
        print(f"4. Resulting SNDR: 10 * log10(Signal_Power / Noise_Power) = {r['sndr']:.2f} dB")
        print(f"5. Resulting ENOB: (SNDR - 1.76) / 6.02 = ({r['sndr']:.2f} - 1.76) / 6.02 = {r['enob']:.2f} bits")

    @staticmethod
    def run(p):
        """シミュレーション実行コア / Core simulation execution"""
        f_in, M = ADCTheoreticalModel.freq_search(p['N_fft'], p['fs'], p['fin_rate'])
        t = np.arange(p['N_fft']) / p['fs']
        
        # 信号生成 / Signal Generation
        v_in = p['v_amplitude'] * np.sin(2 * np.pi * f_in * t) + p['v_offset']
        v_in += p['noise'] * np.random.normal(0, 1, p['N_fft'])
        
        # 理想ADC量子化 / Ideal ADC Quantization
        lsb = (p['v_max'] - p['v_min']) / (2**p['bit'])
        code = np.floor((v_in - p['v_min']) / lsb).astype(int)
        code = np.clip(code, 0, 2**p['bit'] - 1)
        
        # 周波数解析 / Frequency Analysis
        norm_code = (code - np.mean(code)) / (2**(p['bit']-1))
        window = np.hanning(p['N_fft'])
        f_raw = np.fft.fft(norm_code * window)
        power = np.abs(f_raw[:p['N_fft']//2])**2
        # 窓関数の電力補正
        power = power / (np.sum(window)**2)
        power_db = 10 * np.log10(power + 1e-20)
        
        # 指標算出 / Metrics calculation
        bin_in = int(round(f_in / p['fs'] * p['N_fft']))
        sig_bins = 5 # 信号エネルギーと見なす幅
        p_sig = np.sum(power[max(0, bin_in-sig_bins):min(len(power), bin_in+sig_bins)])
        p_total = np.sum(power[1:])
        p_noise = p_total - p_sig
        
        sndr = 10 * np.log10(p_sig / p_noise)
        enob = (sndr - 1.76) / 6.02
        
        # SFDR (単純化のため最大スプリアスを探す)
        power_no_sig = np.copy(power)
        power_no_sig[max(0, bin_in-sig_bins):min(len(power), bin_in+sig_bins)] = 0
        p_spura = np.max(power_no_sig[1:])
        sfdr = 10 * np.log10(p_sig / p_spura)

        # 時間軸ソート (Sorted Wave) / Time sorting
        t_cycle = (t % (1/f_in)) * f_in
        sort_idx = np.argsort(t_cycle)
        
        results = {
            't': t, 'code': code, 't_sort': t_cycle[sort_idx], 'code_sort': code[sort_idx],
            'freqs': np.fft.fftfreq(p['N_fft'], 1/p['fs'])[:p['N_fft']//2], 'psd': power_db,
            'f_in': f_in, 'M': M, 'lsb': lsb, 'sndr': sndr, 'enob': enob, 'sfdr': sfdr
        }
        return results

# =============================================================================
# 4. Visualizer
# =============================================================================
class ADCVisualizer:
    @staticmethod
    def plot(params, current_res, ghost_res=None):
        fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
        
        # Plot 1: Sorted Waveform
        if ghost_res:
            ax1.plot(ghost_res['t_sort'], ghost_res['code_sort'], color='gray', alpha=0.3, label='Previous (Ghost)')
        ax1.plot(current_res['t_sort'], current_res['code_sort'], '.', markersize=2, label='Current Run')
        ax1.set_title('Sorted Time Domain Waveform (1 Cycle Equivalent)')
        ax1.set_xlabel('Phase [0 to 1]')
        ax1.set_ylabel('ADC Output Code [LSB]')
        ax1.set_xlim(params['p1_xmin'], params['p1_xmax'])
        ax1.set_ylim(params['p1_ymin'], params['p1_ymax'])
        ax1.grid(True)
        ax1.legend()
        
        # Plot 2: FFT Power Spectrum
        if ghost_res:
            ax2.plot(ghost_res['freqs']/1000, ghost_res['psd'], color='gray', alpha=0.3)
        ax2.plot(current_res['freqs']/1000, current_res['psd'], linewidth=1)
        ax2.set_title(f'Power Spectrum (SNDR: {current_res["sndr"]:.2f}dB)')
        ax2.set_xlabel('Frequency [kHz]')
        ax2.set_ylabel('Power [dBFS]')
        ax2.set_xlim(params['p2_xmin']/1000, params['p2_xmax']/1000)
        ax2.set_ylim(params['p2_ymin'], params['p2_ymax'])
        ax2.grid(True)
        
        plt.tight_layout()
        return fig

# =============================================================================
# 5. UIBuilder
# =============================================================================
# --- Widgets ---
w_bit = widgets.IntText(value=8, description='Resolution [bit]')
w_v_min = widgets.FloatText(value=0.0, description='Vmin [V]')
w_v_max = widgets.FloatText(value=5.0, description='Vmax [V]')
w_v_amp = widgets.FloatText(value=2.4, description='Amplitude [V]')
w_v_off = widgets.FloatText(value=2.5, description='Offset [V]')

w_fs = widgets.FloatText(value=1000000.0, description='fs [Hz]')
w_Nfft = widgets.IntText(value=4096, description='FFT Points')
w_fin_rate = widgets.FloatText(value=8.1, description='fs/fin Ratio')
w_noise = widgets.FloatText(value=0.01, description='Noise [Vrms]')

w_p1_xmin = widgets.FloatText(value=0.0, description='P1 X-min')
w_p1_xmax = widgets.FloatText(value=1.0, description='P1 X-max')
w_p1_ymin = widgets.FloatText(value=0, description='P1 Y-min')
w_p1_ymax = widgets.FloatText(value=256, description='P1 Y-max')
w_p2_xmin = widgets.FloatText(value=0, description='P2 X-min [Hz]')
w_p2_xmax = widgets.FloatText(value=500000, description='P2 X-max [Hz]')
w_p2_ymin = widgets.FloatText(value=-120, description='P2 Y-min [dB]')
w_p2_ymax = widgets.FloatText(value=10, description='P2 Y-max [dB]')

btn_run = widgets.Button(description="Run Simulation / 実行", button_style='primary')
btn_reset = widgets.Button(description="Reset / リセット", button_style='warning')
btn_save = widgets.Button(description="Save PDF / PDF保存", button_style='success')
btn_log = widgets.Button(description="Show Log / 履歴表示")
btn_clear = widgets.Button(description="Clear Log / 履歴消去")

output = widgets.Output()

# --- Layout ---
ui_model = widgets.VBox([widgets.Label('[Model Parameters]'), w_bit, w_v_min, w_v_max, w_v_amp, w_v_off])
ui_input = widgets.VBox([widgets.Label('[Observed / Input Data]'), w_fs, w_Nfft, w_fin_rate, w_noise])
ui_env = widgets.VBox([widgets.Label('[Environment / Options]'), 
                       widgets.HBox([w_p1_xmin, w_p1_xmax]), widgets.HBox([w_p1_ymin, w_p1_ymax]),
                       widgets.HBox([w_p2_xmin, w_p2_xmax]), widgets.HBox([w_p2_ymin, w_p2_ymax])])

ui_btns = widgets.HBox([btn_run, btn_reset, btn_save, btn_log, btn_clear])
ui_main = widgets.VBox([widgets.HBox([ui_model, ui_input, ui_env]), ui_btns, output])

# --- Logic ---
current_figure = None
last_results = None

def get_params():
    return {
        'bit': w_bit.value, 'v_min': w_v_min.value, 'v_max': w_v_max.value,
        'v_amplitude': w_v_amp.value, 'v_offset': w_v_off.value,
        'fs': w_fs.value, 'N_fft': w_Nfft.value, 'fin_rate': w_fin_rate.value, 'noise': w_noise.value,
        'p1_xmin': w_p1_xmin.value, 'p1_xmax': w_p1_xmax.value, 'p1_ymin': w_p1_ymin.value, 'p1_ymax': w_p1_ymax.value,
        'p2_xmin': w_p2_xmin.value, 'p2_xmax': w_p2_xmax.value, 'p2_ymin': w_p2_ymin.value, 'p2_ymax': w_p2_ymax.value
    }

def on_run_clicked(b):
    global ghost_data, last_results, current_figure
    with output:
        clear_output(wait=True)
        try:
            params = get_params()
            ADCTheoreticalModel.validate(params)
            
            print(f"[{datetime.now().isoformat()}] Starting Simulation... / シミュレーション開始...")
            
            results = ADCTheoreticalModel.run(params)
            ADCTheoreticalModel.print_calculation_steps(params, results)
            
            current_figure = ADCVisualizer.plot(params, results, ghost_data)
            plt.show()
            
            # Save to log
            log_entry = {
                'timestamp': datetime.now().isoformat(),
                'params': {'bit': params['bit'], 'noise': params['noise'], 'N_fft': params['N_fft']},
                'metrics': {'SNDR': results['sndr'], 'ENOB': results['enob'], 'SFDR': results['sfdr']}
            }
            simulation_log.append(log_entry)
            
            # Update ghost and last results
            ghost_data = results
            last_results = results
            
            print("Simulation Finished Successfully. / シミュレーション完了。")
            
        except Exception as e:
            print(f"Error / エラー: {str(e)}")

def on_reset_clicked(b):
    global ghost_data, simulation_log, current_figure
    ghost_data = None
    simulation_log = []
    current_figure = None
    w_bit.value = 8
    w_noise.value = 0.01
    w_p1_ymax.value = 256
    with output:
        clear_output()
        print_documentation()
        print("Reset Completed. / リセット完了。")

def on_save_clicked(b):
    global current_figure
    with output:
        if current_figure is None:
            print("No figure to save. Please run simulation first. / 保存する図がありません。先に実行してください。")
            return
        
        timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
        filename = f"adc_sim_result_{timestamp}.pdf"
        current_figure.savefig(filename, format='pdf')
        print(f"File '{filename}' generated. Starting download... / ファイル生成完了。ダウンロードを開始します...")
        files.download(filename)

def on_log_clicked(b):
    with output:
        if not simulation_log:
            print("Log is empty. / 履歴は空です。")
            return
        print("\n--- Simulation History / 実行履歴 ---")
        for entry in simulation_log:
            print(f"{entry['timestamp']} | Params: {entry['params']} | Metrics: {entry['metrics']}")

def on_clear_clicked(b):
    global simulation_log
    simulation_log = []
    with output:
        print("Log Cleared. / 履歴を消去しました。")

btn_run.on_click(on_run_clicked)
btn_reset.on_click(on_reset_clicked)
btn_save.on_click(on_save_clicked)
btn_log.on_click(on_log_clicked)
btn_clear.on_click(on_clear_clicked)

# =============================================================================
# 6. Initialization
# =============================================================================
with output:
    print_documentation()

display(ui_main)

結果

追加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")