Perform Convolution in Frequency Domain Using FFTW?

Question

I'm trying to convolve two signals $x(n)$ and $h(n)$ in C by using the FFTW library's functions to perform a Fourier transform on each, multiply the appropriate complex components together, and take the IFFT of the resultant product. Unfortunately, my IFFT'd result isn't correct, and I'm not sure what I'm doing incorrectly.

I'm testing using a modified version of code I found online (.c and Makefile here), which populates two arrays with random numbers between 0 and 1. My process is as follows:

1. Create double buffers for x (audio signal) and h (impulse response) of lengths n and nIR respectively.
2. Populate each buffer with random numbers between 0 and 1.
3. Create two complex buffers X and H, each of length nOut = n + nIR - 1
4. Create and then execute two real-to-complex FFT plans, one for each of the signals, both using nOut elements (i.e. plan_forward = fftw_plan_dft_r2c_1d ( nOut, x, X, FFTW_ESTIMATE ); and plan_forwardIR = fftw_plan_dft_r2c_1d (nOut, h, H, FFTW_ESTIMATE);
5. Create complex array fftMulti to hold the frequency products of the two arrays. Its size is nc = nOut / 2 + 1 because FFTW doesn't store the redundant half of the FFT
6. Loop for ( i = 0; i < nc; i++ ) and perform complex multiplication, i.e. fftMulti[i][0] = X[i][0] * H[i][0] - X[i][1] * H[i][1]; for the real component and fftMulti[i][1] = X[i][0] * H[i][1] + X[i][1] * H[i][0]; for the imaginary component
7. Create double array convolvedSig of size nOut. This array will hold the IFFT of fftMulti
8. Create and execute a complex-to-real FFT plan of size nOut, i.e. plan_backwardConv = fftw_plan_dft_c2r_1d(nOut, fftMulti, convolvedSig, FFTW_ESTIMATE);
9. Normalize convolvedSig by dividing each element by nOut

The values in convolvedSig are not correct (they are way too high), but I'm not sure what I'm doing wrong. I also created arrays and complex-to-real plans for my unmodified signals X and H (also of size nOut), and those IFFTs worked just fine (i.e. values were exactly the same as before I performed FFTs on them).

Using my process and/or my code, can someone please help me identify what I'm doing incorrectly?

Answer 1

You need to zero pad x and h.

nZeros = length(x) + length(h) -1

Example in MATLAB:

clear all; close all; home

xSamples = 10;
hSamples = 7;

x = randn(1,xSamples);

h = randn(1,hSamples);

nZeros = xSamples + hSamples - 1;

X = fft(x,nZeros);
H = fft(h,nZeros);

x_conv_h = ifft(X.*H);


figure(1)
plot(real(x_conv_h))
hold on; grid on
plot(conv(x,h),'.r')
legend('fft based convolution','convolution')

Answer 2

FFT convolution of real signals is very easy. Multiplication in the frequency domain is equivalent to convolution in the time domain. The important thing to remember however, is that you are multiplying complex numbers, and therefore, you must to a "complex multiplication".

I will provide some C source for this below.

So the steps are:

1. Do an FFT of your filter kernel,

2. Do an FFT of your "dry" signal.

3. do a complex multiply of the two spectra

4. Perform the inverse FFT of this new spectrum.

Of course if you want to do continuous processing of lenghty signals, then you will need to use the overlap-add or overlap-save method.

If you are using real signals only, on an Intel format (little endian) machine, you can use Surreall FFT plus the multiply function I give below that is in the correct format for the data order.

Here is the code:

#define F_TYPE float;

void CompMulR(F_TYPE a[],F_TYPE b[],F_TYPE res[],long n)
  {
  long i;
  for(i=0;i<n;i++)
    {
    if(i<2)
      res[i]=a[i]*b[i]; // DC & NQ are not complex
    else
      {
      if((i&1)==0)
        {//real
        res[i]= a[i]*b[i] - a[i+1]*b[i+1];
        }
      else
        {//img.
        res[i]= a[i]*b[i-1] + a[i-1]*b[i];
        }
      }
    }
  }



#define N 1024
F_TYPE TimeDom [N];
F_TYPE FrqDom [N];
F_TYPE FrqDom2 [N];
F_TYPE Twidds[N/4];

InitFFT(N, Twidds); //Initialse the twiddle factors
FFT(TimeDom,FrqDom,N,Twidds); // do a forward transform
CompMulR(TimeDom,FrqDom,FrqDom2[],N);//do a complex multiply 
IFFT(TimeDom,FrqDom2,N,Twidds); // do an inverse transform 
//TimeDom now holds the convolution of the 2 signals

The code for the FFT is here: http://ravellescientific.co.uk/sureal-fft/ Hope this helps