CCONV
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 (2a) * (3β) * (5y) 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 |
Description |
|
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) = 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.
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
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
