Note
このページは Jupyter Notebook から生成されました。 ノートブックをダウンロード (.ipynb)
チュートリアル: 伝達関数の測定
このチュートリアルでは、gwexpy を使用して測定データからシステムの伝達関数（Transfer Function）を推定し、 物理モデルへのフィッティングを行ってパラメータを同定する流れを解説します。
シナリオ: ある機械振動系（または電気回路）に白色雑音を入力し、その応答を測定しました。 この入出力データから、システムの共振周波数 \(f_0\) と Q値 \(Q\) を求めます。
ワークフロー:
データ生成: 入力信号（白色雑音）と、それに対するシステムの出力（共振系応答 + 測定ノイズ）をシミュレーションします。
伝達関数の推定: 入出力の時系列データから伝達関数とコヒーレンスを計算し、ボード線図（Bode Plot）を描画します。
モデルフィッティング: 推定した伝達関数に理論モデル（ローレンツ関数など）を当てはめ、パラメータを推定します。
[1]:
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
from gwexpy import TimeSeries
from gwexpy.plot import Plot
1. 実験データの生成（シミュレーション）
仮想的な実験を行います。
システム: 2次共振系（単振動）
共振周波数 \(f_0 = 300\) Hz
Q値 \(Q = 50\)
入力: 白色雑音（White Noise）
測定: サンプリング周波数 2048 Hz, 継続時間 60秒
出力信号には、測定に伴うノイズも付加します。
[2]:
# --- Set ---
fs = 2048.0 # Sample Rate [Hz]
duration = 60.0 # time / duration [s]
f0_true = 300.0 # ResonanceFrequency [Hz]
Q_true = 50.0 # Q
# --- time / durationAxis Signal ---
t = np.linspace(0, duration, int(duration * fs), endpoint=False)
input_data = np.random.normal(0, 1, size=len(t)) # Color
# --- (scipy.signal ) ---
# Transfer Function H(s) = w0^2 / (s^2 + (w0/Q)s + w0^2)
w0 = 2 * np.pi * f0_true
num = [w0**2]
den = [1, w0 / Q_true, w0**2]
system = signal.TransferFunction(num, den)
# time / durationResponse Calculate (lsim)
# Time steps for simulation
_, output_clean, _ = signal.lsim(system, U=input_data, T=t)
# --- Noise ---
# Noise（ Noise ）
measurement_noise = np.random.normal(0, 0.1, size=len(t))
output_data = output_clean + measurement_noise
# --- gwexpy TimeSeries Create ---
ts_input = TimeSeries(input_data, t0=0, sample_rate=fs, name="Input", unit="V")
ts_output = TimeSeries(output_data, t0=0, sample_rate=fs, name="Output", unit="V")
print("Input data shape:", ts_input.shape)
print("Output data shape:", ts_output.shape)
Plot(ts_input, ts_output, separate=True);
Input data shape: (122880,)
Output data shape: (122880,)
2. 伝達関数の推定
入出力の時系列データ ts_input と ts_output から、周波数応答関数（伝達関数）を推定します。 gwexpy の transfer_function メソッドを使用します。
また、測定の信頼性を確認するためにコヒーレンス（Coherence）も計算します。 コヒーレンスが 1 に近い周波数帯域は、入出力の関係が線形であり、ノイズの影響が少ないことを示します。
[3]:
# FFTSet: Mean Data Length
fftlength = 2.0 # 2seconds FFT Mean
# Transfer Function Calculate (Output / Input)
tf = ts_input.transfer_function(ts_output, fftlength=fftlength)
# Coherence Calculate
coh = ts_input.coherence(ts_output, fftlength=fftlength) ** 0.5
# --- Plot ( & Coherence) ---
fig, axes = plt.subplots(3, 1, figsize=(10, 10), sharex=True)
# Amplitude (Magnitude)
ax = axes[0]
ax.loglog(tf.abs(), label="Measured TF")
ax.set_ylabel("Gain [V/V]")
ax.set_title("Bode Plot: Magnitude")
ax.grid(True, which="both", alpha=0.5)
# Phase (Phase)
ax = axes[1]
ax.plot(tf.degree(), label="Measured Phase")
ax.set_ylabel("Phase [deg]")
ax.set_title("Bode Plot: Phase")
ax.set_yticks(np.arange(-180, 181, 90))
ax.grid(True, which="both", alpha=0.5)
# Coherence (Coherence)
ax = axes[2]
ax.plot(coh, color="green", label="Coherence")
ax.set_ylabel("Coherence")
ax.set_xlabel("Frequency [Hz]")
ax.set_ylim(0, 1.1)
ax.set_xlim(10, 1000)
ax.grid(True, which="both", alpha=0.5)
plt.tight_layout()
plt.show()
3. モデルフィッティング
得られた伝達関数（特に共振付近）に対して、理論モデルをフィッティングします。 ここでは、単振動の伝達関数モデルを定義し、最小二乗法でパラメータ（\(A, f_0, Q\)）を求めます。
モデル式（ゲイン \(A\) を含む）:
\[H(f) = A \cdot \frac{f_0^2}{f_0^2 - f^2 + i (f f_0 / Q)}\]
[4]:
# --- Function Define ---
def resonator_model(f, amp, f0, Q):
# f: FrequencyArray
# : amp(Gain), f0(ResonanceFrequency), Q(Q )
# Calculate (Complex number)
numerator = amp * (f0**2)
denominator = (f0**2) - (f**2) + 1j * (f * f0 / Q)
return numerator / denominator
# --- Run ---
# Data 、
# Resonance (100Hz ~ 500Hz) Data 。
tf_crop = tf.crop(100, 500)
# Estimated ( Peak )
# f0=300 、
p0 = {"amp": 1.0, "f0": 300.0, "Q": 10.0}
# Run
# FrequencySeries .fit() method Complex numberData
result = tf_crop.fit(resonator_model, p0=p0)
print(result)
┌─────────────────────────────────────────────────────────────────────────┐
│ Migrad │
├──────────────────────────────────┬──────────────────────────────────────┤
│ FCN = 6.175 (χ²/ndof = 0.0) │ Nfcn = 118 │
│ EDM = 4.72e-06 (Goal: 0.0002) │ │
├──────────────────────────────────┼──────────────────────────────────────┤
│ Valid Minimum │ Below EDM threshold (goal x 10) │
├──────────────────────────────────┼──────────────────────────────────────┤
│ No parameters at limit │ Below call limit │
├──────────────────────────────────┼──────────────────────────────────────┤
│ Hesse ok │ Covariance accurate │
└──────────────────────────────────┴──────────────────────────────────────┘
┌───┬──────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐
│ │ Name │ Value │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit- │ Limit+ │ Fixed │
├───┼──────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤
│ 0 │ amp │ 0.934 │ 0.007 │ │ │ │ │ │
│ 1 │ f0 │ 299.994 │ 0.021 │ │ │ │ │ │
│ 2 │ Q │ 49.6 │ 0.5 │ │ │ │ │ │
└───┴──────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘
┌─────┬────────────────────────────┐
│ │ amp f0 Q │
├─────┼────────────────────────────┤
│ amp │ 4.33e-05 0 -2.30e-3 │
│ f0 │ 0 0.000447 -0.1e-3 │
│ Q │ -2.30e-3 -0.1e-3 0.243 │
└─────┴────────────────────────────┘
4. 結果の確認
フィッティング結果を数値で確認し、測定データに重ねてプロットします。 FitResult.plot() メソッドは、複素数データの場合、自動的にボード線図（振幅・位相）を描画します。
[5]:
# --- Result Display ---
print("--- Estimated Parameters ---")
print(f"Resonance Frequency (f0): {result.params['f0']:.4f} Hz (True: {f0_true})")
print(f"Quality Factor (Q): {result.params['Q']:.4f} (True: {Q_true})")
print(f"Gain (Amp): {result.params['amp']:.4f}")
# --- Plot ---
# result.plot() Complex number 、Amplitude Phase 2 Axes
axes = result.plot()
plt.show()
--- Estimated Parameters ---
Resonance Frequency (f0): 299.9943 Hz (True: 300.0)
Quality Factor (Q): 49.6456 (True: 50.0)
Gain (Amp): 0.9340
