RCORL
Computes the correlation of two real vectors.
Required Arguments
X — Real vector of length N. (Input)
Y — Real vector of length N. (Input)
Z — Real vector of length NZ containing the correlation of X and Y. (Output)
ZHAT — Real vector of length NZ containing the discrete Fourier transform of Z. (Output)
Optional Arguments
IDO — Flag indicating the usage of RCORL. (Input)
Default: IDO = 0.
|
IDO |
Usage |
|
0 |
If this is the only call to RCORL. |
|
If RCORL 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)
Default: IPAD = 0.
IPAD Value
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.
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 RCORL.
Default: NZ = size (Z,1).
FORTRAN 90 Interface
Generic: CALL RCORL (X, Y, Z, ZHAT [, …])
Specific: The specific interface names are S_RCORL and D_RCORL.
FORTRAN 77 Interface
Single: CALL RCORL (IDO, N, X, Y, IPAD, NZ, Z, ZHAT)
Double: The double precision name is DRCORL.
Description
The subroutine RCORL computes the discrete correlation of two 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.
We point out that no complex transforms of x or y are taken since both sequences are real, and we can take real transforms and simulate the complex transform above. This can produce a savings of a factor of six in time as well as save space over using the complex transform.
Comments
1. Workspace may be explicitly provided, if desired, by use of R2ORL/DR2ORL. The reference is:
CALL R2ORL (IDO, N, X, Y, IPAD, NZ, Z, ZHAT, XWK, YWK, WK)
The additional arguments are as follows:
XWK — Real work array of length NZ.
YWK — Real work array of length NZ.
WK — Real work arrary of length 2 * 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) = 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) = 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, F1, F2, FLOAT, PI, SIN, X(N), XNORM, &
Y(N), YNORM, Z(4*N), ZHAT(4*N)
INTRINSIC FLOAT, SIN
! Define functions
F1(A) = SIN(A)
F2(A) = SIN(A*A)
!
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 RCORL (X, Y, Z, ZHAT)
! Find the element of Z with the
! largest normalized value.
XNORM = SNRM2(N,X,1)
YNORM = SNRM2(N,Y,1)
DO 20 I=1, N
Z(I) = Z(I)/(XNORM*YNORM)
20 CONTINUE
K = ISMAX(N,Z,1)
! Print results for the periodic
! case.
WRITE (NOUT,99995)
WRITE (NOUT,99994)
WRITE (NOUT,99997)
WRITE (NOUT,99998) K
WRITE (NOUT,99999) K, Z(K)
! Set up the vectors for the
! nonperiodic 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 RCORL (X, Y, Z, ZHAT, IPAD=1)
! Find the element of Z with the
! largest normalized value.
XNORM = SNRM2(N,X,1)
YNORM = SNRM2(N,Y,1)
DO 40 I=1, N
Z(I) = Z(I)/(XNORM*YNORM)
40 CONTINUE
K = ISMAX(N,Z,1)
! Print results for the nonperiodic
! case.
WRITE (NOUT,99996)
WRITE (NOUT,99994)
WRITE (NOUT,99997)
WRITE (NOUT,99998) K
WRITE (NOUT,99999) K, Z(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 (' value is Z(', I2, ').')
99999 FORMAT (' The normalized value of Z(', I2, ') is', F6.3)
END
Example #1: Periodic case
----------------------------
The element of Z with the largest normalized value is Z(26).
The normalized value of Z(26) is 1.000
Example #2: Nonperiodic case
----------------------------
The element of Z with the largest normalized value is Z(26).
The normalized value of Z(26) is 0.661
