========================================================= --- 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