CPFFT
Computes the cross periodogram of two stationary time series using a fast Fourier transform.
Required Arguments
X — Vector of length NOBS containing the first stationary time series. (Input)
Y — Vector of length NOBS containing the second stationary time series. (Input)
CPM — (⌊N/2⌋ + 1) by 10 matrix containing a summarization of the results of the cross periodogram analysis. (Output)
For k = 0, 1, …, ⌊N/2⌋, the (k + 1)-st element of the j‑th column of CPM is defined as
|
Col. |
Description |
|
1 |
Frequency, ωk where ωk = 2πk/N for IFSCAL = 0 or ωk = k/N for IFSCAL = 1. |
|
2 |
Period, pk where pk = 2π/ ωk for IFSCAL = 0 and pk = 1/ ωk for IFSCAL = 1. If ωk = 0, pk is set to missing. |
|
3 |
X periodogram ordinate, IX(ωk) |
|
4 |
X cosine transformation coefficient, AX(ωk) |
|
5 |
X sine transformation coefficient, BX(ωk) |
|
6 |
Y periodogram ordinate, IY(ωk) |
|
7 |
Y cosine transformation coefficient, AY(ωk) |
|
8 |
Y sine transformation coefficient, BY(ωk) |
|
9 |
Real part of the XY cross periodogram ordinate IXY(ωk). |
|
10 |
Imaginary part of the XY cross periodogram ordinate IXY(ωk). |
Optional Arguments
NOBS — Number of observations in each stationary time series X and Y. (Input)
NOBS must be greater than or equal to two.
Default: NOBS = size (X,1).
IPRINT — Printing option. (Input)
Default: IPRINT = 0.
|
IPRINT |
Action |
|
0 |
No printing is performed. |
|
1 |
Prints the periodogram, cosine and sine series, and the real and imaginary components of the cross periodogram. |
XCNTR — Constant used to center the time series X. (Input)
Default: XCNTR = the arithmetic mean.
YCNTR — Constant used to center the time series Y. (Input)
Default: YCNTR = the arithmetic mean.
NPAD — Number of zeroes used to pad each centered time series. (Input)
NPAD must be greater than or equal to zero. The length of each centered and padded time series is N = NOBS + NPAD.
Default: NPAD = NOBS – 1.
IFSCAL — Option for frequency scale. (Input)
Default: IFSCAL= 0.
|
IFSCAL |
Action |
|
0 |
Frequency in radians per unit time |
|
1 |
Frequency in cycles per unit time |
IPVER — Option for version of the periodogram. (Input)
Default: IPVER = 0.
|
IPVER |
Action |
|
0 |
Compute usual periodogram. |
|
1 |
Compute modified periodogram. |
Refer to the “Description” section for further details.
LDCPM — Leading dimension of CPM exactly as specified in the dimension statement of the calling program. (Input)
LDCPM must be greater than or equal to ⌊N/2⌋ + 1.
Default: LDCPM = size (CPM,1).
FORTRAN 90 Interface
Generic: CALL CPFFT (X, Y, CPM [, …])
Specific: The specific interface names are S_CPFFT and D_CPFFT.
FORTRAN 77 Interface
Single: CALL CPFFT (NOBS, X, Y, IPRINT, XCNTR, YCNTR, NPAD, IFSCAL, IPVER, CPM, LDCPM)
Double: The double precision name is DCPFFT.
Description
Routine CPFFT computes the cross periodogram of two jointly stationary time series given a sample of n = NOBS observations {Xt} and {Yt} for t = 1, 2, …, n.
Let
represent the centered and padded data where N = NOBS + NPAD,
and
is determined by
Similarly, let
represent the centered and padded data where
and
is determined by
The periodogram of the sample sequence {Xt}, t = 1, …, n computed with the padded sequence
is defined by
where
and
represent the
cosine and sine transforms, respectively, and K is the scale factor
The periodogram of the sample sequence {Yt}, t = 1, …, n computed with the padded sequence
is defined by
where
and
represent the
cosine and sine transforms, respectively. Since the periodogram is an even function of the frequency, it is sufficient to estimate the periodogram at the discrete set of nonnegative frequencies
(Here, ⌊a⌋ means the greatest integer less than or equal to a). The routine PFFT is used to compute the periodograms of both
according to the version specified by the argument IPVER. The computational formula for the cross periodogram is given by
where
and
The real part of the (modified) cross periodogram represents the ’raw’ sample cospectrum and the negative of the imaginary part of the (modified) cross periodogram represents the ‘raw’ sample quadrature spectrum (Priestley 1981, page 695). The relationship between the cross periodogram and its complex conjugate is given by
and may be used to recover the cross periodogram at negative frequencies.
Comments
1. Workspace may be explicitly provided, if desired, by use of C2FFT/DC2FFT. The reference is:
CALL C2FFT (NOBS, X, Y, IPRINT, XCNTR, YCNTR, NPAD, IFSCAL, IPVER, CPM, LDCPM, CX, COEF, WFFTC, CPY)
The additional arguments are as follows:
CX — Complex work vector of length N.
COEF — Complex work vector of length N.
WFFTC — Work vector of length 4N + 15.
CPY — Work vector of length 2N.
2. The centered and padded time series are defined by
|
I |
AOV(I) |
|
CX(j) = X(j) ‑ XCNTR |
for j = 1, …, NOBS |
|
CX(j) = 0 |
for j = NOBS + 1, …, N |
|
and |
|
|
CY(j) = Y(j) ‑ YCNTR |
for j = 1, …, NOBS |
|
CY (j) = 0 |
for j = NOBS + 1, …, N where N = NOBS + NPAD. |
3. The cross periodogram IXY(ω) is complex valued in general. The relation
IXY(‑ω) = conj(IXY(ω)) for w > 0.0 recovers the cross periodogram for negative frequencies since real(IXY(‑ω)) = real(IXY(ω)) and imag(IXY(‑ω)) = ‑imag(IXY(ω)). The periodogram I(ω) is an even function of the frequency ω. The relation I(‑ω) = I(ω) for ω > 0.0 recovers the periodogram for negative frequencies.
4. Since cos(ω) is an even function of ω and sin(ω) is an odd function of ω, the cosine and sine transformations, respectively, satisfy A(‑ω) = A(ω) and B(‑ω) = ‑B(ω) for ω > 0.0. Similarly, the complex Fourier coefficients, stored in COEF, satisfy COEF(‑ω) = conj(COEF(ω)).
5. Computation of the 2 * NOBS ‑ 1 cross‑covariances of X and Y using the inverse Fourier transform of the cross periodogram requires NPAD = NOBS ‑ 1.
Example
Consider the Robinson Multichannel Time Series Data (Robinson 1967, page 204) where X is the Wölfer sunspot number and Y is the northern light activity for the time period from 1770 through 1869. Application of routine CPFFT to these data produces the following results. Note that CPFFT sets CPM (1, 2) to the missing value code via routine AMACH.The printing of CPM (1, 2) depends on the computer.
USE GDATA_INT
USE CPFFT_INT
USE WRRRL_INT
IMPLICIT NONE
INTEGER LDCPM, LDRDAT, NDRDAT, NOBS, NPAD
PARAMETER (LDRDAT=100, NDRDAT=4, NOBS=100, NPAD=NOBS-1, &
LDCPM=(NOBS+NPAD)/2+1)
!
INTEGER IPVER, NRCOL, NRROW
REAL CPM(LDCPM,10), FLOAT, RDATA(LDRDAT,NDRDAT), &
X(NOBS), Y(NOBS)
CHARACTER CLABEL1(6)*9, CLABEL2(6)*9, FMT*7, RLABEL(1)*6, &
TITLE*41
INTRINSIC FLOAT
!
EQUIVALENCE (X(1), RDATA(1,2)), (Y(1), RDATA(1,3))
!
DATA TITLE/'Results of the Cross Periodogram Analysis'/
DATA FMT/'(F10.3)'/
DATA CLABEL1/'k+1', 'w(k)', 'p(k)', 'IX(w(k))', 'AX(w(k))', &
'BX(w(k))'/
DATA CLABEL2/'k+1', 'IY(w(k))', 'AY(w(k))', 'BY(w(k))', &
'Real IXY', 'Imag. IXY'/
DATA RLABEL/'NUMBER'/
!
! Robinson Data
CALL GDATA (8, RDATA, NRROW, NRCOL)
! Center on arithmetic means
! Frequency in radians per unit time
! Modified periodogram version
IPVER = 1
! Compute the cross periodogram
CALL CPFFT (X, Y, CPM, IPVER=IPVER)
!
! Print results (First 10 rows)
CALL WRRRL (TITLE, CPM, RLABEL, CLABEL1, 10, 5, FMT=FMT)
CALL WRRRL ('%/', CPM(1:,6), RLABEL, CLABEL2, 10, 5, FMT=FMT)
!
END
Output
Results of the Cross Periodogram Analysis
k+1 w(k) p(k) IX(w(k)) AX(w(k)) BX(w(k))
1 0.000 NaN 0.000 0.000 0.000
2 0.032 199.000 184.159 3.742 -13.044
3 0.063 99.500 1364.408 35.457 -10.354
4 0.095 66.333 2433.933 29.411 39.610
5 0.126 49.750 1351.002 -21.749 29.631
6 0.158 39.800 140.421 -11.716 -1.773
7 0.189 33.167 44.117 -4.671 4.722
8 0.221 28.429 121.186 -11.003 -0.343
9 0.253 24.875 176.275 -4.782 -12.386
10 0.284 22.111 144.867 10.038 -6.642
k+1 IY(w(k)) AY(w(k)) BY(w(k)) Real IXY Imag. IXY
1 0.000 0.000 0.000 0.000 0.000
2 1689.212 -37.480 -16.866 79.776 -552.014
3 4113.003 41.232 -49.122 1970.577 -1314.779
4 3255.785 44.214 36.068 2729.031 -690.474
5 1757.663 -8.162 41.122 1396.006 -652.513
6 1002.050 -30.107 9.778 335.410 -167.954
7 62.360 -6.825 3.972 50.636 13.678
8 1481.396 -38.096 5.487 417.288 -73.451
9 1274.161 -17.176 -31.291 469.704 -63.095
10 488.479 -12.442 -18.267 -3.570 -265.992
