PFFT
Computes the periodogram of a stationary time series using a fast Fourier transform.
Required Arguments
X — Vector of length NOBS containing the stationary time series. (Input)
PM — (⌊N/2⌋ + 1) by 5 matrix that contains a summarization of the periodogram analysis. (Output)
For k = 0, 1, …, ⌊N/2⌋, the (k + 1)-st element of the j‑th column of PM is defined as:
|
Col |
Description |
|
1 |
Frequency, ωk where ωk = 2πk/N for IFSCAL = 0 and ω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 |
Periodogram ordinate, I(ωk). |
|
4 |
Cosine transformation coefficient, A(ωk). |
|
5 |
Sine transformation coefficient, B(ωk). |
Optional Arguments
NOBS — Number of observations in the stationary time series X. (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, and the cosine and sine transformations of the centered and padded time series. |
XCNTR — Constant used to center the time series X. (Input)
Default: XCNTR = the arithmetic mean.
NPAD — Number of zeroes used to pad the centered time series. (Input)
NPAD must be greater than or equal to zero. The length of the 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.
LDPM — Leading dimension of PM exactly as specified in the dimension statement of the calling program. (Input)
LDPM must be greater than or equal to ⌊N/2⌋ + 1.
Default: LDPM = size (PM,1).
FORTRAN 90 Interface
Generic: CALL PFFT (X, PM [, …])
Specific: The specific interface names are S_PFFT and D_PFFT.
FORTRAN 77 Interface
Single: CALL PFFT (NOBS, X, IPRINT, XCNTR, NPAD, IFSCAL, IPVER, PM, LDPM)
Double: The double precision name is DPFFT.
Description
Routine PFFT computes the periodogram of a stationary time series given a sample of n = NOBS observations {Xt} for t = 1, 2, …, n.
Let
for t = 1, …, N represent the centered and padded data where N = NOBS + NPAD,
and
is determined by
The discrete Fourier transform of
for t = 1, …, N is defined by
over the discrete set of frequencies
An alternative representation of
in terms of cosine and sine transforms is
where
and
The fast Fourier transform algorithm is used to compute the discrete Fourier transform. The periodogram of the sample sequence {Xt}, t = 1, …, n computed with the centered and padded sequence
t = 1, …, N is defined by
where K is the scale factor
The scale factor of the usual periodogram relates the ordinates to the sum of squares of
(Fuller 1976, pages 276–277). If the first ordinate (corresponding to k = 0) is replaced by one‑half of its value, then if N is odd, the sum of the ⌊N/2⌋ + 1 ordinates corresponding to k = 0, 1, …, ⌊N/2⌋ is
For N even, if the first ordinate (corresponding to k = 0) and the last ordinate (corresponding to k = N/2) are each replaced by one‑half of their values, then the same relationship holds. The modified periodogram is an asymptotically unbiased estimate of the nonnormalized spectral density function at each frequency ωk (Priestley 1981, page 417). The argument IPVER is used to specify the version of the periodogram.
The alternative representation of the discrete Fourier transform implies
where
and
represent the (scaled) 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
Use of the centered data
(without padding) instead of the original data {Xt} for t = 1, …, n does not affect the asymptotic sampling properties of the periodogram. In fact,
For ωk = 0, both
and
reflect the mean of the data. See Priestley (1981, page 417) for further discussion.
Comments
1. Workspace may be explicitly provided, if desired, by use of P2FT/DP2FT. The reference is:
CALL P2FT (NOBS, X, IPRINT, XCNTR, NPAD, IFSCAL, IPVER, PM, LDPM, 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 is defined by
|
CX(j) = X(j) ‑ XCNTR |
for j = 1, …, NOBS |
|
CX(j) = 0 |
for j = NOBS + 1, …, N |
|
where N = NOBS + NPAD. |
|
3. 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 autocovariances of X using the inverse Fourier transform of the periodogram requires NPAD = NOBS ‑ 1.
Example
Consider the Wölfer Sunspot Data (Anderson 1971, page 660) consisting of the number of sunspots observed each year from 1749 through 1924. The data set for this example consists of the number of sunspots observed from 1770 through 1869. Application of routine PFFT to these data produces the following results.
USE GDATA_INT
USE PFFT_INT
USE WRRRL_INT
IMPLICIT NONE
INTEGER LDPM, NOBS
PARAMETER (LDPM=100, NOBS=100)
!
INTEGER IPVER, NCOL, NPAD, NROW
REAL PM(LDPM,5), RDATA(176,2), REAL, X(NOBS)
CHARACTER CLABEL(6)*9, FMT*20, RLABEL(1)*6
INTRINSIC REAL
!
EQUIVALENCE (X(1), RDATA(22,2))
!
DATA RLABEL/'NONE'/, CLABEL/' ', 'Frequency', 'Period', &
'I(w(k))', 'A(w(k))', 'B(w(k))'/
! Wolfer Sunspot Data for
! years 1770 through 1869
CALL GDATA (2, RDATA, NROW, NCOL)
! Center on arithmetic mean
! Pad standard amount
! Frequency in radians per unit time
! Modified periodogram version
IPVER = 1
! Compute the periodogram
CALL PFFT (X, PM, IPVER=IPVER)
! Print results
FMT = '(F9.4, F6.2, 3F10.2)'
CALL WRRRL (' ', PM, RLABEL, CLABEL, 20, 5, FMT=FMT)
!
END
Output
Frequency Period I(w(k)) A(w(k)) B(w(k))
0.0000 NaN 0.00 0.00 0.00
0.0316 199.00 183.97 3.72 -13.04
0.0631 99.50 1363.37 35.45 -10.32
0.0947 66.33 2427.09 29.31 39.60
0.1263 49.75 1346.64 -21.74 29.56
0.1579 39.80 139.74 -11.69 -1.79
0.1894 33.17 44.67 -4.65 4.80
0.2210 28.43 123.47 -11.11 -0.33
0.2526 24.88 176.04 -4.79 -12.37
0.2842 22.11 143.06 9.92 -6.69
0.3157 19.90 44.17 6.43 1.68
0.3473 18.09 38.95 5.40 3.13
0.3789 16.58 63.20 7.14 3.49
0.4105 15.31 537.64 0.89 23.17
0.4420 14.21 944.68 -30.73 -0.75
0.4736 13.27 162.02 -0.95 -12.69
0.5052 12.44 908.09 -24.51 -17.53
0.5368 11.71 3197.84 34.84 -44.54
0.5683 11.06 1253.82 19.69 29.43
0.5999 10.47 850.45 -8.75 -27.82
