tools.jl

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

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

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

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

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

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

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