Autocorrelation functions, was Re: [PD-dev] using SIMD instructions in externals

Mathieu Bouchard matju at artengine.ca
Fri Jan 20 00:20:45 CET 2006


On Tue, 17 Jan 2006, Ed Kelly wrote:

> You need to use some kind of FFT in order to benefit from the fantastic
> speed increase possible with the Convolution Theorem. Right. A*B=
> fourier AxB, or cepstral A+B. I know a bit about Fourier theory, real
> and imaginary and radix-2 etc optimisations although I have never
> implemented one myself (no need).

I'd rather not use the letter "x" to mean multiplication and not use "*" 
to mean convolution. I'd use "*" and "conv" respectively.

The radix-2 optimisations of FFT are generalizable to any radix. When it
is said that FFT runs in O(n log n) time, what's really meant by log n is
more like, the sum of all prime factors of n. I'm currently thinking of
non-power-of-2 FFT's because I'd like to apply it in the spatial domain
for images that don't have power-of-2 sizes, without having to crop or
scale each picture. I'm not sure how much it would be worth it.

> What I am wondering is how to turn any form of autocorrelation matrix
> into a radix-type algorithm, so that code such as the two examples
> enclosed can be optimized.

>   temp0 += i == 0 ? 0.0 : fabs(in[j] - temp[j]);

Note that the value of i can't be 0, so it's simplifiable to:

>   temp0 += fabs(in[j] - temp[j]);

Then what's the purpose of the temp array? if each value is set only once 
and read immediately, you don't need it.

In any case, AMDF doesn't seem to be FFT-optimisable. Especially, the
absolute value is difficult to get rid of. That doesn't mean that there
isn't a O(n log n) trick to compute it though.

I don't understand your AMDF algorithm. Why do you overwrite temp0 in the 
inner loop, discarding all previous accumulations of fabs ???

In the other algorithm, the autocorrelation section is highly optimisable.
If I'm not too confused, correlation is like convolution with with a
mirror image (in the time domain) of its right-hand function. Note that
the double-fft of a signalblock is the mirror image of the signalblock,
and that the triple-fft is the same as the inverse-fft (supposing a
normalized fft: Pd's fft is not, so extra constants have to be introduced
in order to compensate). Thus the autocorrelation of x(t) and y(t) is
ifft(fft(x)*ifft(y)) where * is pointwise complex multiplication.

However I don't understand the "predictor" portion of the algorithm and I
know that it is running in O(n*n) time. This doesn't mean that optimising
the autocorrelation is useless: you could still get a (maybe) 50% speedup 
of the whole.

>  ...maybe I should take another degree? if only!

A degree in what exactly???... I'd rather use math books than math profs.

 _ _ __ ___ _____ ________ _____________ _____________________ ...
| Mathieu Bouchard - tél:+1.514.383.3801 - http://artengine.ca/matju
| Freelance Digital Arts Engineer, Montréal QC Canada




More information about the Pd-dev mailing list