Core-audio, el algoritmo Goertzel no funciona

Question

Actualmente estoy creando una aplicación que calcula la magnitud a una frecuencia predefinida (16780Hz) en tiempo real desde el micrófono del iPhone.

Tengo los datos de sonido en un búfer e intento procesarlos usando Goertzel, un algoritmo diseñado para esta tarea. Goertzel info. Aquí es donde comienza el problema.

El algoritmo responde con resultados muy positivos cuando se graba un sonido de una frecuencia mucho más baja (5000 Hz) que la definida (16780Hz). De hecho, el resultado es mucho más positivo que el producido cuando se graba un sonido de la frecuencia correcta.

Aquí es mi implementación de Goertzel:

double goertzel(unsigned short *sample, int sampleRate, double Freq, int len) 
{ 

double realW = 2.0 * cos(2.0 * M_PI * Freq/sampleRate); 
double imagW = 2.0 * sin(2.0 * M_PI * Freq/sampleRate); 
double d1 = 0; 
double d2 = 0; 
int z; 
double y; 
for (int i = 0; i < len; i++) { 
    y=(double)(signed short)sample[i] +realW * d1 - d2; 
    d2 = d1; 
    d1 = y; 
} 
double rR = 0.5 * realW *d1-d2; 
double rI = 0.5 * imagW *d1-d2; 

return (sqrt(pow(rR, 2)+pow(rI,2)))/len; 
} /* end function goertzel */ 

Aquí es cómo puedo recuperar el audio si es en absoluto relevante

-(void)startListeningWithFrequency:(float)frequency; 
{ 
OSStatus status; 
//AudioComponentInstance audioUnit; 
AudioComponentDescription desc; 
desc.componentType = kAudioUnitType_Output; 
desc.componentSubType = kAudioUnitSubType_RemoteIO; 
desc.componentFlags = 0; 
desc.componentFlagsMask = 0; 
desc.componentManufacturer = kAudioUnitManufacturer_Apple; 

AudioComponent inputComponent = AudioComponentFindNext(NULL, &desc); 
status = AudioComponentInstanceNew(inputComponent, &audioUnit); 
checkStatus(status); 

UInt32 flag = 1; 
status = AudioUnitSetProperty(audioUnit, kAudioOutputUnitProperty_EnableIO, kAudioUnitScope_Input,kInputBus, &flag, sizeof(flag)); 
checkStatus(status); 

AudioStreamBasicDescription audioFormat; 
audioFormat.mSampleRate   = 44100.00;//44100.00; 
audioFormat.mFormatID   = kAudioFormatLinearPCM; 
audioFormat.mFormatFlags  = kAudioFormatFlagIsPacked | kAudioFormatFlagIsSignedInteger; 
audioFormat.mFramesPerPacket = 1; 
audioFormat.mChannelsPerFrame = 1; 
audioFormat.mBitsPerChannel  = 16; 
// float 
audioFormat.mBytesPerPacket  = 2; 
audioFormat.mBytesPerFrame  = 2; 

status = AudioUnitSetProperty(audioUnit, 
           kAudioUnitProperty_StreamFormat, 
           kAudioUnitScope_Output, 
           kInputBus, 
           &audioFormat, 
           sizeof(audioFormat)); 
checkStatus(status); 
//status = AudioUnitSetProperty(audioUnit, 
//       kAudioUnitProperty_StreamFormat, 
//       kAudioUnitScope_Input, 
//       kOutputBus, 
//       &audioFormat, 
//       sizeof(audioFormat)); 
checkStatus(status); 
AURenderCallbackStruct callbackStruct; 
callbackStruct.inputProc = recordingCallback; 
callbackStruct.inputProcRefCon = self; 
status = AudioUnitSetProperty(audioUnit, 
           kAudioOutputUnitProperty_SetInputCallback, 
           kAudioUnitScope_Global, 
           kInputBus, &callbackStruct, sizeof(callbackStruct)); 
checkStatus(status); 
/* UInt32 shouldAllocateBuffer = 1; 
AudioUnitSetProperty(audioUnit, kAudioUnitProperty_ShouldAllocateBuffer, kAudioUnitScope_Global, 1, &shouldAllocateBuffer, sizeof(shouldAllocateBuffer)); 
*/ 
status = AudioOutputUnitStart(audioUnit); 

} 
static OSStatus recordingCallback(void *inRefCon, 
           AudioUnitRenderActionFlags *ioActionFlags, 
           const AudioTimeStamp *inTimeStamp, 
           UInt32 inBusNumber, 
           UInt32 inNumberFrames, 
           AudioBufferList *ioData) { 
AudioBuffer buffer; 

buffer.mNumberChannels = 1; 
buffer.mDataByteSize = inNumberFrames * 2; 
//NSLog(@"%d",inNumberFrames); 
buffer.mData = malloc(inNumberFrames * 2); 

// Put buffer in a AudioBufferList 
AudioBufferList bufferList; 
bufferList.mNumberBuffers = 1; 
bufferList.mBuffers[0] = buffer; 



OSStatus status; 
status = AudioUnitRender(audioUnit, 
         ioActionFlags, 
         inTimeStamp, 
         inBusNumber, 
         inNumberFrames, 
         &bufferList); 
checkStatus(status); 
//double g = calculateGoertzel((const char *)(&bufferList)->mBuffers[0].mData,16789.0,96000.0); 
UInt16 *q = (UInt16 *)(&bufferList)->mBuffers[0].mData; 
int N = sizeof(q)/sizeof(UInt16); 
double Qr,Qi; 
double theta = 2.0*M_PI*16780/44100; 
double g = goertzel(q,44100,16780,N); 

NSLog(@"goertzel:%f", g); 
} 

Esto devuelve números en los cientos de frecuencia mucho más bajos que 16780Hz , y para frecuencias de 16780Hz devuelve números mucho más pequeños.

Estoy muy frustrado y la ayuda sería muy apreciada.

Answer 1

Sólo una conjetura:

De acuerdo con el teorema de muestreo de Nyquist-Shannon, la frecuencia de muestreo debe ser al menos el doble de la frecuencia que se está tratando de medir. Y el tuyo es, pero apenas. Una velocidad de muestreo de 44.1kHz es el borde externo para medir señales de 22kHz. Una señal de 16 kHz está lo suficientemente cerca del límite que el aliasing podría causar problemas con su análisis de onda. Aquí hay una imagen para ilustrar mi punto:

Por lo tanto, supongo que necesita una mayor frecuencia de muestreo. ¿Por qué no intentas ejecutar una onda sinusoidal pura de 16kHz a través del algoritmo, para ver si mejora con eso? Aliasing será un problema menor si solo tiene una frecuencia única en los datos de prueba. Si obtienes una respuesta más alta de la onda sinusoidal, entonces probablemente solo necesites una frecuencia de muestreo más alta.

Answer 2

Parece que el resonador utilizado en su filtro Goertzel es una aproximación de 1er grado al resonador de 1 polo. Esto disminuirá en gran medida en la precisión y la estabilidad en ángulos de fase altos por paso. Un DFT de 1 compartimiento que use una mejor aproximación a las funciones trigonométricas podría funcionar mejor a tales altas frecuencias.

Y la respuesta de frecuencia del micrófono de iPhone probablemente se dispara a frecuencias tan altas.

añadido:

Para una DFT de 1 separador, a probar esto en su bucle interno:

d1 += (double)sample[i] * cos(2.0*M_PI*i*Freq/sampleRate); 
d2 += (double)sample[i] * sin(2.0*M_PI*i*Freq/sampleRate); 

A continuación, regrese:

dR = d1; 
dI = d2; 
magnitude = sqrt(dR*dR + dI*dI)/(double)len; 

Tenga en cuenta que para una frecuencia fija y frecuencia de muestreo las funciones trigonométricas se pueden calcular previamente fuera de la devolución de llamada de audio y se pueden guardar en una matriz o tabla de búsqueda. Si no hace una optimización como esta, llamar a múltiples funciones trascendentales de doble precisión dentro de su devolución de llamada de audio puede ser demasiado lenta y/o perder mucha energía de la batería, pero puede simular OK en una PC rápida típica.

La DFT se define para una longitud que es un número integral exacto de periodos de la frecuencia Frec Freq, pero otras longitudes funcionarán para aproximaciones que contienen cantidades variables de los llamados errores espectrales de "fuga" y/o festoneado. El ancho de la respuesta de frecuencia del filtro será aproximadamente inversamente proporcional a la longitud DFT. Además, cuanto más cerca esté la frecuencia de Fs/2, más tiempo necesitará el DFT para evitar alias de imágenes complejas, tal vez múltiples períodos de longitud N * Fs/(Fs/2 - Freq) sean de mejor longitud. Puede que necesite guardar o poner en cola las muestras para obtener la longitud adecuada (y no solo usar la duración del búfer que le proporcionó la devolución de llamada de audio).