Computes the convolution of two complex vectors.

The subroutine CCONV computes the discrete convolution of two complex sequences x and y. More precisely, let nx be the length of x and ny denote the length of y. If a circular convolution is desired, then IPAD must be set to zero. We set

and we pad out the shorter vector with zeroes. Then, we compute

where the index on x is interpreted as a positive number between 1 and nz, modulo nz.

The technique used to compute the zi’s is based on the fact that the (complex discrete) Fourier transform maps convolution into multiplication. Thus, the Fourier transform of z is given by

where

The technique used here to compute the convolution is to take the discrete Fourier transform of x and y, multiply the results together component-wise, and then take the inverse transform of this product. It is very important to make sure that nz is a product of small primes if IPAD is set to zero. If nz is a product of small primes, then the computational effort will be proportional to nz log(nz). If IPAD is one, then a a good value is chosen for nz so that the Fourier transforms are efficient and nz ≥ nx + ny - 1. This will mean that both vectors will be padded with zeroes.

In this example, we compute both a periodic and a non-periodic convolution. The idea here is that one can compute a moving average of the type found in digital filtering using this routine. The averaging operator in this case is especially simple and is given by averaging five consecutive points in the sequence. The periodic case tries to recover a noisy function f1 (x) = cos(x) + i sin(x) by averaging five nearby values. The nonperiodic case tries to recover the values of the function f2(x) = exf1 (x) contaminated by noise. The large error for the first and last value printed has to do with the fact that the convolution is averaging the zeroes in the “pad” rather than function values. Notice that the signal size is 100, but we only report the errors at ten points.