reimMatrixFFT

src/__tests__/__snapshots__/index.test.ts.snap

@@ -8,6 +8,7 @@ exports[`existence of exported functions 1`] = `
   "reimPhaseCorrection",
   "reimZeroFilling",
   "reimArrayFFT",
+  "reimMatrixFFT",
   "getOutputArray",
   "xAbsolute",
   "xAbsoluteMedian",
@@ -155,6 +156,7 @@ exports[`existence of exported functions 1`] = `
   "matrixCreateEmpty",
   "matrixCuthillMckee",
   "matrixGetSubMatrix",
+  "matrixHilbertTransform",
   "matrixHistogram",
   "matrixMaxAbsoluteZ",
   "matrixMedian",

src/index.ts

@@ -1,5 +1,6 @@
 export * from './reim/index.ts';
 export * from './reimArray/index.ts';
+export * from './reimMatrix/index.ts';
 
 export * from './x/index.ts';
 

src/matrix/__tests__/matrixHilbertTransform.test.ts

@@ -0,0 +1,61 @@
+import { expect, test } from 'vitest';
+
+import { xHilbertTransform } from '../../x/xHilbertTransform.ts';
+import { matrixHilbertTransform } from '../matrixHilbertTransform.ts';
+
+const row0 = Float64Array.from([1, 0, -1, 0, 1, 0, -1, 0]);
+const row1 = Float64Array.from([0, 1, 0, -1, 0, 1, 0, -1]);
+const row2 = Float64Array.from([1, 2, 3, 4, 3, 2, 1, 0]);
+
+test('matrixHilbertTransform: matches xHilbertTransform row by row', () => {
+  const rows = [row0, row1, row2];
+  const result = matrixHilbertTransform(rows);
+
+  for (let i = 0; i < rows.length; i++) {
+    expect(result[i]).toStrictEqual(xHilbertTransform(rows[i]));
+  }
+});
+
+test('matrixHilbertTransform: returns empty array for empty input', () => {
+  expect(matrixHilbertTransform([])).toStrictEqual([]);
+});
+
+test('matrixHilbertTransform: output arrays are independent (no shared buffers)', () => {
+  const rows = [row0, row1];
+  const result = matrixHilbertTransform(rows);
+
+  expect(result[0]).not.toBe(result[1]);
+  expect(result[0]).not.toBe(rows[0]);
+  expect(result[1]).not.toBe(rows[1]);
+});
+
+test('matrixHilbertTransform: throws RangeError when length is not a power of two', () => {
+  const rows = [Float64Array.from([1, 2, 3, 4, 5, 6])];
+
+  expect(() => matrixHilbertTransform(rows)).toThrowError(RangeError);
+  expect(() => matrixHilbertTransform(rows)).toThrowError(/power of two/);
+});
+
+test('matrixHilbertTransform: throws RangeError when rows have different lengths', () => {
+  const rows = [row0, Float64Array.from([1, 2, 3, 4])];
+
+  expect(() => matrixHilbertTransform(rows)).toThrowError(RangeError);
+  expect(() => matrixHilbertTransform(rows)).toThrowError(/row 1/);
+});
+
+test('matrixHilbertTransform inPlace: result shares references with input', () => {
+  const rows = [Float64Array.from(row0), Float64Array.from(row1)];
+  const result = matrixHilbertTransform(rows, { output: rows });
+
+  expect(result[0]).toBe(rows[0]);
+  expect(result[1]).toBe(rows[1]);
+});
+
+test('matrixHilbertTransform inPlace: produces same values as out-of-place', () => {
+  const rowsCopy = [Float64Array.from(row0), Float64Array.from(row1)];
+  const outOfPlace = matrixHilbertTransform([row0, row1]);
+  matrixHilbertTransform(rowsCopy, { output: rowsCopy });
+
+  expect(rowsCopy[0]).toStrictEqual(outOfPlace[0]);
+  expect(rowsCopy[1]).toStrictEqual(outOfPlace[1]);
+});

src/matrix/index.ts

@@ -10,6 +10,7 @@ export * from './matrixColumnsCorrelation.ts';
 export * from './matrixCreateEmpty.ts';
 export * from './matrixCuthillMckee.ts';
 export * from './matrixGetSubMatrix.ts';
+export * from './matrixHilbertTransform.ts';
 export * from './matrixHistogram.ts';
 export * from './matrixMaxAbsoluteZ.ts';
 export * from './matrixMedian.ts';

src/matrix/matrixHilbertTransform.ts

@@ -0,0 +1,85 @@
+import FFT from 'fft.js';
+
+import { isPowerOfTwo } from '../utils/index.ts';
+
+import { matrixCreateEmpty } from './matrixCreateEmpty.ts';
+
+export interface MatrixHilbertTransformOptions {
+  /**
+   * Float64Array[] variable to store the result of hilbert transform
+   */
+  output?: Float64Array[];
+}
+
+/**
+ * Apply the Hilbert transform to each row of a real-valued matrix, reusing a
+ * single FFT instance for all rows to avoid repeated instantiation overhead.
+ * All rows must have the same length, which must be a power of two.
+ * @param rows - array of real-valued Float64Array rows
+ * @param options - options
+ * @returns array of Hilbert-transformed Float64Array rows
+ * @see DOI: 10.1109/TAU.1970.1162139 "Discrete Hilbert transform"
+ */
+export function matrixHilbertTransform(
+  rows: Float64Array[],
+  options: MatrixHilbertTransformOptions = {},
+): Float64Array[] {
+  if (rows.length === 0) return [];
+
+  const size = rows[0].length;
+
+  if (!isPowerOfTwo(size)) {
+    throw new RangeError(`Row length must be a power of two. Got ${size}.`);
+  }
+
+  for (let j = 1; j < rows.length; j++) {
+    if (rows[j].length !== size) {
+      throw new RangeError(
+        `All rows must have the same length. Expected ${size} but row ${j} has length ${rows[j].length}.`,
+      );
+    }
+  }
+
+  // Single FFT instance reused across all rows
+  const fft = new FFT(size);
+
+  // Multiplier computed once — identical for every row of the same length
+  const half = size >> 1;
+
+  // Shared working buffers reused across all rows
+  const fftResult = new Float64Array(size * 2);
+  const hilbertSignal = new Float64Array(size * 2);
+
+  const {
+    output = matrixCreateEmpty({
+      nbRows: rows.length,
+      nbColumns: size,
+    }),
+  } = options;
+
+  for (let j = 0; j < rows.length; j++) {
+    const row = rows[j];
+
+    fft.realTransform(fftResult, row);
+    fft.completeSpectrum(fftResult);
+
+    const idx = half << 1;
+    fftResult[idx] = 0;
+    fftResult[idx + 1] = 0;
+
+    // Negate negative frequencies
+    for (let c = (half + 1) << 1; c < fftResult.length; c += 2) {
+      fftResult[c] = -fftResult[c];
+      fftResult[c + 1] = -fftResult[c + 1];
+    }
+
+    fft.inverseTransform(hilbertSignal, fftResult);
+
+    const result = output[j];
+    for (let i = 0; i < size; i++) {
+      result[i] = hilbertSignal[i * 2 + 1];
+    }
+  }
+
+  return output;
+}

src/reimMatrix/__tests__/reimMatrixFFT.test.ts

@@ -0,0 +1,143 @@
+import { expect, test } from 'vitest';
+
+import { reimMatrixFFT } from '../reimMatrixFFT.ts';
+
+test('reimMatrixFFT: round-trip on a single row', () => {
+  const data = {
+    re: [Float64Array.from([0, 3, 6, 5])],
+    im: [Float64Array.from([0, 4, 8, 3])],
+  };
+
+  const transformed = reimMatrixFFT(data, { applyZeroShift: true });
+  const restored = reimMatrixFFT(transformed, {
+    inverse: true,
+    applyZeroShift: true,
+  });
+
+  expect(restored.re[0]).toStrictEqual(data.re[0]);
+  expect(restored.im[0]).toStrictEqual(data.im[0]);
+});
+
+test('reimMatrixFFT: round-trip on multiple rows', () => {
+  const data = {
+    re: [
+      Float64Array.from([1, 0, 0, 0]),
+      Float64Array.from([0, 1, 0, 0]),
+      Float64Array.from([0, 3, 6, 5]),
+    ],
+    im: [
+      Float64Array.from([0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0]),
+      Float64Array.from([0, 4, 8, 3]),
+    ],
+  };
+
+  const transformed = reimMatrixFFT(data);
+  const restored = reimMatrixFFT(transformed, { inverse: true });
+
+  for (let i = 0; i < data.re.length; i++) {
+    expect(restored.re[i]).toStrictEqual(data.re[i]);
+    expect(restored.im[i]).toStrictEqual(data.im[i]);
+  }
+});
+
+test('reimMatrixFFT: applyZeroShift round-trip on multiple rows', () => {
+  const data = {
+    re: [Float64Array.from([0, 3, 6, 5]), Float64Array.from([1, 2, 3, 4])],
+    im: [Float64Array.from([0, 4, 8, 3]), Float64Array.from([0, 1, 0, 1])],
+  };
+
+  const transformed = reimMatrixFFT(data, { applyZeroShift: true });
+  const restored = reimMatrixFFT(transformed, {
+    inverse: true,
+    applyZeroShift: true,
+  });
+
+  for (let i = 0; i < data.re.length; i++) {
+    expect(restored.re[i]).toStrictEqual(data.re[i]);
+    expect(restored.im[i]).toStrictEqual(data.im[i]);
+  }
+});
+
+test('reimMatrixFFT: output arrays are independent (no shared buffers)', () => {
+  const data = {
+    re: [Float64Array.from([1, 0, 0, 0]), Float64Array.from([0, 1, 0, 0])],
+    im: [Float64Array.from([0, 0, 0, 0]), Float64Array.from([0, 0, 0, 0])],
+  };
+
+  const result = reimMatrixFFT(data);
+
+  expect(result.re[0]).not.toBe(result.re[1]);
+  expect(result.im[0]).not.toBe(result.im[1]);
+  expect(result.re[0]).not.toBe(data.re[0]);
+  expect(result.im[0]).not.toBe(data.im[0]);
+});
+
+test('reimMatrixFFT: returns empty matrix for empty input', () => {
+  expect(reimMatrixFFT({ re: [], im: [] })).toStrictEqual({ re: [], im: [] });
+});
+
+test('reimMatrixFFT: throws RangeError when rows have different lengths', () => {
+  const data = {
+    re: [
+      Float64Array.from([1, 0, 0, 0]),
+      Float64Array.from([0, 1, 0, 0, 0, 0, 0, 0]),
+    ],
+    im: [
+      Float64Array.from([0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0, 0, 0, 0, 0]),
+    ],
+  };
+
+  expect(() => reimMatrixFFT(data)).toThrowError(RangeError);
+});
+
+test('reimMatrixFFT: throws RangeError indicating which row has the wrong length', () => {
+  const data = {
+    re: [
+      Float64Array.from([1, 0, 0, 0]),
+      Float64Array.from([0, 1, 0, 0]),
+      Float64Array.from([0, 0, 1, 0, 0, 0, 0, 0]),
+    ],
+    im: [
+      Float64Array.from([0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0]),
+      Float64Array.from([0, 0, 0, 0, 0, 0, 0, 0]),
+    ],
+  };
+
+  expect(() => reimMatrixFFT(data)).toThrowError(/row 2/);
+});
+
+test('reimMatrixFFT inPlace: result shares re/im references with input', () => {
+  const data = {
+    re: [Float64Array.from([1, 0, 0, 0]), Float64Array.from([0, 3, 6, 5])],
+    im: [Float64Array.from([0, 0, 0, 0]), Float64Array.from([0, 4, 8, 3])],
+  };
+
+  const result = reimMatrixFFT(data, { inPlace: true });
+
+  expect(result.re[0]).toBe(data.re[0]);
+  expect(result.im[0]).toBe(data.im[0]);
+  expect(result.re[1]).toBe(data.re[1]);
+  expect(result.im[1]).toBe(data.im[1]);
+});
+
+test('reimMatrixFFT inPlace: round-trip restores original values', () => {
+  const data = {
+    re: [Float64Array.from([0, 3, 6, 5]), Float64Array.from([1, 2, 3, 4])],
+    im: [Float64Array.from([0, 4, 8, 3]), Float64Array.from([0, 1, 0, 1])],
+  };
+  const originals = {
+    re: data.re.map((row) => Float64Array.from(row)),
+    im: data.im.map((row) => Float64Array.from(row)),
+  };
+
+  reimMatrixFFT(data, { inPlace: true, applyZeroShift: true });
+  reimMatrixFFT(data, { inPlace: true, inverse: true, applyZeroShift: true });
+
+  for (let i = 0; i < data.re.length; i++) {
+    expect(data.re[i]).toStrictEqual(originals.re[i]);
+    expect(data.im[i]).toStrictEqual(originals.im[i]);
+  }
+});

src/reimMatrix/index.ts

@@ -0,0 +1 @@
+export * from './reimMatrixFFT.ts';