-
Notifications
You must be signed in to change notification settings - Fork 12
/
xHilbertTransform.ts
128 lines (117 loc) · 4.02 KB
/
xHilbertTransform.ts
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
import { NumberArray } from 'cheminfo-types';
import FFT from 'fft.js';
import { nextPowerOfTwo, isPowerOfTwo } from '../utils';
import { xCheck } from './xCheck';
export interface XHilbertTransformOptions {
forceFFT?: boolean;
}
/**
* Performs the Hilbert transform
* @link https://en.wikipedia.org/wiki/Hilbert_transform
* @param array - Array containing values
* @param options
* @returns A new vector with 90 degree shift regarding the phase of the original function
*/
export function xHilbertTransform(
array: NumberArray,
options: XHilbertTransformOptions = {},
): Float64Array {
xCheck(array);
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);
}
}
/**
* Performs the discrete Hilbert transform using fast Fourier transform
* @param array - Array containing values
* @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 {
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(length * 2);
fft.transform(fftResult, complexSignal);
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];
}
const hilbertSignal = new Float64Array(length * 2);
fft.inverseTransform(hilbertSignal, fftResult);
const result = new Float64Array(length);
for (let i = 0; i < length; i++) {
result[i] = hilbertSignal[i * 2 + 1];
}
return result;
}
/**
* Performs the discrete Hilbert transform
* @param array - Array containing values
* @param options
* @returns A new vector with 90 degree shift regarding the phase of the original function
*/
function hilbertTransform(
array: NumberArray,
options: { inClockwise?: boolean } = {},
): Float64Array {
const { inClockwise = true } = options;
const input = [0, ...array, 0];
const result = new Float64Array(array.length);
for (let k = 1; k < input.length - 1; k++) {
let aSum = 0;
for (let i = 0; i < k - 1; i++) {
const log = Math.log((k - i) / (k - i - 1));
aSum += input[i] * log + (input[i + 1] - input[i]) * (-1 + (k - i) * log);
}
const b = input[k - 1] - input[k + 1];
let cSum = 0;
for (let i = k + 1; i < input.length - 1; i++) {
const log = Math.log((i - k) / (i - k + 1));
cSum += input[i] * log + (input[i - 1] - input[i]) * (1 + (i - k) * log);
}
result[k - 1] = ((inClockwise ? 1 : -1) * (aSum + b + cSum)) / Math.PI;
}
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: NumberArray, length: number): Float64Array {
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;
}