fix

series/spectrum.go

@@ -0,0 +1,95 @@
+package series
+
+import (
+	"math"
+
+	"gonum.org/v1/gonum/dsp/fourier"
+)
+
+type SpectrumPoint struct {
+	Frequency float64 `json:"frequency"`
+	Amplitude float64 `json:"amplitude"`
+}
+
+func AverageDelta(times []float64) float64 {
+	if len(times) < 2 {
+		return 0
+	}
+	total := 0.0
+	count := 0
+	for i := 1; i < len(times); i++ {
+		delta := times[i] - times[i-1]
+		if delta <= 0 {
+			continue
+		}
+		total += delta
+		count++
+	}
+	if count == 0 {
+		return 0
+	}
+	return total / float64(count)
+}
+
+func SpectrumPoints(values []float64, dt float64) []SpectrumPoint {
+	if len(values) < 2 || dt <= 0 {
+		return []SpectrumPoint{}
+	}
+
+	centered := Centered(values)
+	fft := fourier.NewFFT(len(centered))
+	coeffs := fft.Coefficients(nil, centered)
+	limit := len(coeffs) / 2
+	points := make([]SpectrumPoint, 0, limit+1)
+	for i := 0; i <= limit; i++ {
+		points = append(points, SpectrumPoint{
+			Frequency: float64(i) / (float64(len(centered)) * dt),
+			Amplitude: math.Hypot(real(coeffs[i]), imag(coeffs[i])) / float64(len(centered)),
+		})
+	}
+	return points
+}
+
+func MainAngularFrequency(times, values []float64) float64 {
+	if len(values) < 2 || len(times) < 2 {
+		return 0
+	}
+
+	dt := AverageDelta(times)
+	if dt <= 0 {
+		return 0
+	}
+
+	centered := Centered(values)
+	fft := fourier.NewFFT(len(centered))
+	coeffs := fft.Coefficients(nil, centered)
+	if len(coeffs) < 2 {
+		return 0
+	}
+
+	bestIndex := 1
+	bestAmp := 0.0
+	for i := 1; i < len(coeffs); i++ {
+		amp := math.Hypot(real(coeffs[i]), imag(coeffs[i]))
+		if amp > bestAmp {
+			bestAmp = amp
+			bestIndex = i
+		}
+	}
+
+	frequencyHz := fft.Freq(bestIndex) / dt
+	return 2 * math.Pi * frequencyHz
+}
+
+func Centered(values []float64) []float64 {
+	centered := make([]float64, len(values))
+	mean := 0.0
+	for _, value := range values {
+		mean += value
+	}
+	mean /= float64(len(values))
+	for i, value := range values {
+		centered[i] = value - mean
+	}
+	return centered
+}