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

Document coefficient vector layout

Open mschauer opened this issue 6 years ago • 4 comments

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.

mschauer avatar May 30 '18 13:05 mschauer

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.

gummif avatar May 30 '18 21:05 gummif

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/

gugushvili avatar May 31 '18 08:05 gugushvili

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?

mschauer avatar May 31 '18 08:05 mschauer

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

mschauer avatar Jun 01 '18 09:06 mschauer