Computes the convolution of two complex vectors.

Required Arguments

X — Complex vector of length NX.   (Input)

Y — Complex vector of length NY.   (Input)

Z — Complex vector of length NZ containing the convolution of X and Y.   (Output)

ZHAT — Complex vector of length NZ containing the discrete complex Fourier transform of Z.   (Output)

Optional Arguments

IDO — Flag indicating the usage of CCONV.   (Input)
   Default: IDO = 0.  

IDO     Usage

0          If this is the only call to CCONV.

If CCONV is called multiple times in sequence with the same NX, NY, and IPAD, IDO should be set to:

1          on the first call

2          on the intermediate calls

            3                              on the final call.

NX — Length of the vector X.   (Input)
Default: NX = size (X,1).

NY — Length of the vector Y.   (Input)
Default: NY = size (Y,1).

IPAD — IPAD should be set to zero for periodic data or to one for nonperiodic data.   (Input)
Default: IPAD =0.

NZ — Length of the vector Z.   (Input/Output)
Upon input: When IPAD is zero, NZ must be at least MAX(NX, NY). When IPAD is one, NZ must be greater than or equal to the smallest integer greater than or equal to
(NX + NY - 1) of the form (2α) * (3β) * (5γ) where alpha, beta, and gamma are nonnegative integers. Upon output, the value for NZ that was used by CCONV.
Default: NZ = size (Z,1).

FORTRAN 90 Interface

Generic:          CALL CCONV (X, Y, Z, ZHAT [,…])

Specific:         The specific interface names are S_CCONV and D_CCONV.

FORTRAN 77 Interface

Single:            CALL CCONV (IDO, NX, X, NY, Y, IPAD, NZ, Z, ZHAT)

Double:          The double precision name is DCCONV.

Description

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

nz := max{nx, ny}

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.

Comments

1.         Workspace may be explicitly provided, if desired, by use of C2ONV/DC2ONV. The reference is:

CALL C2ONV (IDO, NX, X, NY, Y, IPAD, NZ, Z, ZHAT, XWK, YWK, WK)

The additional arguments are as follows:

XWK — Complex work array of length NZ.

YWK — Complex work array of length NZ.

WK — Real work array of length 6 * NZ + 15.

2.         Informational error
Type                                   Code

   4                  1       The length of the vector Z must be large enough to hold the results. An acceptable length is returned in NZ.

Example

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) = ex f1 (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.

 

      USE IMSL_LIBRARIES

 

      IMPLICIT   NONE

      INTEGER    NFLTR, NY

      PARAMETER  (NFLTR=5, NY=100)

!

      INTEGER    I, IPAD, K, MOD, NOUT, NZ

      REAL       CABS, COS, EXP, FLOAT, FLTRER, ORIGER,  &

                SIN, TOTAL1, TOTAL2, TWOPI, X, T1, T2

      COMPLEX    CMPLX, F1, F2, FLTR(NFLTR), Y(NY), Z(2*(NFLTR+NY-1)), &

                ZHAT(2*(NFLTR+NY-1))

      INTRINSIC  CABS, CMPLX, COS, EXP, FLOAT, MOD, SIN

!                                DEFINE FUNCTIONS

      F1(X) = CMPLX(COS(X),SIN(X))

      F2(X) = EXP(X)*CMPLX(COS(X),SIN(X))

!

      CALL RNSET (1234579)

      CALL UMACH (2, NOUT)

      TWOPI = CONST('PI')

      TWOPI = 2.0*TWOPI

!                                 SET UP THE FILTER

      CALL CSET(NFLTR,(0.2,0.0),FLTR,1)

!                                 SET UP Y-VECTOR FOR THE PERIODIC

!                                 CASE.

      DO 20  I=1, NY

         X    = TWOPI*FLOAT(I-1)/FLOAT(NY-1)

         T1   = RNUNF()

         T2   = RNUNF()

         Y(I) = F1(X) + CMPLX(0.5*T1-0.25,0.5*T2-0.25)

   20 CONTINUE

!                                 CALL THE CONVOLUTION ROUTINE FOR THE

!                                 PERIODIC CASE.

      NZ = 2*(NFLTR+NY-1)

      CALL CCONV (FLTR, Y, Z, ZHAT)

!                                 PRINT RESULTS

      WRITE (NOUT,99993)

      WRITE (NOUT,99995)

      TOTAL1 = 0.0

      TOTAL2 = 0.0

      DO 30  I=1, NY

!                                 COMPUTE THE OFFSET FOR THE Z-VECTOR

         IF (I .GE. NY-1) THEN

            K = I - NY + 2

         ELSE

            K = I + 2

         END IF

!

         X      = TWOPI*FLOAT(I-1)/FLOAT(NY-1)

         ORIGER = CABS(Y(I)-F1(X))

         FLTRER = CABS(Z(K)-F1(X))

         IF (MOD(I,11) .EQ. 1) WRITE (NOUT,99997) X, F1(X), ORIGER, &

            FLTRER

         TOTAL1 = TOTAL1 + ORIGER

         TOTAL2 = TOTAL2 + FLTRER

   30 CONTINUE

      WRITE (NOUT,99998) TOTAL1/FLOAT(NY)

      WRITE (NOUT,99999) TOTAL2/FLOAT(NY)

!                                 SET UP Y-VECTOR FOR THE NONPERIODIC

!                                 CASE.

      DO 40  I=1, NY

         X    = FLOAT(I-1)/FLOAT(NY-1)

         T1   = RNUNF()

         T2   = RNUNF()

         Y(I) = F2(X) + CMPLX(0.5*T1-0.25,0.5*T2-0.25)

   40 CONTINUE

!                                 CALL THE CONVOLUTION ROUTINE FOR THE

!                                 NONPERIODIC CASE.

      NZ = 2*(NFLTR+NY-1)

      CALL CCONV (FLTR, Y, Z, ZHAT, IPAD=1)

!                                 PRINT RESULTS

      WRITE (NOUT,99994)

      WRITE (NOUT,99996)

      TOTAL1 = 0.0

      TOTAL2 = 0.0

      DO 50  I=1, NY

         X      = FLOAT(I-1)/FLOAT(NY-1)

         ORIGER = CABS(Y(I)-F2(X))

         FLTRER = CABS(Z(I+2)-F2(X))

         IF (MOD(I,11) .EQ. 1) WRITE (NOUT,99997) X, F2(X), ORIGER, &

            FLTRER

         TOTAL1 = TOTAL1 + ORIGER

         TOTAL2 = TOTAL2 + FLTRER

   50 CONTINUE

      WRITE (NOUT,99998) TOTAL1/FLOAT(NY)

      WRITE (NOUT,99999) TOTAL2/FLOAT(NY)

99993 FORMAT (' Periodic Case')

99994 FORMAT (/, ' Nonperiodic Case')

99995 FORMAT (8X, 'x', 15X, 'f1(x)', 8X, 'Original Error', 5X, &

            'Filtered Error')

99996 FORMAT (8X, 'x', 15X, 'f2(x)', 8X, 'Original Error', 5X, &

            'Filtered Error')

99997 FORMAT (1X, F10.4, 5X, '(', F7.4, ',', F8.4, ' )', 5X, F8.4, &

            10X, F8.4)

99998 FORMAT (' Average absolute error before filter:', F11.5)

99999 FORMAT (' Average absolute error after filter:', F12.5)

      END

Output

 

Periodic Case
 x               f1(x)        Original Error     Filtered Error
 0.0000     ( 1.0000,  0.0000 )       0.1666            0.0773
 0.6981     ( 0.7660,  0.6428 )       0.1685            0.1399
 1.3963     ( 0.1736,  0.9848 )       0.1756            0.0368
 2.0944     (-0.5000,  0.8660 )       0.2171            0.0142
 2.7925     (-0.9397,  0.3420 )       0.1147            0.0200
 3.4907     (-0.9397, -0.3420 )       0.0998            0.0331
 4.1888     (-0.5000, -0.8660 )       0.1137            0.0586
 4.8869     ( 0.1736, -0.9848 )       0.2217            0.0843
 5.5851     ( 0.7660, -0.6428 )       0.1831            0.0744
 6.2832     ( 1.0000,  0.0000 )       0.3234            0.0893
 Average absolute error before filter:    0.19315
 Average absolute error after filter:     0.08296

Nonperiodic Case
 x               f2(x)        Original Error     Filtered Error
 0.0000     ( 1.0000,  0.0000 )       0.0783            0.4336
 0.1111     ( 1.1106,  0.1239 )       0.2434            0.0477
 0.2222     ( 1.2181,  0.2752 )       0.1819            0.0584
 0.3333     ( 1.3188,  0.4566 )       0.0703            0.1267
 0.4444     ( 1.4081,  0.6706 )       0.1458            0.0868
 0.5556     ( 1.4808,  0.9192 )       0.1946            0.0930
 0.6667     ( 1.5307,  1.2044 )       0.1458            0.0734
 0.7778     ( 1.5508,  1.5273 )       0.1815            0.0690
 0.8889     ( 1.5331,  1.8885 )       0.0805            0.0193
 1.0000     ( 1.4687,  2.2874 )       0.2396            1.1708
 Average absolute error before filter:    0.18549
 Average absolute error after filter:     0.09636

Visual Numerics, Inc. Visual Numerics - Developers of IMSL and PV-WAVE http://www.vni.com/ PHONE: 713.784.3131 FAX:713.781.9260