コヒーレントサンプリングに必要な入力周波数を求める

# 正弦波信号の並び替えと可視化

## はじめに
本記事では、サンプリングされた正弦波信号の並び替え（wave sorting）を行い、
その結果を可視化するシミュレーションをPythonで実装します。

また、コヒーレントサンプリングに必要な入力周波数を求める関数 `freq_search` を定義し、
その計算結果を出力するコードについても解説します。

### 設定するパラメータ
以下のパラメータを使用してシミュレーションを行います。

- **サンプリング周波数 (fs):** 1 MHz (1,000,000 サンプル/秒)
- **入力信号周波数 (fin):** `freq_search` 関数で算出
- **サンプル数 (N):** 1024
- **シミュレーション時間:** \(N / fs\)

## コヒーレントサンプリングの周波数算出
```python
import numpy as np

def freq_search(N: int, fs: float, finRate=8):
    """
    コヒーレントサンプリングに必要な入力周波数を求める関数
    Function to calculate the input frequency for coherent sampling.

    Parameters
    ----------
    N : int
        サンプル数(N=log2(num)) | Log2 of sample count
    fs : float
        サンプリング周波数 | Sampling frequency
    finRate : int, optional
        入力周波数の位置 (デフォルト: 8) | Desired input freq position (default: 8)

    Returns
    -------
    fin : float
        算出された入力周波数 | Calculated input frequency
    """
    freq = []
    FFT_points = 2**N  # FFTポイント数 | Number of FFT points
    for n in range(2, N):
        prime = np.exp(n * np.log(2)) - 1  # 2^n-1を算出（オーバーフロー対策） | Compute 2^n - 1
        fin_per_fs = prime / FFT_points    # 正規化周波数 | Normalized frequency
        freq.append(fin_per_fs)

    array = np.array(freq)
    index = (np.abs(array - (1 / finRate))).argmin()  # 最も近いインデックスを取得 | Find nearest index

    fin = freq[index] * fs  # 入力周波数を計算 | Calculate input frequency
    return fin

# ==== パラメータ設定 | Parameter settings ====
N = 10          # 例: 2^10 = 1024 サンプル | Sample count
fs = 1e6        # サンプリング周波数 1 MHz | Sampling frequency
finRate = 8     # 入力周波数の位置 | Target fin/fs ratio

# ==== 関数実行 | Run function ====
calculated_fin = freq_search(N, fs, finRate)

# ==== 2^N を表示 | Print 2^N ====
FFT_points = 2**N
print(f"FFT points (2^{N}) = {FFT_points}")

# ==== 結果表示 | Show results ====
print(f"Calculated Input Frequency: {calculated_fin:.3f} Hz")

# ==== シミュレーション時間を計算 | Compute simulation time ====
sim_time = (FFT_points / fs) * 1.000000000000001  # 少し余裕を持たせる | Add slight margin

# ==== シミュレーション時間を表示 | Show simulation time ====
print(f"Recommended simulation time: {sim_time:.6f} seconds")


```

## シミュレーションコード
```python
import matplotlib.pyplot as plt

def wave_sort(data, fs, fin):
    """
    信号の並び替えを行う関数
    :param data: 入力信号
    :param fs: サンプリング周波数
    :param fin: 入力信号周波数
    :return: 並び替え後の信号
    """
    indices = np.argsort(np.angle(np.exp(1j * 2 * np.pi * fin * np.arange(len(data)) / fs)))
    return data[indices]

# シミュレーションパラメータ
simulation_time = N / fs  # シミュレーション時間

# 時間軸
t = np.arange(0, simulation_time, 1 / fs)
original_data = np.sin(2 * np.pi * calculated_fin * t)  # 正弦波信号

# 並び替え
sorted_data = wave_sort(original_data, fs, calculated_fin)

# プロット
plt.figure(figsize=(10, 6))

# 並び替え前の波形
plt.subplot(2, 1, 1)
plt.plot(original_data, marker='o', linestyle='-', markersize=3, label="Original Data")
plt.title("Original Data")
plt.xlabel("Sample Index")
plt.ylabel("Amplitude")
plt.legend()
plt.grid()

# 並び替え後の波形
plt.subplot(2, 1, 2)
plt.plot(sorted_data, marker='o', linestyle='-', markersize=3, label="Sorted Data", color='r')
plt.title("Sorted Data")
plt.xlabel("Sample Index")
plt.ylabel("Amplitude")
plt.legend()
plt.grid()

plt.tight_layout()
plt.show()
```

## 実行結果と考察

### コヒーレントサンプリングの周波数選定
- `freq_search` 関数を用いて、サンプリングと整合する入力周波数を計算。
- 計算結果を `print()` により確認可能。

### 並び替え前の波形
もとの波形は通常の正弦波信号であり、サンプリング間隔ごとにデータが取得されています。

### 並び替え後の波形
並び替え後の波形は、入力周波数に基づいた角度の昇順でソートされています。
これにより、波形が特定のパターンで整列されることが確認できます。

## まとめ
本記事では、
- コヒーレントサンプリングに必要な入力周波数の算出
- サンプリング信号の並び替え
- ソートされた波形の可視化

を実装しました。この技術は、信号処理やデータ解析の一環として利用可能です。




おまけ　活用事例
```
import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft

# 条件設定
f_in = 124.511e3  # 入力サイン周波数 (Hz)
f_s = 1e6  # サンプリング周波数 (Hz)
N = 2048  # FFT点数
A = 1  # スケール係数

# 計算: サンプリング周期とシミュレーション時間
T_s = 1 / f_s  # サンプリング周期
T_in = 1 / f_in  # 入力サイン周期
T_sim = T_s * N * A  # シミュレーション時間

# 時間軸の生成
t = np.arange(0, T_sim, T_s)  # 0からシミュレーション時間までの時間配列

# サイン波の生成
signal = np.sin(2 * np.pi * f_in * t)

# FFT解析
spectrum = fft(signal, N)  # FFT変換
spectrum_magnitude = np.abs(spectrum[:N//2])  # FFTの前半分（正の周波数成分）
freqs = np.fft.fftfreq(N, T_s)[:N//2]  # 周波数軸

# SNDR（信号対ノイズ＋歪比）計算
signal_power = np.max(spectrum_magnitude) ** 2  # 信号のピークパワー
noise_power = np.sum(spectrum_magnitude ** 2) - signal_power  # ノイズ＋歪成分のパワー
SNDR = 10 * np.log10(signal_power / noise_power)  # SNDR計算

# ENOB（有効ビット数）の計算
ENOB = (SNDR - 1.76) / 6.02  # ENOB計算

# SFDR（スプリアスフリーダイナミックレンジ）の計算
spurious_spectrum = spectrum_magnitude.copy()
spurious_spectrum[np.argmax(spectrum_magnitude)] = 0  # 最大信号成分を除外
SFDR = 20 * np.log10(np.max(spectrum_magnitude) / np.max(spurious_spectrum))  # SFDR計算

# FFTスペクトルの正規化
spectrum_magnitude_db = 20 * np.log10(spectrum_magnitude / np.max(spectrum_magnitude))

# ------------------------------------
# **プロット**
# (1) 入力サイン波
# (2) FFTスペクトル
# ------------------------------------

plt.figure(figsize=(10, 6))

# ---- (1) 入力サイン波のプロット ----
plt.subplot(2, 1, 1)
plt.plot(t[:200], signal[:200], 'b-', label='Input Sine Wave')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title(f'Sampled Sine Wave (f_in = {f_in / 1e3:.3f} kHz, f_s = {f_s / 1e6:.1f} MHz)')
plt.grid(True)
plt.legend()

# ---- (2) FFTスペクトルのプロット ----
plt.subplot(2, 1, 2)
plt.plot(freqs / f_s, spectrum_magnitude_db, 'r-', label='FFT Spectrum (Normalized)')
plt.axhline(y=0, color='k', linestyle='--', label='0 dB Reference')
plt.xlabel('Normalized Frequency (f/f_s)')
plt.ylabel('Magnitude [dB]')
plt.title('FFT Spectrum (0 dB = max amplitude)')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

# ------------------------------------
# **計算結果の表示**
# ------------------------------------
print(f"入力サイン周期 (T_in): {T_in:.9f} s")
print(f"サンプリング周期 (T_s): {T_s:.9f} s")
print(f"シミュレーション時間 (T_sim): {T_sim:.9f} s")
print(f"SNDR: {SNDR:.2f} dB")
print(f"ENOB: {ENOB:.2f} bits")
print(f"SFDR: {SFDR:.2f} dB")
```