Estimates the nonnormalized cross‑spectral density of two stationary time series using a spectral window given the time series data.

Routine CSSWD estimates the nonnormalized cross‑spectral density function of two jointly stationary time series using a spectral window given a sample of n = NOBS observations {Xt} and {Yt} for t = 1, 2, …, n.

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 modified periodogram of

is estimated by

where

and

represent the

cosine and sine transforms, respectively, and K is the scale factor equal to 1/(2πn). The modified periodogram of

is estimated 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

The routine PFFT is used to compute the modified periodograms of both

The computational formula for the cross periodogram is given by

where

and

The routine CPFFT is used to compute the modified cross periodogram between

The nonnormalized spectral density of Xt is estimated by

and the nonnormalized spectral density of Yt is estimated by

where the spectral window Wn(θ) is specified by argument ISWVER. The following spectral windows Wn(θ) are available.

where FM(θ) corresponds to the Fejér kernel of order M.

where DM(θ) represents the Dirichlet kernel. The Tukey‑Hamming window is obtained when
a = 0.23 and the Tukey‑Hanning window is obtained when a = 0.25.

where M is even. If M is odd, then M + 1 is used instead of M in the above formula.

The argument NM specifies the number of window parameters M and, hence, corresponds to the number of spectral density estimates to be computed for a given spectral window. Note that the same spectral window Wn(θ) and set of parameters M are used to obtain both

The above spectral density formulas assume the data {Xt} and {Yt} correspond to a realization of a bivariate discrete‑parameter stationary process observed consecutively in time. In this case, the observations are equally spaced in time with interval Δt = TINT equal to one. However, if the data correspond to a realization of a bivariate continuous‑parameter stationary process recorded at equal time intervals, then the spectral density estimates must be adjusted for the effect of aliasing. In general, the estimate of hX(ω) is given by

and the estimate of hY(ω) is given by

The nonnormalized spectral density is estimated over the set of frequencies

where nƒ = NF. These frequencies are in the scale specified by the argument IFSCAL but are transformed to the scale of radians per unit time for computational purposes. The frequency ω of the desired spectral estimate is assumed to be input in a form already adjusted for the time interval Δt.

The cross‑spectral density function is complex‑valued in general and may be written in the following form:

The cospectrum is estimated by

and the quadrature spectrum is estimated by

Note that the same spectral window Wn(θ) and window parameter M used to derive

are also used to compute

The nonnormalized cross‑spectral density estimate is computed over the same set of frequencies as the nonnormalized spectral density estimates with a similar adjustment for Δt.

An equivalent representation of hXY(ω) is the polar form defined by

The cross-amplitude spectrum is estimated by

and the phase spectrum is estimated by

Finally, the coherency spectrum is estimated by

The coherence or squared coherency is output.

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 CSSWD to these data produces the following results: