NAG Library Function Document

1Purpose

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

2Specification

 #include #include
 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)

3Description

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:
 $1 2 1 - cos π t - 1 2 T , 1 ≤ t ≤ T 1 2 1 - cos π n - t + 1 2 T , n + 1 - T ≤ t ≤ n 1 , otherwise$
where $T=\left[\frac{np}{2}\right]$and $p$ is the tapering proportion.
The unsmoothed sample spectrum
 $f * ω = 1 2 π ∑ t=1 n x t expiωt 2$
is then calculated for frequency values
 $ω k = 2πk K , k = 0 , 1 , … , 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.,
 $ω rl = ν l = 2πl L , l = 0 , 1 , … , L/2$
where $K=r×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\ge 2×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
 $W α = 1 , α ≤ p W α = 1 - α 1-p , p < α ≤ 1$
the proportion $p$ being supplied by you.
The width of the window is fixed as $2\pi /M$ by supplying $M$. A set of averaging weights are constructed:
 $W k = g × W ω k M π , 0 ≤ ω k ≤ π M$
where $g$ is a normalizing constant, and the smoothed spectrum obtained is
 $f ^ ν l = ∑ ω k < π M W k f * ν l + ω k .$
If no smoothing is required $M$ should be set to $n$, in which case the values returned are $\stackrel{^}{f}\left({\nu }_{l}\right)={f}^{*}\left({\nu }_{l}\right)$. Otherwise, in order that the smoothing approximates well to an integration, it is essential that $K\gg 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\simeq n$ and $K\simeq 4n$ should be adequate for usual smoothing situations when $M.
The sampling distribution of $\stackrel{^}{f}\left(\omega \right)$ is approximately that of a scaled ${\chi }_{d}^{2}$ variate, whose degrees of freedom $d$ is provided by the function, together with multiplying limits $mu$, $ml$ from which approximate 95% confidence intervals for the true spectrum $f\left(\omega \right)$ may be constructed as $\left[ml×\stackrel{^}{f},\left(\omega \right),mu,×,\stackrel{^}{f},\left(\omega \right)\right]$. Alternatively, log $\stackrel{^}{f}\left(\omega \right)$ 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.
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

5Arguments

1:    $\mathbf{nx}$IntegerInput
On entry: the length of the time series, $n$.
Constraint: ${\mathbf{nx}}\ge 1$.
2:    $\mathbf{mt_correction}$NagMeanOrTrendInput
On entry: whether the data are to be initially mean or trend corrected. ${\mathbf{mt_correction}}=\mathrm{Nag_NoCorrection}$ for no correction, ${\mathbf{mt_correction}}=\mathrm{Nag_Mean}$ for mean correction, ${\mathbf{mt_correction}}=\mathrm{Nag_Trend}$ for trend correction.
Constraint: ${\mathbf{mt_correction}}=\mathrm{Nag_NoCorrection}$, $\mathrm{Nag_Mean}$ or $\mathrm{Nag_Trend}$.
3:    $\mathbf{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\le {\mathbf{px}}\le 1.0$.
4:    $\mathbf{mw}$IntegerInput
On entry: the value of $M$ which determines the frequency width of the smoothing window as $2\pi /M$. A value of $n$ implies no smoothing is to be carried out.
Constraint: $1\le {\mathbf{mw}}\le {\mathbf{nx}}$.
5:    $\mathbf{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.
If ${\mathbf{mw}}={\mathbf{nx}}$ (i.e., no smoothing is carried out), then pw is not used.
Constraint: $0.0\le {\mathbf{pw}}\le 1.0$. if ${\mathbf{mw}}\ne {\mathbf{nx}}$.
6:    $\mathbf{l}$IntegerInput
On entry: the frequency division, $L$, of smoothed spectral estimates as $2\pi /L$.
Constraints:
• ${\mathbf{l}}\ge 1$;
• l must be a factor of kc (see below).
7:    $\mathbf{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}\ge 2n$, provided $m\le 20$.
Constraints:
• ${\mathbf{kc}}\ge 2×{\mathbf{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:    $\mathbf{lg_spect}$Nag_LoggedSpectraInput
On entry: indicates whether unlogged or logged spectral estimates and confidence limits are required. ${\mathbf{lg_spect}}=\mathrm{Nag_Unlogged}$ for unlogged. ${\mathbf{lg_spect}}=\mathrm{Nag_Logged}$ for logged.
Constraint: ${\mathbf{lg_spect}}=\mathrm{Nag_Unlogged}$ or $\mathrm{Nag_Logged}$.
9:    $\mathbf{x}\left[{\mathbf{kc}}\right]$const doubleInput
On entry: the $n$ data points.
10:  $\mathbf{g}$double **Output
On exit: vector which contains the ng spectral estimates $\stackrel{^}{f}\left({\omega }_{\mathit{i}}\right)$, for $\mathit{i}=0,1,\dots ,\left[L/2\right]$, in ${\mathbf{g}}\left[0\right]$ to ${\mathbf{g}}\left[{\mathbf{ng}}-1\right]$ (logged if ${\mathbf{lg_spect}}=\mathrm{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:  $\mathbf{ng}$Integer *Output
On exit: the number of spectral estimates, $\left[L/2\right]+1$, in g.
12:  $\mathbf{stats}\left[4\right]$doubleOutput
On exit: four associated statistics. These are the degrees of freedom in ${\mathbf{stats}}\left[0\right]$, the lower and upper 95% confidence limit factors in ${\mathbf{stats}}\left[1\right]$ and ${\mathbf{stats}}\left[2\right]$ respectively (logged if ${\mathbf{lg_spect}}=\mathrm{Nag_Logged}$), and the bandwidth in ${\mathbf{stats}}\left[3\right]$.
13:  $\mathbf{fail}$NagError *Input/Output
The NAG error argument (see Section 3.7 in How to Use the NAG Library and its Documentation).

6Error Indicators and Warnings

NE_2_INT_ARG_CONS
On entry, ${\mathbf{kc}}=〈\mathit{\text{value}}〉$ while ${\mathbf{l}}=〈\mathit{\text{value}}〉$. These arguments must satisfy kc%${\mathbf{l}}=0$ when ${\mathbf{l}}>0$.
On entry, ${\mathbf{kc}}=〈\mathit{\text{value}}〉$ while ${\mathbf{nx}}=〈\mathit{\text{value}}〉$. These arguments must satisfy ${\mathbf{kc}}\ge 2×{\mathbf{nx}}$ when ${\mathbf{nx}}>0$.
NE_2_INT_ARG_GT
On entry, ${\mathbf{mw}}=〈\mathit{\text{value}}〉$ while ${\mathbf{nx}}=〈\mathit{\text{value}}〉$. These arguments must satisfy ${\mathbf{mw}}\le {\mathbf{nx}}$.
NE_ALLOC_FAIL
Dynamic memory allocation failed.
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, ${\mathbf{l}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{l}}\ge 1$.
On entry, ${\mathbf{mw}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{mw}}\ge 1$.
On entry, ${\mathbf{nx}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{nx}}\ge 1$.
On entry, pw must not be less than 0.0: ${\mathbf{pw}}=〈\mathit{\text{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: ${\mathbf{pw}}=〈\mathit{\text{value}}〉$.
On entry, px must not be greater than 1.0: ${\mathbf{px}}=〈\mathit{\text{value}}〉$.
NE_REAL_ARG_LT
On entry, px must not be less than 0.0: ${\mathbf{px}}=〈\mathit{\text{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.

7Accuracy

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

8Parallelism and Performance

nag_tsa_spectrum_univar (g13cbc) is not threaded in any implementation.

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 ${\mathbf{kc}}×\mathrm{log}\left({\mathbf{kc}}\right)$ (but see Section 9 in nag_sum_fft_realherm_1d (c06pac) for further details).

10Example

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\pi /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\pi /30$ and shape argument 0.5 and recalculates and prints the logged spectrum and 95% confidence limits.

10.1Program Text

Program Text (g13cbce.c)

10.2Program Data

Program Data (g13cbce.d)

10.3Program Results

Program Results (g13cbce.r)

© The Numerical Algorithms Group Ltd, Oxford, UK. 2017