NAG Library Function Document

nag_tsa_spectrum_univar (g13cbc)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

+− 10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

1 Purpose

nag_tsa_spectrum_univar (g13cbc) calculates the smoothed sample spectrum of a univariate time series using spectral smoothing by the trapezium frequency (Daniell) window.

2 Specification

#include <nag.h>

#include <nagg13.h>

void	nag_tsa_spectrum_univar (Integer nx, NagMeanOrTrend mt_correction, double px, Integer mw, double pw, Integer l, Integer kc, Nag_LoggedSpectra lg_spect, const double x[], double *g, Integer ng, double stats[], NagError *fail)

3 Description

The supplied time series may be mean or trend corrected (by least squares), and tapered, the tapering factors being those of the split cosine bell:

\begin{matrix} \frac{1}{2} (1 - \cos (\frac{π (t - \frac{1}{2})}{T})), & 1 \leq t \leq T \\ \frac{1}{2} (1 - \cos (\frac{π (n - t + \frac{1}{2})}{T})), & n + 1 - T \leq t \leq n \\ 1, & otherwise \end{matrix}

where

T = [\frac{n p}{2}]

and

p

is the tapering proportion.

The unsmoothed sample spectrum

f^{*} (ω) = \frac{1}{2} π {|\sum_{t = 1}^{n} x_{t} \exp (i ω t)|}^{2}

is then calculated for frequency values

ω_{k} = \frac{2 π k}{K}, k = 0, 1, \dots, [K / 2]

where [ ] denotes the integer part.

The smoothed spectrum is returned as a subset of these frequencies for which

K

is a multiple of a chosen value

r

, i.e.,

ω_{r l} = ν_{l} = \frac{2 π l}{L}, l = 0, 1, \dots, [L / 2]

where

K = r \times L

. You will normally fix

L

first, then choose

r

so that

K

is sufficiently large to provide an adequate representation for the unsmoothed spectrum, i.e.,

K \geq 2 \times n

. It is possible to take

L = K

, i.e.,

r = 1

The smoothing is defined by a trapezium window whose shape is supplied by the function

\begin{matrix} W (α) = 1, & |α| \leq p \\ W (α) = \frac{1 - |α|}{1 - p}, & p < |α| \leq 1 \end{matrix}

the proportion

p

being supplied by you.

The width of the window is fixed as

2 π / M

by supplying

M

. A set of averaging weights are constructed:

W_{k} = g \times W (\frac{ω_{k} M}{π}), 0 \leq ω_{k} \leq \frac{π}{M}

where

g

is a normalizing constant, and the smoothed spectrum obtained is

\hat{f} (ν_{l}) = \sum_{|ω_{k}| < \frac{π}{M}} W_{k} f^{*} (ν_{l} + ω_{k}) .

If no smoothing is required

M

should be set to

n

, in which case the values returned are

\hat{f} (ν_{l}) = f^{*} (ν_{l})

. Otherwise, in order that the smoothing approximates well to an integration, it is essential that

K ≫ M

, and preferable, but not essential, that

K

be a multiple of

M

. A choice of

L > M

would normally be required to supply an adequate description of the smoothed spectrum. Typical choices of

L ≃ n

and

K ≃ 4 n

should be adequate for usual smoothing situations when

M < n / 5

The sampling distribution of

\hat{f} (ω)

is approximately that of a scaled

χ_{d}^{2}

variate, whose degrees of freedom

d

is provided by the function, together with multiplying limits

m u

m l

from which approximate 95% confidence intervals for the true spectrum

f (ω)

may be constructed as

[m l \times \hat{f}, (ω), m u, \times, \hat{f}, (ω)]

. Alternatively, log

\hat{f} (ω)

may be returned, with additive limits.

The bandwidth

b

of the corresponding smoothing window in the frequency domain is also provided. Spectrum estimates separated by (angular) frequencies much greater than

b

may be assumed to be independent.

4 References

Bloomfield P (1976) Fourier Analysis of Time Series: An Introduction Wiley

Jenkins G M and Watts D G (1968) Spectral Analysis and its Applications Holden–Day

5 Arguments

1: nx – IntegerInput

On entry: the length of the time series,

n

Constraint:

nx \geq 1

2: mt_correction – NagMeanOrTrendInput

On entry: whether the data are to be initially mean or trend corrected.

mt_correction = Nag_NoCorrection

for no correction,

mt_correction = Nag_Mean

for mean correction,

mt_correction = Nag_Trend

for trend correction.

Constraint:

mt_correction = Nag_NoCorrection

Nag_Mean

Nag_Trend

3: px – doubleInput

On entry: the proportion of the data (totalled over both ends) to be initially tapered by the split cosine bell taper. (A value of 0.0 implies no tapering).

Constraint:

0.0 \leq px \leq 1.0

4: mw – IntegerInput

On entry: the value of

M

which determines the frequency width of the smoothing window as

2 π / M

. A value of

n

implies no smoothing is to be carried out.

Constraint:

1 \leq mw \leq nx

5: pw – doubleInput

On entry: the shape argument,

p

, of the trapezium frequency window.

A value of 0.0 gives a triangular window, and a value of 1.0 a rectangular window.

mw = nx

(i.e., no smoothing is carried out), then pw is not used.

Constraint:

0.0 \leq pw \leq 1.0

. if

mw \neq nx

6: l – IntegerInput

On entry: the frequency division,

L

, of smoothed spectral estimates as

2 π / L

Constraints:

$l \geq 1$ ;
l must be a factor of kc (see below).

7: kc – IntegerInput

On entry: the order of the fast Fourier transform (FFT),

K

, used to calculate the spectral estimates. kc should be a multiple of small primes such as

2^{m}

where

m

is the smallest integer such that

2^{m} \geq 2 n

, provided

m \leq 20

Constraints:

$kc \geq 2 \times nx$ ;
kc must be a multiple of l. The largest prime factor of kc must not exceed $19$ , and the total number of prime factors of kc, counting repetitions, must not exceed $20$ . These two restrictions are imposed by the internal FFT algorithm used.

8: lg_spect – Nag_LoggedSpectraInput

On entry: indicates whether unlogged or logged spectral estimates and confidence limits are required.

lg_spect = Nag_Unlogged

for unlogged.

lg_spect = Nag_Logged

for logged.

Constraint:

lg_spect = Nag_Unlogged

Nag_Logged

9: x[kc] – const doubleInput

On entry: the

n

data points.

10: g – double **Output

On exit: vector which contains the ng spectral estimates

\hat{f} (ω_{i})

, for

i = 0, 1, \dots, [L / 2]

, in

g [0]

g [ng - 1]

(logged if

lg_spect = Nag_Logged

). The memory for this vector is allocated internally. If no memory is allocated to g (e.g., when an input error is detected) then g will be NULL on return. If repeated calls to this function are required then NAG_FREE should be used to free the memory in between calls.

11: ng – Integer *Output

On exit: the number of spectral estimates,

[L / 2] + 1

, in g.

12: stats[ $4$ ] – doubleOutput

On exit: four associated statistics. These are the degrees of freedom in

stats [0]

, the lower and upper 95% confidence limit factors in

stats [1]

and

stats [2]

respectively (logged if

lg_spect = Nag_Logged

), and the bandwidth in

stats [3]

13: fail – NagError *Input/Output

The NAG error argument (see Section 3.6 in the Essential Introduction).

6 Error Indicators and Warnings

NE_2_INT_ARG_CONS: On entry, $kc = ⟨value⟩$ while $l = ⟨value⟩$ . These arguments must satisfy kc% $l = 0$ when $l > 0$ .
On entry, $kc = ⟨value⟩$ while $nx = ⟨value⟩$ . These arguments must satisfy $kc \geq 2 \times nx$ when $nx > 0$ .
NE_2_INT_ARG_GT: On entry, $mw = ⟨value⟩$ while $nx = ⟨value⟩$ . These arguments must satisfy $mw \leq nx$ .
NE_ALLOC_FAIL: Dynamic memory allocation failed.
NE_BAD_PARAM: On entry, argument lg_spect had an illegal value.
On entry, argument mt_correction had an illegal value.
NE_CONFID_LIMIT_FACT: The calculation of confidence limit factors has failed. Spectral estimates (logged if requested) are returned in g, and degrees of freedom and bandwith in stats.
NE_FACTOR_GT: At least one of the prime factors of kc is greater than 19.
NE_INT_ARG_LT: On entry, $l = ⟨value⟩$ .
Constraint: $l \geq 1$ .
On entry, $mw = ⟨value⟩$ .
Constraint: $mw \geq 1$ .
On entry, $nx = ⟨value⟩$ .
Constraint: $nx \geq 1$ .
On entry, pw must not be less than 0.0: $pw = ⟨value⟩$ .
NE_INTERNAL_ERROR: An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
NE_REAL_ARG_GT: On entry, pw must not be greater than 1.0: $pw = ⟨value⟩$ .
On entry, px must not be greater than 1.0: $px = ⟨value⟩$ .
NE_REAL_ARG_LT: On entry, px must not be less than 0.0: $px = ⟨value⟩$ .
NE_SPECTRAL_ESTIM_NEG: One or more spectral estimates are negative. Unlogged spectral estimates are returned in g and the degrees of freedom, unlogged confidence limit factors and bandwith in stats.
NE_TOO_MANY_FACTORS: kc has more than 20 prime factors.

7 Accuracy

The FFT is a numerically stable process, and any errors introduced during the computation will normally be insignificant compared with uncertainty in the data.

8 Parallelism and Performance

Not applicable.

9 Further Comments

nag_tsa_spectrum_univar (g13cbc) carries out a FFT of length kc to calculate the sample spectrum. The time taken by the function for this is approximately proportional to

kc \times \log (kc)

(but see Section 9 in nag_sum_fft_realherm_1d (c06pac) for further details).

10 Example

The example program reads a time series of length 131. It selects the mean correction option, a tapering proportion of 0.2, the option of no smoothing and a frequency division for logged spectral estimates of

2 π / 100

. It then calls nag_tsa_spectrum_univar (g13cbc) to calculate the univariate spectrum and prints the logged spectrum together with 95% confidence limits. The program then selects a smoothing window with frequency width

2 π / 30

and shape argument 0.5 and recalculates and prints the logged spectrum and 95% confidence limits.

NAG Library Function Documentnag_tsa_spectrum_univar (g13cbc)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG Library Function Document

nag_tsa_spectrum_univar (g13cbc)