Computes the correlation of two complex vectors.
Required Arguments
X — Complex vector of length N. (Input)
Y — Complex vector of length N. (Input)
Z — Complex vector of length NZ containing the correlation of X and Y. (Output)
ZHAT — Complex vector of length NZ containing the inverse discrete complex Fourier transform of Z. (Output)
Optional Arguments
IDO — Flag
indicating the usage of CCORL.
(Input)
Default: IDO = 0.
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:
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
(2α) * (3β) * (5γ) 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 α, β, g are nonnegative integers yielding the smallest number of the type 2α3β5γ 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
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
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.
REAL A, COS, F1, F2, FLOAT, PI, SIN, &
COMPLEX CMPLX, X(N), Y(N), Z(4*N), ZHAT(4*N)
INTRINSIC CMPLX, COS, FLOAT, SIN
F2(A) = CMPLX(COS(A*A),SIN(A*A))
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)
! Call the correlation routine for the
CALL CCORL (X, Y, Z, ZHAT, IPAD=0, NZ=NZ)
! Find the element of Z with the
! largest normalized real part.
ZREAL1(I) = REAL(Z(I))/(XNORM*YNORM)
! Print results for the periodic
WRITE (NOUT,99999) K, ZREAL1(K)
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)
! Call the correlation routine for the
CALL CCORL (X, Y, Z, ZHAT, IPAD=1, NZ=NZ)
! Find the element of z with the
! largest normalized real part.
ZREAL1(I) = REAL(Z(I))/(XNORM*YNORM)
! Print results for the nonperiodic
WRITE (NOUT,99999) K, ZREAL1(K)
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)
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
|
Visual Numerics, Inc. PHONE: 713.784.3131 FAX:713.781.9260 |