FFTW 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)

(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

(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))
...
...

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

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

(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:

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.

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

(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
𝛷

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

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)

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

bin to frequency. 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

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)

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]

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

(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)

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

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

See also: spectra, coherence.

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

(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)

(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:

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)

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

frequency to bin. 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

function fbands(sr :: Int,
                wl :: Int,
         bandwidth :: IntOrReal;
      DC :: Bool = false)

Return a vector of Frequencies (in Hz) to which the bins created by a call to function bbands with the same arguments correspond.

See: bbands.

See also: bands.

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

Return a vector with all Fourier discrete frequencies 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]

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

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)

function fres(sr :: Int,
              wl :: Int)

FFT frequency resolution given sampling rate sr and window length wl.

See also: f2b, b2f, fdf, brange.

Examples:

using FourierAnalysis
fres(1024, 2048) # return 0.5

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

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

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.

See also: goertzel, goertzel2.

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.

See:FDobjectsVector, TFobjectsVector.

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.

See:FDobjectsVector, TFobjectsVector.

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

(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)

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

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

See: unwrapPhase, TFAnalyticSignal.

Examples: see examples of amplitude.

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

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.

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.

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

Construct 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)

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

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

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)

(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))

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.

See also: slepians, taperinfo

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

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)

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

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

(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:

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)