CoreAudio code to sample the first channel of the OS X mic
As in the book example, I used the first mic found on the system.
// generate description that will match audio HAL
AudioComponentDescription inputcd = {0};
inputcd.componentType = kAudioUnitType_Output;
inputcd.componentSubType = kAudioUnitSubType_HALOutput;
inputcd.componentManufacturer = kAudioUnitManufacturer_Apple;
AudioComponent comp = AudioComponentFindNext(NULL, &inputcd);
CheckError(AudioComponentInstanceNew(comp, &player->inputUnit),
"Couldn't open component for inputUnit");
propertySize = sizeof (AudioStreamBasicDescription);
CheckError(AudioUnitGetProperty(player->inputUnit,
kAudioUnitProperty_StreamFormat,
kAudioUnitScope_Output, //to other units
inputBus,
&player->streamFormat,
&propertySize),
"Couldn't get ASBD from output unit");
//Filter designed for this rate, so fix at expected frequency
player->streamFormat.mSampleRate = 44100;
CheckError(AudioUnitSetProperty(player->inputUnit,
kAudioUnitProperty_StreamFormat,
kAudioUnitScope_Output,
inputBus,
&player->streamFormat,
propertySize),
"Couldn't set output unit Fs");
AudioStreamBasicDescription deviceFormat;
CheckError(AudioUnitGetProperty(player->inputUnit,
kAudioUnitProperty_StreamFormat,
kAudioUnitScope_Input, //from HW to IO unit
inputBus,
&deviceFormat,
&propertySize),
"Couldn't get ASBD from input unit");
deviceFormat.mSampleRate = 44100;
CheckError(AudioUnitSetProperty(player->inputUnit,
kAudioUnitProperty_StreamFormat,
kAudioUnitScope_Input,
inputBus,
&deviceFormat,
propertySize),
"Couldn't set Fs input unit");
UInt32 bufferSizeFrames = 1<<LOG2M;
propertySize = sizeof(UInt32);
CheckError(AudioUnitSetProperty(player->inputUnit,
kAudioDevicePropertyBufferFrameSize,
kAudioUnitScope_Global,
0,
&bufferSizeFrames,
propertySize),
"Couldn't set buffer frame size to input unit");
Chapter 08 example uses AUGraph to connect the mic HAL unit to default output unit (speaker), but I just want the raw signal from the mic, so the program is significantly simplified: I just have to start the AUHAL unit explicitly:
CheckError(AudioUnitInitialize(player->inputUnit),
"Couldn't initialize input unit");
// Input unit is now ready for use
CheckError(AudioOutputUnitStart(player->inputUnit)
, "AudioOutputUnitStart");
Later, when I want to stop receiving the mic input, I call AudioOutputUnitStop instead. But in the continuous chirp detection scenario, I won't stop the input. The CoreAudio interface for sucking out the sound samples is kAudioOutputUnitProperty_SetInputCallback.
AURenderCallbackStruct callbackStruct;
callbackStruct.inputProc = InputRenderProc;
callbackStruct.inputProcRefCon = player;
CheckError(AudioUnitSetProperty(player->inputUnit,
kAudioOutputUnitProperty_SetInputCallback,
kAudioUnitScope_Global,
0,
&callbackStruct,
sizeof(callbackStruct)),
"Couldn't set input callback");
The callback routine will be called (semi) real-time and therefore is time critical; similar to interrupt handlers, the callback routine must return as quickly as possible. This means that I need to keep pushing the (zero padded) received signal from the audio reception callback to the signal processing algorithm running in another thread.
Memory pool and lockless queue
In the CoreAudio book, data is shuttled from the callback to the downstream through the CARingBuffer, which is specialized for the CA input buffer (this is what CoreAudio returns to me when query the HW for new sample block). Instead of the CARingBuffer the book uses, I wrote my own memory pool and queue for 2 reasons:- I need to send meta data (like the host timestamp for the input block) along with the raw samples to the algorithm.
- I am only interested in 1 channel (the first one) of the potentially stereo mic, so copying the both channels through the CARingBuffer would be a waster of memory.
The interface to the (aligned) memory pool was copied from another blog I wrote more than 5 years ago:
typedef struct llsMP {
size_t _capacity, _available;
void* _pool;/* the memory itself */
void* *_book;/* Array for book keeping, but don't know the size till ctor */
} llsMP;
unsigned char llsMP_alloc(struct llsMP* me,
size_t capacity, size_t memsize, size_t alignment
, unsigned char zero);
void llsMP_free(struct llsMP* me);
struct llsMP* llsMP_new(size_t capacity, size_t memsize, size_t alignment
, unsigned char zero);
void llsMP_delete(struct llsMP* me);
void* llsMP_get(struct llsMP* me);
unsigned char llsMP_return(struct llsMP* me, void* node);
In this primitive, cross-platform API, I don't use more convenient types like uint8_t or bool, because I did not want to require stdint.h. I prefer using the non-memory alloc version if possible, to keep the object as part of the container's object, like this:
typedef struct MyMic {
AudioStreamBasicDescription streamFormat;
AudioUnit inputUnit;
AudioBufferList *inputBuffer;
FILE* t_csv, *x_csv, *f_csv, *c_csv;
struct llsMP mpool;
struct llsQ padded_xQ;
int8_t Eblock, iBlock;
bool firstBlock;
} MyMic;
if (!llsMP_alloc(&player->mpool
, 3 // capacity
, sizeof(struct SampleBlock) // memsize
, 16 // alignment: recommended by vDSP programming guide
, 1)//memset zero, since this pool is used for zero padding
) {
fprintf (stderr, "Can't allocate zero padded memory pool");
exit (errno);
}
struct SampleBlock {
UInt64 mHostTime;//The timestamp for these frames
//UInt64 mWordClockTime //don't bother; always zero
//I can calculate the time of the 1st sample with this just as well
uint32_t iBlock;
UInt32 nFrames;
float sample[1<<LOG2_2M];//2x sample block length, to zero pad
};
struct SampleBlock* block = (struct SampleBlock*)llsMP_get(&player->mpool);
if (!block) {
fprintf(stderr, "Memory pool exhausted\n");
return errno;
}
player->inputBuffer->mBuffers[0].mData = (void*)(&block->sample[1<<LOG2M]);
CheckError(AudioUnitRender(player->inputUnit, ioActionFlags,
inTimeStamp, inBusNumber, inNumberFrames,
player->inputBuffer),
"AudioUnitRender failed");
block->iBlock = player->iBlock;
block->mHostTime = inTimeStamp->mHostTime;
block->nFrames = inNumberFrames;
if (!llsQ_push(&player->padded_xQ, block)) {
fprintf(stderr, "llsQ_push failed\n");
return errno;
}
As you can see, the pointer to the memory pool node then can be pushed into a lockless queue, which works well in this case because there is only 1 writer thread, and 1 reader thread. Note also that the reception time on the host (in nanoseconds) and the block number (starts at 0) are attached to the sample block. If the correlation processor decides there is a positive event, it can then place the event within the audio sampling period (1/44100 Hz here, or < 23 usec).
/* SINGLE writer, SINGLE reader queue */
typedef struct llsQ {
size_t _head, _tail, _mask;
void* *_q;/* array of pointers */
} llsQ;
unsigned char llsQ_alloc(struct llsQ* me, unsigned char exponent);
void llsQ_free(struct llsQ* me);
struct llsQ* llsQ_new(unsigned char exponent);
void llsQ_delete(struct llsQ* me);
unsigned char llsQ_push(struct llsQ* me, void* node);
unsigned char llsQ_pop(struct llsQ* me, void** node);
The main algorithm loop sleeps until the lockless queue becomes not-empty.
for (int8_t iBlock = 0; iBlock < player->Eblock; ) {
struct SampleBlock* block;
if (!llsQ_pop(&player->padded_xQ, (void**)&block)) {
usleep(1000);
continue;
}
// Have a new block to process
The correlation processor now has work to do.
Fast correlation processor
I discussed the central algorithm for the fast correlation using FFT in my previous blog entry. Repeating for convenience:
correlation(x, filter) = iFFT(FFT(x) • FFT*(filter))
where • is the element wise product of 2 vectors of equal length, and * is the complex conjugate operator. Because the filter coefficients are constant, I can pre-calculate FFT*(filter) and store in the DATA segment of the program, as I've done here:
#define LOG2M 10
#define EXP_LOG2M (1<<LOG2M)
#define LOG2_HALFM (LOG2M - 1)
#define LOG2_2M (LOG2M + 1)
// The matched filter from Matlab
float FFT_pulse_conj_vDSPz_flattened[1<<LOG2_2M] __attribute__((aligned(16)))
= {
#include "FFT_pulse_conj_vDSPz_flattened.h"
};
Given the wealth of DSP libraries and algorithms available for free these days, I expected to find a DSP library that is optimized for OS X/iOS, and I did: Apple's Accelerate framework, which seems to take advantage of Intel's SSE. But typical of Apple documentation, the documentation is lacking (to put it mildly). One must begin with the vDSP's packing of the FFT of a real sequence, given here. A few crucial concepts are:
- The result of FFT is a complex number in general, which they keep in 2 DIFFERENT arrays (lots of other implementations, including the venerable Numerical Recipes) keeps them in alternating fashion.
- FFT of a real sequence x is anti-symmetrical about frequency=0; i.e. X(f) = X*(-f), where X = FFT(x).
- X(f=0) and X(f=Fs/2)--FFT at the DC and the Nyquist frequency--are pure real.
Using all of the above properties, vDSP uses N numbers to represent the FFT of a real sequence of length N (which would otherwise require 2N numbers). But because they keep the real and the imaginary parts separately, an in-place FFT requires a prior step of separating out the real sequence into odd and even sequences, as I've done below with vDSP_ctoz.
float creal[1<<LOG2M] __attribute__((aligned(16))) //input/output buffer
, cimag[1<<LOG2M] __attribute__((aligned(16)))
, ftemp[1<<LOG2_2M] __attribute__((aligned(16))) //temporary scratch pad
;
COMPLEX_SPLIT splitc = { .realp = creal, .imagp = cimag }
, split_temp = { .realp = ftemp, .imagp = &ftemp[1<<LOG2M] }
, split_filter = { .realp = FFT_pulse_conj_vDSPz_flattened
, .imagp = FFT_pulse_conj_vDSPz_flattened + (1<<LOG2M) }
;
float FFT_filter_nyq = *split_filter.imagp;
*split_filter.imagp = 0; //see explanation below
...
// When I get a block to process, split out a time series into odd and even
vDSP_ctoz((COMPLEX*)block->sample, 2, &splitc, 1, 1<<LOG2M);
Note the arrangement of the filter's FFT to be consistent with vDSP's representation of FFT of a real sequence. My Matlab code (explained in my previous blog post) anticipates the C code will do this, and writes out the header file correctly. I can then FFT the zero padded sample block with vDSP API:
vDSP_fft_zript(fftSetup, &splitc, 1, &split_temp, LOG2_2M, FFT_FORWARD);
The result is still in the "splitc"--because this is an inplace FFT. The temporary buffer is supposed to speed up the FFT (have not yet measured how much). If I multiply this FFT with the filter's FFT (actually the conjugate thereof), the result would be wrong, because the Nyquist portion masquerading as the imaginary value of FFT(x) @ DC distorts the result. Because vDSP does not take care of this, I have to jump through some hoops to save and then restore. At this this is straight-forward in the 1D case; the 2D case requires setting up a separate vector (huge inconvenience). I think the vDSP folks dropped the ball on this.
float FFT_c_nyq = *splitc.imagp; *splitc.imagp = 0;
vDSP_zvmul(&splitc, 1, &split_filter, 1, &splitc, 1
, 1<<LOG2M//The multiplication is over vector of length M
, 1 // don't conjugate (even though we are correlating);
); // because the filter is already in conjugate form
// Restore the Nyquist frequency portion, which is pure real
*splitc.imagp = FFT_c_nyq * FFT_filter_nyq;
I can then inverse FFT to obtain the fast correlation..
vDSP_fft_zript(fftSetup, &splitc, 1, &split_temp, LOG2_2M, FFT_INVERSE);
vDSP_ztoc(&splitc, 1, (COMPLEX*)ftemp, 2, 1<<LOG2M);
But as I explained in previous blog post, the correlation of the zero-padded block is not the final running correlation unless the result from the previous block is added to the overlap (the overlap save method). That is why the overlap_save (initially zero) is maintained from 1 block to the next:
float overlap_save[1<<LOG2M] __attribute__((aligned(16))) //final correlation
;
...
//Correlation to output: overlap_save + ftemp[0:M-1];
vDSP_vadd(ftemp, 1, overlap_save, 1, ftemp, 1, 1<<LOG2M);
//Correlation to save for next round: overlap_save = ftemp[M:2M-1]
memcpy(overlap_save, &ftemp[1<<LOG2M], sizeof(overlap_save));
Result: clear difference between silence and chirp
To validate the output of the fast correlator, I ran the program twice: the first without any signal (just room noise) and the 2nd while playing the chirp I designed in Matlab on repeat. The Matlab validation script's orange line shows the difference between Matlab's calculation and the output from my CoreAudio program. The 1st panel is the raw sample value, and the middle panel is the magnitude of the element-wise multiplication of the 2 FFTs. The bottom panel is the correlator output.
I am not sure what the mic is picking up during the first few milliseconds of the recording, but there is no obvious hit in the correlator output, which is clearly visible when I play the chirp, as you can see below.
Next steps
When I can further process the correlator output to a yes/no decision between chirp, picking out the sample at which the chirp was received, and mapping that back to a monotonically increasing timestamp on the receiver side should be straight-forward. And thanks to the use of the memory pool and circular queue, this processing pipeline can run ad-infinitum.
