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

Question

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

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

Результаты тестового приложения приведены ниже. Сначала изображены исходные данные сэмплирования, затем пример ожидаемых отфильтрованных результатов (сигнал фильтрации выше 15 Гц), затем, наконец, результаты моего текущего кода 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 я заполняю вход до ближайшей мощности двух. В случае теста образца ниже исходные данные образца составляют 1000 выборок, поэтому они дополняются до 1024. Я не думаю, что это повлияло бы на результаты, но я включаю это ради возможных различий.

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

FFT Example Project code

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

Answer 1

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

Ускорить система также поддерживает больше длины FFT, чем просто полномочия 2.