IPhone FFT с ускорением рамки vDSP

Ответ 1

• Вы помещаете свои данные выборок аудио в реальную часть ввода и нулевую мнимую часть.

• Если вас просто интересует величина каждого бункера в частотной области, вы вычисляете sqrt(re*re + im*im) для каждого выходного лотка. Если вас интересует только относительная величина, вы можете отказаться от sqrt и просто вычислить квадрат величины, (re*re + im*im).

• Вы бы посмотрели величины бина или бункера (см. (2)), которые соответствуют вашей частоте или частотам, представляющим интерес. Если ваша частота дискретизации равна Fs, а ваш размер FFT равен N, то соответствующая частота для выходного лотка i задается f = i * Fs / N. И наоборот, если вас интересует определенная частота f, то интересующий нас бункер i задается i = N * f / Fs.

Дополнительное примечание: вам нужно будет применить подходящую функцию окна (например, Hann aka Hanning) к вашим входным данным FFT, до вычисления самого FFT.

Ответ 2

Вы можете проверить документацию Apple и позаботиться об упаковке данных.

Вот мой пример:

//  main.cpp
//  FFTTest
//
//  Created by Harry-Chris Stamatopoulos on 11/23/12.
//  

/* 
 This is an example of a hilbert transformer using 
 Apple VDSP fft/ifft & other VDSP calls.
 Output signal has a PI/2 phase shift.
 COMPLEX_SPLIT vector "B" was used to cross-check
 real and imaginary parts coherence with the original vector "A"
 that is obtained straight from the fft.
 Tested and working. 
 Cheers!
*/

#include <iostream>
#include <Accelerate/Accelerate.h>
#define PI 3.14159265
#define DEBUG_PRINT 1

int main(int argc, const char * argv[])
{
    float fs = 44100;           //sample rate
    float f0 = 440;             //sine frequency
    uint32_t i = 0;

    uint32_t L = 1024;

    /* vector allocations*/
    float *input = new float [L];
    float *output = new float[L];
    float *mag = new float[L/2];
    float *phase = new float[L/2];

    for (i = 0 ; i < L; i++)
    {
        input[i] = cos(2*PI*f0*i/fs);
    }

    uint32_t log2n = log2f((float)L);
    uint32_t n = 1 << log2n;
    //printf("FFT LENGTH = %lu\n", n);

    FFTSetup fftSetup;
    COMPLEX_SPLIT A;
    COMPLEX_SPLIT B;
    A.realp = (float*) malloc(sizeof(float) * L/2);
    A.imagp = (float*) malloc(sizeof(float) * L/2);

    B.realp = (float*) malloc(sizeof(float) * L/2);
    B.imagp = (float*) malloc(sizeof(float) * L/2);

    fftSetup = vDSP_create_fftsetup(log2n, FFT_RADIX2);

    /* Carry out a Forward and Inverse FFT transform. */
    vDSP_ctoz((COMPLEX *) input, 2, &A, 1, L/2);
    vDSP_fft_zrip(fftSetup, &A, 1, log2n, FFT_FORWARD);

    mag[0] = sqrtf(A.realp[0]*A.realp[0]);

    //get phase
    vDSP_zvphas (&A, 1, phase, 1, L/2);
    phase[0] = 0;

    //get magnitude;
    for(i = 1; i < L/2; i++){
        mag[i] = sqrtf(A.realp[i]*A.realp[i] + A.imagp[i] * A.imagp[i]);
    }

    //after done with possible phase and mag processing re-pack the vectors in VDSP format
    B.realp[0] = mag[0];
    B.imagp[0] = mag[L/2 - 1];;

    //unwrap, process & re-wrap phase
    for(i = 1; i < L/2; i++){
        phase[i] -= 2*PI*i * fs/L;
        phase[i] -= PI / 2 ;
        phase[i] += 2*PI*i * fs/L;
    }

    //construct real & imaginary part of the output packed vector (input to ifft)
    for(i = 1; i < L/2; i++){
        B.realp[i] = mag[i] * cosf(phase[i]);
        B.imagp[i] = mag[i] * sinf(phase[i]);
    }

#if DEBUG_PRINT
    for (i = 0 ; i < L/2; i++)
    {
       printf("A REAL = %f \t A IMAG = %f \n", A.realp[i], A.imagp[i]);
       printf("B REAL = %f \t B IMAG = %f \n", B.realp[i], B.imagp[i]);
    }
#endif
    //ifft
    vDSP_fft_zrip(fftSetup, &B, 1, log2n, FFT_INVERSE);

    //scale factor
    float scale = (float) 1.0 / (2*L);

    //scale values
    vDSP_vsmul(B.realp, 1, &scale, B.realp, 1, L/2);
    vDSP_vsmul(B.imagp, 1, &scale, B.imagp, 1, L/2);

    //unpack B to real interleaved output
    vDSP_ztoc(&B, 1, (COMPLEX *) output, 2, L/2);

    // print output signal values to console
    printf("Shifted signal x = \n");
    for (i = 0 ; i < L/2; i++)
        printf("%f\n", output[i]);

    //release resources
    free(input);
    free(output);
    free(A.realp);
    free(A.imagp);
    free(B.imagp);
    free(B.realp);
    free(mag);
    free(phase);
}

Ответ 3

Одна вещь, с которой вам нужно быть осторожным, это постоянная составляющая рассчитанного БПФ. Я сравнил свои результаты с FFT библиотеки fftw, и воображаемая часть преобразования, вычисленная с помощью библиотеки vDSP, всегда имела другое значение с индексом 0 (что означает частоту 0, то есть DC). Другая мера, которую я применил, состояла в том, чтобы разделить как действительные, так и мнимые части в 2 раза. Я думаю, это связано с алгоритмом, используемым в функции. Кроме того, обе эти проблемы возникли в процессе FFT, но не в процессе IFFT.

Я использовал vDSP_fft_zrip.