Python中的希尔伯特变换？

Posted 2023-03-12

【中文标题】Python中的希尔伯特变换？

【英文标题】：Hilbert Transform in Python?

【发布时间】：2019-10-16 05:52:21

【问题描述】：

我正在尝试从头开始编写希尔伯特变换，但不使用除 fft 和 ifft 之外的任何内置库。我不是专业的数学家，但我在网上找到了这两种用于希尔伯特变换的算法，一种在 C 中，一种在 MATLAB 中。我试图同时实现它们，但它们都没有给我与 SciPy 的 Hilbert 相同的结果。我肯定在我的实施中犯了一些错误。任何见解将不胜感激。

---

第一个实现： （来自 MATLAB 网站） 希尔伯特使用四步算法：

计算输入序列的FFT，将结果存储在向量x中。

创建一个向量 h，其元素 h(i) 具有以下值：

1 为i = 1, (n/2)+1

2 为i = 2, 3, ... , (n/2)

0 为i = (n/2)+2, ... , n

计算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() */

我的尝试：

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 返回的不是信号的希尔伯特变换，而是通过希尔伯特变换 (z = x + j * y) 计算的解析信号。输入信号是 x，y 是 x 的希尔伯特变换。

@peterfranciscook 嗨，彼得，Matlab 页面说它还返回分析信号“希尔伯特返回一个复杂的螺旋序列，有时称为分析信号”。这与它说“计算相应的分析序列”的c代码相同

那么......我会运行你的代码，看看我能找到什么

@peterfranciscook 非常感谢！

【参考方案1】：

我查看了您的代码，进行了一些编辑，并将其与 scipy 和 MATLAB Hilbert 变换进行了比较。函数hilbert_from_scratch返回一个复数序列；实分量是原始信号，复分量是希尔伯特变换。如果您只需要希尔伯特变换，请在返回的数组上使用 np.imag。

import math
from scipy.fftpack import *

def hilbert_from_scratch(u):
    # N : fft length
    # M : number of elements to zero out
    # U : DFT of u
    # v : IDFT of H(U)

    N = len(u)
    # take forward Fourier transform
    U = fft(u)
    M = N - N//2 - 1
    # zero out negative frequency components
    U[N//2+1:] = [0] * M
    # double fft energy except @ DC0
    U[1:N//2] = 2 * U[1:N//2]
    # take inverse Fourier transform
    v = ifft(U)
    return v


if __name__ == '__main__':
    N = 32
    f = 1
    dt = 1.0 / N
    y = []
    for n in range(N):
        x = 2*math.pi*f*dt*n
        y.append(math.sin(x))
    z1 = hilbert_from_scratch(y)
    z2 = hilbert(y)
    print(" n     y fromscratch scipy")
    for n in range(N):
        print(':2d :+5.2f :+10.2f :+5.2f'.format(n, y[n], z1[n], z2[n]))

输出：

 n     y fromscratch scipy
 0 +0.00 +0.00-1.00j +1.00
 1 +0.38 +0.38-0.92j +0.92
 2 +0.71 +0.71-0.71j +0.71
 3 +0.92 +0.92-0.38j +0.38
 4 +1.00 +1.00-0.00j +0.00
 5 +0.92 +0.92+0.38j -0.38
 6 +0.71 +0.71+0.71j -0.71
 7 +0.38 +0.38+0.92j -0.92
 8 +0.00 +0.00+1.00j -1.00
 9 -0.38 -0.38+0.92j -0.92
10 -0.71 -0.71+0.71j -0.71
11 -0.92 -0.92+0.38j -0.38
12 -1.00 -1.00+0.00j -0.00
13 -0.92 -0.92-0.38j +0.38
14 -0.71 -0.71-0.71j +0.71
15 -0.38 -0.38-0.92j +0.92

MATLAB：

>> y = sin(2*pi*linspace(0,1,17)); z = hilbert(y(1:end-1));
>> fprintf('%+2.2f %+2.2f\n',[real(z);imag(z)])
+0.00 -1.00
+0.38 -0.92
+0.71 -0.71
+0.92 -0.38
+1.00 -0.00
+0.92 +0.38
+0.71 +0.71
+0.38 +0.92
+0.00 +1.00
-0.38 +0.92
-0.71 +0.71
-0.92 +0.38
-1.00 +0.00
-0.92 -0.38
-0.71 -0.71
-0.38 -0.92

【讨论】：

非常感谢！