Core-audio, Goertzel algorithm not working

Question

I am currently creating an application which works out the magnitude at a predefined frequency (16780Hz) in realtime from the iPhone's microphone.

I have the sound data in a buffer and I attempt to process it using Goertzel, an algorithm designed for this task. Goertzel info. This is where the problem begins.

The algorithm responds with very positive results when a sound is recorded of a frequency which is much lower (5000Hz) than the defined one (16780Hz). In fact the result is far more positive than that produced when a sound of the correct frequency is recorded.

Here is my implementation of 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 */

Here is how I retrieve the audio if it is at all relevant

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

This returns numbers in the hundreds for frequency much lower than 16780Hz, and for frequencies of 16780Hz returns much smaller numbers.

I am very frustrated and help would be greatly appreciated.

Answer 1

Just a guess:

According to the Nyquist–Shannon sampling theorem, the sampling rate should be at least twice the frequency that you are trying to measure. And yours is, but just barely. A sampling rate of 44.1kHz is the outer edge for measuring 22kHz signals. A signal of 16kHz is close enough to the limit that aliasing might cause problems with your wave analysis. Here's a picture to illustrate my point:

So, I would guess that you need a higher sample rate. Why don't you try running a pure 16kHz sine wave through the algorithm, to see if it does better with that? Aliasing will be less of an issue if you only have a single frequency in the test data. If you get a higher response from the sine wave, then you probably just need a higher sampling rate.

Answer 2

It looks like the resonator used in your Goertzel filter is a 1st degree approximation to 1-pole resonator. This will greatly decrease in accuracy and stability at high phase angles per step. A 1-bin DFT using a better approximation to the trig functions might work better at such high frequencies.

And the iPhone microphone frequency response likely rolls off at such high frequencies.

ADDED:

For a 1-bin DFT, try this in your inner loop:

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

Then return:

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

Note that for a fixed frequency and sample rate the trig functions can be pre-calculated outside of the audio callback and saved in an array or lookup table. If you don't do some optimization like this, calling multiple double precision transcendental functions inside your audio callback may be way too slow and/or waste a lot of battery power, but may simulate OK on a typical fast PC.

The DFT is defined for a length which is an exact integral number of periods of the bin frequency Freq, but other lengths will work for approximations containing varying amounts of so-called spectral "leakage" and/or scalloping errors. The width of the filter frequency response will be roughly inversely proportional to the DFT length. Also the closer the frequency is to Fs/2, the longer the DFT will need to be to avoid complex image aliasing, maybe multiple periods of length N*Fs/(Fs/2 - Freq) would be a better length. You may need to save or queue up samples to get an appropriate length (and not just use the buffer length given you by the audio callback).