CCORL
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.
|
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
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
