Search code examples
iosiphonesignal-processingfftaccelerate-framework

How to get iOS Accelerate FFT results in frequency magnitudes?


I use this code (based on Apple's audioRouch sample):

void FFTHelper::ComputeFFT(Float32* inAudioData, Float32* outFFTData)
{
    if (inAudioData == NULL || outFFTData == NULL) return;

    // Generate a split complex vector from the real data
    vDSP_ctoz((COMPLEX *)inAudioData, 2, &mDspSplitComplex, 1, mFFTLength);

    // Take the fft and scale appropriately
    vDSP_fft_zrip(mSpectrumAnalysis, &mDspSplitComplex, 1, mLog2N, kFFTDirection_Forward);
    vDSP_vsmul(mDspSplitComplex.realp, 1, &mFFTNormFactor, mDspSplitComplex.realp, 1, mFFTLength);
    vDSP_vsmul(mDspSplitComplex.imagp, 1, &mFFTNormFactor, mDspSplitComplex.imagp, 1, mFFTLength);

    // Zero out the nyquist value
    mDspSplitComplex.imagp[0] = 0.0;

    // Complex vector magnitudes squared; single precision.
    // Calculates the squared magnitudes of complex vector A.
    vDSP_zvmags(&mDspSplitComplex, 1, outFFTData, 1, mFFTLength);

}

To compute FFT on the simplest possible - 1Hz sinus wave (shifted up by 1 unit):

    Float32 waveFreq            = 1.0;
    int     samplesCount        = 1024;
    Float32 samplesPerSecond    = 1000;        //sample rate
    Float32 dt = 1 / samplesPerSecond;
    Float32 sd = M_PI * 2.0 * waveFreq;

    FFTHelper *mFFTHelper = new FFTHelper(samplesCount);

    Float32 NyquistMaxFreq  = samplesPerSecond/2.0;
    Float32 fftDataSize     = samplesCount/2.0;

    Float32 *sinusoidOriginal = (Float32 *)malloc(sizeof(Float32) * samplesCount);
    Float32 *outFFTData = (Float32 *)malloc(sizeof(Float32) * fftDataSize);

    // 2. Generate sin samples:
    for (int i = 0; i < samplesCount; i++) {

        Float32 x = dt * i;
        sinusoidOriginal[i] = sin(sd * x) + 1;
        [originalPlot addVector2D:GLVector2DMake(x, sinusoidOriginal[i])];
    }

    mFFTHelper->ComputeFFT(sinusoidOriginal, outFFTData);

    for (int i = 0; i < fftDataSize; i++) {

            Float32 hz = ((Float32)i / (Float32)fftDataSize) * NyquistMaxFreq;
            GLfloat mag = outFFTData[i];
            [fftPlot addVector2D:GLVector2DMake(hz, 0)];
            [fftPlot addVector2D:GLVector2DMake(hz, mag)];

    }

The result I get is:

enter image description here

Black lines are plotter results from FTT, positioned horizontally at their frequencies. DC value (1st black line from the left) looks OK, correctly represents y = sin(x) + 1 vertical offset.

But why 2nd black line, which represents the only frequency present in sinus equation, does not have magnitude = 1 and does not exactly stay at 1Hz?

Can anyone point me to vDSP function to convert FFT results to magnitude units from input signal?


Solution

  • Short answer:

    You seem to be using the code from my answer from this place: https://stackoverflow.com/a/19966776/468812 And it's great, but:

    You can not avoid these extra frequencies in your signal.

    The signal you generate (sin wave) is an infinite signal. And if you "crop it" and "only use a piece of the signal" it introduces "clipping noise" at both ends.

    But you can minimize the noise by taking larger input chunks and using windowing before the FFT. The Accelerate Framework provides some good and simple windowing functions. (Ex. Hann Function, vDSP_hann_window ) Also - use larger input chunks. The larger input the more precise is frequency detection.

    See this article, google: Spectral Leakage, window function.