Note
このページは Jupyter Notebook から生成されました。 ノートブックをダウンロード (.ipynb)
FrequencySeries の新機能と高度な操作
このノートブックでは、gwexpy で拡張された FrequencySeries クラスの新しいメソッドと機能について紹介します。 主に複素スペクトルの扱い、微積分、フィルタリング（スムージング）、および他ライブラリとの連携機能に焦点を当てます。
[1]:
import matplotlib.pyplot as plt
import numpy as np
from gwexpy.frequencyseries import FrequencySeries
from gwexpy.plot import Plot
from gwexpy.timeseries import TimeSeries
plt.rcParams["figure.figsize"] = (10, 6)
1. データの準備
まずは TimeSeries から FFT を用いて FrequencySeries を作成します。 ここでは、特定の周波数成分を持つテスト信号を生成します。
[2]:
fs = 1024
t = np.arange(0, 4, 1 / fs)
exp = np.exp(-t / 1.5)
exp[: int(exp.size / 4)] = 0
data = (
np.sin(2 * np.pi * 10.1 * t)
+ 5 * exp * np.sin(2 * np.pi * 100.1 * t)
+ np.random.normal(scale=0.3, size=len(t))
)
ts = TimeSeries(data, dt=1 / fs, unit="um", name="Test Signal")
print(ts)
ts.plot(title=ts.name)
# FFT Run FrequencySeries (transient Mode )
spec = ts.fft(mode="transient", pad_left=1.0, pad_right=1.0, nfft_mode="next_fast_len")
print(f"Type: {type(spec)}")
print(f"Length: {len(spec)}")
print(f"df: {spec.df}")
TimeSeries([-0.19228838, 0.29767237, -0.04502141, ...,
0.75723094, 0.88698536, 0.8157283 ],
unit: um,
t0: 0.0 s,
dt: 0.0009765625 s,
name: Test Signal,
channel: None)
Type: <class 'gwexpy.frequencyseries.frequencyseries.FrequencySeries'>
Length: 3073
df: 0.16666666666666666 Hz
2. 複素スペクトルの可視化と変換
位相と振幅
phase(), degree(), to_db() メソッドを使用すると、複素スペクトルを直感的な単位に変換できます。
[3]:
# Amplitude dB Transform (ref=1.0, 20*log10)
spec_db = spec.to_db()
# Phase ( 、unwrap=True Continuous )
spec_phase = spec.degree(unwrap=True)
plot = Plot(spec_db, spec_phase, separate=True, sharex=True, xscale="log")
ax = plot.axes
ax[0].set_ylabel("Magnitude [dB (m)]")
ax[0].grid(True, which="both")
ax[1].set_ylabel("Phase [deg]")
ax[1].set_xlabel("Frequency [Hz]")
ax[1].grid(True, which="both")
plot.figure.suptitle("Complex Spectrum Analysis")
plot.show()
3. 周波数ドメインでの微積分
differentiate() および integrate() メソッドにより、周波数ドメインで微分・積分を行うことができます。 引数 order で階数を指定できます（デフォルトは1）。 これは「変位・速度・加速度」の変換（\((2 \pi i f)^n\) の乗算・除算）を簡単に行うための機能です。
[4]:
# (m) -> Velocity (m/s) Differentiate (order=1)
vel_spec = spec.differentiate()
# (m) -> Acceleration (m/s^2) 2 Differentiate (order=2)
accel_spec = spec.differentiate(order=2)
# Integrate : Acceleration -> Velocity
vel_from_accel = accel_spec.integrate()
plot = Plot(
spec.abs(), vel_spec.abs(), accel_spec.abs(), xscale="log", yscale="log", alpha=0.8
)
ax = plot.gca()
ax.get_lines()[0].set_label("Displacement [m]")
ax.get_lines()[1].set_label("Velocity [m/s]")
ax.get_lines()[2].set_label("Acceleration [m/s^2]")
ax.legend()
ax.grid(True, which="both")
ax.set_title("Calculus in Frequency Domain")
ax.set_ylabel("Magnitude")
plot.show()
4. スペクトルのスムージングとピーク検出
スムージング
smooth() メソッドを使用すると、移動平均などによるスペクトルの平滑化が可能です。
[5]:
# Amplitude 11 Sample Smoothing
spec_smooth = spec.smooth(width=11)
plot = Plot(spec.abs(), spec_smooth.abs(), xscale="log", yscale="log")
ax = plot.gca()
ax.get_lines()[0].set_label("Original")
ax.get_lines()[0].set_alpha(0.6)
ax.get_lines()[1].set_label("Smoothed (width=11)")
ax.get_lines()[1].set_color("red")
ax.legend()
ax.grid(True, which="both")
ax.set_title("Spectrum Smoothing")
plot.show()
ピーク検出
find_peaks() メソッドは scipy.signal.find_peaks をラップしており、特定の閾値を超えるピークを簡単に抽出できます。
[6]:
# Amplitude 0.2 Peak
peaks, props = spec.find_peaks(threshold=0.2)
plot = Plot(spec.abs())
ax = plot.gca()
ax.plot(
peaks.abs(), color="red", marker=".", ms=15, lw=0, zorder=3, label="Detected Peaks"
)
ax.set_xlim(1, 150)
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_title("Peak Detection")
ax.legend()
plot.show()
5. 高度な解析機能
群遅延 (Group Delay)
group_delay() メソッドは、位相の周波数微分から群遅延（信号のエンベロープの遅延）を計算します。
[7]:
gd = spec.group_delay()
plot = Plot(gd)
ax = plot.gca()
ax.set_ylabel("Group Delay [s]")
ax.set_xlabel("Frequency [Hz]")
ax.set_xlim(0, 200)
ax.set_title("Group Delay Calculation")
plot.show()
/home/washimi/work/gwexpy/.venv-docs-exec/lib/python3.12/site-packages/gwpy/plot/axes.py:235: UserWarning: Attempt to set non-positive xlim on a log-scaled axis will be ignored.
return super().set_xlim(left=left, right=right, **kwargs)
逆FFT (ifft)
ifft() メソッドは、TimeSeries を返します。mode="transient" で FFT した結果であっても、情報を引き継いで元の長さに戻す (trim=True) などの制御が可能です。
[8]:
# FFT TimeSeries
# mode="auto" 、 FrequencySeries transient
inv_ts = spec.ifft(mode="auto")
red_ts = inv_ts - ts
plot = Plot(ts, inv_ts, red_ts)
ax = plot.gca()
ax.get_lines()[0].set_label("Original")
ax.get_lines()[1].set_label("IFFT Result")
ax.get_lines()[2].set_label("Residual")
ax.get_lines()[1].set_linestyle("--")
ax.get_lines()[1].set_alpha(0.8)
ax.get_lines()[2].set_color("black")
ax.legend()
ax.set_title("Time Domain Round-trip (FFT -> IFFT)")
plot.show()
red_ts.plot(title="Residual Time Series after IFFT");
6. 他ライブラリとの連携
Pandas, xarray, control ライブラリとの相互変換が追加されています。
[9]:
# Pandas Series Transform
pd_series = spec.to_pandas()
print("Pandas index sample:", pd_series.index[:5])
display(pd_series)
# xarray DataArray Transform
try:
da = spec.to_xarray()
print("xarray coord name:", list(da.coords))
display(da)
except ImportError:
print("xarray not installed, skipping DataArray conversion.")
# control.FRD Transform (Control )
try:
from control import FRD
_ = FRD
frd_obj = spec.to_control_frd()
print("Successfully converted to control.FRD")
display(frd_obj)
except ImportError:
print("python-control library not installed")
Pandas index sample: Index([0.0, 0.16666666666666666, 0.3333333333333333, 0.5, 0.6666666666666666], dtype='float64', name='frequency')
frequency
0.000000 0.007069+0.000000j
0.166667 0.000065-0.006216j
0.333333 -0.001391+0.002865j
0.500000 -0.010751-0.001078j
0.666667 -0.002985-0.001148j
...
511.333333 0.001287+0.000084j
511.500000 -0.003222+0.000466j
511.666667 0.001020-0.004634j
511.833333 -0.001818+0.007399j
512.000000 0.001873+0.000000j
Name: Test Signal, Length: 3073, dtype: complex128
xarray coord name: ['frequency']
<xarray.DataArray 'Test Signal' (frequency: 3073)> Size: 49kB
array([ 7.06885196e-03+0.j , 6.52541232e-05-0.00621646j,
-1.39066175e-03+0.00286501j, ...,
1.01996891e-03-0.00463402j, -1.81789723e-03+0.00739938j,
1.87312900e-03+0.j ], shape=(3073,))
Coordinates:
* frequency (frequency) float64 25kB 0.0 0.1667 0.3333 ... 511.7 511.8 512.0
Attributes:
unit: um
channel: None
epoch: 0.0
Successfully converted to control.FRD
FrequencyResponseData(
array([[[ 7.06885196e-03+0.j ,
6.52541232e-05-0.00621646j,
-1.39066175e-03+0.00286501j, ...,
1.01996891e-03-0.00463402j,
-1.81789723e-03+0.00739938j,
1.87312900e-03+0.j ]]], shape=(1, 1, 3073)),
array([0.00000000e+00, 1.66666667e-01, 3.33333333e-01, ...,
5.11666667e+02, 5.11833333e+02, 5.12000000e+02],
shape=(3073,)),
outputs=1, inputs=1)
7. Python Control Library との連携
制御工学の分野で標準的な control ライブラリの Frequency Response Data (FRD) オブジェクトと相互変換が可能です。 これにより、GWExPyで計測した伝達関数を、制御系の設計や解析に直接利用できます。
[10]:
try:
import control
# FrequencySeries -> control.FRD Transform
# frequency_unit="Hz" do 、 rad/s Transform
frd_sys = spec.to_control_frd(frequency_unit="Hz")
print("\n--- Converted to Control FRD ---")
print(str(frd_sys)[:1000] + "\n... (truncated) ...")
# Plot (control )
control.bode(frd_sys) # (Plot Run )
# control.FRD -> FrequencySeries
fs_restored = FrequencySeries.from_control_frd(frd_sys, frequency_unit="Hz")
print("\n--- Restored FrequencySeries ---")
print(fs_restored)
except ImportError:
print("Python Control Systems Library is not installed.")
--- Converted to Control FRD ---
<FrequencyResponseData>: sys[1]
Inputs (1): ['u[0]']
Outputs (1): ['y[0]']
Freq [rad/s] Response
------------ ---------------------
0.000 0.007069 +0j
0.167 6.525e-05 -0.006216j
0.333 -0.001391 +0.002865j
0.500 -0.01075 -0.001078j
0.667 -0.002985 -0.001148j
0.833 0.001477 +0.007722j
1.000 0.01269 -0.004297j
1.167 0.0007285 -0.004045j
1.333 -0.001489-3.472e-05j
1.500 -0.006353 +0.006026j
1.667 -0.01012-0.0009878j
1.833 0.003597 -0.006665j
2.000 0.01067 +0.006731j
2.167 0.01137 +0.001391j
2.333 -0.01063 -0.004407j
2.500 -0.01006 -0.004614j
2.667 -0.001953 +0.004815j
2.833 -0.0009382 +0.006631j
3.000 0.01713 -0.004492j
3.167 0.00139 -0.00142j
3.333 -0.006361 -0.00638j
3.500 -0.003776 +0.005349j
3.667 -0.00991 +0.003714j
3.833 0.0002422 +0.004172j
... (truncated) ...
--- Restored FrequencySeries ---
FrequencySeries([ 7.06885196e-03+0.j ,
6.52541232e-05-0.00621646j,
-1.39066175e-03+0.00286501j, ...,
1.01996891e-03-0.00463402j,
-1.81789723e-03+0.00739938j,
1.87312900e-03+0.j ],
unit: dimensionless,
f0: 0.0 Hz,
df: 0.02652582384864922 Hz,
epoch: None,
name: None,
channel: None)
8. 求積和 (Quadrature Sum)
直交位相の和を計算する機能です。
[11]:
# Generate noisy data
np.random.seed(42)
f = spec.frequencies.value
noise = np.abs(np.random.randn(f.size))
peak = 10.0 * np.exp(-((f - 200) ** 2) / 50.0)
data = noise + peak
raw = FrequencySeries(data, frequencies=f, unit="V", name="Raw Data")
# Smooth
smoothed = raw.smooth(width=10, method="amplitude")
smoothed.name = "Smoothed"
# Convert to dB
raw_db = raw.to_db()
smoothed_db = smoothed.to_db()
raw_db.name = "Raw (dB)"
smoothed_db.name = "Smoothed (dB)"
plot = raw_db.plot(label="Raw", title="Smoothing & dB")
ax = plot.gca()
ax.plot(smoothed_db, label="Smoothed", linewidth=2)
ax.legend()
plot.show()
plt.close()
# Quadrature Sum (Noise Budget example)
noise_a = FrequencySeries(np.ones_like(f), frequencies=f, unit="V", name="Noise A")
noise_b = FrequencySeries(np.ones_like(f) * 2, frequencies=f, unit="V", name="Noise B")
total = noise_a.quadrature_sum(noise_b)
print(f"Noise A: {noise_a.value[0]}, Noise B: {noise_b.value[0]}")
print(f"Total (Sqrt Sum): {total.value[0]}")
Plot(total, noise_a, noise_b, alpha=0.8)
plt.legend(["Total", "Noise A", "Noise B"])
plt.xscale("log")
plt.yscale("log")
plt.show()
Noise A: 1.0, Noise B: 2.0
Total (Sqrt Sum): 2.23606797749979
