From 9d3687f1419552cae10f33ae19c96ca765dc5e29 Mon Sep 17 00:00:00 2001 From: Chris Foster Date: Fri, 17 Feb 2017 17:02:45 +1000 Subject: [PATCH] Window Welch periodogram with Hann window by default. This choice is used by scipy, and seems less likely than the square window to produce completely wrong results for novice users. matlab and octave use the Hamming window, but this works very poorly for signals with a nontrivial DC component when `n != nfft`. --- doc/periodograms.rst | 27 +++++++++++++++++++++------ src/periodograms.jl | 2 +- 2 files changed, 22 insertions(+), 7 deletions(-) diff --git a/doc/periodograms.rst b/doc/periodograms.rst index 19fdd2d3d..1c1b46648 100644 --- a/doc/periodograms.rst +++ b/doc/periodograms.rst @@ -28,12 +28,27 @@ equal to the uncentered variance (or average power) of the original signal. -.. function:: welch_pgram(s, 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. For a Bartlett periodogram, set ``noverlap=0``. See - :func:`periodogram` for description of optional keyword arguments. +.. function:: welch_pgram(s, n=div(length(s), 8), noverlap=div(n, 2); onesided=eltype(s)<:Real, nfft=nextfastfft(n), fs=1, window=hanning) + + Computes the Welch periodogram of a signal ``s`` by dividing the signal + into overlapping segments of length ``n``, averaging the periodograms of + each segment, and returning a ``Periodogram`` object. ``noverlap`` is the + number of samples of overlap, with a new segment starting every + ``n-noverlap`` samples of ``s``. Each segment is padded with zeros to a + length of ``nfft`` for efficiency. + + By default, the signal is windowed using the hanning window as a moderate + tradeoff between high spectral resolution (the square window) and + sensitivity to signals near the noise floor. A window which falls to zero + also reduces the spurious ringing which occurs when ``n!=nfft`` for signals + with large mean. + + Compared to ``periodogram()``, ``welch_pgram()`` reduces noise in the + estimated power spectral density as ``length(s)`` increases. In the basic + periodogram, we increase the frequency resolution instead. For a Bartlett + periodogram, set ``noverlap=0``. + + See :func:`periodogram` for a description of the other keyword arguments. .. function:: mt_pgram(s; onesided=eltype(s)<:Real, nfft=nextfastfft(n), fs=1, nw=4, ntapers=iceil(2nw)-1, window=dpss(length(s), nw, ntapers)) diff --git a/src/periodograms.jl b/src/periodograms.jl index 0875c8e57..84fdd2184 100644 --- a/src/periodograms.jl +++ b/src/periodograms.jl @@ -281,7 +281,7 @@ forward_plan{T<:Union{Complex64, Complex128}}(X::AbstractArray{T}, Y::AbstractAr function welch_pgram{T<:Number}(s::AbstractVector{T}, n::Int=length(s)>>3, noverlap::Int=n>>1; onesided::Bool=eltype(s)<:Real, nfft::Int=nextfastfft(n), fs::Real=1, - window::Union{Function,AbstractVector,Void}=nothing) + window::Union{Function,AbstractVector,Void}=hanning) onesided && T <: Complex && error("cannot compute one-sided FFT of a complex signal") nfft >= n || error("nfft must be >= n")