Estimates the nonnormalized spectral density of a stationary time series using a spectral window given the periodogram.
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)
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.
Generic: CALL SSWP (N, PX, F, M, SX [, ])
Specific: The specific interface names are S_SSWP and D_SSWP.
Single: CALL SSWP (N, PX, NF, F, ISWVER, M, SX)
Double: The double precision name is DSSWP.
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.
where FM(θ) corresponds to the Fejer 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.
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.
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π).
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
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
PHONE: 713.784.3131 FAX:713.781.9260 |