Fast Convolution with vDSP/Accelerate Framework

[hugoderwolf]

I'm a bit stuck with doing fast convolution with Apple's vDSP Framework. I got everything quite running, using an overlap-save algorithm. But the actual multiplication of two FFTed buffers just won't work.

Here's a code excerpt:

Code: Select all

	vDSP_ctoz((COMPLEX *)inBuffer, 2, &sigFFT, 1, bufferSize / 2);
	
	vDSP_fft_zrip(fftSetup, &sigFFT, 1, log2n, FFT_FORWARD);
	vDSP_zvmul (&sigFFT, 1, &irFFT, 1, &sigFFT, 1, bufferSize, 1); // this goes wrong
	vDSP_fft_zrip(fftSetup, &sigFFT, 1, log2n, FFT_INVERSE);
	
	vDSP_ztoc(&sigFFT, 1, (COMPLEX *)resultBuffer, 2, bufferSize / 2);
	float scale = 1.0 / (2*bufferSize);
	vDSP_vsmul(resultBuffer, 1, &scale, resultBuffer, 1, bufferSize);

irFFT holds FFT Data of my Impulse response and is computed during initialization. So what this code does is convert the input buffer to the right form needed for FFT computation, compute FFT, complex multiply with FFTed impulse response, inverse FFT the result, convert it to a linear buffer and scale it correctly.

I'm pretty sure the problem is in this zvmul() operation, because I get the orgiginal inBuffer back when commenting the line out. I also ensured that irFFT ist OK by inverse transforming and verifying the result. Anybody have a clue? What is it I'm not seeing?

[Snapey]

Yes there is a bug in zvmul, so you will have to write a routine that does the split complex multiply manually. There is some twiddling involved when dealing with DSPSplitComplex, so you will need to refer to the vDSP programming guide. Even though you can't use zvmul, you should still be able to vectorise with zvmul etc.

I'm not getting the correct output values either. I've checked the fft/ ifft routines they seem to work OK. I think I have got the problem down to the scaling factor. Putting it as 1/ 4 * bufferSize yields a closer result but its still not right!

I'll let you know if I make any progress!

L

[hugoderwolf]

This was exactly one year ago, and the solution to my problem is nearly as old.

The Problem is the way, vDSP stores the results of real FFTs. For real signals, the DC and Nyquist bins are always real, so to save memory and for convenience reasons, vDSP stores the Nyquist bin as the imaginary part of the DC bin, so there's indeed some twiddling involved when multiplying the two spectra.

[Snapey]

I take it you got it working in the end. Do you have any example code that works? I cannot for the life of me work out what is going wrong!

Cheers,
L

[hugoderwolf]

Here's the snippet that's working here:

Code: Select all

	vDSP_ctoz((DSPComplex *)inBuffer, 2, sigFFT, 1, N);

	vDSP_fft_zrip(*(FFTSetup*)fftSetup, sigFFT, 1, log2n, FFT_FORWARD);

	float preserveIRNyq = irFFT->imagp[0];
	irFFT->imagp[0] = 0;
	float preserveSigNyq = sigFFT->imagp[0];
	sigFFT->imagp[0] = 0;
	vDSP_zvmul(sigFFT, 1, irFFT, 1, sigFFT, 1, N, 1);
	sigFFT->imagp[0] = preserveIRNyq * preserveSigNyq;
	irFFT->imagp[0] = preserveIRNyq;
	vDSP_fft_zrip(*(FFTSetup*)fftSetup, sigFFT, 1, log2n, FFT_INVERSE);

	vDSP_ztoc(sigFFT, 1, (DSPComplex *)resultBuffer, 2, N);

	float scale = 1.0 / (8*N);
	vDSP_vsmul(resultBuffer, 1, &scale, resultBuffer, 1, 2*N);

[Snapey]

thanks!

[to-pse]

I was wondering if anyone here had a pointer to some
sample-code on how to vectorize filter-algorithms to
use them with the vDSP/Accelerate framwork.

So far I have only seen fft-usage...

Best regards,
Tobias

[Snapey]

Are you referring to calculating the filter coefficients on the fly for FIR filters?

[to-pse]

Are you referring to calculating the filter coefficients on the fly for FIR filters?

No - I meant a vectorized version of the actual FIR or IIR implementation to
reduce load instead of using scalar fpu code...

Tobias

[hugoderwolf]

You'd do long FIR filtering (more than about 30 taps) using fast FFT convolution, which is what is done in the above code. If you do it straightforward however, you will have a latency equal to the FFT length.

There are some hybrid methods like this one without latency:
http://www.cs.ust.hk/mjg_lib/bibs/DPSu/ ... s/Ga95.PDF

You can't parralelize IIR filters afaik due to their feedback structure.

[mjdav]

Nice nifty implementation!