Wavelets.jl icon indicating copy to clipboard operation
Wavelets.jl copied to clipboard

multilevel wavelet coefficients

Open tapir27 opened this issue 3 years ago • 3 comments

I am sure I am missing something but I when I run dwt I can only see the first level wavelet detail coefficients. How can I see the other levels?

tapir27 avatar Sep 28 '20 11:09 tapir27

The return value is a vector of scaling and detail coefficients. Where they are depends on the number of levels specified. There are util functions that can be used to determine the ranges of coefficients.

gummif avatar Sep 28 '20 16:09 gummif

I've been using this utility wrapper:

# Define DyadicLevels; a useful way to iterate over a dyadic signal representation,
# such as that returned from `dwt()`.
struct DyadicLevels
    data::Vector
    max_level::Int

    function DyadicLevels(data::Vector, max_level::Union{Int64,Nothing} = nothing)
        max_possible_level = Wavelets.maxtransformlevels(length(data))
        if max_level === nothing
            max_level = max_possible_level
        elseif max_level > max_possible_level
            throw(ArgumentError("Cannot decompose $(max_level) levels with an input of length $(length(data))"))
        end
        return new(data, max_level)
    end
end

function scalingrange(l::Int, max_level::Int)
    return 1:(l - sum(length(detailrange(l, level)) for level in 1:max_level))
end
scalingrange(x::Vector, max_level::Int) = scalingrange(length(x), max_level)

# Allow easy iteration over the wavelet levels
Base.length(ds::DyadicLevels) = min(Wavelets.maxtransformlevels(length(ds.data)), ds.max_level) + 1
function Base.iterate(ds::DyadicLevels, state=(1,))
    scale_idx = state[1]
    if scale_idx > length(ds)
        return nothing
    end

    next_state = (scale_idx + 1,)
    # Special-case the scaling coefficients as the lowest level:
    if scale_idx == 1
        return ds.data[scalingrange(ds.data, ds.max_level)], next_state
    end

    dyadic_slice = ds.data[detailrange(ds.data, length(ds) - (scale_idx - 1))]
    return dyadic_slice, next_state
end

You can use it with a function like the following, to do a quick-n-dirty plot of your wavelet levels:

function plot_wavelet_levels(x, max_level)
    fig = WGLMakie.Figure()
    ax = WGLMakie.Axis(fig[1,1])

    lines!(ax, x; label="Signal", color="black")
    X = dwt(x, wavelet(WT.db2, WT.Filter), max_level)
    for (level, X_range) in enumerate(DyadicLevels(X, max_level))
        t = (0:(length(X_range)-1))*(length(X)-1)/length(X_range)
        if level == 1
            lines!(ax, t, X_range./4; label="Scaling")
        else
            lines!(ax, t, X_range; label="Detail $(level)")
        end
    end
    
    fig[1,2] = Legend(fig, ax, "Wavelet levels")
    display(fig)
    return nothing
end

staticfloat avatar Sep 16 '21 22:09 staticfloat

See also utility functions in https://github.com/JuliaDSP/Wavelets.jl/blob/master/src/mod/Util.jl that may be useful.

gummif avatar Sep 17 '21 10:09 gummif