JuliaDSP / DSP.jl

Filter design, periodograms, window functions, and other digital signal processing functionality
https://docs.juliadsp.org/stable/contents/
Other
381 stars 109 forks source link

Fix `welch_pgram!`, update periodogram docs #557

Closed kcin96 closed 4 months ago

kcin96 commented 5 months ago

This PR

  1. Update docs to mention that power spectral densities are computed for periodograms as per #480
  2. Update periodogram docs with examples.
  3. Fixes welch_pgram! so Floats can work as well as Ints.
  4. Add test to welch_pgram!
codecov[bot] commented 5 months ago

Codecov Report

All modified and coverable lines are covered by tests :white_check_mark:

Project coverage is 97.74%. Comparing base (67d62c9) to head (374fc1d).

Additional details and impacted files ```diff @@ Coverage Diff @@ ## master #557 +/- ## ======================================= Coverage 97.74% 97.74% ======================================= Files 18 18 Lines 3197 3199 +2 ======================================= + Hits 3125 3127 +2 Misses 72 72 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

wheeheee commented 5 months ago

Sorry, I don't really see how floats are not currently accepted as inputs to welch_pgram!. AFAICT your change relaxes the type constraints on s, which seems wrong to me since I feel if someone's using a config they should take care to match it to their input. Care to explain?

kcin96 commented 5 months ago

Currently, when I run this block of code to test for Floats, I get a method error.

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);
ERROR: MethodError: no method matching welch_pgram!(::Vector{Float64}, ::Vector{Float64}, ::WelchConfig{Int64, AbstractFFTs.Frequencies{Float64}, Vector{Float64}, FFTW.rFFTWPlan{Float64, -1, false, 1, Tuple{Int64}}, Vector{Float64}, Vector{ComplexF64}, Float64})
Closest candidates are:
  welch_pgram!(::AbstractVector, ::AbstractVector, ::Int64) 
  welch_pgram!(::AbstractVector, ::AbstractVector, ::Int64, ::Int64; kwargs...) 
  welch_pgram!(::AbstractVector, ::AbstractVector{T}, ::WelchConfig{T}) where T<:Number
  ...
Stacktrace:
 [1] top-level scope
   @ REPL[24]:1

Whereas, if I next run the welch_pgram(s::AbstractVector, config::WelchConfig) function, I get a result.

julia> welch_pgram(x, wconfig)
DSP.Periodograms.Periodogram{Float64, AbstractFFTs.Frequencies{Float64}, Vector{Float64}}([2.656596353078254e-7, 2.8336305486805697e-7, 3.3679303813800914e-7, 4.269139053267828e-7, 5.553500731514562e-7, 7.244120931838832e-7, 9.371335271002283e-7, 1.1973190664269316e-6, 1.5096044010914354e-6, 1.8795284233756564e-6  …  2.313618418942876e-6, 1.8795284233756564e-6, 1.5096044010914354e-6, 1.1973190664269316e-6, 9.371335271002283e-7, 7.244120931838832e-7, 5.553500731514562e-7, 4.269139053267828e-7, 3.3679303813800914e-7, 2.8336305486805697e-7], [0.0, 0.1953125, 0.390625, 0.5859375, 0.78125, 0.9765625, 1.171875, 1.3671875, 1.5625, 1.7578125  …  -1.953125, -1.7578125, -1.5625, -1.3671875, -1.171875, -0.9765625, -0.78125, -0.5859375, -0.390625, -0.1953125])

Not too sure if this is the intent?

wheeheee commented 5 months ago

Oh right, yeah, seems like the type parameters are wrong, but I think the fix should also include another check that float(eltype(data)) == T1 or something along those lines

kcin96 commented 5 months ago

That is a good point! So something like this which checks that the type of input s matches that of config.inbuf

function welch_pgram!(out::AbstractVector, s::AbstractVector, config::WelchConfig{T}) where T<:Number
     if length(out) != length(config.freq)
         throw(DimensionMismatch("""Expected `output` to be of length `length(config.freq)`;
             got `length(output)` = $(length(out)) and `length(config.freq)` = $(length(config.freq))"""))
    elseif eltype(out) != fftabs2type(T)
        throw(ArgumentError("Eltype of output ($(eltype(out))) doesn't match the expected "*
                            "type: $(fftabs2type(T))."))
    elseif eltype(s) != eltype(config.inbuf)
        throw(ArgumentError("Eltype of input `s` ($(eltype(s))) doesn't match the expected `config` "*
                            "type: $(eltype(config.inbuf))."))
end
wheeheee commented 5 months ago

No, the problem is also that the first type parameter of the WelchConfig is the eltype of fs, which isn't what should be checked. The T in this function should instead just be the type parameter in s::AbstractVector{T}. Then eltype(s) can be replaced by T etc

kcin96 commented 5 months ago

I see, you're right, and it is much cleaner now.