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

# Syntax

C# |
---|

public static void g13cb( int nx, int mtx, double px, int mw, double pw, int l, int kc, int lg, double[] xg, out int ng, double[] stats, out int ifail ) |

Visual Basic |
---|

Public Shared Sub g13cb ( _ nx As Integer, _ mtx As Integer, _ px As Double, _ mw As Integer, _ pw As Double, _ l As Integer, _ kc 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 g13cb( int nx, int mtx, double px, int mw, double pw, int l, int kc, int lg, array<double>^ xg, [OutAttribute] int% ng, array<double>^ stats, [OutAttribute] int% ifail ) |

F# |
---|

static member g13cb : nx : int * mtx : int * px : float * mw : int * pw : float * l : int * kc : 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*: 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*: $0\le {\mathbf{mtx}}\le 2$.

- px
- Type: System..::..Double
*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$.

- mw
- Type: System..::..Int32
*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}}$.

- pw
- Type: System..::..Double
*On entry*: $p$, the shape parameter 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), pw is not used.*Constraint*: $0.0\le {\mathbf{pw}}\le 1.0$.

- l
- Type: System..::..Int32
*On entry*: $L$, the frequency division of smoothed spectral estimates as $2\pi /L$.

- kc
- Type: System..::..Int32
*On entry*: $K$, the order of the fast Fourier transform (FFT) 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$.

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

- xg
- Type: array<System..::..Double>[]()[][]An array of size [kc]
*On entry*: the $n$ data points.*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]$ (logged if ${\mathbf{lg}}\ne 0$). The elements ${\mathbf{xg}}\left[\mathit{i}-1\right]$, for $\mathit{i}={\mathbf{ng}}+1,\dots ,{\mathbf{kc}}$, 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}}\ne 0$), 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 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.

$$\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 unsmoothed sample spectrum

is then calculated for frequency values

where [ ] denotes the integer part.

$${f}^{*}\left(\omega \right)=\frac{1}{2\pi}{\left|\sum _{t=1}^{n}{x}_{t}\mathrm{exp}\left(i\omega t\right)\right|}^{2}$$ |

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

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

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

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

the proportion $p$ being supplied by you.

$$\begin{array}{ll}W\left(\alpha \right)=1\text{,}& \left|\alpha \right|\le p\\ W\left(\alpha \right)=\frac{1-\left|\alpha \right|}{1-p}\text{,}& p<\left|\alpha \right|\le 1\end{array}$$ |

The width of the window is fixed as $2\pi /M$ by you 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$.

$${W}_{k}=g\times W\left(\frac{{\omega}_{k}M}{\pi}\right)\text{, \hspace{1em}}0\le {\omega}_{k}\le \frac{\pi}{M}\text{,}$$ |

$$\hat{f}\left({\nu}_{l}\right)={\sum}_{\left|{\omega}_{k}\right|<\frac{\pi}{M}}{W}_{k}{f}^{*}\left({\nu}_{l}+{\omega}_{k}\right)\text{.}$$ |

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.

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

**Note:**g13cb may return useful information for one or more of the following detected errors or warnings.

Errors or warnings detected by the method:

- ${\mathbf{ifail}}=1$
On entry, ${\mathbf{nx}}<1$, or ${\mathbf{mtx}}<0$, or ${\mathbf{mtx}}>2$, or ${\mathbf{px}}<0.0$, or ${\mathbf{px}}>1.0$, or ${\mathbf{mw}}<1$, or ${\mathbf{mw}}>{\mathbf{nx}}$, or ${\mathbf{pw}}<0.0$ and ${\mathbf{mw}}\ne {\mathbf{nx}}$, or ${\mathbf{pw}}>1.0$ and ${\mathbf{mw}}\ne {\mathbf{nx}}$, or ${\mathbf{l}}<1$.

- ${\mathbf{ifail}}=2$
On entry, ${\mathbf{kc}}<2\times {\mathbf{nx}}$, or kc is not a multiple of l, or kc has a prime factor exceeding $19$, or kc has more than $20$ prime factors, counting repetitions.

- ${\mathbf{ifail}}=3$
- This indicates that a serious error has occurred. Check all array subscripts and method parameter lists in calls to g13cb. Seek expert help.

- ${\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

g13cb carries out a FFT of length kc to calculate the sample spectrum. The time taken by the method for this is approximately proportional to ${\mathbf{kc}}\times \mathrm{log}\left({\mathbf{kc}}\right)$ (but see [Further Comments] in c06pa for further details).

# Example

This example reads a time series of length $131$. It then calls g13cb to calculate the univariate spectrum and prints the logged spectrum together with $95\%$ confidence limits.

Example program (C#): g13cbe.cs