diff --git a/src/__tests__/__snapshots__/index.test.ts.snap b/src/__tests__/__snapshots__/index.test.ts.snap index 7330ba0b..d3b02335 100644 --- a/src/__tests__/__snapshots__/index.test.ts.snap +++ b/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", diff --git a/src/index.ts b/src/index.ts index abb3c006..113946a4 100644 --- a/src/index.ts +++ b/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'; diff --git a/src/matrix/__tests__/matrixHilbertTransform.test.ts b/src/matrix/__tests__/matrixHilbertTransform.test.ts new file mode 100644 index 00000000..5a3d2d9a --- /dev/null +++ b/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]); +}); diff --git a/src/matrix/index.ts b/src/matrix/index.ts index 30c7e11d..3c37bb06 100644 --- a/src/matrix/index.ts +++ b/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'; diff --git a/src/matrix/matrixHilbertTransform.ts b/src/matrix/matrixHilbertTransform.ts new file mode 100644 index 00000000..ac1b1460 --- /dev/null +++ b/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; +} diff --git a/src/reimMatrix/__tests__/reimMatrixFFT.test.ts b/src/reimMatrix/__tests__/reimMatrixFFT.test.ts new file mode 100644 index 00000000..ec4995a5 --- /dev/null +++ b/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]); + } +}); diff --git a/src/reimMatrix/index.ts b/src/reimMatrix/index.ts new file mode 100644 index 00000000..5ab93eda --- /dev/null +++ b/src/reimMatrix/index.ts @@ -0,0 +1 @@ +export * from './reimMatrixFFT.ts'; diff --git a/src/reimMatrix/reimMatrixFFT.ts b/src/reimMatrix/reimMatrixFFT.ts new file mode 100644 index 00000000..7ae3d5d1 --- /dev/null +++ b/src/reimMatrix/reimMatrixFFT.ts @@ -0,0 +1,89 @@ +import FFT from 'fft.js'; + +import { zeroShift } from '../reim/zeroShift.ts'; +import type { DataReImMatrix } from '../types/index.ts'; + +export interface ReimMatrixFFTOptions { + inverse?: boolean; + applyZeroShift?: boolean; + /** + * Write the result back into the input arrays instead of allocating new ones. + * @default false + */ + inPlace?: boolean; +} + +/** + * Apply FFT to each row of a complex matrix, reusing a single FFT instance for + * all rows to avoid repeated instantiation overhead. + * All rows must have the same length. + * @param data - complex matrix with re and im arrays of Float64Array rows + * @param options - options + * @returns complex matrix with transformed rows + */ +export function reimMatrixFFT( + data: DataReImMatrix, + options: ReimMatrixFFTOptions = {}, +): DataReImMatrix { + const { re, im } = data; + const numRows = re.length; + + if (numRows === 0) return { re: [], im: [] }; + + const { inverse = false, applyZeroShift = false, inPlace = false } = options; + + const size = re[0].length; + const csize = size << 1; + + for (let j = 0; j < numRows; j++) { + if (re[j].length !== size || im[j].length !== size) { + throw new RangeError( + `All rows must have the same length. Expected ${size} but row ${j} has length ${re[j].length}.`, + ); + } + } + + // Single FFT instance and working buffers reused across all rows + const fft = new FFT(size); + const complexArray = new Float64Array(csize); + const output = new Float64Array(csize); + + const resultRe = new Array(numRows); + const resultIm = new Array(numRows); + + for (let j = 0; j < numRows; j++) { + const reRow = re[j]; + const imRow = im[j]; + + for (let i = 0; i < csize; i += 2) { + complexArray[i] = reRow[i >>> 1]; + complexArray[i + 1] = imRow[i >>> 1]; + } + + const outRe = inPlace ? reRow : new Float64Array(size); + const outIm = inPlace ? imRow : new Float64Array(size); + + if (inverse) { + const input = applyZeroShift + ? zeroShift(complexArray, true) + : complexArray; + fft.inverseTransform(output, input); + for (let i = 0; i < csize; i += 2) { + outRe[i >>> 1] = output[i]; + outIm[i >>> 1] = output[i + 1]; + } + } else { + fft.transform(output, complexArray); + const source = applyZeroShift ? zeroShift(output) : output; + for (let i = 0; i < csize; i += 2) { + outRe[i >>> 1] = source[i]; + outIm[i >>> 1] = source[i + 1]; + } + } + + resultRe[j] = outRe; + resultIm[j] = outIm; + } + + return { re: resultRe, im: resultIm }; +} diff --git a/src/types/DataReImMatrix.ts b/src/types/DataReImMatrix.ts new file mode 100644 index 00000000..e80b2581 --- /dev/null +++ b/src/types/DataReImMatrix.ts @@ -0,0 +1,7 @@ +export interface DataReImMatrix { + /** Array of re rows */ + re: Float64Array[]; + + /** Array of im rows */ + im: Float64Array[]; +} diff --git a/src/types/index.ts b/src/types/index.ts index 7cf638b2..d2947713 100644 --- a/src/types/index.ts +++ b/src/types/index.ts @@ -1,3 +1,4 @@ export * from './DataReIm.ts'; +export * from './DataReImMatrix.ts'; export * from './DataXReIm.ts'; export * from './Point.ts'; diff --git a/src/x/xHilbertTransform.ts b/src/x/xHilbertTransform.ts index c5cb631e..11a79e36 100644 --- a/src/x/xHilbertTransform.ts +++ b/src/x/xHilbertTransform.ts @@ -42,29 +42,39 @@ export function xHilbertTransform( * @returns A new vector with 90 degree shift regarding the phase of the original function * @see DOI: 10.1109/TAU.1970.1162139 "Discrete Hilbert transform" */ -function hilbertTransformWithFFT( - array: NumberArray, -): Float64Array { +function hilbertTransformWithFFT(array: NumberArray) { const length = array.length; const fft = new FFT(length); - const fftResult = new Float64Array(length * 2); - fft.realTransform(fftResult, array); - fft.completeSpectrum(fftResult); - const multiplier = new Float64Array(length); - for (let i = 1; i < length; i++) { - multiplier[i] = Math.sign(length / 2 - i); - } - for (let i = 0; i < length; i++) { - fftResult[i * 2] *= multiplier[i]; - fftResult[i * 2 + 1] *= multiplier[i]; + // Single reusable buffer for FFT spectrum + const spectrum = new Float64Array(length * 2); + + // Forward FFT + fft.realTransform(spectrum, array); + fft.completeSpectrum(spectrum); + + const half = length >> 1; + + // Zero Nyquist + const j = half << 1; + spectrum[j] = 0; + spectrum[j + 1] = 0; + + // Negate negative frequencies + for (let j = (half + 1) << 1; j < spectrum.length; j += 2) { + spectrum[j] = -spectrum[j]; + spectrum[j + 1] = -spectrum[j + 1]; } + const hilbertSignal = new Float64Array(length * 2); - fft.inverseTransform(hilbertSignal, fftResult); + fft.inverseTransform(hilbertSignal, spectrum); + + // Extract imaginary part directly into output const result = new Float64Array(length); for (let i = 0; i < length; i++) { result[i] = hilbertSignal[i * 2 + 1]; } + return result; }