feat(dsp): add power & integral fns

- add power functions: - powerSumSquared() - powerMeanSquared() - powerTimeIntegral() - add integralF/T() - add isComplex() check
thi-ng · Dec 17, 2020 · 88edaac · 88edaac
1 parent a72dfc5
commit 88edaac
Show file tree

Hide file tree

Showing 3 changed files with 84 additions and 0 deletions.
diff --git a/packages/dsp/src/fft/power.ts b/packages/dsp/src/fft/power.ts
@@ -0,0 +1,75 @@
+import type { NumericArray } from "@thi.ng/api";
+import type { ComplexArray } from "../api";
+import { isComplex } from "../util/complex";
+
+/**
+ * Computes the `sum(f(i)^2)` for given array.
+ *
+ * @param window
+ */
+export const integralT = (window: NumericArray) => {
+    let sum = 0;
+    for (let i = window.length; --i >= 0; ) {
+        // prettier-ignore
+        sum += (<number>window[i]) ** 2;
+    }
+    return sum;
+};
+
+/**
+ * Computes the `sum(|c(i)|^2)` for given complex array.
+ *
+ * @param window
+ */
+export const integralF = ([real, img]: ComplexArray) => {
+    let sum = 0;
+    for (let i = real.length; --i >= 0; ) {
+        sum += real[i] ** 2 + img[i] ** 2;
+    }
+    return sum;
+};
+
+/**
+ * Computes sum squared power of given time or frequency domain window.
+ *
+ * @remarks
+ * References:
+ * - http://www.it.uom.gr/teaching/linearalgebra/NumericalRecipiesInC/c13-4.pdf
+ * - http://www.hep.ucl.ac.uk/~rjn/saltStuff/fftNormalisation.pdf
+ *
+ * @param window
+ */
+export const powerSumSquared = (window: NumericArray | ComplexArray) =>
+    isComplex(window)
+        ? integralF(window) / window[0].length
+        : integralT(window);
+
+/**
+ * Computes mean squared power of given time or frequency domain window.
+ *
+ * @remarks
+ * References:
+ * - http://www.it.uom.gr/teaching/linearalgebra/NumericalRecipiesInC/c13-4.pdf
+ * - http://www.hep.ucl.ac.uk/~rjn/saltStuff/fftNormalisation.pdf
+ *
+ * @param window
+ */
+export const powerMeanSquared = (window: NumericArray | ComplexArray) =>
+    powerSumSquared(window) /
+    (isComplex(window) ? window[0].length : window.length);
+
+/**
+ * Computes time-integral squared power of given time or frequency domain
+ * window.
+ *
+ * @remarks
+ * References:
+ * - http://www.it.uom.gr/teaching/linearalgebra/NumericalRecipiesInC/c13-4.pdf
+ * - http://www.hep.ucl.ac.uk/~rjn/saltStuff/fftNormalisation.pdf
+ *
+ * @param window
+ */
+export const powerTimeIntegral = (
+    window: NumericArray | ComplexArray,
+    fs: number
+) => (isComplex(window) ? integralF(window) : integralT(window)) / fs;
diff --git a/packages/dsp/src/index.ts b/packages/dsp/src/index.ts
@@ -53,8 +53,10 @@ export * from "./osc/tri";
 export * from "./osc/wavetable";
 
 export * from "./fft/fft";
+export * from "./fft/power";
 export * from "./fft/window";
 
 export * from "./util/anti-alias";
+export * from "./util/complex";
 export * from "./util/convert";
 export * from "./util/filter-response";
diff --git a/packages/dsp/src/util/complex.ts b/packages/dsp/src/util/complex.ts
@@ -0,0 +1,7 @@
+import type { NumericArray } from "@thi.ng/api";
+import { isNumber } from "@thi.ng/checks";
+import type { ComplexArray } from "../api";
+
+export const isComplex = (
+    buf: NumericArray | ComplexArray
+): buf is ComplexArray => !isNumber(buf[0]);