Periodograms - periodogram estimation

arraysplit(s, n, noverlap, nfft=n, window=nothing; buffer=zeros(eltype(s), nfft))

Split an array s into arrays of length n with overlapping regions of length noverlap.

Arguments

• s: Input array.
• n: Specifies the required subarray length (n ≤ length(s)).
• noverlap: Number of overlapping elements between subarrays.
• nfft: Specifies the length of the split arrays. If length(s) < nfft, then the input is padded with zeros.
• window: An optional scaling vector to be applied to the split arrays.
• buffer: An optional buffer of length nfft. Iterating or indexing the returned AbstractVector always yields the same Vector with different contents. The last result after iteration or indexing calls is stored into buffer.

Returns

An ArraySplit object with split subarrays.

Examples

julia> arraysplit([0.1, 0.2, 0.3, 0.4, 0.5], 3, 1)
2-element DSP.Periodograms.ArraySplit{Vector{Float64}, Float64, Nothing}:
 [0.1, 0.2, 0.3]
 [0.3, 0.4, 0.5]

julia> arraysplit([0.1, 0.2, 0.3, 0.4, 0.5], 3, 2, 8)
3-element DSP.Periodograms.ArraySplit{Vector{Float64}, Float64, Nothing}:
 [0.1, 0.2, 0.3, 0.0, 0.0, 0.0, 0.0, 0.0]
 [0.2, 0.3, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0]
 [0.3, 0.4, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0]

julia> arraysplit([0.1, 0.2, 0.3, 0.4, 0.5], 3, 1, 3, [1, 2, 1])
2-element DSP.Periodograms.ArraySplit{Vector{Float64}, Float64, Vector{Int64}}:
 [0.1, 0.4, 0.3]
 [0.3, 0.8, 0.5]

arraysplit function with buffer

julia> x = [1, 2, 3, 4, 5];

julia> sub_arr, n_overlap, nfft = 3, 1, 3;

julia> x_split = arraysplit(x, sub_arr, n_overlap, nfft, nothing; buffer=zeros(nfft));

julia> x_split[2]   # Returns AbstractVector result and stores it in ArraySplit.buf and buffer
3-element Vector{Float64}:
 3.0
 4.0
 5.0

julia> x_split.buf  # Returns stored results from previous x_split[2] call
3-element Vector{Float64}:
 3.0
 4.0
 5.0

ArraySplit{T<:AbstractVector,S,W} <: AbstractVector{Vector{S}}

ArraySplit object with fields

• s::T
• buf::Vector{S}
• n::Int
• noverlap::Int
• window::W
• k::Int

where buf: buffer, k: number of split arrays, and s, n, noverlap, window are as per arguments in arraysplit.

periodogram(s::AbstractVector; onesided=eltype(s)<:Real, nfft=nextfastfft(length(s)), fs=1, window=nothing)

Computes periodogram of a 1-d signal s by FFT and returns a Periodogram object.

Arguments

• onesided: For real signals, the two-sided periodogram is symmetric and this function returns a one-sided (real only) periodogram by default. A two-sided periodogram can be obtained by setting onesided=false.
• nfft: Specifies the number of points to use for the Fourier transform. If length(s) < nfft, then the input is padded with zeros. By default, nfft is the closest size for which the Fourier transform can be computed with maximal efficiency.
• fs: The sample rate of the original signal.
• window: An optional window function or vector to be applied to the original signal before computing the Fourier transform. The computed periodogram is normalized so that the area under the periodogram is equal to the uncentered variance (or average power) of the original signal.

Returns

A Periodogram object with the 2 computed fields: power, freq.

Examples

Frequency estimate of cos(2π(25)t) with a 1-sided periodogram.

julia> Fs = 100;                     # Sampling frequency

julia> t = (1:Fs)/Fs;                # 100 time samples

julia> x = cos.(2π*25*t);            # Cosine signal with 25Hz frequency

julia> pxx = periodogram(x; fs=Fs);                 # Periodogram

julia> _, max_index = findmax(pxx.power);           # Find index of max power

julia> pxx.power[max_index], pxx.freq[max_index]    # (power, frequency)
(0.5, 25.0)

2-sided periodogram of a rectangle function with Hamming window.

julia> x = rect(50; padding=50);     # Rectangular pulse function

julia> pxx = periodogram(x; onesided=false, nfft=512, fs=1000, window=hamming);

julia> round(maximum(pxx.power), digits=4)  # Expect max power to be close to 0Hz for a DC signal
0.0182

periodogram(s::AbstractMatrix; nfft::NTuple{2,Int}=nextfastfft(size(s)), fs=1, radialsum=false, radialavg=false)

Computes periodogram of a 2-d signal using the 2-d FFT and returns a Periodogram2 or Periodogram object.

Arguments

• nfft: A 2-tuple specifying the number of points to use for the Fourier transform for each respective dimension. If size(s) < nfft, then the input is padded with zeros. By default, nfft is the closest size for which the Fourier transform can be computed with maximal efficiency.
• fs: The sample rate of the original signal in both directions.
• radialsum: For radialsum=true, the value of power[k] is proportional to $\frac{1}{N}\sum_{k\leq |k'|<k+1} |X[k']|^2$.
• radialavg: For radialavg=true, the value of power[k] is proportional to $\frac{1}{N \#\{k\leq |k'|<k+1\}} \sum_{k\leq |k'|<k+1} |X[k']|^2$. The computation of |k'| takes into account non-square signals by scaling the coordinates of the wavevector accordingly.

Returns

• A Periodogram2 object by default with the 3 computed fields: power, freq1, freq2.
• A Periodogram object is returned for a radially summed or averaged periodogram (if radialsum or radialavg is true, respectively). Only one of radialsum or radialavg can be set to true in the function.

Examples

julia> x = [1 1; 0 1; 0 0];

julia> pxx = periodogram(x; nfft=(3, 2));   # Returns Periodogram2

julia> power(pxx)
3×2 Matrix{Float64}:
 1.5  0.166667
 0.5  0.166667
 0.5  0.166667

julia> freq(pxx)
([0.0, 0.3333333333333333, -0.3333333333333333], [0.0, -0.5])

Radial sum and radial average periodograms

julia> x = [1 3; 0 1];

julia> periodogram(x; radialsum=true)  # Returns Periodogram
DSP.Periodograms.Periodogram{Float64, AbstractFFTs.Frequencies{Float64}, Vector{Float64}}([6.25, 4.75], [0.0, 0.5])

julia> periodogram(x; radialavg=true)  # Returns Periodogram
DSP.Periodograms.Periodogram{Float64, AbstractFFTs.Frequencies{Float64}, Vector{Float64}}([6.25, 1.5833333333333333], [0.0, 0.5])

Periodogram{T, F<:Union{Frequencies,AbstractRange}, V<:AbstractVector{T}} <: TFR{T}

A Periodogram object with fields:

• power::V
• freq::F.

See power and freq for further details.

Periodogram2{T, F1<:Union{Frequencies,AbstractRange}, F2<:Union{Frequencies,AbstractRange}, M<:AbstractMatrix{T}} <: TFR{T}

A Periodogram2 object with fields:

• power::M
• freq1::F1
• freq2::F2

See power and freq for further details.

welch_pgram(s::AbstractVector, n=div(length(s), 8), noverlap=div(n, 2); onesided=eltype(s)<:Real,
            nfft=nextfastfft(n), fs=1, window=nothing)

Computes the Welch periodogram of a signal s based on segments with n samples with overlap of noverlap samples, and returns a Periodogram object.

welch_pgram by default divides s into 8 segments with 50% overlap between them. For a Bartlett periodogram, set noverlap=0.

See periodogram for description of optional keyword arguments.

Examples

julia> x = rect(10; padding=20);

julia> power(welch_pgram(x))                    # 1-sided periodogram
2-element Vector{Float64}:
 0.9523809523809523
 0.04761904761904761

julia> power(welch_pgram(x; onesided=false))    # 2-sided periodogram
3-element Vector{Float64}:
 0.9523809523809523
 0.023809523809523805
 0.023809523809523805

julia> power(welch_pgram(x, 5; onesided=false)) # 5 samples segment
5-element Vector{Float64}:
 1.488888888888889
 0.04444444444444444
 0.044444444444444446
 0.044444444444444446
 0.04444444444444444

welch_pgram with $x = sin(2π(100)t) + 2sin(2π(150)t) + noise$

julia> Fs = 1000;                         # Sampling frequency

julia> t = (1:Fs)/Fs;                     # 1000 time samples

julia> f = [100 150];                     # 100Hz & 150Hz frequencies

julia> A = [1; 2];                        # Amplitudes

julia> x = sin.(2π*f.*t)*A + randn(1000); # Noise corrupted x signal

julia> pxx = welch_pgram(x, 100; fs=Fs, nfft=512, window=hamming);

welch_pgram(s::AbstractVector, config::WelchConfig)

Computes the Welch periodogram of the given signal s using a predefined WelchConfig object.

Examples

julia> x = rect(10; padding=20);

julia> wconfig = WelchConfig(x; fs=1000, window=hamming);

julia> welch_pgram(x, wconfig);

welch_pgram!(out::AbstractVector, s::AbstractVector, n=div(length(s), 8),
             noverlap=div(n, 2); onesided=eltype(s)<:Real, nfft=nextfastfft(n),
             fs=1, window=nothing)

Computes the Welch periodogram of a signal s, based on segments with n samples with overlap of noverlap samples, and returns a Periodogram object. The Power Spectral Density (PSD) of Welch periodogram is stored in out.

welch_pgram! by default divides s into 8 segments with 50% overlap between them. For a Bartlett periodogram, set noverlap=0.

See periodogram for description of optional keyword arguments.

Examples

julia> Fs = 100;                     # Sampling frequency

julia> t = (1:Fs)/Fs;                # 100 time samples

julia> x = cos.(2π*25*t);            # Cosine signal with 25Hz frequency

julia> nfft = 512;

julia> y = vec(zeros(1, nfft));

julia> pxx = welch_pgram!(y, x, 10; onesided=false, nfft=nfft, fs=Fs);

julia> y == power(pxx)
true

welch_pgram!(out::AbstractVector, s::AbstractVector, config::WelchConfig)

Computes the Welch periodogram of the given signal s, using a predefined WelchConfig object. The Power Spectral Density (PSD) of Welch periodogram is stored in out.

Examples

julia> Fs = 100;                     # Sampling frequency

julia> t = (1:Fs)/Fs;                # 100 time samples

julia> x = cos.(2π*25*t);            # Cosine signal with 25Hz frequency

julia> nfft = 512;

julia> y = vec(zeros(1, nfft));

julia> wconfig = WelchConfig(x; n=10, onesided=false, nfft=nfft, fs=Fs, window=hamming);

julia> pxx = welch_pgram!(y, x, wconfig);

julia> y == power(pxx)
true

spectrogram(s, n=div(length(s), 8), noverlap=div(n, 2); onesided=eltype(s)<:Real, nfft=nextfastfft(n), fs=1, window=nothing)

Computes the spectrogram of a signal s based on segments with n samples with overlap of noverlap samples, and returns a Spectrogram object.

spectrogram by default divides s into 8 segments with 50% overlap between them. See periodogram for description of optional keyword arguments.

The returned Spectrogram object stores the 3 computed fields: power, freq and time. See power, freq and time for further details.

Examples

julia> Fs = 1000;

julia> t = 0:1/Fs:1-1/Fs;

julia> x = sin.(2π*100*t.*t);

julia> spec = spectrogram(x; fs=Fs);

julia> size(power(spec))
(63, 14)

julia> size(freq(spec))
(63,)

julia> time(spec)
0.0625:0.063:0.8815

Spectrogram{T, F<:Union{Frequencies,AbstractRange}, M<:AbstractMatrix{T}} <: TFR{T}

A Spectrogram object with fields:

• power::M
• freq::F
• time::Float64Range

See power, freq and time for further details.

stft(s, n=div(length(s), 8), noverlap=div(n, 2); onesided=eltype(s)<:Real, nfft=nextfastfft(n), fs=1, window=nothing)

Computes the Short Time Fourier Transform (STFT) of a signal s based on segments with n samples with overlap of noverlap samples, and returns a matrix containing the STFT coefficients.

stft by default divides s into 8 segments with 50% overlap between them. See periodogram for description of optional keyword arguments.

The STFT computes the DFT over K sliding windows (segments) of the signal s. This returns a J x K matrix where J is the number of DFT coefficients and K the number of windowed segments. The kth column of the returned matrix contains the DFT coefficients for the kth segment.

Examples

julia> Fs = 1000;

julia> t = 0:1/Fs:5-1/Fs;

julia> x = sin.(2π*t.*t);

julia> size(stft(x; window=hamming))
(313, 14)

julia> size(stft(x, 500, 250; onesided=false, window=hanning))
(500, 19)

freq(p)

Returns the frequency bin centers for a given Periodogram, Spectrogram, CrossPowerSpectra, or Coherence object.

Returns a tuple of frequency bin centers for a given Periodogram2 object.

Examples

julia> freq(periodogram([1, 2, 3]; fs=210))
2-element AbstractFFTs.Frequencies{Float64}:
  0.0
 70.0

power(p)

For a Periodogram, returns the computed power spectral density (PSD) as a Vector.

For a Spectrogram, returns the computed power spectral density (PSD) at each frequency and time bin as a Matrix. Dimensions are frequency × time.

For a CrossPowerSpectra, returns the pairwise cross power spectral density (CPSD) between each pair of channels at each frequency. Dimensions are channel x channel x frequency.

Examples

julia> power(periodogram([1, 2, 3]; fs=210))
2-element Vector{Float64}:
 0.05714285714285714
 0.009523809523809525

time(p)

Returns the time bin centers for a given Spectrogram object.

Examples

julia> time(spectrogram(0:1/20:1; fs=8000))
0.000125:0.000125:0.0025

coherence(c::Coherence)

Given an Coherence object, returns an n_channels x n_channels x length(freq(c)) array consisting of the pairwise coherences between each channel for each frequency.

Coherence{T, F, A<:AbstractArray{T, 3}}

Holds an n_channels x n_channels x length(freq) array consisting of the pairwise coherences between each channel for each frequency which is accessed by coherence, as well as the frequencies accessed by freq.

See also mt_coherence and mt_coherence!.

mt_pgram(s; onesided=eltype(s)<:Real, nfft=nextfastfft(n), fs=1, nw=4, ntapers=iceil(2nw)-1, window=dpss(length(s), nw, ntapers))
mt_pgram(signal::AbstractVector, config::MTConfig)

Computes the multitaper periodogram of a signal s.

If window is not specified, the signal is tapered with ntapers discrete prolate spheroidal sequences with time-bandwidth product nw. Each sequence is equally weighted; adaptive multitaper is not (yet) supported.

If window is specified, each column is applied as a taper. The sum of periodograms is normalized by the total sum of squares of window.

Returns a Periodogram.

See also mt_pgram! and MTConfig.

mt_pgram!(output, s::AbstractVector{T}; onesided::Bool=eltype(s)<:Real,
    nfft::Int=nextfastfft(length(s)), fs::Real=1,
    nw::Real=4, ntapers::Int=ceil(Int, 2nw)-1,
    window::Union{AbstractMatrix,Nothing}=nothing) where T<:Number
mt_pgram!(output::AbstractVector, signal::AbstractVector, config::MTConfig) -> Periodogram

Computes a multitapered periodogram, storing the output in output. Arguments:

• signal::AbstractVector: should be of length config.n_samples
• output::AbstractVector: should be of length length(config.freq)

Optionally pass an MTConfig object to preallocate temporary variables and choose configuration settings; otherwise, keyword arguments may be passed to choose those settings.

Returns a Periodogram.

See also mt_pgram and MTConfig.

mt_spectrogram(signal::AbstractVector{T}, n::Int=length(s) >> 3,
                              n_overlap::Int=n >> 1; fs=1,
                              onesided::Bool=T <: Real, kwargs...) where {T}
mt_spectrogram(signal::AbstractVector, config::MTSpectrogramConfig)

Compute a multitaper spectrogram, returning a Spectrogram object. Optionally pass a MTSpectrogramConfig object; otherwise, any additional keyword arguments accepted by MTConfig may be passed to configure the tapering.

Returns a Spectrogram.

See also mt_spectrogram!.

mt_spectrogram!(output, signal::AbstractVector{T}, n::Int=length(signal) >> 3,
    n_overlap::Int=n >> 1; fs=1, onesided::Bool=T <: Real, kwargs...) where {T}
mt_spectrogram!(destination::AbstractMatrix, signal::AbstractVector, config::MTSpectrogramConfig)

Computes a multitaper spectrogram using the parameters specified in config. Arguments:

• destination: length(config.mt_config.freq) x length(config.time) matrix. This can be created by DSP.allocate_output(config).
• signal: vector of length config.n_samples
• config: optionally, pass an MTSpectrogramConfig object to hold temporary variables and configuration settings. Otherwise, settings arguments may be passed directly.

Returns a Spectrogram.

See also mt_spectrogram.

mt_cross_power_spectra(signal::AbstractMatrix{T}; fs=1, kwargs...) where {T}
mt_cross_power_spectra(signal::AbstractMatrix, config::MTCrossSpectraConfig)

Computes multitapered cross power spectra between channels of a signal. Arguments:

• signal: n_channels x n_samples
• Optionally pass an MTCrossSpectraConfig object to preallocate temporary variables

and choose configuration settings. Otherwise, any keyword arguments accepted by MTCrossSpectraConfig may be passed here.

Produces a CrossPowerSpectra object holding the n_channels x n_channels x n_frequencies output array (accessed by power) and the corresponding frequencies (accessed by freq).

See also mt_cross_power_spectra! and MTCrossSpectraConfig.

mt_cross_power_spectra!(output, signal::AbstractMatrix; fs=1, kwargs...)
mt_cross_power_spectra!(output, signal::AbstractMatrix, config::MTCrossSpectraConfig)

Computes multitapered cross power spectra between channels of a signal. Arguments:

• output: n_channels x n_channels x length(config.freq). Can be created by DSP.allocate_output(config).
• signal: n_channels x n_samples
• config: MTCrossSpectraConfig{T}: optionally pass a MTCrossSpectraConfig to preallocate temporary and choose configuration settings. Otherwise, one may pass any keyword arguments accepted by this object.

Produces a CrossPowerSpectra object holding the n_channels x n_channels x n_frequencies output array and the corresponding frequencies (accessed by freq).

See also mt_cross_power_spectra and MTCrossSpectraConfig.

CrossPowerSpectra{T,F,A<:AbstractArray{T, 3}}

Access the cross power spectral density (an n_channels x n_channels x length(freq) array) via the function power, and the frequencies by the function freq.

See also mt_cross_power_spectra and mt_cross_power_spectra!.

mt_coherence(signal::AbstractMatrix{T}; fs=1, freq_range = nothing, demean=false, kwargs...) where T
mt_coherence(signal::AbstractMatrix, config::MTCoherenceConfig)

Arguments:

• signal: n_channels x n_samples matrix
• Optionally pass an MTCoherenceConfig to pre-allocate temporary variables and choose configuration settings, otherwise, see MTCrossSpectraConfig for the meaning of the keyword arguments.

Returns a Coherence object.

See also mt_coherence and MTCoherenceConfig.

mt_coherence!(output, signal::AbstractMatrix; fs=1, freq_range=nothing, demean=false, kwargs...)
mt_coherence!(output, signal::AbstractMatrix, config::MTCoherenceConfig)

Computes the pairwise coherences between channels.

• output: n_channels x n_channels matrix
• signal: n_samples x n_channels matrix
• config: optional configuration object that pre-allocates temporary variables and choose settings.

Returns a Coherence object.

See also mt_coherence and MTCoherenceConfig.

WelchConfig(s::AbstractArray; n=size(s, ndims(s))>>3, noverlap=n>>1,
         onesided=eltype(s)<:Real, nfft=nextfastfft(n),
         fs=1, window=nothing)

WelchConfig(nsamples, eltype; n=nsamples>>3, noverlap=n>>1,
         onesided=eltype<:Real, nfft=nextfastfft(n),
         fs=1, window=nothing)

Captures all configuration options for welch_pgram in a single struct (akin to MTConfig). When passed on the second argument of welch_pgram, computes the periodogram based on segments with n samples with overlap of noverlap samples, and returns a Periodogram object. For a Bartlett periodogram, set noverlap=0. See periodogram for description of optional keyword arguments.

Examples

julia> x = 1:8;

julia> wconfig1 = WelchConfig(x; n=length(x), noverlap=0, window=hanning, nfft=16);

julia> wconfig2 = WelchConfig(8, Float64; n=length(x), noverlap=0, window=hanning, nfft=16);

MTConfig{T}(n_samples; fs=1,
        nfft = nextpow(2, n_samples),
        window = nothing,
        nw = 4,
        ntapers = 2 * nw - 1,
        taper_weights = fill(1/ntapers, ntapers),
        onesided::Bool=T<:Real,
        fft_flags = FFTW.MEASURE)

Creates a config object which holds the configuration state and temporary variables used in multitaper computations, e.g. mt_pgram!, mt_spectrogram, MTSpectrogramConfig, MTCrossSpectraConfig, and MTCoherenceConfig.

An MTConfig can be re-used between computations as long as none of the input arguments change.

• n_samples: the number of samples to be used as input when computing multitaper periodograms with this configuration. Used for pre-allocating temporary buffers.
• fs: the number of samples per second of the input signal
• nfft: length of input vector to the FFT; if nfft > n_samples, then the input signal will be zero-padded until it is of length nfft.
• window: window function to use for tapering. If left at the default of nothing, window will be set to dpss(n_samples, nw, ntapers).
• ntapers: the number of tapers to use.
• taper_weights = fill(1/ntapers, ntapers): how to weight the contribution of each taper. The default setting is to simply average them.
• onesided: whether or not to compute a "one-sided" FFT by using that real signal data yields conjugate-symmetry in Fourier space.
• fft_flags: flags to control how the FFT plan is generated.

MTSpectrogramConfig(n_samples, mt_config::MTConfig{T}, n_overlap_samples) where {T}
MTSpectrogramConfig{T}(n_samples, samples_per_window, n_overlap_samples; fs=1, kwargs...) where {T}

Creates a MTSpectrogramConfig which holds configuration and temporary variables for mt_spectrogram and mt_spectrogram!. Any keyword arguments accepted by MTConfig may be passed here, or an MTConfig object itself.

MTCrossSpectraConfig{T}(n_channels, n_samples; fs=1, demean=false, freq_range=nothing,
                        ensure_aligned = T == Float32 || T == Complex{Float32}, kwargs...) where {T}
MTCrossSpectraConfig(n_channels, mt_config::MTConfig{T}; demean=false, freq_range=nothing,
                     ensure_aligned = T == Float32 || T == Complex{Float32})

Creates a configuration object used for mt_cross_power_spectra and mt_cross_power_spectra!.

• n_channels: the number of channels of the input.
• n_samples: the number of samples for each channel of the input.
• demean: if true, the channelwise mean will be subtracted from the input signals before the cross spectral powers are computed.
• freq_range: if nothing, all frequencies are retained. Otherwise, only frequencies between first(freq_range) and last(freq_range) are retained.
• ensure_aligned = T == Float32 || T == Complex{Float32}: perform an extra copy to ensure that the FFT output is memory-aligned
• Additionally, either pass an MTConfig object, or keyword arguments such as fs accepted by MTConfig.

Returns a CrossPowerSpectra object.

MTCoherenceConfig{T}(n_channels, n_samples; fs=1, demean=false, freq_range=nothing, kwargs...) where T
MTCoherenceConfig(cs_config::MTCrossSpectraConfig{T}) where {T}
MTCoherenceConfig(n_channels, mt_config::MTConfig{T}; demean=false, freq_range=nothing,
    ensure_aligned = T == Float32 || T == Complex{Float32}) where {T}

Creates a configuration object for coherences from a MTCrossSpectraConfig. Provides a helper method with the same arugments as MTCrossSpectraConfig to construct the MTCrossSpectraConfig object.

See also mt_coherence and mt_coherence!.