我正在尝试从头开始编写 Hilbert 变换，但不使用任何内置库，除了用于 fft 和 ifft 的库。我不是数学家，但我在网上找到了两个 Hilbert 变换的算法，一个是用 C 语言编写的，另一个是用 MATLAB 编写的。我尝试实现它们，但都没有像 SciPy 的 Hilbert 变换那样给出相同的结果。我肯定在实现中犯了某些错误。非常感谢您的任何见解。
第一个实现: (来自MATLAB网站) Hilbert 使用四步算法：
计算输入序列的 FFT，将结果存储在向量
x中。创建一个向量
h，其元素h(i)具有以下值：当
i = 1, (n/2)+1时，h(i) = 1当
i = 2, 3, ... , (n/2)时，h(i) = 2当
i = (n/2)+2, ... , n时，h(i) = 0
计算
x和h的逐元素乘积。计算步骤 3 中得到的序列的逆 FFT，并返回结果的前
n个元素。
我的尝试：
def generate_array(n):
a = np.hstack((np.full(n//2+1, 2), np.zeros(n//2-1)))
a[[0, n//2]] = 1
return a
def hilbert_from_scratch_2(u):
fft_result = fft(u) #scipy fft
n = len(u)
to_multiply = generate_array(n)
result = np.multiply(n,to_multiply)
return ifft(result) #scipy ifft
第二个实现： （https://www.cfa.harvard.edu/~spaine/am/download/src/transform.c）
void hilbert(double *z, unsigned long n)
{
double x;
unsigned long i, n2;
n2 = n << 1;
/*
* Compute the (bit-reversed) Fourier transform of z.
*/
fft_dif(z, n);
/*
* Form the transform of the analytic sequence by zeroing
* the transform for negative time, except for the (N/2)th.
* element. Since z is now in bit-reversed order, this means
* zeroing every other complex element. The array indices of
* the elements to be zeroed are 6,7,10,11...etc. (The real
* and imaginary parts of the (N/2)th element are in z[2] and
* z[3], respectively.)
*/
for (i = 6; i < n2; i += 4) {
z[i] = 0.;
z[i+1] = 0.;
}
/*
* The 0th and (N/2)th elements get multiplied by 0.5. Test
* for the trivial 1-point transform, just in case.
*/
z[0] *= 0.5;
z[1] *= 0.5;
if (n > 1) {
z[2] *= 0.5;
z[3] *= 0.5;
}
/*
* Compute the inverse transform.
*/
ifft_dit(z, n);
/*
* Normalize the array. The factor of 2 is left over from
* forming the transform in the time domain.
*/
x = 2. / (double)n;
for (i = 0; i < n2; ++i)
z[i] *= x;
return;
} /* hilbert() */
My attempt:
def hilbert_from_scratch(signal):
fast_ft = fft(signal) #scipy fft
for i in range(6,len(signal),4):
fast_ft[i] = 0
fast_ft[i+1] = 0
fast_ft[0] = fast_ft[0]*.5
fast_ft[1] = fast_ft[1]*.5
if(len(fast_ft) > 1):
fast_ft[2] = fast_ft[2]*.5
fast_ft[3] = fast_ft[3]*.5
inverse_fft = ifft(fast_ft) #scipy ifft
x = 2 / len(signal)
for i in range(0,len(signal),1):
inverse_fft[i] = inverse_fft[i]*x
return inverse_fft
欢迎提供任何关于为什么它们都不能像SciPy的hilbert函数一样给出相同结果的见解。
scipy.signal.hilbert返回的不是一个信号的 Hilbert 变换，而是通过 Hilbert 变换计算出来的解析信号 (z = x + j * y)。输入信号是 x，y 是 x 的 Hilbert 变换。 - peterfranciscook