Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Deprecate welch_pgram without explicitly given window #581

Merged
merged 2 commits into from
Nov 14, 2024
Merged

Deprecate welch_pgram without explicitly given window #581

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6284117
Deprecate `welch_pgram` without explicitly given `window`
martinholters Nov 7, 2024
eeb56d5
Test `welch_pgram` deprecation
martinholters Nov 7, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 20 additions & 12 deletions src/periodograms.jl
View file
Open in desktop
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

Comment on lines +610 to 620
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Due to the short length of x, the default n==3 and

julia> hanning(3)
3-element Vector{Float64}:
 0.0
 1.0
 0.0

is not at all better than the default rect window. In the next example, where n=5 is explicitly set, I've changed the window to hanning as it somewhat reasonable there. I think from a pedagogical viewpoint this is ok, and I'd like to keep a complete rewrite of the example out of the scope of this PR.

Sidenote: a variant of the windows without the zeros at the edges might a useful here, e.g.

julia> hanning(5)[begin+1:end-1]
3-element Vector{Float64}:
 0.5
 1.0
 0.5

But that's also beyond the scope of this PR.

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
Open in desktop
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
Loading