iPhone FFT con Accelerate framework vDSP

Question

Tengo dificultades para implementar una FFT usando vDSP. Entiendo la teoría, pero estoy buscando un ejemplo de código específico, por favor.

que tienen datos de un archivo wav de la siguiente manera:

Pregunta 1. ¿Cómo puedo poner los datos de audio en la FFT?

Pregunta 2. ¿Cómo obtengo los datos de salida de la FFT?

Pregunta 3. El objetivo final es comprobar si hay sonidos de baja frecuencia. ¿Cómo haría esto?

-(OSStatus)open:(CFURLRef)inputURL{ 
OSStatus result = -1; 

result = AudioFileOpenURL (inputURL, kAudioFileReadPermission, 0, &mAudioFile); 
if (result == noErr) { 
    //get format info 
    UInt32 size = sizeof(mASBD); 

    result = AudioFileGetProperty(mAudioFile, kAudioFilePropertyDataFormat, &size, &mASBD); 

    UInt32 dataSize = sizeof packetCount; 
    result = AudioFileGetProperty(mAudioFile, kAudioFilePropertyAudioDataPacketCount, &dataSize, &packetCount); 
    NSLog([NSString stringWithFormat:@"File Opened, packet Count: %d", packetCount]); 

    UInt32 packetsRead = packetCount; 
    UInt32 numBytesRead = -1; 
    if (packetCount > 0) { 
     //allocate buffer 
     audioData = (SInt16*)malloc(2 *packetCount); 
     //read the packets 
     result = AudioFileReadPackets (mAudioFile, false, &numBytesRead, NULL, 0, &packetsRead, audioData); 
     NSLog([NSString stringWithFormat:@"Read %d bytes, %d packets", numBytesRead, packetsRead]); 
    } 
} 
return result; 
} 

código FFT a continuación:

log2n = N; 
n = 1 << log2n; 
stride = 1; 
nOver2 = n/2; 

printf("1D real FFT of length log2 (%d) = %d\n\n", n, log2n); 

/* Allocate memory for the input operands and check its availability, 
* use the vector version to get 16-byte alignment. */ 

A.realp = (float *) malloc(nOver2 * sizeof(float)); 
A.imagp = (float *) malloc(nOver2 * sizeof(float)); 
originalReal = (float *) malloc(n * sizeof(float)); 
obtainedReal = (float *) malloc(n * sizeof(float)); 

if (originalReal == NULL || A.realp == NULL || A.imagp == NULL) { 
printf("\nmalloc failed to allocate memory for the real FFT" 
"section of the sample.\n"); 
exit(0); 
} 

/* Generate an input signal in the real domain. */ 
for (i = 0; i < n; i++) 

    originalReal[i] = (float) (i + 1); 

/* Look at the real signal as an interleaved complex vector by 
* casting it. Then call the transformation function vDSP_ctoz to 
* get a split complex vector, which for a real signal, divides into 
* an even-odd configuration. */ 

vDSP_ctoz((COMPLEX *) originalReal, 2, &A, 1, nOver2); 

/* Set up the required memory for the FFT routines and check its 
* availability. */ 

setupReal = vDSP_create_fftsetup(log2n, FFT_RADIX2); 

if (setupReal == NULL) { 

printf("\nFFT_Setup failed to allocate enough memory for" 
"the real FFT.\n"); 

exit(0); 
} 

/* Carry out a Forward and Inverse FFT transform. */ 
vDSP_fft_zrip(setupReal, &A, stride, log2n, FFT_FORWARD); 
vDSP_fft_zrip(setupReal, &A, stride, log2n, FFT_INVERSE); 

/* Verify correctness of the results, but first scale it by 2n. */ 
scale = (float) 1.0/(2 * n); 
vDSP_vsmul(A.realp, 1, &scale, A.realp, 1, nOver2); 
vDSP_vsmul(A.imagp, 1, &scale, A.imagp, 1, nOver2); 

/* The output signal is now in a split real form. Use the function 
* vDSP_ztoc to get a split real vector. */ 
vDSP_ztoc(&A, 1, (COMPLEX *) obtainedReal, 2, nOver2); 

/* Check for accuracy by looking at the inverse transform results. */ 
Compare(originalReal, obtainedReal, n); 

Gracias

Answer 1

1. usted pone sus datos de muestras de audio en la parte real de la entrada, y cero en la parte imaginaria.

2. Si solo está interesado en la magnitud de cada contenedor en el dominio de frecuencia, entonces calcula sqrt(re*re + im*im) para cada bandeja de salida. Si solo está interesado en magnitud relativa, puede soltar el sqrt y calcular la magnitud cuadrada, (re*re + im*im).

3. Miraría las magnitudes de la papelera o los contenedores (ver (2)) que corresponden a su frecuencia o frecuencias de interés. Si su frecuencia de muestreo es Fs, y su tamaño de FFT es N, entonces la frecuencia correspondiente para la bandeja de salida i viene dada por f = i * Fs/N. Por el contrario, si está interesado en una frecuencia específica f, entonces el bin de interés, i, está dado por i = N * f/Fs.

Nota adicional: tendrá que aplicar una adecuada window function (por ejemplo Hann aka Hanning) a los datos de entrada FFT, antes de calcular la FFT en sí.

Answer 2

Puede consultar la documentación de Apple sobre eso y cuidar bien el embalaje de datos. Aquí es mi Ejemplo

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

/* 
This is an example of a hilbert transformer using 
Apple's 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); 

} 

Answer 3

Una cosa que hay que tener cuidado de que es la componente continua de la FFT calculada. Comparé mis resultados con la biblioteca fftw FFT y la parte imaginaria de la transformación calculada con la biblioteca vDSP siempre tuvo un valor diferente en el índice 0 (lo que significa 0 frecuencia, entonces DC). Otra medida que apliqué fue dividir partes reales e imaginarias por un factor de 2. Supongo que esto se debe al algoritmo utilizado en la función. Además, ambos problemas ocurrieron en el proceso de FFT pero no en el proceso de IFFT.

Utilicé vDSP_fft_zrip. Espero que esto pueda ayudar.

Paolo