feat: Allow xHilbertTransform to use FFT for array lengths non power …

…of 2 through and option (#203) feat: Allow xHilbertTransform to use FFT for array lengths non power of 2 through and option (#203) closes: #202
mljs · Dec 6, 2023 · 14c0dde · 14c0dde
1 parent 8153f7c
commit 14c0dde
Show file tree

Hide file tree

Showing 2 changed files with 161 additions and 87 deletions.
diff --git a/src/x/__tests__/xHilbertTransform.test.ts b/src/x/__tests__/xHilbertTransform.test.ts
@@ -1,79 +1,112 @@
-import { xHilbertTransform, xMaxValue } from '../../index';
+import { xHilbertTransform, xMaxValue } from '../..';
 
-describe('test discrete hilbert transform', () => {
+test('test discrete hilbert transform', () => {
   const length = 50;
-  it('test hilbert transform of cos -> sin function', () => {
-    const cos = new Array(length)
-      .fill(0)
-      .map((_, i) => Math.cos((i * Math.PI) / 10));
-    const sin = new Array(length)
-      .fill(0)
-      .map((_, i) => Math.sin((i * Math.PI) / 10));
-    const trs = xHilbertTransform(cos);
-    const result = Array.from(trs);
-    // Excluding some points due to the edge effects
-    for (let i = 5; i < 45; i++) {
-      expect(result[i]).toBeCloseTo(sin[i], 1);
-    }
-  });
+  const cos = new Float64Array(length);
+  const sin = new Float64Array(length);
+  const minusCos = new Float64Array(length);
+  const t = 10;
+  for (let i = 0; i < length; i++) {
+    cos[i] = Math.cos((i * Math.PI) / t);
+    sin[i] = Math.sin((i * Math.PI) / t);
+    minusCos[i] = -Math.cos((i * Math.PI) / t);
+  }
+  const tcos = xHilbertTransform(cos);
+  const tsin = xHilbertTransform(sin);
+  // cos -> sin (Excluding some points due to the edge effects)
+  for (let i = 5; i < 45; i++) {
+    expect(tcos[i]).toBeCloseTo(sin[i], 1);
+  }
+  // sin -> -cos (Excluding some points due to the edge effects)
+  for (let i = 15; i < 35; i++) {
+    expect(tsin[i]).toBeCloseTo(minusCos[i], 1);
+  }
+});
+
+test('test fast hilbert transform', () => {
+  const length = 64;
+  const cos = new Float64Array(length);
+  const sin = new Float64Array(length);
+  const minusCos = new Float64Array(length);
+  const t = 32;
+  for (let i = 0; i < length; i++) {
+    cos[i] = Math.cos((i * Math.PI) / t);
+    sin[i] = Math.sin((i * Math.PI) / t);
+    minusCos[i] = -Math.cos((i * Math.PI) / t);
+  }
+  const tcos = xHilbertTransform(cos);
+  const tsin = xHilbertTransform(sin);
+  // cos -> sin
+  for (let i = 0; i < 64; i++) {
+    expect(tcos[i]).toBeCloseTo(sin[i], 6);
+  }
+  // sin -> -cos
+  for (let i = 0; i < 64; i++) {
+    expect(tsin[i]).toBeCloseTo(minusCos[i], 6);
+  }
+});
 
-  it('test hilbert transform of sin -> -cos function', () => {
-    const minusCos = new Array(length)
-      .fill(0)
-      .map((_, i) => -Math.cos((i * Math.PI) / 10));
-    const sin = new Array(length)
-      .fill(0)
-      .map((_, i) => Math.sin((i * Math.PI) / 10));
-    const trs = xHilbertTransform(sin);
-    const result = Array.from(trs);
-    // Excluding some points due to the edge effects
-    for (let i = 15; i < 35; i++) {
-      expect(result[i]).toBeCloseTo(minusCos[i], 1);
-    }
-  });
+test('test fast hilbert transform of squareWave function', () => {
+  const length = 64;
+  const p = 16;
+  const squareWave = new Float64Array(length);
+  for (let i = 0; i < length; i++) {
+    squareWave[i] = i % p < p / 2 ? 1 : -1;
+  }
+  const result = xHilbertTransform(squareWave);
+  const maxValue = xMaxValue(result);
+  for (let i = 0; i < length / p; i++) {
+    expect(result[i * p]).toStrictEqual(-maxValue);
+    expect(result[i * p + p * 0.5]).toStrictEqual(maxValue);
+    expect(result[i * p + p * 0.25]).toBeCloseTo(0, 10);
+    expect(result[i * p + p * 0.75]).toBeCloseTo(0, 10);
+  }
 });
 
-describe('test fast hilbert transform', () => {
-  const length = 2 ** 6;
-  it('test hilbert transform of cos -> sin function', () => {
-    const cos = new Float64Array(length).map((_, i) =>
-      Math.cos((i * Math.PI) / 32),
-    );
-    const sin = new Float64Array(length).map((_, i) =>
-      Math.sin((i * Math.PI) / 32),
-    );
-    const result = xHilbertTransform(cos);
-    for (let i = 0; i < 64; i++) {
-      expect(result[i]).toBeCloseTo(sin[i], 6);
-    }
-  });
+test('test fast hilbert transform with forceFFT (array length greater than a power of 2)', () => {
+  const length = 74;
+  const cos = new Float64Array(length);
+  const sin = new Float64Array(length);
+  const minusCos = new Float64Array(length);
+  for (let i = 0; i < length; i++) {
+    cos[i] = Math.cos((i * Math.PI) / (length / 2));
+    sin[i] = Math.sin((i * Math.PI) / (length / 2));
+    minusCos[i] = -Math.cos((i * Math.PI) / (length / 2));
+  }
+  const tcos = xHilbertTransform(cos, { forceFFT: true });
+  const tsin = xHilbertTransform(sin, { forceFFT: true });
+  // test hilbert transform of cos -> sin function
+  for (let i = 0; i < length; i++) {
+    expect(tcos[i]).toBeCloseTo(sin[i], 1);
+  }
+  // test hilbert transform of sin -> -cos function
+  for (let i = 0; i < length; i++) {
+    expect(tsin[i]).toBeCloseTo(minusCos[i], 1);
+  }
+});
 
-  it('test hilbert transform of sin -> -cos function', () => {
-    const minusCos = new Float64Array(length).map(
-      (_, i) => -Math.cos((i * Math.PI) / 32),
-    );
-    const sin = new Float64Array(length).map((_, i) =>
-      Math.sin((i * Math.PI) / 32),
-    );
-    const result = xHilbertTransform(sin);
-    for (let i = 0; i < 64; i++) {
-      expect(result[i]).toBeCloseTo(minusCos[i], 6);
-    }
-  });
+test('test fast hilbert transform with forceFFT (array length less than a power of 2)', () => {
+  const length = 54;
+  const x = new Float64Array(length);
+  const cos = new Float64Array(length);
+  const sin = new Float64Array(length);
+  const minusCos = new Float64Array(length);
+  for (let i = 0; i < length; i++) {
+    x[i] = i;
+    cos[i] = Math.cos((i * Math.PI) / (length / 2));
+    sin[i] = Math.sin((i * Math.PI) / (length / 2));
+    minusCos[i] = -Math.cos((i * Math.PI) / (length / 2));
+  }
 
-  it('test hilbert transform of squareWave function', () => {
-    const p = 2 ** 4;
-    const squareWave = new Float64Array(length);
-    for (let i = 0; i < length; i++) {
-      squareWave[i] = i % p < p / 2 ? 1 : -1;
-    }
-    const result = xHilbertTransform(squareWave);
-    const maxValue = xMaxValue(result);
-    for (let i = 0; i < length / p; i++) {
-      expect(result[i * p]).toStrictEqual(-maxValue);
-      expect(result[i * p + p * 0.5]).toStrictEqual(maxValue);
-      expect(result[i * p + p * 0.25]).toBeCloseTo(0, 10);
-      expect(result[i * p + p * 0.75]).toBeCloseTo(0, 10);
-    }
-  });
+  // force hilbert transform to use fft
+  const tcos = xHilbertTransform(cos, { forceFFT: true });
+  const tsin = xHilbertTransform(sin, { forceFFT: true });
+  // test hilbert transform of cos -> sin function
+  for (let i = 0; i < length; i++) {
+    expect(tcos[i]).toBeCloseTo(sin[i], 1);
+  }
+  // test hilbert transform of sin -> -cos function
+  for (let i = 0; i < length; i++) {
+    expect(tsin[i]).toBeCloseTo(minusCos[i], 1);
+  }
 });
diff --git a/src/x/xHilbertTransform.ts b/src/x/xHilbertTransform.ts
@@ -1,6 +1,8 @@
 import { DoubleArray } from 'cheminfo-types';
 import FFT from 'fft.js';
 
+import { nextPowerOfTwo, isPowerOfTwo } from '../utils';
+
 import { xCheck } from './xCheck';
 
 /**
@@ -10,10 +12,20 @@ import { xCheck } from './xCheck';
  * @returns A new vector with 90 degree shift regarding the phase of the original function
  */
 
-export function xHilbertTransform(array: DoubleArray) {
+export function xHilbertTransform(
+  array: DoubleArray,
+  options: { forceFFT?: boolean } = {},
+) {
   xCheck(array);
-  if (Math.log2(array.length) % 1 === 0) {
+  const { forceFFT = false } = options;
+  const length = array.length;
+  if (isPowerOfTwo(length)) {
     return hilbertTransformWithFFT(array);
+  } else if (forceFFT) {
+    return resampling(
+      hilbertTransformWithFFT(resampling(array, nextPowerOfTwo(length))),
+      length,
+    );
   } else {
     return hilbertTransform(array);
   }
@@ -25,27 +37,27 @@ export function xHilbertTransform(array: DoubleArray) {
  * @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"
  */
-export function hilbertTransformWithFFT(array: DoubleArray) {
-  const n = array.length;
-  const fft = new FFT(n);
-  const complexSignal = new Float64Array(n * 2);
-  for (let i = 0; i < n; i++) {
+function hilbertTransformWithFFT(array: DoubleArray) {
+  const length = array.length;
+  const fft = new FFT(length);
+  const complexSignal = new Float64Array(length * 2);
+  for (let i = 0; i < length; i++) {
     complexSignal[i * 2] = array[i];
   }
-  const fftResult = new Float64Array(n * 2);
+  const fftResult = new Float64Array(length * 2);
   fft.transform(fftResult, complexSignal);
-  const multiplier = new Float64Array(n);
-  for (let i = 1; i < n; i++) {
-    multiplier[i] = Math.sign(n / 2 - i);
+  const multiplier = new Float64Array(length);
+  for (let i = 1; i < length; i++) {
+    multiplier[i] = Math.sign(length / 2 - i);
   }
-  for (let i = 0; i < n; i++) {
+  for (let i = 0; i < length; i++) {
     fftResult[i * 2] *= multiplier[i];
     fftResult[i * 2 + 1] *= multiplier[i];
   }
-  const hilbertSignal = new Float64Array(n * 2);
+  const hilbertSignal = new Float64Array(length * 2);
   fft.inverseTransform(hilbertSignal, fftResult);
-  const result = new Float64Array(n);
-  for (let i = 0; i < n; i++) {
+  const result = new Float64Array(length);
+  for (let i = 0; i < length; i++) {
     result[i] = hilbertSignal[i * 2 + 1];
   }
   return result;
@@ -56,7 +68,7 @@ export function hilbertTransformWithFFT(array: DoubleArray) {
  * @param array - Array containing values
  * @returns A new vector with 90 degree shift regarding the phase of the original function
  */
-export function hilbertTransform(
+function hilbertTransform(
   array: DoubleArray,
   options: { inClockwise?: boolean } = {},
 ) {
@@ -79,3 +91,32 @@ export function hilbertTransform(
   }
   return result;
 }
+
+/**
+ * Performs resampling of an input array to the desired length employing linear interpolation.
+ * @param array - Array containing values.
+ * @param length - The length of the resulting array.
+ * @returns It returns a new array of the desired length.
+ * @link https://en.wikipedia.org/wiki/Sample-rate_conversion
+ */
+function resampling(array: DoubleArray, length: number) {
+  xCheck(array);
+  const oldLength = array.length;
+  const ratio = (oldLength - 1) / (length - 1);
+  const result = new Float64Array(length);
+
+  let currentIndex = 0;
+  let floor = Math.floor(currentIndex);
+  let ceil = Math.min(Math.ceil(currentIndex), oldLength - 1);
+  let diff = currentIndex - floor;
+
+  for (let i = 0; i < length; i++) {
+    result[i] = array[floor] * (1 - diff) + array[ceil] * diff;
+    currentIndex += ratio;
+    floor = Math.floor(currentIndex);
+    ceil = Math.min(Math.ceil(currentIndex), oldLength - 1);
+    diff = currentIndex - floor;
+  }
+
+  return result;
+}