godosa-engine: initial public history

1 files changed, 51 insertions, 0 deletions

diff --git a/src/Godosa.Core/Audio/SignalStats.cs b/src/Godosa.Core/Audio/SignalStats.cs
new file mode 100644
index 0000000..34f30df
--- /dev/null
+++ b/src/Godosa.Core/Audio/SignalStats.cs
@@ -0,0 +1,51 @@
+namespace Godosa.Core.Audio;
+
+/// <summary>Signal measurements for audio tests and reports: zero crossings, RMS, spectral centroid.</summary>
+public static class SignalStats
+{
+    public static int ZeroCrossings(float[] x, int from = 0, int to = -1)
+    {
+        if (to < 0) to = x.Length;
+        int n = 0;
+        for (int i = from + 1; i < to; i++) if ((x[i - 1] < 0) != (x[i] < 0)) n++;
+        return n;
+    }
+
+    public static float Rms(float[] x, int from = 0)
+    {
+        double sum = 0;
+        for (int i = from; i < x.Length; i++) sum += x[i] * x[i];
+        return (float)Math.Sqrt(sum / (x.Length - from));
+    }
+
+    /// <summary>Power-weighted mean frequency (where the energy is) of the loudest 2048-sample
+    /// window, naive DFT. Magnitude weighting over-counts faint broadband tails.</summary>
+    public static float Centroid(short[] pcm, int rate = Pcm.SampleRate)
+    {
+        const int n = 2048;
+        int best = 0;
+        double bestEnergy = -1;
+        for (int start = 0; start + n <= Math.Max(pcm.Length, n); start += n / 2)
+        {
+            double e = 0;
+            for (int i = start; i < Math.Min(pcm.Length, start + n); i++) e += (double)pcm[i] * pcm[i];
+            if (e > bestEnergy) (bestEnergy, best) = (e, start);
+        }
+        double weighted = 0, total = 0;
+        for (int k = 1; k < n / 2; k++)
+        {
+            double re = 0, im = 0;
+            for (int i = 0; i < n; i++)
+            {
+                double v = best + i < pcm.Length ? pcm[best + i] : 0;
+                double a = 2 * Math.PI * k * i / n;
+                re += v * Math.Cos(a);
+                im -= v * Math.Sin(a);
+            }
+            double power = re * re + im * im;
+            weighted += power * k * rate / (double)n;
+            total += power;
+        }
+        return (float)(weighted / total);
+    }
+}