CCORL

IDO	Usage
0	If this is the only call to CCORL.
If CCORL 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

N — Length of the X and Y vectors. (Input)
Default: N = size (X,1).

IPAD — IPAD should be set as follows. (Input)
IPAD = 0 for periodic data with X and Y different. IPAD = 1 for nonperiodic data with X and Y different. IPAD = 2 for periodic data with X and Y identical. IPAD = 3 for nonperiodic data with X and Y identical.
Default: IPAD = 0.

NZ — Length of the vector Z. (Input/Output)
Upon input: When IPAD is zero or two, NZ must be at least (2 * N ‑ 1). When IPAD is one or three, NZ must be greater than or equal to the smallest integer greater than or equal to (2 * N ‑ 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 CCORL.
Default: NZ = size (Z,1).

FORTRAN 90 Interface

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

Specific: The specific interface names are S_CCORL and D_CCORL.

FORTRAN 77 Interface

Single: CALL CCORL (IDO, N, X, Y, IPAD, NZ, Z, ZHAT)

Double: The double precision name is DCCORL.

Description

The subroutine CCORL computes the discrete correlation of two complex sequences x and y. More precisely, let n be the length of x and y. If a circular correlation is desired, then IPAD must be set to zero (for x and y distinct) and two (for x = y). We set (on output)

where α,β,y are nonnegative integers yielding the smallest number of the type 2α3β5y satisfying the inequality. Once nz is determined, we pad out the vectors with zeroes. Then, we compute

where the index on x is interpreted as a positive number between one and nz, modulo nz. Note that this means that

contains the correlation of x(⋅- k - 1) with y as k = 0, 1, …, nz /2. Thus, if x(k - 1) = y(k) for all k, then we would expect

to be the largest component of ℜz.

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

where

Thus, the technique used here to compute the correlation is to take the discrete Fourier transform of x and the conjugate of the discrete Fourier transform of 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 or two. If nz is a product of small primes, then the computational effort will be proportional to nz log(nz). If IPAD is one or three, then a good value is chosen for nz so that the Fourier transforms are efficient and nz ≥ 2n - 1. This will mean that both vectors will be padded with zeroes.

Comments

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

CALL C2ORL (IDO, N, X, 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 arrary 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 correlation between two distinct signals x and y. In the first case, we have 100 equally spaced points on the interval [0, 2π] and f1 (x) = cos(x) + i sin(x). We define x and y as follows

Note that the maximum value of z (the correlation of x with y) occurs at i = 26, which corresponds to the offset.

The nonperiodic case uses the function f2 (x) = cos(x2) + i sin(x2). The two input signals are on the interval [0, 4π].

The offset of x to y is again (roughly) 26 and this is where z has its maximum value.

USE IMSL_LIBRARIES

IMPLICIT NONE

INTEGER N

PARAMETER (N=100)

INTEGER I, IPAD, K, NOUT, NZ

REAL A, COS, F1, F2, FLOAT, PI, SIN, &

XNORM, YNORM, ZREAL1(4*N)

COMPLEX CMPLX, X(N), Y(N), Z(4*N), ZHAT(4*N)

INTRINSIC CMPLX, COS, FLOAT, SIN

! Define functions

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

F2(A) = CMPLX(COS(A*A),SIN(A*A))

CALL RNSET (1234579)

CALL UMACH (2, NOUT)

PI = CONST('pi')

! Set up the vectors for the

! periodic case.

DO 10 I=1, N

X(I) = F1(2.0*PI*FLOAT(I-1)/FLOAT(N-1))

Y(I) = F1(2.0*PI*FLOAT(I-1)/FLOAT(N-1)+PI/2.0)

10 CONTINUE

! Call the correlation routine for the

! periodic case.

NZ = 2*N

CALL CCORL (X, Y, Z, ZHAT, IPAD=0, NZ=NZ)

! Find the element of Z with the

! largest normalized real part.

XNORM = SCNRM2(N,X,1)

YNORM = SCNRM2(N,Y,1)

DO 20 I=1, N

ZREAL1(I) = REAL(Z(I))/(XNORM*YNORM)

20 CONTINUE

K = ISMAX(N,ZREAL1,1)

! Print results for the periodic

! case.

WRITE (NOUT,99995)

WRITE (NOUT,99994)

WRITE (NOUT,99997)

WRITE (NOUT,99998) K

WRITE (NOUT,99999) K, ZREAL1(K)

! Set up the vectors for the

! nonperioddic case.

DO 30 I=1, N

X(I) = F2(4.0*PI*FLOAT(I-1)/FLOAT(N-1))

Y(I) = F2(4.0*PI*FLOAT(I-1)/FLOAT(N-1)+PI)

30 CONTINUE

! Call the correlation routine for the

! nonperiodic case.

NZ = 4*N

CALL CCORL (X, Y, Z, ZHAT, IPAD=1, NZ=NZ)

! Find the element of z with the

! largest normalized real part.

XNORM = SCNRM2(N,X,1)

YNORM = SCNRM2(N,Y,1)

DO 40 I=1, N

ZREAL1(I) = REAL(Z(I))/(XNORM*YNORM)

40 CONTINUE

K = ISMAX(N,ZREAL1,1)

! Print results for the nonperiodic

! case.

WRITE (NOUT,99996)

WRITE (NOUT,99994)

WRITE (NOUT,99997)

WRITE (NOUT,99998) K

WRITE (NOUT,99999) K, ZREAL1(K)

99994 FORMAT (1X, 28('-'))

99995 FORMAT (' Case #1: periodic data')

99996 FORMAT (/, ' Case #2: nonperiodic data')

99997 FORMAT (' The element of Z with the largest normalized ')

99998 FORMAT (' real part is Z(', I2, ').')

99999 FORMAT (' The normalized value of real(Z(', I2, ')) is', F6.3)

END

Output

Example #1: periodic case

----------------------------

The element of Z with the largest normalized real part is Z(26).

The normalized value of real(Z(26)) is 1.000

Example #2: nonperiodic case

----------------------------

The element of Z with the largest normalized real part is Z(26).

The normalized value of real(Z(26)) is 0.638