# spectra.jl

Spectra objects created by *FourierAnalysis* are incapsulated in the following structure:

## Spectra

**Categories**: data objects, FDobjects

```
struct Spectra
y :: AbstractArray{T} where T<:Real
sr :: Int
wl :: Int
DC :: Bool
taper :: String
flabels :: Vector{T} where T<:Union{Real, Int}
func :: Function
smoothing :: Smoother
end
```

**Fields**:

`y`

: a real vector holding one *spectrum*, or a real matrix holding several *spectra* in its columns.

`sr`

: the *sampling rate* of the data on which the spectra have been estimated.

`wl`

: the FFT *window length* used for estimating the spectra.

`DC`

: if true, the first row of `y`

holds the *DC level*, otherwise it holds the first positive frequency. Thus, if `DC`

is false, the first dimension of `y`

is equal to $wl÷2$ (integer division), otherwise it is equal to $(wl÷2)+1$ (see Overview). In all constructors it is false by default.

`taper`

: the time-domain *tapering window* used for FFT computation, as a string, with parameters in parentheses for Slepian's *dpss*. See tapers.jl.

`flabels`

: a vector holding all Fourier discrete frequencies in Hz. Those are the *frequency labels* for the rows of `y`

. If `DC`

is true, the first label is $0$, otherwise it is the first positive frequency, which is equal to the frequency resolution $sr/wl$.

`func`

: a *function* applied element-wise to the spectra. In all constructors the default is the `identity`

(do-nothing) function.

`smoothing`

: a Smoother flag indicating whether the spectra have been smoothed across adjacent frequencies. If no smoothing has been applied, it is equal to `noSmoother`

, which is the default in all constructors.

**Note**: In Julia the fields are accessed by the usual dot notation, e.g., you may verify that for *Spectra* object `S`

,

`length(S.flabels) == size(S.y, 1) == (S.wl/2)+S.DC`

A vector of *Spectra* objects is of type SpectraVector.

**Methods for Spectra and SpectraVector objects**

method | Spectra | SpectraVector |
---|---|---|

`bands` | ✔ | ✔ |

`extract` | ✔ | ✔ |

`mean` | ✔ | ✔ |

`smooth` | ✔ | ✔ |

`sameParams` | ✔ |

**Generic Constructors**:

In order to construct a *Spectra* object from univariate and multivariate data, *FourierAnalysis* provides two `spectra`

constuctors, which is what you will use in practice most of the time.

Manual constructors are also possible, for which you have to provide appropriate arguments. The default manual constructor of *Spectra* objects is

`Spectra(y, sr, wl, DC, taper, flabels, func, smoothing).`

Other constructors are also provided:

`Spectra(y, sr, wl, DC, tapers)`

generates the appropriate `flabels`

, set `func`

to `identity`

(do-nothing) and `smoothing`

to `noSmoother`

;

`Spectra(y, sr, wl, taper)`

is like the constructor here above, but also set `DC`

to false;

```
Spectra(𝙎::CrossSpectra;
func::Function=identity)
```

creates a *Spectra* object extracting the spectra from a CrossSpectra object. If a function is provided with the `func`

argument, this function is applied element-wise to the spectra. For instance,

`func=sqrt`

will extract amplitude spectra,`func=log`

will extract log-spectra,`func=decibel`

will extract spectra in deciBels (see`decibel`

).

By default the `identity`

(do-nothing) function is applied, thus (power) spectra are extracted;

If the CrossSpectra object is non-linear, the spectra are uniformly equal to 1.0. See `crossSpectra`

.

```
Spectra(𝙎::CrossSpectraVector;
func::Function=identity)
```

create a SpectraVector object from a CrossSpectraVector object, calling the constructor here above for all `crossSpectra`

objects hold by `𝙎`

.

**Constructors from data**:

When the spectrum is recquired only at one or a few discrete Fourier coefficients, Goertzel's algorithhms may be a more efficient option. See unit goertzel.jl.

`FourierAnalysis.spectra`

— Function```
(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
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))
```