JuliaDSP / Wavelets.jl

A Julia package for fast discrete wavelet transforms and utilities
Other
191 stars 29 forks source link

Document coefficient vector layout #37

Open mschauer opened 6 years ago

mschauer commented 6 years ago

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.

gummif commented 6 years ago

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.

gugushvili commented 6 years ago

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/

mschauer commented 6 years ago

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 commented 6 years ago

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