`FourierAnalysis.Planner`

β TypeFFTW plans to be used by all functions in **Fourier Analysis** are incapsulated in this structure. This is the only *operator object* created by this package, all others being *data objects*.

Spectral and cross-spectral computations (frequency domain objects) need only a forward plan (real FFT), while objects based on the analytic signal (time-frequency domain objects) need both a forward (`.p`

) and a backward (`.ip`

) plan (complex iFFT).

Argument `flags`

must be one of the constants here above. Basic usage of the flags involves a trade-off between the time needed to compute the planner and the efficacy for computing the FFTs. The following flags are sorted in ascending order of time needed to be computed and efficacy: `plan_estimate`

, `plan_measure`

, `plan_patient`

, `plan_exhaustive`

. By default *FourierAnalysis* adopts `plan_estimate`

, which allows the quickest search, but the least optimal choice. Other flags may be worth if several FFT computations are to be done with the same setting. See here for details.

Argument `timelimit`

is the maximum time (in seconds) allowed to FFTW for optimizing the plans. Setting this to `-1.0`

amounts to imposing no time limits.

FFTW plans are computed for a given window length `wl`

and data type `type`

.

`Planners`

objects may be passed as an argument to constructors of FDobjects and TFobjects by using one of the `Planner`

constuctor here below.

**Constructors**

```
Planner(flags :: UInt32,
timelimit :: Union{Int, Float64},
wl :: Int,
type :: Type,
bw :: Bool = false)
```

Use this to create a Planner object, passing as argument a FFTW `flags`

constant and `timelimit`

, the window length `wl`

and the type of the input data `type`

, which must be real, e.g., Float64.

If `bw`

is false (default), a dummy backward plan is created, otherwise a backward plan is created for the complex data type corresponding to `type`

(e.g., if `type`

is Float64, it will created for ComplexF64 data.)

For example, suppose π is a vector of many matrices of multivariate time-series sampled at 128 samples per second and that we want to compute the spectra for all of them using a 256-point FFT. We first create a plan by

`p=Planner(plan_exhaustive, 10.0, 256, eltype(π[1]))`

Then we invoke the Spectra function passing the plan as argument:

`π=spectra(π, sr, t; planner=p)`

A shorter construction syntax is available when only the forward plan is needed and the type of the data is Float64:

```
Planner(flags :: UInt32,
timelimit :: Union{Int, Float64},
wl :: Int)
```

For example, the line above could have been written more shortly as

`p=Planner(plan_exhaustive, 10.0, 256)`

In order to create also a backward plan you would use instead

`p=Planner(plan_exhaustive, 10.0, 256, eltype(π[1]), true)`

`DSP.Periodograms.coherence`

β Method```
(1)
function coherence(π :: CrossSpectra;
allkinds :: Bool = false)
(2)
function coherence(π’ :: CrossSpectraVector;
allkinds :: Bool = false,
check :: Bool = true)
(3)
function coherence(X :: Matrix{T},
sr :: Int,
wl :: Int;
tapering :: Union{Taper, TaperKind} = harris4,
planner :: Planner = getplanner,
slide :: Int = 0,
DC :: Bool = false,
nonlinear :: Bool = false,
smoothing :: Smoother = noSmoother,
tril :: Bool = false,
allkinds :: Bool = false,
β© :: Bool = true) where T<:Real
(4)
function coherence(π :: Vector{Matrix{T}},
< same argument sr, ..., β© of method (3) > where T<:Real
```

**alias**: coh

This identifier may conflict with the DSP package. Invoke it as `FourierAnalysis.coherence`

.

(1)

Construct a Coherence object from a CrossSpectra object, allowing coherence estimates using the Welch method. All non-data fields are copied from the cross-spectra object, i.e., all fields but `y`

, which holds the coherence matrices that are computed by this function.

If `π.tril`

is true, only the lower triangular part of the coherence matrices is computed and the matrices in `y`

are LowerTriangular matrices, otherwise the full matrices are computed and `y`

will hold Hermitian matrices (see LinearAlgebra).

If optional keyword argument `allkinds`

is true, all five kinds of coherence are returned. In this case the output is a 5-tuple of Coherence objects, in the order:

*total*coherence,*real*coherence,*instantaneous*coherence*imaginary*coherence,*lagged*coherence.

If `allkinds`

is false (default), only the *total* (classical) coherence is returned as a single `Coherence`

object.

(2)

Construct a CoherenceVector object from the input CrossSpectraVector object `π’`

, allowing coherence estimates using the Welch method. This method calls method (1) for all objects in `π’`

.

The `allkinds`

optional keyword parameter has the same meaning as in method (1). In this method if the argument is passed as true, the output is a 5-tuple of `CoherenceVector`

objects.

If optional keyword argument `check`

is true (default), it is checked that all `CrossSpectra`

objects in `π’`

have the same value in the ``.sr`

, `.wl`

, `.DC`

, `.taper`

, `.nonlinear`

and `.smoothing`

field. If this is not the case, an error message is printed pointing to the first field that is not identical in all objects and `Nothing`

is returned.

(3)

Given a multivariate time series data matrix `X`

of dimension $t$ x $n$, where $t$ is the number of samples (rows) and $n$ the number of time-series (columns), sampling rate `sr`

and epoch length `wl`

, compute the squared coherence matrices of `X`

, that is, the coherence matrices at all Fourier discrete frequencies obtained from the Welch (sliding window) average cross-spectra.

**optioal keyword arguments**:

`sr`

, `wl`

, `tapering`

, `planner`

and `slide`

have the same meaning as for the `spectra`

function.

`DC`

, `nonlinear`

, `smoothing`

, `tril`

and `β©`

have the same meaning as for the `crossSpectra`

function, to which they apply for estimating the cross-spectra.

The `allkinds`

optional keyword parameter has the same meaning as in method (1).

(4)

Construct a CoherenceVector object from a vector of real multivariate data matrices. Compute the coherence matrices from the cross-spectral matrices estimated using the Welch method as per method (3) for all $k$ data matrices in `π`

.

The $k$ matrices in `π`

must have the same number of columns (*i.e.*, the same number of time-series), but may have any number of (at least `wl`

) rows (samples). All other arguments have the same meaning as in method (3), with the following difference:

`β©`

: if true (default), the method is run in multi-threaded mode across the $k$ data matrices if $k$ is at least twice the number of threads Julia is instructed to use, otherwise this method attempts to run each coherence estimation in multi-threaded mode across series as per method (3). See Threads.

If a `Planner`

is not explicitly passed as an argument, the FFTW plan is computed once and applied for all coherence estimations.

**note for methods (1) and (3)**:

If `tril`

is false (default), the output coherence object(s) is of type `Array{Hermitian,1}`

, which is the `βVector`

type used in package PosDefManifold. Since coherence estimates are symmetric positive definite, they can be straightaway used as argument to PosDefManifold's functions, e.g., for computing matrix moves on geodesics, matrix distances, etc. and the the whole vector output to compute matrix means, spectral embedding and more.

If `tril`

is true, the output is of type `Array{LowerTriangular,1}`

, which is the `πVector`

type used in PosDefManifold, that is, only the lower triangle of the coherencee is computed in order to save time and memory.

**note for methods (2) and (4)**:

If `tril`

is false, the output is of type `Array{Array{Hermitian,1},1}`

, which is the `βVectorβ`

type used in PosDefManifold.

If `tril`

is true, the output is of type `Array{Array{LowerTriangular,1},1}`

, which is the `πVectorβ`

type used in PosDefManifold.

**See**: crossspectra.jl, Spectra, Coherence.

**Examples**:

```
## common code for methods (1)-(4)
using FourierAnalysis, LinearAlgebra
function generateSomeData(sr::Int, t::Int; noise::Real=1.)
# four sinusoids of length t samples and sr sampling rate
# peak amplitude: 0.7, 0.6, 0.5, 0.4
# frequency: 5, 7, 13, 27
# phase: 0, Ο/4, Ο/2, Ο
v1=sinusoidal(0.7, 5, sr, t, 0)
v2=sinusoidal(0.6, 7, sr, t, Ο/4)
v3=sinusoidal(0.5, 13, sr, t, Ο/2)
v4=sinusoidal(0.4, 27, sr, t, Ο)
return hcat(v1, v2, v3, v4) + (randn(t, 4)*noise)
end
sr, wl, t = 128, 512, 8192
# (1)
X=generateSomeData(sr, t) # multivariate data matrix 8192x4
# cross-spectra using default harris4 tapering window
S=crossSpectra(X, sr, wl)
# Only classical coherence
C=FourierAnalysis.coherence(S)
# All 5 kinds of coherence
Ctot, C2real, C3inst, C4imag, C5lag=FourierAnalysis.coherence(S, allkinds=true);
Ctot
# (2)
# generate 3 multivariate data matrices 8192x4
X=[generateSomeData(sr, t) for i=1:3]
# cross-spectra using default harris4 tapering window
S=crossSpectra(X, sr, wl)
# Only classical coherence
C=FourierAnalysis.coherence(S)
# All 5 kinds of coherence
Ctot, C2real, C3inst, C4imag, C5lag=FourierAnalysis.coherence(S, allkinds=true);
Ctot
# (3)
X=generateSomeData(sr, t) # multivariate data matrix 8192x4
# coherence using default harris4 tapering window
C=FourierAnalysis.coherence(X, sr, wl)
# check the coherence matrix at frequency 5Hz
C.y[f2b(5, sr, wl)]
# coherence using hann tapering window
C=FourierAnalysis.coherence(X, sr, wl; tapering=hann)
# using Slepian's multi-tapering
C=FourierAnalysis.coherence(X, sr, wl; tapering=slepians(sr, wl))
# compute non-linear coherence (phase-locking value)
C=FourierAnalysis.coherence(X, sr, wl; tapering=slepians(sr, wl), nonlinear=true)
# compute only the lower triangle of coherence matrices
C=FourierAnalysis.coherence(X, sr, wl; tapering=slepians(sr, wl), tril=true)
# compute all five kinds of coherence
Ctot, Creal, Cinst, Cimag, Clag=FourierAnalysis.coherence(X, sr, wl;
tapering=slepians(sr, wl), tril=true, allkinds=true);
Ctot
# smooth a-posteriori the coherence
C2=smooth(blackmanSmoother, C)
# or compute coherence already smoothed
C=FourierAnalysis.coherence(X, sr, wl;
tapering=slepians(sr, wl), tril=true, smoothing=blackmanSmoother)
# mean coherence matrix in 8Hz-12Hz range
M=mean(C, (8, 12)) # or also M=mean(C, (8.0, 12.0))
# extract all coherence matrices in 8Hz-12Hz range
E=extract(C, (8, 12))
# coherence matrices averaged in 2Hz band-pass regions
B=bands(C, 2)
# (4)
# generate 3 multivariate data matrices 8192x4
X=[generateSomeData(sr, t) for i=1:3]
# The examples here below use exactly the same syntax as method (3).
# However, since the input here is a vector of data matrices
# and not a single data matrix, the examples here below create a vector
# of the object created by the examples of method (3).
# For example:
# coherence using the default harris4 tapering window
# this creates a CoherenceVector object
C=FourierAnalysis.coherence(X, sr, wl)
# check the first Coherence object
C[1]
# check the coherence matrix at fr. 5Hz for the first Coherence object
C[1].y[f2b(5, sr, wl)]
# coherence using Hann tapering window
C=FourierAnalysis.coherence(X, sr, wl; tapering=hann)
# using Slepian's multi-tapering
C=coherence(X, sr, wl; tapering=slepians(sr, wl))
# compute non-linear coherence
C=FourierAnalysis.coherence(X, sr, wl; tapering=slepians(sr, wl), nonlinear=true)
# compute only the lower triangle of coherence matrices
C=FourierAnalysis.coherence(X, sr, wl; tapering=slepians(sr, wl), tril=true)
# compute all five kinds of coherence
Ctot, Creal, Cinst, Cimag, Clag=FourierAnalysis.coherence(X, sr, wl;
tapering=slepians(sr, wl), tril=true, allkinds=true);
Ctot
# smooth a-posteriori all coherence objects in S
C2=smooth(blackmanSmoother, C)
# or compute them all already smoothed
C=FourierAnalysis.coherence(X, sr, wl; tapering=parzen, smoothing=blackmanSmoother)
# mean coherence matrix in 8Hz-12Hz range for all coherence objects
M=mean(C, (8, 12)) # or also M=mean(C, (8.0, 12.0))
# grand-average mean of the above across all Coherence objects
meanM=mean(mean(C, (8, 12)))
# extract all coherence matrices in 8Hz-12Hz range for all coherence objects
E=extract(C, (8, 12))
# grand average of coherence matrices in 8Hz-12Hz range for all coh. objects
meanE=mean(extract(C, (8, 12)))
# coherence matrices averaged in 2Hz band-pass regions for all coh. objects
B=bands(C, 2)
# Pre-compute a FFT planner and pass it as argument
# (this interesting if the function is to be called repeatedly).
plan=Planner(plan_exhaustive, 10.0, wl, eltype(X[1])) # wait 10s
C=FourierAnalysis.coherence(X, sr, wl; planner=plan)
# how faster is this?
using BenchmarkTools
@benchmark(FourierAnalysis.coherence(X, sr, wl))
@benchmark(FourierAnalysis.coherence(X, sr, wl; planner=plan))
...
...
```

`DSP.Periodograms.coherence`

β Method```
(5)
function coherence(πβ :: TFAnalyticSignalVector,
πβ :: TFAnalyticSignalVector,
frange :: fInterval,
trange :: tInterval;
nonlinear :: Bool = false,
allkinds :: Bool = false,
mode :: Function = extract,
func :: Function = identity,
w :: Vector = [],
check :: Bool = true)
(6)
function coherence(π±β :: Vector{Vector{T}},
π±β :: Vector{Vector{T}},
sr :: Int,
wl :: Int,
frange :: fInterval,
trange :: tInterval,
bandwidth :: IntOrReal = 2;
nonlinear :: Bool = false,
allkinds :: Bool = false,
mode :: Function = extract,
func :: Function = identity,
w :: Vector = [],
filtkind :: FilterDesign = Butterworth(2),
fmin :: IntOrReal = bandwidth,
fmax :: IntOrReal = srΓ·2,
fsmoothing :: Smoother = noSmoother,
tsmoothing :: Smoother = noSmoother,
planner :: Planner = getplanner,
β© :: Bool = true) where T<:Real
```

**alias**: coh

If optional keyword parameter `nonlinear`

is false (default), estimate linear (weighted) coherence measure, otherwise estimate the the (weighted) phase concentration measure.

If optional keyword argument `allkinds`

is true all five kinds of coherence are returned. In this case the output is a 5-tuple of `Coherence`

matrices, in the order:

*total*coherence,*real*coherence,*instantaneous*coherence*imaginary*coherence,*lagged*coherence.

If `allkinds`

is false (default) only the *total* (classical) coherence is returned as a single `Coherence`

matrix.

(5) The desired measure is obtained averaging across the TFAnalyticSignal objects in `π`

. Since this method uses pre-computed analytic signal objects, their `.nonlinear`

field must agree with the `nonlinear`

argument passed to this function.

`frange`

, `trange`

, `w`

, `mode`

and `func`

have the same meaning as in the `meanAmplitude`

function, however keep in mind that the two possible `mode`

functions, i.e., `extract`

and `mean`

, in this function operate on complex numbers.

The checks performed in the `meanAmplitude`

function are performed here too. In addition, if `check`

is true, also check that

- if
`nonlinear`

is true, all objects in`π`

are nonlinear; - if
`nonlinear`

is false, all objects in`π`

are linear.

If either check fails, print an error message and return `Nothing`

.

(6) Estimate the analytic signal of all data vectors in `π±`

calling the `TFanalyticsignal`

constructor and then use method (5) to obtained the desired measure.

`frange`

, `trange`

, `mode`

, `func`

, `w`

and `check`

have the same meaning as in the `meanAmplitude`

function. The other arguments are passed to the `TFanalyticsignal`

constructor, to which the reader is referred for understanding their action.

**See also**: `meanAmplitude`

, `meanDirection`

, timefrequencybi.jl.

**Examples**: see the examples of `comodulation`

function.

`FourierAnalysis.TFamplitude`

β Method```
(1)
function TFamplitude(Z::TFAnalyticSignal;
func::Function=identity)
(2)
function TFamplitude(π::TFAnalyticSignalVector;
func::Function=identity)
(3)
function TFamplitude(x :: Vector{T},
sr :: Int,
wl :: Int,
bandwidth :: IntOrReal = 2;
func :: Function = identity,
filtkind :: FilterDesign = Butterworth(2),
fmin :: IntOrReal = bandwidth,
fmax :: IntOrReal = srΓ·2,
fsmoothing :: Smoother = noSmoother,
tsmoothing :: Smoother = noSmoother,
planner :: Planner = getplanner,
β© :: Bool = true) where T<:Real =
(4)
function TFamplitude(π± :: Vector{Vector{T}},
< same arguments as method (3) >
```

(1)

Construct a TFAmplitude object computing the amplitude of TFAnalyticSignal object `Z`

. Optional keyword argument `func`

is a function to be applied element-wise to the data matrix of the output. By default, the `identity`

(do nothing) function is applied.

(2)

Construct a TFAmplitudeVector object from a TFAnalyticSignalVector object executing method (1) for all TFAnalyticSignal objects in `π`

(3)

Call `TFanalyticsignal`

to obtain the time-frequency analytic signal of real signal vector `x`

and construct a TFAmplitude object holding the time-frequency amplitude (the mudulus, often referred to as the *envelope*) of `x`

.

All arguments are used for regulating the estimation of the analytic signal, with the exception of `func`

, `fsmoothing`

and `fsmoothing`

:

`func`

is an optional function to be applied to the amplitude data matrix output.

In order to estimate the analytic signal in the time-frequency domain this function calls the `TFanalyticsignal`

constructor (method (1) therein), with both `fsmoothing`

and `tsmoothing`

arguments set to `noSmoother`

. Arguments `fsmoothing`

and `fsmoothing`

are then used to smooth the amplitude.

In order to obtain amplitude estimations on smoothed analytic signal instead, create a TFAnalyticSignal object passing a Smoother to the `TFanalyticsignal`

constructor and then use method (1) to obtain the amplitude. Such amplitude estimation can be further smoothed using the `smooth`

function, as shown in the examples.

For the meaning of all other arguments, which are passed to function `TFanalyticsignal`

, see the documentation therein.

(4)

Construct a TFAmplitudeVector object from a vector of real signal vectors `π±`

, executing method (3) for all of them. In order to estimate the time-frequency analytic signal for a vector of signals, method (2) of `TFanalyticsignal`

is called.

**See**: `TFanalyticsignal`

, TFAmplitude.

**See also**: `amplitude`

.

**Examples**: see the examples of `TFanalyticsignal`

.

`FourierAnalysis.TFanalyticsignal`

β Method```
(1)
function TFanalyticsignal(x :: Vector{T},
sr :: Int,
wl :: Int = 0,
bandwidth :: IntOrReal = 2;
fmin :: IntOrReal = bandwidth,
fmax :: IntOrReal = srΓ·2,
filtkind :: FilterDesign = Butterworth(2),
nonlinear :: Bool = false,
fsmoothing :: Smoother = noSmoother,
tsmoothing :: Smoother = noSmoother,
planner :: Planner = getplanner,
β© :: Bool = true) where T<:Real
(2)
function TFanalyticsignal(π± :: Vector{Vector{T}},
<same arguments as method (1)>
```

(1)

Given sampling rate `sr`

, construct a TFAnalyticSignal object from univariate data `x`

, that is, a time-frequency representation of real vector `x`

.

Call three functions in series performing the three following operations:

i. pass the data throught a bank of pass-band filters calling function `filterbank`

,

ii. compute the analytic signal calling function `analyticsignal`

(method (1) therein),

iii. If requested, smooth the analytic signal along the time and/or frequency dimension calling function `smooth`

.

The arguments passed to each of these functions are:

filterbank | analyticsignal | smooth |
---|---|---|

`x` | `wl` | `fsmoothing` |

`sr` | `nonlinear` | `tsmoothing` |

`bandwidth` | `planner` | |

`filtkind` | `β©` | |

`fmin` | ||

`fmax` | ||

`β©` |

Refer to the documentation of each function to learn the meaning of each argument.

By default the `wl`

argument is set to `length(x)`

.

If `β©`

is true (default), filtering and computation fo the analytic signal are multi-threaded across band-pass regions as long as the number of the regions is at least twice the number of threads Julia is instructed to use. See Threads.

(2)

Given sampling rate `sr`

, construct a TFAnalyticSignalVector object from a vector of univariate data `π±`

, that is, from the time-frequency representations of the vectors in `π±`

.

This method operates as method (1) with the following exceptions:

By default the `wl`

argument is set to `length(x)`

for each vectors in `π±`

. If another values is given, it will be used for all of them.

If `β©`

is true (default), the method is run in multi-threaded mode across the vectors in `π±`

as long as their number is at least twice the number of threads Julia is instructed to use, otherwise this method attempts to run each analytic signal estimation in multi-threaded mode like in method (1). See Threads.

If a `Planner`

is not explicitly passed as an argument, the FFTW plan is computed once and applied for all analytic signal estimations.

**See**: `filterbank`

, `analyticsignal`

, `smooth`

.

**Examples**:

```
using Plots, FourierAnalysis
# generate some data
sr, t, bandwidth=128, 512, 2
h=taper(harris4, t)
x1=sinusoidal(10, 8, sr, t, 0)
x2=sinusoidal(10, 19, sr, t, 0)
x=Vector((x1+x2).*h.y+randn(t))
y1=sinusoidal(10, 6, sr, t, 0)
y2=sinusoidal(10, 16, sr, t, 0)
y=Vector((y1+y2).*h.y+randn(t))
plot([x, y])
# TFAnalyticSignal object for x (method (1))
Y = TFanalyticsignal(x, sr, sr*4; fmax=32)
# vector of TFAnalyticSignal objects for x and y (method (2))
π = TFanalyticsignal([x, y], sr, sr*4; fmax=32)
# (for shortening comments: TF=time-frequency, AS=analytic signal)
# mean AS in a TF region (8:12Hz and samples 1:128) for the first object
m=mean(π[1], (8, 12), (1, 128)) # Output a complex number
# extract the AS in a TF region (8:12Hz and samples 1:128) for the first object
E=extract(π[1], (8, 12), (1, 128)) # Output a complex matrix
# mean AS in a TF region (8:12Hz and samples 1:128) for the two objects
π¦=mean(π, (8, 12), (1, 128)) # Output a vector of complex numbers
# same computation without checking homogeneity of the two objects in π (faster)
π¦=mean(π, (8, 12), (1, 128); check=false)
# extract the AS in a TF region (8:12Hz and samples 1:128) for the two objects
π=extract(π, (8, 12), (8, 12)) # Output a vector of complex matrices
# plot the real part of the AS of x (see unit recipes.jl)
# gather first useful attributes for the plot
using Plots.Measures
tfArgs=(right_margin = 2mm,
top_margin = 2mm,
xtickfont = font(10, "Times"),
ytickfont = font(10, "Times"))
heatmap(tfAxes(Y)..., real(Y.y);
c=:fire, tfArgs...)
# ...the imaginary part
heatmap(tfAxes(Y)..., imag(Y.y);
c=:bluesreds, tfArgs...)
# ...the amplitude
heatmap(tfAxes(Y)..., amplitude(Y.y);
c=:amp, tfArgs...)
# ...the amplitude of the AS smoothed in frequency and time
heatmap(tfAxes(Y)..., amplitude(smooth(hannSmoother, hannSmoother, Y).y);
c=:fire, tfArgs...)
# ...the phase
heatmap(tfAxes(Y)..., phase(Y.y);
c=:bluesreds, tfArgs...)
# or generate a TFAmplitude object from the AS
A=TFamplitude(Y)
# and plot it (with different colors)
heatmap(tfAxes(A)..., A.y;
c=:fire, tfArgs...)
# generate a TFPhase object
Ο΄=TFphase(Y)
# and plot it (with custom colors)
heatmap(tfAxes(Ο΄)..., Ο΄.y;
c=:fire, tfArgs...)
# compute and plot phase in [0, 2Ο]
heatmap(tfAxes(Y)..., TFphase(Y; func=x->x+Ο).y;
c=:amp, tfArgs...)
# compute and plot unwrapped phase
heatmap(tfAxes(Y)..., TFphase(Y; unwrapped=true).y;
c=:amp, tfArgs...)
# smooth time-frequency AS: smooth frequency
Z=smooth(blackmanSmoother, noSmoother, Y)
# plot amplitude of smoothed analytic signal
heatmap(tfAxes(Z)..., amplitude(Z.y);
c=:amp, tfArgs...)
# not equivalently (!), create an amplitude object and smooth it:
# in this case the amplitude is smoothed, not the AS
A=smooth(blackmanSmoother, noSmoother, TFamplitude(Y))
heatmap(tfAxes(A)..., A.y;
c=:fire, tfArgs...)
# Smoothing raw phase estimates is unappropriate
# since the phase is a discontinous function, however it makes sense
# to smooth phase if the phase is unwrapped.
heatmap(tfAxes(Y)...,
smooth(blackmanSmoother, noSmoother, TFphase(Y; unwrapped=true)).y;
c=:amp, tfArgs...)
# smooth AS: smooth both frequency and time
E=smooth(blackmanSmoother, blackmanSmoother, Y)
# plot amplitude of smoothed analytic signal
heatmap(tfAxes(E)..., amplitude(E.y);
c=:fire, tfArgs...)
# plot phase of smoothed analytic signal
heatmap(tfAxes(E)..., phase(E.y);
c=:bluesreds, tfArgs...)
# not equivalently (!), create amplitude and phase objects and smooth them
A=smooth(blackmanSmoother, blackmanSmoother, TFamplitude(Y))
heatmap(tfAxes(A)..., A.y;
c=:fire, tfArgs...)
Ο΄=smooth(blackmanSmoother, blackmanSmoother, TFphase(Y, unwrapped=true))
heatmap(tfAxes(Ο΄)..., Ο΄.y;
c=:fire, tfArgs...)
# smooth again
Ο΄=smooth(blackmanSmoother, blackmanSmoother, Ο΄)
heatmap(tfAxes(Ο΄)..., Ο΄.y;
c=:fire, tfArgs...)
# and again ...
# create directly smoothed AS
Y=TFanalyticsignal(x, sr, t, bandwidth;
fmax=32,
fsmoothing=hannSmoother,
tsmoothing=hannSmoother)
# plot amplitude of smoothed analytic signal
heatmap(tfAxes(Y)..., amplitude(Y.y);
c=:amp, tfArgs...)
# create directly smoothed Amplitude
A=TFamplitude(x, sr, t, bandwidth;
fmax=32,
fsmoothing=hannSmoother,
tsmoothing=hannSmoother)
# plot smoothed amplitude
heatmap(tfAxes(A)..., A.y;
c=:amp, tfArgs...)
# compute a TFAnalyticSignal object with non-linear AS
Y=TFanalyticsignal(x, sr, t, bandwidth; fmax=32, nonlinear=true)
# check that it is non-linear
Y.nonlinear
# check that the amplitude is now 1.0 everywhere
norm(amplitude(Y.y)-ones(eltype(Y.y), size(Y.y))) # must be zero
# plot non-linear phase
heatmap(tfAxes(Y)..., phase(Y.y);
c=:bkr, tfArgs...)
# get the center frequencies of TFAmplitude object A
A.flabels
# extract the amplitude in a time-frequency region
extract(A, (2, 10), (1, 256)) # output a real matrix
# extract the amplitude in a time-frequency region at only one frequency
extract(A, 10, (1, 256)) # output a row vector
# extract the amplitude at one temporal sample at one frequency
extract(A, 10, 12) # or extract(A, 10.0, 12)
# extract amplitude at one temporal sample in a frequency region
extract(A, (10, 12), 12) # or extract(A, (10.0, 12.0), 12)
# extract amplitude at one temporal sample and all frequencies
extract(A, :, 12) # output a (column) vector
# compute the mean in a time-frequency region:
mean(A, (2, 10), (1, 256)) # output a real number
# is equivalent to (but may be less efficient than)
mean(extract(A, (2, 10), (1, 256)))
# using column sign for extracting all time samples
extract(A, (2, 10), :)
# This:
extract(A, :, :)
# is equivalent to this:
A.y
# but if you don't need to extract all frequencies,
# use the extract function to control what frequencies will be extracted:
# This
extract(A, (4, 5), 10)
# is not equivalent to this
A.y[4:5, 10]
# since the `extract` function finds the correct rows corresponding
# to the sought frequencies (in Hz), while A.y[4:5, 10]
# just returns the elements [4:5, 10] in the TF amplitude object
# Although the first center frequency in A is 2Hz, its
# band-pass region is 1-3Hz, therefore the frequency range 1:10 is accepted
mean(A, (1, 10), (1, 256))
# but this result in an error (see the REPL) since 0.5 is out of range:
mean(A, (0.5, 10), (1, 256))
# using a colon sign for time range
a=mean(A, (1, 10), :)
# using an integer for time range (indicates one specific sample)
a=mean(A, (1, 10), 16)
# using a colon sign for frequency range
a=mean(A, :, (1, 16))
# using a real number for frequency range
a=mean(A, 8.5, (1, 16))
# NB: the `extract` and `mean` functions work with the same syntax
# for objects TFAnayticSignal, TFAmplitude and TFPhase.
```

`FourierAnalysis.TFphase`

β Method```
(1)
function TFphase(Z :: TFAnalyticSignal;
func :: Function = identity,
unwrapped :: Bool = false)
(2)
function TFphase(π :: TFAnalyticSignalVector;
func ::Function = identity,
unwrapped ::Bool = false)
(3)
function TFphase(x :: Vector{T},
sr :: Int,
wl :: Int,
bandwidth :: IntOrReal = 2;
unwrapped :: Bool = false,
func :: Function = identity,
filtkind :: FilterDesign = Butterworth(2),
fmin :: IntOrReal = bandwidth,
fmax :: IntOrReal = srΓ·2,
nonlinear :: Bool = false,
fsmoothing :: Smoother = noSmoother,
tsmoothing :: Smoother = noSmoother,
planner :: Planner = getplanner,
β© :: Bool = true) where T<:Real
(4)
function TFphase(π± :: Vector{Vector{T}},
< same arguments as method (3) >
```

(1)

Construct a TFPhase object computing the phase of TFAnalyticSignal object `Z`

. By default the phase is represented in $[βΟ, Ο]$.

If optional keyword argument `unwrapped`

is true (false by defalut), the phase is unwrapped, that is, it holds the cumulative sum of the phase along the time dimension once this is represented in $[0, 2Ο]$.

Optional keyword argument `func`

is a function to be applied element-wise to the data matrix of the output. By default, the `identity`

(do nothing) function is applied. If `unwrapped`

is true, the function is applied on the unwrapped phase.

(2)

Construct a TFPhaseVector object from a TFAnalyticSignalVector object executing method (1) for all TFAnalyticSignal objects in `π`

(3)

Call `TFanalyticsignal`

to obtain the time-frequency analytic signal of real signal vector `x`

and construct a TFPhase object holding the time-frequency phase (argument) of `x`

.

All arguments are used for regulating the estimation of the analytic signal, with the exception of `unwrapped`

, `func`

, `fsmoothing`

and `fsmoothing`

.

`unwrapped`

has the same meaning as in method (1) and (2).

`func`

is an optional function to be applied to the phase data matrix output. If `unwrapped`

is true, it is applied to the unwrapped phase.

In order to estimate the analytic signal in the time-frequency domain this function calls the `TFanalyticsignal`

constructor (method (1) therein), with both `fsmoothing`

and `tsmoothing`

arguments set to `noSmoother`

. `fsmoothing`

and `fsmoothing`

are then used to smooth the phase if `unwrapped`

is true.

In order to obtain phase estimations on smoothed analytic signal instead, create a TFAnalyticSignal object passing a Smoother to the `TFanalyticsignal`

constructor and then use method (1) to obtain the phase.

For the meaning of all other arguments, which are passed to function `TFanalyticsignal`

, see the documentation therein.

(4)

Construct a TFPhaseVector object from a vector of real signal vectors `π±`

, executing method (3) for all of them. In order to estimate the time-frequency analytic signal for a vector of signals, method (2) of `TFanalyticsignal`

is called.

**See**: `TFanalyticsignal`

, TFPhase.

**See also**: `phase`

, `unwrapPhase`

, `polar`

.

**Examples**: see the examples of `TFanalyticsignal`

.

`FourierAnalysis.amplitude`

β Method```
(1)
function amplitude(c::Complex;
func::Function=identity) = func(abs(c))
(2)
function amplitude(A::AbstractArray{T};
func::Function=identity) where T<:Complex
(3)
function amplitude(A::TFAnalyticSignal;
func::Function=identity)
(4)
function amplitude(π::TFAnalyticSignalVector;
func::Function=identity)
```

(1)

Return the amplitude (modulus) of a complex number. This corresponds to Julia's abs function. It is here provided for syntactic consistency with the following methods.

(2)

Return the amplitude of a complex array `Z`

. Typically, `Z`

holds analytic signal, in which case the output is the analytic (instantaneous) amplitude (also known as envelope). The output is a real array of the same size as `Z`

.

(3)

Return a real matrix with the analytic (instantaneous) amplitude of the TFAnalyticSignal object `Z`

. The output is of the same size as the data field `Z.y`

.

(4)

As (3), but return a vector of amplitude matrices for all TFAnalyticSignal objects in π

~

In all methods if a function is provided by the optional keyword argument `func`

, it is applied element-wise to the output. For example,

- passing
`func=x->x^2`

will return the power, - passing
`func=x->log(x^2)`

will return the log-power, - passing
`func=x->decibel(x^2)`

will return the power in deciBels.

**See**: TFAnalyticSignal.

**Examples**:

```
using FourierAnalysis, Plots
x=sinusoidal(10, 2, 128, t*4, 0).*sinusoidal(10, 1, 128, t*4, 0)
# amplitude and phase of a vector using analytic signal standard method
y=analyticsignal(x)
a=amplitude(y)
Ο=phase(y, func=x->(x+Ο)/2Ο*50)
plot([x, a, Ο]; labels=["signal", "amplitude", "phase"])
# see what happen if `x` contains energy in frequencies below sr/wl Hz
# (see documentation of `analyticSignal` function)
y=analyticsignal(x, 64)
a=amplitude(y)
Ο=phase(y, func=x->(x+Ο)/2Ο*50)
plot([x, a, Ο]; labels=["signal", "amplitude", "phase"])
# unwrapped phase
# the line below will do nothing as argument `unwrapdims` is 0 by default
Ο2=unwrapPhase(phase(y))
# this will do the job
Ο2=unwrapPhase(phase(y); unwrapdims=1)
plot([x, a, Ο2./25]; labels=["signal", "amplitude", "unwr. phase"])
# amplitude from analytic signal of a data matrix holding multiple series
X=randn(t, 4)
Y=analyticsignal(X)
A=amplitude(Y)
plot(A[:, 1:2])
# phase
π·=phase(Y)
plot(π·[:, 1:1])
# unwrapped phase
π·2=unwrapPhase(π·; unwrapdims=1)
plot(π·2)
# phase represented in [-1, 1]
π·=phase(Y, func=x->(x+Ο)/2Ο)
plot(π·[:, 1:1])
# sine of the phase
π·=phase(Y, func=sin)
plot(π·[:, 1:1])
# get Amplitude and phase from analytic Signal
A, π·=polar(Y)
A
π·
```

`FourierAnalysis.analyticsignal`

β Method```
(1)
function analyticsignal( X :: Union{Vector{T}, Matrix{T}},
wl :: Int = size(X, 1);
nonlinear :: Bool = false,
planner :: Planner = getplanner,
β© :: Bool = true) where T<:Real
(2)
function analyticsignal( π :: Vector{Matrix{T}},
wl :: Int;
nonlinear :: Bool = false,
planner :: Planner = getplanner,
β© :: Bool = true) where T<:Real
```

(1)

Compute the analytic signal(AS) of vector `X`

or of all column vectors of matrix `X`

via the FFT and iFFT procedure, as explained in Marple(1999). If `wl`

=size(X, 1) (default), use the standard method passing to the FFT and iFFT all samples in `X`

altogether, whereas if `wl`

<size(X, 1) a sliding-windows method is used (see below).

Return the analytic signal `π`

, a complex vector if `X`

is a vector or a complex matrix holding in its columns the analytic signal of the columns of `X`

if `X`

is a matrix. `π`

has the same number of samples (rows) as `X`

.

Contrarely to what is done in the DSP package, the DC level of the signal is removed, therefore, if the input signal features a non-zero DC level, the real part of the AS will be equal to the input signal with the DC level removed. The imaginary part of `π`

is the Hilbert transform of such no-DC `X`

.

The sliding-windows AS allows an efficient estimation of the AS for vectors and matrices when they are very large; it proceeds computing the AS on 50% sliding overlapping windows and forming the AS by retaining the central half of each window. The number of points effectively used to obtain the final estimation is $wl$Γ·2 (integer division). `wl`

must be even for using this estimation. This procedure produces edge effects, thus the first and last $wlΓ·4$ samples of the AS estimation should be discarded. In practice, one requires the AS of a larger data segment and trims at least $wlΓ·4$ samples at the beginning and end of the estimation. This is done automatically by the `TFanalyticsignal`

function.

Also, the sliding-windows AS method creates small discontinuities at sample $wl$Γ·4 and then every $wl$Γ·2 samples, therefore $wl$ should be chosen as large as possible.

In order to avoid FFT computation of very long epochs, if `wl`

> 2^14, then `wl`

is set to 2^10. Below this limit, as long as the computations are feasable, use the standard method. If you absolutely need to use the sliding-windows method, see window length in FFTW for setting efficiently argument `wl`

.

The input signal should be previously band-pass or high-pass filtered so as not to contain frequency components below the first discrete Fourier frequency obtained with windows of `wl`

samples, that is, below sr/wl Hz, where sr is the sampling rate.

**Optional Keyword Arguments**:

`nonlinear`

, if true, the analytic signal is normalized so that its amplitude is $1.0$ at all points. This allow non-linear univariate and bivariate estimations (see timefrequencyuni.jl and timefrequencybi.jl).

`planner`

is an instance of the `Planner`

object, holding the forward and backward FFTW plans used to compute the FFTs and the iFFTs. By default the planner is computed, but it can be passed as an argumet here if it is pre-computed. This is interesting if the `analyticsignal`

function is to be invoked repeatedly.

if `β©`

is true, the method is run in multi-threaded mode across the series in `X`

if the number of series is at least twice the number of threads Julia is instructed to use. See Threads.

(2)

Compute the analytic signal for all $k$ multivariate data matrices given as a vector of matrices `π`

. Return a vector of matrices hodling the corresponding analytic signals as in method (1). The FFT and iFFT plans are computed only once. The $k$ matrices in `π`

may have different number of columns (i.e., different number of series) and different number of rows (samples). However, the number of rows must be larger than `wl`

for all of them.

If `β©`

is true, this method run in multi-threaded mode across the matrices in `π`

if the number of matrices is at least twice the number of threads Julia is instructed to use, otherwise it tries to run each analytic signal estimation in multi-threaded mode as per method (1). See Threads.

This function is called by the following functions operating on time-frequency reprsentations: `TFanalyticsignal`

, `TFamplitude`

, `TFphase`

, `meanAmplitude`

, `concentration`

, `meanDirection`

, `comodulation`

, `coherence`

.

**References** Marple S.L. (1999) Computing the Discrete-Time Analytic Signal via FFT. IEEE Transactions on Signal Processing 47(9), 2600-3.

**Examples**:

```
using FourierAnalysis, FFTW, LinearAlgebra, Statistics, Plots, DSP
t=128; lab=["x", "real(y)", "imag(y)"]
# Analytic signal of one vector
x=sinusoidal(10, 2, 128, t, Ο/2; DC=10) # cosine
y=analyticsignal(x)
# note that real(y) is x without the DC level, i.e., x=real(y)+DC
plot([x, real(y), imag(y)]; labels=lab)
# make a check
s=sinusoidal(10, 2, 128, t, 0) # sine
norm(s-imag(y)) # should be zero
# Analytic Signal by DSP.jl
y2=hilbert(x)
norm(s-imag(y2)) # should be zero
# DSP.jl does not remove the DC level
# thus x=real(y2) in this case
plot([x, real(y2), imag(y2)]; labels=lab)
# Analytic signal of multiple vectors
x=hcat(x, sinusoidal(10, 3, 128, t, Ο/2; DC=10))
y=analyticsignal(x)
# sliding-windows analytic signal of one vector
# (note edge effects)
x=sinusoidal(10, 2, 128, t*4, Ο/2; DC=0)
y=analyticsignal(x, t)
plot([x, real(y), imag(y)]; labels=lab)
# sliding-windows analytic signal of multiple vectors
x=hcat(x, sinusoidal(10, 3, 128, t*4, Ο/2; DC=0))
y=analyticsignal(x, t)
```

`FourierAnalysis.b2f`

β Method```
function b2f(bin :: Int,
sr :: Int,
wl :: Int;
DC :: Bool = false)
```

**b**in **to f**requency. Return the closest discrete Fourier frequency (in Hz) that corresponds to a bin (position) in a real-FFT vector, given sampling rate `sr`

and window length `wl`

. The FFT vector is assumed to be 1-based, as always in Julia.

If `DC`

is false, the first discrete frequency is assumed to be at bin 1, otherwise the DC level is assumed to be at bin 1 and the first discrete frequency at bin 2.

**See also**: `f2b`

, `fres`

, `fdf`

, `brange`

.

**Examples**:

```
using FourierAnalysis
f2b(20, 512, 1024) # return 40
f2b(10, 128, 128) # return 10
```

`FourierAnalysis.bands`

β Method```
function bands(S :: Union{FDobjects, FDobjectsVector}
bandwidth :: IntOrReal)
```

Return band-pass average of spectral, cross-spectral or coherence estimates in equally spaced band-pass regions with the given `bandwidth`

. `bandwidth`

can be given as an integer or as a real number. See `bbands`

for details on the definition of band-pass regions.

Band-pass average is not supported for time-frequency objects as for those objects a similar averaging is natively avaiable using argument `bandwidth`

in their constructors.

The output of this function is as it follows:

- for univariate Spectra objects (i.e., those hodling one spectrum only), a real column vector,
- for multivariate Spectra objects, a real matrix,
- for SpectraVector objects, a vector of the above,
- for CrossSpectra and Coherence objects, a vector of Hermitian or LowerTriangular matrices, depending on how the object has been cosntructed,
- for CrossSpectraVector and CoherenceVector objects, a vector of the above.

**See**: `bbands`

.

**Examples**:

```
using FourierAnalysis, Plots
# example with univariate Spectra objects (one series -> one spectrum)
sr, t, f, a = 128, 256, 10, 1
# create a sinusoidal superimposed to white noise
v=sinusoidal(a, f, sr, t*16, 0) + randn(t*16)
# compute the spectrum
Ξ£=spectra(v, sr, t)
# mean spectra in 2Hz-band-pass regions
b=bands(Ξ£, 2)
plot(b)
# example with multivariate spectra (several series -> several spectra)
Ξ£=spectra(hcat(v, v+randn(t*16), v+randn(t*16) ), sr, t)
# mean spectra in 2Hz-band-pass regions for all time-series
b=bands(Ξ£, 2)
plot(b)
# plot mean spectra in 2Hz-band-pass regions for time-series 2 and 3 only
plot(bands(Ξ£, 2)[:, 2:3])
# example with CrossSpectra objects (the same goes for Coherence objects)
X=broadcast(+, v, randn(t*16, 3))*randn(3, 3)
Ξ£=crossSpectra(X, sr, t)
# mean cross-spectra in 4Hz-band-pass regions
B=bands(Ξ£, 4)
# example with multiple CrossSpectra objects
X2=broadcast(+, v, randn(t*16, 3))*randn(3, 3)
Ξ£=crossSpectra([X, X2], sr, t) # this is a CrossSpectraVector
# mean cross-spectra in 4Hz-band-pass regions for all cross-spectra objects
B=bands(Ξ£, 4)
```

`FourierAnalysis.bbands`

β Method```
function bbands(sr :: Int,
wl :: Int,
bandwidth :: IntOrReal;
DC :: Bool = false)
```

Return a vector of integers holding the limits of all `bandwidth`

-spaced band-pass regions of a real-FFT, in bins of discrete Fourier frequencies, from one to $wlΓ·2$ (integer division).

This is used by function `bands`

.

To know the frequencies in Hz to which these bins correspond, call `fbands`

.

**See**: `bands`

.

**See also**: `fbands`

.

**Examples**:

```
using FourierAnalysis
bbands(128, 256, 16) # return [1, 32, 64, 96, 128]
fbands(128, 256, 16) # return [0.5, 16.0, 32.0, 48.0, 64.0]
bbands(128, 256, 16; DC=true) # return [2, 33, 65, 97, 129]
fbands(128, 256, 16; DC=true) # return [0.5, 16.0, 32.0, 48.0, 64.0]
bbands(128, 128, 16) # return [1, 16, 32, 48, 64]
fbands(128, 128, 16) # return [1.0, 16.0, 32.0, 48.0, 64.0]
```

`FourierAnalysis.brange`

β Method```
function brange(wl :: Int;
DC :: Bool = false)
```

Return a range of bins for a real-FFT vector covering all Fourier discrete frequencies given window length `wl`

.

If `DC`

is false, the range is $1:(wlΓ·2)$ (integer division), otherwise it is $1:(wlΓ·2)+1$.

**See also**: `f2b`

, `fres`

, `b2f`

, `fdf`

.

**Examples**:

```
using FourierAnalysis
brange(0.5, 8) # return 1:4
```

`FourierAnalysis.comodulation`

β Method```
(1)
function comodulation( πβ :: TFAmplitudeVector,
πβ :: TFAmplitudeVector,
frange :: fInterval,
trange :: tInterval;
mode :: Function = extract,
func :: Function = identity,
w :: Vector = [],
check :: Bool = true) =
(2)
function comodulation( πβ :: TFAnalyticSignalVector,
πβ :: TFAnalyticSignalVector,
< arguments frange, trange, mode, func, w and check as in method (1) >
(3)
function comodulation(π±β :: Vector{Vector{T}},
π±β :: Vector{Vector{T}},
sr :: Int,
wl :: Int,
frange :: fInterval,
trange :: tInterval,
bandwidth :: IntOrReal = 2;
mode :: Function = extract,
func :: Function = identity,
w :: Vector = [],
filtkind :: FilterDesign = Butterworth(2),
fmin :: IntOrReal = bandwidth,
fmax :: IntOrReal = srΓ·2,
fsmoothing :: Smoother = noSmoother,
tsmoothing :: Smoother = noSmoother,
planner :: Planner = getplanner,
β© :: Bool = true) where T<:Real
```

**alias**: com

(1)

Given a pair of TFAmplitudeVector objects, estimate their (weighted) comodulation measure. `πβ`

and `πβ`

must hold the same number of objects and the time-frequency planes of all the objects should be congruent.

**arguments**:

`frange`

and `trange`

have the same meaning as in the `meanAmplitude`

method.

**optional keyword arguments**

`mode`

, `func`

and `check`

have the same meaning as in the `meanAmplitude`

method.

`w`

may be a vector of non-negative real weights associated to each pair of input objects. By default the unweighted version of the measure is computed.

(2)

Given a pair of TFAnalyticSignalVector object, compute the amplitude of all objects and estimate the (weighted) comodulation as per method (1). Like in method (1), `πβ`

and `πβ`

must hold the same number of objects and the time-frequency planes of all the objects should be congruent. In addition, since using this method all TFAnalyticSignal objects in `πβ`

and `πβ`

must be `linear`

, if `check`

is true (default) and this is not the case, print an error and return `Nothing`

. The checks on `frange`

and `trange`

performed by method (1) are also performed by this method.

(3)

Estimate the amplitude of all data vectors in `π±β`

and `π±β`

calling the `TFamplitude`

constructor and then estimate the (weighted) comodulation measure across the constructed amplitude objects as per method (1).

`frange`

, `trange`

, `mode`

, `func`

, `w`

and `check`

have the same meaning as in method (1). The other arguments are passed to the `TFamplitude`

constructors, to which the reader is referred for their meaning.

If a `planner`

for FFT computations is not provided, it is computed once and applied for all amplitude estimations.

**See also**: `coherence`

, timefrequencybi.jl.

**Examples**:

```
using FourierAnalysis
# generate 100 pairs of data vectors
sr, t, bandwidth=128, 512, 2
h=taper(harris4, t)
π±β=[sinusoidal(2, 10, sr, t, 0).*h.y+randn(t) for i=1:100]
π±β=[sinusoidal(2, 10, sr, t, 0).*h.y+randn(t) for i=1:100]
# compute their (linear) analytic signal
πβ=TFanalyticsignal(π±β, sr, wl, bandwidth; fmax=32, nonlinear=false)
πβ=TFanalyticsignal(π±β, sr, wl, bandwidth; fmax=32, nonlinear=false)
# compute their amplitude
πβ=TFamplitude(πβ)
πβ=TFamplitude(πβ)
# compute the Com averaging in a TF region from TFAnalyticSignal objects
# πβ and πβ must be linear
Com=comodulation(πβ, πβ, (8, 12), :; mode=mean)
# compute the Com averaging in a TF region from TFAmplitudeVector objects
# πβ and πβ must be linear
Com=comodulation(πβ, πβ, (8, 12), :; mode=mean)
# compute the Com averaging in a TF region directly from data
# In this case you don't have to worry about linearity
Com=comodulation(π±β, π±β, sr, wl, (8, 12), :, bandwidth; mode=mean)
# compute comodulation from smoothed amplitude:
Com=comodulation(π±β, π±β, sr, wl, (8, 12), :, bandwidth;
mode=mean,
fsmoothing=blackmanSmoother,
tsmoothing=blackmanSmoother)
# you can go faster pre-computing a FFTW plan.
# This is useful when you have to call the comodulation function several times
plan=Planner(plan_patient, 5, wl, Float64, true)
Com=comodulation(π±β, π±β, sr, wl, (8, 12), :, bandwidth; mode=mean, planner=plan)
# compute the Com in a TF region from TFAnalyticSignalVector objects
Com=comodulation(πβ, πβ, (8, 12), :; mode=extract)
# compute the Com in a TF region from TFAmplitudeVector objects
Com=comodulation(πβ, πβ, (8, 12), :; mode=extract)
# compute the Com in a TF region directly from data
Com=comodulation(π±β, π±β, sr, wl, (8, 12), :, bandwidth; mode=extract)
# All these operations can be done also for coherence measures, for example
Coh=coherence(πβ, πβ, (8, 12), :; mode=mean)
Coh=coherence(πβ, πβ, (8, 12), :; mode=extract)
# Compute all 5 coherence types
Coh=coherence(πβ, πβ, (8, 12), :; mode=extract, allkinds=true)
# phase coherence (phase-locking value)
# we obtain this measure from non-linear TFAnalyticSignalVector objects
πβ=TFanalyticsignal(π±β, sr, wl, bandwidth; fmax=32, nonlinear=true)
πβ=TFanalyticsignal(π±β, sr, wl, bandwidth; fmax=32, nonlinear=true)
Coh=coherence(πβ, πβ, (8, 12), :; mode=mean, nonlinear=true)
# or directly from data (no need to worry about non-linearity in this case)
Coh=coherence(π±β, π±β, sr, wl, (8, 12), :, bandwidth; mode=mean, nonlinear=true)
```

`FourierAnalysis.concentration`

β Method```
(1)
function concentration( π :: TFAnalyticSignalVector,
frange :: fInterval,
trange :: tInterval;
nonlinear :: Bool = false,
mode :: Function = extract,
func :: Function = identity,
w :: Vector = [],
check :: Bool = true)
(2)
function concentration( π± :: Vector{Vector{T}},
sr :: Int,
wl :: Int,
frange :: fInterval,
trange :: tInterval,
bandwidth :: IntOrReal = 2;
nonlinear :: Bool = false,
mode :: Function = extract,
func :: Function = identity,
w :: Vector = [],
check :: Bool = true,
filtkind :: FilterDesign = Butterworth(2),
fmin :: IntOrReal = bandwidth,
fmax :: IntOrReal = srΓ·2,
fsmoothing :: Smoother = noSmoother,
tsmoothing :: Smoother = noSmoother,
planner :: Planner = getplanner,
β© :: Bool = true) where T<:Real
```

**alias**: con

If optional keyword parameter `nonlinear`

is false (default), estimate the (weighted) concentration measure, otherwise estimate the (weighted) phase concentration measure.

(1) The desired measure is obtained averaging across the TFAnalyticSignal objects in `π`

. Since this method uses pre-computed analytic signal objects, their `.nonlinear`

field must agree with the `nonlinear`

argument passed to this function.

`frange`

, `trange`

, `w`

, `mode`

and `func`

have the same meaning as in the `meanAmplitude`

function, however keep in mind that the two possible `mode`

functions, i.e., `extract`

and `mean`

, in this function operate on complex numbers.

The checks performed in the `meanAmplitude`

function are performed here too. In addition, if `check`

is true, also check that

- if
`nonlinear`

is true, all objects in`π`

are nonlinear; - if
`nonlinear`

is false, all objects in`π`

are linear.

If either check fails, print an error message and return `Nothing`

.

(2) Estimate the analytic signal of all data vectors in `π±`

calling the `TFanalyticsignal`

constructor and then use method (1) to obtained the desired measure.

`frange`

, `trange`

, `mode`

, `func`

, `w`

and `check`

have the same meaning as in the `meanAmplitude`

function. The other arguments are passed to the `TFanalyticsignal`

constructor, to which the reader is referred for understanding their action.

**See also**: `meanAmplitude`

, `meanDirection`

, timefrequencybi.jl.

**Examples**: see the examples of `meanAmplitude`

function.

`FourierAnalysis.crossSpectra`

β Method```
(1)
function crossSpectra( X :: Matrix{T},
sr :: Int,
wl :: Int;
tapering :: Union{Taper, TaperKind} = harris4,
planner :: Planner = getplanner,
slide :: Int = 0,
DC :: Bool = false,
nonlinear :: Bool = false,
smoothing :: Smoother = noSmoother,
tril :: Bool = false,
β© :: Bool = true) where T<:Real
(2)
function crossSpectra( π :: Vector{Matrix{T}},
< same argument sr, ..., β© of method (1) > where T<:Real
```

(1)

Construct a CrossSpectra objects from real multivariate data using the Welch method.

Given sampling rate `sr`

and epoch length `wl`

, compute the cross-spectral matrices of dimension $n$x$n$ for a multivariate data matrix `X`

of dimension $t$x$n$, where $t$ is the number of samples (rows) and `n>1`

is the number of time-series (columns).

The cross-spectral matrices are hold in the `.y`

vector field of the created object. The length of `.y`

depends upon the `wl`

argument and `DC`

optional keyword argument (see below).

**Optional Keyword Arguments**:

`sr`

, `wl`

, `tapering`

, `planner`

and `slide`

have the same meaning as for the `spectra`

function.

`DC`

: if true the cross-spectral matrix of the DC level is returned in the first position of `y`

(see the fields of the CrossSpectra object), otherwise (default) the matrices in `y`

start with the first positive discrete frequency, that is, $sr/wl$ Hz.

`nonlinear`

: if true, the amplitude information is eliminated from the DFT coefficients, that is, they are normalized by their modulus before being averaged. This leads to non-linear estimates (Congedo, 2018; Pascual-Marqui 2007) where the diagonal elements of the cross-spectral matrices (the spectra) are 1.0 for all frequencies. By default, it is false.

`smoothing`

: apply a smoothing function of type Smoother to the cross-spectral matrices across frequencies. By default no smoothing is applied.

`tril`

: if false (default), the whole cross-spectra matrices will be computed, otherwise only their lower triangular part (see below).

`β©`

: if true (default), the method is run in multi-threaded mode across the series in `X`

if the number of series is at least twice the number of threads Julia is instructed to use. See Threads.

**Return**

If `tril`

is false (default), the output is of type `Array{Hermitian,1}`

, which is the `βVector`

type used in package PosDefManifold. Since cross-spectral estimates are Hermitian positive definite, they can be straightaway used as argument to PosDefManifold's functions, e.g., for computing matrix moves on geodesics, matrix distances, etc. and the the whole vector output to compute matrix means, spectral embedding and more.

If `tril`

is true, the output is of type `Array{LowerTriangular,1}`

, which is the `πVector`

type used in PosDefManifold, that is, only the lower triangle of the cross-spectra is computed in order to save time and memory.

(2)

Construct a CrossSpectraVector object from a vector of real multivariate data matrices. Compute the cross-spectral matrices using the Welch method as per method (1) for all $k$ data matrices in `π`

.

The $k$ matrices in `π`

must have the same number of columns (i.e., the same number of time-series), but may have any number of (at least `wl`

) rows (samples). All other arguments have the same meaning as in method (1), with the following difference:

`β©`

: if true (default), the method is run in multi-threaded mode across the $k$ data matrices if $k$ is at least twice the number of threads Julia is instructed to use, otherwise this method attempts to run each cross-spectral estimation in multi-threaded mode across series as per method (1). See Threads.

If a `planner`

is not explicitly passed as an argument, the FFTW plan is computed once and applied for all cross-spectral estimations.

**Return**

If `tril`

is false, the output is of type `Array{Array{Hermitian,1},1}`

, which is the `βVectorβ`

type used in PosDefManifold.

If `tril`

is true, the output is of type `Array{Array{LowerTriangular,1},1}`

, which is the `πVectorβ`

type used in PosDefManifold.

**See**: CrossSpectra.

**Examples**:

```
using FourierAnalysis, Plots, LinearAlgebra
function generateSomeData(sr::Int, t::Int; noise::Real=1.)
# four sinusoids of length t samples and sr sampling rate
# peak amplitude: 0.7, 0.6, 0.5, 0.4
# frequency: 5, 7, 13, 27
# phase: 0, Ο/4, Ο/2, Ο
v1=sinusoidal(0.7, 5, sr, t, 0)
v2=sinusoidal(0.6, 7, sr, t, Ο/4)
v3=sinusoidal(0.5, 13, sr, t, Ο/2)
v4=sinusoidal(0.4, 27, sr, t, Ο)
return hcat(v1, v2, v3, v4) + (randn(t, 4)*noise)
end
sr, wl, t = 128, 512, 8192
# (1)
X=generateSomeData(sr, t) # multivariate data matrix 8192x4
# cross-spectra using default harris4 tapering window
S=crossSpectra(X, sr, wl)
# check the cross-spectral matrix at frequency 5Hz
S.y[f2b(5, sr, wl)]
# check only the diagonal part of this matrix as a vector
diag(S.y[f2b(5, sr, wl)])
# cross-spectra using hann tapering window
S=crossSpectra(X, sr, wl; tapering=hann)
# using Slepian's multi-tapering
S=crossSpectra(X, sr, wl; tapering=slepians(sr, wl))
# compute non-linear cross-spectra
S=crossSpectra(X, sr, wl; tapering=slepians(sr, wl), nonlinear=true)
# compute only the lower triangle of cross-spectral matrices
S=crossSpectra(X, sr, wl; tapering=slepians(sr, wl), tril=true)
# smooth a-posteriori the cross-spectra
S2=smooth(blackmanSmoother, S)
# or compute cross-spectra already smoothed
S=crossSpectra(X, sr, wl;
tapering=slepians(sr, wl), tril=true, smoothing=blackmanSmoother)
# mean cross-spectral matrix in 8Hz-12Hz range
M=mean(S, (8, 12)) # or also M=mean(S, (8.0, 12.0))
# extract all cross-spectral matrices in 8Hz-12Hz range
E=extract(S, (8, 12))
# cross-spectral matrices averaged in 2Hz band-pass regions
B=bands(S, 2)
# Get the spectra from a CrossSpectra object
PowerSpectra=Spectra(S)
# Get the amplitude spectra from a CrossSpectra object
AmpSpectra=Spectra(S, func=β)
# Get the log10-spectra from a CrossSpectra object
log10Spectra=Spectra(S, func=log10)
# plot the spectra (see recipes.jl)
plot(AmpSpectra; fmax=32, xspace=4, ytitle="Amplitude")
# (2)
# generate 3 multivariate data matrices 8192x4
X=[generateSomeData(sr, t) for i=1:3]
# The examples here below use exactly the same syntax as the previous method.
# However, since the input here is a vector of data matrices
# and not a single data matrix, the examples here below create a vector
# of the object created by the examples above.
# For example:
# cross-spectra using the default harris4 tapering window
# this creates a CrossSpectraVector object
S=crossSpectra(X, sr, wl)
# check the cross-spectral matrix at fr. 5Hz for the first CrossSpectra object
S[1].y[f2b(5, sr, wl)]
# check only the diagonal part of this matrix as a vector
diag(S[1].y[f2b(5, sr, wl)])
# cross-spectra using Hamming's tapering window
S=crossSpectra(X, sr, wl; tapering=hamming)
# using Slepian's multi-tapering
S=crossSpectra(X, sr, wl; tapering=slepians(sr, wl))
# compute non-linear cross-spectra
S=crossSpectra(X, sr, wl; tapering=slepians(sr, wl), nonlinear=true)
# compute only the lower triangle of cross-spectral matrices
S=crossSpectra(X, sr, wl; tapering=slepians(sr, wl), tril=true)
# smooth a-posteriori all CrossSpectra objects in S
S2=smooth(blackmanSmoother, S)
# or compute them all already smoothed
S=crossSpectra(X, sr, wl; tapering=parzen, smoothing=blackmanSmoother)
# mean cross-spectral matrix in 8Hz-12Hz range for all CrossSpectra (CS) objects
M=mean(S, (8, 12)) # or also M=mean(S, (8.0, 12.0))
# grand-average mean of the above across all CS objects
meanM=mean(mean(S, (8, 12)))
# extract all cross-spectral matrices in 8Hz-12Hz range for all CS objects
E=extract(S, (8, 12))
# grand average of cross-spectral matrices in 8Hz-12Hz range for all CS objects
meanE=mean(extract(S, (8, 12)))
# cross-spectral matrices averaged in 2Hz band-pass regions for all CS objects
B=bands(S, 2)
# Get and plot the spectra from a CrossSpectra object
plot(Spectra(S[1]); fmax=32, xspace=4)
# Pre-compute a FFT planner and pass it as argument
# (this interesting if the function is to be called repeatedly).
plan=Planner(plan_exhaustive, 10.0, wl, eltype(X[1])) # wait 10s
S=crossSpectra(X, sr, wl; planner=plan)
# how faster is this?
using BenchmarkTools
@benchmark(crossSpectra(X, sr, wl))
@benchmark(crossSpectra(X, sr, wl; planner=plan))
...
...
```

`FourierAnalysis.decibel`

β Method```
(1)
function decibel(S :: Union{Real, AbstractArray{T}}) where T<:Real
(2)
function decibel(S1 :: Union{Real, AbstractArray{T}},
S2 :: Union{Real, AbstractArray{T}}) where T<:Real
```

Convert (1) a measure `S`

, or (2) a ratio between two measures `S1`

./`S2`

into deciBels.

Input measures can be real numbers or real arrays of any dimensions.

For array input, the ratio and the conversion is computed element-wise.

**Examples**:

```
using FourierAnalysis
v=sinusoidal(3., 1, 128, 256, 0)
s=spectra(v, 128, 256; func=decibel) # compute the spectra in dB
s.y # show the spectra
decibel(s.y)
decibel(10.0)
N=abs.(randn(3, 3))
decibel(N)
```

`FourierAnalysis.extract`

β Method```
(1)
function extract(S :: FDobjects,
frange :: fInterval)
(2)
function extract(π :: FDobjectsVector,
frange :: fInterval;
w :: Vector = [],
check :: Bool = true)
(3)
function extract(Y :: TFobjects,
frange :: fInterval,
trange :: tInterval)
(4)
function extract(π :: TFobjectsVector,
frange :: fInterval,
trange :: tInterval;
w :: Vector = [],
check :: Bool = true)
```

**alias**: `extr`

Extract data in a frequency region from FDobjects and data in a time-frequency region from TFobjects. The frequency and time region are indicated by `frange`

and `trange`

, which are of type fInterval and tInterval, respectively.

The input/output types of this function for a region with more then one frequency and more than one sample is reported in the following table:

method | input object | output |
---|---|---|

(1.1) | Spectra | a real matrix with spectra in `frange` arranged in columnsΒΉ |

(1.2) | CrossSpectra | a vector of complex matrices holding the cross-spectra in `frange` Β² |

(1.3) | Coherence | a vector of real matrices holding the coherence in `frange` Β² |

(2.1) | SpectraVector | a vector of matrices of type (1.1) |

(2.2) | CrossSpectraVector | a vector of vectors of type (1.2) |

(2.3) | CoherenceVector | a vector of vectors of type (1.3) |

(3.1) | TFAnalyticSignal | a complex matrix holding the analytic signal in [`frange` , `trange` ] |

(3.2) | TFAmplitude | a real matrix holding the amplitude in [`frange` , `trange` ] |

(3.3) | TFPhase | a real matrices holding the phase in [`frange` , `trange` ] |

(4.1) | TFAnalyticSignalVector | a vector of matrices of type (3.1) |

(4.2) | TFAmplitudeVector | a vector of matrices of type (3.2) |

(4.3) | TFPhaseVector | a vector of matrices of type (3.3) |

Legend: ΒΉ *each column refers to a time-series on which the spectra have been computed.* Β² *depending on how the objects has been created, the matrices may be either Hermitian or LowerTriangular. See the documentation of CrossSpectra and Coherence.

Note that depending on the arguments the type of the output may loose one or two dimensions. For instance,

- if the Spectra object holds only one spectrum, (1.1) will output a column vector and (2.1) a vector of column vectors.
- if
`frange`

points to a single frequency, (1.1) will output a row vector and (2.1) a vector of row vectors. - if both the above two conditions hold, (1.1) will output a real number and (2.1) a vector.
- if
`frange`

points to a single frequency, (1.2), (1.3) will output a matrix and (2.2), (2.3) a vector of matrices. - If
`frange`

points to a single frequency band, (3.1), (3.2), (3.3) will output a row vector and (4.1), (4.2), (4.3) a vector of row vectors. - If
`trange`

points to a single time sample, (3.1), (3.2), (3.3) will output a column vector and (4.1), (4.2), (4.3) a vector of column vectors. - if both the above two conditions hold, (3.1), (3.2), (3.3) will output a number and (4.1), (4.2), (4.3) a vector.

Method (2) and (4) allows the following *optional keyword arguments*:

`w`

, a $k$-vector of non-negative integers or real numbers, where $k$ is the numbers of objects hold in the input FDobjectsVector or TFobjectsVector. `w`

is a vector of weights for the regions extracted from the input objects. By default, no weights are assigned.

`check`

, a boolean. If it is true (default), it is checked that the non-data fields of the input objects are all the same (for example, sampling rate, bandwidth, etc.). Set it to false to improve speed.

**See also**: `mean`

.

**Examples**:

```
using FourierAnalysis
# example with univariate Spectra objects (one series -> one spectrum)
sr, t, f, a = 128, 256, 10, 1
# create a sinusoidal superimposed to white noise
v=sinusoidal(a, f, sr, t*16, 0) + randn(t*16)
# compute univariate spectra
Ξ£=spectra(v, sr, t)
# spectra in between 8Hz and 12Hz
s=extract(Ξ£, (8, 12))
# spectra in between 8Hz and 12.5Hz
s=extract(Ξ£, (8, 12.5))
# spectra at 10Hz
s=extract(Ξ£, 10) # or s=extract(S, (10, 10))
# these two expressions are equivalent: s=extract(Ξ£, :), s=Ξ£.y
# example with multivariate spectra (several series -> several spectra)
Ξ£=spectra(hcat(v, v+randn(t*16)), sr, t)
# spectra in between 8Hz and 12Hz
S=extract(Ξ£, (8, 12))
# spectra at 10Hz
S=extract(Ξ£, 10)
# example with CrossSpectra objects (the same goes for Coherence objects)
X=broadcast(+, v, randn(t*16, 3))*randn(3, 3)
Ξ£=crossSpectra(X, sr, t)
# cross-spectra in between 8Hz and 12Hz (Hermitian matrices)
S=extract(Ξ£, (8, 12))
Ξ£=crossSpectra(X, sr, t; tril=true)
# cross-spectra in between 8Hz and 12Hz (LowerTriangular matrices)
S=extract(Ξ£, (8, 12))
# example with multiple cross-spectra
X2=broadcast(+, v, randn(t*16, 3))*randn(3, 3)
Ξ£=crossSpectra([X, X2], sr, t) # this is a CrossSpectraVector
S=extract(Ξ£, (8, 12); w=[0.4, 0.6])
# now S[1] holds the cross-spectra in range 8-12Hz for X
# and S[2] holds the cross-spectra in range 8-12Hz for X2
# example with time-frequency objects
# (work in the same way for TFAnalyticSignal, TFAmplitude and TFPhase)
Y = TFanalyticsignal(v, sr, t)
# analytic signal within frequencies 8Hz and 12Hz and time samples 1 to 64.
AS=extract(Y, (8, 12), (1, 64))
# all analytic signal within frequencies 8Hz and 12Hz.
AS=extract(Y, (8.0, 12), :) # accept integers and reals for frequencies
# all analytic signal within time samples 1 to 64.
AS=extract(Y, :, (1, 64))
# example with multiple time-frequency objects
# (notice how the type of the output changes)
Y = TFanalyticsignal([v, v+randn(t*16)], sr, t)
AS=extract(Y, (8, 12), (1, 64))
AS=extract(Y, (8), :)
AS=extract(Y, 8, 2)
```

`FourierAnalysis.f2b`

β Method```
function f2b(f :: IntOrReal,
sr :: Int,
wl :: Int;
DC :: Bool = false)
```

**f**requency **to b**in. Return the bin (position) in a real-FFT vector best matching a frequency `f`

(in Hz), given sampling rate `sr`

and window length `wl`

. The frequency can be given either as an integer or as a real number.

If the requested `f`

is exactly in between two Fourier discrete frequencies, then the smallest of the two equidistant frequencies is returned.

The FFT vector is assumed to be 1-based (as always in Julia). If `DC`

is false, the first discrete frequency is assumed to be at bin 1, otherwise the DC is assumed to be at bin 1 and the first discrete frequency at bin 2.

If `DC`

is false return 0 for frequencies inferior to half the frequency resolution.

**See also**: `fres`

, `b2f`

, `fdf`

, `brange`

.

**Examples**:

```
using FourierAnalysis
f2b(10, 512, 1024) # return 20
```

`FourierAnalysis.fbands`

β Method`FourierAnalysis.fdf`

β Method```
function fdf(sr :: Int,
wl :: Int;
DC :: Bool = false)
```

Return a vector with all **F**ourier **d**iscrete **f**requencies for a real-FFT, given sampling rate `sr`

and window length `wl`

. If `DC`

is false, the first discrete frequency starts at bin (position) 1 and the length of the vector is $wlΓ·2$ (integer division), otherwise the DC level is at position 1. and the length of the vector is $(wlΓ·2)+1$.

**See also**: `f2b`

, `fres`

, `b2f`

, `brange`

.

**Examples**:

```
using FourierAnalysis
fdf(8, 16)
# return the 8-element Array{Float64,1}:
# [0.5, 1.0, 1.5, 2.0, 2.5, 3, 3.5, 4.0]
```

`FourierAnalysis.filterbank`

β Method```
function filterbank(x :: Vector{T},
sr :: Int,
bandwidth :: IntOrReal = 2;
filtkind :: FilterDesign = Butterworth(2),
fmin :: IntOrReal = bandwidth,
fmax :: IntOrReal = srΓ·2,
β© :: Bool = true) where T<:Real
```

Pass signal vector `x`

throught a bank of band-pass filters, given sampling rate `sr`

and `bandwidth`

parameter.

The kind of filters is specified by optional keyword argument `filtkind`

, of type FilterDesign, using the DSP package. By default `filtkind`

is a forward-backward (linear phase response) Butterworth IIR filter of order 2 in each direction (hence, 4th order total). See notes on DSP package useful functions for tips on how to design other IIR and FIR filters.

Return 2-tuple `(f, Y)`

, where `f`

is a vector holding the center frequencies of the filter bank band-pass regions and `Y`

a matrix holding in the $i^{th}$ column the signal `x`

band-pass iltered by the $i^{th}$ band-pass filter. Hence, `size(Y, 1)=length(x)`

and `size(Y, 2)=length(f)`

.

The filter bank is designed by means of argument `bandwidth`

and optional keyword arguments `fmin`

and `fmax`

. These three arguments can be passed either as integers or real numbers. All band-pass regions have bandwidth equal to the `bandwidth`

argument and overlap with adjacent band-pass regions. By default the lower limit of the first band-pass region is set at `bandwidth`

Hz, and successive band-pass regions are defined up to, but excluding, the Nyquist frequency ($sr/2$). If `fmin`

is specified (in Hz), the center frequency of the first band-pass region is set as close as possible, but not below, `fmin`

. If `fmax`

is specified (in Hz), the center frequency of the last band-pass region is set as close as possible, but not above, `fmax`

.

Here are some examples of filter bank definitions given different arguments (`sr`

=16 in all examples).

bandwidth | fmin | fmax | center frequencies | band-pass regions |
---|---|---|---|---|

4 | - | - | 4 | 2-6 |

2 | - | - | 2, 3, 4, 5, 6 | 1-3, 2-4, 3-5, 4-6, 5-7 |

2 | 3 | 7 | 3, 4, 5, 6 | 2-4, 3-5, 4-6, 5-7 |

1 | 3 | 5 | 3, 3.5, 4, 4.5, 5 | 2.5-3.5, 3-4, 3.5-4.5, 4-5, 4.5-5.5 |

1.1 | 3 | 5 | 2.75, 3.3, 3.85, 4.4, 4.95 | 2.2-3.3, 2.75-8.85, 3.3-4.4, 3.85-4.95, 4.4-5.5 |

1.9 | 3 | 5 | 2.85, 3.8, 4.75 | 1.9-3.8, 2.85-4.75, 3.8-5.7 |

At least `sr`

samples should be trimmed at the beginning and end of the output signal `Y`

, as those samples are severely distorted by the filtering process.

If keyword optional argument β© is true (default), the filtering is multi-threaded across band-pass filters. See Threads.

This function is called by the following functions operating on time-frequency reprsentations: `TFanalyticsignal`

, `TFamplitude`

, `TFphase`

, `meanAmplitude`

, `concentration`

, `meanDirection`

, `comodulation`

, `coherence`

.

**Examples**:

```
using FourierAnalysis, DSP, Plots
# generate a sinusoidal + noise
f, sr, t = 8, 128, 512
v=sinusoidal(1., f, sr, t, 0)
x=v+randn(t)
flabels, Y=filterbank(x, 128)
flabels, Y=filterbank(x, 128; fmin=4, fmax=32)
flabels, Y=filterbank(x, 128, 4; fmin=4, fmax=32)
flabels, Y=filterbank(x, 128, 4;
filtkind=Chebyshev2(8, 10),
fmin=4,
fmax=16)
# trick for plotting the signal filtered in the band-pass regions
for i=1:size(Y, 2) Y[:, i].+=convert(eltype(Y), i)*1.5 end
mylabels=Array{String}(undef, 1, length(flabels))
for i=1:length(flabels) mylabels[1, i]=string(flabels[i])*" Hz" end
plot(Y; c=:grey, labels=mylabels)
```

`FourierAnalysis.fres`

β Method`FourierAnalysis.goertzel`

β Method```
function goertzel( x :: AbstractVector{T},
f :: IntOrReal,
sr :: Int,
wl :: Int,
startAT :: Int = 1) where T<:Union{Real, Complex}
```

Given a time series as input vector `x`

sampled at `sr`

sampling rate, return the DFT complex coefficient at the discrete Fourier frequency which is the closest to `f`

. The coefficient is computed for the time window starting at `startAt`

(one-based) and of lenght `wl`

.

If `startAT`

=1 (default) the first `wl`

samples of the vector are considered.

`wl`

does not need to be power of 2.

**See also**: `goertzel_fast`

, `goertzel2`

.

**Examples**:

```
using FourierAnalysis
sr, t, f, a = 128, 128, 5, 10
v=sinusoidal(a, f, sr, t, 0)
c=goertzel(v, f, sr, t) # should output 0+aim
```

`FourierAnalysis.goertzel2`

β Method```
function goertzel2( x :: AbstractVector{T},
f :: IntOrReal,
sr :: Int,
wl :: Int,
startAT:: Int = 1) where T<:Union{Real, Complex}
```

Like the `goertzel`

function, but allows estimating the DFT coefficient in the whole positive real line up to the Nyquist frequency. This is useful when the DFT coefficient is sought for a frequency which is not a discrete Fourier Frequency.

Using an appropriate taper this function allows almost exact recovery of both amplitude and phase at all frequencies in case of one sinusoid, whereas the `goertzel`

function can do so only at exact Fourier discrete frequencies.

**See also**: `goertzel_fast`

, `goertzel`

.

**Examples**:

```
using FourierAnalysis
sr, t, f, a = 128, 128, 5, 10
ftrue=f+0.15
v=sinusoidal(a, ftrue, sr, t, 0 )
c=goertzel(v, f, sr, t) # should be 0+aim, but clearly fails
d=goertzel2(v, ftrue, sr, t) # this get closer
```

`FourierAnalysis.goertzel_fast`

β Method```
function goertzel_fast( x :: AbstractVector{T},
wl :: Int,
a :: Real,
c :: Real,
s :: Real,
d :: Real,
startAT:: Int = 1) where T<:Union{Real, Complex}
```

Fast version of the `goertzel`

function to be preferred if the function is invoked repeatedly and speed is of concern. The user provides as arguments:

- a time series as input vector
`x`

, `wl`

, the number of samples in`x`

(or the window length),`a`

= $2/wl$,`c`

= $cos(p)$ where $p$=`2*Ο*round(UInt16, (f*wl)/sr)/wl`

, $f$ is the desired frequency and $sr$ the sampling rate,`s`

= $sin(p)$`d`

= 2*`c`

`startAt`

as in the`goertzel`

function.

`FourierAnalysis.isLinear`

β Method`function isLinear(π::Union{FDobjectsVector, TFobjectsVector})`

Return true if all objects in `π`

are linear. By definition, Spectra and TFAmplitude objects are linear. CrossSpectra, Coherence, TFAnalyticSignal and TFPhase objects may be linear or non-linear.

`FourierAnalysis.isNonLinear`

β Method`function isNonLinear(π::Union{FDobjectsVector, TFobjectsVector})`

Return true if all objects in `π`

are non-linear. By definition, Spectra and TFAmplitude objects are linear. CrossSpectra, Coherence, TFAnalyticSignal and TFPhase objects may be linear or non-linear.

`FourierAnalysis.isUnwrapped`

β Method```
(1)
function isUnwrapped(Ο΄::TFPhase)
(2)
function isUnwrapped(π―::TFPhaseVector)
```

(1) Return true if the TFPhase objects Ο΄ have the phase unwrapped.

(2) Return true if all TFPhase objects in π― have the phase unwrapped.

**See**: `unwrapPhase`

, TFPhase, TFPhaseVector.

`FourierAnalysis.meanAmplitude`

β Method```
(1)
function meanAmplitude( π :: TFAmplitudeVector,
frange :: fInterval,
trange :: tInterval;
mode :: Function = extract,
func :: Function = identity,
w :: Vector = [],
check :: Bool = true)
(2)
function meanAmplitude( π :: TFAnalyticSignalVector,
< same arguments as method (1) >
(3)
function meanAmplitude( π± :: Vector{Vector{T}},
sr :: Int,
wl :: Int,
frange :: fInterval,
trange :: tInterval,
bandwidth :: IntOrReal = 2;
mode :: Function = extract,
func :: Function = identity,
w :: Vector = [],
check :: Bool = true,
filtkind :: FilterDesign = Butterworth(2),
fmin :: IntOrReal = bandwidth,
fmax :: IntOrReal = srΓ·2,
fsmoothing :: Smoother = noSmoother,
tsmoothing :: Smoother = noSmoother,
planner :: Planner = getplanner,
β© :: Bool = true) where T<:Real
```

**alias**: mamp

(1)

Given a TFAmplitudeVector object, estimate the (weighted) mean amplitude measure across those objects. The time-frequency planes of all the objects in `π`

should be congruent.

**arguments**:

`frange`

and `trange`

define a time-frequency region on which the estimation is sought. `frange`

is a fInterval type and delimits center frequencies of the filter bank. `trange`

is a tInterval type and delimits time samples. To obtain the estimation in the whole time-frequency plane use the column (:) sign for both arguments.

**optional keyword arguments**

If `mode=extract`

is passed (default), the measure will be computed for all points in the chosen time-frequency region. If `mode=mean`

is passed, it will be computed on the mean of these points (grand-average). The `extract`

and `mean`

functions are generic methods of *FourierAnalysis*.

**Note**: with `mode=mean`

the output of the function is always a real number, whereas with `mode=extract`

the output may be a real number, a real row or column vector or a real matrix, depending on the shape of the chosen time-frequency region.

Passing a function with the `func`

argument you can derive your own time-frequency measures.

`w`

may be a vector of non-negative real weights associated to each input object. By default the unweighted version of the measure is computed.

If `check`

is true (default), check that if the column sign is passed

- as
`frange`

argument, all input objects have the same number of rows (center frequencies); - as
`trange`

argument, all input objects have the same number of columns (time samples).

If either check fails, print an error message and return `Nothing`

. No other range checks are performed.

(2)

Given a TFAnalyticSignalVector object, compute the amplitude of all objects in `π`

and estimate the (weighted) mean amplitude measure across those objects as per method (1). In addition, since using this method all TFAnalyticSignal in `π`

must be `linear`

, if `check`

is true (default) and this is not the case, print an error and return `Nothing`

. The checks on `frange`

and `trange`

performed by method (1) are also performed by this method.

(3)

Estimate the amplitude of all data vectors in `π±`

calling the `TFamplitude`

constructor and then estimate the (weighted) mean amplitude measure across the constructed amplitude objects as per method (1).

`frange`

, `trange`

, `mode`

, `func`

, `w`

and `check`

have the same meaning as in method (1). The other arguments are passed to the `TFamplitude`

constructor, to which the reader is referred for their meaning.

**See also**: `concentration`

, `meanDirection`

, timefrequencybi.jl.

**Examples**:

```
using FourierAnalysis
# generate 100 vectors of data
sr, t, bandwidth=128, 512, 2
h=taper(harris4, t)
π±=[sinusoidal(2, 10, sr, t, 0).*h.y+randn(t) for i=1:100]
π=TFanalyticsignal(π±, sr, t, bandwidth; fmax=32)
π=TFamplitude(π)
# mean amplitude in a TF region from a TFAnalyticSignalVector object
MAmp=meanAmplitude(π, (4, 16), :)
heatmap(MAmp; c=:fire) # y axis labels are not correct
# mean amplitude in a TF region from a TFAmplitudeVector object
MAmp=meanAmplitude(π, (4, 16), :)
# mean amplitude in a TF region directly from data
MAmp=meanAmplitude(π±, sr, t, (4, 16), :, bandwidth)
# NB: in the above, the analytic signal objects must all
# be linear, since meanAmplitude is computed from amplitude
# and the amplitude of non-linear analytic signal is uniformy equal to 1.
# All these computations can be obtained averaging in a TF region, e.g.,
MAmp=meanAmplitude(π, (4, 16), :; mode=mean) # output a real number
# and can be obtained on smoothed Amplitude, e.g.,
MAmp=meanAmplitude(π±, sr, t, (4, 16), :;
fsmoothing=blackmanSmoother,
tsmoothing=blackmanSmoother)
# or, equivalently, and using the alias `mamp`,
MAmp=mamp(smooth(blackmanSmoother, blackmanSmoother, π), (4, 16), :)
# A similar syntax is used for the other univariate measures, e.g.,
# concentration averaging in a TF region from a TFAnalyticSignalVector object
ConM=concentration(π, (4, 16), (128, 384); mode=mean)
# concentration in a TF region directly from data (using the alias `con`)
ConE=con(π±, sr, t, (4, 16), (128, 384), bandwidth; mode=extract)
heatmap(Con; c=:fire) # y axis labels are not correct
NB: ConM is not at all equivalent to mean(ConE) !
# mean direction averaging in a TF region directly from data
MDir=meanDirection(π±, sr, t, (4, 16), :, bandwidth; mode=mean)
# mean direction in a TF region from a TFAnalyticSignalVector object
MDir=meanDirection(π, (4, 16), :)
# and for the non-linear counterpart:
# phase concentration in a TF region directly from data
Con=concentration(π±, sr, t, (8, 12), :; nonlinear=true)
# phase concentration at a single TF point
Con=concentration(π±, sr, t, 10, 256; nonlinear=true)
# phase mean direction averaging in a TF region directly from data
# and using the alias `mdir`
MDir=mdir(π±, sr, t, (8, 12), :; mode=mean, nonlinear=true)
# If you try to compute a non-linear measure from a linear
# AnalyticSignal object you will get en error (see the REPL), e.g.,
Con=con(π, (8, 12), (1, 512); mode=mean, nonlinear=true)
# In order to compute non-linear measures from analytic signal objects
# first we need to compute non-linear analytic signal objects:
π=TFanalyticsignal(π±, sr, t, bandwidth; fmax=32, nonlinear=true)
# then, we can obtain for example the phase concentration
Con=con(π, (8, 12), :; mode=mean, nonlinear=true)
# and the phase mean direction
MDir=meanDirection(π, (8, 12), :; nonlinear=true)
```

`FourierAnalysis.meanDirection`

β Method```
(1)
function meanDirection( π :: TFAnalyticSignalVector,
frange :: fInterval,
trange :: tInterval;
nonlinear :: Bool = false,
mode :: Function = extract,
func :: Function = identity,
w :: Vector = [],
check :: Bool = true)
(2)
function meanDirection( π± :: Vector{Vector{T}},
sr :: Int,
wl :: Int,
frange :: fInterval,
trange :: tInterval,
bandwidth :: IntOrReal = 2;
nonlinear :: Bool = false,
mode :: Function = extract,
func :: Function = identity,
w :: Vector = [],
check :: Bool = true,
filtkind :: FilterDesign = Butterworth(2),
fmin :: IntOrReal = bandwidth,
fmax :: IntOrReal = srΓ·2,
fsmoothing :: Smoother = noSmoother,
tsmoothing :: Smoother = noSmoother,
planner :: Planner = getplanner,
β© :: Bool = true) where T<:Real
```

This function features two methods that use exactly the same syntax as the two corresponding methods of the `concentration`

function. All arguements have exactly the same meaning as therein. Only the output differs:

if optional keyword parameter `nonlinear`

is false (default), estimate the (weighted) mean direction measure, otherwise estimate the (weighted) phase mean direction measure.

**alias**: mdir

**See also**: `meanAmplitude`

, `concentration`

, timefrequencybi.jl.

**Examples**: see the examples of `meanAmplitude`

.

`FourierAnalysis.phase`

β Method```
(1)
function phase(z::Complex; func::Function=identity)
(2)
function phase(Z::AbstractArray{T};
unwrapdims::Int=0,
func::Function=identity) where T<:Complex
(3)
function phase(Z::TFAnalyticSignal;
unwrapped::Bool=false,
func::Function=identity)
(4)
function phase(π::TFAnalyticSignalVector;
unwrapped::Bool=false,
func::Function=identity)
```

(1)

Return the phase (argument) of a complex number. This corresponds to a standard atan2 function. It is here provided for syntactic consistency with the following methods.

(2)

Return the phase of a complex array `Z`

. Typically, `Z`

holds analytic signal, in which case the output is the analytic (instantaneous) phase. The output is a real array of the same size as `Z`

.

If optional keyword argument `unwrapdims`

is > 0, return the unwrapped phase along the `unwrapdims`

dimension of the array. For example, if `Z`

is a matrix, passing `unwrapdims=1`

unwrap the phase indipendently along its columns.

(3)

Return a real matrix with the analytic (instantaneous) phase of the TFAnalyticSignal object `Z`

. The output is of the same size as the data field `Z.y`

.

If optional keyword argument `unwrapped`

is true, return the unwrapped phase along the time dimension of the analytic signal (dims=2).

(4)

As (3), but return a vector of phase matrices for all TFAnalyticSignal objects in π―.

~

In all methods by default the phase is returned in [βΟ, Ο]. If a function is provided by the optional keyword argument `func`

, it is applied to the phase. For example

- passing
`func=x->x+Ο`

will return the phase in [0, 2Ο], - passing
`func=x->x/Ο`

will return the phase in [-1, 1], - passing
`func=sin`

will return the sine of the phase.

If in method (2) `unwrapdims`

is >0 or in method (3) and (4) `unwrapped`

is true, the function `func`

is applied to the unwrapped phase.

**See**: `unwrapPhase`

, TFAnalyticSignal.

**Examples**: see examples of `amplitude`

.

`FourierAnalysis.polar`

β Method```
(1)
function polar(c::Complex)
(2)
function polar(Z::AbstractArray{T}) where T<:Complex
(3)
function polar(Z::TFAnalyticSignal)
```

(1)

Return the amplitude (modulus) and phase (argument) of a complex number as a 2-tuple.

(2)

Return the amplitude and phase of a complex array `Z`

. Typically, `Z`

holds analytic signal, in which case return the analytic (instantaneous) amplitude and phase. The output is a tuple of two real arrays of the same size as data field `Z.y`

.

(3)

Return the analytic (instantaneous) amplitude and phase of the TFAnalyticSignal object `Z`

. The output is a tuple of two real arrays of the same size as data field `Z.y`

.

~

In all methods the phase is returned in [βΟ, Ο].

**See**: `amplitude`

, `phase`

, TFAnalyticSignal.

**Examples**: see examples of `amplitude`

.

`FourierAnalysis.sameParams`

β Function```
function sameParams(π :: TFobjectsVector,
funcname :: String) =
```

Return true if all objects in π have the same `bandwidth`

, `nonlinear`

, `fsmoothing`

and `tsmoothing`

field, otherwise print an error message pointing to the first field that is not identical in all objects and return `Nothing`

. This method applies to all TFobjectsVector types, that is, TFAnalyticSignalVector, TFAmplitudeVector and TFPhaseVector.

`funcname`

has the same meaning as in the previous method.

`FourierAnalysis.sameParams`

β Function```
function sameParams(π :: FDobjectsVector,
funcname :: String)
```

Return true if all objects in π have the same `sr`

, `wl`

, `DC`

, `taper`

, `func`

(only for SpectraVector objects), `nonlinear`

(only for CrossSpectraVector and CoherenceVector) and `smoothing`

fields, otherwise print an error message pointing to the first field that is not identical in all objects and return `Nothing`

. This method applies to all FDobjectsVector types, that is, to SpectraVector, CrossSpectraVector and CoherenceVector.

`funcname`

is an optional string that the user can provide. It is inserted into the error message to locate the part of the code that generated the error. By defalut, "unknown" is used.

`FourierAnalysis.sinusoidal`

β Function```
function sinusoidal(a :: IntOrReal,
f :: IntOrReal,
sr :: Int,
t :: Int,
ΞΈ :: IntOrReal = 0.;
DClevel = 0.)
```

Generate a sinusoidal wave with peak amplitude `a`

, frequency `f`

, sampling rate `sr`

, duration (in samples) `t`

, angle `ΞΈ`

(ΞΈ=0 makes a sine, ΞΈ=Ο/2 makes a cosine) and optional keyword argument `DC`

(float), the DC level defaulting to zero. It is adopted the convention that a sine wave starts at zero.

**Examples**:

```
using FourierAnalysis, Plots
# create and plot a sinusoidal wave of 128 samples with
# peak amplitude 1, frequency 12Hz, sr=64, phase=Ο/2
v=sinusoidal(1., 12, 64, 128, Ο/2)
plot(v)
# estimate amplitude of a sinusoidal wave using Goertzel algorithm
f, sr, t = 32, 128, 128
v=sinusoidal(3., f, sr, t, 0)
c=goertzel(v, f, sr, t) # c should be equal to 0+3.0im
```

`FourierAnalysis.slepians`

β FunctionConstruct a Taper objects holding Slepian's multi-tapering *discrete prolate spheroidal sequences*, given sampling rate `sr`

, window length `wl`

and the `bandwidth`

argument in Hz. For EEG data, 1<=bandwidth<=2 is an adequate choice.

The 'half-bandwidth' parameter `Ξ±`

used in the DSP package and in the universal Taper constructor is set as

` `Ξ±=(bandwidth/2)*wl/sr`.`

The optimal number of dpss is heuristically set to

` `n=max(1, trunc(Int, 2*Ξ±)-trunc(Int, log(2*Ξ±)))`.`

The created object can be passed as argument in constructors `spectra`

, `crossSpectra`

and `coherence`

.

**See**: plot tapering windows.

**Examples**:

```
using FourierAnalysis
sr, t, f, a = 128, 128, 10, 0.5
# create a sinusoidal superimposed to white noise
v=sinusoidal(a, f, sr, t*16, 0) + randn(t*16)
# create a data matrix
X=broadcast(+, v, randn(t*16, 3))*randn(3, 3)
# compute spectra using slepian multi-tapering with bandwidth 1.5
H=slepians(sr, t, 2)
S=spectra(X, sr, t; tapering=H)
using Plots
plot(H)
plot(S)
```

`FourierAnalysis.smooth`

β Method```
(1)
function smooth(smoothing::Smoother,
v::Vector{R}) where R<:RealOrComplex
(2)
function smooth(smoothing::Smoother,
X::Matrix{R}; dims::Int=1) where R<:RealOrComplex
(2)
function smooth(smoother :: Smoother,
S :: Union{FDobjects, FDobjectsVector})
(3)
function smooth(fsmoothing :: Smoother,
tsmoothing :: Smoother,
Y :: Union{TFobjects, TFobjectsVector})
```

Apply a smoothing function of type Smoother to

- (1) a vector of real or complex numbers,
- (2) a real of complex matrix along dimension
`dims`

(default=1), - (3) a FDobjects or all objects in a FDobjectsVector,
- (4) a TFobjects or all objects in a TFobjectsVector.

Methods (1) and (2) are provided for low-level computations. Methods (3) and (4) are constructors; for all methods the output is always of the same type as the input.

Method (3) smooths across the frequency dimension:

- for Spectra objects this amounts to smoothing the column vectors in their
`.y`

field, - for CrossSpectra and Coherence objects this amounts to smoothing adjacent matrices in their .y field.

Method (4) smooths across the frequency dimension, time dimension or both. This amounts to smoothing across the column vectors (frequency) and/or row vectors (time) in the `.y`

field of the object. A smoother must be specified for the frequency dimension (`fsmoothing`

) and for the time dimension (`tsmoothing`

). Either one may be `noSmoother`

, but if the two are different from `noSmoother`

, then they must be the same. If smoothing is requested in both the frequency and time dimension, then the data is smoothed first in the time then in the frequency dimension. For TFPhase objects, smoothing is allowed only if the phase is unwrapped.

This function allow smoothing frequency domain and time-frequency domain objects after they have been created, however, smoothing can also be requested upon creation. For example, see the documentation of Spectra.

For methods (1), (2) and (3), if `Smoother`

is `noSmoother`

, then the input is returned unchanged. For method (4) this is the case if both `fsmoother`

and `tsmoother`

are `noSmoother`

.

The data input must hold in the concerned dimension at least three elements for applying an Hann or Hamming smoother and at least five elements for applying the Blackman smoother.

**Maths**

Smoothing of a series $x$ composed of $k$ elements is carried out at element $i$ such as

$x_{i}=ax_{i-2}+bx_{i-1}+cx_{i}+bx_{i+1}+ax_{i+2}$.

The coefficients are

smoothing window | a | b | c |
---|---|---|---|

Hann | 0 | 0.25 | 0.50 |

Hamming | 0 | 0.23 | 0.54 |

Blackman | 0.04 | 0.25 | 0.42 |

For 3-point smoothers, the first point is smoothed as

$x_{1}=\frac{c}{b+c}x_{1} + \frac{b}{b+c}x_{2}$

and the last (the $k^{th}$) as

$x_{k}=\frac{c}{b+c}x_{k} + \frac{b}{b+c}x_{k-1}$.

For 5-point smoothers, the first point is smoothed as

$x_{1}=\frac{c}{a+b+c}x_{1} + \frac{b}{a+b+c}x_{2} + \frac{a}{a+b+c}x_{3}$,

the second as

$x_{2}=\frac{b}{a+2b+c}x_{1} + \frac{c}{a+2b+c}x_{2} + \frac{b}{a+2b+c}x_{3} + \frac{a}{a+2b+c}x_{4}$,

the second to last as

$x_{k-1}=\frac{a}{a+2b+c}x_{k-3} + \frac{b}{a+2b+c}x_{k-2} + \frac{c}{a+2b+c}x_{k-1} + \frac{b}{a+2b+c}x_{k}$

and the last as

$x_{k}=\frac{a}{a+b+c}x_{k-2} + \frac{b}{a+b+c}x_{k-1} + \frac{c}{a+b+c}x_{k}$.

**See**: Smoother

**Examples**:

```
using FourierAnalysis, Plots
sr, t, f, a = 128, 128, 10, 0.5
# create a sinusoidal superimposed to white noise
v=sinusoidal(a, f, sr, t*16, 0) + randn(t*16)
# compute Amplitude Spectra
Ξ£=spectra(v, sr, t; func=β)
bar(Ξ£.y, labels="raw amplitude spectra")
#smooth spectra
Ξ£2=smooth(blackmanSmoother, Ξ£)
bar!(Ξ£2.y, labels="smoothed amplitude spectra")
# smooth cross-spectra (or coherence) matrices
X=broadcast(+, v, randn(t*16, 3))*randn(3, 3)
S=crossSpectra(X, sr, t) # or coherence (X, sr, t)
# smooth the cross-spectra # or coherence
S2=smooth(blackmanSmoother, S)
# smooth time-frequency object
Y = TFanalyticsignal(v, sr, sr*4)
# smooth frequency
Z=smooth(blackmanSmoother, noSmoother, Y)
# plot amplitude of smoothed analytic signal
heatmap(Z, amplitude)
# smooth AS: smooth both frequency and time
E=smooth(blackmanSmoother, blackmanSmoother, Y)
# plot real part of smoothed analytic signal
heatmap(Z, real)
```

`FourierAnalysis.spectra`

β Method```
(1)
function spectra( X :: Union{Vector{T}, Matrix{T}},
sr :: Int,
wl :: Int;
tapering :: Union{Taper, TaperKind} = harris4,
planner :: Planner = getplanner,
slide :: Int = 0,
DC :: Bool = false,
func :: Function = identity, # equal to x->x
smoothing :: Smoother = noSmoother,
β© :: Bool = true) where T<:Real
(2)
function spectra( π :: Vector{Matrix{T}},
< same argument sr, ..., β© of method (1) > where T<:Real
```

(1)

Construct a Spectra object from real univariate or multivariate data. Given sampling rate `sr`

and epoch length `wl`

, compute the Welch power spectra of a vector (univariate) or of a data matrix (multivariate) `X`

of dimension $t$x$n$, where $t$ is the number of samples (rows) and $n$ is the number of time-series (columns).

The spectra are hold in the `.y`

field of the created object. If `X`

is a vector, `.y`

is a vector, whereas if `X`

is a matrix, `.y`

is a matrix holding in its columns the spectra of the signals in the columns of `X`

. The size of the spectra depends on the `DC`

optional keyword argument (see below), as reported in the documentation of the Spectra structure.

**Optional Keyword Arguments**:

`tapering`

: this is a tapering object of type Taper or a tapering window kind of type TaperKind. By default the *harris4* tapering window is used. If no tapering is sought, pass as argument `tapering=rectangular`

. This same syntax is the most convenient way to specify all simple tapering window, e.g., `tapering=hann`

, `tapering=hamming`

, etc. For *discrete prolate spheroidal sequences (dpss)* multi-tapering, use instead the `slepians`

constructor, e.g., pass as argument something like `tapering=slepians(sr, wl, 2)`

.

`planner`

: this is an instance of the `Planner`

object, holding the forward FFTW plan used to compute the FFTs. By default the planner is computed by this method, but it can be passed as an argumet here if it is pre-computed. This is interesting if this function is to be invoked repeatedly.

`slide`

: this is the number of samples the windows slide on (Welch method). By default the number of samples is chosen to allow 50% overlap.

`DC`

: if true, the spectrum/a of the DC level is returned in the first row of `y`

(see the fields of the Spectra object), otherwise (default) the rows in `y`

start with the first positive discrete frequency, that is, $sr/wl$ Hz.

`func`

: this is a function that operate element-wise to transfrom the power spectra before returning them, including anonymous functions. Common choices are:

`func=sqrt`

return the amplitude spectra,`func=log`

return the log-power spectra,`func=decibel`

return the power spectra in deciBels (see`decibel`

),

By default no function is applied and the power spectra are returned. If smoothing has been requested (see below), it is applied after the function.

`smoothing`

: it applies a smoothing function of type Smoother to the spectra across adjacent frequencies. By default no smoothing is applied.

`β©`

: if true (default), the method run in multi-threaded mode across the series in `X`

if the number of series is at least twice the number of threads Julia is instructed to use. See Threads.

(2)

Construct a SpectraVector object from a vector of real univariate (vectors) or multivariate data (matrices). Compute the spectra as per method (1) for all $k$ data vectors or data matrices in `π`

.

If `π`

holds data matrices, the $k$ matrices in `π`

must have the same number of columns (i.e., the same number of time series), but may have any number of (at least `wl`

) rows (samples). All other arguments have the same meaning as in method (1), with the following differences:

`β©`

: if true (default), the method run in multi-threaded mode across the $k$ data matrices if $k$ is at least twice the number of threads Julia is instructed to use, otherwise this method attempts to run each spectra estimation in multi-threaded mode across series as per method (1). See Threads.If a

`planner`

is not explicitly passed as an argument, the FFTW plan is computed once and applied for all spectral estimations.

**See**: Spectra, plot spectra.

**See also**: `crossSpectra`

, `coherence`

, goertzel.jl.

**Examples**:

```
using FourierAnalysis
###################################################################
# (1)
# Check that correct amplitude spectra is returned at all discrete
# Fourier Frequency (using a rectangular taper).
# Sinusoids are created at all frequencies with amplitude 10 at the
# first frequency and then incrementing by 10 units along frequencies.
# NB when t is even, correct amplitude for the last frequency is obtained
# only if the sinusoidal has a particular phase.
sr, t, wl= 16, 32, 16
V=Matrix{Float64}(undef, t, wl)
for i=1:wl V[:, i]=sinusoidal(10*i, b2f(i, sr, t), sr, t, Ο/6) end
# using FFTW.jl only
using FFTW
P=plan_rfft(V, 1)*(2/t);
Ξ£=abs.(P*V)
using Plots
bar(Ξ£[brange(t, DC=true), :], labels="")
# using FourierAnalysis.jl
Ξ£2=spectra(V, sr, t; tapering=rectangular, func=β, DC=true)
using Plots
bar(Ξ£2.y[brange(t, DC=true), :], labels="")
#############################################################################
### Check amplitude spectra on long signals obtained by welch methods
# one sinusoidal is at an exact discrete Fourier Frequency and the other not
# Rectangular window
sr, t, f, a = 128, 128, 10, 0.5
v=sinusoidal(a, f, sr, t*16)+sinusoidal(a, f*3.5+0.5, sr, t*16)+randn(t*16);
Ξ£=spectra(v, sr, t; tapering=rectangular, func=β)
bar(Ξ£.y, labels="rectangular")
# harris4 window (default)
Ξ£2=spectra(v, sr, t; func=β)
bar!(Ξ£2.y, labels="harris4")
#smooth spectra
Ξ£3=smooth(blackmanSmoother, Ξ£2)
bar!(Ξ£3.y, labels="smoothed")
#############################################################################
function generateSomeData(sr::Int, t::Int; noise::Real=1.)
# four sinusoids of length t samples and sr sampling rate
# peak amplitude: 0.7, 0.6, 0.5, 0.4
# frequency: 5, 7, 13, 27
# phase: 0, Ο/4, Ο/2, Ο
v1=sinusoidal(0.7, 5, sr, t, 0)
v2=sinusoidal(0.6, 7, sr, t, Ο/4)
v3=sinusoidal(0.5, 13, sr, t, Ο/2)
v4=sinusoidal(0.4, 27, sr, t, Ο)
return hcat(v1, v2, v3, v4) + (randn(t, 4)*noise)
end
sr, wl, t = 128, 512, 8192
X=generateSomeData(sr, t)
# multivariate data matrix 8192x4
# compute spectra
S=spectra(X, sr, wl)
# check the spectrum of first series
S.y[:, 1]
# gather some plot attributes to get nice plots
using Plots.Measures
spectraArgs=(fmax = 32,
left_margin = 2mm,
bottom_margin = 2mm,
xtickfont = font(11, "Times"),
ytickfont = font(11, "Times"))
plot(S; spectraArgs...)
plot(S; xspace=2, spectraArgs...)
# use a custom simple taperig window
S=spectra(X, sr, wl; tapering=riesz)
# use Slepian's multi-tapering
S=spectra(X, sr, wl; tapering=slepians(sr, wl, 1.5))
# construct with smoothing
S=spectra(X, sr, wl; tapering=slepians(sr, wl, 1.5), smoothing=hannSmoother)
# compute Amplitude Spectra instead
S=spectra(X, sr, wl; tapering=slepians(sr, wl, 1.5), func=β)
# plot Aplitude spectra
plot(S; ytitle="Amplitude", spectraArgs...)
# smooth the spectra a-posteriori
S=smooth(blackmanSmoother, S)
# extract spectra in range (8Hz to 12Hz)
e=extract(S, (8, 12))
# extract spectra in range (8Hz to 12Hz) for series 1 and 2
e=extract(S, (8, 12))[:, 1:2]
# extract the spectra at 10Hz only (1 bin per series)
e=extract(S, 10)
# average spectra in the 8Hz-12Hz range
bar(mean(S, (8.0, 12.0)))
# average across series of the average spectra in the 8Hz-12Hz range
mean(mean(S, (8.0, 12.0)))
# average spectrum across all frequencies for each series
bar(mean(S, :))
# average spectra in equally-spaced 2Hz-band-pass regions for all series
plot(bands(S, 2))
# average spectra in equally-spaced 2Hz-band-pass regions for series 1 and 2
plot(bands(S, 2)[:, 1:2])
# (2)
# generate 3 multivariate data matrices 8192x4
X=[generateSomeData(sr, t) for i=1:3]
# Now the call to the spectra function will generate 3 Spectra objects
S=spectra(X, sr, wl)
plot(S[1]; spectraArgs...)
plot(S[2]; spectraArgs...)
plot(S[3]; spectraArgs...)
# when you want to compute the spectra of many data matrices you may want
# to do it using a fast FFTW plan (wait up to 10s for computing the plan)
plan=Planner(plan_exhaustive, 10.0, wl)
S=spectra(X, sr, wl; planner=plan)
# how faster is this?
using BenchmarkTools
@benchmark(spectra(X, sr, wl))
@benchmark(spectra(X, sr, wl; planner=plan))
# average spectra in range (8Hz-12Hz) for all series of all objects
M=mean(S, (8, 12))
# plot the average spectrum across all series for the three objects
# using Julia standard mean function
plot(mean(S[1].y[:, i] for i=1:size(S[1].y, 2)))
plot!(mean(S[2].y[:, i] for i=1:size(S[2].y, 2)))
plot!(mean(S[3].y[:, i] for i=1:size(S[3].y, 2)))
# average spectra in range (4Hz-32.5Hz) across objects for the 4 series
plot(mean(mean(S, (4, 32.5))))
# extract spectra in range (8Hz-12Hz) for all series and all subjects
extract(S, (8, 12))
# if you enter en illegal range, nothing will be done and you will get
# an error in the REPL explaining what is wrong in your argument
extract(S, (0, 12))
mean(S, (1, 128))
# extract 4Hz-band-pass average spectra for all series and all objects
bands(S, 4)
# Apply smoothing in the spectra computations
S=spectra(X, sr, t; smoothing=blackmanSmoother)
plot(S[1]; spectraArgs...)
plot(S[2]; spectraArgs...)
plot(S[3]; spectraArgs...)
# plot spectra in in 1Hz band-pass regions for all series in S[1]
plot(bands(S[1], 1))
# use slepian multi-tapering
S=spectra(X, sr, wl; tapering=slepians(sr, wl, 1.))
plot(S[1]; spectraArgs...)
# average spectra across objects
plot(mean(s.y for s β S))
```

`FourierAnalysis.taper`

β Method```
function taper( kind :: TaperKind,
wl :: Int;
Ξ± :: Real = 2.,
n :: Int = ceil(Int, 2*Ξ±)-1,
padding :: Int = 0,
type :: Type{T} = Float64) where T<:Union{Real, Complex}
```

Universal constructor of Taper objects, given a tapering window `kind`

, of type TaperKind and the window length `wl`

. NB: for the `hamming`

and `blackman`

type use `FourierAnalysis.hamming`

and `FourierAnalysis.blackman`

to avoid conflicts with the DSP package.

Return a vector of length `wl`

for all types of tapers, but for the dpss (Slepian multi-tapers), for which return a matrix of size `wl`

x `n`

.

If optional keyword argument `padding`

is >0, then the actual window length will be `wl`

+ `padding`

.

The `type`

optional keyword argument can be use to specify the the type of elements of the typering window. By default, this is the `Float64`

type.

Optional keywords arguments `Ξ±`

and `n`

currently apply only for slepian multi-tapering (*discrete prolate spheroidal sequences*):

`Ξ±`

is the *half bandwidth* (hbw) parameter as per the DSP.jl package. This unit is used in many dpss implementations and it is often reported that "typical values" for `Ξ±`

are 2, 2.5, 3, 3.5 or 4. However, the optimal smoothing increases with the ratio between window length and sampling rate, thus these values in absolute terms are not useful. In fact, the larger the *hbw* and the higher `n`

, the smoother the spectra will be (variance reduction). In order to overcome these difficulties, for Slepian's multitapering *FourirAnalysis* implements the `slepians`

constructor, which allows the bandwidth parameter to be given in Hz and the number of tapering windows to be chosen automatically.

`n`

is the number of tapering windows. For slepian tapers this is the number of the discrete prolate spheroidal sequences. As in the DSP package, by default it is set to `ceil(Int, 2*Ξ±)-1`

, however, depending on how large this number is, low eigenvalues may correspond to the last sequences, therefore those should be discarded.

**See**: plot tapering windows.

**Examples**:

```
using FourierAnalysis
## Use the constructor
sr, t, f, a = 128, 128, 10, 0.5
# create a sinusoidal superimposed to white noise
v=sinusoidal(a, f, sr, t*16, 0) + randn(t*16)
# create a data matrix
X=broadcast(+, v, randn(t*16, 3))*randn(3, 3)
# compute spectra using hamming tapering window
# we need to prepend 'FourierAnalysis.' since `hamming`
# is declared also in DSP.jl
H=taper(FourierAnalysis.hamming, t)
S=spectra(X, sr, t; tapering=H)
# you can obtain the same thing with
S=spectra(X, sr, t; tapering=H)
# which will create the hamming tapering window on the fly,
# thus calling explicitly the constructor is interesting
# only if you need to reuse the same tapering window many times.
## Plot tapering windows using the standard plot function
using Plots
tapers=[TaperKind(i) for i=1:8]
X=zeros(t, 8)
for i=1:8 X[:, i] = taper(tapers[i], t).y end
mylabels=Array{String}(undef, 1, 8)
for i=1:8 mylabels[1, i]=string(tapers[i]) end
plot(X; labels=mylabels)
## using the recipe declared in recipes.jl
plot(taper(parzen, 256))
plot(taper(slepian, 256, Ξ±=4, n=7))
```

`FourierAnalysis.taperinfo`

β Method`function taperinfo(taper::Taper)`

Return the name of the tapering window(s) encapsulated in the Taper object as a string.

Only for Slepian's discrete prolate spheroidal sequences (dpss), their parameters, namely, $Ξ±$ (half-bandwidth) and $n$ (number of windows), are reported within parentheses as well.

**Examples**:

```
H=taper(hamming, 128*8)
taperinfo(H)
H=slepians(128, 128*8, 2)
taperinfo(H)
```

`FourierAnalysis.tfAxes`

β MethodGenerate labels for the axes to be used in heatmap plots of TFobjects. On the y-axis the labels correspond to the center frequencies of the band-pass filter bank used to obtain the time-frequency representation, that is, to the `.flabels`

field of the object.

See plot time-frequency objects for examples.

`FourierAnalysis.unwrapPhase`

β Method```
(1)
function unwrapPhase(Z::AbstractArray{T};
unwrapdims::Int=0) where T<:Complex
(2)
function unwrapPhase(Ο΄::AbstractArray{T};
unwrapdims::Int=0) where T<:Real
(3)
unwrapPhase(Ο΄::TFPhase) [constructor of a TFPhase object]
(4)
unwrapPhase(π―::TFPhaseVector) [constructor of a TFPhaseVector object]
```

(1)

If optional keyword argument `unwrapdims`

is > 0, compute the phase (argument) from a *complex* array and unwrap it along the `unwrapdims`

dimension, otherwise (default) return `Z`

. Typically, `Z`

holds analytic signal.

(2)

If optional keyword argument `unwrapdims`

is > 0, unwrap along the `unwrapdims`

dimension a *real* array holding phase data in [βΟ, Ο], otherwise return `Ο΄`

.

(3)

Construct a TFPhase object by unwrapping its phase along the time dimension and copying all other fields from the `Ο΄`

object. If `Ο΄.func`

is different from the `identity`

(do nothing) function, return instead an error message.

(4)

As (3), but conctruct a TFPhaseVector holding TFPhase objects in π― with the phase unwrapped. `Ο΄.func`

must be the identity function for all Ο΄ β π―.

The unwrapped phase is defined as the cumulative sum of the phase (along the relevant dimension) once this is represented in [0, 2Ο].

**Examples**: see examples of `amplitude`

.

`Statistics.mean`

β Method```
(1)
function mean(S :: FDobjects,
frange :: fInterval)
(2)
function mean(π :: FDobjectsVector,
frange :: fInterval;
w :: Vector = [],
check :: Bool = true)
(3)
function mean(Y :: TFobjects,
frange :: fInterval,
trange :: tInterval)
(4)
function mean(π :: TFobjectsVector,
frange :: fInterval,
trange :: tInterval;
w :: Vector = [],
check :: Bool = true)
```

Return the mean of data in a frequency region from FDobjects and data in a time-frequency region from TFobjects. The frequency and time region are indicated by `frange`

and `trange`

, which are of type fInterval and tInterval, respectively.

The complete input/output types for this function is reported in the following table:

method | input object | output |
---|---|---|

(1.1) | Spectra | a vector holding the mean spectra in `frange` ΒΉ |

(1.2) | CrossSpectra | a complex matrix holding the mean cross-spectra in `frange` Β² |

(1.3) | Coherence | a real matrix holding the mean coherence in `frange` Β² |

(2.1) | SpectraVector | a vector of vectors of type (1.1) |

(2.2) | CrossSpectraVector | a vector of matrices of type (1.2) |

(2.3) | CoherenceVector | a vector of matrices of type (1.3) |

(3.1) | TFAnalyticSignal | a complex number holding the mean analytic signal in [`frange` , `trange` ] |

(3.2) | TFAmplitude | a real number holding the mean amplitude in [`frange` , `trange` ] |

(3.3) | TFPhase | a real number holding the mean phase in [`frange` , `trange` ] |

(4.1) | TFAnalyticSignalVector | a vector of numbers of type (3.1) |

(4.2) | TFAmplitudeVector | a vector of numbers of type (3.2) |

(4.3) | TFPhaseVector | a vector of numbers of type (3.3) |

legend: ΒΉ*each element of the vector refers to a time-series on which the spectra have been computed.* Β² *depending on how the objects has been created, the matrices may be either Hermitian or LowerTriangular.*

Method (2) and (4) allows the following *optional keyword arguments*:

`w`

, a $k$-vector of non-negative integers or real numbers, where $k$ is the numbers of objects hold in the input FDobjectsVector or TFobjectsVector. `w`

is a vector of weights for the means extracted from the input objects. By default, no weights are assigned.

`check`

, a boolean. If it is true (default), it is checked that the non-data fields of the input objects are all the same (for example, sampling rate, bandwidth, etc.).

**See also**: `extract`

.

**Examples**:

```
using FourierAnalysis, Plots
# example with univariate Spectra objects (one series -> one spectrum)
sr, t, f, a = 128, 256, 10, 1
# create a sinusoidal superimposed to white noise
v=sinusoidal(a, f, sr, t*16, 0) + randn(t*16)
# compute the spectrum
Ξ£=spectra(v, sr, t)
# mean spectrum in between 8Hz and 12Hz
s=mean(Ξ£, (8, 12))
# mean spectrum in between 8Hz and 12.5Hz
s=mean(Ξ£, (8, 12.5))
# example with multivariate spectra (several series -> several spectra)
Ξ£=spectra(hcat(v, v+randn(t*16)), sr, t)
# mean spectra in between 8Hz and 12Hz
S=mean(Ξ£, (8, 12))
# mean spectra at 10Hz, i.e., the spectra at 10Hz
S=mean(Ξ£, 10)
# example with CrossSpectra objects (the same goes for Coherence objects)
X=broadcast(+, v, randn(t*16, 3))*randn(3, 3)
Ξ£=crossSpectra(X, sr, t)
# mean cross-spectra in between 8Hz and 12Hz (an Hermitian matrix)
S=mean(Ξ£, (8, 12))
Ξ£=crossSpectra(X, sr, t; tril=true)
# mean cross-spectra in between 8Hz and 12Hz (a LowerTriangular matrix)
S=mean(Ξ£, (8.0, 12.0)) # accept integers and reals for frequencies
# example with multiple CrossSpectra objects
X2=broadcast(+, v, randn(t*16, 3))*randn(3, 3)
Ξ£=crossSpectra([X, X2], sr, t) # this is a CrossSpectraVector
S=mean(Ξ£, (8, 12); w=[0.4, 0.6])
# now S[1] will hold the mean cross-spectrum in range 8-12Hz for X
# and S[2] will hold the mean cross-spectrum in range 8-12Hz for X2
# example with time-frequency objects
# (work in the same way for TFAnalyticSignal, TFAmplitude and TFPhase)
Y = TFanalyticsignal(v, sr, t)
# mean analytic signal within frequencies 8Hz and 12Hz and time samples 1 to 64.
as=mean(Ξ£, (8, 12), (1, 64))
# mean analytic signal within frequencies 8Hz and 12Hz.
as=mean(Ξ£, (8, 12), :)
# mean analytic signal within time samples 1 to 64.
as=mean(Ξ£, :, (1, 64))
# example with multiple time-frequency objects
Y = TFanalyticsignal([v, v+randn(t*16)], sr, t)
AS=mean(Y, (8, 12), (1, 64))
# get the mean across TFobjects of those means
m=mean(mean(Y, (8, 12), (1, 64)))
AS=mean(Y, (8), :)
AS=mean(Y, 8, 2)
```