From 4789ae85443613274e66d5cb0a3ad27c62d5fdd7 Mon Sep 17 00:00:00 2001 From: vincentcp Date: Mon, 9 Apr 2018 14:48:50 +0200 Subject: [PATCH 1/4] add dual dwt functionality --- notebooks/README.ipynb | 2 +- src/WaveletsCopy.jl | 2 +- src/dwt/discretewavelets.jl | 3 ++ src/dwt/dwttransform.jl | 36 ++++++++++++++-------- src/dwt/scaling_coefficients.jl | 54 ++++++++++++++++++++++----------- test/runtests.jl | 2 +- test/suite_dwtstep.jl | 4 +++ test/suite_evaluation.jl | 2 +- 8 files changed, 72 insertions(+), 33 deletions(-) diff --git a/notebooks/README.ipynb b/notebooks/README.ipynb index a3ab2a8..16d85bc 100644 --- a/notebooks/README.ipynb +++ b/notebooks/README.ipynb @@ -30,7 +30,7 @@ }, "outputs": [], "source": [ - "using WaveletsCopys \n", + "using WaveletsCopy \n", "using Plots" ] }, diff --git a/src/WaveletsCopy.jl b/src/WaveletsCopy.jl index 855fa65..8d61a02 100644 --- a/src/WaveletsCopy.jl +++ b/src/WaveletsCopy.jl @@ -1,5 +1,5 @@ # disable precompilation during development -#__precompile__() +__precompile__(true) module WaveletsCopy diff --git a/src/dwt/discretewavelets.jl b/src/dwt/discretewavelets.jl index 666887c..4b93892 100644 --- a/src/dwt/discretewavelets.jl +++ b/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 diff --git a/src/dwt/dwttransform.jl b/src/dwt/dwttransform.jl index 302c6fc..f432ac2 100644 --- a/src/dwt/dwttransform.jl +++ b/src/dwt/dwttransform.jl @@ -12,18 +12,24 @@ 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 +70,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 diff --git a/src/dwt/scaling_coefficients.jl b/src/dwt/scaling_coefficients.jl index b6bbecb..eae58c7 100644 --- a/src/dwt/scaling_coefficients.jl +++ b/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) @@ -58,13 +61,21 @@ function scaling_coefficients{T}(f::Function, s::CompactSequence{T}, L::Int, fem scaling_coefficients(fcoefs, filter, fembedding; options...) end -# Convenience function: wavelet to filter +# Convenience function: default scaling_coefficients{T}(f::AbstractArray, w::DiscreteWavelet{T}, bnd::PeriodicBoundary; options...) = - scaling_coefficients(f, _scalingcoefficient_filter(filter(Dul(), Scl(), w)), PeriodicEmbedding(); options...) + scaling_coefficients(f, Primal, w, bnd; options...) -# Convenience function: wavelet to filter +# Convenience function: default 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!(c, f, Primal, w, bnd; options...) + +# Convenience function: wavelet to filter +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...) + +# Convenience function: wavelet to filter +scaling_coefficients!{T}(c::AbstractArray, f::AbstractArray, s::Side, w::DiscreteWavelet{T}, bnd::PeriodicBoundary; options...) = + scaling_coefficients!(c, f, _scalingcoefficient_filter(filter(inv(s), Scl(), w)), PeriodicEmbedding(); options...) # function samples to scaling coeffients @@ -92,10 +103,14 @@ _scalingcoefficient_filter(f::CompactSequence) = reverse(CompactSequence(recursion_algorithm(f, 0), f.offset)) "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)) +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, 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 +120,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 +144,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<x^p, w, d, nothing, side=Primal) + scaling_coefficients = DWT.scaling_coefficients(x->x^p, Primal, w, d, nothing) for k in 0:(1< Date: Mon, 9 Apr 2018 15:46:31 +0200 Subject: [PATCH 2/4] added tests --- notebooks/README.ipynb | 4 +- src/WaveletsCopy.jl | 2 +- src/dwt/dwtstep.jl | 79 +++++------ src/dwt/dwttransform.jl | 1 - src/dwt/scaling_coefficients.jl | 20 +-- src/dwt/util/recipes.jl | 4 +- src/dwt/util/util_functions.jl | 3 +- src/dwt/wvlt.jl | 3 +- src/filterbank.jl | 223 -------------------------------- test/suite_evaluation.jl | 39 ++++-- 10 files changed, 78 insertions(+), 300 deletions(-) delete mode 100644 src/filterbank.jl diff --git a/notebooks/README.ipynb b/notebooks/README.ipynb index 16d85bc..6340e80 100644 --- a/notebooks/README.ipynb +++ b/notebooks/README.ipynb @@ -4,8 +4,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "[![Build Status](https://travis-ci.org/vincentcp/Wavelets.jl.svg?branch=master)](https://travis-ci.org/vincentcp/Wavelets.jl)\n", - "[![Coverage Status](https://coveralls.io/repos/github/vincentcp/Wavelets.jl/badge.svg?branch=master)](https://coveralls.io/github/vincentcp/Wavelets.jl?branch=master)\n", + "[![Build Status](https://travis-ci.org/vincentcp/WaveletsCopy.jl.svg?branch=master)](https://travis-ci.org/vincentcp/WaveletsCopy.jl)\n", + "[![Coverage Status](https://coveralls.io/repos/github/vincentcp/WaveletsCopy.jl/badge.svg?branch=master)](https://coveralls.io/github/vincentcp/WaveletsCopy.jl?branch=master)\n", "# Wavelets\n", "A julia package for a fast Discrete Wavelet Transform, and plotting and evaluation of wavelets " ] diff --git a/src/WaveletsCopy.jl b/src/WaveletsCopy.jl index 8d61a02..da66893 100644 --- a/src/WaveletsCopy.jl +++ b/src/WaveletsCopy.jl @@ -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 diff --git a/src/dwt/dwtstep.jl b/src/dwt/dwtstep.jl index f2bb4d7..2ac78b9 100644 --- a/src/dwt/dwtstep.jl +++ b/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 diff --git a/src/dwt/dwttransform.jl b/src/dwt/dwttransform.jl index f432ac2..4e6a148 100644 --- a/src/dwt/dwttransform.jl +++ b/src/dwt/dwttransform.jl @@ -9,7 +9,6 @@ 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(x, w::DiscreteWavelet, bnd::WaveletBoundary, L::Int=maxtransformlevels(x)) = diff --git a/src/dwt/scaling_coefficients.jl b/src/dwt/scaling_coefficients.jl index eae58c7..0b4653f 100644 --- a/src/dwt/scaling_coefficients.jl +++ b/src/dwt/scaling_coefficients.jl @@ -61,23 +61,10 @@ function scaling_coefficients{T}(f::Function, s::CompactSequence{T}, L::Int, fem scaling_coefficients(fcoefs, filter, fembedding; options...) end -# Convenience function: default -scaling_coefficients{T}(f::AbstractArray, w::DiscreteWavelet{T}, bnd::PeriodicBoundary; options...) = - scaling_coefficients(f, Primal, w, bnd; options...) - -# Convenience function: default -scaling_coefficients!{T}(c::AbstractArray, f::AbstractArray, w::DiscreteWavelet{T}, bnd::PeriodicBoundary; options...) = - scaling_coefficients!(c, f, Primal, w, bnd; options...) - # Convenience function: wavelet to filter 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...) -# Convenience function: wavelet to filter -scaling_coefficients!{T}(c::AbstractArray, f::AbstractArray, s::Side, w::DiscreteWavelet{T}, bnd::PeriodicBoundary; options...) = - scaling_coefficients!(c, 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...) @assert isdyadic(f) @@ -102,9 +89,10 @@ end _scalingcoefficient_filter(f::CompactSequence) = reverse(CompactSequence(recursion_algorithm(f, 0), f.offset)) -"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...) +# 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, s::Side, w::DWT.DiscreteWavelet{T}, bnd::WaveletBoundary, d=ndyadicscales(scaling_coefficients); grid=false, options...) diff --git a/src/dwt/util/recipes.jl b/src/dwt/util/recipes.jl index 7de1b4d..d308884 100644 --- a/src/dwt/util/recipes.jl +++ b/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 diff --git a/src/dwt/util/util_functions.jl b/src/dwt/util/util_functions.jl index 661c3f7..e716aab 100644 --- a/src/dwt/util/util_functions.jl +++ b/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))) diff --git a/src/dwt/wvlt.jl b/src/dwt/wvlt.jl index 347940e..f690c35 100644 --- a/src/dwt/wvlt.jl +++ b/src/dwt/wvlt.jl @@ -30,10 +30,11 @@ scaling = Scl(); DWT.name(::Scl) = "scaling" wavelet = Wvl(); DWT.name(::Wvl) = "wavelet" Base.inv(::Prl) = Dul() Base.inv(::Dul) = Prl() +Base.inv(::Type{Prl}) = Dul +Base.inv(::Type{Dul}) = Prl ############################################################################### # vanishingmoments ############################################################################### -vanishingmoments{WT<:DiscreteWavelet}(side::Side, kind::Kind, ::Type{WT}) = vanishingmoments(s, WT) vanishingmoments{WT<:DiscreteWavelet}(::Prl, ::Type{WT}) = throw("unimplemented") vanishingmoments{WT<:DiscreteWavelet}(::Dul, W::Type{WT}) = _vanishingmoments(Prl(), W, is_orthogonal(W)) _vanishingmoments(::Prl, W, is_orthogonal::Type{True}) = vanishingmoments(Prl(), W) diff --git a/src/filterbank.jl b/src/filterbank.jl deleted file mode 100644 index b3a6257..0000000 --- a/src/filterbank.jl +++ /dev/null @@ -1,223 +0,0 @@ -# filterbank.jl - -module Filterbanks - -using WaveletsCopy.Sequences - -import Base: transpose, ctranspose, eltype, getindex - -export FilterPair, FilterMatrix, PolyphaseMatrix, Filterbank - -export polyphasematrix, modulationmatrix, polyphase_analysis!, polyphase_synthesis!, - lowpassfilter, highpassfilter, primal_lowpassfilter, primal_highpassfilter, - dual_lowpassfilter, dual_highpassfilter - - -"A pair of two filters." -struct FilterPair{F1 <: Sequence,F2 <: Sequence} - f1 :: F1 - f2 :: F2 -end - -eltype{F1,F2}(::Type{FilterPair{F1,F2}}) = eltype(F1) - -lowpassfilter(fp::FilterPair) = fp.f1 - -highpassfilter(fp::FilterPair) = fp.f2 - - -function getindex(fp::FilterPair, k::Int) - if k == 1 - fp.f1 - elseif k == 2 - fp.f2 - else - throw(BoundsError()) - end -end - -ctranspose(fb::FilterPair) = FilterPair(fb.f1', fb.f2') - - -"A 2x2 matrix of filters." -struct FilterMatrix{F11,F12,F21,F22} - a11 :: F11 - a12 :: F12 - a21 :: F21 - a22 :: F22 -end - -PolyphaseMatrix = FilterMatrix - -eltype{F11,F12,F21,F22}(::Type{FilterMatrix{F11,F12,F21,F22}}) = eltype(F11) - -transpose(m::FilterMatrix) = FilterMatrix(m.a11, m.a21, m.a12, m.a22) - -ctranspose(m::FilterMatrix) = FilterMatrix(m.a11', m.a21', m.a12', m.a22') - -# Done removed: call(m::FilterMatrix, z) = [ztransform(m.a11, z) ztransform(m.a12, z); ztransform(m.a21, z) ztransform(m.a22, z)] -(m::FilterMatrix)(z) = [ztransform(m.a11, z) ztransform(m.a12, z); ztransform(m.a21, z) ztransform(m.a22, z)] - -polyphasematrix(fp::FilterPair) = FilterMatrix(evenpart(fp[1]), evenpart(fp[2]), oddpart(fp[1]), oddpart(fp[2])) - -modulationmatrix(fp::FilterPair) = FilterMatrix(fp[1], alternating(fp[1]), fp[2], alternating(fp[2])) - - - -"A Filterbank groups several objects related to a two-phase filterbank." -struct Filterbank{P1 <: FilterMatrix, P2 <: FilterMatrix, FP1 <: FilterPair, FP2 <: FilterPair} - "The primal filter pair, to be used on the synthesis side." - primal_pair :: FP1 - "The dual filter pair, to be used on the analysis side." - dual_pair :: FP2 - "The polyphase matrix for the synthesis side." - pm_synthesis :: P2 - "The polyphase matrix for the analysis side." - pm_analysis :: P1 -end - -Filterbank(lowpass::Sequence) = Filterbank(FilterPair(lowpass, alternating_flip(lowpass))) - -Filterbank(primal_lowpass::Sequence, dual_lowpass::Sequence) = - Filterbank(FilterPair(primal_lowpass, alternating_flip(dual_lowpass)), - FilterPair(dual_lowpass, alternating_flip(primal_lowpass))) - -# If no dual pair is given, assume orthogonality. -# TODO: verify orthogonality -Filterbank(primal_pair::FilterPair) = Filterbank(primal_pair, primal_pair) - -Filterbank(primal_pair::FilterPair, dual_pair::FilterPair) = - Filterbank(primal_pair, dual_pair, polyphasematrix(primal_pair), polyphasematrix(dual_pair)') - -eltype{P1,P2,FP1,FP2}(::Type{Filterbank{P1,P2,FP1,FP2}}) = eltype(P1) - -primal_lowpassfilter(fb::Filterbank) = lowpassfilter(fb.primal_pair) -primal_highpassfilter(fb::Filterbank) = highpassfilter(fb.primal_pair) - -dual_lowpassfilter(fb::Filterbank) = lowpassfilter(fb.dual_pair) -dual_highpassfilter(fb::Filterbank) = highpassfilter(fb.dual_pair) - -# "Return a temporary array with elements of the given type and of length L. Memory is allocated only -# on the first call of this function and reused afterwards." -# @generated temparray{T,L}(::Type{T}, ::Type{Val{L}}) = quote -# a = $(Array(T,L-1)) -# return a -# end - -# Split the (finite) signal x embedded in (infinite) embedding into a low pass signal y1 -# and a high pass signal y2. The filters are given by the polyphasematrix f. -function polyphase_analysis!(y1, y2, x, f::PolyphaseMatrix, embedding) - Heven = f.a11 - Hodd = f.a12 - Geven = f.a21 - Godd = f.a22 - - T = eltype(y1) - # Determine for which range the convolution of the filters with x results - # in evaluations of x that not occur outside of the embedding. - # Remark (TODO?): this could be even more efficient. - L = length(x) - lower_i = min(length(y1)-1, maximum(map(lastindex,(Heven,Geven,Godd,Hodd)))) - L%2 == 0? - upper_i = L>>1 - 1 + minimum(map(firstindex,(Heven,Geven,Godd,Hodd))) : - upper_i = L>>1 - L%2 + minimum(map(firstindex,(Heven,Geven,Godd,Hodd))) - upper_i = max(0, upper_i) - # Lower boundary region - for i in 0:lower_i - y1i = zero(T) - y2i = zero(T) - for l = firstindex(Heven):lastindex(Heven) - y1i += Heven[l] * embedding[x, 2*i-2*l] - end - for l = firstindex(Geven):lastindex(Geven) - y2i += Geven[l] * embedding[x, 2*i-2*l] - end - for l = firstindex(Hodd):lastindex(Hodd) - y1i += Hodd[l] * embedding[x, 2*i-2*l+1] - end - for l = firstindex(Godd):lastindex(Godd) - y2i += Godd[l] * embedding[x, 2*i-2*l+1] - end - y1[i+1] = y1i - y2[i+1] = y2i - end - - # Middle region - @inbounds for i in lower_i+1:upper_i-1 - # for l in minimum(map(firstindex,(Heven,Geven,Hodd,Godd))):maximum(map(lastindex,(Heven,Geven,Hodd,Godd))) - # @assert (0<=2(i-l)>1 - 1 - xj_e = zero(T) - xj_o = zero(T) - for l = firstindex(Heven):lastindex(Heven) - xj_e += Heven[l] * embedding[y1,j-l] - end - for l = firstindex(Geven):lastindex(Geven) - xj_e += Geven[l] * embedding[y2,j-l] - end - for l = firstindex(Hodd):lastindex(Hodd) - xj_o += Hodd[l] * embedding[y1,j-l] - end - for l = firstindex(Godd):lastindex(Godd) - xj_o += Godd[l] * embedding[y2,j-l] - end - x[2*j+1] = xj_e - x[2*j+2] = xj_o - end -end - - - -end # module diff --git a/test/suite_evaluation.jl b/test/suite_evaluation.jl index f799469..d8c3236 100644 --- a/test/suite_evaluation.jl +++ b/test/suite_evaluation.jl @@ -209,7 +209,7 @@ function coefficient_util_test() @testset "$(rpad("coefficient util tests",P))" begin levels = [0,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3] for i in 1:16 - @test DWT.level(16,i) == levels[i] + @test DWT.level(16,i) == levels[i] end @test DWT.wavelet_index(4,1,0) == (scaling, 2, 0) @test DWT.wavelet_index(4,2,0) == (scaling, 2, 1) @@ -223,6 +223,10 @@ function coefficient_util_test() @test DWT.wavelet_index(4,2,2) == (wavelet, 0, 0) @test DWT.wavelet_index(4,3,2) == (wavelet, 1, 0) @test DWT.wavelet_index(4,4,2) == (wavelet, 1, 1) + @test DWT.wavelet_indices(3)[1] == (scaling, 0, 0) + @test DWT.wavelet_indices(3)[2] == (wavelet, 0, 0) + @test DWT.wavelet_indices(3)[3] == (wavelet, 1, 0) + @test DWT.wavelet_indices(3)[4] == (wavelet, 1, 1) @test DWT.support(Primal, 8, 1, 3, DWT.db1) == (0,1) @test DWT.support(Primal, 8, 2, 3, DWT.db1) == (0,1) @test DWT.support(Primal, 8, 3, 3, DWT.db1) == (.0,.5) @@ -240,6 +244,8 @@ function implementation_test() @test DWT.name(DWT.wavelet) == "wavelet" @test DWT.name(DWT.Primal) == "primal" @test DWT.name(DWT.Dual) == "dual" + print_implemented_wavelets() + print_all_implemented_wavelets() @test DWT.is_symmetric(DWT.TestWavelet{Float16}) == False @test DWT.is_biorthogonal(DWT.TestWavelet{Float16}) == False @test DWT.is_orthogonal(DWT.TestWavelet{Float16}) == False @@ -340,20 +346,26 @@ function implementation_test() end function eval_wavelet_test() - @testset "$(rpad("evaluate wavelet in point",P))" begin - @test DWT._periodize((.25,.5),-1,1)==((.25,.5),) - @test DWT._periodize((-1.25,1.5),-1,1)==((-1,1),) - @test DWT._periodize((-1.25,.5),-1,1)==((.75,1.),(-1.,.5)) - @test DWT._periodize((-1.5,0.),-1,1)==((0.5,1.),(-1,0)) - for side in (Primal, Dual) - for kind in (scaling, DWT.wavelet) - @test DWT.in_periodic_support(1,DWT.periodic_support(side,kind, DWT.cdf11, 0,0)...) - @test DWT.in_periodic_support(0,DWT.periodic_support(side,kind, DWT.cdf11, 3,0)...) - @test !DWT.in_periodic_support(0,DWT.periodic_support(side,kind, DWT.cdf11, 1,1)...) - end + @testset "$(rpad("evaluate wavelet in point",P))" begin + @test DWT._periodize((.25,.5),-1,1)==((.25,.5),) + @test DWT._periodize((-1.25,1.5),-1,1)==((-1,1),) + @test DWT._periodize((-1.25,.5),-1,1)==((.75,1.),(-1.,.5)) + @test DWT._periodize((-1.5,0.),-1,1)==((0.5,1.),(-1,0)) + + t = linspace(0,1,10) + for side in (Primal, Dual) + for kind in (scaling, DWT.wavelet) + @test DWT.in_periodic_support(1,DWT.periodic_support(side,kind, DWT.cdf11, 0,0)...) + @test DWT.in_periodic_support(0,DWT.periodic_support(side,kind, DWT.cdf11, 3,0)...) + @test !DWT.in_periodic_support(0,DWT.periodic_support(side,kind, DWT.cdf11, 1,1)...) + for w in [db2, cdf15] + @test 1+evaluate_periodic.(Primal, scaling, w, 3, 0, t-1) ≈ 1+evaluate.(Primal, scaling, w, 3, 0, t)+evaluate.(Primal, scaling, w, 3, 0, t-1)+evaluate.(Primal, scaling, w, 3, 0, t+1) + end + end + end end - end end + implementation_test() recursiontest() primalfunctiontest() @@ -369,6 +381,7 @@ coefficient_util_test() using Plots plot(DWT.cdf11) plot(DWT.db1) +plot(DWT.db2,periodic=true) # using Plots # plot(layout=(2,2)) From 5cd403f79d15a190652790584d6e300d42d60328 Mon Sep 17 00:00:00 2001 From: vincentcp Date: Tue, 10 Apr 2018 10:07:09 +0200 Subject: [PATCH 3/4] bugfixes in pointswise evaluation --- README.md | 10 ++++------ notebooks/README.ipynb | 2 +- src/dwt/util/periodize.jl | 22 +++++++++------------- src/dwt/util/recursion.jl | 3 ++- src/dwt/wvlt.jl | 36 +++++++++++++++++++++++++++--------- test/suite_dwtstep.jl | 6 ++---- test/suite_evaluation.jl | 13 +++++++++---- 7 files changed, 54 insertions(+), 38 deletions(-) diff --git a/README.md b/README.md index 2611cd7..42327d3 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ [![Build Status](https://travis-ci.org/vincentcp/Wavelets.jl.svg?branch=master)](https://travis-ci.org/vincentcp/Wavelets.jl) [![Coverage Status](https://coveralls.io/repos/github/vincentcp/Wavelets.jl/badge.svg?branch=master)](https://coveralls.io/github/vincentcp/Wavelets.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) - - diff --git a/notebooks/README.ipynb b/notebooks/README.ipynb index 6340e80..902bf78 100644 --- a/notebooks/README.ipynb +++ b/notebooks/README.ipynb @@ -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)`" diff --git a/src/dwt/util/periodize.jl b/src/dwt/util/periodize.jl index 4d73eee..b994709 100644 --- a/src/dwt/util/periodize.jl +++ b/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 diff --git a/src/dwt/util/recursion.jl b/src/dwt/util/recursion.jl index 147021e..c899d44 100644 --- a/src/dwt/util/recursion.jl +++ b/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< Date: Tue, 10 Apr 2018 10:44:21 +0200 Subject: [PATCH 4/4] more bugfixes and testing on pointwise evaluation --- src/dwt/wvlt.jl | 22 ++++++++++++---------- test/suite_evaluation.jl | 12 ++++-------- 2 files changed, 16 insertions(+), 18 deletions(-) diff --git a/src/dwt/wvlt.jl b/src/dwt/wvlt.jl index c477500..d9ac81f 100644 --- a/src/dwt/wvlt.jl +++ b/src/dwt/wvlt.jl @@ -229,11 +229,12 @@ function evaluate_periodic_in_dyadic_points!{T}(f::AbstractArray{T}, side::DWT.S # Wavelet evaluation uses more schratch than scaling evaluation DWT.evaluate_in_dyadic_points!(scratch, side, kind, w, j ,k ,d, scratch2, scratch3; options...) - try - DWT._periodize!(f, scratch, -Int((1<