g13ca calculates the smoothed sample spectrum of a univariate time series using one of four lag windows – rectangular, Bartlett, Tukey or Parzen window.

# Syntax

C# |
---|

public static void g13ca( int nx, int mtx, double px, int iw, int mw, int ic, int nc, double[] c, int kc, int l, int lg, double[] xg, out int ng, double[] stats, out int ifail ) |

Visual Basic |
---|

Public Shared Sub g13ca ( _ nx As Integer, _ mtx As Integer, _ px As Double, _ iw As Integer, _ mw As Integer, _ ic As Integer, _ nc As Integer, _ c As Double(), _ kc As Integer, _ l As Integer, _ lg As Integer, _ xg As Double(), _ <OutAttribute> ByRef ng As Integer, _ stats As Double(), _ <OutAttribute> ByRef ifail As Integer _ ) |

Visual C++ |
---|

public: static void g13ca( int nx, int mtx, double px, int iw, int mw, int ic, int nc, array<double>^ c, int kc, int l, int lg, array<double>^ xg, [OutAttribute] int% ng, array<double>^ stats, [OutAttribute] int% ifail ) |

F# |
---|

static member g13ca : nx : int * mtx : int * px : float * iw : int * mw : int * ic : int * nc : int * c : float[] * kc : int * l : int * lg : int * xg : float[] * ng : int byref * stats : float[] * ifail : int byref -> unit |

#### Parameters

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

- mtx
- Type: System..::..Int32
*On entry*: if covariances are to be calculated by the method (${\mathbf{ic}}=0$), mtx must specify whether the data are to be initially mean or trend corrected.- ${\mathbf{mtx}}=0$
- For no correction.
- ${\mathbf{mtx}}=1$
- For mean correction.
- ${\mathbf{mtx}}=2$
- For trend correction.

*Constraint*: if ${\mathbf{ic}}=0$, $0\le {\mathbf{mtx}}\le 2$If covariances are supplied (${\mathbf{ic}}\ne 0$), mtx is not used.

- px
- Type: System..::..Double
*On entry*: if covariances are to be calculated by the method (${\mathbf{ic}}=0$), px must specify the proportion of the data (totalled over both ends) to be initially tapered by the split cosine bell taper.*Constraint*: $0.0\le {\mathbf{px}}\le 1.0$.

- iw
- Type: System..::..Int32
*On entry*: the choice of lag window.- ${\mathbf{iw}}=1$
- Rectangular.
- ${\mathbf{iw}}=2$
- Bartlett.
- ${\mathbf{iw}}=3$
- Tukey.
- ${\mathbf{iw}}=4$
- Parzen.

*Constraint*: $1\le {\mathbf{iw}}\le 4$.

- mw
- Type: System..::..Int32
*On entry*: $M$, the ‘cut-off’ point of the lag window. Windowed covariances at lag $M$ or greater are zero.*Constraint*: $1\le {\mathbf{mw}}\le {\mathbf{nx}}$.

- ic
- Type: System..::..Int32
*On entry*: indicates whether covariances are to be calculated in the method or supplied in the call to the method.- ${\mathbf{ic}}=0$
- Covariances are to be calculated.
- ${\mathbf{ic}}\ne 0$
- Covariances are to be supplied.

- nc
- Type: System..::..Int32
*On entry*: the number of covariances to be calculated in the method or supplied in the call to the method.*Constraint*: ${\mathbf{mw}}\le {\mathbf{nc}}\le {\mathbf{nx}}$.

- c
- Type: array<System..::..Double>[]()[][]An array of size [nc]

- kc
- Type: System..::..Int32
*On entry*: if ${\mathbf{ic}}=0$, kc must specify the order of the fast Fourier transform (FFT) used to calculate the covariances. kc should be a product of small primes such as ${2}^{m}$ where $m$ is the smallest integer such that ${2}^{m}\ge {\mathbf{nx}}+{\mathbf{nc}}$, provided $m\le 20$.If ${\mathbf{ic}}\ne 0$, that is covariances are supplied, kc is not used.

- l
- Type: System..::..Int32
*On entry*: $L$, the frequency division of the spectral estimates as $\frac{2\pi}{L}$. Therefore it is also the order of the FFT used to construct the sample spectrum from the covariances. l should be a product of small primes such as ${2}^{m}$ where $m$ is the smallest integer such that ${2}^{m}\ge 2M-1$, provided $m\le 20$.

- lg
- Type: System..::..Int32
*On entry*: indicates whether unlogged or logged spectral estimates and confidence limits are required.- ${\mathbf{lg}}=0$
- Unlogged.
- ${\mathbf{lg}}\ne 0$
- Logged.

- xg
- Type: array<System..::..Double>[]()[][]An array of size [nxg]
*On entry*: if the covariances are to be calculated, then xg must contain the nx data points. If covariances are supplied, xg may contain any values.*On exit*: contains the ng spectral estimates, $\hat{f}\left({\omega}_{\mathit{i}}\right)$, for $\mathit{i}=0,1,\dots ,\left[L/2\right]$ in ${\mathbf{xg}}\left[0\right]$ to ${\mathbf{xg}}\left[{\mathbf{ng}}-1\right]$ respectively (logged if ${\mathbf{lg}}=1$). The elements ${\mathbf{xg}}\left[\mathit{i}-1\right]$, for $\mathit{i}={\mathbf{ng}}+1,\dots ,{\mathbf{nxg}}$ contain $0.0$.

- ng
- Type: System..::..Int32%

- stats
- Type: array<System..::..Double>[]()[][]An array of size [$4$]
*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}}=1$), and the bandwidth in ${\mathbf{stats}}\left[3\right]$.

- ifail
- Type: System..::..Int32%
*On exit*: ${\mathbf{ifail}}={0}$ unless the method detects an error or a warning has been flagged (see [Error Indicators and Warnings]).

# Description

The smoothed sample spectrum is defined as

where $M$ is the window width, and is calculated for frequency values

where $\left[\right]$ denotes the integer part.

$$\hat{f}\left(\omega \right)=\frac{1}{2\pi}\left({C}_{0}+2\sum _{k=1}^{M-1}{w}_{k}{C}_{k}\mathrm{cos}\left(\omega k\right)\right)\text{,}$$ |

$${\omega}_{i}=\frac{2\pi i}{L}\text{, \hspace{1em}}i=0,1,\dots ,\left[L/2\right]\text{,}$$ |

The autocovariances ${C}_{k}$ may be supplied by you, or constructed from a time series ${x}_{1},{x}_{2},\dots ,{x}_{n}$, as

the fast Fourier transform (FFT) being used to carry out the convolution in this formula.

$${C}_{k}=\frac{1}{n}\sum _{t=1}^{n-k}{x}_{t}{x}_{t+k}\text{,}$$ |

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

where $T=\left[\frac{np}{2}\right]$ and $p$ is the tapering proportion.

$$\begin{array}{ll}\frac{1}{2}\left(1-\mathrm{cos}\left(\pi \left(t-\frac{1}{2}\right)/T\right)\right)\text{,}& 1\le t\le T\\ & \\ \frac{1}{2}\left(1-\mathrm{cos}\left(\pi \left(n-t+\frac{1}{2}\right)/T\right)\right)\text{,}& n+1-T\le t\le n\\ & \\ 1\text{,}& \text{otherwise,}\end{array}$$ |

The smoothing window is defined by

which for the various windows is defined over $0\le \alpha <1$ by

$${w}_{k}=W\left(\frac{k}{M}\right)\text{, \hspace{1em}}k\le M-1\text{,}$$ |

rectangular:

Bartlett:

Tukey:

Parzen:

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 method, 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.

$$W\left(\alpha \right)=1$$ |

$$W\left(\alpha \right)=1-\alpha $$ |

$$W\left(\alpha \right)=\frac{1}{2}\left(1+\mathrm{cos}\left(\pi \alpha \right)\right)$$ |

$$\begin{array}{ll}W\left(\alpha \right)=1-6{\alpha}^{2}+6{\alpha}^{3}\text{,}& 0\le \alpha \le \frac{1}{2}\\ & \\ W\left(\alpha \right)=2{\left(1-\alpha \right)}^{3}\text{,}& \frac{1}{2}<\alpha <1\text{.}\end{array}$$ |

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.

# References

Bloomfield P (1976)

*Fourier Analysis of Time Series: An Introduction*WileyJenkins G M and Watts D G (1968)

*Spectral Analysis and its Applications*Holden–Day# Error Indicators and Warnings

Errors or warnings detected by the method:

- ${\mathbf{ifail}}=1$
On entry, ${\mathbf{nx}}<1$, or ${\mathbf{mtx}}<0$ and ${\mathbf{ic}}=0$, or ${\mathbf{mtx}}>2$ and ${\mathbf{ic}}=0$, or ${\mathbf{px}}<0.0$, or ${\mathbf{px}}>1.0$, or ${\mathbf{iw}}<1$, or ${\mathbf{iw}}>4$, or ${\mathbf{mw}}<1$, or ${\mathbf{mw}}>{\mathbf{nx}}$, or ${\mathbf{nc}}<{\mathbf{mw}}$, or ${\mathbf{nc}}>{\mathbf{nx}}$, or ${\mathbf{nxg}}<\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{kc}},{\mathbf{l}}\right)$ and ${\mathbf{ic}}=0$, or ${\mathbf{nxg}}<{\mathbf{l}}$ and ${\mathbf{ic}}\ne 0$.

- ${\mathbf{ifail}}=2$
On entry, ${\mathbf{kc}}<{\mathbf{nx}}+{\mathbf{nc}}$, or kc has a prime factor exceeding $19$, or kc has more than $20$ prime factors, counting repetitions. This error only occurs when ${\mathbf{ic}}=0$.

- ${\mathbf{ifail}}=3$
On entry, ${\mathbf{l}}<2\times {\mathbf{mw}}-1$, or l has a prime factor exceeding $19$, or l has more than $20$ prime factors, counting repetitions.

- ${\mathbf{ifail}}=4$

- ${\mathbf{ifail}}=5$

# 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.

# Parallelism and Performance

None.

# Further Comments

g13ca carries out two FFTs of length kc to calculate the covariances and one FFT of length l to calculate the sample spectrum. The time taken by the method for an FFT of length $n$ is approximately proportional to $n\mathrm{log}\left(n\right)$ (but see [Further Comments] in c06pa for further details).

# Example

This example reads a time series of length $256$. It selects the mean correction option, a tapering proportion of $0.1$, the Parzen smoothing window and a cut-off point for the window at lag $100$. It chooses to have $100$ auto-covariances calculated and unlogged spectral estimates at a frequency division of $2\pi /200$. It then calls g13ca to calculate the univariate spectrum and statistics and prints the autocovariances and the spectrum together with its $95\%$ confidence multiplying limits.

Example program (C#): g13cae.cs