The IMSL_CONVOL1D function computes the discrete convolution of two one-dimensional arrays.
This routine requires an IDL Advanced Math and Stats license. For more information, contact your sales or technical support representative.
Let nx be the length of x, and ny denote the length of y. If the keyword PERIODIC is set, then nz = max{nx, ny}, otherwise nz is set to the smallest whole number, nz ≥ nx + ny – 1, of the form:
The arrays x and y are then zero-padded to a length nz. Then, we compute:
where the index on x is interpreted as a nonnegative number between 0 and nz – 1.
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 following equation is true:
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 the product of small primes if PERIODIC is set. If nz is a product of small primes, then the computational effort will be proportional to nz log (nz). If PERIODIC is not set, then nz is chosen to be a product of small primes.
We point out that if x and y are not complex, then no complex transforms of x or y are taken, since a real transforms can simulate the complex transform above. Such a strategy is six times faster and requires less space than when using the complex transform.
Example
This example computes simple moving-average digital filter plots of 5-point and 25- point moving-average filters of noisy data. Results are shown in the figures that follow the example.
PRO Convol1d_ex1
IMSL_RANDOMOPT, SET = 1234579L
ny = 100
t = FINDGEN(ny)/(ny-1)
y = SIN(2*!PI*t) + .5*IMSL_RANDOM(ny, /Uniform) -.25
win=0
FOR nfltr = 5, 25, 20 DO BEGIN
nfltr_str = strcompress(nfltr,/Remove_All)
fltr = fltarr(nfltr)
fltr(*) = 1./nfltr
z = IMSL_CONVOL1D(fltr, y, /Periodic)
WINDOW, win++
PLOT, y, LINESTYLE = 1, TITLE = nfltr_str + $
'-point Moving Average'
OPLOT, shift(z, -nfltr/2)
ENDFOR
END
Syntax
Result = IMSL_CONVOL1D(X, Y [, DIRECT=value] [, PERIODIC=value])
Return Value
A one-dimensional array containing the discrete convolution of x and y.
Arguments
X
One-dimensional array.
Y
One-dimensional array.
Keywords
DIRECT (optional)
If present and nonzero, causes the computations to be done by the direct method instead of the FFT method regardless of the size of the vectors passed in.
PERIODIC (optional)
If present and nonzero, then a circular convolution is computed.
Version History