NAG CL Interface
g13cbc (uni_spectrum_daniell)
1
Purpose
g13cbc calculates the smoothed sample spectrum of a univariate time series using spectral smoothing by the trapezium frequency (Daniell) window.
2
Specification
void 
g13cbc (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) 

The function may be called by the names: g13cbc, nag_tsa_uni_spectrum_daniell or nag_tsa_spectrum_univar.
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:
where
$T=\left[\frac{np}{2}\right]$and
$p$ is the tapering proportion.
The unsmoothed sample spectrum
is then calculated for frequency values
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.,
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\ge 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
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:
where
$g$ is a normalizing constant, and the smoothed spectrum obtained is
If no smoothing is required
$M$ should be set to
$n$, in which case the values returned are
$\hat{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<n/5$.
The sampling distribution of $\hat{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\times \hat{f},\left(\omega \right),mu,\times ,\hat{f},\left(\omega \right)\right]$. Alternatively, log $\hat{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.
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:
$\mathbf{nx}$ – Integer
Input

On entry: the length of the time series, $n$.
Constraint:
${\mathbf{nx}}\ge 1$.

2:
$\mathbf{mt\_correction}$ – NagMeanOrTrend
Input

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}$ – double
Input

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}$ – Integer
Input

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}$ – double
Input

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}$ – Integer
Input

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}$ – Integer
Input

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\times {\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_LoggedSpectra
Input

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 double
Input

On entry: the $n$ data points.

10:
$\mathbf{g}$ – double **
Output

On exit: vector which contains the
ng spectral estimates
$\hat{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]$ – double
Output

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 7 in the Introduction to the NAG Library CL Interface).
6
Error Indicators and Warnings
 NE_2_INT_ARG_CONS

On entry,
${\mathbf{kc}}=\u2329\mathit{\text{value}}\u232a$ while
${\mathbf{l}}=\u2329\mathit{\text{value}}\u232a$. These arguments must satisfy
kc%
${\mathbf{l}}=0$ when
${\mathbf{l}}>0$.
On entry, ${\mathbf{kc}}=\u2329\mathit{\text{value}}\u232a$ while ${\mathbf{nx}}=\u2329\mathit{\text{value}}\u232a$. These arguments must satisfy ${\mathbf{kc}}\ge 2\times {\mathbf{nx}}$ when ${\mathbf{nx}}>0$.
 NE_2_INT_ARG_GT

On entry, ${\mathbf{mw}}=\u2329\mathit{\text{value}}\u232a$ while ${\mathbf{nx}}=\u2329\mathit{\text{value}}\u232a$. These arguments must satisfy ${\mathbf{mw}}\le {\mathbf{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, ${\mathbf{l}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{l}}\ge 1$.
On entry, ${\mathbf{mw}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{mw}}\ge 1$.
On entry, ${\mathbf{nx}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{nx}}\ge 1$.
On entry,
pw must not be less than 0.0:
${\mathbf{pw}}=\u2329\mathit{\text{value}}\u232a$.
 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}}=\u2329\mathit{\text{value}}\u232a$.
On entry,
px must not be greater than 1.0:
${\mathbf{px}}=\u2329\mathit{\text{value}}\u232a$.
 NE_REAL_ARG_LT

On entry,
px must not be less than 0.0:
${\mathbf{px}}=\u2329\mathit{\text{value}}\u232a$.
 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
g13cbc is not threaded in any implementation.
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}}\times \mathrm{log}\left({\mathbf{kc}}\right)$ (but see
Section 9 in
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\pi /100$. It then calls 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.1
Program Text
10.2
Program Data
10.3
Program Results