正弦波信号の並び替えと可視化
はじめに
本記事では、サンプリングされた正弦波信号の並び替え（wave sorting）を行い、
その結果を可視化するシミュレーションをPythonで実装します。
また、コヒーレントサンプリングに必要な入力周波数を求める関数 freq_search を定義し、
その計算結果を出力するコードについても解説します。
設定するパラメータ
以下のパラメータを使用してシミュレーションを行います。
- サンプリング周波数 (fs): 1 MHz (1,000,000 サンプル/秒)
-
入力信号周波数 (fin):
freq_search関数で算出 - サンプル数 (N): 1024
- シミュレーション時間: (N / fs)
コヒーレントサンプリングの周波数算出
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")
シミュレーションコード
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")