Python实现Hilbert变换遇问题：自研代码与SciPy结果不符求解析

Let's break down exactly where each of your implementations went off-track compared to SciPy's hilbert function. I'll go through them one by one—small missteps in FFT-based transforms can throw off results completely, so we'll focus on the key mistakes first.

Issue 1: MATLAB-Inspired Implementation #

Your core logic for the frequency-domain multiplier (generate_array) is actually correct for MATLAB's Hilbert transform approach—but you made a critical mistake in combining it with the FFT result.

Looking at your code:

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)  # THIS IS THE PROBLEM
    return ifft(result) #scipy ifft

You're multiplying the signal length n by to_multiply instead of multiplying the actual FFT result by to_multiply. This completely ignores the frequency content of your input signal! No wonder it doesn't match SciPy's output.

Fixed Version: #

import numpy as np
from scipy.fft import fft, ifft

def generate_array(n):
    a = np.zeros(n)
    a[0] = 1  # DC component
    if n % 2 == 0:
        a[1:n//2] = 2  # Positive frequencies (excluding Nyquist)
        a[n//2] = 1    # Nyquist frequency component
    else:
        a[1:(n+1)//2] = 2  # Positive frequencies for odd-length signals
    return a

def hilbert_from_scratch_2(u):
    fft_result = fft(u)
    n = len(u)
    to_multiply = generate_array(n)
    # Multiply FFT result by the frequency-domain mask (not n!)
    transformed_fft = fft_result * to_multiply
    return ifft(transformed_fft)

Issue 2: C-Language-Inspired Implementation #

The C code you referenced uses a real/imaginary interleaved array format (e.g., z[0] = real part of first bin, z[1] = imaginary part, z[2] = real part of second bin, etc.) for its FFT output. But SciPy's fft returns a complex-valued array (each element is a single complex number with real and imaginary parts). This mismatch in data structures is why your index-based zeroing is wrong.

Other problems in this implementation:

1. You're using bit-reversed FFT indices from the C code, but SciPy's fft returns frequencies in standard order (not bit-reversed).
2. The normalization step (x = 2 / len(signal)) conflicts with SciPy's ifft, which already applies a 1/n normalization factor automatically.

Fixed Version: #

def hilbert_from_scratch(signal):
    n = len(signal)
    fast_ft = fft(signal)
    
    # Zero out all negative frequency bins (standard FFT order)
    if n % 2 == 0:
        # For even-length signals: keep DC and Nyquist, zero rest of negative frequencies
        fast_ft[n//2 + 1 :] = 0
    else:
        # For odd-length signals: zero all bins starting from (n+1)//2
        fast_ft[(n+1)//2 :] = 0
    
    # Scale positive frequencies (excluding DC/Nyquist) by 2, matching Hilbert transform requirements
    fast_ft[1 : (n//2 if n%2==0 else (n+1)//2)] *= 2
    
    # No extra normalization needed—Scipy's ifft handles 1/n automatically
    return ifft(fast_ft)

Why This Works #

SciPy's hilbert function generates the analytic signal by:

1. Computing the FFT of the input.
2. Zeroing all negative frequency components.
3. Scaling positive frequency components (excluding DC and Nyquist) by 2.
4. Applying the inverse FFT.

Both fixed implementations align with this logic, so they should now produce results nearly identical to scipy.signal.hilbert (minor floating-point differences aside).

内容的提问来源于stack exchange，提问作者Brandan B