Autocorrelations using vDSP functions

Question

Given a 1D vector of floats or doubles, how can the autocorrelations for that vector be calculated using functions from the vDSP library in the Accelerate Framework?

One would suspect the vDSP_acor() and vDSP_acorD() functions would perform this calculation, but the documentation vDSP_Library.pdf (available here) doesn't do a very good job explaining how the function arguments are used.

Similarly, the vDSP_conv() and vDSP_convD() functions mention the ability to perform correlations and convolutions between two vectors, but don't provide enough explanation or example code for me to be able to use them successfully. For example, if a filter kernel is used to convolve a 2D matrix, I would imagine two calls to vDSP_convD() would be needed, with different values of signalStride, but this is omitted from the documentation. Another omission is how the data in the filter must be packed. If padded with zeros, does it matter if the zeros come first, last, or do they need to be evenly distributed on either side of the non-zero entries? Are there requirements for filter length, result length and input length?

Suggestions for a useful examples: Implementation of autocorrelation of a vector with itself using both vDSP_acor() and vDSP_conv(). Dyadic multiplication of two arrays in the frequency domain that are packed as real data that has been forward FFT'd using vDSP_fft2d_zrip() that would be used in the calculation of an autocorrelation function before the IFT returns the un-normalized answer. Implementation of a gaussian kernel convolutionon a 1D and 2D array. Generally this is a fantastic library (can you say FAST?!), but I've found these particular functions a bit hard to understand, and the aforementioned examples would probably be widely used because they are so common in signal processing and image analysis.

Suggestions for maintainers of the vDSP_Library reference document: I assume "spatial domain" and "time domain" are equivalent throughout the document. If not, please do make that distinction. Also, please check that any formulas have well-defined parameters that match the declared names of arguments in the functions being discussed.

Footnote: here the autocorrelation I refer to is defined by: A[T] = <(X[t]-m)(X[t-T]-m)>/v, where A[T] is the autocorrelation at lag T, t is the index of the signal X, m is the average of X over all t, v is the variance of X over all t, and the angle brackets <> indicate an average over all available pairs of X's that are T apart.

Answer 1

The documentation for vDSP_acorD seems reasonably clear to me, if a little sparse, and there does appear to be a typo in the parameter descriptions.

void vDSP_acorD (double * A,
    double * C, int N, int M);

A is the input signal, C is the autocorrelation output, N is the number of samples in A, and M is the number of output values that you require in C (i.e. C will hold lag values from 0 to M - 1).

If vDSP_acorD is not available then you may be able to use vDSP_conv since convolution is the same operation as correlation with one of the input signals reversed.

Alternatively you can roll your own fast autocorrelation using the fact that autocorrelation is equivalent to the inverse FFT of the power spectrum, so:

auto_correlation = IFFT(MAG(FFT(x)))

where

FFT = forward FFT
IFFT = inverse FFT
MAG = magnitude of complex spectrum (sqrt(re * re + im * im))

and then use the FFT routines from vDSP.