iOS Accelerate low-pass FFT filter mirroring result

Question

I am trying to port an existing FFT based low-pass filter to iOS using the Accelerate vDSP framework.

It seems like the FFT works as expected for about the first 1/4 of the sample. But then after that the results seem wrong, and even more odd are mirrored (with the last half of the signal mirroring most of the first half).

You can see the results from a test application below. First is plotted the original sampled data, then an example of the expected filtered results (filtering out signal higher than 15Hz), then finally the results of my current FFT code (note that the desired results and example FFT result are at a different scale than the original data):

FFT Results

The actual code for my low-pass filter is as follows:

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

In the original code the library was handling non power-of-two sizes of input data; in my Accelerate code I am padding out the input to the nearest power of two. In the case of the sample test below the original sample data is 1000 samples so it's padded to 1024. I don't think that would affect results but I include that for the sake of possible differences.

If you want to experiment with a solution, you can download the sample project that generates the graphs here (in the FFTTest folder):

FFT Example Project code

Thanks for any insight, I've not worked with FFT's before so I feel like I am missing something critical.

One big problem with your method (rather than whatever implementation problems you might be having) is that trying to apply a brick wall filter in the frequency domain like this will generate huge artefacts in the time domain. You need to use a windowing method to avoid this. — Paul R

hotpaw2 hotpaw2 · Accepted Answer · 2013-05-28T05:12:38

If you want a strictly real (not complex) result, then the data before the IFFT must be conjugate symmetric. If you don't want the result to be mirror symmetric, then don't zero the imaginary component before the IFFT. Merely zeroing bins before the IFFT creates a filter with a huge amount of ripple in the passband.

The Accelerate framework also supports more FFT lengths than just powers of 2.

iOS Accelerate low-pass FFT filter mirroring result

1 Answers