IDL Analyst Reference Guide: Transforms

IMSL_CONVOL1D

The IMSL_CONVOL1D function computes the discrete convolution of two one-dimensional arrays.

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

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

If present and nonzero, then a circular convolution is computed.

Discussion

The IMSL_CONVOL1D function computes the discrete convolution of two sequences x and y.

Let n_x be the length of x, and n_y denote the length of y. If the keyword PERIODIC is set, then n_z = max{n_x, n_y}, otherwise n_z is set to the smallest whole number, n_z ³ n_x + n_y – 1, of the form:

The arrays x and y are then zero-padded to a length n_z. Then, we compute:

where the index on x is interpreted as a nonnegative number between 0 and n_z – 1.

The technique used to compute the z_i'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 n_z is a product of small primes, then the computational effort will be proportional to n_z log  (n_z). If PERIODIC is not set, then n_z 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 figures Figure 9-1 and Figure 9-2.

PRO Convol1d_ex1  
   IMSL_RANDOMOPT, SET = 1234579L  
   ; Set the random number seed.  
   ny = 100  
   t = FINDGEN(ny)/(ny-1)  
   y = SIN(2*!PI*t) + .5*IMSL_RANDOM(ny, /Uniform) -.25  
   ; Define a 1-period sine wave with added noise.  
   win=0  
   FOR nfltr = 5, 25, 20 DO BEGIN  
      nfltr_str = strcompress(nfltr,/Remove_All)  
      fltr = fltarr(nfltr)  
      fltr(*) = 1./nfltr  
      ; Define the NFLTR-point moving average array.  
      z = IMSL_CONVOL1D(fltr, y, /Periodic)  
      ; Convolve filter and signal, using keyword Periodic.  
      WINDOW, win++  
      PLOT, y, LINESTYLE = 1, TITLE = nfltr_str + $  
         '-point Moving Average'  
      OPLOT, shift(z, -nfltr/2)  
   ENDFOR  
END  

Figure 9-1: 5 Point Moving Average

Figure 9-1: 5 Point Moving Average

Figure 9-2: 25 Point Moving Average

Figure 9-2: 25 Point Moving Average

Version History

  IDL Online Help (March 06, 2007)