Wavelet filter definitions

[ElOceanografo]

Implementing #42, I noticed filter definitions are different in this package than in at least one other wavelet package I've used (wmtsa in R). The code below plots the equivalently-named wavelets from these two packages.

using Wavelets
using Plots
using RCall; R"library(wmtsa)"

julia_wavelets = [:coif6, :db2, :db4, :db6, :db8, :db10, :haar, :sym4, :sym6, :sym8]
r_wavelets = ["c6", "d2", "d4", "d6", "d8", "d10", "haar", "s4", "s6", "s8"]

subplots = []
for (wj, wr) in zip(julia_wavelets, r_wavelets)
    hr = Array(reval("wavDaubechies('$wr')")[:wavelet])
    gr = Array(reval("wavDaubechies('$wr')")[:scaling])
    p = plot(hr, color=:pink, label="Wavelet (R)", 
        title="Julia: $(String(wj)), R: $wr", 
        background_color_legend=:transparent, size=(600, 1000))
    plot!(p, gr, color=:red, label="Scaling (R)")

    g, h = @eval WT.makeqmfpair(wavelet(WT.$wj))
    plot!(p, h, color=:lightblue, label="Wavelet (Julia)")
    plot!(p, g, color=:blue, label="Scaling (Julia)")
    push!(subplots, p)
end
plot(subplots..., layout=(5,2))

This may just be a Pepsi/Coke problem with slightly different definitions, but it seemed worth double-checking.

[gugushvili]

I guess it's okay. But it's worth indicating the source of Julia definitions.

P.S. It looks like the wavethresh package in R uses slightly different definitions than wmtsa (and other packages depending on Percival and Walden, like waveslim and wavelets).

[gugushvili]

Figures indicate that filters with the same names in Julia and wmtsa don't have the same width (they differ by factor 2, e.g. 8 vs. 16). Thus, for instance, sym4 in Julia must correspond to s8 in R. Also, by default, wavDaubechies() scales filters with sqrt(2) (can be turned off by including the option normalized=FALSE).

[gummif]

The numbers of the wavelet types correspond to the number of vanishing moments of the wavelets. I think it is consistent among all the orthogonal and biorthogonal wavelet types that are parametrized. The filters are also all normalized in l2-norm, defining orthonormal bases.

[gugushvili]

That numbers were referring to vanishing moments was clear, there are authors who use that convention.

I checked it for the sym4/s8 filter, and absolute magnitudes of the filter coefficients do agree in Julia and R, but signs don't. It would be useful to document all such differences, and to actually write down maths behind the Julia implementation (for R it's Percival & Walden).

[gummif]

I have begun writing down the mathematics (based on Mallat). I then need to verify that it corresponds to the actual implementation.

Something to think about is whether it makes sense to allow the user to specify the filter convention to use, for possible consistency and interoperability with other packages.