Merge branch 'development'

README.md

@@ -2,7 +2,7 @@
 [![Build Status](https://travis-ci.org/vincentcp/WaveletsCopy.jl.svg?branch=master)](https://travis-ci.org/vincentcp/WaveletsCopy.jl)
 [![Coverage Status](https://coveralls.io/repos/github/vincentcp/WaveletsCopy.jl/badge.svg?branch=master)](https://coveralls.io/github/vincentcp/WaveletsCopy.jl?branch=master)
 # Wavelets
-A julia package for a fast Discrete Wavelet Transform, and plotting and evaluation of wavelets 
+A julia package for a fast Discrete Wavelet Transform, and plotting and evaluation of wavelets
 
 
 ```julia
@@ -12,7 +12,7 @@ A julia package for a fast Discrete Wavelet Transform, and plotting and evaluati
 
 
 ```julia
-using WaveletsCopy 
+using WaveletsCopy
 using Plots
 ```
 
@@ -50,7 +50,7 @@ a = rand(1<<8);
 dwt(a, db3, perbound);
 ```
 
-Only periodic boundaries boundaries and dyadic length are supported. 
+Only periodic boundaries boundaries and dyadic length are supported.
 
 
 ```julia
@@ -166,7 +166,7 @@ plot!()
 
 
 
-Also pointwise evaluation of scaling functions is possible up to some given precision. (There are still errors in the code depending on the point evaluation)
+Also pointwise evaluation of scaling functions is possible up to some given precision.
 `evaluate{T, S<:Real}(side::Side, kind::Scl, w::DiscreteWavelet{T}, j::Int, k::Int, x::S; xtol::S=1e-5)`
 
 `evaluate_periodic{T, S<:Real}(side::Side, kind::Kind, w::DiscreteWavelet{T}, j::Int, k::Int, x::S; xtol::S=1e-5)`
@@ -264,5 +264,3 @@ plot(Dual, wavelet,cdf13)
 
 
 ![svg](README_files/README_33_0.svg)
-
-

notebooks/README.ipynb

@@ -30,7 +30,7 @@
    },
    "outputs": [],
    "source": [
-    "using WaveletsCopys \n",
+    "using WaveletsCopy \n",
     "using Plots"
    ]
   },
@@ -6109,7 +6109,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Also pointwise evaluation of scaling functions is possible up to some given precision. (There are still errors in the code depending on the point evaluation)\n",
+    "Also pointwise evaluation of scaling functions is possible up to some given precision.\n",
     "`evaluate{T, S<:Real}(side::Side, kind::Scl, w::DiscreteWavelet{T}, j::Int, k::Int, x::S; xtol::S=1e-5)`\n",
     "\n",
     "`evaluate_periodic{T, S<:Real}(side::Side, kind::Kind, w::DiscreteWavelet{T}, j::Int, k::Int, x::S; xtol::S=1e-5)`"

src/WaveletsCopy.jl

@@ -1,5 +1,5 @@
 # disable precompilation during development
-#__precompile__()
+__precompile__(true)
 
 
 module WaveletsCopy
@@ -8,7 +8,7 @@ using Reexport
 using CardinalBSplines
 
 include("sequences/sequence.jl")
-include("filterbank.jl")
+include("filterbanks/filterbank.jl")
 include("dwt/discretewavelets.jl")
 
 @reexport using .DWT

src/dwt/discretewavelets.jl

@@ -36,6 +36,9 @@ include("wvlt.jl")
 Filterbank(w::DiscreteWavelet) =
     Filterbank( FilterPair(filter(Prl(), Scl(), w), filter(Prl(), Wvl(), w)),
                 FilterPair(filter(Dul(), Scl(), w), filter(Dul(), Wvl(), w)) )
+DualFilterbank(w::DiscreteWavelet) =
+    Filterbank( FilterPair(filter(Dul(), Scl(), w), filter(Dul(), Wvl(), w)),
+                FilterPair(filter(Prl(), Scl(), w), filter(Prl(), Wvl(), w)) )
 
 "DWT groups the data that fully characterize a discrete wavelet transform."
 struct DWT_Data

src/dwt/dwtstep.jl

@@ -47,7 +47,7 @@ idwtstep_size(sc, dc, fb::Filterbank, bnd::WaveletBoundary) = idwtstep_size(leng
 
 
 function dwtstep(x, fb::Filterbank, bnd::WaveletBoundary)
-    sc_len, dc_len = dwtstep_size(length(x), fb, bnd)
+    sc_len, dc_len = dwtstep_size(x, fb, bnd)
     T = promote_type(eltype(x), eltype(fb))
     sc = zeros(T, sc_len)
     dc = zeros(T, dc_len)
@@ -73,41 +73,42 @@ idwtstep!(x, sc, dc, fb::Filterbank, bnd::PeriodicBoundary) =
 
 
 ## Symmetric boundaries
-
-
-function dwtstep!(sc, dc, x, fb::Filterbank, bnd::SymmetricBoundary)
-    if iseven(length(x))
-        if isodd(sublength(dual_lowpassfilter(fb)))
-            polyphase_analysis!(sc, dc, x, fb.pm_analysis, SymmetricEmbedding{:wp,:wp,:even,:even}())
-        end
-    end
-end
-
-function idwtstep!(x, sc, dc, fb::Filterbank, bnd::SymmetricBoundary)
-    if iseven(length(x))
-        if isodd(sublength(primal_lowpassfilter(fb)))
-            polyphase_synthesis!(x, sc, dc, fb.pm_synthesis, SymmetricEmbedding{:hp,:hp,:even,:even}())
-        end
-    end
-end
-
-
-"Compute a matrix representation of one DWT step for a signal of size n."
-function dwtstep_matrix(n, fb::Filterbank, bnd::WaveletBoundary)
-    T = eltype(fb)
-    A = zeros(n,n)
-    sc_len, dc_len = dwtstep_size(n, fb, bnd)
-    sc = zeros(T, sc_len)
-    dc = zeros(T, dc_len)
-    x = zeros(T, n)
-    for i = 1:n
-        if i > 1
-            x[i-1] = 0
-        end
-        x[i] = 1
-        dwtstep!(sc, dc, x, fb, bnd)
-        A[1:sc_len,i] = sc
-        A[sc_len+1:n,i] = dc
-    end
-    A
-end
+# TODO implement
+#
+#
+# function dwtstep!(sc, dc, x, fb::Filterbank, bnd::SymmetricBoundary)
+#     if iseven(length(x))
+#         if isodd(sublength(dual_lowpassfilter(fb)))
+#             polyphase_analysis!(sc, dc, x, fb.pm_analysis, SymmetricEmbedding{:wp,:wp,:even,:even}())
+#         end
+#     end
+# end
+#
+# function idwtstep!(x, sc, dc, fb::Filterbank, bnd::SymmetricBoundary)
+#     if iseven(length(x))
+#         if isodd(sublength(primal_lowpassfilter(fb)))
+#             polyphase_synthesis!(x, sc, dc, fb.pm_synthesis, SymmetricEmbedding{:hp,:hp,:even,:even}())
+#         end
+#     end
+# end
+#
+#
+# "Compute a matrix representation of one DWT step for a signal of size n."
+# function dwtstep_matrix(n, fb::Filterbank, bnd::WaveletBoundary)
+#     T = eltype(fb)
+#     A = zeros(n,n)
+#     sc_len, dc_len = dwtstep_size(n, fb, bnd)
+#     sc = zeros(T, sc_len)
+#     dc = zeros(T, dc_len)
+#     x = zeros(T, n)
+#     for i = 1:n
+#         if i > 1
+#             x[i-1] = 0
+#         end
+#         x[i] = 1
+#         dwtstep!(sc, dc, x, fb, bnd)
+#         A[1:sc_len,i] = sc
+#         A[sc_len+1:n,i] = dc
+#     end
+#     A
+# end

src/dwt/dwttransform.jl

@@ -9,21 +9,26 @@ end
 function idwt!
 end
 
-dwt_size(x, w::DiscreteWavelet, bnd::WaveletBoundary) = dwt_size(x, Filterbank(w), bnd)
 dwt_size(x, fb::Filterbank, bnd::WaveletBoundary) = length(x)
 
-dwt_dual(x, w::DiscreteWavelet, bnd::WaveletBoundary, L::Int=maxtransformlevels(x)) =
-    dwt(x, dual(Filterbank(w)), bnd, L)
+dwt(x, w::DiscreteWavelet, bnd::WaveletBoundary, L::Int=maxtransformlevels(x)) =
+    dwt(x, Primal, w, bnd, L)
 
-idwt_dual(x, w::DiscreteWavelet, bnd::WaveletBoundary, L::Int=maxtransformlevels(x)) =
-    idwt(x, dual(Filterbank(w)), bnd, L)
+idwt(x, w::DiscreteWavelet, bnd::WaveletBoundary, L::Int=maxtransformlevels(x)) =
+    idwt(x, Primal, w, bnd, L)
 
-dwt(x, w::DiscreteWavelet, bnd::WaveletBoundary, L::Int=maxtransformlevels(x)) =
+dwt(x, s::Prl, w::DiscreteWavelet, bnd::WaveletBoundary, L::Int=maxtransformlevels(x)) =
     dwt(x, Filterbank(w), bnd, L)
 
-idwt(x, w::DiscreteWavelet, bnd::WaveletBoundary, L::Int=maxtransformlevels(x)) =
+idwt(x, s::Prl, w::DiscreteWavelet, bnd::WaveletBoundary, L::Int=maxtransformlevels(x)) =
     idwt(x, Filterbank(w), bnd, L)
 
+dwt(x, s::Dul, w::DiscreteWavelet, bnd::WaveletBoundary, L::Int=maxtransformlevels(x)) =
+    dwt(x, DualFilterbank(w), bnd, L)
+
+idwt(x, s::Dul, w::DiscreteWavelet, bnd::WaveletBoundary, L::Int=maxtransformlevels(x)) =
+    idwt(x, DualFilterbank(w), bnd, L)
+
 function dwt!(y, x, fb::Filterbank, bnd::WaveletBoundary, L::Int=maxtransformlevels(x))
     lx = length(x)
     @assert 2^L <= lx
@@ -64,12 +69,18 @@ function idwt(x, fb::Filterbank, bnd::WaveletBoundary, L::Int=maxtransformlevels
     y
 end
 
-function full_dwt{T}(x, w::DiscreteWavelet{T}, bnd::WaveletBoundary)
-    coefs = scaling_coefficients(x, w, bnd)
-    dwt(coefs, Filterbank(w), bnd)
+full_dwt{T}(x, w::DiscreteWavelet{T}, bnd::WaveletBoundary) =
+    full_dwt(x, Primal, w, bnd)
+
+full_idwt{T}(x, w::DiscreteWavelet{T}, bnd::WaveletBoundary) =
+    full_idwt(x, Primal, w, bnd)
+
+function full_dwt{T}(x, s::Side, w::DiscreteWavelet{T}, bnd::WaveletBoundary)
+    coefs = scaling_coefficients(x, s, w, bnd)
+    dwt(coefs, s, w, bnd)
 end
 
-function full_idwt{T}(x, w::DiscreteWavelet{T}, bnd::WaveletBoundary)
-    scalingcoefs = idwt(x, Filterbank(w), bnd)
-    scaling_coefficients_to_dyadic_grid(scalingcoefs, w, bnd)
+function full_idwt{T}(x, s::Side, w::DiscreteWavelet{T}, bnd::WaveletBoundary)
+    scalingcoefs = idwt(x, s, w, bnd)
+    scaling_coefficients_to_dyadic_grid(scalingcoefs, s, w, bnd)
 end

src/dwt/scaling_coefficients.jl

@@ -39,9 +39,12 @@
 """
 function scaling_coefficients end
 
+scaling_coefficients(f::Function, w::DWT.DiscreteWavelet, L::Int, fembedding, a::Real=0, b::Real=1; options...) =
+    scaling_coefficients(f, Dual, w, L, fembedding, a, b; options...)
+
 
 # Function on the interval [a,b] to function on [0,1]
-function scaling_coefficients(f::Function, w::DWT.DiscreteWavelet, L::Int, fembedding, a::Real=0, b::Real=1; side::Side=Dul(), options...)
+function scaling_coefficients(f::Function, side::DWT.Side, w::DWT.DiscreteWavelet, L::Int, fembedding, a::Real=0, b::Real=1; options...)
       T = promote_type(eltype(w), eltype(a), eltype(b))
       a = T(a); b = T(b)
       flt = filter(side, Scl(), w)
@@ -59,13 +62,8 @@ function scaling_coefficients{T}(f::Function, s::CompactSequence{T}, L::Int, fem
 end
 
 # Convenience function: wavelet to filter
-scaling_coefficients{T}(f::AbstractArray, w::DiscreteWavelet{T}, bnd::PeriodicBoundary; options...) =
-    scaling_coefficients(f, _scalingcoefficient_filter(filter(Dul(), Scl(), w)), PeriodicEmbedding(); options...)
-
-# Convenience function: wavelet to filter
-scaling_coefficients!{T}(c::AbstractArray, f::AbstractArray, w::DiscreteWavelet{T}, bnd::PeriodicBoundary; options...) =
-    scaling_coefficients!(c, f, _scalingcoefficient_filter(filter(Dul(), Scl(), w)), PeriodicEmbedding(); options...)
-
+scaling_coefficients{T}(f::AbstractArray, s::Side, w::DiscreteWavelet{T}, bnd::PeriodicBoundary; options...) =
+    scaling_coefficients(f, _scalingcoefficient_filter(filter(inv(s), Scl(), w)), PeriodicEmbedding(); options...)
 
 # function samples to scaling coeffients
 function scaling_coefficients{T}(f::AbstractArray, filter::CompactSequence{T}, fembedding; n::Int=length(f), options...)
@@ -91,11 +89,16 @@ end
 _scalingcoefficient_filter(f::CompactSequence) =
     reverse(CompactSequence(recursion_algorithm(f, 0), f.offset))
 
+# Remove
+# "Transforms scaling coeffients back to function evaluations on the dyadic grid."
+# scaling_coefficients_to_dyadic_grid{T}(scaling_coefficients::AbstractArray, w::DWT.DiscreteWavelet{T}, bnd::WaveletBoundary, d=ndyadicscales(scaling_coefficients); options...) =
+    # scaling_coefficients_to_dyadic_grid(scaling_coefficients, Primal, w, bmd, d; options...)
+
 "Transforms scaling coeffients back to function evaluations on the dyadic grid."
-function scaling_coefficients_to_dyadic_grid{T}(scaling_coefficients::AbstractArray, w::DWT.DiscreteWavelet{T}, bnd::WaveletBoundary, d=ndyadicscales(scaling_coefficients); grid=false, options...)
-    function_evals = zeros(T,DWT.scaling_coefficients_to_dyadic_grid_length(w,d))
-    scratch = zeros(DWT.scaling_coefficients_to_dyadic_grid_scratch_length(w,d))
-    scratch2 = zeros(DWT.scaling_coefficients_to_dyadic_grid_scratch2_length(w,d))
+function scaling_coefficients_to_dyadic_grid{T}(scaling_coefficients::AbstractArray, s::Side, w::DWT.DiscreteWavelet{T}, bnd::WaveletBoundary, d=ndyadicscales(scaling_coefficients); grid=false, options...)
+    function_evals = zeros(T,DWT.scaling_coefficients_to_dyadic_grid_length(s,w,d))
+    scratch = zeros(DWT.scaling_coefficients_to_dyadic_grid_scratch_length(s,w,d))
+    scratch2 = zeros(DWT.scaling_coefficients_to_dyadic_grid_scratch2_length(s,w,d))
 
     scaling_coefficients_to_dyadic_grid!(function_evals, scaling_coefficients, w, bnd, scratch, scratch2; options...)
     grid ?
@@ -105,16 +108,21 @@ end
 
 # Scaling coefficients to function evaluations on dyadic grid (assumes periodic extension)
 "In place method of scaling_coefficients_to_dyadic_grid."
-function scaling_coefficients_to_dyadic_grid!{T}(function_evals::AbstractArray, scaling_coefficients::AbstractArray, w::DWT.DiscreteWavelet{T}, ::DWT.PeriodicBoundary, scratch=nothing, scratch2=nothing; grid=false, options...)
+scaling_coefficients_to_dyadic_grid!{T}(function_evals::AbstractArray, scaling_coefficients::AbstractArray, w::DWT.DiscreteWavelet{T}, bnd::DWT.PeriodicBoundary, scratch=nothing, scratch2=nothing; options...) =
+    scaling_coefficients_to_dyadic_grid!(function_evals, scaling_coefficients, Primal, w, bnd, scratch, scratch2; options...)
+
+
+"In place method of scaling_coefficients_to_dyadic_grid."
+function scaling_coefficients_to_dyadic_grid!{T}(function_evals::AbstractArray, scaling_coefficients::AbstractArray, s::Side, w::DWT.DiscreteWavelet{T}, ::DWT.PeriodicBoundary, scratch=nothing, scratch2=nothing; grid=false, options...)
     d = Int(log2(length(function_evals)))
     @assert length(function_evals) == length(scratch)
-    @assert DWT.scaling_coefficients_to_dyadic_grid_scratch2_length(w, d) == length(scratch2)
+    @assert DWT.scaling_coefficients_to_dyadic_grid_scratch2_length(s, w, d) == length(scratch2)
     @assert isdyadic(scaling_coefficients)
     j = ndyadicscales(scaling_coefficients)
     function_evals[:] = T(0)
     for (c_i,c) in enumerate(scaling_coefficients)
         k = c_i - 1
-        DWT.evaluate_periodic_in_dyadic_points!(scratch, DWT.Prl(), DWT.Scl(), w, j, k, d, scratch2, nothing)
+        DWT.evaluate_periodic_in_dyadic_points!(scratch, s, DWT.Scl(), w, j, k, d, scratch2, nothing)
         scale!(scratch, c)
         for i in 1:length(function_evals)
             function_evals[i] += scratch[i]
@@ -124,12 +132,12 @@ function scaling_coefficients_to_dyadic_grid!{T}(function_evals::AbstractArray,
     nothing
 end
 
-scaling_coefficients_to_dyadic_grid_length(w::DWT.DiscreteWavelet, d::Int) = 1<<d
+scaling_coefficients_to_dyadic_grid_length(s::Side, w::DWT.DiscreteWavelet, d::Int) = 1<<d
 
-scaling_coefficients_to_dyadic_grid_scratch_length(w::DWT.DiscreteWavelet, d::Int) =
-    scaling_coefficients_to_dyadic_grid_length(w::DWT.DiscreteWavelet, d::Int)
-scaling_coefficients_to_dyadic_grid_scratch2_length(w::DWT.DiscreteWavelet, d::Int) =
-    DWT.evaluate_periodic_in_dyadic_points_scratch_length(DWT.Prl(), DWT.Scl(), w, d, 0, d)
+scaling_coefficients_to_dyadic_grid_scratch_length(s::Side, w::DWT.DiscreteWavelet, d::Int) =
+    scaling_coefficients_to_dyadic_grid_length(s, w, d)
+scaling_coefficients_to_dyadic_grid_scratch2_length(s::Side, w::DWT.DiscreteWavelet, d::Int) =
+    DWT.evaluate_periodic_in_dyadic_points_scratch_length(s, DWT.Scl(), w, d, 0, d)
 
 function DWT.support(side::Side, n::Int, i::Int, l::Int, w::DiscreteWavelet)
     kind, j, k = wavelet_index(n,i,l)

src/dwt/util/periodize.jl

@@ -33,18 +33,14 @@ end
 #          src1
 #        + ________________
 #  =>      res1, res2, res3
-function _periodize!{T}(dest::AbstractArray{T}, src::AbstractArray{T}, istart)
-  L = length(dest)
-  j = mod(istart, L)
-  srclength = length(src)
-  (j == 0) && (j = L)
-  for i in 1:length(dest)
-    t = T(0)
-    for m in i:L:srclength
-      t += src[m]
+function _periodize!{T}(dest::AbstractArray{T}, src::AbstractArray{T}, istart::Int, step::Int=1)
+    L = length(dest)
+    srclength = length(src)
+    for i in 1:L
+        t = T(0)
+        for m in mod(istart-1+step*(i-1),step*L)+1:step*L:srclength
+            t += src[m]
+        end
+        dest[i] = t
     end
-    dest[j] = t
-    j += 1
-    (j > L) && (j = 1)
-  end
 end

src/dwt/util/recipes.jl

@@ -40,5 +40,5 @@ end
     Dual, w
   end
 end
-
-@recipe f(x::LinSpace, f::AbstractVector) = collect(x), f
+# REMOVE
+# @recipe f(x::LinSpace, f::AbstractVector) = collect(x), f

src/dwt/util/recursion.jl

@@ -88,7 +88,8 @@ function dyadicpointsofrecursion{T}(side::Side, kind::Kind, w::DiscreteWavelet{T
     if L >= 0
         linspace(T(s[1]), T(s[2]), (1<<L)*H+1)
     else
-        linspace(T(s[1]), T(s[2]), H+1)[1:1<<-L:end]
+        # include zero, therefore, do some rounding
+        (1<<-L)*(cld(s[1],1<<-L):fld(s[2],1<<-L))
     end
 end
 

src/dwt/util/util_functions.jl

@@ -1,5 +1,4 @@
-# Convenience method
-eltype(x, y) = promote_type(eltype(x), eltype(y))
+# Convenience methods
 # WAVELET INDEXING AND SIZES
 
 maxtransformlevels(x::AbstractArray) = maxtransformlevels(minimum(size(x)))