Convolution using real-valued FFT

I have implemented the basic overlap-add convolution using a standard complex FFT routine.
This has lots of zeroes, first from padding data and IR (impulse response) to FFT size,
then I have to copy data and IR to complex numbers (with zero imaginary part).

Now I tried instead to replace the complex FFT routine with a real-valued FFT
so I can skip copying all data and IR to complex numbers.
The real-FFT is Sorensen's old code, but now the final convolved output gets all screwed up.

The output from Sorensen's real-FFT is said to be:
[ re(0), re(1), ... re(N/2), im(N/2-1), im(N/2-2), ... im(1) ]

which I don't quite understand.

Anyway, I can't figure out how to do the convolution.
I try the regular complex number mult, using the addressing given above.

But I'm misunderstanding the basic data format - output gets shuffled around.
It's not noise, but like echoes front and back.

So, can anyone please point me to a place that explains how to convolve
IR and data in real-value FFT's? (or just tell me ;-)

I haven't found this specific case when trawling all the usual papers and web-articles
on convolution.


Btw, I have successfully replaced the complex FFT with a Hartley transform, which works on real values only.
This lets me skip the complex numbers.
Of course I have to do the right Hartley multiplication for the convolution of data and IR,
but that was easy to get from some article. (It's fairly similar to the Sorensen FT)

So I do already have a real-value algorithm that works, but I'd like to be able
to switch around to get the best FFT, since there are a lot more of those than fHT's.

Thanks for any tips,

	Joakim

This reference looks pretty good and describes the case where you are using a 1/2 length fft and putting real data in both the real and imaginary buffers. I assume this is what you are talking about.



http://www.mpi-hd.mpg.de/astrophysik/HEA/internal/Numerical_Recipes/f12-3.pdf

Equation 12.3.5 shows how to convert the fft output back to the conventional full-length fft output in normal order.


Bob

Thanks for your response.
The NR comments alerted me to check the order of inputs;
that was probably the problem.

Meanwhile I had found another FFT method that takes straight input and returns straight output.
That one worked right for me.
Unfortunately it runs slowest of all I tested, so need to fiddle some more...

Still, most FFT's I've looked at are implemented like

FFT(data, n) {
   hartley(data, n)
   shuffle(data, n)
}

IFFT(data, n) {
   unshuffle(data, n)
   hartley(data, n)
}

- it seems I can extract the hartley from the best FFT and use that for convolution.
The extra iterations over the full data/FT array is useless in this case.

The Hartley convolution means accessing the arrays like
y(0) ops ir(0)
y(1), ir(1) ops y(N-1), ir(N-1)
y(2), ir(2) ops y(N-2), ir(N-2)
etc
(where "ops" are some muls and adds)
I hope the backward iteration is not too cache-thrashing?

NB: haven't tested FFTW or djbfft, I want single file to drop into project.

Cheers,
	Joakim


On 2013-03-14 12:31, radams2000@gmail.com wrote:
> This reference looks pretty good and describes the case where you are using a 1/2 length fft and putting real data in both the real and imaginary buffers. I assume this is what you are talking about.
>
>
>
> http://www.mpi-hd.mpg.de/astrophysik/HEA/internal/Numerical_Recipes/f12-3.pdf
>
> Equation 12.3.5 shows how to convert the fft output back to the conventional full-length fft output in normal order.
>
>
> Bob