Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_tsa_uni_spectrum_daniell (g13cb)

## Purpose

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

## Syntax

[xg, ng, stats, ifail] = g13cb(nx, mtx, px, mw, pw, l, lg, xg, 'kc', kc)
[xg, ng, stats, ifail] = nag_tsa_uni_spectrum_daniell(nx, mtx, px, mw, pw, l, lg, xg, 'kc', kc)

## 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:
 $12 1-cosπ t-12 /T , 1≤t≤T 12 1-cosπ n-t+12 /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*ω=12π ∑t=1nxtexpiω t 2$
is then calculated for frequency values
 $ωk=2 π kK, 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πlL, 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 you supplying $M$. A set of averaging weights are constructed:
 $Wk=g×W ωkM π , 0≤ωk≤πM,$
where $g$ is a normalizing constant, and the smoothed spectrum obtained is
 $f^νl=∑ωk< πMWkf*ν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.

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

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{nx}$int64int32nag_int scalar
$n$, the length of the time series.
Constraint: ${\mathbf{nx}}\ge 1$.
2:     $\mathrm{mtx}$int64int32nag_int scalar
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$.
3:     $\mathrm{px}$ – double scalar
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:     $\mathrm{mw}$int64int32nag_int scalar
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:     $\mathrm{pw}$ – double scalar
$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$.
6:     $\mathrm{l}$int64int32nag_int scalar
$L$, the frequency division of smoothed spectral estimates as $2\pi /L$.
Constraints:
• ${\mathbf{l}}\ge 1$;
• l must be a factor of kc.
7:     $\mathrm{lg}$int64int32nag_int scalar
Indicates whether unlogged or logged spectral estimates and confidence limits are required.
${\mathbf{lg}}=0$
For unlogged.
${\mathbf{lg}}\ne 0$
For logged.
8:     $\mathrm{xg}\left({\mathbf{kc}}\right)$ – double array
The $n$ data points.

### Optional Input Parameters

1:     $\mathrm{kc}$int64int32nag_int scalar
Default: the dimension of the array xg.
$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$.
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.

### Output Parameters

1:     $\mathrm{xg}\left({\mathbf{kc}}\right)$ – double array
Contains the ng spectral estimates $\stackrel{^}{f}\left({\omega }_{\mathit{i}}\right)$, for $\mathit{i}=0,1,\dots ,\left[L/2\right]$, in ${\mathbf{xg}}\left(1\right)$ to ${\mathbf{xg}}\left({\mathbf{ng}}\right)$ (logged if ${\mathbf{lg}}\ne 0$). The elements ${\mathbf{xg}}\left(\mathit{i}\right)$, for $\mathit{i}={\mathbf{ng}}+1,\dots ,{\mathbf{kc}}$, contain $0.0$.
2:     $\mathrm{ng}$int64int32nag_int scalar
The number of spectral estimates, $\left[L/2\right]+1$, in xg.
3:     $\mathrm{stats}\left(4\right)$ – double array
Four associated statistics. These are the degrees of freedom in ${\mathbf{stats}}\left(1\right)$, the lower and upper $95%$ confidence limit factors in ${\mathbf{stats}}\left(2\right)$ and ${\mathbf{stats}}\left(3\right)$ respectively (logged if ${\mathbf{lg}}\ne 0$), and the bandwidth in ${\mathbf{stats}}\left(4\right)$.
4:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{0}}$ unless the function detects an error (see Error Indicators and Warnings).

## Error Indicators and Warnings

Note: nag_tsa_uni_spectrum_daniell (g13cb) may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

${\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×{\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 function argument lists in calls to nag_tsa_uni_spectrum_daniell (g13cb). Seek expert help.
W  ${\mathbf{ifail}}=4$
One or more spectral estimates are negative. Unlogged spectral estimates are returned in xg, and the degrees of freedom, unlogged confidence limit factors and bandwidth in stats.
W  ${\mathbf{ifail}}=5$
The calculation of confidence limit factors has failed. This error will not normally occur. Spectral estimates (logged if requested) are returned in xg, and degrees of freedom and bandwidth in stats.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

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

nag_tsa_uni_spectrum_daniell (g13cb) 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 Further Comments in nag_sum_fft_realherm_1d (c06pa) for further details).

## Example

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

fprintf('g13cb example results\n\n');

% Smoothing parameters
nx  = int64(131);
mtx = int64(1);
px  = 0.2;
mw  = int64(30);
pw = 0.5;
l  = int64(100);
lg = int64(1);

% Series
xg = zeros(400, 1);
xg(1:131) = ...
[11.500;  9.890;  8.728;  8.400;  8.230;  8.365;  8.383;  8.243;  8.080;
8.244;  8.490;  8.867;  9.469;  9.786; 10.100; 10.714; 11.320; 11.900;
12.390; 12.095; 11.800; 12.400; 11.833; 12.200; 12.242; 11.687; 10.883;
10.138;  8.952;  8.443;  8.231;  8.067;  7.871;  7.962;  8.217;  8.689;
8.989;  9.450;  9.883; 10.150; 10.787; 11.000; 11.133; 11.100; 11.800;
12.250; 11.350; 11.575; 11.800; 11.100; 10.300;  9.725;  9.025;  8.048;
7.294;  7.070;  6.933;  7.208;  7.617;  7.867;  8.309;  8.640;  9.179;
9.570; 10.063; 10.803; 11.547; 11.550; 11.800; 12.200; 12.400; 12.367;
12.350; 12.400; 12.270; 12.300; 11.800; 10.794;  9.675;  8.900;  8.208;
8.087;  7.763;  7.917;  8.030;  8.212;  8.669;  9.175;  9.683; 10.290;
10.400; 10.850; 11.700; 11.900; 12.500; 12.500; 12.800; 12.950; 13.050;
12.800; 12.800; 12.800; 12.600; 11.917; 10.805;  9.240;  8.777;  8.683;
8.649;  8.547;  8.625;  8.750;  9.110;  9.392;  9.787; 10.340; 10.500;
11.233; 12.033; 12.200; 12.300; 12.600; 12.800; 12.650; 12.733; 12.700;
12.259; 11.817; 10.767;  9.825;  9.150];

% Calculate smooth spectrum
[xg, ng, stats, ifail] = ...
g13cb( ...
nx, mtx, px, mw, pw, l, lg, xg);

% Display results
fprintf('Frequency width of smoothing window = 1/%d\n',mw);
fprintf('Degrees of freedom = %4.1f      Bandwidth = %7.4f\n\n', ...
stats(1), stats(4));
fprintf('95 percent confidence limits -  Lower = %7.4f  Upper = %7.4f\n', ...
stats(2:3));
fprintf('\n      Spectrum       Spectrum       Spectrum       Spectrum\n');
fprintf('      estimate       estimate       estimate       estimate\n');
result = [double([1:ng]); xg(1:ng)'];
for j = 1:4:ng
fprintf('%4d%10.4f', result(:,j:min(j+3,ng)));
fprintf('\n');
end

```
```g13cb example results

Frequency width of smoothing window = 1/30
Degrees of freedom =  7.0      Bandwidth =  0.1767

95 percent confidence limits -  Lower = -0.8275  Upper =  1.4213

Spectrum       Spectrum       Spectrum       Spectrum
estimate       estimate       estimate       estimate
1   -0.1776   2   -0.4561   3   -0.1784   4    1.9042
5    2.1094   6    1.7061   7   -0.7659   8   -1.4734
9   -1.5939  10   -2.1157  11   -2.9151  12   -2.7055
13   -2.8200  14   -3.4077  15   -3.8813  16   -3.6607
17   -4.0601  18   -4.4756  19   -4.2700  20   -4.3092
21   -4.5711  22   -4.8111  23   -4.5658  24   -4.7285
25   -5.4386  26   -5.5081  27   -5.2325  28   -5.0262
29   -4.4539  30   -4.4764  31   -4.9152  32   -5.8492
33   -5.5872  34   -4.9804  35   -4.8904  36   -5.2666
37   -5.7643  38   -5.8620  39   -5.5011  40   -5.7129
41   -6.3894  42   -6.4027  43   -6.1352  44   -6.5766
45   -7.3676  46   -7.1405  47   -6.1674  48   -5.8600
49   -6.1036  50   -6.2673  51   -6.4321
```