Module_Time_Series_Procedures¶

This file is part of statpack.

statpack is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

statpack is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You can find a copy of the GNU Lesser General Public License in the statpack/doc directory.

Authors: Pascal Terray (LOCEAN/IPSL, Paris, France)

MODULE EXPORTING SUBROUTINES AND FUNCTIONS FOR TIME SERIES ANALYSIS

LATEST REVISION : 16/09/2020

`subroutine comp_smooth ( x, smooth_factor )`¶

Purpose¶

Smooth a time series.

Arguments¶

X (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, input vector containing size(X) observations which must be smoothed with a smoothing factor of SMOOTH_FACTOR. On exit, the smoothed vector.

SMOOTH_FACTOR (INPUT) integer(i4b)

On entry, the smoothing factor. The smoothing factor must be greater than 0 and less than size(X).

Further Details¶

The input vector is smoothed with a moving average of, approximately, (2 * SMOOTH_FACTOR) + 1 terms.

For further details, see:

Olagnon, M., 1996:

Traitement de donnees numeriques avec Fortran 90. Masson, 264 pages, Chapter 11.1.2, ISBN 2-225-85259-6.

`subroutine comp_smooth ( x, smooth_factor, dimvar )`¶

Purpose¶

Smooth the rows or the columns of a matrix.

Arguments¶

X (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, input matrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables which must be smoothed with a smoothing factor of SMOOTH_FACTOR. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details.

SMOOTH_FACTOR (INPUT) integer(i4b)

On entry, the smoothing factor. The smoothing factor must be greater than 0 and less than size(X,3-DIMVAR).

DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows:

DIMVAR = 1, the input matrix X contains size(X,2) observations on size(X,1) variables and the rows of X will be smoothed.

DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables and the columns of X will be smoothed.

The default is DIMVAR = 1.

Further Details¶

The input matris is smoothed along the specified dimension with a moving average of, approximately, (2 * SMOOTH_FACTOR) + 1 terms.

For further details, see:

Olagnon, M., 1996:

Traitement de donnees numeriques avec Fortran 90. Masson, 264 pages, Chapter 11.1.2, ISBN 2-225-85259-6.

`subroutine comp_smooth ( x, smooth_factor )`¶

Purpose¶

Smooth a tridimensional array along the third dimension.

Arguments¶

X (INPUT/OUTPUT) real(stnd), dimension(:,:,:)

On entry, input tridimensional array containing size(X,3) observations on size(X,1) by size(X,2) variables which must be smoothed with a smoothing factor of SMOOTH_FACTOR.

SMOOTH_FACTOR (INPUT) integer(i4b)

On entry, the smoothing factor. The smoothing factor must be greater than 0 and less than size(X,3).

Further Details¶

The input tridimensional array is smoothed along the third dimension with a moving average of, approximately, (2 * SMOOTH_FACTOR) + 1 terms.

For further details, see:

Olagnon, M., 1996:

Traitement de donnees numeriques avec Fortran 90. Masson, 264 pages, Chapter 11.1.2, ISBN 2-225-85259-6.

`subroutine comp_trend ( y, nt, itdeg, robust, trend, ntjump, maxiter, rw, no, ok )`¶

Purpose¶

COMP_TREND extracts a smoothed component from a time series using a LOESS method. It returns the smoothed component (e.g. the trend) and, optionally, the robustness weights.

Arguments¶

Y (INPUT) real(stnd), dimension(:)

On entry, the time series to be decomposed.

NT (INPUT) integer(i4b)

On entry, the length of the trend smoother. The value of NT should be an odd integer greater than or equal to 3. As NT increases the values of the trend component become smoother.

ITDEG (INPUT) integer(i4b)

On entry, the degree of locally-fitted polynomial in trend smoothing. The value must be 0, 1 or 2.

ROBUST (INPUT) logical(lgl)

On entry, TRUE if robustness iterations are to be used, FALSE otherwise.

Robustness iterations are carried out until convergence of the trend component, with MAXITER iterations maximum. Convergence occurs if the maximum changes in the trend fit is less than 1% of the component’s range after the previous iteration.

TREND (OUTPUT) real(stnd), dimension(:)

On output, the smoothed (e.g. trend) component.

TREND must verify: size(TREND) = size(Y).

NTJUMP (INPUT, OPTIONAL) integer(i4b)

On entry, the skipping value for trend smoothing. By default, NTJUMP is set to NT/10.

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, the maximum number of robustness iterations.

The default is 15. This argument is not used if ROBUST=FALSE.

RW (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On output, final robustness weights. All RW elements are 1 if ROBUST=FALSE .

RW must verify: size(RW) = size(Y).

NO (OUTPUT, OPTIONAL) integer(i4b)

On output, if:

ROBUST=TRUE : the number of robustness iterations. The iterations end if a convergence criterion is met or if the number is MAXITER.

ROBUST=FALSE : NO is set to 0.

OK (OUTPUT, OPTIONAL) logical(lgl)

On output, if:

ROBUST=TRUE : OK is set to TRUE if the convergence criterion is met and to FALSE otherwise.

ROBUST=FALSE : OK is set to TRUE.

Further Details¶

This subroutine is adapted from subroutine STL developped by Cleveland and coworkers at AT&T Bell Laboratories.

This subroutine decomposes a time series into trend and residual components, assuming that the time series has no seasonal cycle or other harmonic components. The algorithm uses LOESS interpolation to smooth the time series and find the trend.

The LOESS smoother for estimating the trend is specified with three parameters: a width (e.g. NT), a degree (e.g. ITDEG) and a jump (e.g. NTJUMP). The width specifies the number of data points that the local interpolation uses to smooth each point, the degree specifies the degree of the local polynomial that is fit to the data, and the jump specifies how many points are skipped between Loess interpolations, with linear interpolation being done between these points.

If the optional ROBUST argument is set to true, the process is iterative and includes robustness iterations that take advandages of the weighted-least-squares underpinnings of LOESS to remove the effects of outliers.

Note that, finally, that this subroutine expects equally spaced data with no missing values.

For further details, see:

Cleveland, R.B., Cleveland, W.S., McRae, J.E., and Terpenning, I.:

STL: A Seasonal-Trend Decomposition Procedure Based on Loess. Statistics Research Report, AT&T Bell Laboratories.

Cleveland, R.B., Cleveland, W.S., McRae, J.E., and Terpenning, I., 1990:

STL: A Seasonal-Trend Decomposition Procedure Based on Loess. J. Official Stat., 6, 3-73.

Crotinger, J., 2017:

Java implementation of Seasonal-Trend-Loess time-series decomposition algorithm. https://github.com/ServiceNow/stl-decomp-4j

`subroutine comp_trend ( y, nt, itdeg, robust, trend, ntjump, maxiter, rw, no, ok )`¶

Purpose¶

COMP_TREND extracts smoothed components from the (time series) columns of a matrix using a LOESS method. It returns the smoothed components (e.g. the trends) and, optionally, the robustness weights.

Arguments¶

Y (INPUT) real(stnd), dimension(:,:)

On entry, the matrix to be decomposed.

NT (INPUT) integer(i4b)

On entry, the length of the trend smoother. The value of NT should be an odd integer greater than or equal to 3. As NT increases the values of the trend component become smoother.

ITDEG (INPUT) integer(i4b)

On entry, the degree of locally-fitted polynomial in trend smoothing. The value must be 0, 1 or 2.

ROBUST (INPUT) logical(lgl)

On entry, TRUE if robustness iterations are to be used, FALSE otherwise.

Robustness iterations are carried out until convergence of the trend component, with MAXITER iterations maximum. Convergence occurs if the maximum changes in the trend fit is less than 1% of the component’s range after the previous iteration.

TREND (OUTPUT) real(stnd), dimension(:,:)

On output, the smoothed (e.g. trend) components.

TREND must verify: size(TREND,1) = size(Y,1) and size(TREND,2) = size(Y,2).

NTJUMP (INPUT, OPTIONAL) integer(i4b)

On entry, the skipping value for trend smoothing.

By default, NTJUMP is set to NT/10.

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, the maximum number of robustness iterations.

The default is 15. This argument is not used if ROBUST=FALSE.

RW (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On output, final robustness weights. All RW elements are 1 if ROBUST=FALSE .

RW must verify: size(RW,1) = size(Y,1) and size(RW,2) = size(Y,2).

NO (OUTPUT, OPTIONAL) integer(i4b), dimension(:)

On output, if

ROBUST=TRUE : NO(i) is the number of robustness iterations for each time series Y(:,i). The iterations end if a convergence criterion is met or if the number is MAXITER for each time series Y(:,i).

ROBUST=FALSE : NO(:) is set to 0.

NO must verify: size(NO) = size(Y,2).

OK (OUTPUT, OPTIONAL) logical(lgl), dimension(:)

On output, if

ROBUST=TRUE : OK(i) is set to TRUE if the convergence criterion is met for time series Y(:,i) and to FALSE otherwise.

ROBUST=FALSE : OK(:) is set to TRUE.

OK must verify: size(OK) = size(Y,2).

Further Details¶

This subroutine is adapted from subroutine STL developped by Cleveland and coworkers at AT&T Bell Laboratories.

This subroutine decomposes a multi-channel time series into trend and residual components, assuming that the multi-channel time series has no seasonal cycle or other harmonic components. The algorithm uses LOESS interpolation to smooth the multi-channel time series and find the trends.

The LOESS smoother for estimating the trends is specified with three parameters: a width (e.g. NT), a degree (e.g. ITDEG) and a jump (e.g. NTJUMP). The width specifies the number of data points that the local interpolation uses to smooth each point, the degree specifies the degree of the local polynomial that is fit to the data, and the jump specifies how many points are skipped between Loess interpolations, with linear interpolation being done between these points.

If the optional ROBUST argument is set to true, the process is iterative and includes robustness iterations that take advandages of the weighted-least-squares underpinnings of LOESS to remove the effects of outliers.

Note that, finally, that this subroutine expects equally spaced data with no missing values.

For further details, see description of COMP_STL and :

Cleveland, R.B., Cleveland, W.S., McRae, J.E., and Terpenning, I.:

STL: A Seasonal-Trend Decomposition Procedure Based on Loess. Statistics Research Report, AT&T Bell Laboratories.

Cleveland, R.B., Cleveland, W.S., McRae, J.E., and Terpenning, I., 1990:

STL: A Seasonal-Trend Decomposition Procedure Based on Loess. J. Official Stat., 6, 3-73.

Crotinger, J., 2017:

Java implementation of Seasonal-Trend-Loess time-series decomposition algorithm. https://github.com/ServiceNow/stl-decomp-4j

`subroutine comp_stlez ( y, np, ns, isdeg, itdeg, robust, season, trend, ni, nt, nl, ildeg, nsjump, ntjump, nljump, maxiter, rw, no, ok )`¶

Purpose¶

COMP_STLEZ decomposes a time series into seasonal and trend components. It returns the components and, optionally, the robustness weights.

COMP_STLEZ offers an easy to use version of COMP_STL subroutine, also included in STATPACK, by defaulting most parameters values associated with the three LOESS smoothers used in COMP_STL.

Arguments¶

Y (INPUT) real(stnd), dimension(:)

On entry, the time series to be decomposed.

NP (INPUT) integer(i4b)

On entry, the period of the seasonal component. For example, if the time series is monthly with a yearly cycle, then NP=12 should be used. NP must be greater than 1.

NS (INPUT) integer(i4b)

On entry, the length of the seasonal smoother. The value of NS should be an odd integer greater than or equal to 3; NS>6 is recommended. As NS increases the values of the seasonal component at a given point in the seasonal cycle (e.g., January values of a monthly series with a yearly cycle) become smoother.

ISDEG (INPUT) integer(i4b)

On entry, the degree of locally-fitted polynomial in seasonal smoothing. The value must be 0 or 1.

ITDEG (INPUT) integer(i4b)

On entry, the degree of locally-fitted polynomial in trend smoothing. The value must be 0, 1 or 2.

ROBUST (INPUT) logical(lgl)

On entry, TRUE if robustness iterations are to be used, FALSE otherwise.

Robustness iterations are carried out until convergence of both seasonal and trend components, with MAXITER iterations maximum. Convergence occurs if the maximum changes in individual seasonal and trend fits are less than 1% of the component’s range after the previous iteration.

SEASON (OUTPUT) real(stnd), dimension(:)

On output, the seasonal component.

SEASON must verify: size(SEASON) = size(Y).

TREND (OUTPUT) real(stnd), dimension(:)

On output, the trend component.

TREND must verify: size(TREND) = size(Y).

NI (INPUT, OPTIONAL) integer(i4b)

On entry, the number of loops for updating the seasonal and trend components. The value of NI should be a positive integer.

By default, NI=2 if ROBUST=FALSE and NI=1 if ROBUST=TRUE.

NT (INPUT, OPTIONAL) integer(i4b)

On entry, the length of the trend smoother. The value of NT should be an odd integer greater than or equal to 3. A value of NT between 1.5 * NP and 2 * NP is recommended. As NT increases the values of the trend component become smoother.

By default, NT is set to the smallest odd integer greater than or equal to (1.5 * NP) / (1-(1.5/NS)).

NL (INPUT, OPTIONAL) integer(i4b)

On entry, the length of the low-pass filter. The value of NL should be an odd integer greater than or equal to 3. The smallest odd integer greater than or equal to NP is recommended.

By default, NL is set to the smallest odd integer greater than or equal to NP.

ILDEG (INPUT, OPTIONAL) integer(i4b)

On entry, the degree of locally-fitted polynomial in low-pass smoothing.

By default, ILDEG is set to ITDEG.

NSJUMP (INPUT, OPTIONAL) integer(i4b)

On entry, the skipping value for seasonal smoothing. The seasonal smoother skips ahead NSJUMP points and then linearly interpolates in between. The value of NSJUMP should be a positive integer; if NSJUMP=1, a seasonal smooth is calculated at all size(Y) points. To make the procedure run faster, a reasonable choice for NSJUMP is 10% or 20% of NS.

By default, NSJUMP is set to NS/10.

NTJUMP (INPUT, OPTIONAL) integer(i4b)

On entry, the skipping value for trend smoothing.

By default, NTJUMP is set to NT/10.

NLJUMP (INPUT, OPTIONAL) integer(i4b)

On entry, the skipping value for the low-pass filter.

By default, NLJUMP is set to NL/10.

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, the maximum number of robustness iterations.

The default is 15. This argument is not used if ROBUST=FALSE.

RW (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On output, final robustness weights. All RW elements are 1 if ROBUST=FALSE .

RW must verify: size(RW) = size(Y).

NO (OUTPUT, OPTIONAL) integer(i4b)

On output, if

ROBUST=TRUE : the number of robustness iterations. The iterations end if a convergence criterion is met or if the number is MAXITER.

ROBUST=FALSE : NO is set to 0.

OK (OUTPUT, OPTIONAL) logical(lgl)

On output, if

ROBUST=TRUE : OK is set to TRUE if the convergence criterion is met and to FALSE otherwise.

ROBUST=FALSE : OK is set to TRUE.

Further Details¶

This subroutine is a FORTRAN90 implementation of subroutine STLEZ developped by Cleveland and coworkers at AT&T Bell Laboratories.

At a minimum, COMP_STLEZ requires specifying the periodicity of the data (e.g. NP, 12 for monthly), the width of the LOESS smoother used to smooth the cyclic seasonal sub-series (e.g. NS) and the degree of the locally-fitted polynomial in seasonal (e.g. ISDEG) and trend (e.g. ITDEG) smoothing.

COMP_STLEZ sets, by default, others parameters of the STL procedure to the values recommended in Cleveland et al. (1990). It also includes tests of convergence if robust iterations are carried out. Otherwise, COMP_STLEZ is similar to COMP_STL.

For further details, see description of COMP_STL and:

Cleveland, R.B., Cleveland, W.S., McRae, J.E., and Terpenning, I.:

STL: A Seasonal-Trend Decomposition Procedure Based on Loess. Statistics Research Report, AT&T Bell Laboratories.

Cleveland, R.B., Cleveland, W.S., McRae, J.E., and Terpenning, I., 1990:

STL: A Seasonal-Trend Decomposition Procedure Based on Loess. J. Official Stat., 6, 3-73.

Crotinger, J., 2017:

Java implementation of Seasonal-Trend-Loess time-series decomposition algorithm. https://github.com/ServiceNow/stl-decomp-4j

`subroutine comp_stlez ( y, np, ns, isdeg, itdeg, robust, season, trend, ni, nt, nl, ildeg, nsjump, ntjump, nljump, maxiter, rw, no, ok )`¶

Purpose¶

COMP_STLEZ decomposes the (time series) columns of a matrix into seasonal and trend components. It returns the components and, optionally, the robustness weights.

COMP_STLEZ offers an easy to use version of COMP_STL subroutine, also included in STATPACK, by defaulting most parameters values associated with the three LOESS smoothers used in COMP_STL.

Arguments¶

Y (INPUT) real(stnd), dimension(:,:)

On entry, the matrix to be decomposed.

NP (INPUT) integer(i4b)

On entry, the period of the seasonal component. For example, if the time series is monthly with a yearly cycle, then NP=12 should be used. NP must be greater than 1.

NS (INPUT) integer(i4b)

On entry, the length of the seasonal smoother. The value of NS should be an odd integer greater than or equal to 3; NS>6 is recommended. As NS increases the values of the seasonal component at a given point in the seasonal cycle (e.g., January values of a monthly series with a yearly cycle) become smoother.

ISDEG (INPUT) integer(i4b)

On entry, the degree of locally-fitted polynomial in seasonal smoothing. The value must be 0 or 1.

ITDEG (INPUT) integer(i4b)

On entry, the degree of locally-fitted polynomial in trend smoothing. The value must be 0, 1 or 2.

ROBUST (INPUT) logical(lgl)

On entry, TRUE if robustness iterations are to be used, FALSE otherwise. Robustness iterations are carried out until convergence of both seasonal and trend components, with MAXITER iterations maximum. Convergence occurs if the maximum changes in individual seasonal and trend fits are less than 1% of the component’s range after the previous iteration.

SEASON (OUTPUT) real(stnd), dimension(:,:)

On output, the seasonal components.

SEASON must verify: size(SEASON,1) = size(Y,1) and size(SEASON,2) = size(Y,2).

TREND (OUTPUT) real(stnd), dimension(:,:)

On output, the trend components.

TREND must verify: size(TREND,1) = size(Y,1) and size(TREND,2) = size(Y,2).

NI (INPUT, OPTIONAL) integer(i4b)

On entry, the number of loops for updating the seasonal and trend components. The value of NI should be a positive integer.

By default, NI=2 if ROBUST=FALSE and NI=1 if ROBUST=TRUE.

NT (INPUT, OPTIONAL) integer(i4b)

On entry, the length of the trend smoother. The value of NT should be an odd integer greater than or equal to 3.

A value of NT between 1.5 * NP and 2 * NP is recommended. As NT increases the values of the trend component become smoother.

By default, NT is set to the smallest odd integer greater than or equal to (1.5 * NP) / (1-(1.5/NS)).

NL (INPUT, OPTIONAL) integer(i4b)

On entry, the length of the low-pass filter. The value of NL should be an odd integer greater than or equal to 3.

The smallest odd integer greater than or equal to NP is recommended.

By default, NL is set to the smallest odd integer greater than or equal to NP.

ILDEG (INPUT, OPTIONAL) integer(i4b)

On entry, the degree of locally-fitted polynomial in low-pass smoothing.

By default, ILDEG is set to ITDEG.

NSJUMP (INPUT, OPTIONAL) integer(i4b)

On entry, the skipping value for seasonal smoothing. The seasonal smoother skips ahead NSJUMP points and then linearly interpolates in between. The value of NSJUMP should be a positive integer; if NSJUMP=1, a seasonal smooth is calculated at all size(Y) points. To make the procedure run faster, a reasonable choice for NSJUMP is 10% or 20% of NS.

By default, NSJUMP is set to NS/10.

NTJUMP (INPUT, OPTIONAL) integer(i4b)

On entry, the skipping value for trend smoothing.

By default, NTJUMP is set to NT/10.

NLJUMP (INPUT, OPTIONAL) integer(i4b)

On entry, the skipping value for the low-pass filter.

By default, NLJUMP is set to NL/10.

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, the maximum number of robustness iterations.

The default is 15. This argument is not used if ROBUST=FALSE.

RW (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On output, final robustness weights. All RW elements are 1 if ROBUST=FALSE .

RW must verify: size(RW,1) = size(Y,1) and size(RW,2) = size(Y,2).

NO (OUTPUT, OPTIONAL) integer(i4b), dimension(:)

On output, if

ROBUST=TRUE : NO(i) is the number of robustness iterations for each time series Y(:,i). The iterations end if a convergence criterion is met or if the number is MAXITER for each time series Y(:,i).

ROBUST=FALSE : NO(:) is set to 0.

NO must verify: size(NO) = size(Y,2).

OK (OUTPUT, OPTIONAL) logical(lgl), dimension(:)

On output, if

ROBUST=TRUE : OK(i) is set to TRUE if the convergence criterion is met for time series Y(:,i) and to FALSE otherwise.

ROBUST=FALSE : OK(:) is set to TRUE.

OK must verify: size(OK) = size(Y,2).

Further Details¶

This subroutine is a FORTRAN90 implementation of subroutine STLEZ developped by Cleveland and coworkers at AT&T Bell Laboratories.

At a minimum, COMP_STLEZ requires specifying the periodicity of the data (e.g. NP, 12 for monthly), the width of the LOESS smoother used to smooth the cyclic seasonal sub-series (e.g. NS) and the degree of the locally-fitted polynomial in seasonal (e.g. ISDEG) and trend (e.g. ITDEG) smoothing.

COMP_STLEZ sets, by default, others parameters of the STL procedure to the values recommended in Cleveland et al. (1990). It also includes tests of convergence if robust iterations are carried out. Otherwise, COMP_STLEZ is similar to COMP_STL.

For further details, see:

Cleveland, R.B., Cleveland, W.S., McRae, J.E., and Terpenning, I.:

STL: A Seasonal-Trend Decomposition Procedure Based on Loess. Statistics Research Report, AT&T Bell Laboratories.

Cleveland, R.B., Cleveland, W.S., McRae, J.E., and Terpenning, I., 1990:

STL: A Seasonal-Trend Decomposition Procedure Based on Loess. J. Official Stat., 6, 3-73.

Crotinger, J., 2017:

Java implementation of Seasonal-Trend-Loess time-series decomposition algorithm. https://github.com/ServiceNow/stl-decomp-4j

`subroutine comp_stl ( y, np, ni, no, isdeg, itdeg, ildeg, nsjump, ntjump, nljump, ns, nt, nl, rw, season, trend )`¶

Purpose¶

COMP_STL decomposes a time series into seasonal and trend components using LOESS smoothers. It returns the components and robustness weights.

Arguments¶

Y (INPUT) real(stnd), dimension(:)

On entry, the time series to be decomposed.

NP (INPUT) integer(i4b)

On entry, the period of the seasonal component. For example, if the time series is monthly with a yearly cycle, then NP=12 should be used. NP must be greater than 1.

NI (INPUT) integer(i4b)

On entry, the number of loops for updating the seasonal and trend components. The value of NI should be a strictly positive integer.

NO (INPUT) integer(i4b)

On entry, the number of robustness iterations. The value of NO should be a positive integer.

ISDEG (INPUT) integer(i4b)

On entry, the degree of locally-fitted polynomial in seasonal smoothing. The value must be 0 or 1.

ITDEG (INPUT) integer(i4b)

On entry, the degree of locally-fitted polynomial in trend smoothing. The value must be 0, 1 or 2.

ILDEG (INPUT) integer(i4b)

On entry, the degree of locally-fitted polynomial in low-pass smoothing. The value must be 0, 1 or 2.

NSJUMP (INPUT/OUTPUT) integer(i4b)

On entry, the skipping value for seasonal smoothing. The seasonal smoother skips ahead NSJUMP points and then linearly interpolates in between. The value of NSJUMP should be a positive integer; if NSJUMP=1, a seasonal smooth is calculated at all size(Y) points. To make the procedure run faster, a reasonable choice for NSJUMP is 10% or 20% of NS.

NTJUMP (INPUT/OUTPUT) integer(i4b)

On entry, the skipping value for trend smoothing.

NLJUMP (INPUT/OUTPUT) integer(i4b)

On entry, the skipping value for the low-pass filter.

NS (INPUT/OUTPUT) integer(i4b)

On entry, the length of the seasonal smoother. The value of NS should be an odd integer greater than or equal to 3; NS>6 is recommended. As NS increases the values of the seasonal component at a given point in the seasonal cycle (e.g., January values of a monthly series with a yearly cycle) become smoother.

NT (INPUT/OUTPUT) integer(i4b)

On entry, the length of the trend smoother. The value of NT should be an odd integer greater than or equal to 3. A value of NT between 1.5 * NP and 2 * NP is recommended. As NT increases the values of the trend component become smoother.

NL (INPUT/OUTPUT) integer(i4b)

On entry, the length of the low-pass filter. The value of NL should be an odd integer greater than or equal to 3. The smallest odd integer greater than or equal to NP is recommended.

RW (OUTPUT) real(stnd), dimension(:)

On output, final robustness weights. All RW elements are 1 if NO=0 .

RW must verify: size(RW) = size(Y).

SEASON (OUTPUT) real(stnd), dimension(:)

On output, the seasonal component.

SEASON must verify: size(SEASON) = size(Y).

TREND (OUTPUT) real(stnd), dimension(:)

On output, the trend component.

TREND must verify: size(TREND) = size(Y).

Further Details¶

This subroutine is a FORTRAN90 implementation of subroutine STL developped by Cleveland and coworkers at AT&T Bell Laboratories.

This subroutine decomposes a time series into seasonal, trend and residual components. The algorithm uses LOESS interpolation and smoothers to smooth the time series and estimate the seasonal (or harmonic) component and the trend. This process is iterative with many steps and may include robustness iterations that take advantage of the weighted-least-squares underpinnings of LOESS to remove the effects of outliers.

There are three LOESS smoothers in COMP_STL and each require three parameters: a width, a degree, and a jump. The width specifies the number of data points that the local interpolation uses to smooth each point, the degree specifies the degree of the local polynomial that is fit to the data, and the jump specifies how many points are skipped between LOESS interpolations, with linear interpolation being done between these points.

The LOESS smoother for estimating the trend is specified with the following parameters: a width (e.g. NT), a degree (e.g. ITDEG) and a jump (e.g. NTJUMP).

The LOESS smoother for estimating the seasonal component is specified with the following parameters: a width (e.g. NS), a degree (e.g. ISDEG) and a jump (e.g. NSJUMP).

The LOESS smoother for low-pass filtering is specified with the following parameters: a width (e.g. NL), a degree (e.g. ILDEG) and a jump (e.g. NLJUMP).

If the NO argument is set to an integer value greater than 0, the process includes also robustness iterations that take advandages of the weighted-least-squares underpinnings of LOESS to remove the effects of outliers.

Note that, finally, that this subroutine expects equally spaced data with no missing values.

For further details, see:

Cleveland, R.B., Cleveland, W.S., McRae, J.E., and Terpenning, I.:

STL: A Seasonal-Trend Decomposition Procedure Based on Loess. Statistics Research Report, AT&T Bell Laboratories.

Cleveland, R.B., Cleveland, W.S., McRae, J.E., and Terpenning, I., 1990:

STL: A Seasonal-Trend Decomposition Procedure Based on Loess. J. Official Stat., 6, 3-73.

Crotinger, J., 2017:

Java implementation of Seasonal-Trend-Loess time-series decomposition algorithm. https://github.com/ServiceNow/stl-decomp-4j

`subroutine comp_stl ( y, np, ni, no, isdeg, itdeg, ildeg, nsjump, ntjump, nljump, ns, nt, nl, rw, season, trend )`¶

Purpose¶

COMP_STL decomposes the (time series) columns of a matrix into seasonal and trend components using LOESS smoothers. It returns the components and robustness weights.

Arguments¶

Y (INPUT) real(stnd), dimension(:,:)

On entry, the time series to be decomposed.

NP (INPUT) integer(i4b)

On entry, the period of the seasonal component. For example, if the time series is monthly with a yearly cycle, then NP=12 should be used. NP must be greater than 1.

NI (INPUT) integer(i4b)

On entry, the number of loops for updating the seasonal and trend components. The value of NI should be a strictly positive integer.

NO (INPUT) integer(i4b)

On entry, the number of robustness iterations. The value of NO should be a positive integer.

ISDEG (INPUT) integer(i4b)

On entry, the degree of locally-fitted polynomial in seasonal smoothing. The value must be 0 or 1.

ITDEG (INPUT) integer(i4b)

On entry, the degree of locally-fitted polynomial in trend smoothing. The value must be 0, 1 or 2.

ILDEG (INPUT) integer(i4b)

On entry, the degree of locally-fitted polynomial in low-pass smoothing. The value must be 0, 1 or 2.

NSJUMP (INPUT/OUTPUT) integer(i4b)

On entry, the skipping value for seasonal smoothing. The seasonal smoother skips ahead NSJUMP points and then linearly interpolates in between. The value of NSJUMP should be a positive integer; if NSJUMP=1, a seasonal smooth is calculated at all size(Y) points. To make the procedure run faster, a reasonable choice for NSJUMP is 10% or 20% of NS.

NTJUMP (INPUT/OUTPUT) integer(i4b)

On entry, the skipping value for trend smoothing.

NLJUMP (INPUT/OUTPUT) integer(i4b)

On entry, the skipping value for the low-pass filter.

NS (INPUT/OUTPUT) integer(i4b)

On entry, the length of the seasonal smoother. The value of NS should be an odd integer greater than or equal to 3; NS>6 is recommended. As NS increases the values of the seasonal component at a given point in the seasonal cycle (e.g., January values of a monthly series with a yearly cycle) become smoother.

NT (INPUT/OUTPUT) integer(i4b)

On entry, the length of the trend smoother. The value of NT should be an odd integer greater than or equal to 3. A value of NT between 1.5 * NP and 2 * NP is recommended. As NT increases the values of the trend component become smoother.

NL (INPUT/OUTPUT) integer(i4b)

On entry, the length of the low-pass filter. The value of NL should be an odd integer greater than or equal to 3. The smallest odd integer greater than or equal to NP is recommended.

RW (OUTPUT) real(stnd), dimension(:,:)

On output, final robustness weights. All RW elemets are 1 if NO=0 .

RW must verify: size(RW,1) = size(Y,1) and size(RW,2) = size(Y,2).

SEASON (OUTPUT) real(stnd), dimension(:,:)

On output, the seasonal components.

SEASON must verify: size(SEASON,1) = size(Y,1) and size(SEASON,2) = size(Y,2).

TREND (OUTPUT) real(stnd), dimension(:,:)

On output, the trend components.

TREND must verify: size(TREND,1) = size(Y,1) and size(TREND,2) = size(Y,2).

Further Details¶

This subroutine is a FORTRAN90 implementation of subroutine STL developped by Cleveland and coworkers at AT&T Bell Laboratories.

This subroutine decomposes a multi-channel time series into seasonal, trend and residual components. The algorithm uses LOESS interpolation and smoothers to smooth the multi-channel time series and estimate the seasonal (or harmonic) components and the trends. This process is iterative with many steps and may include robustness iterations that take advantage of the weighted-least-squares underpinnings of LOESS to remove the effects of outliers.

There are three LOESS smoothers in COMP_STL and each require three parameters: a width, a degree, and a jump. The width specifies the number of data points that the local interpolation uses to smooth each point, the degree specifies the degree of the local polynomial that is fit to the data, and the jump specifies how many points are skipped between LOESS interpolations, with linear interpolation being done between these points.

The LOESS smoother for estimating the trend is specified with the following parameters: a width (e.g. NT), a degree (e.g. ITDEG) and a jump (e.g. NTJUMP).

The LOESS smoother for estimating the seasonal component is specified with the following parameters: a width (e.g. NS), a degree (e.g. ISDEG) and a jump (e.g. NSJUMP).

The LOESS smoother for low-pass filtering is specified with the following parameters: a width (e.g. NL), a degree (e.g. ILDEG) and a jump (e.g. NLJUMP).

If the NO argument is set to an integer value greater than 0, the process includes also robustness iterations that take advandages of the weighted-least-squares underpinnings of LOESS to remove the effects of outliers.

Note that, finally, that this subroutine expects equally spaced data with no missing values.

For further details, see:

Cleveland, R.B., Cleveland, W.S., McRae, J.E., and Terpenning, I.:

STL: A Seasonal-Trend Decomposition Procedure Based on Loess. Statistics Research Report, AT&T Bell Laboratories.

Cleveland, R.B., Cleveland, W.S., McRae, J.E., and Terpenning, I., 1990:

STL: A Seasonal-Trend Decomposition Procedure Based on Loess. J. Official Stat., 6, 3-73.

Crotinger, J., 2017:

Java implementation of Seasonal-Trend-Loess time-series decomposition algorithm. https://github.com/ServiceNow/stl-decomp-4j

`subroutine ma ( x, len, ave)`¶

Purpose¶

Smooth the vector X with a moving average of length LEN and output the result in the vector AVE.

Arguments¶

X (INPUT) real(stnd), dimension(:)

On entry, the vector to smooth.

LEN (INPUT) integer(i4b)

On entry, the length of the moving average. The argument LEN must be >=1 and < size(x).

AVE (OUTPUT) real(stnd), dimension(size(x))

On output, AVE(1:size(X)-LEN+1) contains the smoothed values and AVE(size(X)-LEN+2:size(X)) is unchanged.

Further Details¶

This subroutine is a low-level subroutine used by subroutines COMP_STLEZ and COMP_STL.

`subroutine detrend ( vec, trend, orig, slope )`¶

Purpose¶

Subroutine DETREND detrends a time series (e.g. the argument VEC).

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

The time series vector to be detrended.

TREND (INPUT) integer(i4b)

If:

TREND=1 The mean of the time series is removed

TREND=2 The drift from the time series is removed by using the formula:

drift = (VEC(size(VEC)) - VEC(1))/(size(VEC) - 1)

TREND=3 The least-squares line from the time series is removed.

For other values of TREND nothing is done.

ORIG (OUTPUT, OPTIONAL) real(stnd)

On exit, the constant term if TREND=1 or 3.

SLOPE (OUTPUT, OPTIONAL) real(stnd)

On exit, the linear term if TREND=2 or 3.

Further Details¶

On exit, the original time series may be recovered with the formula

VEC(i) = VEC(i) + ORIG + SLOPE * real(i-1,stnd)

for i=1, size(vec), in all the cases.

`subroutine detrend ( mat, trend, orig, slope )`¶

Purpose¶

Subroutine DETREND detrends a multi-channel time series (e.g. the argument MAT). Each row of matrix MAT is a real time series

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

The multi-channel time series matrix to be detrended.

TREND (INPUT) integer(i4b)

If:

TREND=1 The means of the time series are removed

TREND=2 The drifts from the time series are removed by using the formula:

drift(:) = (MAT(:,size(MAT,2)) - MAT(:,1))/(size(MAT,2) - 1)

TREND=3 The least-squares lines from the time series are removed.

For other values of TREND nothing is done.

ORIG (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the constant terms if TREND=1 or 3.

The size of ORIG must verify: size(ORIG) = size(MAT,1) .

SLOPE (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the linear terms if TREND=2 or 3.

The size of SLOPE must verify: size(SLOPE) = size(MAT,1) .

Further Details¶

On exit, the original time series may be recovered with the formula

MAT(j,i) = MAT(j,i) + ORIG(j) + SLOPE(j) * real(i-1,stnd)

for i=1, size(MAT,2) and j=1, size(MAT,1), in all the cases.

`subroutine hwfilter ( vec, pl, ph, initfft, trend, win )`¶

Purpose¶

Subroutine HWFILTER filters a time series (e.g. the argument VEC) in the frequency band limited by periods PL and PH by windowed filtering (PL and PH are expressed in number of points, i.e. PL=6(18) and PH=32(96) selects periods between 1.5 yrs and 8 yrs for quarterly (monthly) data).

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

The time series vector to be filtered.

Size(VEC) must be greater or equal to 4.

PL (INPUT) integer(i4b)

Minimum period of oscillation of desired component. Use PL=0 for high-pass filtering frequencies corresponding to periods shorter than PH, PL must be equal to 0 or greater or equal to 2. Moreover, PL must be less or equal to size(VEC).

PH (INPUT) integer(i4b)

Maximum period of oscillation of desired component. USE PH=0 for low-pass filtering frequencies corresponding to periods longer than PL. PH must be equal to 0 or greater or equal to 2. Moreover, PH must be less or equal to size(VEC).

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if INITFFT is set to false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine HWFILTER in order to sets up constants and functions for use by subroutine FFT_ROW which is called inside subroutine HWFILTER (the call to INIT_FFT must have the following form:

call init_fft( size(VEC) )

If INITFFT is set to true, the call to INIT_FFT is done inside subroutine HWFILTER and a call to END_FFT is also done before leaving subroutine HWFILTER.

The default is INITFFT=true.

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+/-1 The mean of the time series is removed before time filtering

TREND=+/-2 The drift from the time series is removed before time filtering by using the formula: drift = (VEC(size(VEC)) - VEC(1))/(size(VEC) - 1)

TREND=+/-3 The least-squares line from the time series is removed before time filtering.

IF TREND=-1,-2 or -3, the mean, drift or least-squares line is reintroduced post-filtering, respectively. For other values of TREND nothing is done before or after filtering.

WIN (INPUT, OPTIONAL) real(stnd)

By default, Hamming window filtering is used (i.e. WIN=0.54). SET WIN=0.5 for Hanning window or WIN=1 for rectangular window.

WIN must be greater or equal to O.5 and less or equal to 1.

Further Details¶

Use PL=0 for high-pass filtering frequencies corresponding to periods shorter than PH, or PH=0 for low-pass filtering frequencies corresponding to periods longer than PL.

Setting PH<PL is also allowed and performs band rejection of periods between PH and PL (i.e. in that case the meaning of the PL and PH arguments are reversed).

Examples:

For quarterly data:

call hwfilter( vec, pl=6, ph=32) returns component with periods between 1.5 and 8 yrs.

For monthly data:

call hwfilter( vec, pl=0, ph=24) returns component with all periods less than 2 yrs.

For more details and algorithm, see:

Iacobucci, A., and Noullez, A., 2005:

A Frequency Selective Filter for Short-Length Time Series. Computational Economics, 25,75-102.

`subroutine hwfilter ( mat, pl, ph, initfft, trend, win, max_alloc )`¶

Purpose¶

Subroutine HWFILTER filters a multi-channel time series (e.g. the argument MAT) in the frequency band limited by periods PL and PH by windowed filtering (PL and PH are expressed in number of points, i.e. PL=6(18) and PH=32(96) selects periods between 1.5 yrs and 8 yrs for quarterly (monthly) data).

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

The multi-channel time series matrix to be filtered. Each column of MAT corresponds to one observation.

Size(MAT,2) must be greater or equal to 4.

PL (INPUT) integer(i4b)

Minimum period of oscillation of desired component. Use PL=0 for high-pass filtering frequencies corresponding to periods shorter than PH, PL must be equal to 0 or greater or equal to 2. Moreover, PL must be less or equal to size(MAT,2).

PH (INPUT) integer(i4b)

Maximum period of oscillation of desired component. USE PH=0 for low-pass filtering frequencies corresponding to periods longer than PL. PH must be equal to 0 or greater or equal to 2. Moreover, PH must be less or equal to size(MAT,2).

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if INITFFT is set to false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine HWFILTER in order to sets up constants and functions for use by subroutine FFT_ROW which is called inside subroutine HWFILTER (the call to INIT_FFT must have the following form:

call init_fft( shape(MAT), dim=2_i4b )

If INITFFT is set to true, the call to INIT_FFT is done inside subroutine HWFILTER and a call to END_FFT is also done before leaving subroutine HWFILTER.

The default is INITFFT=true.

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+/-1 The means of the time series are removed before time filtering

TREND=+/-2 The drifts from the time series are removed before time filtering by using the formula: drift(:) = (MAT(:,size(MAT,2)) - MAT(:,1))/(size(MAT,2) - 1)

TREND=+/-3 The least-squares lines from the time series are removed before time filtering.

IF TREND=-1,-2 or -3, the means, drifts or least-squares lines are reintroduced post-filtering, respectively. For other values of TREND nothing is done before or after filtering.

WIN (INPUT, OPTIONAL) real(stnd)

By default, Hamming window filtering is used (i.e. WIN=0.54). SET WIN=0.5 for Hanning window or WIN=1 for rectangular window.

WIN must be greater or equal to O.5 and less or equal to 1.

MAX_ALLOC (INPUT, OPTIONAL) integer(i4b)

MAX_ALLOC is a factor which allows to reduce the workspace used to compute the Fourier transform of the data if necessary at the expense of increasing the computing time. MAX_ALLOC must be greater or equal to 1 and less or equal to size(MAT,1).

The default is MAX_ALLOC= size(MAT,1).

Further Details¶

Use PL=0 for high-pass filtering frequencies corresponding to periods shorter than PH, or PH=0 for low-pass filtering frequencies corresponding to periods longer than PL.

Setting PH<PL is also allowed and performs band rejection of periods between PH and PL (i.e. in that case the meaning of the PL and PH arguments are reversed).

Examples:

For quarterly data:

call hwfilter( mat, pl=6, ph=32) returns components with periods between 1.5 and 8 yrs.

For monthly data:

call hwfilter( mat, pl=0, ph=24) returns components with all periods less than 2 yrs.

For more details and algorithm, see :

Iacobucci, A., and Noullez, A., 2005:

A Frequency Selective Filter for Short-Length Time Series. Computational Economics, 25,75-102.

`subroutine hwfilter2 ( vec, pl, ph, trend, win )`¶

Purpose¶

Subroutine HWFILTER2 filters a time series (e.g. the argument VEC) in the frequency band limited by periods PL and PH by windowed filtering (PL and PH are expressed in number of points, i.e. PL=6(18) and PH=32(96) selects periods between 1.5 yrs and 8 yrs for quarterly (monthly) data).

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

The time series vector to be filtered. Size(VEC) must be greater or equal to 4.

PL (INPUT) integer(i4b)

Minimum period of oscillation of desired component. Use PL=0 for high-pass filtering frequencies corresponding to periods shorter than PH, PL must be equal to 0 or greater or equal to 2. Moreover, PL must be less or equal to size(VEC).

PH (INPUT) integer(i4b)

Maximum period of oscillation of desired component. USE PH=0 for low-pass filtering frequencies corresponding to periods longer than PL. PH must be equal to 0 or greater or equal to 2. Moreover, PH must be less or equal to size(VEC).

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+/-1 The mean of the time series is removed before time filtering

TREND=+/-2 The drift from the time series is removed before time filtering by using the formula: drift = (VEC(size(VEC)) - VEC(1))/(size(VEC) - 1)

TREND=+/-3 The least-squares line from the time series is removed before time filtering.

IF TREND=-1,-2 or -3, the mean, drift or least-squares line is reintroduced post-filtering, respectively. For other values of TREND nothing is done before or after filtering.

WIN (INPUT, OPTIONAL) real(stnd)

By default, Hamming window filtering is used (i.e. WIN=0.54). SET WIN=0.5 for Hanning window or WIN=1 for rectangular window.

WIN must be greater or equal to O.5 and less or equal to 1.

Further Details¶

Use PL=0 for high-pass filtering frequencies corresponding to periods shorter than PH, or PH=0 for low-pass filtering frequencies corresponding to periods longer than PL.

Setting PH<PL is also allowed and performs band rejection of periods between PH and PL (i.e. in that case the meaning of the PL and PH arguments are reversed).

The unique difference between HWFILTER2 and HWFILTER is the use of the Goertzel method for computing the Fourier transform of the data instead of a Fast Fourier Transform algorithm.

Examples:

For quarterly data:

call hwfilter2( vec, pl=6, ph=32) returns component with periods between 1.5 and 8 yrs.

For monthly data:

call hwfilter2( vec, pl=0, ph=24) returns component with all periods less than 2 yrs.

For more details and algorithm, see :

Iacobucci, A., and Noullez, A., 2005:

A Frequency Selective Filter for Short-Length Time Series. Computational Economics, 25,75-102.

Goertzel, G., 1958:

An Algorithm for the Evaluation of Finite Trigonometric Series. The American Mathematical Monthly, Vol. 65, No. 1, pp. 34-35

`subroutine hwfilter2 ( mat, pl, ph, trend, win )`¶

Purpose¶

Subroutine HWFILTER2 filters a multi-channel time series (e.g. the argument MAT) in the frequency band limited by periods PL and PH by windowed filtering (PL and PH are expressed in number of points, i.e. PL=6(18) and PH=32(96) selects periods between 1.5 yrs and 8 yrs for quarterly (monthly) data).

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

The multi-channel time series matrix to be filtered. Each column of MAT corresponds to one observation.

Size(MAT,2) must be greater or equal to 4.

PL (INPUT) integer(i4b)

Minimum period of oscillation of desired component. Use PL=0 for high-pass filtering frequencies corresponding to periods shorter than PH, PL must be equal to 0 or greater or equal to 2. Moreover, PL must be less or equal to size(MAT,2).

PH (INPUT) integer(i4b)

Maximum period of oscillation of desired component. USE PH=0 for low-pass filtering frequencies corresponding to periods longer than PL. PH must be equal to 0 or greater or equal to 2. Moreover, PH must be less or equal to size(MAT,2).

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+/-1 The means of the time series are removed before time filtering

TREND=+/-2 The drifts from the time series are removed before time filtering by using the formula: drift(:) = (MAT(:,size(MAT,2)) - MAT(:,1))/(size(MAT,2) - 1)

TREND=+/-3 The least-squares lines from the time series are removed before time filtering.

IF TREND=-1,-2 or -3, the means, drifts or least-squares lines are reintroduced post-filtering, respectively. For other values of TREND nothing is done before or after filtering.

WIN (INPUT, OPTIONAL) real(stnd)

By default, Hamming window filtering is used (i.e. WIN=0.54). SET WIN=0.5 for Hanning window or WIN=1 for rectangular window.

WIN must be greater or equal to O.5 and less or equal to 1.

Further Details¶

Use PL=0 for high-pass filtering frequencies corresponding to periods shorter than PH, or PH=0 for low-pass filtering frequencies corresponding to periods longer than PL.

Setting PH<PL is also allowed and performs band rejection of periods between PH and PL (i.e. in that case the meaning of the PL and PH arguments are reversed).

The unique difference between HWFILTER2 and HWFILTER is the use of the Goertzel method for computing the Fourier transform of the data instead of a Fast Fourier Transform algorithm.

Examples:

For quarterly data:

call hwfilter2( mat, pl=6, ph=32) returns components with periods between 1.5 and 8 yrs.

For monthly data:

call hwfilter2( mat, pl=0, ph=24) returns components with all periods less than 2 yrs.

For more details and algorithm, see :

Iacobucci, A., and Noullez, A., 2005:

A Frequency Selective Filter for Short-Length Time Series. Computational Economics, 25,75-102.

Goertzel, G., 1958:

An Algorithm for the Evaluation of Finite Trigonometric Series. The American Mathematical Monthly, Vol. 65, No. 1, pp. 34-35

`function lp_coef ( pl, k, fc, notest_fc )`¶

Purpose¶

Function LP_COEF computes the K-term least squares approximation to an -ideal- low pass filter with cutoff period PL (e.g. cutoff frequency FC = 1/PL).

This filter has a transfer function with a transition band of width delta surrounding FC, where delta = 4 * pi/K when FC is expressed in radians.

Arguments¶

PL (INPUT) integer(i4b)

Minimum period of oscillation of desired component. The corresponding cutoff frequency is FC=1/PL (i.e. filter has zero response in the interval [FC+1/K,Nyquist] and one response in the interval [0,FC-1/K].

PL must be greater than 2 and FC must verify the following inequalities:

FC - 1/K >= 0

FC + 1/K < 0.5

K (INPUT) integer(i4b)

The number of filter terms to be computed. K must be greater or equal to 3 and odd.

FC (INPUT, OPTIONAL) real(stnd)

The user chosen cutoff frequency in cycles per sample interval. If the optional argument FC is used, the PL argument is not used to determine the cutoff frequency.

FC must verify the following inequalities:

FC - 1/K >= 0

FC + 1/K < 0.5

NOTEST_FC (INPUT, OPTIONAL) logical(lgl)

On input, if this optional logical argument is set to true, the two tests on the cutoff frequency (e.g. FC - 1/K >= 0 and FC + 1/K < 0.5) are bypassed. However, in that case, the cutoff frequency FC must still verify the inequalities 0 < FC < 0.5.

Further Details¶

Function LP_COEF computes symmetric linear low-pass filter coefficients using a least squares approximation to an ideal low-pass filter with convergence factors (i.e. Lanczos window) which reduce overshoot and ripple (Bloomfield, 1976).

This low-pass filter has a transfer function which changes from approximately one to zero in a transition band about the ideal cutoff frequency FC (FC=1/PL), that is from (FC - 1/K) to (FC + 1/K), as discussed in section 6.4 of Bloomfield (1976). The user must specify the cutoff period (or the cutoff frequency) and the number of filter coefficients, which must be odd.

The user must also choose the number of filter terms, K, so that (FC - 1/K) >= 0 and (FC + 1/K) < 0.5 if the optional logical argument NOTEST_FC is not used or is not set to true.

In addition, K must be chosen as a compromise between:

A sharp cutoff, that is, 1/K small; and

Minimizing the number of data points lost by the filtering operations (e.g. (K-1)/2 data points will be lost from each end of the series).

The subroutine returns the normalized low-pass filter coefficients.

This function is adapted from the STARPAC software developed by the National Institute of Standards and Technology (NIST). For more details and algorithm, see

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York, Chapter 6.

`function lp_coef2 ( pl, k, fc, win, notest_fc )`¶

Purpose¶

Function LP_COEF2 computes the K-term least squares approximation to an -ideal- low pass filter with cutoff period PL (e.g. cutoff frequency FC = 1/PL) by windowed filtering (e.g. Hamming window is used).

This filter has a transfer function with a transition band of width delta surrounding FC, where delta = 4 * pi/K when FC is expressed in radians.

Arguments¶

PL (INPUT) integer(i4b)

Minimum period of oscillation of desired component. The corresponding cutoff frequency is FC=1/PL (i.e. filter has zero response in the interval [FC+1/K,Nyquist] and one response in the interval [0,FC-1/K].

PL must be greater than two and FC must verify the following inequalities:

FC - 1/K >= 0

FC + 1/K < 0.5

K (INPUT) integer(i4b)

The number of filter terms to be computed. K must be greater or equal to 3 and odd.

FC (INPUT, OPTIONAL) real(stnd)

The user chosen cutoff frequency in cycles per sample interval. If the optional argument FC is used, the PL argument is not used to determine the cutoff frequency.

FC must verify the following inequalities:

FC - 1/K >= 0

FC + 1/K < 0.5

WIN (INPUT, OPTIONAL) real(stnd)

By default, Hamming window filtering is used (i.e. WIN=0.54). Set WIN=0.5 for Hanning window or WIN=1 for rectangular window.

WIN must be greater or equal to O.5 and less or equal to 1, otherwise WIN is reset to 0.54.

NOTEST_FC (INPUT, OPTIONAL) logical(lgl)

On input, if this optional logical argument is set to true, the two tests on the cutoff frequency (e.g. FC - 1/K >= 0 and FC + 1/K < 0.5) are bypassed. However, in that case, the cutoff frequency FC must still verify the inequalities 0 < FC < 0.5.

Further Details¶

Function LP_COEF2 computes symmetric linear low-pass filter coefficients using a least squares approximation to an ideal low-pass filter. The Hamming window is used to reduce overshoot and ripple in the transfert function of the ideal low-pass filter.

This low-pass filter has a transfer function which changes from approximately one to zero in a transition band about the ideal cutoff frequency FC (FC=1/PL), that is from (FC - 1/K) to (FC + 1/K), as discussed in section 6.4 of Bloomfield (1976). The user must specify the cutoff period (or the cutoff frequency) and the number of filter coefficients, which must be odd.

The user must also choose the number of filter terms, K, so that (FC - 1/K) >= 0 and (FC + 1/K) < 0.5 if the optional logical argument NOTEST_FC is not used or is not set to true.

The overshoot and the associated ripples in the ideal transfert function are reduced by the use of the Hamming window.

In addition, K must be chosen as a compromise between:

A sharp cutoff, that is, 1/K small; and

Minimizing the number of data points lost by the filtering operations (e.g. (K-1)/2 data points will be lost from each end of the series).

The subroutine returns the normalized low-pass filter coefficients.

For more details and algorithm, see

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York, Chapter 6.

`function hp_coef ( ph, k, fc, notest_fc )`¶

Purpose¶

Function HP_COEF computes the K-term least squares approximation to an -ideal- high pass filter with cutoff period PH (e.g. cutoff frequency FC = 1/PH).

This filter has a transfer function with a transition band of width delta surrounding FC, where delta = 4 * pi/K when FC is expressed in radians.

Arguments¶

PH (INPUT) integer(i4b)

Maximum period of oscillation of desired component. The corresponding cutoff frequency is FC=1/PH (i.e. filter has one response in the interval [FC+1/K,Nyquist] and zero response in the interval [0,FC-1/K].

PH must be greater than two and FC must verify the following inequalities:

FC - 1/K >= 0

FC + 1/K < 0.5

K (INPUT) integer(i4b)

The number of filter terms to be computed. K must be greater or equal to 3 and odd.

FC (INPUT, OPTIONAL) real(stnd)

The user chosen cutoff frequency in cycles per sample interval. If the optional argument FC is used, the PH argument is not used to determine the cutoff frequency.

FC must verify the following inequalities:

FC - 1/K >= 0

FC + 1/K < 0.5

NOTEST_FC (INPUT, OPTIONAL) logical(lgl)

On input, if this optional logical argument is set to true, the two tests on the cutoff frequency (e.g. FC - 1/K >= 0 and FC + 1/K < 0.5) are bypassed. However, in that case, the cutoff frequency FC must still verify the inequalities 0 < FC < 0.5.

Further Details¶

Function HP_COEF computes symmetric linear high-pass filter coefficients from the corresponding low-pass filter as given by function LP_COEF. This is equivalent to subtracting the low-pass filtered series from the original time series.

This high-pass filter has a transfer function which changes from approximately zero to one in a transition band about the ideal cutoff frequency FC (FC=1/PH), that is from (FC - 1/K) to (FC + 1/K), as discussed in section 6.4 of Bloomfield (1976). The user must specify the cutoff period (or the cutoff frequency) and the number of filter coefficients, which must be odd.

The user must also choose the number of filter terms, K, so that (FC - 1/K) >= 0 and (FC + 1/K) < 0.5 if the optional logical argument NOTEST_FC is not used or is not set to true.

In addition, K must be chosen as a compromise between:

A sharp cutoff, that is, 1/K small; and

Minimizing the number of data points lost by the filtering operations (e.g. (K-1)/2 data points will be lost from each end of the series).

The subroutine returns the high-pass filter coefficients.

This function is adapted from the STARPAC software developed by the National Institute of Standards and Technology (NIST).

For more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York, Chapter 6.

`function hp_coef2 ( ph, k, fc, win, notest_fc )`¶

Purpose¶

Function HP_COEF2 computes the K-term least squares approximation to an -ideal- high pass filter with cutoff period PH (e.g. cutoff frequency FC = 1/PH) by windowed filtering (e.g. Hamming window is used).

This filter has a transfer function with a transition band of width delta surrounding FC, where delta = 4 * pi/K when FC is expressed in radians.

Arguments¶

PH (INPUT) integer(i4b)

Maximum period of oscillation of desired component. The corresponding cutoff frequency is FC=1/PH (i.e. filter has one response in the interval [FC+1/K,Nyquist] and zero response in the interval [0,FC-1/K].

PH must be greater than two and FC must verify the following inequalities:

FC - 1/K >= 0

FC + 1/K < 0.5

K (INPUT) integer(i4b)

The number of filter terms to be computed. K must be greater or equal to 3 and odd.

FC (INPUT, OPTIONAL) real(stnd)

The user chosen cutoff frequency in cycles per sample interval. If the optional argument FC is used, the PH argument is not used to determine the cutoff frequency.

FC must verify the following inequalities:

FC - 1/K >= 0

FC + 1/K < 0.5

WIN (INPUT, OPTIONAL) real(stnd)

By default, Hamming window filtering is used (i.e. WIN=0.54). Set WIN=0.5 for Hanning window or WIN=1 for rectangular window.

WIN must be greater or equal to O.5 and less or equal to 1, otherwise WIN is reset to 0.54.

NOTEST_FC (INPUT, OPTIONAL) logical(lgl)

On input, if this optional logical argument is set to true, the two tests on the cutoff frequency (e.g. FC - 1/K >= 0 and FC + 1/K < 0.5) are bypassed. However, in that case, the cutoff frequency FC must still verify the inequalities 0 < FC < 0.5.

Further Details¶

Function HP_COEF2 computes symmetric linear high-pass filter coefficients from the corresponding low-pass filter as given by function LP_COEF2. This is equivalent to subtracting the low-pass filtered series from the original time series.

This high-pass filter has a transfer function which changes from approximately zero to one in a transition band about the ideal cutoff frequency FC (FC=1/PH), that is from (FC - 1/K) to (FC + 1/K), as discussed in section 6.4 of Bloomfield (1976). The user must specify the cutoff period (or the cutoff frequency) and the number of filter coefficients, which must be odd.

The user must also choose the number of filter terms, K, so that (FC - 1/K) >= 0 and (FC + 1/K) < 0.5 if the optional logical argument NOTEST_FC is not used or is not set to true.

The overshoot and the associated ripples in the ideal transfert function are reduced by the use of the Hamming window.

In addition, K must be chosen as a compromise between:

A sharp cutoff, that is, 1/K small; and

Minimizing the number of data points lost by the filtering operations (e.g. (K-1)/2 data points will be lost from each end of the series).

The subroutine returns the high-pass filter coefficients.

For more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York, Chapter 6.

`function bd_coef ( pl, ph, k, fch, fcl, notest_fc )`¶

Purpose¶

Function BD_COEF computes the K-term least squares approximation to an -ideal- band pass filter with cutoff periods PL and PH (e.g. cutoff frequencies 1/PL and 1/PH).

PL and PH are expressed in number of points, i.e. PL=6(18) and PH=32(96) selects periods between 1.5 yrs and 8 yrs for quarterly(monthly) data).

Alternatively, the user can directly specify the two cutoff frequencies, FCL and FCH, corresponding to PL and PH.

Arguments¶

PL (INPUT) integer(i4b)

Minimum period of oscillation of desired component. The corresponding cutoff frequency is 1/PL. PL must be greater than two and must verify the following inequalities:

1/PH + 1.3/(K+1) <= 1/PL - 1.3/(K+1)

1/PL + 1/K < 0.5

PH (INPUT) integer(i4b)

Maximum period of oscillation of desired component. The corresponding cutoff frequency is 1/PH. PH must be greater than two and 1/PH must verify the following inequalities:

0 <= 1/PH - 1/K

1/PH + 1.3/(K+1) <= 1/PL - 1.3/(K+1)

K (INPUT) integer(i4b)

The number of filter terms to be computed. K must be greater or equal to 3 and odd.

FCH (INPUT, OPTIONAL) real(stnd)

The user chosen (low) cutoff frequency in cycles per sample interval. If the optional argument FCH is used, the PH argument is not used to determine the (low) cutoff frequency.

FCH must verify the following inequalities:

0 <= FCH - 1/K

FCH + 1.3/(K+1) <= FCL - 1.3/(K+1)

FCL (INPUT, OPTIONAL) real(stnd)

The user chosen (high) cutoff frequency in cycles per sample interval. If the optional argument FCL is used, the PL argument is not used to determine the cutoff (high) frequency.

FCL must verify the following inequalities:

FCH + 1.3/(K+1) <= FCL - 1.3/(K+1)

FCL + 1/K < 0.5

NOTEST_FC (INPUT, OPTIONAL) logical(lgl)

On input, if this optional logical argument is set to true, the tests on the cutoff frequencies FCH and FCL (e.g. FCH - 1/K >= 0 and FCL + 1/K < 0.5) are bypassed. However, in that case FCH and FCL must still verify the inequalities FCH > 0, FCL < 0.5 and FCH + 1.3/(K+1) <= FCL - 1.3/(K+1) .

Further Details¶

Function BD_COEF computes symmetric linear band-pass filter coefficients using a least squares approximation to an ideal band-pass filter that has convergence factors which reduce overshoot and ripple (Bloomfield, 1976).

This band-pass filter is computed as the difference between two low-pass filters with cutoff frequencies 1/PH and 1/PL, respectively (or FCH and FCL).

This band-pass filter has a transfer function which changes from approximately zero to one and one to zero in the transition bands about the ideal cutoff frequencies 1/PH and 1/PL), that is from (1/PH - 1/K) to (1/PH + 1/K) and (1/PL - 1/K) to (1/PL + 1/K), respectively. The user must specify the two cutoff periods and the number of filter coefficients, which must be odd.

The user must also choose the number of filter terms, K, so that:

0<=(1/PH - 1/K)

(1/PH + 1.3/(K+1)) <= (1/PL - 1.3/(K+1))

(1/PL + 1/K) < 0.5

However, if the optional logical argument NOTEST_FC is used and is set to true, the two tests

0<=(1/PH - 1/K)

(1/PL + 1/K) < 0.5

are bypassed.

In addition, K must be chosen as a compromise between:

A sharp cutoff, that is, 1/K small; and

Minimizing the number of data points lost by the filtering operations (e.g. (K-1)/2 data points will be lost from each end of the series).

The subroutine returns the difference between the two corresponding normalized low-pass filter coefficients as computed by function LP_COEF.

For more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York, Chapter 6.

Duchon, C., 1979:

Lanczos filtering in one and two dimensions. Journal of applied meteorology, vol. 18, 1016-1022.

`function bd_coef2 ( pl, ph, k, fch, fcl, win, notest_fc )`¶

Purpose¶

Function BD_COEF2 computes the K-term least squares approximation to an -ideal- band pass filter with cutoff periods PL and PH (e.g. cutoff frequencies 1/PL and 1/PH) by windowed filtering (e.g. Hamming window is used).

PL and PH are expressed in number of points, i.e. PL=6(18) and PH=32(96) selects periods between 1.5 yrs and 8 yrs for quarterly(monthly) data).

Alternatively, the user can directly specify the two cutoff frequencies, FCL and FCH, corresponding to PL and PH.

Arguments¶

PL (INPUT) integer(i4b)

Minimum period of oscillation of desired component. The corresponding cutoff frequency is 1/PL. PL must be greater than two and must verify the following inequalities:

1/PH < 1/PL

1/PL + 1/K < 0.5

PH (INPUT) integer(i4b)

Maximum period of oscillation of desired component. The corresponding cutoff frequency is 1/PH. PH must be greater than two and 1/PH must verify the following inequalities:

0 <= 1/PH - 1/K

1/PH < 1/PL

K (INPUT) integer(i4b)

The number of filter terms to be computed. K must be greater or equal to 3 and odd.

FCH (INPUT, OPTIONAL) real(stnd)

The user chosen (low) cutoff frequency in cycles per sample interval. If the optional argument FCH is used, the PH argument is not used to determine the (low) cutoff frequency.

FCH must verify the following inequalities:

0 <= FCH - 1/K

FCH < FCL

FCL (INPUT, OPTIONAL) real(stnd)

The user chosen (high) cutoff frequency in cycles per sample interval. If the optional argument FCL is used, the PL argument is not used to determine the cutoff (high) frequency.

FCL must verify the following inequalities:

FCH < FCL

FCL + 1/K < 0.5

WIN (INPUT, OPTIONAL) real(stnd)

By default, Hamming window filtering is used (i.e. WIN=0.54). Set WIN=0.5 for Hanning window or WIN=1 for rectangular window.

WIN must be greater or equal to O.5 and less or equal to 1, otherwise WIN is reset to 0.54.

NOTEST_FC (INPUT, OPTIONAL) logical(lgl)

On input, if this optional logical argument is set to true, the tests on the cutoff frequencies FCH and FCL (e.g. FCH - 1/K >= 0 and FCL + 1/K < 0.5) are bypassed. However, in that case FCH and FCL must still verify the inequalities FCH > 0, FCL < 0.5 and FCH < FCL .

Further Details¶

Function BD_COEF2 computes symmetric linear band-pass filter coefficients using a least squares approximation to an ideal band-pass filter. The Hamming window is used to reduce overshoot and ripple in the transfert function of the ideal low-pass filter.

This band-pass filter is computed as the difference between two low-pass filters with cutoff frequencies 1/PH and 1/PL, respectively (or FCH and FCL).

This band-pass filter has a transfer function which changes from approximately zero to one and one to zero in the transition bands about the ideal cutoff frequencies 1/PH and 1/PL), that is from (1/PH - 1/K) to (1/PH + 1/K) and (1/PL - 1/K) to (1/PL + 1/K), respectively. The user must specify the two cutoff periods and the number of filter coefficients, which must be odd. The user must also choose the number of filter terms, K, so that:

0 <= (1/PH - 1/K)

1/PH < 1/PL

(1/PL + 1/K) < 0.5

However, if the optional logical argument NOTEST_FC is used and is set to true, the two tests

0 <= (1/PH - 1/K)

(1/PL + 1/K) < 0.5

are bypassed.

The overshoot and the associated ripples in the ideal transfert function are reduced by the use of the Hamming window.

In addition, K must be chosen as a compromise between:

A sharp cutoff, that is, 1/K small; and

Minimizing the number of data points lost by the filtering operations (e.g. (K-1)/2 data points will be lost from each end of the series).

The subroutine returns the difference between the two corresponding normalized low-pass filter coefficients as computed by function LP_COEF2.

For more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York, Chapter 6.

`function pk_coef ( freq, k, notest_freq )`¶

Purpose¶

Function PK_COEF computes the K-term least squares approximation to an -ideal- band pass filter with peak response near one at the single frequency FREQ (e.g. the peak response is at period=1/FREQ).

Arguments¶

FREQ (INPUT) real(stnd)

The band pass filter will have unit response at the single frequency FREQ. FREQ is expressed in cycles per sample interval.

The frequency FREQ must also be greater or equal to ( 1.3/(K+1) + 1/K ) and less than 0.5 - ( 1.3/(K+1) + 1/K ).

K (INPUT) integer(i4b)

The number of filter terms to be computed. K must be greater or equal to 3 and odd.

NOTEST_FREQ (INPUT, OPTIONAL) logical(lgl)

On input, if this optional logical argument is set to true, the frequency FREQ must only be greater or equal to 1.3/(K+1) and less than 0.5 - 1.3/(K+1).

Further Details¶

Function PK_COEF computes symmetric linear band-pass filter coefficients using a least squares approximation to an ideal band-pass filter that has convergence factors which reduce overshoot and ripple (Bloomfield, 1976).

This band-pass filter is computed as the difference between two low-pass filters with cutoff frequencies FCL and FCH, respectively (Duchon, 1979).

This band-pass filter has a transfer function which changes from approximately zero to one and one to zero in the transition bands about the cutoff frequencies FCH and FCL, that is from (FCH - 1/K) to FREQ and FREQ to (FCL + 1/K), respectively. The user must specify the frequency FREQ with unit response and the number of filter coefficients, which must be odd. The user must also choose the number of filter terms, K, as a compromise between:

A sharp cutoff, that is, 1/K small; and

Minimizing the number of data points lost by the filtering operations (e.g. (K-1)/2 data points will be lost from each end of the series).

The subroutine computes the two cutoff frequencies FCL and FCH as described by Duchon (1979) and returns the difference between the two corresponding normalized low-pass filter coefficients as computed by function LP_COEF.

For more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York, Chapter 6.

Duchon, C., 1979:

Lanczos filtering in one and two dimensions. Journal of applied meteorology, vol. 18, 1016-1022.

`function moddan_coef ( k, smooth_param )`¶

Purpose¶

This function computes the impulse response function (e.g. weights) corresponding to a number of applications of modified Daniell filters as done in subroutine MODDAN_FILTER.

Arguments¶

K (INPUT) integer(i4b)

The number of filter weights to be computed. K must be equal to 2 * (2+sum(SMOOTH_PARAM(:)))- 1

SMOOTH_PARAM (INPUT) integer(i4b), dimension(:)

The array of the half-lengths of the modified Daniell filters to be applied. All the values in SMOOTH_PARAM(:) must be greater than 0.

Size(SMOOTH_PARAM) must be greater or equal to 1.

Further Details¶

For definition, more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York, Chapter 6.

`subroutine freq_func ( nfreq, coef, freqr, four_freq, freq )`¶

Purpose¶

Subroutine FREQ_FUNC computes the frequency response function (e.g. the transfer function) of the symmetric linear filter given by the argument COEF(:).

The frequency response function is computed at NFREQ frequencies regularly sampled between 0 and the Nyquist frequency if the optional logical argument FOUR_FREQ is not used or at the NFREQ Fourier frequencies 2 * pi * j/nfreq for j=0 to NFREQ-1 if this argument is used and set to true.

Arguments¶

NFREQ (INPUT) integer(i4b)

The number of frequencies at which the frequency response function must be evaluated.

COEF (INPUT) real(stnd), dimension(:)

The array of symmetric linear filter coefficients.

Size(COEF) must be greater or equal to 3 and odd.

FREQR (OUTPUT) real(stnd), dimension(NFREQ)

On output, the frequency response function.

FOUR_FREQ (INPUT, OPTIONAL) logical(lgl)

On input, if this argument is set to true the frequency response function is evaluated at the Fourier frequencies 2 * pi * j/nfreq for j=0 to NFREQ-1.

FREQ (OUTPUT, OPTIONAL) real(stnd), dimension(NFREQ)

The NFREQ frequencies, in cycles per sample interval, at which the frequency response function are evaluated.

Further Details¶

For more details, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York, Chapter 6.

Oppenheim, A.V., and Schafer, R.W., 1999:

Discrete-Time Signal Processing. Second Edition. Prentice-Hall, New Jersey.

`subroutine symlin_filter ( vec, coef, trend, nfilt )`¶

Purpose¶

Subroutine SYMLIN_FILTER performs a symmetric filtering operation on an input time series (e.g. the argument VEC).

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

On input, the vector containing the time series to be filtered. On output, the filtered time series is returned in VEC(:NFILT). Note that (size(COEF)-1)/2 data points will be lost from each end of the series, so that NFILT (NFILT= size(VEC) - size(COEF) + 1) time observations are returned and the remainig and ending part of VEC(:) is set to zero.

Size(VEC) must be greater or equal to 4.

COEF (INPUT) real(stnd), dimension(:)

The array of symmetric linear filter coefficients.

Size(COEF) must be odd, greater or equal to 3 and less or equal to size(VEC).

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+/-1 The mean of the time series is removed before time filtering

TREND=+/-2 The drift from the time series is removed before time filtering by using the formula: drift = (VEC(size(VEC)) - VEC(1))/(size(VEC) - 1)

TREND=+/-3 The least-squares line from the time series is removed before time filtering.

IF TREND=-1,-2 or -3, the mean, drift or least-squares line is reintroduced post-filtering, respectively. For other values of TREND nothing is done before or after filtering.

NFILT (OUTPUT, OPTIONAL) integer(i4b)

The number of time observations in the filtered time series. On output, NFILT= size(VEC) - size(COEF) + 1.

Further Details¶

The filtering is done in place and (size(COEF)-1)/2 observations will be lost from each end of the time series.

Note, also, that the filtered time series is shifted in time and is stored in VEC(1:NFILT) on output, with NFILT= size(VEC) - size(COEF) + 1.

The symmetric linear filter coefficients (e.g. the array COEF) can be computed with the help of functions LP_COEF, LP_COEF2, HP_COEF, HP_COEF2, BD_COEF and BD_COEF2.

`subroutine symlin_filter ( mat, coef, trend, nfilt )`¶

Purpose¶

Subroutine SYMLIN_FILTER performs a symmetric filtering operation on an input multi-channel time series (e.g. the argument MAT).

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

The multi-channel time series matrix to be filtered. Each column of MAT corresponds to one observation. On output, the multi-channel filtered time series are returned in MAT(:,:NFILT). Note that (size(COEF)-1)/2 observations will be lost from each end of the multi-channel series, so that NFILT (NFILT= size(MAT,2) - size(COEF) + 1) time observations are returned and the remainig part of MAT(:,:) is set to zero.

Size(MAT,2) must be greater or equal to 4.

COEF (INPUT) real(stnd), dimension(:)

The array of symmetric linear filter coefficients.

Size(COEF) must be odd, greater or equal to 3 and less or equal to size(MAT,2).

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+/-1 The means of the time series are removed before time filtering

TREND=+/-2 The drifts from the time series are removed before time filtering by using the formula: drift(:) = (MAT(:,size(MAT,2)) - MAT(:,1))/(size(MAT,2) - 1)

TREND=+/-3 The least-squares lines from the time series are removed before time filtering.

IF TREND=-1,-2 or -3, the means, drifts or least-squares lines are reintroduced post-filtering, respectively. For other values of TREND nothing is done before or after filtering.

NFILT (OUTPUT, OPTIONAL) integer(i4b)

The number of time observations in the filtered multi-channel time series. On output, NFILT= size(MAT,2) - size(COEF) + 1.

Further Details¶

The filtering is done in place and (size(COEF)-1)/2 observations will be lost from each end of the multi-channel series.

Note, also, that the filtered multi-channel time series is shifted in time and is stored in MAT(:,1:NFILT) on output, with NFILT= size(MAT,2) - size(COEF) + 1.

The symmetric linear filter coefficients (e.g. the array COEF) can be computed with the help of functions LP_COEF, LP_COEF2, HP_COEF, HP_COEF2, BD_COEF and BD_COEF2.

`subroutine symlin_filter2 ( vec, coef, trend, usefft, initfft )`¶

Purpose¶

Subroutine SYMLIN_FILTER2 performs a symmetric filtering operation on an input time series (e.g. the argument VEC).

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

On input, the vector containing the time series to be filtered. On output, the filtered time series is returned.

Size(VEC) must be greater or equal to 4.

COEF (INPUT) real(stnd), dimension(:)

The array of symmetric linear filter coefficients.

Size(COEF) must be odd, greater or equal to 3 and less or equal to size(VEC).

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+/-1 The mean of the time series is removed before time filtering

TREND=+/-2 The drift from the time series is removed before time filtering by using the formula: drift = (VEC(size(VEC)) - VEC(1))/(size(VEC) - 1)

TREND=+/-3 The least-squares line from the time series is removed before time filtering.

IF TREND=-1,-2 or -3, the mean, drift or least-squares line is reintroduced post-filtering, respectively. For other values of TREND nothing is done before or after filtering.

USEFFT (INPUT, OPTIONAL) logical(lgl)

On input, if USEFFT is used and is set to true, the symmetric linear filter is applied to the argument VEC by using a Fast Fourier Transform and the convolution theorem.

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if INITFFT is set to false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine SYMLIN_FILTER2 in order to sets up constants and functions for use by subroutine FFT_ROW which is called inside subroutine SYMLIN_FILTER2 (the call to INIT_FFT must have the following form:

call init_fft( shape(VEC) )

If INITFFT is set to true, the call to INIT_FFT is done inside subroutine SYMLIN_FILTER2 and a call to END_FFT is also done before leaving subroutine SYMLIN_FILTER2. This optional argument has an effect only if argument USEFFT is used with the value true.

The default is INITFFT=true .

Further Details¶

No time observations will be lost, however the first and last (size(COEF)-1)/2 time observations are affected by end effects.

If USEFFT is used with the value true, the values at both ends of the output series are computed by assuming that the input series is part of a periodic sequence of period size(VEC). Otherwise, each end of the filtered time series is estimated by truncated the symmetric linear filter coefficients array.

The symmetric linear filter coefficients (e.g. the array COEF) can be computed with the help of functions LP_COEF, LP_COEF2, HP_COEF, HP_COEF2, BD_COEF and BD_COEF2.

`subroutine symlin_filter2 ( mat, coef, trend, usefft, initfft )`¶

Purpose¶

Subroutine SYMLIN_FILTER2 performs a symmetric filtering operation on an input multi-channel time series (e.g. the argument MAT).

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

The multi-channel time series matrix to be filtered. Each column of MAT corresponds to one observation. On output, the multi-channel filtered time series are returned.

Size(MAT,2) must be greater or equal to 4.

COEF (INPUT) real(stnd), dimension(:)

The array of symmetric linear filter coefficients.

Size(COEF) must be odd, greater or equal to 3 and less or equal to size(MAT,2).

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+/-1 The means of the time series are removed before time filtering

TREND=+/-2 The drifts from the time series are removed before time filtering by using the formula: drift(:) = (MAT(:,size(MAT,2)) - MAT(:,1))/(size(MAT,2) - 1)

TREND=+/-3 The least-squares lines from the time series are removed before time filtering.

IF TREND=-1,-2 or -3, the means, drifts or least-squares lines are reintroduced post-filtering, respectively. For other values of TREND nothing is done before or after filtering.

USEFFT (INPUT, OPTIONAL) logical(lgl)

On input, if USEFFT is used and is set to true, the symmetric linear filter is applied to the argument VEC by using a Fast Fourier Transform and the convolution theorem.

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if INITFFT is set to false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine SYMLIN_FILTER2 in order to sets up constants and functions for use by subroutine FFT_ROW which is called inside subroutine SYMLIN_FILTER2 (the call to INIT_FFT must have the following form:

call init_fft( shape(MAT), dim=2_i4b )

If INITFFT is set to true, the call to INIT_FFT is done inside subroutine SYMLIN_FILTER2 and a call to END_FFT is also done before leaving subroutine SYMLIN_FILTER2. This optional argument has an effect only if argument USEFFT is used with the value true.

The default is INITFFT=true .

Further Details¶

No time observations will be lost, however the first and last (size(COEF)-1)/2 time observations are affected by end effects.

If USEFFT is used with the value true, the values at both ends of the output multi-channel time series are computed by assuming that the input multi-channel series is part of a periodic sequence of period size(VEC). Otherwise, each end of the filtered multi-channel time series is estimated by truncated the symmetric linear filter coefficients array.

The symmetric linear filter coefficients (e.g. the array COEF) can be computed with the help of functions LP_COEF, LP_COEF2, HP_COEF, HP_COEF2, BD_COEF and BD_COEF2.

`subroutine dan_filter ( vec, nsmooth, sym, trend )`¶

Purpose¶

Subroutine DAN_FILTER smooths an input time series (e.g. the argument VEC) by applying a Daniell filter (e.g. a simple moving average) of length NSMOOTH.

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

On input, the vector containing the time series to be filtered. On output, the filtered time series is returned.

Size(VEC) must be greater or equal to 4.

NSMOOTH (INPUT) integer(i4b)

The length of the Daniell filter to be applied to the time series. NSMOOTH must be odd. Moreover, NSMOOTH must be greater or equal to 3 and less or equal to size(VEC).

SYM (INPUT, OPTIONAL) real(stnd)

An optional indictor variable used to designate wether the series has an even symmetry (SYM = one), an odd symmetry (SYM = -one) or no symmetry (SYM = zero). Other values than -one, one or zero are not allowed for the optional argument SYM.

The default value for SYM is one.

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+/-1 The mean of the time series is removed before time filtering

TREND=+/-2 The drift from the time series is removed before time filtering by using the formula: drift = (VEC(size(VEC)) - VEC(1))/(size(VEC) - 1)

TREND=+/-3 The least-squares line from the time series is removed before time filtering.

IF TREND=-1,-2 or -3, the mean, drift or least-squares line is reintroduced post-filtering, respectively. For other values of TREND nothing is done before or after filtering.

Further Details¶

Subroutine DAN_FILTER smooths an input time series by applying a Daniell filter as discussed in chapter 7 of Bloomfield (1976).

This subroutine use the hypothesis of the (even or odd) symmetry of the input time series to avoid losing values from the ends of the series.

For more details and algorithm, see:

Bloomfield, P.,1976:

Fourier analysis of time series- An introduction, John Wiley and Sons, New York, Chapter 7.

`subroutine dan_filter ( mat, nsmooth, sym, trend )`¶

Purpose¶

Subroutine DAN_FILTER smooths an input multi-channel time series (the argument MAT) by applying a Daniell filter (e.g. a simple moving average) of length NSMOOTH.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

The multi-channel time series matrix to be filtered. Each column of MAT corresponds to one observation. On output, the multi-channel filtered time series are returned.

Size(MAT,2) must be greater or equal to 4.

NSMOOTH (INPUT) integer(i4b)

The length of the Daniell filter to be applied to the time series. NSMOOTH must be odd. Moreover, NSMOOTH must be greater or equal to 3 and less or equal to size(MAT,2).

SYM (INPUT, OPTIONAL) real(stnd)

An optional indictor variable used to designate wether the series has an even symmetry (SYM = one), an odd symmetry (SYM = -one) or no symmetry (SYM = zero). Other values than -one, one or zero are not allowed for the optional argument SYM.

The default value for SYM is one.

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+/-1 The means of the time series are removed before time filtering

TREND=+/-2 The drifts from the time series are removed before time filtering by using the formula: drift(:) = (MAT(:,size(MAT,2)) - MAT(:,1))/(size(MAT,2) - 1)

TREND=+/-3 The least-squares lines from the time series are removed before time filtering.

IF TREND=-1,-2 or -3, the means, drifts or least-squares lines are reintroduced post-filtering, respectively. For other values of TREND nothing is done before or after filtering.

Further Details¶

Subroutine DAN_FILTER smooths an input multi-channel time series by applying a Daniell filter as discussed in chapter 7 of Bloomfield (1976).

This subroutine may use the hypothesis of the (even or odd) symmetry of the input time series to avoid losing values from the ends of the series.

For more details and algorithm, see

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York, Chapter 7.

`subroutine moddan_filter ( vec, smooth_param, sym, trend )`¶

Purpose¶

Subroutine MODDAN_FILTER smooths an input time series (e.g. the argument VEC) by applying a sequence of modified Daniell filters.

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

On input, the vector containing the time series to be filtered. On output, the filtered time series is returned. Size(VEC) must be greater or equal to 4.

SMOOTH_PARAM (INPUT) integer(i4b), dimension(:)

The array of the half-lengths of the modified Daniell filters to be applied to the time series. All the values in SMOOTH_PARAM(:) must be greater than 0 and less than size(VEC) .

Size(SMOOTH_PARAM) must be greater or equal to 1.

SYM (INPUT, OPTIONAL) real(stnd)

An optional indictor variable used to designate wether the series has an even symmetry (SYM = one), an odd symmetry (SYM = -one) or no symmetry (SYM = zero). Other values than -one, one or zero are not allowed for the optional argument SYM.

The default value for SYM is one.

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+/-1 The mean of the time series is removed before time filtering

TREND=+/-2 The drift from the time series is removed before time filtering by using the formula: drift = (VEC(size(VEC)) - VEC(1))/(size(VEC) - 1)

TREND=+/-3 The least-squares line from the time series is removed before time filtering.

IF TREND=-1,-2 or -3, the mean, drift or least-squares line is reintroduced post-filtering, respectively. For other values of TREND nothing is done before or after filtering.

Further Details¶

Subroutine MODDAN_FILTER smooths an input time series by applying a sequence of modified Daniell filters as discussed in chapter 7 of Bloomfield (1976). This subroutine use the hypothesis of the (even or odd) symmetry of the input time series to avoid losing values from the ends of the series.

For more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York, Chapter 7.

`subroutine moddan_filter ( mat, smooth_param, sym, trend )`¶

Purpose¶

Subroutine MODDAN_FILTER smooths an input multi-channel time series (the argument MAT) by applying a sequence of modified Daniell filters.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

The multi-channel time series matrix to be filtered. Each column of MAT corresponds to one observation. On output, the multi-channel filtered time series are returned.

Size(MAT,2) must be greater or equal to 4.

SMOOTH_PARAM (INPUT) integer(i4b), dimension(:)

The array of the half-lengths of the modified Daniell filters to be applied to the time series. All the values in SMOOTH_PARAM(:) must be greater than 0 and less than size(MAT,2) .

Size(SMOOTH_PARAM) must be greater or equal to 1.

SYM (INPUT, OPTIONAL) real(stnd)

An optional indictor variable used to designate wether the series has an even symmetry (SYM = one), an odd symmetry (SYM = -one) or no symmetry (SYM = zero). Other values than -one, one or zero are not allowed for the optional argument SYM.

The default value for SYM is one.

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+/-1 The means of the time series are removed before time filtering

TREND=+/-2 The drifts from the time series are removed before time filtering by using the formula: drift(:) = (MAT(:,size(MAT,2)) - MAT(:,1))/(size(MAT,2) - 1)

TREND=+/-3 The least-squares lines from the time series are removed before time filtering.

IF TREND=-1,-2 or -3, the means, drifts or least-squares lines are reintroduced post-filtering, respectively. For other values of TREND nothing is done before or after filtering.

Further Details¶

Subroutine MODDAN_FILTER smooths an input multi-channel time series by applying a sequence of modified Daniell filters as discussed in chapter 7 of Bloomfield (1976). This subroutine may use the hypothesis of the (even or odd) symmetry of the input time series to avoid losing values from the ends of the series.

For more details and algorithm see

Bloomfield, P., 1976: Fourier analysis of time series- An introduction,

John Wiley and Sons, New York, Chapter 7.

`function extend ( vec, index, sym )`¶

Purpose¶

This function returns the INDEX-th term in the series VEC, extending it if necessary with even or odd symmetry according to the sign of SYM, which should be either plus or minus one. (Note: the value zero will result in the extended value being zero).

Arguments¶

VEC (INPUT) real(stnd), dimension(:)

On input, the vector containing the time series. If size(VEC) is zero, the extended value returned is zero.

INDEX (INPUT) integer(i4b)

On input, the index of the desired term in the time series. INDEX may be any integer.

SYM (INPUT) real(stnd)

An indictor variable used to designate wether the series has an even symmetry (SYM = one), an odd symmetry (SYM = -one) or no symmetry (SYM = zero). Other values than -one, one or zero are not allowed, however no checking is done on the SYM argument.

Further Details¶

For more details and algorithm, see

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York, Chapter 6.

`function extend ( mat, index, sym )`¶

Purpose¶

This function returns the INDEX-th term in the multi-channel series MAT, extending it if necessary with even or odd symmetry according to the sign of SYM, which should be either plus or minus one. Note: the value zero will result in the extended value being zero.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On input, the matrix containing the multi-channel time series. Each column of MAT corresponds to one observation. If size(MAT,2) is zero, the extended vector (which is dimensionned as size(MAT,1)) returned is zero.

INDEX (INPUT) integer(i4b)

On input, the index of the desired term in the multi-channel time series.

SYM (INPUT) real(stnd)

An indictor variable used to designate wether the series has an even symmetry (SYM = one), an odd symmetry (SYM = -one) or no symmetry (SYM = zero). Other values than -one, one or zero are not allowed, however no checking is done on the SYM argument.

Further Details¶

For more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York, Chapter 6.

`subroutine taper ( vec, taperp )`¶

Purpose¶

Subroutine TAPER applies a split-cosine-bell taper on an input time series (e.g. the argument VEC).

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

On input, the vector containing the time series to be tapered. On output, the tapered time series is returned.

TAPERP (INPUT) real(stnd)

The total percentage of the data to be tapered. TAPERP must be greater than zero and less or equal to one, otherwise the series is not tapered.

Further Details¶

This subroutine is adapted from Bloomfield (1976).

For more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York, Chapter 5.

`subroutine taper ( mat, taperp )`¶

Purpose¶

Subroutine TAPER applies a split-cosine-bell taper on an input multi-channel time series (the argument MAT).

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On input, the matrix containing the multi-channel time series. Each column of MAT corresponds to one observation.

TAPERP (INPUT) real(stnd)

The total percentage of the data to be tapered. TAPERP must be greater than zero and less or equal to one, otherwise the series is not tapered.

Further Details¶

This subroutine is adapted from Bloomfield (1976).

For more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York, Chapter 5.

`function data_window ( n, win, taperp )`¶

Purpose¶

Function DATA_WINDOW computes data windows used in spectral computations.

Arguments¶

N (INPUT) integer(i4b)

The size of the data window. N must be an even positive integer.

WIN (INPUT) integer(i4b)

On entry, this argument specify the form of the data window. If:

WIN=+1 The Bartlett window is used

WIN=+2 The square window is used

WIN=+3 The Welch window is used

WIN=+4 The Hann window is used

WIN=+5 The Hamming window is used

WIN=+6 A split-cosine-bell window is used

For other values of WIN, a square window is returned.

TAPERP (INPUT, OPTIONAL) real(stnd)

The total percentage of the data to be tapered if WIN=6. TAPERP must be greater than zero and less or equal to one, otherwise the default value is used.

The default is 0.2 .

Further Details¶

For more details, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York, Chapter 5.

`function estim_dof ( wk, win, smooth_param, l0, nseg, overlap )`¶

Purpose¶

Function ESTIM_DOF computes “the equivalent number of degrees of freedom” of power and cross spectrum estimates as calculated by subroutines POWER_SPECTRUM, CROSS_SPECTRUM, POWER_SPECTRUM2 and CROSS_SPECTRUM2.

Arguments¶

WK (INPUT) real(stnd), dimension(:)

On entry, this argument specify the data window used in the computations of the power and/or cross spectra.

Spectral computations are at (Size(WK)/2)+1 frequencies if the optional argument L0 is absent and are at ((Size(WK)+L0)/2)+1 frequencies if L0 is present (L0 is the number of zeros added to each segment).

Size(WK) must be greater or equal to 4 and Size(WK)+L0 must be even.

WIN (INPUT, OPTIONAL) integer(i4b)

On entry, this argument specify the form of the data window given in argument WK. If:

WIN=+1 The Bartlett window is used

WIN=+2 The square window is used

WIN=+3 The Welch window is used

WIN=+4 The Hann window is used

WIN=+5 The Hamming window is used

WIN=+6 A split-cosine-bell window is used

For other values of WIN, a message error is issued and the program is stopped.

The default is WIN=+3, e.g. the Welch window.

SMOOTH_PARAM (INPUT, OPTIONAL) integer(i4b), dimension(:)

if SMOOTH_PARAM is used, the power and/or cross spectrum have been estimated by repeated smoothing of the periodogram with modified Daniell weights.

On entry, SMOOTH_PARAM(:) gives the array of the half-lengths of the modified Daniell filters that have been applied.

All the values in SMOOTH_PARAM(:) must be greater than 0 and less than ((size(WK)+L0)/2) + 1 .

Size(SMOOTH_PARAM) must be greater or equal to 1.

L0 (INPUT, OPTIONAL) integer(i4b)

The number of zeros added to the time series (or segment) in order to obtain more finely spaced spectral estimates. L0 must be a positive integer. Moreover, Size(VEC)+L0 must be even.

The default is L0=0, e.g. no zeros are added to the time series.

NSEG (INPUT, OPTIONAL) integer(i4b)

The number of segments if the spectra have been computed by POWER_SPECTRUM2 and CROSS_SPECTRUM2 . NSEG must be a positive integer.

The segments are assumed to be independent or to overlap by one half of their length if the optional argument OVERLAP is used and is set to true. Let L = size(WK). Then, the number of segments may be computed as follows:

N/L if OVERLAP=false

(2N/L)-1 if OVERLAP=true

where N is equal to:

the length of the original time series (call it M) if this length is evenly divisible by L,

M+L-mod(M,L) if M is not evenly divisible L.

The default is NSEG=1, e.g. the time series is not segmented.

OVERLAP (INPUT, OPTIONAL) logical(lgl)

If OVERLAP is set to false, the spectrum estimates have been computed from nonoverlapping segments.

If OVERLAP is set to true, the spectrum estimates have been computed from overlapped segments (subroutines POWER_SPECTRUM2 and CROSS_SPECTRUM2 may overlap the segments by one half of their length.

The default is OVERLAP=false .

Further Details¶

The computed equivalent number of degrees of freedom must be divided by two for the zero and Nyquist frequencies.

Furthermore, the computed equivalent number of degrees of freedom is not right near the zero and Nyquist frequencies if the PSD estimates have been smoothed by modified Daniell filters. The reason is that ESTIM_DOF assumes that smoothing involves averaging independent frequency ordinates. This is true except near the zero and Nyquist frequencies where an average may contain contributions from negative frequencies, which are identical to and hence not independent of positive frequency spectral values. Thus, the number of degrees of freedom in PSD estimates near the 0 and Nyquist frequencies are as little as half the number of degrees of freedom of the spectral estimates away from these frequency extremes if the optional argument SMOOTH_PARAM is used.

For more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York, Chapter 8.

Welch, P.D., 1967:

The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE trans. on audio and electroacoustics, Vol. Au-15, 2, 70-73.

`function estim_dof2 ( wk, l0, win, nsmooth, nseg, overlap )`¶

Purpose¶

Function ESTIM_DOF2 computes “the equivalent number of degrees of freedom” of power and cross spectrum estimates as calculated by subroutines POWER_SPCTRM, CROSS_SPCTRM, POWER_SPCTRM2 and CROSS_SPCTRM2.

Arguments¶

WK (INPUT) real(stnd), dimension(:)

On entry, this argument specifies the data window used in the computations of the power and/or cross spectra.

Spectral computations are at ((Size(WK)+L0)/2)+1 frequencies (L0 is the number of zeros added to each segment).

Size(WK) must be greater or equal to 4 and Size(WK)+L0 must be even.

L0 (INPUT) integer(i4b)

The number of zeros added to the time series (or segment) in order to obtain more finely spaced spectral estimates. L0 must be a positive integer. Moreover, Size(VEC)+L0 must be even.

WIN (INPUT, OPTIONAL) integer(i4b)

On entry, this argument specify the form of the data window given in argument WK. If:

WIN=+1 The Bartlett window is used

WIN=+2 The square window is used

WIN=+3 The Welch window is used

WIN=+4 The Hann window is used

WIN=+5 The Hamming window is used

WIN=+6 A split-cosine-bell window is used

For other values of WIN, a message error is issued and the program is stopped.

The default is WIN=+3, e.g. the Welch window.

NSMOOTH (INPUT, OPTIONAL) integer(i4b)

if NSMOOTH is used, the power and/or cross spectra have been estimated by smoothing the periodogram with Daniell weights.

On entry, NSMOOTH gives the length of the Daniell filter that has been applied.

Setting NSMOOTH=0 on entry is equivalent to omit the optional argument NSMOOTH. Otherwise, NSMOOTH must be odd, greater than 2 and less or equal to ((size(WK)+L0)/2) + 1 .

NSEG (INPUT, OPTIONAL) integer(i4b)

The number of segments if the spectra have been computed by POWER_SPCTRM2 and CROSS_SPCTRM2 . NSEG must be a positive integer.

The segments are assumed to be independent or to overlap by one half of their length if the optional argument OVERLAP is used and is set to true. Let L = size(WK). Then, the number of segments may be computed as follows:

N/L if OVERLAP=false

(2N/L)-1 if OVERLAP=true

where N is equal to:

the length of the original time series (call it M) if this length is evenly divisible by L,

M+L-mod(M,L) if M is not evenly divisible L.

The default is NSEG=1, e.g. the time series is not segmented.

OVERLAP (INPUT, OPTIONAL) logical(lgl)

If OVERLAP is set to false, the spectrum estimates have been computed from nonoverlapping segments.

If OVERLAP is set to true, the spectrum estimates have been computed from overlapped segments (subroutines POWER_SPCTRUM2 and CROSS_SPCTRUM2 may overlap the segments by one half of their length.

The default is OVERLAP=false .

Further Details¶

For more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York, Chapter 8.

Welch, P.D., 1967:

The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE trans. on audio and electroacoustics, Vol. Au-15, 2, 70-73.

`subroutine comp_conflim ( edof, probtest, conlwr, conupr, testcoher )`¶

Purpose¶

Subroutine COMP_CONFLIM estimates confidence limit factors for spectral estimates and, optionally, critical value for testing the null hypothesis that squared coherency is zero.

Arguments¶

EDOF (INPUT) real(stnd)

On entry, the equivalent number of degrees of freedom of the power spectrum estimates.

PROBTEST (INPUT, OPTIONAL) real(stnd)

On entry, a probability. PROBTEST is the critical probability which is used to determine the lower and upper confidence limit factors (e.g. the optional arguments CONLWR and CONUPR ) and the critical value for testing the null hypothesis that the squared coherency is zero (e.g. the TESTCOHER optional argument). PROBTEST must verify 0. < P < 1.

The default is 0.05 .

CONLWR (OUTPUT, OPTIONAL) real(stnd)

CONUPR (OUTPUT, OPTIONAL) real(stnd)

On output, these arguments specify the lower and upper (1-PROBTEST) * 100% confidence limit factors, respectively. Multiply the PSD estimates by these constants to get the lower and upper limits of a (1-PROBTEST) * 100% confidence interval for the PSD estimates.

TESTCOHER (OUTPUT, OPTIONAL) real(stnd)

On output, this argument specifies the critical value for testing the null hypothesis that the squared coherency is zero at the PROBTEST * 100% significance level (e.g. squared coherencies less than TESTCOHER should be regarded as not significantly different from zero at the PROBTEST * 100% significance level).

`subroutine comp_conflim ( edof, probtest, conlwr, conupr, testcoher )`¶

Purpose¶

Subroutine COMP_CONFLIM estimates confidence limit factors for spectral estimates and, optionally, critical values for testing the null hypothesis that squared coherencies are zero.

Arguments¶

EDOF (INPUT) real(stnd), dimension(:)

On entry, the equivalent number of degrees of freedom of the power spectrum estimates.

PROBTEST (INPUT, OPTIONAL) real(stnd)

On entry, a probability. PROBTEST is the critical probability which is used to determine the lower and upper confidence limit factors (e.g. the optional arguments CONLWR and CONUPR ) and the critical value for testing the null hypothesis that the squared coherency is zero (e.g. the TESTCOHER optional argument). PROBTEST must verify 0. < P < 1.

The default is 0.05 .

CONLWR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

CONUPR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On output, these arguments specify the lower and upper (1-PROBTEST) * 100% confidence limit factors, respectively. Multiply the PSD estimates by these constants to get the lower and upper limits of a (1-PROBTEST) * 100% confidence interval for the PSD estimates.

CONLWR must verify: size(CONLWR) = size(EDOF) .

CONUPR must verify: size(CONUPR) = size(EDOF) .

TESTCOHER (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On output, this argument specifies the critical values for testing the null hypothesis that the squared coherencies are zero at the PROBTEST * 100% significance level (e.g. squared coherencies less than TESTCOHER(:) should be regarded as not significantly different from zero at the PROBTEST * 100% significance level).

TESTCOHER must verify: size(TESTCOHER) = size(EDOF) .

`subroutine spctrm_ratio ( edofn, edofd, lwr_ratio, upr_ratio, pinterval )`¶

Purpose¶

Subroutine SPCTRM_RATIO calculates a pointwise tolerance interval for the ratio of two estimated spectra under the assumption that the two “true” underlying spectra are the same.

Arguments¶

EDOFN (INPUT) real(stnd)

On exit, the equivalent number of degrees of freedom of the first estimated spectrum (e.g. the numerator of the ratio of the two estimated spectra).

EDOFN must be greater than zero.

EDOFD (INPUT) real(stnd)

On exit, the equivalent number of degrees of freedom of the second estimated spectrum (e.g. the denominator of the ratio of the two estimated spectra).

EDOFD must be greater than zero.

LWR_RATIO (OUTPUT) real(stnd)

UPR_RATIO (OUTPUT) real(stnd)

On output, these arguments specify the lower and upper critical ratios of the computed PINTERVAL * 100% tolerance interval for the ratio of the power spectral density estimates.

PINTERVAL (INPUT, OPTIONAL) real(stnd)

On entry, a probability. This probability is used to determine the upper and lower critical ratios of the computed tolerance interval. A PINTERVAL * 100% tolerance interval is computed and output in the two arguments LWR_RATIO and UPR_RATIO. PINTERVAL must verify: 0. < PINTERVAL < 1.

The default value is 0.90, e.g. a 90% tolerance interval is computed.

Further Details¶

For more details, see:

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford, Chapter 4.

`subroutine spctrm_ratio ( edofn, edofd, lwr_ratio, upr_ratio, pinterval )`¶

Purpose¶

Subroutine SPCTRM_RATIO calculates pointwise tolerance intervals for the ratio of two estimated spectra under the assumption that the two “true” underlying spectra are the same.

Arguments¶

EDOFN (INPUT) real(stnd), dimension(:)

On exit, the equivalent number of degrees of freedom of the first estimated spectrum (e.g. the numerator of the ratio of the two estimated spectra).

Elements of EDOFN(:) must be greater than zero.

EDOFD (INPUT) real(stnd), dimension(:)

On exit, the equivalent number of degrees of freedom of the second estimated spectrum (e.g. the denominator of the ratio of the two estimated spectra). Elements of EDOFD(:) must be greater than zero.

EDOFD must verify: size(EDOFD) = size(EDOFN) .

LWR_RATIO (OUTPUT) real(stnd), dimension(:)

UPR_RATIO (OUTPUT) real(stnd), dimension(:)

On output, these arguments specify the lower and upper critical ratios of the computed PINTERVAL * 100% tolerance interval for the ratio of the power spectral density estimates.

LWR_RATIO must verify: size(LWR_RATIO) = size(EDOFN) .

UPR_RATIO must verify: size(UPR_RATIO) = size(EDOFN) .

PINTERVAL (INPUT, OPTIONAL) real(stnd)

On entry, a probability. This probability is used to determine the upper and lower critical ratios of the computed tolerance interval. A PINTERVAL * 100% tolerance interval is computed and output in the two arguments LWR_RATIO and UPR_RATIO. PINTERVAL must verify: 0. < PINTERVAL < 1.

The default value is 0.90, e.g. a 90% tolerance interval is computed .

Further Details¶

For more details, see:

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford, Chapter 4.

`subroutine spctrm_ratio2 ( psvecn, psvecd, edofn, edofd, prob, min_ratio, max_ratio, prob_min_ratio, prob_max_ratio )`¶

Purpose¶

Subroutine SPCTRM_RATIO2 calculates a conservative critical probability value (e.g. p-value) for testing the hypothesis of a common spectrum for two estimated spectra (e.g. the arguments PSVECN, PSVECD). This conservative critical probability value is computed from the minimim and maximum values of the ratio of the two estimated spectra and the associated probabilities of obtaining, respectively, a value less (for the minimum ratio) and higher (for the maximum ratio) than attained under the null hypothesis of a common spectrum for the two time series.

Arguments¶

PSVECN (INPUT) real(stnd), dimension(:)

On entry, a real vector containing the Power Spectral Density (PSD) estimates of the first time series (e.g. the numerator of the ratio of the two estimated spectra).

All elements in PSVECN(:) must be greater or equal to zero and size(PSVECN) must be greater or equal to 2.

PSVECD (INPUT) real(stnd), dimension(:)

On entry, a real vector containing the Power Spectral Density (PSD) estimates of the second time series (e.g. the denominator of the ratio of the two estimated spectra). All elements in PSVECD(:) must be greater than zero and size(PSVECD) must be greater or equal to 2.

PSVECD must also verify: size(PSVECD) = size(PSVECN) .

EDOFN (INPUT) real(stnd)

On exit, the equivalent number of degrees of freedom of the first estimated spectrum (e.g. the numerator of the ratio of the two estimated spectra).

EDOFN must be greater than zero.

EDOFD (INPUT) real(stnd)

On exit, the equivalent number of degrees of freedom of the second estimated spectrum (e.g. the denominator of the ratio of the two estimated spectra).

EDOFD must be greater than zero.

PROB (OUTPUT) real(stnd)

On exit, the conservative critical probability value (e.g. p-value) computed under the hypothesis that the two “true” underlying spectra are the same. See the description of the PROB_MIN_RATIO and PROB_MAX_RATIO optional arguments for more details.

MIN_RATIO (OUTPUT, OPTIONAL) real(stnd)

MAX_RATIO (OUTPUT, OPTIONAL) real(stnd)

On output, these arguments give, respectively, the minimum and maximum values of the ratio of the two PSD estimates.

PROB_MIN_RATIO (OUTPUT, OPTIONAL) real(stnd)

PROB_MAX_RATIO (OUTPUT, OPTIONAL) real(stnd)

On output, these arguments give, respectively, the probabilities of obtaining a smaller value of the minimum ratio (e.g. the argument MIN_RATIO) and a greater value of the maximum ratio (e.g. the argument MAX_RATIO) under the null hypothesis that the two “true” underlying spectra are the same.

The PROB argument is computed as 2 * min(PROB_MIN_RATIO,PROB_MAX_RATIO).

Further Details¶

This statistical test relies on the assumptions that the different spectral ordinates have the same sampling distribution and are independent of each other for each series. This means, in particular, that the spectral ordinates corresponding to the zero and Nyquist frequencies must be excluded from the PSVECN and PSVECD vectors before calling SPCTRM_RATIO2 and that the two estimated spectra have not been obtained by smoothing the periodogram in the frequency domain.

It is also assumed that the PSVECN and PSVECD realizations are independent.

For more details, see:

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford, Chapter 4.

`subroutine spctrm_ratio2 ( psmatn, psmatd, edofn, edofd, prob, min_ratio, max_ratio, prob_min_ratio, prob_max_ratio )`¶

Purpose¶

Subroutine SPCTRM_RATIO2 calculates conservative critical probability values (e.g. p-values) for testing the hypothesis of a common spectrum for the elements of two estimated multi-channel spectra (e.g. the arguments PSMATN, PSMATD).

These conservative critical probability values are computed from the minimim and maximum values of the ratio of the two estimated multi-channel spectra and the associated probabilities of obtaining, respectively, a value less (for the minimum ratio) and higher (for the maximum ratio) than attained under the null hypothesis of a common spectrum for the two multi-channel time series.

Arguments¶

PSMATN (INPUT) real(stnd), dimension(:,:)

On entry, a real matrix containing the Power Spectral Density (PSD) estimates of the first multi-channel time series (e.g. the numerator of the ratio of the two estimated multi-channel spectra). Each row of the real matrix PSMATN contains the estimated spectrum of the corresponding “row” of the first multi-channel times series.

All elements in PSMATN(:,:) must be greater or equal to zero and size(PSMATN,2) must be greater or equal to 2.

PSMATD (INPUT) real(stnd), dimension(:,:)

On entry, a real matrix containing the Power Spectral Density (PSD) estimates of the second multi-channel time series (e.g. the denominator of the ratio of the two estimated multi-channel spectra). Each row of the real matrix PSMATD contains the estimated spectrum of the corresponding “row” of the second multi-channel times series. All elements in PSMATD(:,:) must be greater than zero and size(PSMATD,2) must be greater or equal to 2.

PSMATD must also verify:

size(PSMATD,1) = size(PSMATN,1) ,

size(PSMATD,2) = size(PSMATN,2) .

EDOFN (INPUT) real(stnd)

On exit, the equivalent number of degrees of freedom of the first estimated spectrum (e.g. the numerator of the ratio of the two estimated spectra).

EDOFN must be greater than zero.

EDOFD (INPUT) real(stnd)

On exit, the equivalent number of degrees of freedom of the second estimated spectrum (e.g. the denominator of the ratio of the two estimated spectra).

EDOFD must be greater than zero.

PROB (OUTPUT) real(stnd), dimension(:)

On exit, the conservative critical probability values (e.g. p-values) computed under the hypothesis that the two “true” underlying multi-channel spectra are the same. See the description of the PROB_MIN_RATIO and PROB_MAX_RATIO optional arguments for more details.

PROB must verify: size(PROB) = size(PSMATN,1) .

MIN_RATIO (OUTPUT, OPTIONAL) real(stnd), dimension(:)

MAX_RATIO (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On output, these arguments give, respectively, the minimum and maximum values of the ratio of the two multi-channel PSD estimates.

MIN_RATIO must verify: size(MIN_RATIO) = size(PSMATN,1) .

MAX_RATIO must verify: size(MAX_RATIO) = size(PSMATN,1) .

PROB_MIN_RATIO (OUTPUT, OPTIONAL) real(stnd), dimension(:)

PROB_MAX_RATIO (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On output, these arguments give, respectively, the probabilities of obtaining a smaller value of the minimum ratio (e.g. the argument MIN_RATIO) and a greater value of the maximum ratio (e.g. the argument MAX_RATIO) under the null hypothesis that the two “true” underlying multi-channel spectra are the same. The PROB(:) argument is calculated as 2 * min(PROB_MIN_RATIO(:),PROB_MAX_RATIO(:)).

PROB_MIN_RATIO must verify: size(PROB_MIN_RATIO) = size(PSMATN,1) .

PROB_MAX_RATIO must verify: size(PROB_MAX_RATIO) = size(PSMATN,1) .

Further Details¶

This statistical test relies on the assumptions that the different spectral ordinates have the same sampling distribution and are independent of each other for each series. This means, in particular, that the spectral ordinates corresponding to the zero and Nyquist frequencies must be excluded from the PSMATN and PSMATD matrices before calling SPCTRM_RATIO2 and that the two estimated multi-channel spectra have not been obtained by smoothing the periodogram in the frequency domain.

It is also assumed that the PSMATN and PSMATD realizations are independent.

For more details, see:

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford, Chapter 4.

`subroutine spctrm_ratio3 ( psvecn, psvecd, edofn, edofd, chi2_stat, prob )`¶

Purpose¶

Subroutine SPCTRM_RATIO3 calculates an approximate critical probability value (e.g. p-value) for testing the hypothesis of a common spectrum for two estimated spectra (e.g. the arguments PSVECN, PSVECD). This approximate critical probability value is derived from the following CHI2 statistic :

CHI2_STAT = ( 2/EDOFN + 2/EDOFD )**(-1) [ sum k=1 to nf ] log( PSVECN(k) / PSVECD(k) )**(2)

where nf = size(PSVECN) = size(PSVECD). In order to derive an approximate critical probability value, it is assumed that CHI2_STAT has an approximate CHI2 distribution with nf degrees of freedom.

Arguments¶

PSVECN (INPUT) real(stnd), dimension(:)

On entry, a real vector containing the Power Spectral Density (PSD) estimates of the first time series (e.g. the numerator of the ratio of the two estimated spectra).

All elements in PSVECN(:) must be greater than zero and size(PSVECN) must be greater or equal to 2.

PSVECD (INPUT) real(stnd), dimension(:)

On entry, a real vector containing the Power Spectral Density (PSD) estimates of the second time series (e.g. the denominator of the ratio of the two estimated spectra).

All elements in PSVECD(:) must be greater than zero and size(PSVECD) must be greater or equal to 2.

PSVECD must verify: size(PSVECD) = size(PSVECN) .

EDOFN (INPUT) real(stnd)

On exit, the equivalent number of degrees of freedom of the first estimated spectrum (e.g. the numerator of the ratio of the two estimated spectra).

EDOFN must be greater than one.

EDOFD (INPUT) real(stnd)

On exit, the equivalent number of degrees of freedom of the second estimated spectrum (e.g. the denominator of the ratio of the two estimated spectra).

EDOFD must be greater than one.

CHI2_STAT (OUTPUT) real(stnd)

On output, the CHI2 statistic which is assumed to follow a CHI2 distribution with size(PSVECN) degrees of freedom under the null hypothesis of a common spectrum.

PROB (OUTPUT) real(stnd)

On exit, the aproximate critical probability value (e.g. p-value) computed under the hypothesis that the two “true” underlying spectra are the same. PROB is calculated as the probability of obtaining a value greater or equal to CHI2_STAT under the hypothesis of a common spectrum for the two series.

Further Details¶

This statistical test relies on the assumptions that the different spectral ordinates have the same sampling distribution and are independent of each other for each time series. This means, in particular, that the spectral ordinates corresponding to the zero and Nyquist frequencies must be excluded from the PSVECN and PSVECD vectors before calling SPCTRM_RATIO3 and that the two estimated spectra have not been obtained by smoothing the periodogram in the frequency domain.

It is also assumed that the PSVECN and PSVECD realizations are independent.

`subroutine spctrm_ratio3 ( psmatn, psmatd, edofn, edofd, chi2_stat, prob )`¶

Purpose¶

Subroutine SPCTRM_RATIO3 calculates approximate critical probability values (e.g. p-values) for testing the hypothesis of a common spectrum for two estimated multi-channel spectra (eg the arguments PSMATN, PSMATD). Thess approximate critical probability values are derived from the following CHI2 statistics :

CHI2_STAT(:n) = ( 2/EDOFN + 2/EDOFD )**(-1) [ sum k=1 to nf ] log( PSMATN(:n,k) / PSMATD(:n,k) )**(2)

where n = size(PSMATN,1) = size(PSMATD,1) = size(CHI2_STAT) and nf = size(PSMATN,2) = size(PSMATD,2). In order to derive approximate critical probability values, it is assumed that each element of CHI2_STAT(:n) has an approximate CHI2 distribution with nf degrees of freedom.

Arguments¶

PSMATN (INPUT) real(stnd), dimension(:,:)

On entry, a real matrix containing the Power Spectral Density (PSD) estimates of the first multi-channel time series (e.g. the numerator of the ratio of the two estimated multi-channel spectra). Each row of the real matrix PSMATN contains the estimated spectrum of the corresponding “row” of the first multi-channel times series.

All elements in PSMATN(:,:) must be greater than zero and size(PSMATN,2) must be greater or equal to 2.

PSMATD (INPUT) real(stnd), dimension(:,:)

On entry, a real matrix containing the Power Spectral Density (PSD) estimates of the second multi-channel time series (e.g. the denominator of the ratio of the two estimated multi-channel spectra). Each row of the real matrix PSMATD contains the estimated spectrum of the corresponding “row” of the second multi-channel times series. All elements in PSMATD(:,:) must be greater than zero and size(PSMATD,2) must be greater or equal to 2.

PSMATD must also verify:

size(PSMATD,1) = size(PSMATN,1) ,

size(PSMATD,2) = size(PSMATN,2) .

EDOFN (INPUT) real(stnd)

On exit, the equivalent number of degrees of freedom of the first estimated spectrum (e.g. the numerator of the ratio of the two estimated spectra).

EDOFN must be greater than one.

EDOFD (INPUT) real(stnd)

On exit, the equivalent number of degrees of freedom of the second estimated spectrum (e.g. the denominator of the ratio of the two estimated spectra).

EDOFD must be greater than one.

CHI2_STAT (OUTPUT) real(stnd), dimension(:)

On output, the CHI2 statistics which are assumed to follow a CHI2 distribution with size(PSMATN,2) degrees of freedom under the null hypothesis of a common spectrum.

CHI2_STAT must verify: size(CHI2_STAT) = size(PSMATN,1) .

PROB (OUTPUT) real(stnd), dimension(:)

On exit, the aproximate critical probability values (e.g. p-values) computed under the hypothesis that the two “true” underlying multi-channel spectra are the same. Each element of PROB(:) is calculated as the probability of obtaining a value greater or equal to the corresponding element of CHI2_STAT(:) under the hypothesis of a common spectrum for the two (single) series.

PROB must verify: size(PROB) = size(PSMATN,1) .

Further Details¶

This statistical test relies on the assumptions that the different spectral ordinates have the same sampling distribution and are independent of each other for each time series. This means, in particular, that the spectral ordinates corresponding to the zero and Nyquist frequencies must be excluded from the PSMATN and PSMATD matrices before calling SPCTRM_RATIO3 and that the two estimated multi-channel spectra have not been obtained by smoothing the periodogram in the frequency domain.

It is also assumed that the PSMATN and PSMATD realizations are independent.

`subroutine spctrm_ratio4 ( psvecn, psvecd, edofn, edofd, range_stat, prob )`¶

Purpose¶

Subroutine SPCTRM_RATIO4 calculates an approximate critical probability value (e.g. p-value) for testing the hypothesis of a common shape for two estimated spectra (e.g. the arguments PSVECN, PSVECD). This approximate critical probability value is derived from the following RANGE statistic :

RANGE_STAT = ( 2/EDOFN + 2/EDOFD )**(-1/2) * ( maxval( logratio(:nf) ) - minval( logratio(:nf) ) )

where nf = size(PSVECN) = size(PSVECD) and logratio(:nf) is given as

logratio(:nf) = log( PSVECN(:nf) / PSVECD(:nf) )

In order to derive an approximate critical probability value, it is assumed that the elements of the vector logratio(:nf) are independent and follow approximately a normal distribution with mean (1/EDOFN) - (1/EDOFD) and variance (2/EDOFN) + (2/EDOFD). Than, the distribution of the statistic RANGE_STAT may be approximated by the distribution function of the range of nf independent normal random variables with mean and variance as specified above.

Arguments¶

PSVECN (INPUT) real(stnd), dimension(:)

On entry, a real vector containing the Power Spectral Density (PSD) estimates of the first time series (e.g. the numerator of the ratio of the two estimated spectra).

All elements in PSVECN(:) must be greater than zero and size(PSVECN) must be greater or equal to 2.

PSVECD (INPUT) real(stnd), dimension(:)

On entry, a real vector containing the Power Spectral Density (PSD) estimates of the second time series (e.g. the denominator of the ratio of the two estimated spectra).

All elements in PSVECD(:) must be greater than zero and size(PSVECD) must be greater or equal to 2.

PSVECD must also verify: size(PSVECD) = size(PSVECN) .

EDOFN (INPUT) real(stnd)

On exit, the equivalent number of degrees of freedom of the first estimated spectrum (e.g. the numerator of the ratio of the two estimated spectra).

EDOFN must be greater than one.

EDOFD (INPUT) real(stnd)

On exit, the equivalent number of degrees of freedom of the second estimated spectrum (e.g. the denominator of the ratio of the two estimated spectra).

EDOFD must be greater than one.

RANGE_STAT (OUTPUT) real(stnd)

On output, the range statistic which is assumed to follow the distribution of the range of nf=size(PSVECN) independent standard normal variates under the null hypothesis of a common shape spectrum.

PROB (OUTPUT) real(stnd)

On exit, the aproximate critical probability value (e.g. p-value) computed under the hypothesis that the two “true” underlying spectra have the same shape. PROB is calculated as the probability of obtaining a value greater or equal to RANGE_STAT under the hypothesis of a common shape spectrum for the two series.

Further Details¶

This statistical test relies on the assumptions that the different spectral ordinates have the same sampling distribution and are independent of each other for each time series. This means, in particular, that the spectral ordinates corresponding to the zero and Nyquist frequencies must be excluded from the PSVECN and PSVECD vectors before calling SPCTRM_RATIO4 and that the two estimated spectra have not been obtained by smoothing the periodogram in the frequency domain.

It is also assumed that the PSVECN and PSVECD realizations are independent.

For more details, see:

Coates, D.S., and Diggle, P.J., 1986:

Tests for comparing two estimated spectral densities. J. Time series Analysis, Vol. 7, pp. 7-20 .

Potscher, B.,M., and Reschenhofer, E., 1988:

Discriminating between two spectral densities in case of replicated observations. J. Time series Analysis, Vol. 9, pp. 221-224 .

Potscher, B.,M., and Reschenhofer, E., 1989:

Distribution of the Coates-Diggle test statistic in case of replicated observations. Statistics, Vol. 20, pp. 417-421 .

`subroutine spctrm_ratio4 ( psmatn, psmatd, edofn, edofd, range_stat, prob )`¶

Purpose¶

Subroutine SPCTRM_RATIO4 calculates approximate critical probability values (e.g. p-values) for testing the hypothesis of a common shape for two estimated multi-channel spectra (e.g. the arguments PSMATN, PSMATD). These approximate critical probability values are derived from the following range statistics :

RANGE_STAT(:n) = ( 2/EDOFN + 2/EDOFD )**(-1/2) * ( maxval( logratio(:n,:nf), dim=2 ) - minval( logratio(:n,:nf), dim=2 ) )

where n = size(PSMATN,1) = size(PSMATD,1), nf = size(PSMATN,2) = size(PSMATD,2) and logratio(:n,:nf) is given as

logratio(:n,:nf) = log( PSMATN(:n,:nf) / PSMATD(:n,:nf) )

In order to derive approximate critical probability values, it is assumed that the elements of the vectors logratio(i,:nf), for i=1 to n, are independent and follow approximately a normal distribution with mean (1/EDOFN) - (1/EDOFD) and variance (2/EDOFN) + (2/EDOFD). Than, the distribution of the statistics RANGE_STAT(:n) may be approximated by the distribution function of the range of nf independent normal random variables with mean and variance as specified above.

Arguments¶

PSMATN (INPUT) real(stnd), dimension(:,:)

On entry, a real matrix containing the Power Spectral Density (PSD) estimates of the first multi-channel time series (e.g. the numerator of the ratio of the two estimated multi-channel spectra). Each row of the real matrix PSMATN contains the estimated spectrum of the corresponding “row” of the first multi-channel times series.

All elements in PSMATN(:,:) must be greater than zero and size(PSMATN,2) must be greater or equal to 2.

PSMATD (INPUT) real(stnd), dimension(:,:)

On entry, a real matrix containing the Power Spectral Density (PSD) estimates of the second multi-channel time series (e.g. the denominator of the ratio of the two estimated multi-channel spectra). Each row of the real matrix PSMATD contains the estimated spectrum of the corresponding “row” of the second multi-channel times series.

All elements in PSMATD(:,:) must be greater than zero and size(PSMATD,2) must be greater or equal to 2.

PSMATD must also verify:

size(PSMATD,1) = size(PSMATN,1) ,

size(PSMATD,2) = size(PSMATN,2) .

EDOFN (INPUT) real(stnd)

On exit, the equivalent number of degrees of freedom of the first estimated spectrum (e.g. the numerator of the ratio of the two multi-channel estimated spectra).

EDOFN must be greater than one.

EDOFD (INPUT) real(stnd)

On exit, the equivalent number of degrees of freedom of the second estimated spectrum (e.g. the denominator of the ratio of the two multi-channel estimated spectra).

EDOFD must be greater than one.

RANGE_STAT (OUTPUT) real(stnd), dimension(:)

On output, the range statistics which are assumed to follow the distribution of the range of nf=size(PSMATN,2) independent standard normal variates under the null hypothesis of a common shape spectrum.

RANGE_STAT must verify: size(RANGE_STAT) = size(PSMATN,1) .

PROB (OUTPUT) real(stnd), dimension(:)

On exit, the aproximate critical probability values (e.g. p-values) computed under the hypothesis that the two “true” underlying multi-channel spectra have the same shape.

Each element of PROB(:) is calculated as the probability of obtaining a value greater or equal to the corresponding element of RANGE_STAT(:) under the hypothesis of a common shape spectrum for the two multi-channel series.

Further Details¶

This statistical test relies on the assumptions that the different spectral ordinates have the same sampling distribution and are independent of each other for each multi-channel time series. This means, in particular, that the spectral ordinates corresponding to the zero and Nyquist frequencies must be excluded from the PSMATN and PSMATD matrices before calling SPCTRM_RATIO4 and that the two estimated multi-channel spectra have not been obtained by smoothing the periodogram in the frequency domain.

It is also assumed that the PSMATN and PSMATD multi-channel realizations are independent.

For more details, see:

Coates, D.S., and Diggle, P.J., 1986:

Tests for comparing two estimated spectral densities. J. Time series Analysis, Vol. 7, pp. 7-20 .

Potscher, B.,M., and Reschenhofer, E., 1988:

Discriminating between two spectral densities in case of replicated observations. J. Time series Analysis, Vol. 9, pp. 221-224 .

Potscher, B.,M., and Reschenhofer, E., 1989:

Distribution of the Coates-Diggle test statistic in case of replicated observations. Statistics, Vol. 20, pp. 417-421 .

`subroutine spctrm_diff ( psvec1, psvec2, ks_stat, prob, nrep, norm, initseed )`¶

Purpose¶

Subroutine SPCTRM_DIFF calculates an approximate critical probability value (e.g. p-value) for testing the hypothesis of a common shape for two estimated spectra (e.g. the arguments PSVEC1 and PSVEC2). This approximate critical probability value is derived from the following Kolmogorov-Smirnov statistic (e.g. the KS_STAT argument) :

D = [ sup m=1 to nf ] | F1(m) - F2(m) |

where nf = size(PSVEC1) = size(PSVEC2), F1(:) and F2(:) are the normalized cumulative periodograms computed from the estimated spectra PSVEC1(:) and PSVEC2(:). The distribution of D under the null hypothesis of a common shape for the spectra of the two series is approximated by calculating D for some large number (e.g. the NREP argument) of random interchanges of periodogram ordinates at each frequency for the two estimated spectra (e.g. the arguments PSVEC1(:) and PSVEC2(:)).

Arguments¶

PSVEC1 (INPUT) real(stnd), dimension(:)

On entry, a real vector containing the Power Spectral Density (PSD) estimates of the first time series.

size(PSVECN) must be greater or equal to 10.

PSVEC2 (INPUT) real(stnd), dimension(:)

On entry, a real vector containing the Power Spectral Density (PSD) estimates of the second time series.

PSVEC2 must verify: size(PSVEC2) = size(PSVEC1) .

KS_STAT (OUTPUT) real(stnd)

On output, the Kolmogorov-Smirnov statistic.

PROB (OUTPUT) real(stnd)

On exit, the aproximate critical probability value (e.g. p-value) computed under the hypothesis that the two “true” underlying spectra have the same shape. PROB is calculated as the probability of obtaining a value greater or equal to KS_STAT under the hypothesis of a common shape for the spectra of the two series.

NREP (INPUT, OPTIONAL) integer(i4b)

On entry, when argument NREP is present, NREP specifies the number of random interchanges of periodogram ordinates at each frequency in order to approximate the randomization distribution of KS_STAT under the null hypothesis.

The default is 99.

NORM (INPUT, OPTIONAL) logical(lgl)

On entry, if:

NORM=true, KS_STAT is calculated from normalized cumulative periodograms.

NORM=false KS_STAT is calculated from non-normalized cumulative periodograms. In that case the null hypothesis is that the spectra of the two time series is the same.

The default is NORM=true.

INITSEED (INPUT, OPTIONAL) logical(lgl)

On entry, if INITSEED=true, a call to RANDOM_SEED_() without arguments is done in the subroutine, in order to initiate a non-repeatable reset of the seed used by the STATPACK random generator.

The default is INITSEED=false.

Further Details¶

This statistical randomization test relies on the assumptions that the different spectral ordinates have the same sampling distribution and are independent of each other. This means, in particular, that the spectral ordinates corresponding to the zero and Nyquist frequencies must be excluded from the PSVEC1 and PSVEC2 vectors before calling SPCTRM_DIFF and that the two estimated spectra have not been obtained by smoothing the periodogram in the frequency domain.

This randomization test could only be used to compare two periodograms or two spectral estimates computed as the the average of, say, r periodograms for each time series.

For more details, see:

Diggle, P.J., and Fisher, N.I., 1991:

Nonparametric comparison of cumulative periodograms, Applied Statistics, Vol. 40, No 3, pp. 423-434.

`subroutine spctrm_diff ( psmat1, psmat2, ks_stat, prob, nrep, norm, initseed )`¶

Purpose¶

Subroutine SPCTRM_DIFF calculates approximate critical probability values (e.g. p-value) for testing the hypothesis of a common shape for two estimated multi-channel spectra (e.g. the arguments PSMAT1 and PSMAT2). These approximate critical probability values are derived from the following Kolmogorov-Smirnov statistics (e.g. the KS_STAT(:) argument) :

D(j) = [ sup m=1 to nf ] | F1(j,m) - F2(j,m) | for j=1, … , size(PSMAT1,1)

where nf = size(PSMAT1,2) = size(PSMAT2,2), F1(:,:) and F2(:,:) are the normalized cumulative periodograms computed from the estimated spectra PSMAT1(:,:) and PSMAT2(:,:). The distribution of D under the null hypothesis of a common shape for the spectra of the two multi-channel series is approximated by calculating D for some large number (e.g. the NREP argument) of random interchanges of periodogram ordinates at each frequency for the two estimated multi-channel spectra (e.g. the arguments PSMAT1(:,:) and PSMAT2(:,:)).

Arguments¶

PSMAT1 (INPUT) real(stnd), dimension(:,:)

On entry, a real matrix containing the Power Spectral Density (PSD) estimates of the first multi-channel time series. Each row of the real matrix PSMAT1 contains the estimated spectrum of the corresponding “row” of the first multi-channel times series.

size(PSMATN,2) must be greater or equal to 10.

PSMAT2 (INPUT) real(stnd), dimension(:,:)

On entry, a real matrix containing the Power Spectral Density (PSD) estimates of the second multi-channel time series. Each row of the real matrix PSMAT2 contains the estimated spectrum of the corresponding “row” of the second multi-channel times series.

PSMAT2 must verify:

size(PSMAT2,1) = size(PSMAT1,1) ,

size(PSMAT2,2) = size(PSMAT1,2) .

KS_STAT (OUTPUT) real(stnd), dimension(:)

On output, the Kolmogorov-Smirnov statistics for the multi-channel times series .

KS_STAT must verify: size(KS_STAT) = size(PSMAT1,1) .

PROB (OUTPUT) real(stnd), dimension(:)

On exit, the aproximate critical probability values (e.g. p-values) computed under the hypothesis that the two “true” underlying multi-channel spectra have the same shape.

PROB is calculated as the probability of obtaining a value greater or equal to KS_STAT under the hypothesis of a common shape for the spectra of the two series.

PROB must verify: size(PROB) = size(PSMAT1,1) .

NREP (INPUT, OPTIONAL) integer(i4b)

On entry, when argument NREP is present, NREP specifies the number of random interchanges of periodogram ordinates at each frequency in order to approximate the randomization distribution of KS_STAT under the null hypothesis.

The default is 99.

NORM (INPUT, OPTIONAL) logical(lgl)

On entry, if:

NORM=true, KS_STAT is calculated from normalized cumulative periodograms.

NORM=false KS_STAT is calculated from non-normalized cumulative periodograms. In that case the null hypothesis is that the spectra of the two multi-channel time series is the same.

The default is NORM=true.

INITSEED (INPUT, OPTIONAL) logical(lgl)

On entry, if INITSEED=true, a call to RANDOM_SEED_() without arguments is done in the subroutine, in order to initiate a non-repeatable reset of the seed used by the STATPACK random generator.

The default is INITSEED=false.

Further Details¶

This statistical randomization test relies on the assumptions that the different spectral ordinates have the same sampling distribution and are independent of each other. This means, in particular, that the spectral ordinates corresponding to the zero and Nyquist frequencies must be excluded from the PSMAT1 and PSMAT2 matrices before calling SPCTRM_DIFF and that the two estimated multi-channel spectra have not been obtained by smoothing the periodograms in the frequency domain.

This randomization test could only be used to compare two periodograms or two spectral estimates computed as the the average of, say, r periodograms for each time series.

For more details see:

Diggle, P.J., and Fisher, N.I., 1991:

Nonparametric comparison of cumulative periodograms, Applied Statistics, Vol. 40, No 3, pp. 423-434.

`subroutine spctrm_diff2 ( psvec1, psvec2, chi2_stat, prob, nrep, initseed )`¶

Purpose¶

Subroutine SPCTRM_DIFF2 calculates an approximate critical probability value (e.g. p-value) for testing the hypothesis of a common underlying spectrum for the two estimated spectra (e.g. the arguments PSVEC1 and PSVEC2). This approximate critical probability value is derived from the following CHI2 statistic (e.g. the CHI2_STAT argument) :

CHI2_STAT = (1/nf) [ sum k=1 to nf ] log( PSVEC1(k) / PSVEC2(k) )**(2)

where nf = size(PSVEC1) = size(PSVEC2). The distribution of CHI2_STAT under the null hypothesis of a common spectrum for the two series is approximated by calculating CHI2_STAT for some large number (e.g. the NREP argument) of random interchanges of periodogram ordinates at each frequency for the two estimated spectra (e.g. the arguments PSVEC1(:) and PSVEC2(:)).

Arguments¶

PSVEC1 (INPUT) real(stnd), dimension(:)

On entry, a real vector containing the Power Spectral Density (PSD) estimates of the first time series.

All elements in PSVEC1(:) must be greater than zero. size(PSVECN) must be greater or equal to 10.

PSVEC2 (INPUT) real(stnd), dimension(:)

On entry, a real vector containing the Power Spectral Density (PSD) estimates of the second time series.

All elements in PSVEC2(:) must be greater than zero.

PSVEC2 must verify: size(PSVEC2) = size(PSVEC1) .

CHI2_STAT (OUTPUT) real(stnd)

On output, the CHI2 statistic.

PROB (OUTPUT) real(stnd)

On exit, the aproximate critical probability value (e.g. p-value) computed under the hypothesis that the two “true” underlying spectra are the same.

PROB is calculated as the probability of obtaining a value greater or equal to CHI2_STAT under the hypothesis of a common spectrum for the two series.

NREP (INPUT, OPTIONAL) integer(i4b)

On entry, when argument NREP is present, NREP specifies the number of random interchanges of periodogram ordinates at each frequency in order to approximate the randomization distribution of CHI2_STAT under the null hypothesis.

The default is 99.

INITSEED (INPUT, OPTIONAL) logical(lgl)

On entry, if INITSEED=true, a call to RANDOM_SEED_() without arguments is done in the subroutine, in order to initiates a non-repeatable reset of the seed used by the STATPACK random generator.

The default is INITSEED=false.

Further Details¶

This statistical randomization test relies on the assumptions that the different spectral ordinates have the same sampling distribution and are independent of each other. This means, in particular, that the spectral ordinates corresponding to the zero and Nyquist frequencies must be excluded from the PSVEC1 and PSVEC2 vectors before calling SPCTRM_DIFF2 and that the two estimated spectra have not been obtained by smoothing the periodograms in the frequency domain.

This randomization test could only be used to compare two periodograms or two spectral estimates computed as the the average of, say, r periodograms for each time series.

Finally, none of the spectral estimates must be zero.

For more details see:

Diggle, P.J., and Fisher, N.I., 1991:

Nonparametric comparison of cumulative periodograms, Applied Statistics, Vol. 40, No 3, pp. 423-434.

`subroutine spctrm_diff2 ( psmat1, psmat2, chi2_stat, prob, nrep, initseed )`¶

Purpose¶

Subroutine SPCTRM_DIFF2 calculates approximate critical probability values (e.g. p-value) for testing the hypothesis of a common spectrum for two estimated multi-channel spectra (e.g. the arguments PSMAT1 and PSMAT2). These approximate critical probability values are derived from the following CHI2 statistics (e.g. the CHI2_STAT(:) argument) :

CHI2_STAT(:n) = ( 1/nf ) [ sum k=1 to nf ] log( PSMAT1(:n,k) / PSMAT2(:n,k) )**(2)

where n = size(PSMAT1,1) = size(PSMAT2,1) = size(CHI2_STAT) and nf = size(PSMAT2,2) = size(PSMAT2,2).

The distribution of CHI2_STAT under the null hypothesis of a common spectrum for the spectra of the two multi-channel series is approximated by calculating CHI2_STAT for some large number (e.g. the NREP argument) of random interchanges of periodogram ordinates at each frequency for the two estimated multi-channel spectra (e.g. the arguments PSMAT1(:,:) and PSMAT2(:,:)).

Arguments¶

PSMAT1 (INPUT) real(stnd), dimension(:,:)

On entry, a real matrix containing the Power Spectral Density (PSD) estimates of the first multi-channel time series. Each row of the real matrix PSMAT1 contains the estimated spectrum of the corresponding “row” of the first multi-channel times series.

All elements in PSMAT1(:,:) must be greater than zero. size(PSMATN,2) must be greater or equal to 10.

PSMAT2 (INPUT) real(stnd), dimension(:,:)

On entry, a real matrix containing the Power Spectral Density (PSD) estimates of the second multi-channel time series. Each row of the real matrix PSMAT2 contains the estimated spectrum of the corresponding “row” of the second multi-channel times series.

All elements in PSMAT2(:,:) must be greater than zero.

PSMAT2 must verify:

size(PSMAT2,1) = size(PSMAT1,1) ,

size(PSMAT2,2) = size(PSMAT2,2) .

CHI2_STAT (OUTPUT) real(stnd), dimension(:)

On output, the CHI2 statistics computed from the multi-channel times series .

CHI2_STAT must verify: size(CHI2_STAT) = size(PSMAT1,1) .

PROB (OUTPUT) real(stnd), dimension(:)

On exit, the aproximate critical probability values (e.g. p-values) computed under the hypothesis that the two “true” underlying multi-channel spectra are the same.

PROB is calculated as the probability of obtaining a value greater or equal to CHI2_STAT under the hypothesis of a common spectrum for the two multi-channel series.

PROB must verify: size(PROB) = size(PSMAT1,1) .

NREP (INPUT, OPTIONAL) integer(i4b)

On entry, when argument NREP is present, NREP specifies the number of random interchanges of periodogram ordinates at each frequency in order to approximate the randomization distribution of CHI2_STAT under the null hypothesis.

The default is 99.

INITSEED (INPUT, OPTIONAL) logical(lgl)

On entry, if INITSEED=true, a call to RANDOM_SEED_() without arguments is done in the subroutine, in order to initiate a non-repeatable reset of the seed used by the STATPACK random generator.

The default is INITSEED=false.

Further Details¶

This statistical randomization test relies on the assumptions that the different spectral ordinates have the same sampling distribution and are independent of each other. This means, in particular, that the spectral ordinates corresponding to the zero and Nyquist frequencies must be excluded from the PSMAT1 and PSMAT2 matrices before calling SPCTRM_DIFF2 and that the two estimated multi-channel spectra have not been obtained by smoothing the periodograms in the frequency domain.

This randomization test could only be used to compare two periodograms or two spectral estimates computed as the the average of, say, r periodograms for each time series.

Finally, none of the spectral estimates must be zero.

For more details see:

Diggle, P.J., and Fisher, N.I., 1991:

Nonparametric comparison of cumulative periodograms, Applied Statistics, Vol. 40, No 3, pp. 423-434.

`subroutine power_spctrm ( vec, psvec, freq, fftvec, edof, bandwidth, conlwr, conupr, initfft, normpsd, nsmooth, trend, win, taperp, probtest )`¶

Purpose¶

Subroutine POWER_SPCTRM computes a Fast Fourier Transform (FFT) estimate of the power spectrum of a real time series, VEC. The real valued sequence VEC must be of even length.

The Power Spectral Density (PSD) estimates are returned in units which are the square of the data (if NORMPSD=false) or in spectral density units (if NORMPSD=true).

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the real time series for which the power spectrum must be estimated. If WIN/=2 or TREND=1, 2 or 3, VEC is used as workspace and is transformed.

Size(VEC) must be an even (positive) integer greater or equal to 4.

PSVEC (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the Power Spectral Density (PSD) estimates of VEC.

PSVEC must verify: size(PSVEC) = size(VEC)/2 + 1 .

FREQ (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the frequencies at which the spectral quantities are calculated in cycles per unit of time.

The spectral estimates are taken at frequencies (i-1)/size(VEC) for i=1,2, … , (size(VEC)/2 + 1).

FREQ must verify: size(FREQ) = size(VEC)/2 + 1 .

FFTVEC (OUTPUT, OPTIONAL) complex(stnd), dimension(:)

On exit, a complex vector of length (size(VEC)/2)+1 containing the Fast Fourier Transform of the product of the (detrended, e.g. the TREND argument) real time series VEC with the choosen window function (e.g. The WIN argument).

FFTVEC must verify: size(FFTVEC) = size(VEC)/2 + 1 .

EDOF (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the equivalent number of degrees of freedom of the power spectrum estimates.

EDOF must verify: size(EDOF) = size(VEC)/2 + 1 .

BANDWIDTH (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the bandwidth of the power spectrum estimates.

BANDWIDTH must verify: size(BANDWIDTH) = size(VEC)/2 + 1 .

CONLWR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

CONUPR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On output, these arguments specify the lower and upper (1-PROBTEST) * 100% confidence limit factors, respectively. Multiply the PSD estimates (e.g. the PSVEC(:) argument) by these constants to get the lower and upper limits of a (1-PROBTEST) * 100% confidence interval for the PSD estimates.

CONLWR must verify: size(CONLWR) = size(VEC)/2 + 1 .

CONUPR must verify: size(CONUPR) = size(VEC)/2 + 1 .

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if:

INITFFT = false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine POWER_SPCTRM in order to sets up constants and functions for use by subroutine FFT which is called inside subroutine POWER_SPCTRM. This call to INITFFT must have the following form:

call init_fft( size(VEC)/2 )

INITFFT = true, the call to INIT_FFT is done inside subroutine POWER_SPCTRM and a call to END_FFT is also done before leaving subroutine POWER_SPCTRM.

The default is INITFFT=true .

NORMPSD (INPUT, OPTIONAL) logical(lgl)

On entry, if:

NORMPSD = true, the PSD estimates are normalized in such a way that the total area under the power spectrum is equal to the variance of the time series VEC.

NORMPSD = false, the sum of the PSD estimates (e.g. sum( PSVEC(2:) ) is equal to the variance of the time series.

The default is NORMPSD=true .

NSMOOTH (INPUT, OPTIONAL) integer(i4b)

if NSMOOTH is used, the PSD estimates are computed by smoothing the periodogram with Daniell weights (e.g. a simple moving average).

On entry, NSMOOTH gives the length of the Daniell filter to be applied.

Setting NSMOOTH=0 on entry is equivalent to omit the optional argument NSMOOTH. Otherwise, NSMOOTH must be odd, greater than 2 and less or equal to size(VEC)/2+1 .

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+1 The mean of the time series is removed before computing the spectrum

TREND=+2 The drift from the time series is removed before computing the spectrum by using the formula: drift = (VEC(size(VEC)) - VEC(1))/(size(VEC) - 1)

TREND=+3 The least-squares line from the time series is removed before computing the spectrum.

For other values of TREND nothing is done before estimating the power spectrum.

The default is TREND=1, e.g. the mean of the time series is removed before the computations.

WIN (INPUT, OPTIONAL) integer(i4b)

On entry, this argument specify the data window used in the computations of the power spectrum. If:

WIN=+1 The Bartlett window is used

WIN=+2 The square window is used

WIN=+3 The Welch window is used

WIN=+4 The Hann window is used

WIN=+5 The Hamming window is used

WIN=+6 A split-cosine-bell window is used

The default is WIN=3, e.g. the Welch window is used.

TAPERP (INPUT, OPTIONAL) real(stnd)

The total percentage of the data to be tapered if WIN=6. TAPERP must be greater than zero and less or equal to one, otherwise the default value is used.

The default is 0.2 .

PROBTEST (INPUT, OPTIONAL) real(stnd)

On entry, a probability. PROBTEST is the critical probability which is used to determine the lower and upper confidence limit factors (e.g. the optional arguments CONLWR and CONUPR ).

PROBTEST must verify: 0. < P < 1.

The default is 0.05 .

Further Details¶

After removing the mean or the trend from the time series (e.g. TREND=1,2,3), the selected data window (e.g. WIN=1,2,3,4,5,6) is applied to the time series and the PSD estimates are computed by the FFT of this transformed time series. Optionally, theses PSD estimates may then be smoothed in the frequency domain by a Daniell filter (e.g. if NSMOOTH is used).

For definitions, more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York.

Welch, P.D., 1967:

The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE trans. on audio and electroacoustics, Vol. Au-15, 2, 70-73.

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford.

`subroutine power_spctrm ( mat, psmat, freq, fftmat, edof, bandwidth, conlwr, conupr, initfft, normpsd, nsmooth, trend, win, taperp, probtest )`¶

Purpose¶

Subroutine POWER_SPCTRM computes a Fast Fourier Transform (FFT) estimate of the power spectra of the rows of the real matrix, MAT. size(MAT,2) must be of even length.

The Power Spectral Density (PSD) estimates are returned in units which are the square of the data (if NORMPSD=false) or in spectral density units (if NORMPSD=true).

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the real time series for which power spectra must be estimated. Each row of MAT is a real time series.

If WIN/=2 or TREND=1, 2 or 3, MAT is used as workspace and is transformed.

Size(MAT,2) must be an even (positive) integer greater or equal to 4.

PSMAT (OUTPUT) real(stnd), dimension(:,:)

On exit, a real matrix containing the Power Spectral Density (PSD) estimates for each row of the real matrix MAT.

The shape of PSMAT must verify:

size(PSMAT,1) = size(MAT,1) ;

size(PSMAT,2) = size(MAT,2)/2 + 1 .

FREQ (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length (size(MAT,2)/2)+1 containing the frequencies at which the spectral quantities are calculated in cycles per unit of time.

The spectral estimates are taken at frequencies (i-1)/size(VEC) for i=1,2, … , (size(MAT,2)/2 + 1).

FREQ must verify: size(FREQ) = size(MAT,2)/2 + 1 .

FFTMAT (OUTPUT, OPTIONAL) complex(stnd), dimension(:,:)

On exit, a complex matrix containing the Fast Fourier Transform of the product of the (detrended, e.g. the TREND argument) real time series in each row of MAT with the choosen window function (e.g. The WIN argument).

The shape of FFTMAT must verify:

size(FFTMAT,1) = size(MAT,1) ;

size(FFTMAT,2) = size(MAT,2)/2 + 1 .

EDOF (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the equivalent number of degrees of freedom of the power spectrum estimates.

EDOF must verify: size(EDOF) = size(MAT,2)/2 + 1 .

BANDWIDTH (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the bandwidth of the power spectrum estimates.

BANDWIDTH must verify: size(BANDWIDTH) = size(MAT,2)/2 + 1 .

CONLWR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

CONUPR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On output, these arguments specify the lower and upper (1-PROBTEST) * 100% confidence limit factors, respectively. Multiply the PSD estimates (e.g. the PSMAT(:,:) argument) by these constants to get the lower and upper limits of a (1-PROBTEST) * 100% confidence interval for the PSD estimates.

CONLWR must verify: size(CONLWR) = size(MAT,2)/2 + 1 .

CONUPR must verify: size(CONUPR) = size(MAT,2)/2 + 1 .

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if:

INITFFT = false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine POWER_SPCTRM in order to sets up constants and functions for use by subroutine FFT which is called inside subroutine POWER_SPCTRM. This call to INITFFT must have the following form:

call init_fft( (/ size(MAT,1), size(MAT,2)/2 /), dim=2_i4b )

INITFFT = true, the call to INIT_FFT is done inside subroutine POWER_SPCTRM and a call to END_FFT is also done before leaving subroutine POWER_SPCTRM

The default is INITFFT=true .

NORMPSD (INPUT, OPTIONAL) logical(lgl)

On entry, if:

NORMPSD is set to true, the PSD estimates are normalized in such a way that the total area under the power spectrum is equal to the variance of the time series MAT.

NORMPSD = false, the sum of the PSD estimates for each row of MAT (e.g. sum( PSMAT(:,2:), dim=2 ) is equal to the variance of the corresponding time series.

The default is NORMPSD=true .

NSMOOTH (INPUT, OPTIONAL) integer(i4b)

if NSMOOTH is used, the PSD estimates are computed by smoothing the periodogram with Daniell weights (e.g. a simple moving average).

On entry, NSMOOTH gives the length of the Daniell filter to be applied.

Setting NSMOOTH=0 on entry is equivalent to omit the optional argument NSMOOTH. Otherwise, NSMOOTH must be odd, greater than 2 and less or equal to size(MAT,2)/2+1 .

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+1 The means of the time series are removed before computing the spectra

TREND=+2 The drifts from the time series are removed before computing the spectra by using the formula: drift(i) = (MAT(i,size(MAT,2)) - MAT(i,1))/(size(MAT,2) - 1)

TREND=+3 The least-squares lines from the time series are removed before computing the spectra.

For other values of TREND nothing is done before estimating the power spectra.

The default is TREND=1, e.g. the means of the time series are removed before the computations.

WIN (INPUT, OPTIONAL) integer(i4b)

On entry, this argument specify the data window used in the computations of the power spectrum. If:

WIN=+1 The Bartlett window is used

WIN=+2 The square window is used

WIN=+3 The Welch window is used

WIN=+4 The Hann window is used

WIN=+5 The Hamming window is used

WIN=+6 A split-cosine-bell window is used WIN=+1 The Bartlett window is used

The default is WIN=3, e.g. the Welch window is used.

TAPERP (INPUT, OPTIONAL) real(stnd)

The total percentage of the data to be tapered if WIN=6. TAPERP must be greater than zero and less or equal to one, otherwise the default value is used.

The default is 0.2 .

PROBTEST (INPUT, OPTIONAL) real(stnd)

On entry, a probability. PROBTEST is the critical probability which is used to determine the lower and upper confidence limit factors (e.g. the optional arguments CONLWR and CONUPR ).

PROBTEST must verify: 0. < P < 1.

The default is 0.05 .

Further Details¶

After removing the mean or the trend from the time series (e.g. TREND=1,2,3), the selected data window (e.g. WIN=1,2,3,4,5,6) is applied to the time series and the PSD estimates are computed by the FFT of these transformed time series. Optionally, theses PSD estimates may then be smoothed in the frequency domain by modified Daniell filters (e.g. if SMOOTH_PARAM is used).

For definitions, more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York.

Welch, P.D., 1967:

The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE trans. on audio and electroacoustics, Vol. Au-15, 2, 70-73.

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford.

`subroutine cross_spctrm ( vec, vec2, psvec, psvec2, phase, coher, freq, edof, bandwidth, conlwr, conupr, testcoher, ampli, co_spect, quad_spect, prob_coher, initfft, normpsd, nsmooth, trend, win, taperp, probtest )`¶

Purpose¶

Subroutine CROSS_SPCTRM computes Fast Fourier Transform (FFT) estimates of the power and cross spectra of two real time series, VEC and VEC2. The real valued sequences VEC and VEC2 must be of even length.

The Power Spectral Density (PSD) and Cross Spectral Density (CSD) estimates are returned in units which are the square of the data (if NORMPSD=false) or in spectral density units (if NORMPSD=true).

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the first real time series for which the power and cross spectra must be estimated. If WIN/=2 or TREND=1, 2 or 3, VEC is used as workspace and is transformed.

Size(VEC) must be an even (positive) integer greater or equal to 4.

VEC2 (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the second real time series for which the power and cross spectra must be estimated. If WIN/=2 or TREND=1, 2 or 3, VEC2 is used as workspace and is transformed.

VEC2 must verify: size(VEC2) = size(VEC).

PSVEC (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the Power Spectral Density (PSD) estimates of VEC.

PSVEC must verify: size(PSVEC) = size(VEC)/2 + 1 .

PSVEC2 (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC2)/2)+1 containing the Power Spectral Density (PSD) estimates of VEC2.

PSVEC2 must verify: size(PSVEC2) = size(VEC)/2 + 1 .

PHASE (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the phase of the cross spectrum, given in fractions of a circle (e.g. on the closed interval (0,1) ).

PHASE must verify: size(PHASE) = size(VEC)/2 + 1 .

COHER (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the squared coherency estimates for all frequencies.

COHER must verify: size(COHER) = size(VEC)/2 + 1 .

FREQ (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the frequencies at which the spectral quantities are calculated in cycles per unit of time.

The spectral estimates are taken at frequencies (i-1)/size(VEC) for i=1,2, … , (size(VEC)/2 + 1).

FREQ must verify: size(FREQ) = size(VEC)/2 + 1 .

EDOF (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the equivalent number of degrees of freedom of the power spectrum estimates.

EDOF must verify: size(EDOF) = size(VEC)/2 + 1 .

BANDWIDTH (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the bandwidth of the power spectrum estimates.

BANDWIDTH must verify: size(BANDWIDTH) = size(VEC)/2 + 1 .

CONLWR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

CONUPR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On output, these arguments specify the lower and upper (1-PROBTEST) * 100% confidence limit factors, respectively. Multiply the PSD estimates (e.g. the PSVEC(:) and PSVEC2(:) arguments) by these constants to get the lower and upper limits of a (1-PROBTEST) * 100% confidence interval for the PSD estimates.

CONLWR must verify: size(CONLWR) = size(VEC)/2 + 1 .

CONUPR must verify: size(CONUPR) = size(VEC)/2 + 1 .

TESTCOHER (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On output, this argument specifies the critical value for testing the null hypothesis that the squared coherency is zero at the PROBTEST * 100% significance level (e.g. elements of COHER(:) less than TESTCOHER(:) should be regarded as not significantly different from zero at the PROBTEST * 100% significance level).

TESTCOHER must verify: size(TESTCOHER) = size(VEC)/2 + 1 .

AMPLI (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the cross-amplitude spectrum.

AMPLI must verify: size(AMPLI) = (size(VEC)/2) + 1 .

CO_SPECT (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the co-spectrum (e.g. the real part of cross-spectrum).

CO_SPECT must verify: size(CO_SPECT) = (size(VEC)/2) + 1 .

QUAD_SPECT (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the quadrature spectrum (e.g. the imaginary part of cross-spectrum with a minus sign).

QUAD_SPECT must verify: size(QUAD_SPECT) = (size(VEC)/2) + 1 .

PROB_COHER (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the probabilities that the computed sample squared coherencies came from an ergodic stationary bivariate process with (corresponding) squared coherencies equal to zero.

PROB_COHER must verify: size(PROB_COHER) = (size(VEC)/2) + 1 .

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if:

INITFFT = false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine CROSS_SPCTRM in order to sets up constants and functions for use by subroutine FFT which is called inside subroutine CROSS_SPCTRM. This call to INITFFT must have the following form:

call init_fft( size(VEC)/2 )

INITFFT is set to true, the call to INIT_FFT is done inside subroutine CROSS_SPCTRM and a call to END_FFT is also done before leaving subroutine CROSS_SPCTRM.

The default is INITFFT=true .

NORMPSD (INPUT, OPTIONAL) logical(lgl)

On entry, if:

NORMPSD = true, the power and cross spectra estimates are normalized in such a way that the total area under the power spectrum is equal to the variance of the time series VEC and VEC2.

NORMPSD = false, the sum of the PSD estimates (e.g. sum(PSVEC(2:)) and sum(PSVEC2(2:)) ) is equal to the variance of the corresponding time series.

The default is NORMPSD=true .

NSMOOTH (INPUT, OPTIONAL) integer(i4b)

If NSMOOTH is used, the PSD and CSD estimates are computed by smoothing the periodogram with Daniell weights (e.g. a simple moving average).

On entry, NSMOOTH gives the length of the Daniell filter to be applied.

Setting NSMOOTH=0 on entry is equivalent to omit the optional argument NSMOOTH. Otherwise, NSMOOTH must be odd, greater than 2 and less or equal to size(VEC)/2+1 .

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+1 The mean of the two time series is removed before computing the power and cross spectra.

TREND=+2 The drift from the two time series is removed before computing the power and cross spectra.

TREND=+3 The least-squares line from the two time series is removed before computing the power and cross spectra.

For other values of TREND nothing is done before estimating the power and cross spectra.

The default is TREND=1, e.g. the means of the time series are removed before the computations.

WIN (INPUT, OPTIONAL) integer(i4b)

On entry, this argument specify the data window used in the computations of the power and cross spectra. If:

WIN=+1 The Bartlett window is used

WIN=+2 The square window is used

WIN=+3 The Welch window is used

WIN=+4 The Hann window is used

WIN=+5 The Hamming window is used

WIN=+6 A split-cosine-bell window is used

The default is WIN=3, e.g. the Welch window is used.

TAPERP (INPUT, OPTIONAL) real(stnd)

The total percentage of the data to be tapered if WIN=6. TAPERP must be greater than zero and less or equal to one, otherwise the default value is used.

The default is 0.2 .

PROBTEST (INPUT, OPTIONAL) real(stnd)

On entry, a probability. PROBTEST is the critical probability which is used to determine the lower and upper confidence limit factors (e.g. the optional arguments CONLWR and CONUPR ) and the critical value for testing the null hypothesis that the squared coherency is zero (e.g. the TESTCOHER optional argument).

PROBTEST must verify: 0. < P < 1.

The default is 0.05 .

Further Details¶

After removing the mean or the trend from the time series (e.g. TREND=1,2,3), the selected data window (e.g. WIN=1,2,3,4,5,6) is applied to the time series and the PSD and CSD estimates are computed by the FFT of these transformed time series. Optionally, theses PSD and CSD estimates may then be smoothed in the frequency domain by modified Daniell filters (e.g. if argument NSMOOTH is used).

For definitions, more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York.

Welch, P.D., 1967:

The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE trans. on audio and electroacoustics, Vol. Au-15, 2, 70-73.

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford.

`subroutine cross_spctrm ( vec, mat, psvec, psmat, phase, coher, freq, edof, bandwidth, conlwr, conupr, testcoher, ampli, co_spect, quad_spect, prob_coher, initfft, normpsd, nsmooth, trend, win, taperp, probtest )`¶

Purpose¶

Subroutine CROSS_SPCTRM computes Fast Fourier Transform (FFT) estimates of the power and cross spectra of the real time series, VEC, and the multi-channel real time series MAT.

The Power Spectral Density (PSD) and Cross Spectral Density (CSD) estimates are returned in units which are the square of the data (if NORMPSD=false) or in spectral density units (if NORMPSD=true).

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the real time series for which the power and cross spectra must be estimated. If WIN/=2 or TREND=1, 2 or 3, VEC is used as workspace and is transformed.

Size(VEC) must be an even (positive) integer greater or equal to 4.

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the multi-channel real time series for which the power and cross spectra must be estimated. Each row of MAT is a real time series. If WIN/=2 or TREND=1, 2 or 3, MAR is used as workspace and is transformed.

The shape of MAT must verify: size(MAT,2) = size(VEC).

PSVEC (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the Power Spectral Density (PSD) estimates of VEC.

PSVEC must verify: size(PSVEC) = size(VEC)/2 + 1 .

PSMAT (OUTPUT) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and (size(VEC)/2)+1 columns containing the Power Spectral Density (PSD) estimates of each row of MAT.

The shape of PSMAT must verify:

size(PSMAT,1) = size(MAT,1) ;

size(PSMAT,2) = size(VEC)/2 + 1 .

PHASE (OUTPUT) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and (size(VEC)/2)+1 columns containing the phase of the cross spectrum, given in fractions of a circle (e.g. on the closed interval (0,1) ).

The shape of PHASE must verify:

size(PHASE,1) = size(MAT,1) ;

size(PHASE,2) = size(VEC)/2 + 1 .

COHER (OUTPUT) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and (size(VEC)/2)+1 columns containing the squared coherency estimates for all frequencies.

The shape of COHER must verify:

size(COHER,1) = size(MAT,1) ;

size(COHER,2) = size(VEC)/2 + 1 .

FREQ (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the frequencies at which the spectral quantities are calculated in cycles per unit of time.

The spectral estimates are taken at frequencies (i-1)/size(VEC) for i=1,2, … , (size(VEC)/2 + 1).

FREQ must verify: size(FREQ) = size(VEC)/2 + 1 .

EDOF (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the equivalent number of degrees of freedom of the power spectrum estimates.

EDOF must verify: size(EDOF) = size(VEC)/2 + 1 .

BANDWIDTH (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the bandwidth of the power spectrum estimates.

BANDWIDTH must verify: size(BANDWIDTH) = size(VEC)/2 + 1 .

CONLWR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

CONUPR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On output, these arguments specify the lower and upper (1-PROBTEST) * 100% confidence limit factors, respectively. Multiply the PSD estimates (e.g. the PSVEC(:) and PSMAT(:,:) arguments) by these constants to get the lower and upper limits of a (1-PROBTEST) * 100% confidence interval for the PSD estimates.

CONLWR must verify: size(CONLWR) = size(VEC)/2 + 1 .

CONUPR must verify: size(CONUPR) = size(VEC)/2 + 1 .

TESTCOHER (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On output, this argument specifies the critical value for testing the null hypothesis that the squared coherency is zero at the PROBTEST * 100% significance level (e.g. elements of COHER(:,:) less than TESTCOHER(:) should be regarded as not significantly different from zero at the PROBTEST * 100% significance level).

TESTCOHER must verify: size(TESTCOHER) = size(VEC)/2 + 1 .

AMPLI (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and (size(VEC)/2)+1 columns containing the cross-amplitude spectra.

The shape of AMPLI must verify:

size(AMPLI,1) = size(MAT,1) ;

size(AMPLI,2) = (size(VEC)/2) + 1 .

CO_SPECT (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and (size(VEC)/2)+1 columns containing the co-spectra (e.g. the real part of cross-spectra).

The shape of CO_SPECT must verify:

size(CO_SPECT,1) = size(MAT,1) ;

size(CO_SPECT,2) = (size(VEC)/2) + 1 .

QUAD_SPECT (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and (size(VEC)/2)+1 columns containing the quadrature spectrum (e.g. the imaginary part of cross-spectrum with a minus sign).

The shape of QUAD_SPECT must verify:

size(QUAD_SPECT,1) = size(MAT,1) ;

size(QUAD_SPECT,2) = (size(VEC)/2) + 1 .

PROB_COHER (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and (size(VEC)/2)+1 columns containing the probabilities that the computed sample squared coherencies came from an ergodic stationary bivariate process with (corresponding) squared coherencies equal to zero.

The shape of PROB_COHER must verify:

size(PROB_COHER,1) = size(MAT,1) ;

size(PROB_COHER,2) = (size(VEC)/2) + 1 .

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if:

INITFFT = false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine CROSS_SPCTRM in order to sets up constants and functions for use by subroutine FFT which is called inside subroutine CROSS_SPCTRM. This call to INITFFT must have the following form:

call init_fft( (/ size(MAT,1), size(MAT,2)/2 /), dim=2_i4b )

INITFFT = true, the call to INIT_FFT is done inside subroutine CROSS_SPCTRM and a call to END_FFT is also done before leaving subroutine CROSS_SPCTRM.

The default is INITFFT=true .

NORMPSD (INPUT, OPTIONAL) logical(lgl)

On entry, if:

NORMPSD = true, the power and cross spectra estimates are normalized in such a way that the total area under the power spectra is equal to the variance of the time series contained in VEC and in each row of MAT.

NORMPSD = false, the sum of the PSD estimates (e.g. sum(PSVEC(2:)) and sum(PSMAT(:,2:),dim=2) ) is equal to the variance of the corresponding time series.

The default is NORMPSD=true .

NSMOOTH (INPUT, OPTIONAL) integer(i4b)

If NSMOOTH is used, the PSD and CSD estimates are computed by smoothing the periodogram with Daniell weights (e.g. a simple moving average).

On entry, NSMOOTH gives the length of the Daniell filter to be applied.

Setting NSMOOTH=0 on entry is equivalent to omit the optional argument NSMOOTH. Otherwise, NSMOOTH must be odd, greater than 2 and less or equal to size(VEC)/2+1 .

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+1 The means of the time series are removed before computing the power and cross spectra

TREND=+2 The drifts from time series are removed before computing the power and cross spectra

TREND=+3 The least-squares lines from time series are removed before computing the power and cross spectra.

For other values of TREND nothing is done before estimating the power and cross spectra.

The default is TREND=1, e.g. the means of the time series are removed before the computations.

WIN (INPUT, OPTIONAL) integer(i4b)

On entry, this argument specify the data window used in the computations of the power and cross spectra. If:

WIN=+1 The Bartlett window is used

WIN=+2 The square window is used

WIN=+3 The Welch window is used

WIN=+4 The Hann window is used

WIN=+5 The Hamming window is used

WIN=+6 A split-cosine-bell window is used

The default is WIN=3, e.g. the Welch window is used.

TAPERP (INPUT, OPTIONAL) real(stnd)

The total percentage of the data to be tapered if WIN=6. TAPERP must be greater than zero and less or equal to one, otherwise the default value is used.

The default is 0.2 .

PROBTEST (INPUT, OPTIONAL) real(stnd)

On entry, a probability. PROBTEST is the critical probability which is used to determine the lower and upper confidence limit factors (e.g. the optional arguments CONLWR and CONUPR ) and the critical value for testing the null hypothesis that the squared coherency is zero (e.g. the TESTCOHER optional argument).

PROBTEST must verify: 0. < P < 1.

The default is 0.05 .

Further Details¶

After removing the mean or the trend from the time series (e.g. TREND=1,2,3), the selected data window (e.g. WIN=1,2,3,4,5,6) is applied to the time series and the PSD and CSD estimates are computed by the FFT of these transformed time series. Optionally, theses PSD and CSD estimates may then be smoothed in the frequency domain by modified Daniell filters (e.g. if argument NSMOOTH is used).

For definitions, more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York.

Welch, P.D., 1967:

The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE trans. on audio and electroacoustics, Vol. Au-15, 2, 70-73.

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford.

`subroutine power_spctrm2 ( vec, l, psvec, freq, edof, bandwidth, conlwr, conupr, initfft, overlap, normpsd, nsmooth, trend, trend2, win, taperp, l0, probtest )`¶

Purpose¶

Subroutine POWER_SPCTRM2 computes a Fast Fourier Transform (FFT) estimate of the power spectrum of a real time series.

The Power Spectral Density (PSD) estimates are returned in units which are the square of the data (if NORMPSD=false) or in spectral density units (if NORMPSD=true).

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the real time series for which the power spectrum must be estimated. If TREND=1, 2 or 3, VEC is used as workspace and is transformed.

Size(VEC) must be greater or equal to 4.

L (INPUT) integer(i4b)

On entry, an integer used to segment the time series. L is the length of the segments. L must be a positive even integer, less or equal to size(VEC), but greater or equal to 4.

Spectral computations are at (L/2)+1 frequencies if the optional argument L0 is absent and are at ((L+L0)/2)+1 frequencies if L0 is present (L0 is the number of zeros added to each segment).

Suggested values for L+L0 are 16, 32, 64 or 128 (e.g. an integer power of two, in order to speed the computations).

PSVEC (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the Power Spectral Density (PSD) estimates of VEC.

PSVEC must verify: size(PSVEC) = ((L+L0)/2) + 1 .

FREQ (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the frequencies at which the spectral quantities are calculated in cycles per unit of time.

The spectral estimates are taken at frequencies (i-1)/(L+L0) for i=1,2, … , ((L+L0)/2 + 1).

FREQ must verify: size(FREQ) = (L+L0)/2 + 1 .

EDOF (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the equivalent number of degrees of freedom of the power spectrum estimates.

EDOF must verify: size(EDOF) = ((L+L0)/2) + 1 .

BANDWIDTH (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the bandwidth of the power spectrum estimates.

BANDWIDTH must verify: size(BANDWIDTH) = ((L+L0)/2) + 1 .

CONLWR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

CONUPR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On output, these arguments specify the lower and upper (1-PROBTEST) * 100% confidence limit factors, respectively. Multiply the PSD estimates (e.g. the PSVEC(:) argument) by these constants to get the lower and upper limits of a (1-PROBTEST) * 100% confidence interval for the PSD estimates.

CONLWR must verify: size(CONLWR) = ((L+L0)/2) + 1 .

CONUPR must verify: size(CONUPR) = ((L+L0)/2) + 1 .

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if:

INITFFT = false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine POWER_SPCTRM2 in order to sets up constants and functions for use by subroutine FFT which is called inside subroutine POWER_SPCTRM2. This call to INITFFT must have the following form:

call init_fft( (L+L0)/2 )

INITFFT = true, the call to INIT_FFT is done inside subroutine POWER_SPCTRM2 and a call to END_FFT is also done before leaving subroutine POWER_SPCTRM2.

The default is INITFFT=true .

OVERLAP (INPUT, OPTIONAL) logical(lgl)

If:

OVERLAP = false, the subroutine segments the data without any overlapping.

OVERLAP=true, the subroutine overlaps the segments by one half of their length (which is equal to L).

In both cases, zeros are eventually added to each segment (if argument L0 is present) and each segment will be FFT’d, and the resulting periodograms will averaged together to obtain a Power Spectrum Density estimate at the ((L+L0)/2)+1 frequencies.

The default is OVERLAP=false .

NORMPSD (INPUT, OPTIONAL) logical(lgl)

On entry, if:

NORMPSD = true, the PSD estimates are normalized in such a way that the total area under the power spectrum is equal to the variance of the time series VEC.

NORMPSD is set to false, the sum of the PSD estimates (e.g. sum( PSVEC(2:) ) is equal to the variance of the time series.

The default is NORMPSD=true .

NSMOOTH (INPUT, OPTIONAL) integer(i4b)

If NSMOOTH is used, the PSD estimates are computed by smoothing the periodogram with Daniell weights (e.g. a simple moving average).

On entry, NSMOOTH gives the length of the Daniell filter to be applied.

Setting NSMOOTH=0 on entry is equivalent to omit the optional argument NSMOOTH. Otherwise, NSMOOTH must be odd, greater than 2 and less or equal to (L+L0)/2+1 .

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+1 The mean of the time series is removed before computing the spectrum

TREND=+2 The drift from the time series is removed before computing the spectrum by using the formula: drift = (VEC(size(VEC)) - VEC(1))/(size(VEC) - 1)

TREND=+3 The least-squares line from the time series is removed before computing the spectrum.

For other values of TREND nothing is done before estimating the power spectrum.

The default is TREND=1, e.g. the mean of the time series is removed before the computations.

TREND2 (INPUT, OPTIONAL) integer(i4b)

If:

TREND2=+1 The mean of the time segment is removed before computing the spectrum on this segment.

TREND2=+2 The drift from the time segment is removed before computing the spectrum on this segment.

TREND2=+3 The least-squares line from the time segment is removed before computing the spectrum on this segment.

For other values of TREND2 nothing is done before estimating the power spectrum on each segment.

The default is TREND2=0, e.g. nothing is done before estimating the power spectrum on each segment.

WIN (INPUT, OPTIONAL) integer(i4b)

On entry, this argument specify the data window used in the computations of the power spectrum. If:

WIN=+1 The Bartlett window is used

WIN=+2 The square window is used

WIN=+3 The Welch window is used

WIN=+4 The Hann window is used

WIN=+5 The Hamming window is used

WIN=+6 A split-cosine-bell window is used

The default is WIN=3, e.g. the Welch window is used.

TAPERP (INPUT, OPTIONAL) real(stnd)

The total percentage of the data to be tapered if WIN=6. TAPERP must be greater than zero and less or equal to one, otherwise the default value is used.

The default is 0.2 .

L0 (INPUT, OPTIONAL) integer(i4b)

The number of zeros added to each time segment in order to obtain more finely spaced spectral estimates. L+L0 must be a positive even integer.

The default is L0=0, e.g. no zeros are added to each time segment.

PROBTEST (INPUT, OPTIONAL) real(stnd)

On entry, a probability. PROBTEST is the critical probability which is used to determine the lower and upper confidence limit factors (e.g. the optional arguments CONLWR and CONUPR ).

PROBTEST must verify: 0. < P < 1.

The default is 0.05 .

Further Details¶

After removing the mean or the trend from the time series (e.g. TREND=1,2,3), the series is padded with zero on the right such that the length of the resulting time series is evenly divisible by L (a positive even integer). The length, N, of this resulting time series is the first integer greater than or equal to size(VEC) which is evenly divisible by L. If size(VEC) is not evenly divisible by L, N is equal to size(VEC)+L-mod(size(VEC),L).

Optionally, the mean or the trend may also be removed from each time segment (e.g. TREND2=1,2,3). Optionally, zeros may be added to each time segment (e.g. the optional arguemnt L0) if more finely spaced spectral esimates are desired.

The stability of the PSD estimates depends on the averaging process. That is, the greater the number of segments ( N/L if OVERLAP=false and (2N/L)-1 if OVERLAP=true), the more stable the resulting PSD estimates.

Optionally, theses PSD estimates may then be smoothed again in the frequency domain by a Daniell filter (e.g. if argument NSMOOTH is used).

For definitions, more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York.

Welch, P.D., 1967:

The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE trans. on audio and electroacoustics, Vol. Au-15, 2, 70-73.

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford.

`subroutine power_spctrm2 ( mat, l, psmat, freq, edof, bandwidth, conlwr, conupr, initfft, overlap, normpsd, nsmooth, trend, trend2, win, taperp, l0, probtest )`¶

Purpose¶

Subroutine POWER_SPCTRM2 computes Fast Fourier Transform (FFT) estimates of the power spectra of the multi-channel real time series MAT (e.g. each row of MAT contains a time series).

The Power Spectral Density (PSD) estimates are returned in units which are the square of the data (if NORMPSD=false) or in spectral density units (if NORMPSD=true).

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the multi-channel real time series for which the power spectra must be estimated. Each row of MAT is a real time series. If TREND=1, 2 or 3, MAT is used as workspace and is transformed.

Size(MAT,2) must be greater or equal to 4.

L (INPUT) integer(i4b)

On entry, an integer used to segment the time series. L is the length of the segments. L must be a positive even integer less or equal to size(MAT,2), but greater or equal to 4.

Spectral computations are at (L/2)+1 frequencies if the optional argument L0 is absent and are at ((L+L0)/2)+1 frequencies if L0 is present (L0 is the number of zeros added to each segment).

Suggested values for L+L0 are 16, 32, 64 or 128 (e.g. an integer power of two, in order to speed the computations).

PSMAT (OUTPUT) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and ((L+L0)/2) + 1 columns containing the Power Spectral Density (PSD) estimates of each row of MAT.

The shape of PSMAT must verify:

size(PSMAT,1) = size(MAT,1) ;

size(PSMAT,2) = ((L+L0)/2) + 1 .

FREQ (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the frequencies at which the spectral quantities are calculated in cycles per unit of time.

The spectral estimates are taken at frequencies (i-1)/(L+L0) for i=1,2, … , ((L+L0)/2 + 1).

FREQ must verify: size(FREQ) = (L+L0)/2 + 1 .

EDOF (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the equivalent number of degrees of freedom of the power spectrum estimates.

EDOF must verify: size(EDOF) = ((L+L0)/2) + 1 .

BANDWIDTH (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the bandwidth of the power spectrum estimates.

BANDWIDTH must verify: size(BANDWIDTH) = ((L+L0)/2) + 1 .

CONLWR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

CONUPR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On output, these arguments specify the lower and upper (1-PROBTEST) * 100% confidence limit factors, respectively. Multiply the PSD estimates (e.g. the PSMAT(:,:) argument) by these constants to get the lower and upper limits of a (1-PROBTEST) * 100% confidence interval for the PSD estimates.

CONLWR must verify: size(CONLWR) = ((L+L0)/2) + 1 .

CONUPR must verify: size(CONUPR) = ((L+L0)/2) + 1 .

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if:

INITFFT = false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine POWER_SPCTRM2 in order to sets up constants and functions for use by subroutine FFT which is called inside subroutine POWER_SPCTRM2. This call to INITFFT must have the following form:

call init_fft( (/ size(MAT,1), (L+L0)/2 /), dim=2_i4b )

INITFFT = true, the call to INIT_FFT is done inside subroutine POWER_SPCTRM2 and a call to END_FFT is also done before leaving subroutine POWER_SPCTRM2.

The default is INITFFT=true .

OVERLAP (INPUT, OPTIONAL) logical(lgl)

If:

OVERLAP = false, the subroutine segments the data without any overlapping.

OVERLAP=true, the subroutine overlaps the segments by one half of their length (which is equal to L).

In both cases, zeros are eventually added to each segment (if argument L0 is present) and each segment will be FFT’d, and the resulting periodograms will averaged together to obtain a Power Spectrum Density estimate at the ((L+L0)/2)+1 frequencies.

The default is OVERLAP=false .

NORMPSD (INPUT, OPTIONAL) logical(lgl)

On entry, if:

NORMPSD = true, the PSD estimates are normalized in such a way that the total area under the power spectrum is equal to the variance of the time series VEC.

NORMPSD is set to false, the sum of the PSD estimates (e.g. sum(PSMAT(:,2:),dim=2) ) is equal to the variance of the corresponding time series.

The default is NORMPSD=true .

NSMOOTH (INPUT, OPTIONAL) integer(i4b)

If NSMOOTH is used, the PSD estimates are computed by smoothing the periodogram with Daniell weights (e.g. a simple moving average).

On entry, NSMOOTH gives the length of the Daniell filter to be applied.

Setting NSMOOTH=0 on entry is equivalent to omit the optional argument NSMOOTH. Otherwise, NSMOOTH must be odd, greater than 2 and less or equal to (L+L0)/2+1 .

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+1 The means of the time series are removed before computing the spectra

TREND=+2 The drifts from time series are removed before computing the spectra

TREND=+3 The least-squares lines from time series are removed before computing the spectra.

For other values of TREND nothing is done before estimating the power and cross spectra. The default is TREND=1, e.g. the means of the time series are removed before the computations.

TREND2 (INPUT, OPTIONAL) integer(i4b)

If:

TREND2=+1 The mean of the time segment is removed before computing the spectrum on this segment.

TREND2=+2 The drift from the time segment is removed before computing the spectrum on this segment.

TREND2=+3 The least-squares line from the time segment is removed before computing the spectrum on this segment.

For other values of TREND2 nothing is done before estimating the power spectrum on each segment.

The default is TREND2=0, e.g. nothing is done before estimating the power spectrum on each segment.

WIN (INPUT, OPTIONAL) integer(i4b)

On entry, this argument specify the data window used in the computations of the power spectrum. If:

WIN=+1 The Bartlett window is used

WIN=+2 The square window is used

WIN=+3 The Welch window is used

WIN=+4 The Hann window is used

WIN=+5 The Hamming window is used

WIN=+6 A split-cosine-bell window is used

The default is WIN=3, e.g. the Welch window is used.

TAPERP (INPUT, OPTIONAL) real(stnd)

The total percentage of the data to be tapered if WIN=6. TAPERP must be greater than zero and less or equal to one, otherwise the default value is used.

The default is 0.2 .

L0 (INPUT, OPTIONAL) integer(i4b)

The number of zeros added to each time segment in order to obtain more finely spaced spectral estimates. L+L0 must be a positive even integer.

The default is L0=0, e.g. no zeros are added to each time segment.

PROBTEST (INPUT, OPTIONAL) real(stnd)

On entry, a probability. PROBTEST is the critical probability which is used to determine the lower and upper confidence limit factors (e.g. the optional arguments CONLWR and CONUPR).

PROBTEST must verify: 0. < P < 1.

The default is 0.05 .

Further Details¶

After removing the mean or the trend from the time series (e.g. TREND=1,2,3), the series are padded with zero on the right such that the length of the resulting time series is evenly divisible by L (a positive even integer). The length, N, of this resulting time series is the first integer greater than or equal to size(MAT,2) which is evenly divisible by L. If size(MAT,2) is not evenly divisible by L, N is equal to size(MAT,2)+L-mod(size(MAT,2),L).

Optionally, the mean or the trend may also be removed from each time segment (e.g. TREND2=1,2,3). Optionally, zeros may be added to each time segment (e.g. the optional arguemnt L0) if more finely spaced spectral esimates are desired.

The stability of the PSD estimates depends on the averaging process. That is, the greater the number of segments ( N/L if OVERLAP=false and (2N/L)-1 if OVERLAP=true), the more stable the resulting PSD estimates.

Optionally, theses PSD estimates may then be smoothed again in the frequency domain by modified Daniell filters (e.g. if argument NSMOOTH is used).

For definitions, more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York.

Welch, P.D., 1967:

The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE trans. on audio and electroacoustics, Vol. Au-15, 2, 70-73.

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford.

`subroutine cross_spctrm2 ( vec, vec2, l, psvec, psvec2, phase, coher, freq, edof, bandwidth, conlwr, conupr, testcoher, ampli, co_spect, quad_spect, prob_coher, initfft, overlap, normpsd, nsmooth, trend, trend2, win, taperp, l0, probtest )`¶

Purpose¶

Subroutine CROSS_SPCTRM2 computes Fast Fourier Transform (FFT) estimates of the power and cross spectra of two real time series.

The Power Spectral Density (PSD) and Cross Spectral Density (CSD) estimates are returned in units which are the square of the data (if NORMPSD=false) or in spectral density units (if NORMPSD=true).

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the first real time series for which the power and cross spectra must be estimated. If TREND=1, 2 or 3, VEC is used as workspace and is transformed.

Size(VEC) must be greater or equal to 4.

VEC2 (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the second real time series for which the power and cross spectra must be estimated. If TREND=1, 2 or 3, VEC2 is used as workspace and is transformed.

VEC2 must verify: size(VEC2) = size(VEC).

L (INPUT) integer(i4b)

On entry, an integer used to segment the time series. L is the length of the segments. L must be a positive even integer, less or equal to size(VEC), but greater or equal to 4.

Spectral computations are at (L/2)+1 frequencies if the optional argument L0 is absent and are at ((L+L0)/2)+1 frequencies if L0 is present (L0 is the number of zeros added to each segment).

Suggested values for L+L0 are 16, 32, 64 or 128 (e.g. an integer power of two, in order to speed the computations).

PSVEC (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the Power Spectral Density (PSD) estimates of VEC.

PSVEC must verify: size(PSVEC) = ((L+L0)/2) + 1 .

PSVEC2 (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the Power Spectral Density (PSD) estimates of VEC2.

PSVEC2 must verify: size(PSVEC2) = ((L+L0)/2) + 1 .

PHASE (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the phase of the cross spectrum, given in fractions of a circle (e.g. on the closed interval (0,1) ).

PHASE must verify: size(PHASE) = ((L+L0)/2) + 1 .

COHER (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the squared coherency estimates for all frequencies.

COHER must verify: size(COHER) = ((L+L0)/2) + 1 .

FREQ (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the frequencies at which the spectral quantities are calculated in cycles per unit of time.

The spectral estimates are taken at frequencies (i-1)/(L+L0) for i=1,2, … , ((L+L0)/2 + 1).

FREQ must verify: size(FREQ) = (L+L0)/2 + 1 .

EDOF (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the equivalent number of degrees of freedom of the power spectrum estimates.

EDOF must verify: size(EDOF) = ((L+L0)/2) + 1 .

BANDWIDTH (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the bandwidth of the power spectrum estimates.

BANDWIDTH must verify: size(BANDWIDTH) = ((L+L0)/2) + 1 .

CONLWR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

CONUPR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On output, these arguments specify the lower and upper (1-PROBTEST) * 100% confidence limit factors, respectively. Multiply the PSD estimates (e.g. the PSVEC(:) and PSVEC2(:) arguments) by these constants to get the lower and upper limits of a (1-PROBTEST) * 100% confidence interval for the PSD estimates.

CONLWR must verify: size(CONLWR) = ((L+L0)/2) + 1 .

CONUPR must verify: size(CONUPR) = ((L+L0)/2) + 1 .

TESTCOHER (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On output, this argument specifies the critical value for testing the null hypothesis that the squared coherency is zero at the PROBTEST * 100% significance level (e.g. elements of COHER(:) less than TESTCOHER(:) should be regarded as not significantly different from zero at the PROBTEST * 100% significance level).

TESTCOHER must verify: size(TESTCOHER) = ((L+L0)/2) + 1 .

AMPLI (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the cross-amplitude spectrum.

AMPLI must verify: size(AMPLI) = ((L+L0)/2) + 1 .

CO_SPECT (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the co-spectrum (e.g. the real part of cross-spectrum).

CO_SPECT must verify: size(CO_SPECT) = ((L+L0)/2) + 1 .

QUAD_SPECT (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the quadrature spectrum (e.g. the imaginary part of cross-spectrum with a minus sign).

QUAD_SPECT must verify: size(QUAD_SPECT) = ((L+L0)/2) + 1 .

PROB_COHER (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the probabilities that the computed sample squared coherencies came from an ergodic stationary bivariate process with (corresponding) squared coherencies equal to zero.

PROB_COHER must verify: size(PROB_COHER) = ((L+L0)/2)+1 .

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if:

INITFFT = false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine CROSS_SPCTRM2 in order to sets up constants and functions for use by subroutine FFT which is called inside subroutine CROSS_SPCTRM2. This call to INITFFT must have the following form:

call init_fft( (L+L0)/2 )

INITFFT = true, the call to INIT_FFT is done inside subroutine CROSS_SPCTRM2 and a call to END_FFT is also done before leaving subroutine CROSS_SPCTRM2.

The default is INITFFT=true .

OVERLAP (INPUT, OPTIONAL) logical(lgl)

If:

OVERLAP = false, the subroutine segments the data without any overlapping.

OVERLAP=true, the subroutine overlaps the segments by one half of their length (which is equal to L).

In both cases, zeros are eventually added to each segment (if argument L0 is present) and each segment will be FFT’d, and the resulting periodograms will averaged together to obtain a Power Spectrum Density estimate at the ((L+L0)/2)+1 frequencies.

The default is OVERLAP=false .

NORMPSD (INPUT, OPTIONAL) logical(lgl)

On entry, if:

NORMPSD = true, the power and cross spectra estimates are normalized in such a way that the total area under the power spectrum is equal to the variance of the time series VEC and VEC2.

NORMPSD = false, the sum of the PSD estimates (e.g. sum(PSVEC(2:)) and sum(PSVEC2(2:)) ) is equal to the variance of the corresponding time series.

The default is NORMPSD=true .

NSMOOTH (INPUT, OPTIONAL) integer(i4b)

if NSMOOTH is used, the PSD and CSD estimates are computed by smoothing the periodogram with Daniell weights (e.g. a simple moving average).

On entry, NSMOOTH gives the length of the Daniell filter to be applied.

Setting NSMOOTH=0 on entry is equivalent to omit the optional argument NSMOOTH. Otherwise, NSMOOTH must be odd, greater than 2 and less or equal to (L+L0)/2+1 .

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+1 The mean of the two time series is removed before computing the spectra

TREND=+2 The drift from the two time series is removed before computing the spectra

TREND=+3 The least-squares line from the two time series is removed before computing the spectra.

For other values of TREND nothing is done before estimating the power and cross spectra. The default is TREND=1, e.g. the means of the time series are removed before the computations.

TREND2 (INPUT, OPTIONAL) integer(i4b)

If:

TREND2=+1 The mean of the time segment is removed before computing the cross-spectrum on this segment.

TREND2=+2 The drift from the time segment is removed before computing the cross-spectrum on this segment.

TREND2=+3 The least-squares line from the time segment is removed before computing the cross-spectrum on this segment.

For other values of TREND2 nothing is done before estimating the cross-spectrum on each segment.

The default is TREND2=0, e.g. nothing is done before estimating the power spectrum on each segment.

WIN (INPUT, OPTIONAL) integer(i4b)

On entry, this argument specify the data window used in the computations of the power and cross spectra. If:

WIN=+1 The Bartlett window is used

WIN=+2 The square window is used

WIN=+3 The Welch window is used

WIN=+4 The Hann window is used

WIN=+5 The Hamming window is used

WIN=+6 A split-cosine-bell window is used

The default is WIN=3, e.g. the Welch window is used.

TAPERP (INPUT, OPTIONAL) real(stnd)

The total percentage of the data to be tapered if WIN=6. TAPERP must be greater than zero and less or equal to one, otherwise the default value is used.

The default is 0.2 .

L0 (INPUT, OPTIONAL) integer(i4b)

The number of zeros added to each time segment in order to obtain more finely spaced spectral estimates. L+L0 must be a positive even integer.

The default is L0=0, e.g. no zeros are added to each time segment.

PROBTEST (INPUT, OPTIONAL) real(stnd)

On entry, a probability. PROBTEST is the critical probability which is used to determine the lower and upper confidence limit factors (e.g. the optional arguments CONLWR and CONUPR ) and the critical value for testing the null hypothesis that the squared coherency is zero (e.g. the TESTCOHER optional argument).

PROBTEST must verify: 0. < P < 1.

The default is 0.05 .

Further Details¶

After removing the mean or the trend from the two time series (e.g. TREND=1,2,3), the series are padded with zero on the right such that the length of the resulting two time series is evenly divisible by L (a positive even integer). The length, N, of these resulting time series is the first integer greater than or equal to size(VEC) which is evenly divisible by L. If size(VEC) is not evenly divisible by L, N is equal to size(VEC)+L-mod(size(VEC),L).

Optionally, the mean or the trend may also be removed from each time segment (e.g. TREND2=1,2,3). Optionally, zeros may be added to each time segment (e.g. the optional arguemnt L0) if more finely spaced spectral esimates are desired.

The stability of the power and cross spectra estimates depends on the averaging process. That is, the greater the number of segments ( N/L if OVERLAP=false and (2N/L)-1 if OVERLAP=true), the more stable the resulting power and cross spectra estimates.

Optionally, these power and cross spectra estimates may then be smoothed again in the frequency domain by modified Daniell filters (e.g. if argument NSMOOTH is used).

For definitions, more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York.

Welch, P.D., 1967:

The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE trans. on audio and electroacoustics, Vol. Au-15, 2, 70-73.

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford.

`subroutine cross_spctrm2 ( vec, mat, l, psvec, psmat, phase, coher, freq, edof, bandwidth, conlwr, conupr, testcoher, ampli, co_spect, quad_spect, prob_coher, initfft, overlap, normpsd, nsmooth, trend, trend2, win, taperp, l0, probtest )`¶

Purpose¶

Subroutine CROSS_SPCTRM2 computes Fast Fourier Transform (FFT) estimates of the power and cross spectra of the real time series, VEC, and the multi-channel real time series MAT.

The Power Spectral Density (PSD) and Cross Spectral Density (CSD) estimates are returned in units which are the square of the data (if NORMPSD=false) or in spectral density units (if NORMPSD=true).

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the real time series for which the power and cross spectra must be estimated. If TREND=1, 2 or 3, VEC is used as workspace and is transformed.

Size(VEC) must be greater or equal to 4.

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the multi-channel real time series for which the power and cross spectra must be estimated. Each row of MAT is a real time series. If TREND=1, 2 or 3, MAT is used as workspace and is transformed.

The shape of MAT must verify: size(MAT,2) = size(VEC).

L (INPUT) integer(i4b)

On entry, an integer used to segment the time series. L is the length of the segments. L must be a positive even integer, less or equal to size(VEC), but greater or equal to 4.

Spectral computations are at (L/2)+1 frequencies if the optional argument L0 is absent and are at ((L+L0)/2)+1 frequencies if L0 is present (L0 is the number of zeros added to each segment).

Suggested values for L+L0 are 16, 32, 64 or 128 (e.g. an integer power of two, in order to speed the computations).

PSVEC (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the Power Spectral Density (PSD) estimates of VEC.

PSVEC must verify: size(PSVEC) = ((L+L0)/2) + 1 .

PSMAT (OUTPUT) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and ((L+L0)/2) + 1 columns containing the Power Spectral Density (PSD) estimates of each row of MAT.

The shape of PSMAT must verify:

size(PSMAT,1) = size(MAT,1) ;

size(PSMAT,2) = ((L+L0)/2) + 1 .

PHASE (OUTPUT) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and ((L+L0)/2) + 1 columns containing the phase of the cross spectrum, given in fractions of a circle (e.g. on the closed interval (0,1) ).

The shape of PHASE must verify:

size(PHASE,1) = size(MAT,1) ;

size(PHASE,2) = ((L+L0)/2) + 1 .

COHER (OUTPUT) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and ((L+L0)/2) + 1 columns containing the squared coherency estimates for all frequencies.

The shape of COHER must verify:

size(COHER,1) = size(MAT,1) ;

size(COHER,2) = ((L+L0)/2) + 1 .

FREQ (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the frequencies at which the spectral quantities are calculated in cycles per unit of time.

The spectral estimates are taken at frequencies (i-1)/(L+L0) for i=1,2, … , ((L+L0)/2 + 1).

FREQ must verify: size(FREQ) = (L+L0)/2 + 1 .

EDOF (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the equivalent number of degrees of freedom of the power spectrum estimates.

EDOF must verify: size(EDOF) = ((L+L0)/2) + 1 .

BANDWIDTH (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, the bandwidth of the power spectrum estimates.

BANDWIDTH must verify: size(BANDWIDTH) = ((L+L0)/2) + 1 .

CONLWR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

CONUPR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On output, these arguments specify the lower and upper (1-PROBTEST) * 100% confidence limit factors, respectively. Multiply the PSD estimates (e.g. the PSVEC(:) and PSMAT(:,:) arguments) by these constants to get the lower and upper limits of a (1-PROBTEST) * 100% confidence interval for the PSD estimates.

CONLWR must verify: size(CONLWR) = ((L+L0)/2) + 1 .

CONUPR must verify: size(CONUPR) = ((L+L0)/2) + 1 .

TESTCOHER (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On output, this argument specifies the critical value for testing the null hypothesis that the squared coherency is zero at the PROBTEST * 100% significance level (e.g. elements of COHER(:,:) less than TESTCOHER(:) should be regarded as not significantly different from zero at the PROBTEST * 100% significance level).

TESTCOHER must verify: size(TESTCOHER) = ((L+L0)/2) + 1 .

AMPLI (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and ((L+L0)/2) + 1 columns containing the cross-amplitude spectra.

The shape of AMPLI must verify:

size(AMPLI,1) = size(MAT,1) ;

size(AMPLI,2) = ((L+L0)/2) + 1 .

CO_SPECT (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and ((L+L0)/2) + 1 columns containing the co-spectra (e.g. the real part of cross-spectra).

The shape of CO_SPECT must verify:

size(CO_SPECT,1) = size(MAT,1) ;

size(CO_SPECT,2) = ((L+L0)/2) + 1 .

QUAD_SPECT (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and ((L+L0)/2) + 1 columns containing the quadrature spectrum (e.g. the imaginary part of cross-spectrum with a minus sign).

The shape of QUAD_SPECT must verify:

size(QUAD_SPECT,1) = size(MAT,1) ;

size(QUAD_SPECT,2) = ((L+L0)/2) + 1 .

PROB_COHER (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and ((L+L0)/2) + 1 columns containing the probabilities that the computed sample squared coherencies came from an ergodic stationary bivariate process with (corresponding) squared coherencies equal to zero.

The shape of PROB_COHER must verify:

size(PROB_COHER,1) = size(MAT,1) ; size(PROB_COHER,2) = ((L+L0)/2) + 1 .

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if:

INITFFT = false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine CROSS_SPCTRM2 in order to sets up constants and functions for use by subroutine FFT which is called inside subroutine CROSS_SPCTRM2. This call to INITFFT must have the following form:

call init_fft( (/ size(MAT,1), (L+L0)/2 /), dim=2_i4b )

INITFFT is set to true, the call to INIT_FFT is done inside subroutine CROSS_SPCTRM2 and a call to END_FFT is also done before leaving subroutine CROSS_SPCTRM2.

The default is INITFFT=true .

OVERLAP (INPUT, OPTIONAL) logical(lgl)

If:

OVERLAP = false, the subroutine segments the data without any overlapping.

OVERLAP=true, the subroutine overlaps the segments by one half of their length (which is equal to L).

In both cases, zeros are eventually added to each segment (if argument L0 is present) and each segment will be FFT’d, and the resulting periodograms will averaged together to obtain a Power Spectrum Density estimate at the ((L+L0)/2)+1 frequencies.

The default is OVERLAP=false .

NORMPSD (INPUT, OPTIONAL) logical(lgl)

On entry, if:

NORMPSD = true, the power and cross spectra estimates are normalized in such a way that the total area under the power spectra is equal to the variance of the time series contained in VEC and in each row of MAT.

NORMPSD is set to false, the sum of the PSD estimates (e.g. sum(PSVEC(2:)) and sum(PSMAT(:,2:),dim=2) ) is equal to the variance of the corresponding time series.

The default is NORMPSD=true .

NSMOOTH (INPUT, OPTIONAL) integer(i4b)

If NSMOOTH is used, the PSD estimates are computed by smoothing the periodogram with Daniell weights (e.g. a simple moving average).

On entry, NSMOOTH gives the length of the Daniell filter to be applied.

Setting NSMOOTH=0 on entry is equivalent to omit the optional argument NSMOOTH. Otherwise, NSMOOTH must be odd, greater than 2 and less or equal to (L+L0)/2+1 .

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+1 The means of the time series are removed before computing the spectra

TREND=+2 The drifts from time series are removed before computing the spectra

TREND=+3 The least-squares lines from time series are removed before computing the spectra.

For other values of TREND nothing is done before estimating the power and cross spectra. The default is TREND=1, e.g. the means of the time series are removed before the computations.

TREND2 (INPUT, OPTIONAL) integer(i4b)

If:

TREND2=+1 The mean of the time segment is removed before computing the cross-spectrum on this segment.

TREND2=+2 The drift from the time segment is removed before computing the cross-spectrum on this segment.

TREND2=+3 The least-squares line from the time segment is removed before computing the cross-spectrum on this segment.

For other values of TREND2 nothing is done before estimating the cross-spectrum on each segment.

The default is TREND2=0, e.g. nothing is done before estimating the power spectrum on each segment.

WIN (INPUT, OPTIONAL) integer(i4b)

On entry, this argument specify the data window used in the computations of the power and cross spectra. If:

WIN=+1 The Bartlett window is used

WIN=+2 The square window is used

WIN=+3 The Welch window is used

WIN=+4 The Hann window is used

WIN=+5 The Hamming window is used

WIN=+6 A split-cosine-bell window is used

The default is WIN=3, e.g. the Welch window is used.

TAPERP (INPUT, OPTIONAL) real(stnd)

The total percentage of the data to be tapered if WIN=6. TAPERP must be greater than zero and less or equal to one, otherwise the default value is used.

The default is 0.2 .

L0 (INPUT, OPTIONAL) integer(i4b)

The number of zeros added to each time segment in order to obtain more finely spaced spectral estimates. L+L0 must be a positive even integer.

The default is L0=0, e.g. no zeros are added to each time segment.

PROBTEST (INPUT, OPTIONAL) real(stnd)

On entry, a probability. PROBTEST is the critical probability which is used to determine the lower and upper confidence limit factors (e.g. the optional arguments CONLWR and CONUPR ) and the critical value for testing the null hypothesis that the squared coherency is zero (e.g. the TESTCOHER optional argument).

PROBTEST must verify 0. < P < 1.

The default is 0.05 .

Further Details¶

After removing the mean or the trend from the time series (e.g. TREND=1,2,3), the series are padded with zero on the right such that the length of the resulting time series is evenly divisible by L (a positive even integer). The length, N, of these resulting time series is the first integer greater than or equal to size(VEC) which is evenly divisible by L. If size(VEC) is not evenly divisible by L, N is equal to size(VEC)+L-mod(size(VEC),L).

Optionally, the mean or the trend may also be removed from each time segment (e.g. TREND2=1,2,3). Optionally, zeros may be added to each time segment (e.g. the optional arguemnt L0) if more finely spaced spectral esimates are desired.

The stability of the power and cross spectra estimates depends on the averaging process. That is, the greater the number of segments ( N/L if OVERLAP=false and (2N/L)-1 if OVERLAP=true), the more stable the resulting power and cross spectra estimates.

Optionally, these power and cross spectra estimates may then be smoothed again in the frequency domain by modified Daniell filters (e.g. if argument NSMOOTH is used).

For definitions, more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York.

Welch, P.D., 1967:

The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE trans. on audio and electroacoustics, Vol. Au-15, 2, 70-73.

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford.

`subroutine power_spectrum ( vec, psvec, freq, fftvec, edof, bandwidth, conlwr, conupr, initfft, normpsd, smooth_param, trend, win, taperp, probtest )`¶

Purpose¶

Subroutine POWER_SPECTRUM computes a Fast Fourier Transform (FFT) estimate of the power spectrum of a real time series, VEC. The real valued sequence VEC must be of even length.

The Power Spectral Density (PSD) estimates are returned in units which are the square of the data (if NORMPSD=false) or in spectral density units (if NORMPSD=true).

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the real time series for which the power spectrum must be estimated. If WIN/=2 or TREND=1, 2 or 3, VEC is used as workspace and is transformed.

Size(VEC) must be an even (positive) integer greater or equal to 4.

PSVEC (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the Power Spectral Density (PSD) estimates of VEC.

PSVEC must verify: size(PSVEC) = size(VEC)/2 + 1 .

FREQ (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the frequencies at which the spectral quantities are calculated in cycles per unit of time.

The spectral estimates are taken at frequencies (i-1)/size(VEC) for i=1,2, … , (size(VEC)/2 + 1).

FREQ must verify: size(FREQ) = size(VEC)/2 + 1 .

FFTVEC (OUTPUT, OPTIONAL) complex(stnd), dimension(:)

On exit, a complex vector of length (size(VEC)/2)+1 containing the Fast Fourier Transform of the product of the (detrended, e.g. the TREND argument) real time series VEC with the choosen window function (e.g. The WIN argument).

FFTVEC must verify: size(FFTVEC) = size(VEC)/2 + 1 .

EDOF (OUTPUT, OPTIONAL) real(stnd)

On exit, the equivalent number of degrees of freedom of the power spectrum estimates.

BANDWIDTH (OUTPUT, OPTIONAL) real(stnd)

On exit, the bandwidth of the power spectrum estimates.

CONLWR (OUTPUT, OPTIONAL) real(stnd)

CONUPR (OUTPUT, OPTIONAL) real(stnd)

On output, these arguments specify the lower and upper (1-PROBTEST) * 100% confidence limit factors, respectively. Multiply the PSD estimates (e.g. the PSVEC(:) argument) by these constants to get the lower and upper limits of a (1-PROBTEST) * 100% confidence interval for the PSD estimates.

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if:

INITFFT = false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine POWER_SPECTRUM in order to sets up constants and functions for use by subroutine FFT which is called inside subroutine POWER_SPECTRUM. This call to INITFFT must have the following form:

call init_fft( size(VEC)/2 )

INITFFT = true, the call to INIT_FFT is done inside subroutine POWER_SPECTRUM and a call to END_FFT is also done before leaving subroutine POWER_SPECTRUM.

The default is INITFFT=true .

NORMPSD (INPUT, OPTIONAL) logical(lgl)

On entry, if:

NORMPSD = true, the PSD estimates are normalized in such a way that the total area under the power spectrum is equal to the variance of the time series VEC.

NORMPSD is set to false, the sum of the PSD estimates (e.g. sum( PSVEC(2:) ) is equal to the variance of the time series.

The default is NORMPSD=true .

SMOOTH_PARAM (INPUT, OPTIONAL) integer(i4b), dimension(:)

if SMOOTH_PARAM is used, the PSD estimates are computed by repeated smoothing of the periodogram with modified Daniell weights.

On entry, SMOOTH_PARAM(:) gives the array of the half-lengths of the modified Daniell filters to be applied.

All the values in SMOOTH_PARAM(:) must be greater than 0 and less than size(VEC)/2+1 .

Size(SMOOTH_PARAM) must be greater or equal to 1.

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+1 The mean of the time series is removed before computing the spectrum

TREND=+2 The drift from the time series is removed before computing the spectrum by using the formula: drift = (VEC(size(VEC)) - VEC(1))/(size(VEC) - 1)

TREND=+3 The least-squares line from the time series is removed before computing the spectrum.

For other values of TREND nothing is done before estimating the power spectrum.

The default is TREND=1, e.g. the mean of the time series is removed before the computations.

WIN (INPUT, OPTIONAL) integer(i4b)

On entry, this argument specify the data window used in the computations of the power spectrum. If:

WIN=+1 The Bartlett window is used

WIN=+2 The square window is used

WIN=+3 The Welch window is used

WIN=+4 The Hann window is used

WIN=+5 The Hamming window is used

WIN=+6 A split-cosine-bell window is used

The default is WIN=3, e.g. the Welch window is used.

TAPERP (INPUT, OPTIONAL) real(stnd)

The total percentage of the data to be tapered if WIN=6. TAPERP must be greater than zero and less or equal to one, otherwise the default value is used.

The default is 0.2 .

PROBTEST (INPUT, OPTIONAL) real(stnd)

On entry, a probability. PROBTEST is the critical probability which is used to determine the lower and upper confidence limit factors (e.g. the optional arguments CONLWR and CONUPR).

PROBTEST must verify: 0. < P < 1.

The default is 0.05 .

Further Details¶

After removing the mean or the trend from the time series (e.g. TREND=1,2,3), the selected data window (e.g. WIN=1,2,3,4,5,6) is applied to the time series and the PSD estimates are computed by the FFT of this transformed time series. Optionally, theses PSD estimates may then be smoothed in the frequency domain by modified Daniell filters (e.g. if SMOOTH_PARAM is used).

The computed equivalent number of degrees of freedom and bandwidth must be divided by two for the zero and Nyquist frequencies.

Furthermore, the computed equivalent number of degrees of freedom, bandwidth, lower and upper (1-PROBTEST) * 100% confidence limit factors are not right near the zero and Nyquist frequencies if the PSD estimates have been smoothed by modified Daniell filters. The reason is that POWER_SPECTRUM assumes that smoothing involves averaging independent frequency ordinates. This is true except near the zero and Nyquist frequencies where an average may contain contributions from negative frequencies, which are identical to and hence not independent of positive frequency spectral values. Thus, the number of degrees of freedom in PSD estimates near the 0 and Nyquist frequencies are as little as half the number of degrees of freedom of the spectral estimates away from these frequency extremes if the optional argument SMOOTH_PARAM is used.

If the optional argument SMOOTH_PARAM is used, the computed equivalent number of degrees of freedom, bandwidth, lower and upper (1-PROBTEST) * 100% confidence limit factors are right for PSD estimates at frequencies

(i-1)/Size(VEC) for i= (nparam+1)/2 + 1 to ( Size(VEC) - nparam + 1)/2

where nparam = 2 * (2+sum(SMOOTH_PARAM(:)))- 1, (e.g. for frequencies i/Size(VEC) for i = (nparam+1)/2, … , ( Size(VEC) - nparam - 1)/2 ).

For definitions, more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York.

Welch, P.D., 1967:

The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE trans. on audio and electroacoustics, Vol. Au-15, 2, 70-73.

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford.

`subroutine power_spectrum ( mat, psmat, freq, fftmat, edof, bandwidth, conlwr, conupr, initfft, normpsd, smooth_param, trend, win, taperp, probtest )`¶

Purpose¶

Subroutine POWER_SPECTRUM computes a Fast Fourier Transform (FFT) estimate of the power spectra of the rows of the real matrix, MAT. size(MAT,2) must be of even length.

The Power Spectral Density (PSD) estimates are returned in units which are the square of the data (if NORMPSD=false) or in spectral density units (if NORMPSD=true).

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the real time series for which power spectra must be estimated. Each row of MAT is a real time series. If WIN/=2 or TREND=1, 2 or 3, MAT is used as workspace and is transformed.

Size(MAT,2) must be an even (positive) integer greater or equal to 4.

PSMAT (OUTPUT) real(stnd), dimension(:,:)

On exit, a real matrix containing the Power Spectral Density (PSD) estimates for each row of the real matrix MAT.

The shape of PSMAT must verify:

size(PSMAT,1) = size(MAT,1) ;

size(PSMAT,2) = size(MAT,2)/2 + 1 .

FREQ (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length (size(MAT,2)/2)+1 containing the frequencies at which the spectral quantities are calculated in cycles per unit of time.

The spectral estimates are taken at frequencies (i-1)/size(VEC) for i=1,2, … , (size(MAT,2)/2 + 1).

FREQ must verify: size(FREQ) = size(MAT,2)/2 + 1 .

FFTMAT (OUTPUT, OPTIONAL) complex(stnd), dimension(:,:)

On exit, a complex matrix containing the Fast Fourier Transform of the product of the (detrended, e.g. the TREND argument) real time series in each row of MAT with the choosen window function (e.g. The WIN argument).

The shape of FFTMAT must verify:

size(FFTMAT,1) = size(MAT,1) ;

size(FFTMAT,2) = size(MAT,2)/2 + 1 .

EDOF (OUTPUT, OPTIONAL) real(stnd)

On exit, the equivalent number of degrees of freedom of the power spectrum estimates.

BANDWIDTH (OUTPUT, OPTIONAL) real(stnd)

On exit, the bandwidth of the power spectrum estimates.

CONLWR (OUTPUT, OPTIONAL) real(stnd)

CONUPR (OUTPUT, OPTIONAL) real(stnd)

On output, these arguments specify the lower and upper (1-PROBTEST) * 100% confidence limit factors, respectively. Multiply the PSD estimates (e.g. the PSMAT(:,:) argument) by these constants to get the lower and upper limits of a (1-PROBTEST) * 100% confidence interval for the PSD estimates.

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if:

INITFFT = false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine POWER_SPECTRUM in order to sets up constants and functions for use by subroutine FFT which is called inside subroutine POWER_SPECTRUM. This call to INITFFT must have the following form:

call init_fft( (/ size(MAT,1), size(MAT,2)/2 /), dim=2_i4b )

INITFFT = true, the call to INIT_FFT is done inside subroutine POWER_SPECTRUM and a call to END_FFT is also done before leaving subroutine POWER_SPECTRUM

The default is INITFFT=true .

NORMPSD (INPUT, OPTIONAL) logical(lgl)

On entry, if:

NORMPSD = true, the PSD estimates are normalized in such a way that the total area under the power spectrum is equal to the variance of the time series MAT.

NORMPSD = false, the sum of the PSD estimates for each row of MAT (e.g. sum( PSMAT(:,2:), dim=2 ) is equal to the variance of the corresponding time series.

The default is NORMPSD=true .

SMOOTH_PARAM (INPUT, OPTIONAL) integer(i4b), dimension(:)

if SMOOTH_PARAM is used, the PSD estimates are computed by repeated smoothing of the periodogram with modified Daniell weights.

On entry, SMOOTH_PARAM(:) gives the array of the half-lengths of the modified Daniell filters to be applied.

All the values in SMOOTH_PARAM(:) must be greater than 0 and less than size(MAT,2)/2 + 1 .

Size(SMOOTH_PARAM) must be greater or equal to 1.

TREND (INPUT, OPTIONAL) integer(i4b)

If

TREND=+1 The means of the time series are removed before computing the spectra

TREND=+2 The drifts from the time series are removed before computing the spectra by using the formula: drift(i) = (MAT(i,size(MAT,2)) - MAT(i,1))/(size(MAT,2) - 1)

TREND=+3 The least-squares lines from the time series are removed before computing the spectra.

For other values of TREND nothing is done before estimating the power spectra.

The default is TREND=1, e.g. the means of the time series are removed before the computations.

WIN (INPUT, OPTIONAL) integer(i4b)

On entry, this argument specify the data window used in the computations of the power spectrum. If:

WIN=+1 The Bartlett window is used

WIN=+2 The square window is used

WIN=+3 The Welch window is used

WIN=+4 The Hann window is used

WIN=+5 The Hamming window is used

WIN=+6 A split-cosine-bell window is used

The default is WIN=3, e.g. the Welch window is used.

TAPERP (INPUT, OPTIONAL) real(stnd)

The total percentage of the data to be tapered if WIN=6. TAPERP must be greater than zero and less or equal to one, otherwise the default value is used.

The default is 0.2 .

PROBTEST (INPUT, OPTIONAL) real(stnd)

On entry, a probability. PROBTEST is the critical probability which is used to determine the lower and upper confidence limit factors (e.g. the optional arguments CONLWR and CONUPR).

PROBTEST must verify: 0. < P < 1.

The default is 0.05 .

Further Details¶

After removing the mean or the trend from the time series (e.g. TREND=1,2,3), the selected data window (e.g. WIN=1,2,3,4,5,6) is applied to the time series and the PSD estimates are computed by the FFT of these transformed time series. Optionally, theses PSD estimates may then be smoothed in the frequency domain by modified Daniell filters (e.g. if SMOOTH_PARAM is used).

The computed equivalent number of degrees of freedom and bandwidth must be divided by two for the zero and Nyquist frequencies.

Furthermore, the computed equivalent number of degrees of freedom, bandwidth, lower and upper (1-PROBTEST) * 100% confidence limit factors are not right near the zero and Nyquist frequencies if the PSD estimates have been smoothed by modified Daniell filters. The reason is that POWER_SPECTRUM assumes that smoothing involves averaging independent frequency ordinates. This is true except near the zero and Nyquist frequencies where an average may contain contributions from negative frequencies, which are identical to and hence not independent of positive frequency spectral values. Thus, the number of degrees of freedom in PSD estimates near the 0 and Nyquist frequencies are as little as half the number of degrees of freedom of the spectral estimates away from these frequency extremes if the optional argument SMOOTH_PARAM is used.

If the optional argument SMOOTH_PARAM is used, the computed equivalent number of degrees of freedom, bandwidth, lower and upper (1-PROBTEST) * 100% confidence limit factors are right for PSD estimates at frequencies

(i-1)/Size(MAT,2) for i= (nparam+1)/2 + 1 to ( Size(MAT,2) - nparam + 1)/2

where nparam = 2 * (2+sum(SMOOTH_PARAM(:)))- 1, (e.g. for frequencies i/Size(MAT,2) for i = (nparam+1)/2, … , ( Size(MAT,2) - nparam - 1)/2 ).

For definitions, more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York.

Welch, P.D., 1967:

The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE trans. on audio and electroacoustics, Vol. Au-15, 2, 70-73.

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford.

`subroutine cross_spectrum ( vec, vec2, psvec, psvec2, phase, coher, freq, edof, bandwidth, conlwr, conupr, testcoher, ampli, co_spect, quad_spect, prob_coher, initfft, normpsd, smooth_param, trend, win, taperp, probtest )`¶

Purpose¶

Subroutine CROSS_SPECTRUM computes Fast Fourier Transform (FFT) estimates of the power and cross spectra of two real time series, VEC and VEC2. The real valued sequences VEC and VEC2 must be of even length.

The Power Spectral Density (PSD) and Cross Spectral Density (CSD) estimates are returned in units which are the square of the data (if NORMPSD=false) or in spectral density units (if NORMPSD=true).

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the first real time series for which the power and cross spectra must be estimated. If WIN/=2 or TREND=1, 2 or 3, VEC is used as workspace and is transformed.

Size(VEC) must be an even (positive) integer greater or equal to 4.

VEC2 (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the second real time series for which the power and cross spectra must be estimated. If WIN/=2 or TREND=1, 2 or 3, VEC2 is used as workspace and is transformed.

VEC2 must verify: size(VEC2) = size(VEC).

PSVEC (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the Power Spectral Density (PSD) estimates of VEC.

PSVEC must verify: size(PSVEC) = size(VEC)/2 + 1 .

PSVEC2 (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC2)/2)+1 containing the Power Spectral Density (PSD) estimates of VEC2.

PSVEC2 must verify: size(PSVEC2) = size(VEC)/2 + 1 .

PHASE (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the phase of the cross spectrum, given in fractions of a circle (e.g. on the closed interval (0,1) ).

PHASE must verify: size(PHASE) = size(VEC)/2 + 1 .

COHER (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the squared coherency estimates for all frequencies.

COHER must verify: size(COHER) = size(VEC)/2 + 1 .

FREQ (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the frequencies at which the spectral quantities are calculated in cycles per unit of time.

The spectral estimates are taken at frequencies (i-1)/size(VEC) for i=1,2, … , (size(VEC)/2 + 1).

FREQ must verify: size(FREQ) = size(VEC)/2 + 1 .

EDOF (OUTPUT, OPTIONAL) real(stnd)

On exit, the equivalent number of degrees of freedom of the power and cross spectrum estimates.

BANDWIDTH (OUTPUT, OPTIONAL) real(stnd)

On exit, the bandwidth of the power and cross spectrum estimates.

CONLWR (OUTPUT, OPTIONAL) real(stnd)

CONUPR (OUTPUT, OPTIONAL) real(stnd)

On output, these arguments specify the lower and upper (1-PROBTEST) * 100% confidence limit factors, respectively. Multiply the PSD estimates (e.g. the PSVEC(:) and PSVEC2(:) arguments ) by these constants to get the lower and upper limits of a (1-PROBTEST) * 100% confidence interval for the PSD estimates.

TESTCOHER (OUTPUT, OPTIONAL) real(stnd)

On output, this argument specifies the critical value for testing the null hypothesis that the squared coherency is zero at the PROBTEST * 100% significance level (e.g. elements of COHER(:) less than TESTCOHER should be regarded as not significantly different from zero at the PROBTEST * 100% significance level).

AMPLI (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the cross-amplitude spectrum.

AMPLI must verify: size(AMPLI) = (size(VEC)/2) + 1 .

CO_SPECT (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the co-spectrum (e.g. the real part of cross-spectrum).

CO_SPECT must verify: size(CO_SPECT) = (size(VEC)/2) + 1 .

QUAD_SPECT (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the quadrature spectrum (e.g. the imaginary part of cross-spectrum with a minus sign).

QUAD_SPECT must verify: size(QUAD_SPECT) = (size(VEC)/2) + 1 .

PROB_COHER (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the probabilities that the computed sample squared coherencies came from an ergodic stationary bivariate process with (corresponding) squared coherencies equal to zero.

PROB_COHER must verify: size(PROB_COHER) = (size(VEC)/2) + 1 .

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if:

INITFFT = false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine CROSS_SPECTRUM in order to sets up constants and functions for use by subroutine FFT which is called inside subroutine CROSS_SPECTRUM. This call to INITFFT must have the following form:

call init_fft( size(VEC)/2 )

INITFFT = true, the call to INIT_FFT is done inside subroutine CROSS_SPECTRUM and a call to END_FFT is also done before leaving subroutine CROSS_SPECTRUM.

The default is INITFFT=true .

NORMPSD (INPUT, OPTIONAL) logical(lgl)

On entry, if:

NORMPSD = true, the power and cross spectra estimates are normalized in such a way that the total area under the power spectrum is equal to the variance of the time series VEC and VEC2.

NORMPSD = false, the sum of the PSD estimates (e.g. sum(PSVEC(2:)) and sum(PSVEC2(2:)) ) is equal to the variance of the corresponding time series.

The default is NORMPSD=true .

SMOOTH_PARAM (INPUT, OPTIONAL) integer(i4b), dimension(:)

if SMOOTH_PARAM is used, the power and cross spectra estimates are computed by repeated smoothing of the periodograms and cross-periodogram with modified Daniell weights.

On entry, SMOOTH_PARAM(:) gives the array of the half-lengths of the modified Daniell filters to be applied.

All the values in SMOOTH_PARAM(:) must be greater than 0 and less than size(VEC)/2+1 .

Size(SMOOTH_PARAM) must be greater or equal to 1.

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+1 The mean of the two time series is removed before computing the power and cross spectra.

TREND=+2 The drift from the two time series is removed before computing the power and cross spectra.

TREND=+3 The least-squares line from the two time series is removed before computing the power and cross spectra.

For other values of TREND nothing is done before estimating the power and cross spectra.

The default is TREND=1, e.g. the means of the time series are removed before the computations.

WIN (INPUT, OPTIONAL) integer(i4b)

On entry, this argument specify the data window used in the computations of the power and cross spectra. If

WIN=+1 The Bartlett window is used

WIN=+2 The square window is used

WIN=+3 The Welch window is used

WIN=+4 The Hann window is used

WIN=+5 The Hamming window is used

WIN=+6 A split-cosine-bell window is used

The default is WIN=3, e.g. the Welch window is used.

TAPERP (INPUT, OPTIONAL) real(stnd)

The total percentage of the data to be tapered if WIN=6. TAPERP must be greater than zero and less or equal to one, otherwise the default value is used.

The default is 0.2 .

PROBTEST (INPUT, OPTIONAL) real(stnd)

On entry, a probability. PROBTEST is the critical probability which is used to determine the lower and upper confidence limit factors (e.g. the optional arguments CONLWR and CONUPR ) and the critical value for testing the null hypothesis that the squared coherency is zero (e.g. the TESTCOHER optional argument).

PROBTEST must verify: 0. < P < 1.

The default is 0.05 .

Further Details¶

After removing the mean or the trend from the time series (e.g. TREND=1,2,3), the selected data window (e.g. WIN=1,2,3,4,5,6) is applied to the time series and the PSD and CSD estimates are computed by the FFT of these transformed time series. Optionally, theses PSD and CSD estimates may then be smoothed in the frequency domain by modified Daniell filters (e.g. if argument SMOOTH_PARAM is used).

The computed equivalent number of degrees of freedom and bandwidth must be divided by two for the zero and Nyquist frequencies.

Furthermore, the computed equivalent number of degrees of freedom, bandwidth, lower and upper (1-PROBTEST) * 100% confidence limit factors and critical value for the squared coherency (e.g. arguments EDOF, BANDWIDTH, CONLWR, CONUPR and TESTCOHER) are not right near the zero and Nyquist frequencies if the PSD estimates have been smoothed by modified Daniell filters. The reason is that CROSS_SPECTRUM assumes that smoothing involves averaging independent frequency ordinates. This is true except near the zero and Nyquist frequencies where an average may contain contributions from negative frequencies, which are identical to and hence not independent of positive frequency spectral values. Thus, the number of degrees of freedom in PSD estimates near the 0 and Nyquist frequencies are as little as half the number of degrees of freedom of the spectral estimates away from these frequency extremes if the optional argument SMOOTH_PARAM is used.

If the optional argument SMOOTH_PARAM is used, the computed equivalent number of degrees of freedom, bandwidth, lower and upper (1-PROBTEST) * 100% confidence limit factors and critical value for the squared coherency are right for PSD estimates at frequencies

(i-1)/size(VEC) for i= (nparam+1)/2 + 1 to ( size(VEC) - nparam + 1)/2

where nparam = 2 * (2+sum(SMOOTH_PARAM(:)))- 1, (e.g. for frequencies i/size(VEC) for i = (nparam+1)/2, … , (size(VEC)-nparam-1)/2 ) .

For definitions, more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York.

Welch, P.D., 1967:

The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE trans. on audio and electroacoustics, Vol. Au-15, 2, 70-73.

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford.

`subroutine cross_spectrum ( vec, mat, psvec, psmat, phase, coher, freq, edof, bandwidth, conlwr, conupr, testcoher, ampli, co_spect, quad_spect, prob_coher, initfft, normpsd, smooth_param, trend, win, taperp, probtest )`¶

Purpose¶

Subroutine CROSS_SPECTRUM computes Fast Fourier Transform (FFT) estimates of the power and cross spectra of the real time series, VEC, and the multi-channel real time series MAT.

The Power Spectral Density (PSD) and Cross Spectral Density (CSD) estimates are returned in units which are the square of the data (if NORMPSD=false) or in spectral density units (if NORMPSD=true).

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the real time series for which the power and cross spectra must be estimated. If WIN/=2 or TREND=1, 2 or 3, VEC is used as workspace and is transformed.

Size(VEC) must be an even (positive) integer greater or equal to 4.

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the multi-channel real time series for which the power and cross spectra must be estimated. Each row of MAT is a real time series. If WIN/=2 or TREND=1, 2 or 3, MAR is used as workspace and is transformed.

The shape of MAT must verify: size(MAT,2) = size(VEC).

PSVEC (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the Power Spectral Density (PSD) estimates of VEC.

PSVEC must verify: size(PSVEC) = size(VEC)/2 + 1 .

PSMAT (OUTPUT) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and (size(VEC)/2)+1 columns containing the Power Spectral Density (PSD) estimates of each row of MAT.

The shape of PSMAT must verify:

size(PSMAT,1) = size(MAT,1) ;

size(PSMAT,2) = size(VEC)/2 + 1 .

PHASE (OUTPUT) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and (size(VEC)/2)+1 columns containing the phase of the cross spectrum, given in fractions of a circle (e.g. on the closed interval (0,1) ).

The shape of PHASE must verify:

size(PHASE,1) = size(MAT,1) ;

size(PHASE,2) = size(VEC)/2 + 1 .

COHER (OUTPUT) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and (size(VEC)/2)+1 columns containing the squared coherency estimates for all frequencies.

The shape of COHER must verify:

size(COHER,1) = size(MAT,1) ;

size(COHER,2) = size(VEC)/2 + 1 .

FREQ (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length (size(VEC)/2)+1 containing the frequencies at which the spectral quantities are calculated in cycles per unit of time.

The spectral estimates are taken at frequencies (i-1)/size(VEC) for i=1,2, … , (size(VEC)/2 + 1).

FREQ must verify: size(FREQ) = size(VEC)/2 + 1 .

EDOF (OUTPUT, OPTIONAL) real(stnd)

On exit, the equivalent number of degrees of freedom of the power and cross spectrum estimates.

BANDWIDTH (OUTPUT, OPTIONAL) real(stnd)

On exit, the bandwidth of the power and cross spectrum estimates.

CONLWR (OUTPUT, OPTIONAL) real(stnd)

CONUPR (OUTPUT, OPTIONAL) real(stnd)

On output, these arguments specify the lower and upper (1-PROBTEST) * 100% confidence limit factors, respectively. Multiply the PSD estimates (e.g. the PSVEC(:) and PSMAT(:,:) arguments ) by these constants to get the lower and upper limits of a (1-PROBTEST) * 100% confidence interval for the PSD estimates.

TESTCOHER (OUTPUT, OPTIONAL) real(stnd)

On output, this argument specifies the critical value for testing the null hypothesis that the squared coherency is zero at the PROBTEST * 100% significance level (e.g. elements of COHER(:,:) less than TESTCOHER should be regarded as not significantly different from zero at the PROBTEST * 100% significance level).

AMPLI (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and (size(VEC)/2)+1 columns containing the cross-amplitude spectra.

The shape of AMPLI must verify:

size(AMPLI,1) = size(MAT,1) ;

size(AMPLI,2) = (size(VEC)/2) + 1 .

CO_SPECT (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and (size(VEC)/2)+1 columns containing the co-spectra (e.g. the real part of cross-spectra).

The shape of CO_SPECT must verify:

size(CO_SPECT,1) = size(MAT,1) ;

size(CO_SPECT,2) = (size(VEC)/2) + 1 .

QUAD_SPECT (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and (size(VEC)/2)+1 columns containing the quadrature spectrum (e.g. the imaginary part of cross-spectrum with a minus sign).

The shape of QUAD_SPECT must verify:

size(QUAD_SPECT,1) = size(MAT,1) ;

size(QUAD_SPECT,2) = (size(VEC)/2) + 1 .

PROB_COHER (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and (size(VEC)/2)+1 columns containing the probabilities that the computed sample squared coherencies came from an ergodic stationary bivariate process with (corresponding) squared coherencies equal to zero.

The shape of PROB_COHER must verify:

size(PROB_COHER,1) = size(MAT,1) ;

size(PROB_COHER,2) = (size(VEC)/2) + 1 .

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if:

INITFFT = false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine CROSS_SPECTRUM in order to sets up constants and functions for use by subroutine FFT which is called inside subroutine CROSS_SPECTRUM. This call to INITFFT must have the following form:

call init_fft( (/ size(MAT,1), size(MAT,2)/2 /), dim=2_i4b )

INITFFT = true, the call to INIT_FFT is done inside subroutine CROSS_SPECTRUM and a call to END_FFT is also done before leaving subroutine CROSS_SPECTRUM.

The default is INITFFT=true .

NORMPSD (INPUT, OPTIONAL) logical(lgl)

On entry, if:

NORMPSD = true, the power and cross spectra estimates are normalized in such a way that the total area under the power spectra is equal to the variance of the time series contained in VEC and in each row of MAT.

NORMPSD = false, the sum of the PSD estimates (e.g. sum(PSVEC(2:)) and sum(PSMAT(:,2:),dim=2) ) is equal to the variance of the corresponding time series.

The default is NORMPSD=true .

SMOOTH_PARAM (INPUT, OPTIONAL) integer(i4b), dimension(:)

if SMOOTH_PARAM is used, the power and cross spectra estimates are computed by repeated smoothing of the periodograms and cross-periodogram with modified Daniell weights.

On entry, SMOOTH_PARAM(:) gives the array of the half-lengths of the modified Daniell filters to be applied.

All the values in SMOOTH_PARAM(:) must be greater than 0 and less than (size(VEC)/2)+1 .

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+1 The means of the time series are removed before computing the power and cross spectra

TREND=+2 The drifts from time series are removed before computing the power and cross spectra

TREND=+3 The least-squares lines from time series are removed before computing the power and cross spectra.

For other values of TREND nothing is done before estimating the power and cross spectra.

The default is TREND=1, e.g. the means of the time series are removed before the computations.

WIN (INPUT, OPTIONAL) integer(i4b)

On entry, this argument specify the data window used in the computations of the power and cross spectra. If:

WIN=+1 The Bartlett window is used

WIN=+2 The square window is used

WIN=+3 The Welch window is used

WIN=+4 The Hann window is used

WIN=+5 The Hamming window is used

WIN=+6 A split-cosine-bell window is used

The default is WIN=3, e.g. the Welch window is used.

TAPERP (INPUT, OPTIONAL) real(stnd)

The total percentage of the data to be tapered if WIN=6. TAPERP must be greater than zero and less or equal to one, otherwise the default value is used.

The default is 0.2 .

PROBTEST (INPUT, OPTIONAL) real(stnd)

On entry, a probability. PROBTEST is the critical probability which is used to determine the lower and upper confidence limit factors (e.g. the optional arguments CONLWR and CONUPR ) and the critical value for testing the null hypothesis that the squared coherency is zero (e.g. the TESTCOHER optional argument).

PROBTEST must verify: 0. < P < 1.

The default is 0.05 .

Further Details¶

After removing the mean or the trend from the time series (e.g. TREND=1,2,3), the selected data window (e.g. WIN=1,2,3,4,5,6) is applied to the time series and the PSD and CSD estimates are computed by the FFT of these transformed time series. Optionally, theses PSD and CSD estimates may then be smoothed in the frequency domain by modified Daniell filters (e.g. if argument SMOOTH_PARAM is used).

The computed equivalent number of degrees of freedom and bandwidth must be divided by two for the zero and Nyquist frequencies.

Furthermore, the computed equivalent number of degrees of freedom, bandwidth, lower and upper (1-PROBTEST) * 100% confidence limit factors and critical value for the squared coherency (e.g. arguments EDOF, BANDWIDTH, CONLWR, CONUPR and TESTCOHER) are not right near the zero and Nyquist frequencies if the PSD estimates have been smoothed by modified Daniell filters. The reason is that CROSS_SPECTRUM assumes that smoothing involves averaging independent frequency ordinates. This is true except near the zero and Nyquist frequencies where an average may contain contributions from negative frequencies, which are identical to and hence not independent of positive frequency spectral values. Thus, the number of degrees of freedom in PSD estimates near the 0 and Nyquist frequencies are as little as half the number of degrees of freedom of the spectral estimates away from these frequency extremes if the optional argument SMOOTH_PARAM is used.

If the optional argument SMOOTH_PARAM is used, the computed equivalent number of degrees of freedom, bandwidth, lower and upper (1-PROBTEST) * 100% confidence limit factors and critical value for the squared coherency are right for PSD estimates at frequencies

(i-1)/size(VEC) for i= (nparam+1)/2 + 1 to ( size(VEC) - nparam + 1)/2

where nparam = 2 * (2+sum(SMOOTH_PARAM(:)))- 1, (e.g. for frequencies i/size(VEC) for i = (nparam+1)/2, … , (size(VEC)-nparam-1)/2 ) .

For definitions, more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York.

Welch, P.D., 1967:

The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE trans. on audio and electroacoustics, Vol. Au-15, 2, 70-73.

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford.

`subroutine power_spectrum2 ( vec, l, psvec, freq, edof, bandwidth, conlwr, conupr, initfft, overlap, normpsd, smooth_param, trend, trend2, win, taperp, l0, probtest )`¶

Purpose¶

Subroutine POWER_SPECTRUM2 computes a Fast Fourier Transform (FFT) estimate of the power spectrum of a real time series.

The Power Spectral Density (PSD) estimates are returned in units which are the square of the data (if NORMPSD=false) or in spectral density units (if NORMPSD=true).

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the real time series for which the power spectrum must be estimated. If TREND=1, 2 or 3, VEC is used as workspace and is transformed.

Size(VEC) must be greater or equal to 4.

L (INPUT) integer(i4b)

On entry, an integer used to segment the time series. L is the length of the segments. L must be a positive even integer, less or equal to size(VEC), but greater or equal to 4.

Spectral computations are at (L/2)+1 frequencies if the optional argument L0 is absent and are at ((L+L0)/2)+1 frequencies if L0 is present (L0 is the number of zeros added to each segment).

Suggested values for L+L0 are 16, 32, 64 or 128 (e.g. an integer power of two, in order to speed the computations).

PSVEC (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the Power Spectral Density (PSD) estimates of VEC.

PSVEC must verify: size(PSVEC) = ((L+L0)/2) + 1 .

FREQ (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the frequencies at which the spectral quantities are calculated in cycles per unit of time.

The spectral estimates are taken at frequencies (i-1)/(L+L0) for i=1,2, … , ((L+L0)/2 + 1).

FREQ must verify: size(FREQ) = (L+L0)/2 + 1 .

EDOF (OUTPUT, OPTIONAL) real(stnd)

On exit, the equivalent number of degrees of freedom of the power spectrum estimates.

BANDWIDTH (OUTPUT, OPTIONAL) real(stnd)

On exit, the bandwidth of the power spectrum estimates.

CONLWR (OUTPUT, OPTIONAL) real(stnd)

CONUPR (OUTPUT, OPTIONAL) real(stnd)

On output, these arguments specify the lower and upper (1-PROBTEST) * 100% confidence limit factors, respectively. Multiply the PSD estimates (e.g. the PSVEC(:) argument) by these constants to get the lower and upper limits of a (1-PROBTEST) * 100% confidence interval for the PSD estimates.

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if:

INITFFT = false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine POWER_SPECTRUM2 in order to sets up constants and functions for use by subroutine FFT which is called inside subroutine POWER_SPECTRUM2. This call to INITFFT must have the following form:

call init_fft( (L+L0)/2 )

INITFFT = true, the call to INIT_FFT is done inside subroutine POWER_SPECTRUM2 and a call to END_FFT is also done before leaving subroutine POWER_SPECTRUM2.

The default is INITFFT=true .

OVERLAP (INPUT, OPTIONAL) logical(lgl)

If:

OVERLAP = false, the subroutine segments the data without any overlapping.

OVERLAP = true, the subroutine overlaps the segments by one half of their length (which is equal to L).

In both cases, zeros are eventually added to each segment (if argument L0 is present) and each segment will be FFT’d, and the resulting periodograms will averaged together to obtain a Power Spectrum Density estimate at the ((L+L0)/2)+1 frequencies.

The default is OVERLAP=false .

NORMPSD (INPUT, OPTIONAL) logical(lgl)

On entry, if:

NORMPSD = true, the PSD estimates are normalized in such a way that the total area under the power spectrum is equal to the variance of the time series VEC.

NORMPSD = false, the sum of the PSD estimates (e.g. sum( PSVEC(2:) ) is equal to the variance of the time series.

The default is NORMPSD=true .

SMOOTH_PARAM (INPUT, OPTIONAL) integer(i4b), dimension(:)

If SMOOTH_PARAM is used, the PSD estimates are computed by repeated smoothing of the periodogram with modified Daniell weights.

On entry, SMOOTH_PARAM(:) gives the array of the half-lengths of the modified Daniell filters to be applied.

All the values in SMOOTH_PARAM(:) must be greater than 0 and less than ((L+L0)/2) + 1 .

Size(SMOOTH_PARAM) must be greater or equal to 1.

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+1 The mean of the time series is removed before computing the spectrum

TREND=+2 The drift from the time series is removed before computing the spectrum by using the formula: drift = (VEC(size(VEC)) - VEC(1))/(size(VEC) - 1)

TREND=+3 The least-squares line from the time series is removed before computing the spectrum.

For other values of TREND nothing is done before estimating the power spectrum.

The default is TREND=1, e.g. the mean of the time series is removed before the computations.

TREND2 (INPUT, OPTIONAL) integer(i4b)

If:

TREND2=+1 The mean of the time segment is removed before computing the spectrum on this segment.

TREND2=+2 The drift from the time segment is removed before computing the spectrum on this segment.

TREND2=+3 The least-squares line from the time segment is removed before computing the spectrum on this segment.

For other values of TREND2 nothing is done before estimating the power spectrum on each segment.

The default is TREND2=0, e.g. nothing is done before estimating the power spectrum on each segment.

WIN (INPUT, OPTIONAL) integer(i4b)

On entry, this argument specify the data window used in the computations of the power spectrum. If:

WIN=+1 The Bartlett window is used

WIN=+2 The square window is used

WIN=+3 The Welch window is used

WIN=+4 The Hann window is used

WIN=+5 The Hamming window is used

WIN=+6 A split-cosine-bell window is used

The default is WIN=3, e.g. the Welch window is used.

TAPERP (INPUT, OPTIONAL) real(stnd)

The total percentage of the data to be tapered if WIN=6. TAPERP must be greater than zero and less or equal to one, otherwise the default value is used.

The default is 0.2 .

L0 (INPUT, OPTIONAL) integer(i4b)

The number of zeros added to each time segment in order to obtain more finely spaced spectral estimates. L+L0 must be a positive even integer.

The default is L0=0, e.g. no zeros are added to each time segment.

PROBTEST (INPUT, OPTIONAL) real(stnd)

On entry, a probability. PROBTEST is the critical probability which is used to determine the lower and upper confidence limit factors (e.g. the optional arguments CONLWR and CONUPR).

PROBTEST must verify: 0. < P < 1.

The default is 0.05 .

Further Details¶

After removing the mean or the trend from the time series (e.g. TREND=1,2,3), the series is padded with zero on the right such that the length of the resulting time series is evenly divisible by L (a positive even integer). The length, N, of this resulting time series is the first integer greater than or equal to size(VEC) which is evenly divisible by L. If size(VEC) is not evenly divisible by L, N is equal to size(VEC)+L-mod(size(VEC),L).

Optionally, the mean or the trend may also be removed from each time segment (e.g. TREND2=1,2,3). Optionally, zeros may be added to each time segment (e.g. the optional arguemnt L0) if more finely spaced spectral esimates are desired.

The stability of the PSD estimates depends on the averaging process. That is, the greater the number of segments ( N/L if OVERLAP=false and (2N/L)-1 if OVERLAP=true), the more stable the resulting PSD estimates.

Optionally, theses PSD estimates may then be smoothed again in the frequency domain by modified Daniell filters (e.g. if argument SMOOTH_PARAM is used).

The computed equivalent number of degrees of freedom and bandwidth must be divided by two for the zero and Nyquist frequencies.

Furthermore, the computed equivalent number of degrees of freedom, bandwidth, lower and upper (1-PROBTEST) * 100% confidence limit factors are not right near the zero and Nyquist frequencies if the PSD estimates have been smoothed by modified Daniell filters. The reason is that POWER_SPECTRUM2 assumes that smoothing involves averaging independent frequency ordinates. This is true except near the zero and Nyquist frequencies where an average may contain contributions from negative frequencies, which are identical to and hence not independent of positive frequency spectral values. Thus, the number of degrees of freedom in PSD estimates near the 0 and Nyquist frequencies are as little as half the number of degrees of freedom of the spectral estimates away from these frequency extremes if the optional argument SMOOTH_PARAM is used.

If the optional argument SMOOTH_PARAM is used, the computed equivalent number of degrees of freedom, bandwidth, lower and upper (1-PROBTEST) * 100% confidence limit factors are right for PSD estimates at frequencies

(i-1)/(L+L0) for i= (nparam+1)/2 + 1 to ( (L+L0) - nparam + 1)/2

where nparam = 2 * (2+sum(SMOOTH_PARAM(:)))- 1, (e.g. for frequencies i/(L+L0) for i = (nparam+1)/2, … , ( (L+L0) - nparam - 1)/2 ).

For definitions, more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York.

Welch, P.D., 1967:

The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE trans. on audio and electroacoustics, Vol. Au-15, 2, 70-73.

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford.

`subroutine power_spectrum2 ( mat, l, psmat, freq, edof, bandwidth, conlwr, conupr, initfft, overlap, normpsd, smooth_param, trend, trend2, win, taperp, l0, probtest )`¶

Purpose¶

Subroutine POWER_SPECTRUM2 computes Fast Fourier Transform (FFT) estimates of the power spectra of the multi-channel real time series MAT (e.g. each row of MAT contains a time series).

The Power Spectral Density (PSD) estimates are returned in units which are the square of the data (if NORMPSD=false) or in spectral density units (if NORMPSD=true).

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the multi-channel real time series for which the power spectra must be estimated. Each row of MAT is a real time series. If TREND=1, 2 or 3, MAT is used as workspace and is transformed.

Size(MAT,2) must be greater or equal to 4.

L (INPUT) integer(i4b)

On entry, an integer used to segment the time series. L is the length of the segments. L must be a positive even integer less or equal to size(MAT,2), but greater or equal to 4.

Spectral computations are at (L/2)+1 frequencies if the optional argument L0 is absent and are at ((L+L0)/2)+1 frequencies if L0 is present (L0 is the number of zeros added to each segment).

Suggested values for L+L0 are 16, 32, 64 or 128 (e.g. an integer power of two, in order to speed the computations).

PSMAT (OUTPUT) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and ((L+L0)/2) + 1 columns containing the Power Spectral Density (PSD) estimates of each row of MAT.

The shape of PSMAT must verify:

size(PSMAT,1) = size(MAT,1) ;

size(PSMAT,2) = ((L+L0)/2) + 1 .

FREQ (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the frequencies at which the spectral quantities are calculated in cycles per unit of time.

The spectral estimates are taken at frequencies (i-1)/(L+L0) for i=1,2, … , ((L+L0)/2 + 1).

FREQ must verify: size(FREQ) = (L+L0)/2 + 1 .

EDOF (OUTPUT, OPTIONAL) real(stnd)

On exit, the equivalent number of degrees of freedom of the power spectrum estimates.

BANDWIDTH (OUTPUT, OPTIONAL) real(stnd)

On exit, the bandwidth of the power spectrum estimates.

CONLWR (OUTPUT, OPTIONAL) real(stnd)

CONUPR (OUTPUT, OPTIONAL) real(stnd)

On output, these arguments specify the lower and upper (1-PROBTEST) * 100% confidence limit factors, respectively. Multiply the PSD estimates (e.g. the PSMAT(:,:) argument) by these constants to get the lower and upper limits of a (1-PROBTEST) * 100% confidence interval for the PSD estimates.

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if:

INITFFT = false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine POWER_SPECTRUM2 in order to sets up constants and functions for use by subroutine FFT which is called inside subroutine POWER_SPECTRUM2. This call to INITFFT must have the following form:

call init_fft( (/ size(MAT,1), (L+L0)/2 /), dim=2_i4b )

INITFFT = true, the call to INIT_FFT is done inside subroutine POWER_SPECTRUM2 and a call to END_FFT is also done before leaving subroutine POWER_SPECTRUM2.

The default is INITFFT=true .

OVERLAP (INPUT, OPTIONAL) logical(lgl)

If:

OVERLAP = false, the subroutine segments the data without any overlapping.

OVERLAP = true, the subroutine overlaps the segments by one half of their length (which is equal to L).

In both cases, zeros are eventually added to each segment (if argument L0 is present) and each segment will be FFT’d, and the resulting periodograms will averaged together to obtain a Power Spectrum Density estimate at the ((L+L0)/2)+1 frequencies.

The default is OVERLAP=false .

NORMPSD (INPUT, OPTIONAL) logical(lgl)

On entry, if:

NORMPSD = true, the PSD estimates are normalized in such a way that the total area under the power spectrum is equal to the variance of the corresponding time series in MAT.

NORMPSD = false, the sum of the PSD estimates (e.g. sum(PSMAT(:,2:),dim=2) ) is equal to the variance of the corresponding time series.

The default is NORMPSD=true .

SMOOTH_PARAM (INPUT, OPTIONAL) integer(i4b), dimension(:)

If SMOOTH_PARAM is used, the PSD estimates are computed by repeated smoothing of the periodogram with modified Daniell weights.

On entry, SMOOTH_PARAM(:) gives the array of the half-lengths of the modified Daniell filters to be applied.

All the values in SMOOTH_PARAM(:) must be greater than 0 and less than ((L+L0)/2)+1 .

Size(SMOOTH_PARAM) must be greater or equal to 1.

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+1 The means of the time series are removed before computing the spectra

TREND=+2 The drifts from time series are removed before computing the spectra

TREND=+3 The least-squares lines from time series are removed before computing the spectra.

For other values of TREND nothing is done before estimating the power and cross spectra.

The default is TREND=1, e.g. the means of the time series are removed before the computations.

TREND2 (INPUT, OPTIONAL) integer(i4b)

If:

TREND2=+1 The mean of the time segment is removed before computing the spectrum on this segment.

TREND2=+2 The drift from the time segment is removed before computing the spectrum on this segment.

TREND2=+3 The least-squares line from the time segment is removed before computing the spectrum on this segment.

For other values of TREND2 nothing is done before estimating the power spectrum on each segment.

The default is TREND2=0, e.g. nothing is done before estimating the power spectrum on each segment.

WIN (INPUT, OPTIONAL) integer(i4b)

On entry, this argument specify the data window used in the computations of the power spectrum. If:

WIN=+1 The Bartlett window is used

WIN=+2 The square window is used

WIN=+3 The Welch window is used

WIN=+4 The Hann window is used

WIN=+5 The Hamming window is used

WIN=+6 A split-cosine-bell window is used

The default is WIN=3, e.g. the Welch window is used.

TAPERP (INPUT, OPTIONAL) real(stnd)

The total percentage of the data to be tapered if WIN=6. TAPERP must be greater than zero and less or equal to one, otherwise the default value is used.

The default is 0.2 .

L0 (INPUT, OPTIONAL) integer(i4b)

The number of zeros added to each time segment in order to obtain more finely spaced spectral estimates. L+L0 must be a positive even integer.

The default is L0=0, e.g. no zeros are added to each time segment.

PROBTEST (INPUT, OPTIONAL) real(stnd)

On entry, a probability. PROBTEST is the critical probability which is used to determine the lower and upper confidence limit factors (e.g. the optional arguments CONLWR and CONUPR).

PROBTEST must verify: 0. < P < 1.

The default is 0.05 .

Further Details¶

After removing the mean or the trend from the time series (e.g. TREND=1,2,3), the series are padded with zero on the right such that the length of the resulting time series is evenly divisible by L (a positive even integer). The length, N, of this resulting time series is the first integer greater than or equal to size(MAT,2) which is evenly divisible by L. If size(MAT,2) is not evenly divisible by L, N is equal to size(MAT,2)+L-mod(size(MAT,2),L).

Optionally, the mean or the trend may also be removed from each time segment (e.g. TREND2=1,2,3). Optionally, zeros may be added to each time segment (e.g. the optional arguemnt L0) if more finely spaced spectral esimates are desired.

The stability of the PSD estimates depends on the averaging process. That is, the greater the number of segments ( N/L if OVERLAP=false and (2N/L)-1 if OVERLAP=true), the more stable the resulting PSD estimates.

Optionally, theses PSD estimates may then be smoothed again in the frequency domain by modified Daniell filters (e.g. if argument SMOOTH_PARAM is used).

The computed equivalent number of degrees of freedom and bandwidth must be divided by two for the zero and Nyquist frequencies.

Furthermore, the computed equivalent number of degrees of freedom, bandwidth, lower and upper (1-PROBTEST) * 100% confidence limit factors are not right near the zero and Nyquist frequencies if the PSD estimates have been smoothed by modified Daniell filters. The reason is that POWER_SPECTRUM2 assumes that smoothing involves averaging independent frequency ordinates. This is true except near the zero and Nyquist frequencies where an average may contain contributions from negative frequencies, which are identical to and hence not independent of positive frequency spectral values. Thus, the number of degrees of freedom in PSD estimates near the 0 and Nyquist frequencies are as little as half the number of degrees of freedom of the spectral estimates away from these frequency extremes if the optional argument SMOOTH_PARAM is used.

If the optional argument SMOOTH_PARAM is used, the computed equivalent number of degrees of freedom, bandwidth, lower and upper (1-PROBTEST) * 100% confidence limit factors are right for PSD estimates at frequencies

(i-1)/(L+L0) for i= (nparam+1)/2 + 1 to ( (L+L0) - nparam + 1)/2

where nparam = 2 * (2+sum(SMOOTH_PARAM(:)))- 1, (e.g. for frequencies i/(L+L0) for i = (nparam+1)/2, … , ( (L+L0) - nparam - 1)/2 ).

For definitions, more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York.

Welch, P.D., 1967:

The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE trans. on audio and electroacoustics, Vol. Au-15, 2, 70-73.

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford.

`subroutine cross_spectrum2 ( vec, vec2, l, psvec, psvec2, phase, coher, freq, edof, bandwidth, conlwr, conupr, testcoher, ampli, co_spect, quad_spect, prob_coher, initfft, overlap, normpsd, smooth_param, trend, trend2, win, taperp, l0, probtest )`¶

Purpose¶

Subroutine CROSS_SPECTRUM2 computes Fast Fourier Transform (FFT) estimates of the power and cross spectra of two real time series.

The Power Spectral Density (PSD) and Cross Spectral Density (CSD) estimates are returned in units which are the square of the data (if NORMPSD=false) or in spectral density units (if NORMPSD=true).

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the first real time series for which the power and cross spectra must be estimated. If TREND=1, 2 or 3, VEC is used as workspace and is transformed.

Size(VEC) must be greater or equal to 4.

VEC2 (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the second real time series for which the power and cross spectra must be estimated. If TREND=1, 2 or 3, VEC2 is used as workspace and is transformed.

VEC2 must verify: size(VEC2) = size(VEC).

L (INPUT) integer(i4b)

On entry, an integer used to segment the time series. L is the length of the segments. L must be a positive even integer, less or equal to size(VEC), but greater or equal to 4.

Spectral computations are at (L/2)+1 frequencies if the optional argument L0 is absent and are at ((L+L0)/2)+1 frequencies if L0 is present (L0 is the number of zeros added to each segment).

Suggested values for L+L0 are 16, 32, 64 or 128 (e.g. an integer power of two, in order to speed the computations).

PSVEC (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the Power Spectral Density (PSD) estimates of VEC.

PSVEC must verify: size(PSVEC) = ((L+L0)/2) + 1 .

PSVEC2 (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the Power Spectral Density (PSD) estimates of VEC2.

PSVEC2 must verify: size(PSVEC2) = ((L+L0)/2) + 1 .

PHASE (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the phase of the cross spectrum, given in fractions of a circle (e.g. on the closed interval (0,1) ).

PHASE must verify: size(PHASE) = ((L+L0)/2) + 1 .

COHER (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the squared coherency estimates for all frequencies.

COHER must verify: size(COHER) = ((L+L0)/2) + 1 .

FREQ (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the frequencies at which the spectral quantities are calculated in cycles per unit of time.

The spectral estimates are taken at frequencies (i-1)/(L+L0) for i=1,2, … , ((L+L0)/2 + 1).

FREQ must verify: size(FREQ) = (L+L0)/2 + 1 .

EDOF (OUTPUT, OPTIONAL) real(stnd)

On exit, the equivalent number of degrees of freedom of the power and cross spectrum estimates.

BANDWIDTH (OUTPUT, OPTIONAL) real(stnd)

On exit, the bandwidth of the power and cross spectrum estimates.

CONLWR (OUTPUT, OPTIONAL) real(stnd)

CONUPR (OUTPUT, OPTIONAL) real(stnd)

On output, these arguments specify the lower and upper (1-PROBTEST) * 100% confidence limit factors, respectively. Multiply the PSD estimates (e.g. the PSVEC(:) and PSVEC2(:) arguments ) by these constants to get the lower and upper limits of a (1-PROBTEST) * 100% confidence interval for the PSD estimates.

TESTCOHER (OUTPUT, OPTIONAL) real(stnd)

On output, this argument specifies the critical value for testing the null hypothesis that the squared coherency is zero at the PROBTEST * 100% significance level (e.g. elements of COHER(:) less than TESTCOHER should be regarded as not significantly different from zero at the PROBTEST * 100% significance level).

AMPLI (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the cross-amplitude spectrum.

AMPLI must verify: size(AMPLI) = ((L+L0)/2) + 1 .

CO_SPECT (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the co-spectrum (e.g. the real part of cross-spectrum).

CO_SPECT must verify: size(CO_SPECT) = ((L+L0)/2) + 1 .

QUAD_SPECT (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the quadrature spectrum (e.g. the imaginary part of cross-spectrum with a minus sign).

QUAD_SPECT must verify: size(QUAD_SPECT) = ((L+L0)/2) + 1 .

PROB_COHER (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the probabilities that the computed sample squared coherencies came from an ergodic stationary bivariate process with (corresponding) squared coherencies equal to zero.

PROB_COHER must verify: size(PROB_COHER) = ((L+L0)/2)+1 .

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if INITFFT is set to false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine CROSS_SPECTRUM2 in order to sets up constants and functions for use by subroutine FFT which is called inside subroutine CROSS_SPECTRUM2 (the call to INITFFT must have the following form:

call init_fft( (L+L0)/2 )

If INITFFT is set to true, the call to INIT_FFT is done inside subroutine CROSS_SPECTRUM2 and a call to END_FFT is also done before leaving subroutine CROSS_SPECTRUM2.

The default is INITFFT=true .

OVERLAP (INPUT, OPTIONAL) logical(lgl)

If:

OVERLAP = false, the subroutine segments the data without any overlapping.

OVERLAP = true, the subroutine overlaps the segments by one half of their length (which is equal to L).

In both cases, zeros are eventually added to each segment (if argument L0 is present) and each segment will be FFT’d, and the resulting periodograms will averaged together to obtain a Power Spectrum Density estimate at the ((L+L0)/2)+1 frequencies.

The default is OVERLAP=false .

NORMPSD (INPUT, OPTIONAL) logical(lgl)

On entry, if NORMPSD is set to true, the power and cross spectra estimates are normalized in such a way that the total area under the power spectrum is equal to the variance of the time series VEC and VEC2. If NORMPSD is set to false, the sum of the PSD estimates (e.g. sum(PSVEC(2:)) and sum(PSVEC2(2:)) ) is equal to the variance of the corresponding time series.

The default is NORMPSD=true .

SMOOTH_PARAM (INPUT, OPTIONAL) integer(i4b), dimension(:)

If SMOOTH_PARAM is used, the power and cross spectra estimates are computed by repeated smoothing of the periodograms and cross-periodogram with modified Daniell weights.

On entry, SMOOTH_PARAM(:) gives the array of the half-lengths of the modified Daniell filters to be applied.

All the values in SMOOTH_PARAM(:) must be greater than 0 and less than (L+L0)/2 + 1 .

Size(SMOOTH_PARAM) must be greater or equal to 1.

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+1 The mean of the two time series is removed before computing the spectra

TREND=+2 The drift from the two time series is removed before computing the spectra

TREND=+3 The least-squares line from the two time series is removed before computing the spectra.

For other values of TREND nothing is done before estimating the power and cross spectra.

The default is TREND=1, e.g. the means of the time series are removed before the computations.

TREND2 (INPUT, OPTIONAL) integer(i4b)

If:

TREND2=+1 The mean of the time segment is removed before computing the cross-spectrum on this segment.

TREND2=+2 The drift from the time segment is removed before computing the cross-spectrum on this segment.

TREND2=+3 The least-squares line from the time segment is removed before computing the cross-spectrum on this segment.

For other values of TREND2 nothing is done before estimating the cross-spectrum on each segment.

The default is TREND2=0, e.g. nothing is done before estimating the power spectrum on each segment.

WIN (INPUT, OPTIONAL) integer(i4b)

On entry, this argument specify the data window used in the computations of the power and cross spectra. If:

WIN=+1 The Bartlett window is used

WIN=+2 The square window is used

WIN=+3 The Welch window is used

WIN=+4 The Hann window is used

WIN=+5 The Hamming window is used

WIN=+6 A split-cosine-bell window is used

The default is WIN=3, e.g. the Welch window is used.

TAPERP (INPUT, OPTIONAL) real(stnd)

The total percentage of the data to be tapered if WIN=6. TAPERP must be greater than zero and less or equal to one, otherwise the default value is used.

The default is 0.2 .

L0 (INPUT, OPTIONAL) integer(i4b)

The number of zeros added to each time segment in order to obtain more finely spaced spectral estimates. L+L0 must be a positive even integer.

The default is L0=0, e.g. no zeros are added to each time segment.

PROBTEST (INPUT, OPTIONAL) real(stnd)

On entry, a probability. PROBTEST is the critical probability which is used to determine the lower and upper confidence limit factors (e.g. the optional arguments CONLWR and CONUPR ) and the critical value for testing the null hypothesis that the squared coherency is zero (e.g. the TESTCOHER optional argument).

PROBTEST must verify: 0. < P < 1.

The default is 0.05 .

Further Details¶

After removing the mean or the trend from the two time series (e.g. TREND=1,2,3), the series are padded with zero on the right such that the length of the resulting two time series is evenly divisible by L (a positive even integer). The length, N, of these resulting time series is the first integer greater than or equal to size(VEC) which is evenly divisible by L. If size(VEC) is not evenly divisible by L, N is equal to size(VEC)+L-mod(size(VEC),L).

Optionally, the mean or the trend may also be removed from each time segment (e.g. TREND2=1,2,3). Optionally, zeros may be added to each time segment (e.g. the optional arguemnt L0) if more finely spaced spectral esimates are desired.

The stability of the power and cross spectra estimates depends on the averaging process. That is, the greater the number of segments ( N/L if OVERLAP=false and (2N/L)-1 if OVERLAP=true), the more stable the resulting power and cross spectra estimates.

Optionally, these power and cross spectra estimates may then be smoothed again in the frequency domain by modified Daniell filters (e.g. if argument SMOOTH_PARAM is used).

The computed equivalent number of degrees of freedom and bandwidth must be divided by two for the zero and Nyquist frequencies.

Furthermore, the computed equivalent number of degrees of freedom, bandwidth, lower and upper (1-PROBTEST) * 100% confidence limit factors and critical value for the squared coherency (e.g. arguments EDOF, BANDWIDTH, CONLWR, CONUPR and TESTCOHER) are not right near the zero and Nyquist frequencies if the PSD estimates have been smoothed by modified Daniell filters. The reason is that CROSS_SPECTRUM2 assumes that smoothing involves averaging independent frequency ordinates. This is true except near the zero and Nyquist frequencies where an average may contain contributions from negative frequencies, which are identical to and hence not independent of positive frequency spectral values. Thus, the number of degrees of freedom in PSD estimates near the 0 and Nyquist frequencies are as little as half the number of degrees of freedom of the spectral estimates away from these frequency extremes if the optional argument SMOOTH_PARAM is used.

If the optional argument SMOOTH_PARAM is used, the computed equivalent number of degrees of freedom, bandwidth, lower and upper (1-PROBTEST) * 100% confidence limit factors and critical value for the squared coherency are right for PSD estimates at frequencies

(i-1)/(L+L0) for i= (nparam+1)/2 + 1 to ( (L+L0) - nparam + 1)/2

where nparam = 2 * (2+sum(SMOOTH_PARAM(:)))- 1, (e.g. for frequencies i/(L+L0) for i = (nparam+1)/2, … , ((L+L0)-nparam-1)/2 ) .

For definitions, more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York.

Welch, P.D., 1967:

The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE trans. on audio and electroacoustics, Vol. Au-15, 2, 70-73.

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford.

`subroutine cross_spectrum2 ( vec, mat, l, psvec, psmat, phase, coher, freq, edof, bandwidth, conlwr, conupr, testcoher, ampli, co_spect, quad_spect, prob_coher, initfft, overlap, normpsd, smooth_param, trend, trend2, win, taperp, l0, probtest )`¶

Purpose¶

Subroutine CROSS_SPECTRUM2 computes Fast Fourier Transform (FFT) estimates of the power and cross spectra of the real time series, VEC, and the multi-channel real time series MAT.

The Power Spectral Density (PSD) and Cross Spectral Density (CSD) estimates are returned in units which are the square of the data (if NORMPSD=false) or in spectral density units (if NORMPSD=true).

Arguments¶

VEC (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the real time series for which the power and cross spectra must be estimated. If TREND=1, 2 or 3, VEC is used as workspace and is transformed.

Size(VEC) must be greater or equal to 4.

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the multi-channel real time series for which the power and cross spectra must be estimated. Each row of MAT is a real time series. If TREND=1, 2 or 3, MAT is used as workspace and is transformed.

The shape of MAT must verify: size(MAT,2) = size(VEC).

L (INPUT) integer(i4b)

On entry, an integer used to segment the time series. L is the length of the segments. L must be a positive even integer, less or equal to size(VEC), but greater or equal to 4.

Spectral computations are at (L/2)+1 frequencies if the optional argument L0 is absent and are at ((L+L0)/2)+1 frequencies if L0 is present (L0 is the number of zeros added to each segment).

Suggested values for L+L0 are 16, 32, 64 or 128 (e.g. an integer power of two, in order to speed the computations).

PSVEC (OUTPUT) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the Power Spectral Density (PSD) estimates of VEC.

PSVEC must verify: size(PSVEC) = ((L+L0)/2) + 1 .

PSMAT (OUTPUT) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and ((L+L0)/2) + 1 columns containing the Power Spectral Density (PSD) estimates of each row of MAT.

The shape of PSMAT must verify:

size(PSMAT,1) = size(MAT,1) ;

size(PSMAT,2) = ((L+L0)/2) + 1 .

PHASE (OUTPUT) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and ((L+L0)/2) + 1 columns containing the phase of the cross spectrum, given in fractions of a circle (e.g. on the closed interval (0,1) ).

The shape of PHASE must verify:

size(PHASE,1) = size(MAT,1) ;

size(PHASE,2) = ((L+L0)/2) + 1 .

COHER (OUTPUT) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and ((L+L0)/2) + 1 columns containing the squared coherency estimates for all frequencies.

The shape of COHER must verify:

size(COHER,1) = size(MAT,1) ;

size(COHER,2) = ((L+L0)/2) + 1 .

FREQ (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a real vector of length ((L+L0)/2)+1 containing the frequencies at which the spectral quantities are calculated in cycles per unit of time.

The spectral estimates are taken at frequencies (i-1)/(L+L0) for i=1,2, … , ((L+L0)/2 + 1).

FREQ must verify: size(FREQ) = (L+L0)/2 + 1 .

EDOF (OUTPUT, OPTIONAL) real(stnd)

On exit, the equivalent number of degrees of freedom of the power and cross spectrum estimates.

BANDWIDTH (OUTPUT, OPTIONAL) real(stnd)

On exit, the bandwidth of the power and cross spectrum estimates.

CONLWR (OUTPUT, OPTIONAL) real(stnd)

CONUPR (OUTPUT, OPTIONAL) real(stnd)

On output, these arguments specify the lower and upper (1-PROBTEST) * 100% confidence limit factors, respectively. Multiply the PSD estimates (e.g. the PSVEC(:) and PSMAT(:,:) arguments ) by these constants to get the lower and upper limits of a (1-PROBTEST) * 100% confidence interval for the PSD estimates.

TESTCOHER (OUTPUT, OPTIONAL) real(stnd)

On output, this argument specifies the critical value for testing the null hypothesis that the squared coherency is zero at the PROBTEST * 100% significance level (e.g. elements of COHER(:,:) less than TESTCOHER should be regarded as not significantly different from zero at the PROBTEST * 100% significance level).

AMPLI (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and ((L+L0)/2) + 1 columns containing the cross-amplitude spectra.

The shape of AMPLI must verify:

size(AMPLI,1) = size(MAT,1) ;

size(AMPLI,2) = ((L+L0)/2) + 1 .

CO_SPECT (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and ((L+L0)/2) + 1 columns containing the co-spectra (e.g. the real part of cross-spectra).

The shape of CO_SPECT must verify:

size(CO_SPECT,1) = size(MAT,1) ;

size(CO_SPECT,2) = ((L+L0)/2) + 1 .

QUAD_SPECT (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and ((L+L0)/2) + 1 columns containing the quadrature spectrum (e.g. the imaginary part of cross-spectrum with a minus sign).

The shape of QUAD_SPECT must verify:

size(QUAD_SPECT,1) = size(MAT,1) ;

size(QUAD_SPECT,2) = ((L+L0)/2) + 1 .

PROB_COHER (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, a real matrix with size(MAT,1) rows and ((L+L0)/2) + 1 columns containing the probabilities that the computed sample squared coherencies came from an ergodic stationary bivariate process with (corresponding) squared coherencies equal to zero.

The shape of PROB_COHER must verify:

size(PROB_COHER,1) = size(MAT,1) ;

size(PROB_COHER,2) = ((L+L0)/2) + 1 .

INITFFT (INPUT, OPTIONAL) logical(lgl)

On entry, if:

INITFFT = false, it is assumed that a call to subroutine INIT_FFT has been done before calling subroutine CROSS_SPECTRUM2 in order to sets up constants and functions for use by subroutine FFT which is called inside subroutine CROSS_SPECTRUM2. This call to INITFFT must have the following form:

call init_fft( (/ size(MAT,1), (L+L0)/2 /), dim=2_i4b )

INITFFT = true, the call to INIT_FFT is done inside subroutine CROSS_SPECTRUM2 and a call to END_FFT is also done before leaving subroutine CROSS_SPECTRUM2.

The default is INITFFT=true .

OVERLAP (INPUT, OPTIONAL) logical(lgl)

If:

OVERLAP = false, the subroutine segments the data without any overlapping.

OVERLAP = true, the subroutine overlaps the segments by one half of their length (which is equal to L).

In both cases, zeros are eventually added to each segment (if argument L0 is present) and each segment will be FFT’d, and the resulting periodograms will averaged together to obtain a Power Spectrum Density estimate at the ((L+L0)/2)+1 frequencies.

The default is OVERLAP=false .

NORMPSD (INPUT, OPTIONAL) logical(lgl)

On entry, if:

NORMPSD = true, the power and cross spectra estimates are normalized in such a way that the total area under the power spectra is equal to the variance of the time series contained in VEC and in each row of MAT.

NORMPSD = false, the sum of the PSD estimates (e.g. sum(PSVEC(2:)) and sum(PSMAT(:,2:),dim=2) ) is equal to the variance of the corresponding time series.

The default is NORMPSD=true .

SMOOTH_PARAM (INPUT, OPTIONAL) integer(i4b), dimension(:)

If SMOOTH_PARAM is used, the power and cross spectra estimates are computed by repeated smoothing of the periodograms and cross-periodogram with modified Daniell weights.

On entry, SMOOTH_PARAM(:) gives the array of the half-lengths of the modified Daniell filters to be applied.

All the values in SMOOTH_PARAM(:) must be greater than 0 and less than ((L+L0)/2)+1 .

TREND (INPUT, OPTIONAL) integer(i4b)

If:

TREND=+1 The means of the time series are removed before computing the spectra

TREND=+2 The drifts from time series are removed before computing the spectra

TREND=+3 The least-squares lines from time series are removed before computing the spectra.

For other values of TREND nothing is done before estimating the power and cross spectra.

The default is TREND=1, e.g. the means of the time series are removed before the computations.

TREND2 (INPUT, OPTIONAL) integer(i4b)

If:

TREND2=+1 The mean of the time segment is removed before computing the cross-spectrum on this segment.

TREND2=+2 The drift from the time segment is removed before computing the cross-spectrum on this segment.

TREND2=+3 The least-squares line from the time segment is removed before computing the cross-spectrum on this segment.

For other values of TREND2 nothing is done before estimating the cross-spectrum on each segment.

The default is TREND2=0, e.g. nothing is done before estimating the power spectrum on each segment.

WIN (INPUT, OPTIONAL) integer(i4b)

On entry, this argument specify the data window used in the computations of the power and cross spectra. If:

WIN=+1 The Bartlett window is used

WIN=+2 The square window is used

WIN=+3 The Welch window is used

WIN=+4 The Hann window is used

WIN=+5 The Hamming window is used

WIN=+6 A split-cosine-bell window is used

The default is WIN=3, e.g. the Welch window is used.

TAPERP (INPUT, OPTIONAL) real(stnd)

The total percentage of the data to be tapered if WIN=6. TAPERP must be greater than zero and less or equal to one, otherwise the default value is used.

The default is 0.2 .

L0 (INPUT, OPTIONAL) integer(i4b)

The number of zeros added to each time segment in order to obtain more finely spaced spectral estimates. L+L0 must be a positive even integer.

The default is L0=0, e.g. no zeros are added to each time segment.

PROBTEST (INPUT, OPTIONAL) real(stnd)

On entry, a probability. PROBTEST is the critical probability which is used to determine the lower and upper confidence limit factors (e.g. the optional arguments CONLWR and CONUPR ) and the critical value for testing the null hypothesis that the squared coherency is zero (e.g. the TESTCOHER optional argument).

PROBTEST must verify: 0. < P < 1.

The default is 0.05 .

Further Details¶

After removing the mean or the trend from the time series (e.g. TREND=1,2,3), the series are padded with zero on the right such that the length of the resulting time series is evenly divisible by L (a positive even integer). The length, N, of these resulting time series is the first integer greater than or equal to size(VEC) which is evenly divisible by L. If size(VEC) is not evenly divisible by L, N is equal to size(VEC)+L-mod(size(VEC),L).

Optionally, the mean or the trend may also be removed from each time segment (e.g. TREND2=1,2,3). Optionally, zeros may be added to each time segment (e.g. the optional arguemnt L0) if more finely spaced spectral esimates are desired.

The stability of the power and cross spectra estimates depends on the averaging process. That is, the greater the number of segments ( N/L if OVERLAP=false and (2N/L)-1 if OVERLAP=true), the more stable the resulting power and cross spectra estimates.

Optionally, these power and cross spectra estimates may then be smoothed again in the frequency domain by modified Daniell filters (e.g. if argument SMOOTH_PARAM is used).

The computed equivalent number of degrees of freedom and bandwidth must be divided by two for the zero and Nyquist frequencies.

Furthermore, the computed equivalent number of degrees of freedom, bandwidth, lower and upper (1-PROBTEST) * 100% confidence limit factors and critical value for the squared coherency (e.g. arguments EDOF, BANDWIDTH, CONLWR, CONUPR and TESTCOHER) are not right near the zero and Nyquist frequencies if the PSD estimates have been smoothed by modified Daniell filters. The reason is that CROSS_SPECTRUM2 assumes that smoothing involves averaging independent frequency ordinates. This is true except near the zero and Nyquist frequencies where an average may contain contributions from negative frequencies, which are identical to and hence not independent of positive frequency spectral values. Thus, the number of degrees of freedom in PSD estimates near the 0 and Nyquist frequencies are as little as half the number of degrees of freedom of the spectral estimates away from these frequency extremes if the optional argument SMOOTH_PARAM is used.

If the optional argument SMOOTH_PARAM is used, the computed equivalent number of degrees of freedom, bandwidth, lower and upper (1-PROBTEST) * 100% confidence limit factors and critical value for the squared coherency are right for PSD estimates at frequencies

(i-1)/(L+L0) for i= (nparam+1)/2 + 1 to ( (L+L0) - nparam + 1)/2

where nparam = 2 * (2+sum(SMOOTH_PARAM(:)))- 1, (e.g. for frequencies i/(L+L0) for i = (nparam+1)/2, … , ((L+L0)-nparam-1)/2 ) .

For definitions, more details and algorithm, see:

Bloomfield, P., 1976:

Fourier analysis of time series- An introduction. John Wiley and Sons, New York.

Welch, P.D., 1967:

The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE trans. on audio and electroacoustics, Vol. Au-15, 2, 70-73.

Diggle, P.J., 1990:

Time series: a biostatistical introduction. Clarendon Press, Oxford.

Module_Time_Series_Procedures¶

subroutine comp_smooth ( x, smooth_factor )¶

Purpose¶

Arguments¶

Further Details¶

subroutine comp_smooth ( x, smooth_factor, dimvar )¶

Purpose¶

Arguments¶

Further Details¶

subroutine comp_smooth ( x, smooth_factor )¶

Purpose¶

Arguments¶

Further Details¶

subroutine comp_trend ( y, nt, itdeg, robust, trend, ntjump, maxiter, rw, no, ok )¶

Purpose¶

Arguments¶

Further Details¶

subroutine comp_trend ( y, nt, itdeg, robust, trend, ntjump, maxiter, rw, no, ok )¶

Purpose¶

Arguments¶

Further Details¶

subroutine comp_stlez ( y, np, ns, isdeg, itdeg, robust, season, trend, ni, nt, nl, ildeg, nsjump, ntjump, nljump, maxiter, rw, no, ok )¶

Purpose¶

Arguments¶

Further Details¶

subroutine comp_stlez ( y, np, ns, isdeg, itdeg, robust, season, trend, ni, nt, nl, ildeg, nsjump, ntjump, nljump, maxiter, rw, no, ok )¶

Purpose¶

Arguments¶

Further Details¶

subroutine comp_stl ( y, np, ni, no, isdeg, itdeg, ildeg, nsjump, ntjump, nljump, ns, nt, nl, rw, season, trend )¶

Purpose¶

Arguments¶

Further Details¶

subroutine comp_stl ( y, np, ni, no, isdeg, itdeg, ildeg, nsjump, ntjump, nljump, ns, nt, nl, rw, season, trend )¶

Purpose¶

Arguments¶

Further Details¶

subroutine ma ( x, len, ave)¶

Purpose¶

Arguments¶

Further Details¶

subroutine detrend ( vec, trend, orig, slope )¶

Purpose¶

Arguments¶

Further Details¶

subroutine detrend ( mat, trend, orig, slope )¶

Purpose¶

Arguments¶

Further Details¶

subroutine hwfilter ( vec, pl, ph, initfft, trend, win )¶

Purpose¶

Arguments¶

Further Details¶

subroutine hwfilter ( mat, pl, ph, initfft, trend, win, max_alloc )¶

Purpose¶

Arguments¶

Further Details¶

subroutine hwfilter2 ( vec, pl, ph, trend, win )¶

Purpose¶

Arguments¶

Further Details¶

subroutine hwfilter2 ( mat, pl, ph, trend, win )¶

Purpose¶

Arguments¶

Further Details¶

function lp_coef ( pl, k, fc, notest_fc )¶

Purpose¶

Arguments¶

Further Details¶

function lp_coef2 ( pl, k, fc, win, notest_fc )¶

Purpose¶

Arguments¶

Further Details¶

function hp_coef ( ph, k, fc, notest_fc )¶

Purpose¶

Arguments¶

Further Details¶

function hp_coef2 ( ph, k, fc, win, notest_fc )¶

Purpose¶

Arguments¶

`subroutine comp_smooth ( x, smooth_factor )`¶

`subroutine comp_smooth ( x, smooth_factor, dimvar )`¶

`subroutine comp_smooth ( x, smooth_factor )`¶

`subroutine comp_trend ( y, nt, itdeg, robust, trend, ntjump, maxiter, rw, no, ok )`¶

`subroutine comp_trend ( y, nt, itdeg, robust, trend, ntjump, maxiter, rw, no, ok )`¶

`subroutine comp_stlez ( y, np, ns, isdeg, itdeg, robust, season, trend, ni, nt, nl, ildeg, nsjump, ntjump, nljump, maxiter, rw, no, ok )`¶

`subroutine comp_stlez ( y, np, ns, isdeg, itdeg, robust, season, trend, ni, nt, nl, ildeg, nsjump, ntjump, nljump, maxiter, rw, no, ok )`¶

`subroutine comp_stl ( y, np, ni, no, isdeg, itdeg, ildeg, nsjump, ntjump, nljump, ns, nt, nl, rw, season, trend )`¶

`subroutine comp_stl ( y, np, ni, no, isdeg, itdeg, ildeg, nsjump, ntjump, nljump, ns, nt, nl, rw, season, trend )`¶

`subroutine ma ( x, len, ave)`¶

`subroutine detrend ( vec, trend, orig, slope )`¶

`subroutine detrend ( mat, trend, orig, slope )`¶

`subroutine hwfilter ( vec, pl, ph, initfft, trend, win )`¶

`subroutine hwfilter ( mat, pl, ph, initfft, trend, win, max_alloc )`¶

`subroutine hwfilter2 ( vec, pl, ph, trend, win )`¶

`subroutine hwfilter2 ( mat, pl, ph, trend, win )`¶

`function lp_coef ( pl, k, fc, notest_fc )`¶

`function lp_coef2 ( pl, k, fc, win, notest_fc )`¶

`function hp_coef ( ph, k, fc, notest_fc )`¶

`function hp_coef2 ( ph, k, fc, win, notest_fc )`¶

`function bd_coef ( pl, ph, k, fch, fcl, notest_fc )`¶

`function bd_coef2 ( pl, ph, k, fch, fcl, win, notest_fc )`¶

`function pk_coef ( freq, k, notest_freq )`¶

`function moddan_coef ( k, smooth_param )`¶

`subroutine freq_func ( nfreq, coef, freqr, four_freq, freq )`¶

`subroutine symlin_filter ( vec, coef, trend, nfilt )`¶

`subroutine symlin_filter ( mat, coef, trend, nfilt )`¶

`subroutine symlin_filter2 ( vec, coef, trend, usefft, initfft )`¶

`subroutine symlin_filter2 ( mat, coef, trend, usefft, initfft )`¶

`subroutine dan_filter ( vec, nsmooth, sym, trend )`¶

`subroutine dan_filter ( mat, nsmooth, sym, trend )`¶

`subroutine moddan_filter ( vec, smooth_param, sym, trend )`¶

`subroutine moddan_filter ( mat, smooth_param, sym, trend )`¶

`function extend ( vec, index, sym )`¶

`function extend ( mat, index, sym )`¶

`subroutine taper ( vec, taperp )`¶