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 real-valued
FFTs using only a single complex-valued FFT.

The Trick:

You can take two real-valued FFTs at the same time by making use of
the imaginary part of the transform.

Suppose we have the time-indexed 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 real-valued sequence, a[n], in the real part of
x[n] and a real-valued 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(x1-x1(end:-1:1))*j;
im1=real(x1-x1(end:-1:1))/(j)+imag(x1+x1(end:-1:1));
re2=real(x2+x2(end:-1:1))+imag(x2-x2(end:-1:1))*j;
im2=real(x2-x2(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


home