Wavelets.jl
Wavelets.jl copied to clipboard
Document coefficient vector layout
I a hope I am not missing anything but I could not find the layout of the coefficient vector. Inspecting the size of the coefficient also does not immediately suggest a certain order to me.
You are right, this is not well documented. There are some utiliities in https://github.com/JuliaDSP/Wavelets.jl/blob/master/src/mod/Util.jl , especially the dyadicdetailrange
and dyadicscalingrange
functions (assuming the signal is dyadic, a power of 2). There is also detailrange
for other more general sizes of signals, but no scalingrange
function. These give you the range of the vector that corresponds to each level of the coefficients. I will take a look at adding missing functions and documenting these better.
What would be very helpful is to have equivalents of the functions accessC, accessD, putC and putD from the wavethresh package in R. These allow to access or replace smooth and detail coefficients in entire resolution levels of DWT. See https://cran.r-project.org/web/packages/wavethresh/
Ah, so my first guess for the memory layout was actually right according to
dyadicdetailrange(j::Integer) = (2^j+1):(2^(j+1))
I just did not believe it: i started with smooth signal of length 128 and expected the coefficients xt[65:128]
small and of the same order, but that was not the case for the example I looked at:
julia> t = linspace(0,1,128)
0.0:0.007874015748031496:1.0
julia> x = sin.(2pi*t)
julia> using Wavelets
julia> xt = dwt(x, wavelet(WT.db4, WT.Periodic))
128-element Array{Float64,1}:
-5.48173e-16
5.31657
-4.19097
4.11042
-0.290043
0.489353
0.33211
-0.574409
0.0127308
0.00217123
⋮
-2.50345e-6
-2.30573e-6
-2.08546e-6
-1.84479e-6
-1.58608e-6
-1.31184e-6
-1.02478e-6
-7.27685e-7
-4.23475e-7
julia> xt[65:80]
16-element Array{Float64,1}:
0.0239865
-0.00864036
0.00214933
6.5407e-7
9.53168e-7
1.24294e-6
1.52056e-6
1.78329e-6
2.02859e-6
2.25404e-6
2.45744e-6
2.63679e-6
2.79036e-6
2.91662e-6
3.01436e-6
3.0826e-6
Maybe these are boundary effects?
Ah, I found my mistake, the time points t = linspace(0,1,128)
were not at equal distance on a circle of length 1
, this works.
julia> t = (1:128)./128
julia> x = sin.(2pi*t)
julia> xt = dwt(x, wavelet(WT.db4, WT.Periodic))
julia> xt[65:80]
16-element Array{Float64,1}:
-1.10697e-7
1.86888e-7
4.82674e-7
7.73811e-7
1.0575e-6
1.331e-6
1.59168e-6
1.83703e-6
2.06469e-6
2.27247e-6
2.45836e-6
2.62058e-6
2.75756e-6
2.86798e-6
2.95079e-6
3.00517e-6