iOS Accelerate低通FFT滤波器镜像结果

Posted 2023-02-25

【中文标题】iOS Accelerate低通FFT滤波器镜像结果

【英文标题】：iOS Accelerate low-pass FFT filter mirroring result

【发布时间】：2013-05-27 18:35:38

【问题描述】：

我正在尝试使用 Accelerate vDSP 框架将现有的基于 FFT 的低通滤波器移植到 ios。

在样本的前 1/4 左右，FFT 似乎按预期工作。但在那之后结果似乎是错误的，甚至更奇怪的是镜像（信号的后半部分反映了前半部分的大部分）。

您可以从下面的测试应用程序中看到结果。首先绘制原始采样数据，然后是预期过滤结果的示例（过滤掉高于 15Hz 的信号），最后是我当前 FFT 代码的结果（请注意，预期结果和示例 FFT 结果的比例不同于原始数据）：

我的低通滤波器的实际代码如下：

double *lowpassFilterVector(double *accell, uint32_t sampleCount, double lowPassFreq, double sampleRate )

    double stride = 1;

    int ln = log2f(sampleCount);
    int n = 1 << ln;

    // So that we get an FFT of the whole data set, we pad out the array to the next highest power of 2.
    int fullPadN = n * 2;
    double *padAccell = malloc(sizeof(double) * fullPadN);
    memset(padAccell, 0, sizeof(double) * fullPadN);
    memcpy(padAccell, accell, sizeof(double) * sampleCount);

    ln = log2f(fullPadN);
    n = 1 << ln;

    int nOver2 = n/2;

    DSPDoubleSplitComplex A;
    A.realp = (double *)malloc(sizeof(double) * nOver2);
    A.imagp = (double *)malloc(sizeof(double) * nOver2);

    // This can be reused, just including it here for simplicity.
    FFTSetupD setupReal = vDSP_create_fftsetupD(ln, FFT_RADIX2);

    vDSP_ctozD((DSPDoubleComplex*)padAccell,2,&A,1,nOver2);

    // Use the FFT to get frequency counts
    vDSP_fft_zripD(setupReal, &A, stride, ln, FFT_FORWARD);


    const double factor = 0.5f;
    vDSP_vsmulD(A.realp, 1, &factor, A.realp, 1, nOver2);
    vDSP_vsmulD(A.imagp, 1, &factor, A.imagp, 1, nOver2);
    A.realp[nOver2] = A.imagp[0];
    A.imagp[0] = 0.0f;
    A.imagp[nOver2] = 0.0f;

    // Set frequencies above target to 0.

    // This tells us which bin the frequencies over the minimum desired correspond to
    NSInteger binLocation = (lowPassFreq * n) / sampleRate;

    // We add 2 because bin 0 holds special FFT meta data, so bins really start at "1" - and we want to filter out anything OVER the target frequency
    for ( NSInteger i = binLocation+2; i < nOver2; i++ )
    
        A.realp[i] = 0;
    

    // Clear out all imaginary parts
    bzero(A.imagp, (nOver2) * sizeof(double));
    //A.imagp[0] = A.realp[nOver2];


    // Now shift back all of the values
    vDSP_fft_zripD(setupReal, &A, stride, ln, FFT_INVERSE);

    double *filteredAccell = (double *)malloc(sizeof(double) * fullPadN);

    // Converts complex vector back into 2D array
    vDSP_ztocD(&A, stride, (DSPDoubleComplex*)filteredAccell, 2, nOver2);

    //  Have to scale results to account for Apple's FFT library algorithm, see:
    // http://developer.apple.com/library/ios/#documentation/Performance/Conceptual/vDSP_Programming_Guide/UsingFourierTransforms/UsingFourierTransforms.html#//apple_ref/doc/uid/TP40005147-CH202-15952
    double scale = (float)1.0f / fullPadN;//(2.0f * (float)n);
    vDSP_vsmulD(filteredAccell, 1, &scale, filteredAccell, 1, fullPadN);

    // Tracks results of conversion
    printf("\nInput & output:\n");
    for (int k = 0; k < sampleCount; k++)
    
        printf("%3d\t%6.2f\t%6.2f\t%6.2f\n", k, accell[k], padAccell[k], filteredAccell[k]);
    


    // Acceleration data will be replaced in-place.
    return filteredAccell;

在原始代码中，库处理非二次幂大小的输入数据；在我的 Accelerate 代码中，我将输入填充到最接近的 2 次幂。在示例测试的情况下，原始样本数据是 1000 个样本，因此它被填充到 1024。我认为这不会影响结果，但为了可能的差异，我将其包括在内。

如果您想试验解决方案，可以在此处下载生成图表的示例项目（在 FFTTest 文件夹中）：

FFT Example Project code

感谢您的任何见解，我以前没有使用过 FFT，所以我觉得我错过了一些关键的东西。

【问题讨论】：

您的方法的一个大问题（而不是您可能遇到的任何实现问题）是尝试像这样在频域中应用砖墙滤波器会在时域中产生巨大的伪影。您需要使用窗口化方法来避免这种情况。

你找到解决方案了吗？

【参考方案1】：

如果您想要一个严格实数（不复数）的结果，那么 IFFT 之前的数据必须是共轭对称的。如果您不希望结果镜像对称，则不要在 IFFT 之前将虚部归零。仅在 IFFT 之前将 bin 归零会在通带中创建一个具有大量波纹的滤波器。

Accelerate 框架还支持更多的 FFT 长度，而不仅仅是 2 的幂。