SSWP
Estimates the nonnormalized spectral density of a stationary time series using a spectral window given the periodogram.
Required Arguments
N — Number of observations in the centered and padded time series X. (Input)
N must be greater than or equal to two.
PX — Vector of length ⌊N/2⌋ + 1 containing the (modified) periodogram of X. (Input)
The periodogram ordinate evaluated at (angular) frequency wk = 2πk/N is given by
PX(k + 1), k = 0, 1, …, ⌊N/2⌋.
F — Vector of length NF containing the (angular) frequencies at which the spectral density is estimated. (Input)
M — Spectral window parameter. (Input)
M must be greater than or equal to one and less than N.
SX — Vector of length NF containing the estimate of the spectral density of the time series X. (Output)
Optional Arguments
NF — Number of (angular) frequencies. (Input)
NF must be greater than or equal to one.
Default: NF = size (F,1).
ISWVER — Option for version of the spectral window. (Input)
Default: ISWVER = 1.
ISWVER | Action |
|---|
1 | Modified Bartlett |
2 | Daniell |
3 | Tukey‑Hamming |
4 | Tukey‑Hanning |
5 | Parzen |
6 | Bartlett‑Priestley |
Refer to the “Algorithm” section for further details.
FORTRAN 90 Interface
Generic: CALL SSWP (N, PX, F, M, SX [, …])
Specific: The specific interface names are S_SSWP and D_SSWP.
FORTRAN 77 Interface
Single: CALL SSWP (N, PX, NF, F, ISWVER, M, SX)
Double: The double precision name is DSSWP.
Description
Routine SSWP estimates the nonnormalized spectral density function of a stationary time series using a spectral window given the modified periodogram of the appropriately centered and padded data
The routine
PFFT may be used to obtain the modified periodogram
over the discrete set of nonnegative frequencies
The symmetry of the periodogram is used to recover the ordinates at negative frequencies.
The estimate of the nonnormalized spectral density hX(ω) is computed according to
where the spectral window Wn(θ) is specified by argument ISWVER. The following spectral windows Wn(θ) are available.
Modified Bartlett
where FM(θ) corresponds to the Fejér kernel of order M.
Daniell
Tukey
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.
Parzen
where M is even. If M is odd, then M + 1 is used instead of M in the above formula.
Bartlett-Priestley
Only one window parameter M may be specified so that only one estimate of hX(ω) is computed. The nonnormalized spectral density is estimated over the set of frequencies
ω = ƒi, i = 1, …, nƒ
where nƒ = NF. These frequencies are in the scale of radians per unit time. The time sampling interval Δt is assumed to be equal to one.
Comments
1. The periodogram of
X may be computed using the routine
PFFT. Estimation of the spectral density of
X using the modified periodogram preserves the scale of the spectral density up to adjustment for the time sampling interval.
2. The time sampling interval, TINT, is assumed to be equal to one. This assumption is appropriate for discrete parameter processes. The adjustment for continuous parameter processes (TINT > 0.0) involves multiplication of the frequency vector F by 1/TINT and multiplication of the spectral density estimate by TINT.
3. To convert the frequency scale from radians per unit time to cycles per unit time, multiply F by 1/(2π).
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 SSWP to these data produces the following results:
USE IMSL_LIBRARIES
IMPLICIT NONE
INTEGER LDPM, LDSM, NF, NM, NOBS
REAL PI, NPAD
PARAMETER (NF=20, NM=3, NOBS=100, PI=3.141592654, &
LDPM=NOBS, LDSM=NF)
!
INTEGER I, IPVER, ISWVER, J, M(NM), N, NCOL, NROW
REAL F(NF), PM(LDPM,5), PX(LDPM), RDATA(176,2), FLOAT, &
SM(NF,5), SX(NF), X(NOBS)
CHARACTER CLABEL(6)*9, FMT*20, RLABEL(1)*6, TITLE*60
INTRINSIC FLOAT
!
EQUIVALENCE (PX(1), PM(1,3)), (F(1), SM(1,1))
EQUIVALENCE (X(1), RDATA(22,2))
!
DATA RLABEL/'NONE'/, CLABEL/' ', 'Frequency', 'Period', &
'M = 10', 'M = 20', 'M = 30'/
! Wolfer Sunspot Data for
! years 1770 through 1869
CALL GDATA (2, RDATA, NROW, NCOL)
! Center on arithmetic mean
! Pad standard amount
NPAD = NOBS-1
! Frequency in radians per unit time
! Modified periodogram version
IPVER = 1
! Compute periodogram
CALL PFFT (X, PM, IPVER=IPVER)
! Number of observations used to
! compute the periodogram
N = NOBS + NPAD
! Determine frequency and period
! at which to evaluate the spectral
! density
DO 10 I=1, NF
SM(I,1) = PI*FLOAT(I)/FLOAT(NF)
SM(I,2) = 2.0*FLOAT(NF)/FLOAT(I)
10 CONTINUE
! Spectral window parameters
M(1) = 10
M(2) = 20
M(3) = 30
! Compute spectral density using
! the Parzen window
ISWVER = 5
DO 20 J=1, NM
CALL SSWP (N, PX, F, M(J), SX, ISWVER=ISWVER)
! Copy into SM
CALL SCOPY (NF, SX, 1, SM(1:,2+J), 1)
20 CONTINUE
! Print results
TITLE = 'Spectral Density Using the Parzen Window'
FMT = '(F9.4, F6.2, 3F10.2)'
CALL WRRRL (TITLE, SM, RLABEL, CLABEL, FMT=FMT)
! Compute spectral density using
! the Bartlett-Priestley window
ISWVER = 6
DO 30 J=1, NM
CALL SSWP (N, PX, F, M(J), SX, ISWVER=ISWVER)
! Copy into SM
CALL SCOPY (NF, SX, 1, SM(1:,2+J), 1)
30 CONTINUE
! Print results
TITLE = '%/Spectral Density Using the Bartlett-Priestley '// &
& 'Window'
CALL WRRRL (TITLE, SM, RLABEL, CLABEL, FMT=FMT)
!
END
Output
Spectral Density Using the Parzen Window
Frequency Period M = 10 M = 20 M = 30
0.1571 40.00 659.64 617.42 619.73
0.3142 20.00 666.95 554.70 339.61
0.4712 13.33 653.73 770.64 860.49
0.6283 10.00 598.77 857.80 1046.13
0.7854 8.00 497.47 582.85 550.77
0.9425 6.67 367.72 266.33 186.98
1.0996 5.71 240.65 121.46 104.79
1.2566 5.00 142.41 76.17 76.74
1.4137 4.44 81.28 54.20 47.19
1.5708 4.00 49.13 40.16 41.39
1.7279 3.64 32.57 27.58 26.46
1.8850 3.33 22.44 16.52 14.40
2.0420 3.08 15.53 10.93 9.87
2.1991 2.86 11.19 8.30 8.32
2.3562 2.67 8.66 6.18 5.86
2.5133 2.50 6.93 4.75 4.22
2.6704 2.35 5.51 4.62 4.35
2.8274 2.22 4.47 4.91 5.24
2.9845 2.11 3.61 4.23 4.75
3.1416 2.00 2.62 2.44 2.27
Spectral Density Using the Bartlett-Priestley Window
Frequency Period M = 10 M = 20 M = 30
0.1571 40.00 604.34 712.73 757.61
0.3142 20.00 564.28 176.81 107.08
0.4712 13.33 767.63 927.14 981.10
0.6283 10.00 900.32 1190.30 1172.23
0.7854 8.00 607.45 494.85 571.65
0.9425 6.67 237.16 127.65 87.36
1.0996 5.71 103.34 113.93 135.34
1.2566 5.00 75.74 74.88 57.57
1.4137 4.44 52.64 44.98 38.59
1.5708 4.00 38.50 44.56 50.59
1.7279 3.64 27.35 25.28 21.76
1.8850 3.33 15.68 13.84 13.10
2.0420 3.08 10.33 9.79 7.41
2.1991 2.86 7.95 8.31 8.67
2.3562 2.67 6.04 5.86 7.08
2.5133 2.50 4.56 3.67 2.90
2.6704 2.35 4.44 4.38 4.06
2.8274 2.22 4.99 5.62 5.40
2.9845 2.11 4.31 5.07 5.08
3.1416 2.00 2.43 2.23 2.44
