Two FFTs
=========================================================
 DSP Trick Submission Form 
=========================================================
THIS WORK IS PLACED IN THE PUBLIC DOMAIN
Name: The Two FFTs Trick
Category: Algorithmic Trick.
Application: I used this this trick to make fast convolution even
faster.
Introduction: Use this trick to perform two simultaneous realvalued
FFTs using only a single complexvalued FFT.
The Trick:
You can take two realvalued FFTs at the same time by making use of
the imaginary part of the transform.
Suppose we have the timeindexed series
X[n] = a[n] + jb[n]
The transform is (r is real and i is imaginary)
X[k] = Xr[k] + jXi[j]
By splitting both the real and imaginary parts into even and odd (e
and o) parts, you can retrieve the transforms of a and b.
Xa[k] = Xre[k] + jXio[k]
Xb[k] = (Xro[k] + jXie[k])/j
I divided Xb by j because it's the imaginary part, but we want it to
be real. If you put a realvalued sequence, a[n], in the real part of
x[n] and a realvalued sequence, b[n], in the the real part of x[n],
you'll get the transforms in Xa[k] and Xb[k].
An example of circular convolution using this trick in matlab is
below. It's very easy to see what's going on using this high level
syntax. Remember that the even and odd parts are
Xe[n] = (X[n] + X[n])/2
Xo[n] = (X[n]  X[n])/2
The n is simply an axis reversal. You won't see any divide by twos.
That's because I've multiplied them and moved the result outside the
final inverse transform. The reason for this is that division is
slow. If you remember that the definition of the Inverse DFT has a
division as the last step, this is easy to understand. You can merge
the divide by four into this division. Then, with a little effort, you
can get two fast convolutions for the cost of two FFTs, one IFFT, and
some extra additions.
a=rand(1,6);b=rand(1,6);c=rand(1,6);d=rand(1,6);
ifft(fft(a).*fft(c))
ifft(fft(b).*fft(d))
xtmp1=fft(a+b*j);
xtmp2=fft(c+d*j);
x1=xtmp1(2:end);
x2=xtmp2(2:end);
re1=real(x1+x1(end:1:1))+imag(x1x1(end:1:1))*j;
im1=real(x1x1(end:1:1))/(j)+imag(x1+x1(end:1:1));
re2=real(x2+x2(end:1:1))+imag(x2x2(end:1:1))*j;
im2=real(x2x2(end:1:1))/(j)+imag(x2+x2(end:1:1));
y1=[real(xtmp1(1))*4 re1].*[real(xtmp2(1)) re2];
y2=[imag(xtmp1(1))*4 im1].*[imag(xtmp2(1)) im2];
ifft(y1+y2*j)/4
