FFT で合成振動データの卓越周波数を求める

time_s,accel_m_s2
0.000,0
0.125,1
0.250,0
0.375,-1
0.500,0
0.625,1
0.750,0
0.875,-1

import csv
from pathlib import Path
import numpy as np

SOURCE = Path("vibration_data.csv")
OUTPUT_DIR = Path("outputs") / "fft_result"
OUTPUT = OUTPUT_DIR / "spectrum.csv"
EXPECTED_DOMINANT_HZ = 2.0
DT_TOLERANCE = 1e-12
FREQUENCY_TOLERANCE = 1e-12


def main() -> None:
    if not SOURCE.is_file():
        raise FileNotFoundError(f"Source CSV not found: {SOURCE}")
    if OUTPUT_DIR.exists():
        raise FileExistsError(f"Output folder already exists: {OUTPUT_DIR}")

    times = []
    values = []
    with SOURCE.open("r", encoding="utf-8-sig", newline="") as stream:
        reader = csv.DictReader(stream)
        required = {"time_s", "accel_m_s2"}
        if reader.fieldnames is None or not required.issubset(reader.fieldnames):
            raise ValueError("CSV is missing a required column.")

        for row in reader:
            times.append(float(row["time_s"]))
            values.append(float(row["accel_m_s2"]))

    if len(times) < 2:
        raise ValueError("At least two samples are required.")

    time_array = np.asarray(times, dtype=float)
    signal = np.asarray(values, dtype=float)
    intervals = np.diff(time_array)
    dt = intervals[0]

    if dt <= 0:
        raise ValueError("Sample interval must be positive.")
    if not np.allclose(intervals, dt, rtol=0.0, atol=DT_TOLERANCE):
        raise ValueError("Sampling interval is not constant.")

    sample_rate = 1.0 / dt
    centered = signal - signal.mean()
    spectrum = np.fft.rfft(centered)
    frequencies = np.fft.rfftfreq(len(centered), d=dt)
    magnitudes = np.abs(spectrum)

    if len(frequencies) < 2:
        raise ValueError("No nonzero frequency bin is available.")

    dominant_index = 1 + int(np.argmax(magnitudes[1:]))
    dominant_frequency = float(frequencies[dominant_index])

    if abs(dominant_frequency - EXPECTED_DOMINANT_HZ) > FREQUENCY_TOLERANCE:
        raise RuntimeError(
            f"Expected {EXPECTED_DOMINANT_HZ} Hz, got {dominant_frequency} Hz."
        )

    OUTPUT_DIR.parent.mkdir(parents=True, exist_ok=True)
    OUTPUT_DIR.mkdir()
    with OUTPUT.open("x", encoding="utf-8", newline="") as stream:
        writer = csv.writer(stream)
        writer.writerow(["frequency_hz", "fft_magnitude"])
        for frequency, magnitude in zip(frequencies, magnitudes):
            writer.writerow([f"{frequency:.6f}", f"{magnitude:.6f}"])

    resolution = sample_rate / len(centered)
    print(f"Samples: {len(centered)}.")
    print(f"Sampling frequency: {sample_rate:.1f} Hz.")
    print(f"Frequency resolution: {resolution:.1f} Hz.")
    print(f"Dominant nonzero frequency: {dominant_frequency:.1f} Hz.")
    print(f"Output: {OUTPUT.as_posix()}")


if __name__ == "__main__":
    main()

frequency_hz,fft_magnitude
0.000000,0.000000
1.000000,0.000000
2.000000,4.000000
3.000000,0.000000
4.000000,0.000000

Samples: 8.
Sampling frequency: 8.0 Hz.
Frequency resolution: 1.0 Hz.
Dominant nonzero frequency: 2.0 Hz.
Output: outputs/fft_result/spectrum.csv