使用FFTW在C#中计算Hilbert变换

Question

我希望在C#中实现Hilbert变换。从这篇文章中，我看到最快的FFT开源实现似乎是FFTW，所以我下载了该示例并使用它来学习如何使用C#的fftw包装器。

我有一个包含200,000个点的当前信号用于测试。通过fft获取Hilbert变换相对简单：

1. 计算fft。 2. 除0和Nyquist分量（如果样本大小为偶数，则为0和n / 2 + 1）外，将所有正频率乘以2。 3. 将所有负频率（[n / 2 + 1，n]）乘以0。 4. 计算逆fft。

到目前为止，我已经完成了所有步骤。唯一的问题是逆fft。我无法通过fftw获得与Matlab的ifft相同的结果。

我的代码

    RealArray _input;
    ComplexArray _fft;

    void ComputeFFT()
    {
        _fft = new ComplexArray(_length / 2 + 1);
        _input.Set(Data);
        _plan = Plan.Create1(_length, _input, _fft, Options.Estimate);
        _plan.Execute();
    }

迄今为止，我仅使用正频率进行了fft。 因此，我不需要将负频率乘以零：它们甚至不存在。 使用以下代码，我可以获取我的原始信号：

    double[] ComputeIFFT(ComplexArray input)
    {
        double[] temp = new double[_length];
        RealArray output = new RealArray(_length);

        _plan = Plan.Create1(_length, input, output, Options.Estimate);
        _plan.Execute();
        temp = output.ToArray();

        for (int i = 0; i < _length; ++i)
        {
            temp[i] /= _length;
        }

        return temp;
    }

当我尝试从信号中获取复杂的反演时，问题就出现了。

    void ComputeHilbert()
    {
        double[] fft = FFT.ToArray();
        double[] h = new double[_length / 2 + 1];
        double[] temp = new double[_length * 2];

        bool fftLengthIsOdd = (_length | 1) == 1;

        h[0] = 1;

        for (int i = 1; i < _length / 2; i++) h[i] = 2;

        if (!fftLengthIsOdd) h[(_length / 2)] = 1;

        for (int i = 0; i <= _length / 2; i++)
        {
            temp[2 * i] = fft[2*i] * h[i];
            temp[2 * i + 1] = fft[2*i + 1] * h[i];
        }

        ComplexArray _tempHilbert = new ComplexArray(_length);

        _tempHilbert.Set(temp);

        _hilbert = ComputeIFFT(_tempHilbert);
        _hilbertComputed = true;
    }

重要的是需要注意，当我在一个ComplexArray对象上应用ToArray()方法时，得到的结果是一个长度为原数组两倍的double[]数组，其中实部和虚部是连续的。例如对于一个包含“3+1i”的ComplexArray对象，将得到一个双精度向量[3,1]。

因此，现在我有了以下形式的数据：[DC 频率，2*正频率，奈奎斯特频率，零点]。如果我将这些数据导出到Matlab并进行IFFT计算，我会得到与其hilbert(signal)相同的结果。

但是，如果我尝试使用fftw提供的IFFT，那么从奈奎斯特频率到结尾（也就是说，零点会影响fftw）得到的值会很奇怪。

这是我用来进行此操作的IFFT代码：

double[] ComputeIFFT(ComplexArray input)
{
    double[] temp;
    ComplexArray output = new ComplexArray(_length);

    _plan = Plan.Create1(_length, input, output, Direction.Backward, Options.Estimate);
    _plan.Execute();
    temp = output.ToArray();

    for (int i = 0; i < _length; ++i)
    {
        temp[i] /= _length;
    }

    return temp;
}

所以，简单来说，我的问题在于我使用的ifft计算方式。它似乎无法很好地处理零值。或许Matlab能够理解它必须采用一些不同的方法，而我应该手动执行，但我不知道该怎么做。

非常感谢您的帮助，提前致以感激！

- lyurealm

Answer 1

所以问题出在ComputeIFFT函数上。在for循环中，我使用的是i < _length，但是temp数组的长度为2 * _length，因为它包含实数和虚数值。

这就是为什么我只得到了一半的正确值。

正确的代码如下：

double[] ComputeIFFT(ComplexArray input)
{
    double[] temp;
    ComplexArray output = new ComplexArray(_length);

    _plan = Plan.Create1(_length, input, output, Direction.Backward, Options.Estimate);
    _plan.Execute();
    temp = output.ToArray();

    for (int i = 0; i < temp.Length; ++i)
    {
        temp[i] /= _length;
    }

    return temp;
}

我希望这对于任何试图在C#中通过FFTW实现Hilbert变换的人都有所帮助。

- lyurealm