Skip to content

Commit

Permalink
Deprecate welch_pgram without explicitly given window (#581)
Browse files Browse the repository at this point in the history
  • Loading branch information
@martinholters
martinholters authored Nov 14, 2024
1 parent 3c1b783 commit 8902bef
Show file tree
Hide file tree
Showing 2 changed files with 34 additions and 19 deletions.
32 changes: 20 additions & 12 deletions src/periodograms.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -530,11 +530,11 @@ end
"""
WelchConfig(s::AbstractArray; n=size(s, ndims(s))>>3, noverlap=n>>1,
onesided=eltype(s)<:Real, nfft=nextfastfft(n),
fs=1, window=nothing)
fs=1, window)
WelchConfig(nsamples, eltype; n=nsamples>>3, noverlap=n>>1,
onesided=eltype<:Real, nfft=nextfastfft(n),
fs=1, window=nothing)
fs=1, window)
Captures all configuration options for [`welch_pgram`](@ref) in a single struct (akin to
[`MTConfig`](@ref)). When passed on the second argument of [`welch_pgram`](@ref), computes the
Expand All @@ -559,7 +559,7 @@ julia> wconfig2 = WelchConfig(8, Float64; n=length(x), noverlap=0, window=hannin
"""
function WelchConfig(nsamples, ::Type{T}; n::Int=nsamples >> 3, noverlap::Int=n >> 1,
onesided::Bool=T <: Real, nfft::Int=nextfastfft(n),
fs::Real=1, window::Union{Function,AbstractVector,Nothing}=nothing) where T
fs::Real=1, window::Union{Function,AbstractVector,Nothing}=welch_config_default_window()) where T

onesided && T <: Complex && throw(ArgumentError("cannot compute one-sided FFT of a complex signal"))
nfft >= n || throw(DomainError((; nfft, n), "nfft must be >= n"))
Expand All @@ -579,6 +579,12 @@ function WelchConfig(data::AbstractArray; kwargs...)
return WelchConfig(size(data, ndims(data)), eltype(data); kwargs...)
end

function welch_config_default_window()
# deprecation added in 0.8, TODO for 0.9: replace with just `hanning`
Base.depwarn("Omitting `window` is deprecated; specify `window=nothing` for the old behaviour or `window=hanning` for the future default.", :WelchConfig)
return nothing
end

# Compute an estimate of the power spectral density of a signal s via Welch's
# method. The resulting periodogram has length N and is computed with an overlap
# region of length M. The method is detailed in "The Use of Fast Fourier Transform
Expand All @@ -587,7 +593,7 @@ end
# vol AU-15, pp 70-73, 1967.
"""
welch_pgram(s::AbstractVector, n=div(length(s), 8), noverlap=div(n, 2); onesided=eltype(s)<:Real,
nfft=nextfastfft(n), fs=1, window=nothing)
nfft=nextfastfft(n), fs=1, window)
Computes the Welch periodogram of a signal `s` based on segments with `n` samples
with overlap of `noverlap` samples, and returns a [`Periodogram`](@ref) object.
Expand All @@ -601,24 +607,26 @@ See [`periodogram`](@ref) for description of optional keyword arguments.
```jldoctest
julia> x = rect(10; padding=20);
julia> power(welch_pgram(x)) # 1-sided periodogram
julia> power(welch_pgram(x; window=nothing)) # 1-sided periodogram
2-element Vector{Float64}:
0.9523809523809523
0.04761904761904761
julia> power(welch_pgram(x; onesided=false)) # 2-sided periodogram
julia> power(welch_pgram(x; onesided=false, window=nothing)) # 2-sided periodogram
3-element Vector{Float64}:
0.9523809523809523
0.023809523809523805
0.023809523809523805
julia> power(welch_pgram(x, 5; onesided=false)) # 5 samples segment
julia> power(welch_pgram(x, 5; onesided=false, window=hanning)) # 5 samples segment
5-element Vector{Float64}:
1.488888888888889
0.04444444444444444
0.044444444444444446
0.044444444444444446
0.04444444444444444
0.8888888888888888
0.3807834425694269
0.008105446319461968
0.008105446319461968
0.3807834425694269
```
welch_pgram with ``x = sin(2π(100)t) + 2sin(2π(150)t) + noise ``
```jldoctest
Expand Down
21 changes: 14 additions & 7 deletions test/periodograms.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,31 +106,37 @@ end
4.0,
13.656854249492380]
@test power(periodogram(data, onesided=false)) ≈ data0
@test power(welch_pgram(data, length(data), 0, onesided=false)) ≈ data0
@test power(@test_deprecated welch_pgram(data, length(data), 0, onesided=false)) ≈ data0
@test power(welch_pgram(data, length(data), 0, onesided=false, window=nothing)) ≈ data0
@test power(spectrogram(data, length(data), 0, onesided=false)) ≈ data0
@test power(periodogram(complex.([data;], [data;]), onesided=false)) ≈ data0*2
@test power(welch_pgram(complex.([data;], [data;]), length(data), 0, onesided=false)) ≈ data0*2
@test power(@test_deprecated welch_pgram(complex.([data;], [data;]), length(data), 0, onesided=false)) ≈ data0*2
@test power(welch_pgram(complex.([data;], [data;]), length(data), 0, onesided=false, window=nothing)) ≈ data0*2
@test power(spectrogram(complex.([data;], [data;]), length(data), 0, onesided=false)) ≈ data0*2

# # ~~~~~~~~ Tests with no window ~~~~~~~~~~~~~~~~~~~
# Matlab: p = pwelch(0:7, [1, 1], 0, 2, 1, 'twosided')
expected = Float64[34.5, 0.5]
@test power(welch_pgram(data, 2, 0; onesided=false)) ≈ expected
@test power(@test_deprecated welch_pgram(data, 2, 0; onesided=false)) ≈ expected
@test power(welch_pgram(data, 2, 0; onesided=false, window=nothing)) ≈ expected
@test mean(power(spectrogram(data, 2, 0; onesided=false)), dims=2) ≈ expected

# Matlab: p = pwelch(0:7, [1, 1, 1], 0, 3, 1, 'twosided')
expected = Float64[25.5, 1.0, 1.0]
@test power(welch_pgram(data, 3, 0; onesided=false)) ≈ expected
@test power(@test_deprecated welch_pgram(data, 3, 0; onesided=false)) ≈ expected
@test power(welch_pgram(data, 3, 0; onesided=false, window=nothing)) ≈ expected
@test mean(power(spectrogram(data, 3, 0; onesided=false)), dims=2) ≈ expected

# Matlab: p = pwelch(0:7, [1, 1, 1], 1, 3, 1, 'twosided')
expected = Float64[35.0, 1.0, 1.0]
@test power(welch_pgram(data, 3, 1; onesided=false)) ≈ expected
@test power(@test_deprecated welch_pgram(data, 3, 1; onesided=false)) ≈ expected
@test power(welch_pgram(data, 3, 1; onesided=false, window=nothing)) ≈ expected
@test mean(power(spectrogram(data, 3, 1; onesided=false)), dims=2) ≈ expected

# Matlab: p = pwelch(0:7, [1, 1, 1, 1], 1, 4, 1, 'twosided')
expected = Float64[45, 2, 1, 2]
@test power(welch_pgram(data, 4, 1; onesided=false)) ≈ expected
@test power(@test_deprecated welch_pgram(data, 4, 1; onesided=false)) ≈ expected
@test power(welch_pgram(data, 4, 1; onesided=false, window=nothing)) ≈ expected
@test mean(power(spectrogram(data, 4, 1; onesided=false)), dims=2) ≈ expected

# ~~~~~~~~~~~ This one tests periodogram ~~~~~~~~~~~~
Expand Down Expand Up @@ -189,7 +195,8 @@ end
2
]
@test power(periodogram(data; nfft=32)) ≈ expected
@test power(welch_pgram(data, length(data), 0; nfft=32)) ≈ expected
@test power(@test_deprecated welch_pgram(data, length(data), 0; nfft=32)) ≈ expected
@test power(welch_pgram(data, length(data), 0; nfft=32, window=nothing)) ≈ expected
@test power(spectrogram(data, length(data), 0; nfft=32)) ≈ expected

# Padded periodogram with window
Expand Down

0 comments on commit 8902bef

Please sign in to comment.