iOS ускоряет результат зеркального отображения низкочастотного фильтра FFT

Question

Я пытаюсь перенести существующий низкочастотный фильтр на основе FFT на iOS с использованием инфраструктуры Accelerate vDSP.

Похоже, что БПФ работает, как и ожидалось, для первой четверти образца. Но затем результаты кажутся неправильными, и еще более странные отражаются (причем большая часть сигнала отражает большую часть первой половины).

Вы можете увидеть результаты тестового приложения ниже. Сначала наносятся исходные данные выборки, затем пример ожидаемых отфильтрованных результатов (фильтрация сигнала выше 15 Гц), затем, наконец, результаты моего текущего кода БПФ (обратите внимание, что желаемые результаты и пример результата БПФ имеют разный масштаб, чем исходные данные):

Фактический код для моего фильтра нижних частот выглядит следующим образом:

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;
}

В исходном коде библиотека обрабатывала входные данные не двух степеней; в моем Ускоряющем коде я дополняю вход до ближайшей степени двух. В случае нижеприведенного выборочного теста исходные данные выборки составляют 1000 выборок, поэтому они дополняются до 1024. Я не думаю, что это повлияет на результаты, но я включаю их ради возможных различий.

Если вы хотите поэкспериментировать с решением, вы можете скачать пример проекта, который генерирует графики здесь (в папке FFTTest):

FFT Пример кода проекта

Спасибо за понимание, я раньше не работал с БПФ, поэтому чувствую, что упускаю что-то критическое.

Answer 1

Если вам нужен строго реальный (не сложный) результат, то данные перед IFFT должны быть сопряженно-симметричными. Если вы не хотите, чтобы результат был зеркально симметричным, то не обнуляйте мнимый компонент перед IFFT. Просто обнуление бункеров до того, как IFFT создаст фильтр с огромным количеством пульсаций в полосе пропускания.

Платформа Accelerate также поддерживает больше длин БПФ, чем просто степени 2.