From 9a4337f9c00632c278cb79ce78b27976aca6f037 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thimot=C3=A9=20Dupuch?= Date: Thu, 14 Aug 2025 20:41:47 +0200 Subject: [PATCH 1/3] Only code reorganization (src/ directory) --- src/Plot/Plot.jl | 8 + src/{mod/Plot.jl => Plot/plot.jl} | 55 +- src/Threshold/Threshold.jl | 27 + src/Threshold/basis_functions.jl | 56 ++ src/Threshold/denoising.jl | 122 +++++ src/Threshold/entropy.jl | 129 +++++ src/Threshold/threshold.jl | 129 +++++ src/Transforms/Transforms.jl | 10 + src/Transforms/transforms.jl | 226 ++++++++ src/{mod => Transforms}/transforms_filter.jl | 258 ++++----- src/{mod => Transforms}/transforms_lifting.jl | 208 ++++---- .../transforms_maximal_overlap.jl | 12 +- src/Util/Util.jl | 44 ++ src/Util/dyadic.jl | 20 + src/Util/non_dyadic.jl | 27 + src/{mod/Util.jl => Util/util.jl} | 286 ++++------ src/WT/WT.jl | 13 + src/WT/wt.jl | 454 ++++++++++++++++ src/Wavelets.jl | 71 ++- src/mod/Threshold.jl | 462 ---------------- src/mod/Transforms.jl | 239 --------- src/mod/WT.jl | 493 ------------------ 22 files changed, 1691 insertions(+), 1658 deletions(-) create mode 100644 src/Plot/Plot.jl rename src/{mod/Plot.jl => Plot/plot.jl} (70%) create mode 100644 src/Threshold/Threshold.jl create mode 100644 src/Threshold/basis_functions.jl create mode 100644 src/Threshold/denoising.jl create mode 100644 src/Threshold/entropy.jl create mode 100644 src/Threshold/threshold.jl create mode 100644 src/Transforms/Transforms.jl create mode 100644 src/Transforms/transforms.jl rename src/{mod => Transforms}/transforms_filter.jl (65%) rename src/{mod => Transforms}/transforms_lifting.jl (69%) rename src/{mod => Transforms}/transforms_maximal_overlap.jl (90%) create mode 100644 src/Util/Util.jl create mode 100644 src/Util/dyadic.jl create mode 100644 src/Util/non_dyadic.jl rename src/{mod/Util.jl => Util/util.jl} (53%) create mode 100644 src/WT/WT.jl create mode 100644 src/WT/wt.jl delete mode 100644 src/mod/Threshold.jl delete mode 100644 src/mod/Transforms.jl delete mode 100644 src/mod/WT.jl diff --git a/src/Plot/Plot.jl b/src/Plot/Plot.jl new file mode 100644 index 0000000..cac8d37 --- /dev/null +++ b/src/Plot/Plot.jl @@ -0,0 +1,8 @@ +module Plot + +include("plot.jl") + +export + wplotdots, wplotim + +end diff --git a/src/mod/Plot.jl b/src/Plot/plot.jl similarity index 70% rename from src/mod/Plot.jl rename to src/Plot/plot.jl index 21c30a4..3046a74 100644 --- a/src/mod/Plot.jl +++ b/src/Plot/plot.jl @@ -1,7 +1,6 @@ -module Plot -export wplotdots, wplotim -using ..Util, ..WT, ..Transforms -using LinearAlgebra +using ..Util +using ..WT +using ..Transforms # PLOTTING UTILITIES @@ -9,8 +8,8 @@ using LinearAlgebra # as tuple (d,l) function wplotdots(x::AbstractVector, t::Real=0, r::Real=1) isdyadic(x) || throw(ArgumentError("array must be of dyadic size")) - n = length(x) - c = wcount(x, t, level=0) + n = length(x) + c = wcount(x, t, level=0) d = Vector{Float64}(undef, c) l = Vector{Int}(undef, c) range = 0:1/n:1-eps() @@ -22,7 +21,7 @@ function wplotdots(x::AbstractVector, t::Real=0, r::Real=1) for j = 0:J-1 rind = 2^(J-1-j):2^(J-j):n for i in 1:dyadicdetailn(j) - if abs(x[dyadicdetailindex(j,i)]) >= t + if abs(x[dyadicdetailindex(j, i)]) >= t d[k] = range[rind[i]] l[k] = j k += 1 @@ -30,7 +29,7 @@ function wplotdots(x::AbstractVector, t::Real=0, r::Real=1) end end end - return (d,l) + return (d, l) end # return a matrix of detail coefficient values where row j+1 is level j @@ -42,11 +41,11 @@ function wplotim(x::AbstractVector) @inbounds begin for j = 0:J-1 - dr = dyadicdetailrange(j) - m = 2^(J-j) - for i in eachindex(dr) - A[j+1,1+(i-1)*m:i*m] .= x[dr[i]] - end + dr = dyadicdetailrange(j) + m = 2^(J - j) + for i in eachindex(dr) + A[j+1, 1+(i-1)*m:i*m] .= x[dr[i]] + end end end return A @@ -55,20 +54,20 @@ end # return an array of scaled detail coefficients and unscaled scaling coefficients # ready to be plotted as an image function wplotim(x::AbstractArray, L::Integer, wt::Union{DiscreteWavelet,Nothing}=nothing; - wabs::Bool=true, power::Real=0.7, pnorm::Real=1) + wabs::Bool=true, power::Real=0.7, pnorm::Real=1) isdyadic(x) || throw(ArgumentError("array must be of dyadic size")) dim = ndims(x) (dim == 2 || dim == 3) || throw(ArgumentError("dimension $(dim) not supported")) - n = size(x,1) - cn = size(x,3) # color dimension - n == size(x,2) || throw(ArgumentError("array must be square")) + n = size(x, 1) + cn = size(x, 3) # color dimension + n == size(x, 2) || throw(ArgumentError("array must be square")) (cn == 1 || cn == 3) || throw(ArgumentError("third dimension $(cn) not supported")) J = ndyadicscales(n) - nsc = 2^(J-L) + nsc = 2^(J - L) # do wavelet transform if wt != nothing - if size(x,3)>1 + if size(x, 3) > 1 x = dwtc(x, wt, L) else x = dwt(x, wt, L) @@ -76,21 +75,21 @@ function wplotim(x::AbstractArray, L::Integer, wt::Union{DiscreteWavelet,Nothing end # scaling coefs - scs = x[1:nsc,1:nsc,:] + scs = x[1:nsc, 1:nsc, :] scale01!(scs) # detail coefs xts = copy(x) - wabs && (broadcast!(abs,xts,xts)) - xts[1:nsc,1:nsc,:] .= 0 + wabs && (broadcast!(abs, xts, xts)) + xts[1:nsc, 1:nsc, :] .= 0 scale01!(xts) - for j=1:n, i=1:n - @inbounds xts[i,j,:] .= norm(xts[i,j,:],pnorm).^(power) + for j = 1:n, i = 1:n + @inbounds xts[i, j, :] .= norm(xts[i, j, :], pnorm) .^ (power) end # merge and reshape final image scale01!(xts) - xts[1:nsc,1:nsc,:] = scs + xts[1:nsc, 1:nsc, :] = scs return xts end # scale elements of z to the interval [0,1] @@ -98,11 +97,7 @@ function scale01!(z) mi = minimum(z) ma = maximum(z) for i in eachindex(z) - @inbounds z[i] = (z[i] - mi)/(ma-mi) + @inbounds z[i] = (z[i] - mi) / (ma - mi) end return z end - - - -end diff --git a/src/Threshold/Threshold.jl b/src/Threshold/Threshold.jl new file mode 100644 index 0000000..c96d3fe --- /dev/null +++ b/src/Threshold/Threshold.jl @@ -0,0 +1,27 @@ +module Threshold + +using LinearAlgebra: rmul!, norm +using Statistics: median! + + +include("threshold.jl") +include("basis_functions.jl") +include("denoising.jl") +include("entropy.jl") + + +export + # threshold + threshold!, threshold, HardTH, SoftTH, SemiSoftTH, SteinTH, + BiggestTH, PosTH, NegTH, + + # denoising + DNFT, VisuShrink, denoise, noisest, + + # basis functions + matchingpursuit, bestbasistree, + + # entropy + coefentropy, Entropy, ShannonEntropy, LogEnergyEntropy + +end diff --git a/src/Threshold/basis_functions.jl b/src/Threshold/basis_functions.jl new file mode 100644 index 0000000..e53cbd9 --- /dev/null +++ b/src/Threshold/basis_functions.jl @@ -0,0 +1,56 @@ + +# BASIS FUNCTIONS + +# Matching Pursuit +# see: Mallat (2009) p.642 "A wavelet tour of signal processing" +# find sparse vector y such that ||x - f(y)|| < tol approximately +# f is the operation of a M by N (M= -1 + @assert tol > 0 + r = x + n = 1 + + if !oop + y = zeros(eltype(x), length(ft(x))) + else # out of place functions f and ft + y = zeros(eltype(x), N) + tmp = similar(x, N) + ftr = similar(x, N) + aphi = similar(x, length(x)) + end + spat = zeros(eltype(x), length(y)) # sparse for atom computation + nmax == -1 && (nmax = length(y)) + + while norm(r) > tol && n <= nmax + # find largest inner product + !oop && (ftr = ft(r)) + oop && ft(ftr, r, tmp) + i = findmaxabs(ftr) + + # project on i-th atom + spat[i] = ftr[i] + !oop && (aphi = f(spat)) + oop && f(aphi, spat, tmp) + spat[i] = 0 + + # update residual, r = r - aphi + broadcast!(-, r, r, aphi) + + y[i] += ftr[i] + n += 1 + end + return y +end +function findmaxabs(x::AbstractVector) + m = abs(x[1]) + k = 1 + @inbounds for i in eachindex(x) + if abs(x[i]) > m + k = i + m = abs(x[i]) + end + end + return k +end diff --git a/src/Threshold/denoising.jl b/src/Threshold/denoising.jl new file mode 100644 index 0000000..65ac8c6 --- /dev/null +++ b/src/Threshold/denoising.jl @@ -0,0 +1,122 @@ +using ..Transforms + + + +# DENOISING + +abstract type DNFT end + +struct VisuShrink <: DNFT + th::THType # threshold type + t::Float64 # threshold for noise level sigma=1, use sigma*t in application + VisuShrink(th, t) = new(th, t) +end +# define type for signal length n +function VisuShrink(n::Int) + return VisuShrink(DEFAULT_TH, sqrt(2 * log(n))) +end + +const DEFAULT_WAVELET = wavelet(WT.sym5, WT.Filter) # default wavelet type + +# denoise signal x by thresholding in wavelet space +# estnoise is (x::AbstractArray, wt::Union{DiscreteWavelet,Nothing}) +function denoise(x::AbstractArray, + wt::Union{DiscreteWavelet,Nothing}=DEFAULT_WAVELET; + L::Int=min(maxtransformlevels(x), 6), + dnt::S=VisuShrink(size(x, 1)), + estnoise::Function=noisest, + TI::Bool=false, + nspin::Union{Int,Tuple}=tuple([8 for i = 1:ndims(x)]...)) where {S<:DNFT} + iscube(x) || throw(ArgumentError("array must be square/cube")) + sigma = estnoise(x, wt) + + if TI + wt == nothing && error("TI not supported with wt=nothing") + y = zeros(eltype(x), size(x)) + xt = similar(x) + pns = prod(nspin) + + if ndims(x) == 1 + z = similar(x) + for i = 1:pns + shift = i - 1 + Util.circshift!(z, x, shift) + + Transforms.dwt_oop!(xt, z, wt, L) + threshold!(xt, dnt.th, sigma * dnt.t) + Transforms.idwt_oop!(z, xt, wt, L) + + Util.circshift!(xt, z, -shift) + arrayadd!(y, xt) + end + else # ndims > 1 + for i = 1:pns + shift = nspin2circ(nspin, i) + z = circshift(x, shift) + + Transforms.dwt_oop!(xt, z, wt, L) + threshold!(xt, dnt.th, sigma * dnt.t) + Transforms.idwt_oop!(z, xt, wt, L) + + z = circshift(z, -shift) + arrayadd!(y, z) + end + end + rmul!(y, 1 / pns) + else # !TI + if wt == nothing + y = copy(x) + threshold!(y, dnt.th, sigma * dnt.t) + else + y = dwt(x, wt, L) + threshold!(y, dnt.th, sigma * dnt.t) + if isa(wt, GLS) + idwt!(y, wt, L) + else + y2 = idwt(y, wt, L) + y = y2 + end + end + end + + return y +end +# add z to y +function arrayadd!(y::AbstractArray, z::AbstractArray) + length(y) == length(z) || throw(DimensionMismatch("lengths must be equal")) + for i in eachindex(y) + @inbounds y[i] += z[i] + end + return y +end + + +# estimate the std. dev. of the signal noise, assuming Gaussian distribution +function noisest(x::AbstractArray, wt::Union{DiscreteWavelet,Nothing}=DEFAULT_WAVELET, L::Integer=1) + if wt == nothing + y = x + else + y = dwt(x, wt, L) + end + dr = y[detailrange(y, L)] + return mad!(dr) / 0.6745 +end +# Median absolute deviation +function mad!(y::AbstractArray) + m = median!(y) + for i in eachindex(y) + y[i] = abs(y[i] - m) + end + return median!(y) +end + +# convert index i to a circshift array starting at 0 shift +nspin2circ(nspin::Int, i::Int) = nspin2circ((nspin,), i) +function nspin2circ(nspin::Tuple, i::Int) + c1 = CartesianIndices(nspin)[i].I + c = Vector{Int}(undef, length(c1)) + for k in 1:length(c1) + c[k] = c1[k] - 1 + end + return c +end diff --git a/src/Threshold/entropy.jl b/src/Threshold/entropy.jl new file mode 100644 index 0000000..90b947e --- /dev/null +++ b/src/Threshold/entropy.jl @@ -0,0 +1,129 @@ +using ..Transforms + + +# ENTROPY MEASURES + +abstract type Entropy end +struct ShannonEntropy <: Entropy end #Coifman-Wickerhauser +struct LogEnergyEntropy <: Entropy end + +# Entropy measures: Additive with coefentropy(0) = 0 +# all coefs assumed to be on [-1,1] after normalization with nrm +# given x and y, where x has "more concentrated energy" than y +# then coefentropy(x, et, norm) <= coefentropy(y, et, norm) should be satisfied. + +function coefentropy(x::T, et::ShannonEntropy, nrm::T) where {T<:AbstractFloat} + s = (x / nrm)^2 + if s == 0.0 + return -zero(T) + else + return -s * log(s) + end +end +function coefentropy(x::T, et::LogEnergyEntropy, nrm::T) where {T<:AbstractFloat} + s = (x / nrm)^2 + if s == 0.0 + return -zero(T) + else + return -log(s) + end +end +function coefentropy(x::AbstractArray{T}, et::Entropy, nrm::T=norm(x)) where {T<:AbstractFloat} + @assert nrm >= 0 + sum = zero(T) + nrm == sum && return sum + for i in eachindex(x) + @inbounds sum += coefentropy(x[i], et, nrm) + end + return sum +end + + +# find the best tree that is a subset of the input tree (use :full to find the best tree) +# for wpt +function bestbasistree(y::AbstractVector{T}, wt::DiscreteWavelet, L::Integer=maxtransformlevels(y), et::Entropy=ShannonEntropy()) where {T<:AbstractFloat} + bestbasistree(y, wt, maketree(length(y), L, :full), et) +end +function bestbasistree(y::AbstractVector{T}, wt::DiscreteWavelet, tree::BitVector, et::Entropy=ShannonEntropy()) where {T<:AbstractFloat} + + isvalidtree(y, tree) || throw(ArgumentError("invalid tree")) + + tree[1] || copy(tree) # do nothing + + x = copy(y) + n = length(y) + tmp = Vector{T}(undef, n) + ntree = length(tree) + entr_bf = Vector{T}(undef, ntree) + nrm = norm(y) + + Lmax = maxtransformlevels(n) + L = Lmax + k = 1 + while L > 0 + ix = 1 + Lfw = Lmax - L + nj = detailn(n, Lfw) + + @assert nj <= n + dtmp = Transforms.unsafe_vectorslice(tmp, 1, nj) + while ix <= n + @assert nj + ix - 1 <= n + dx = Transforms.unsafe_vectorslice(x, ix, nj) + + entr_bf[k] = coefentropy(dx, et, nrm) + + dwt!(dtmp, dx, wt, 1) + copyto!(dx, dtmp) + + ix += nj + k += 1 + end + L -= 1 + end + + # entropy of fully transformed signal (end nodes) + n_af = 2^(Lmax - 1) + entr_af = Vector{T}(undef, n_af) + n_coef_af = div(n, n_af) + for i in 1:n_af + range = (i-1)*n_coef_af+1:i*n_coef_af + entr_af[i] = coefentropy(x[range], et, nrm) + end + + # make the best tree + besttree = copy(tree) + for i in 1:ntree + if (i > 1 && !besttree[i>>1]) || !tree[i] # parent is 0 or input tree-node is 0 + besttree[i] = false + else + if entr_bf[i] <= bestsubtree_entropy(entr_bf, entr_af, i) + besttree[i] = false + else + besttree[i] = true + end + end + end + + @assert isvalidtree(y, besttree) + return besttree::BitVector +end + +# the entropy of best subtree +function bestsubtree_entropy(entr_bf::Array, entr_af::Array, i::Int) + n = length(entr_bf) + n_af = length(entr_af) + @assert isdyadic(n + 1) + @assert isdyadic(n_af) + @assert n + 1 == n_af << 1 + + if n < (i << 1) # bottom of tree + sum = entr_af[i-n_af+1] + else + sum = bestsubtree_entropy(entr_bf, entr_af, i << 1) + sum += bestsubtree_entropy(entr_bf, entr_af, i << 1 + 1) + end + + lowestentropy = min(entr_bf[i], sum) + return lowestentropy +end diff --git a/src/Threshold/threshold.jl b/src/Threshold/threshold.jl new file mode 100644 index 0000000..7980263 --- /dev/null +++ b/src/Threshold/threshold.jl @@ -0,0 +1,129 @@ +using ..Util +using ..WT +using ..Transforms + + +# THRESHOLD TYPES AND FUNCTIONS + +abstract type THType end +struct HardTH <: THType end +struct SoftTH <: THType end +struct SemiSoftTH <: THType end +struct SteinTH <: THType end +struct BiggestTH <: THType end +struct PosTH <: THType end +struct NegTH <: THType end + +const DEFAULT_TH = HardTH() + +# biggest m-term approximation (best m-term approximation for orthogonal transforms) +# result is m-sparse +function threshold!(x::AbstractArray{<:Number}, TH::BiggestTH, m::Int) + @assert m >= 0 + n = length(x) + m > n && (m = n) + ind = sortperm(x, alg=QuickSort, by=abs) + @inbounds begin + for i = 1:n-m + x[ind[i]] = 0 + end + end + return x +end + +# hard +function threshold!(x::AbstractArray{<:Number}, TH::HardTH, t::Real) + @assert t >= 0 + @inbounds begin + for i in eachindex(x) + if abs(x[i]) <= t + x[i] = 0 + end + end + end + return x +end + +# soft +function threshold!(x::AbstractArray{<:Number}, TH::SoftTH, t::Real) + @assert t >= 0 + @inbounds begin + for i in eachindex(x) + sh = abs(x[i]) - t + if sh < 0 + x[i] = 0 + else + x[i] = sign(x[i]) * sh + end + end + end + return x +end + +# semisoft +function threshold!(x::AbstractArray{<:Number}, TH::SemiSoftTH, t::Real) + @assert t >= 0 + @inbounds begin + for i in eachindex(x) + if x[i] <= 2 * t + sh = abs(x[i]) - t + if sh < 0 + x[i] = 0 + elseif sh - t < 0 + x[i] = sign(x[i]) * sh * 2 + end + end + end + end + return x +end + +# stein +function threshold!(x::AbstractArray{<:Number}, TH::SteinTH, t::Real) + @assert t >= 0 + @inbounds begin + for i in eachindex(x) + sh = 1 - t * t / (x[i] * x[i]) + if sh < 0 + x[i] = 0 + else + x[i] = x[i] * sh + end + end + end + return x +end + +# shrink negative elements to 0 +function threshold!(x::AbstractArray{<:Number}, TH::NegTH) + @inbounds begin + for i in eachindex(x) + if x[i] < 0 + x[i] = 0 + end + end + end + return x +end + +# shrink positive elements to 0 +function threshold!(x::AbstractArray{<:Number}, TH::PosTH) + @inbounds begin + for i in eachindex(x) + if x[i] > 0 + x[i] = 0 + end + end + end + return x +end + +# the non inplace functions +function threshold(x::AbstractArray{T}, TH::THType, t::Real) where {T<:Number} + y = Vector{T}(undef, size(x)) + return threshold!(copyto!(y, x), TH, t) +end +function threshold(x::AbstractArray{T}, TH::THType) where {T<:Number} + y = Vector{T}(undef, size(x)) + return threshold!(copyto!(y, x), TH) +end diff --git a/src/Transforms/Transforms.jl b/src/Transforms/Transforms.jl new file mode 100644 index 0000000..e17ca42 --- /dev/null +++ b/src/Transforms/Transforms.jl @@ -0,0 +1,10 @@ +module Transforms + +include("transforms.jl") +include("transforms_filter.jl") +include("transforms_lifting.jl") +include("transforms_maximal_overlap.jl") + +export + dwt, idwt, dwt!, idwt!, wpt, iwpt, wpt!, iwpt!, modwt, imodwt +end diff --git a/src/Transforms/transforms.jl b/src/Transforms/transforms.jl new file mode 100644 index 0000000..db51abc --- /dev/null +++ b/src/Transforms/transforms.jl @@ -0,0 +1,226 @@ +using ..Util +using ..WT +using FFTW + +# TODO Use StridedArray instead of AbstractArray where writing to array. +# TODO change integer dependent wavelets to parametric types (see "Value types", https://docs.julialang.org/en/v1/manual/types/index.html#%22Value-types%22-1) +const DWTArray = AbstractArray +const WPTArray = AbstractVector +const ValueType = Union{AbstractFloat,Complex} + +const FVector = StridedVector # e.g., work space vectors + +# DWT + +""" + dwt(x, wt[, L=maxtransformlevels(x)]) + +Perform a discrete wavelet transform of the array `x`. +The wavelet type `wt` determines the transform type +(filter or lifting) and the wavelet class, see `wavelet`. + +The number of transformation levels `L` can be any non-negative +integer such that the size of `x` is divisible by `L`. +Performs the identity transformation if `L==0`. + +# Examples +```julia +dwt(x, wavelet(WT.coif6)) +``` + +**See also:** `idwt`, `dwt!`, `wavelet` +""" +function dwt end + +""" + idwt(x, wt[, L=maxtransformlevels(x)]) + +Perform an inverse discrete wavelet transform of the array `x`, +the inverse of `dwt(x, wt, L)`. + +**See also:** `dwt`, `idwt!` +""" +function idwt end + +""" + dwt!(y, x, wt::OrthoFilter[, L=maxtransformlevels(x)]) + + dwt!(y, wt::GLS[, L=maxtransformlevels(x)]) + +Same as `dwt` but without array allocation. +Perform "out of place" transform with a filter, or +a inplace transform with a lifting scheme. The difference +between the filter and lifting methods is due to the +structure of the transform algorithms. + +**See also:** `idwt!` +""" +function dwt! end + +""" + idwt!(y, x, wt::OrthoFilter[, L=maxtransformlevels(x)]) + + idwt!(y, wt::GLS[, L=maxtransformlevels(x)]) + +The inverse of `dwt!`. + +**See also:** `dwt!` +""" +function idwt! end + + +# WPT + +""" + wpt + +Perform a discrete wavelet packet transform of the array `x`. +**See also:** `dwt`, `wavelet` +""" +function wpt end + +""" + iwpt + +Perform an inverse discrete wavelet packet transform of the array `x`. +**See also:** `idwt`, `wavelet` +""" +function iwpt end + +""" + wpt! + +Same as `wpt` but without array allocation. +""" +function wpt! end + +""" + iwpt! + +Same as `iwpt` but without array allocation. +""" +function iwpt! end + + +# DWT (discrete wavelet transform) + +for (Xwt, Xwt!, _Xwt!, fw) in ((:dwt, :dwt!, :_dwt!, true), + (:idwt, :idwt!, :_dwt!, false)) + @eval begin + # filter + function ($Xwt)(x::DWTArray{T}, filter::OrthoFilter, + L::Integer=maxtransformlevels(x)) where {T<:ValueType} + y = Array{T}(undef, size(x)) + return ($_Xwt!)(y, x, filter, L, $fw) + end + function ($Xwt!)(y::DWTArray{<:ValueType}, x::DWTArray{<:ValueType}, filter::OrthoFilter, + L::Integer=maxtransformlevels(x)) + return ($_Xwt!)(y, x, filter, L, $fw) + end + # lifting + function ($Xwt)(x::DWTArray{T}, scheme::GLS, + L::Integer=maxtransformlevels(x)) where {T<:ValueType} + y = Array{T}(undef, size(x)) + copyto!(y, x) + return ($_Xwt!)(y, scheme, L, $fw) + end + function ($Xwt!)(y::DWTArray{T}, scheme::GLS, + L::Integer=maxtransformlevels(y)) where {T<:ValueType} + return ($_Xwt!)(y, scheme, L, $fw) + end + end # begin +end # for + +# WPT (wavelet packet transform) + +for (Xwt, Xwt!, _Xwt!, fw) in ((:wpt, :wpt!, :_wpt!, true), + (:iwpt, :iwpt!, :_wpt!, false)) + @eval begin + function ($Xwt)(x::WPTArray{T}, wt::DiscreteWavelet, + L::Integer=maxtransformlevels(x)) where {T<:ValueType} + return ($Xwt)(x, wt, maketree(length(x), L, :full)) + end + # filter + function ($Xwt)(x::WPTArray{T}, filter::OrthoFilter, + tree::BitVector=maketree(x, :full)) where {T<:ValueType} + y = Array{T}(undef, size(x)) + return ($_Xwt!)(y, x, filter, tree, $fw) + end + function ($Xwt!)(y::WPTArray{T}, x::WPTArray{T}, + filter::OrthoFilter) where {T<:ValueType} + return ($Xwt!)(y, x, filter, maketree(x, :full)) + end + function ($Xwt!)(y::WPTArray{T}, x::WPTArray{T}, + filter::OrthoFilter, L::Integer) where {T<:ValueType} + return ($Xwt!)(y, x, filter, maketree(length(x), L, :full)) + end + function ($Xwt!)(y::WPTArray{T}, x::WPTArray{T}, + filter::OrthoFilter, tree::BitVector) where {T<:ValueType} + return ($_Xwt!)(y, x, filter, tree, $fw) + end + # lifting + function ($Xwt)(x::WPTArray{T}, scheme::GLS, + tree::BitVector=maketree(x, :full)) where {T<:ValueType} + y = Array{T}(undef, size(x)) + copyto!(y, x) + return ($_Xwt!)(y, scheme, tree, $fw) + end + function ($Xwt!)(y::WPTArray{T}, scheme::GLS) where {T<:ValueType} + return ($Xwt!)(y, scheme, maketree(y, :full)) + end + function ($Xwt!)(y::WPTArray{T}, scheme::GLS, L::Integer) where {T<:ValueType} + return ($Xwt!)(y, scheme, maketree(length(x), L, :full)) + end + function ($Xwt!)(y::WPTArray{T}, scheme::GLS, tree::BitVector) where {T<:ValueType} + return ($_Xwt!)(y, scheme, tree, $fw) + end + end # begin +end # for + + +# DWTC (column-wise discrete wavelet transform) +#dwtc(::AbstractArray, ::DiscreteWavelet) +#idwtc(::AbstractArray, ::DiscreteWavelet) + +# SWT (stationary wavelet transform) +# swt(::AbstractVector, ::DiscreteWavelet) +# iswt(::AbstractVector, ::DiscreteWavelet) + +# Int -> Float +for Xwt in (:dwt, :idwt, :dwtc, :idwtc, :wpt, :iwpt) + @eval $Xwt(x::AbstractArray{<:Integer}, args...) = $Xwt(float(x), args...) +end + +# non-exported "out of place" functions +for (Xwt_oop!, Xwt!) in ((:dwt_oop!, :dwt!), (:idwt_oop!, :idwt!)) + @eval begin + # filter + function ($Xwt_oop!)(y::DWTArray{T}, x::DWTArray{T}, filter::OrthoFilter, + L::Integer=maxtransformlevels(x)) where {T<:ValueType} + return ($Xwt!)(y, x, filter, L) + end + # lifting + function ($Xwt_oop!)(y::DWTArray{T}, x::DWTArray{T}, scheme::GLS, + L::Integer=maxtransformlevels(x)) where {T<:ValueType} + copyto!(y, x) + return ($Xwt!)(y, scheme, L) + end + end # begin +end # for + +# Array with shared memory +function unsafe_vectorslice(A::Array{T}, i::Int, n::Int) where {T} + return unsafe_wrap(Array, pointer(A, i), n)::Vector{T} +end +function unsafe_vectorslice(A::StridedArray{T}, i::Int, n::Int) where {T} + return @view A[i:(i-1+n)] +end + +# linear indices of start of rows/cols/planes +# 2-D, size(A) = (m,n) +row_idx(i, m) = i +col_idx(i, m) = 1 + (i - 1) * m +# 3-D, size(A) = (m,n,d) +row_idx(i, j, m, n=m) = row_idx(i, n) + (j - 1) * n * m +col_idx(i, j, m, n=m) = col_idx(i, m) + (j - 1) * n * m +plane_idx(i, j, m) = i + (j - 1) * m diff --git a/src/mod/transforms_filter.jl b/src/Transforms/transforms_filter.jl similarity index 65% rename from src/mod/transforms_filter.jl rename to src/Transforms/transforms_filter.jl index 1c1c074..75b966e 100644 --- a/src/mod/transforms_filter.jl +++ b/src/Transforms/transforms_filter.jl @@ -11,16 +11,16 @@ # 1-D # writes to y function _dwt!(y::AbstractVector{Ty}, x::AbstractVector{Tx}, - filter::OrthoFilter, L::Integer, fw::Bool) where {Tx<:Number, Ty<:Number} + filter::OrthoFilter, L::Integer, fw::Bool) where {Tx<:Number,Ty<:Number} T = promote_type(Tx, Ty) - si = Vector{T}(undef, length(filter)-1) # tmp filter vector + si = Vector{T}(undef, length(filter) - 1) # tmp filter vector scfilter, dcfilter = WT.makereverseqmfpair(filter, fw, T) return _dwt!(y, x, filter, L, fw, dcfilter, scfilter, si) end function _dwt!(y::AbstractVector{<:Number}, x::AbstractVector{<:Number}, - filter::OrthoFilter, L::Integer, fw::Bool, - dcfilter::Vector{T}, scfilter::Vector{T}, - si::Vector{T}, snew::Vector{T} = Vector{T}(undef, ifelse(L>1, length(x)>>1, 0))) where T<:Number + filter::OrthoFilter, L::Integer, fw::Bool, + dcfilter::Vector{T}, scfilter::Vector{T}, + si::Vector{T}, snew::Vector{T}=Vector{T}(undef, ifelse(L > 1, length(x) >> 1, 0))) where {T<:Number} n = length(x) size(x) == size(y) || throw(DimensionMismatch("in and out array size must match")) @@ -30,11 +30,11 @@ function _dwt!(y::AbstractVector{<:Number}, x::AbstractVector{<:Number}, throw(ArgumentError("size must have a sufficient power of 2 factor")) y === x && throw(ArgumentError("in array is out array")) - L > 1 && length(snew) != n>>1 && + L > 1 && length(snew) != n >> 1 && throw(ArgumentError("length of snew incorrect")) if L == 0 - return copyto!(y,x) + return copyto!(y, x) end s = x # s is current scaling coefs location filtlen = length(filter) @@ -45,49 +45,49 @@ function _dwt!(y::AbstractVector{<:Number}, x::AbstractVector{<:Number}, for l in lrange if fw # detail coefficients - filtdown!(dcfilter, si, y, detailindex(n,l,1), detailn(n,l), s, 1,-filtlen+1, true) + filtdown!(dcfilter, si, y, detailindex(n, l, 1), detailn(n, l), s, 1, -filtlen + 1, true) # scaling coefficients - filtdown!(scfilter, si, y, 1, detailn(n,l), s, 1, 0, false) + filtdown!(scfilter, si, y, 1, detailn(n, l), s, 1, 0, false) else # scaling coefficients - filtup!(false, scfilter, si, y, 1, detailn(n,l-1), s, 1, -filtlen+1, false) + filtup!(false, scfilter, si, y, 1, detailn(n, l - 1), s, 1, -filtlen + 1, false) # detail coefficients - filtup!(true, dcfilter, si, y, 1, detailn(n,l-1), x, detailindex(n,l,1), 0, true) + filtup!(true, dcfilter, si, y, 1, detailn(n, l - 1), x, detailindex(n, l, 1), 0, true) end # if not final iteration: copy to tmp location - l != lrange[end] && copyto!(snew,1,y,1,detailn(n, fw ? l : l-1)) + l != lrange[end] && copyto!(snew, 1, y, 1, detailn(n, fw ? l : l - 1)) L > 1 && (s = snew) end return y end function unsafe_dwt1level!(y::AbstractVector{<:Number}, x::AbstractVector{<:Number}, - filter::OrthoFilter, fw::Bool, - dcfilter::FVector{T}, scfilter::FVector{T}, - si::FVector{T}) where T<:Number + filter::OrthoFilter, fw::Bool, + dcfilter::FVector{T}, scfilter::FVector{T}, + si::FVector{T}) where {T<:Number} n = length(x) l = 1 filtlen = length(filter) if fw # detail coefficients - filtdown!(dcfilter, si, y, detailindex(n,l,1), detailn(n,l), x, 1,-filtlen+1, true) + filtdown!(dcfilter, si, y, detailindex(n, l, 1), detailn(n, l), x, 1, -filtlen + 1, true) # scaling coefficients - filtdown!(scfilter, si, y, 1, detailn(n,l), x, 1, 0, false) + filtdown!(scfilter, si, y, 1, detailn(n, l), x, 1, 0, false) else # scaling coefficients - filtup!(false, scfilter, si, y, 1, detailn(n,l-1), x, 1, -filtlen+1, false) + filtup!(false, scfilter, si, y, 1, detailn(n, l - 1), x, 1, -filtlen + 1, false) # detail coefficients - filtup!(true, dcfilter, si, y, 1, detailn(n,l-1), x, detailindex(n,l,1), 0, true) + filtup!(true, dcfilter, si, y, 1, detailn(n, l - 1), x, detailindex(n, l, 1), 0, true) end return y end function dwt_transform_strided!(y::AbstractArray{<:Number}, x::AbstractArray{<:Number}, - msub::Int, nsub::Int, stride::Int, idx_func::Function, - tmpvec::FVector{T}, tmpvec2::FVector{T}, - filter::OrthoFilter, fw::Bool, - dcfilter::FVector{T}, scfilter::FVector{T}, si::FVector{T}) where T<:Number - for i=1:msub + msub::Int, nsub::Int, stride::Int, idx_func::Function, + tmpvec::FVector{T}, tmpvec2::FVector{T}, + filter::OrthoFilter, fw::Bool, + dcfilter::FVector{T}, scfilter::FVector{T}, si::FVector{T}) where {T<:Number} + for i = 1:msub xi = idx_func(i) stridedcopy!(tmpvec, x, xi, stride, nsub) unsafe_dwt1level!(tmpvec2, tmpvec, filter, fw, dcfilter, scfilter, si) @@ -96,11 +96,11 @@ function dwt_transform_strided!(y::AbstractArray{<:Number}, x::AbstractArray{<:N end function dwt_transform_cols!(y::AbstractArray{<:Number}, x::AbstractArray{<:Number}, - msub::Int, nsub::Int, idx_func::Function, - tmpvec::FVector{T}, - filter::OrthoFilter, fw::Bool, - dcfilter::FVector{T}, scfilter::FVector{T}, si::FVector{T}) where T<:Number - for i=1:nsub + msub::Int, nsub::Int, idx_func::Function, + tmpvec::FVector{T}, + filter::OrthoFilter, fw::Bool, + dcfilter::FVector{T}, scfilter::FVector{T}, si::FVector{T}) where {T<:Number} + for i = 1:nsub xi = idx_func(i) copyto!(tmpvec, 1, x, xi, msub) ya = unsafe_vectorslice(y, xi, msub) @@ -111,19 +111,19 @@ end # 2-D # writes to y function _dwt!(y::AbstractMatrix{Ty}, x::AbstractMatrix{Tx}, - filter::OrthoFilter, L::Integer, fw::Bool) where {Tx<:Number, Ty<:Number} + filter::OrthoFilter, L::Integer, fw::Bool) where {Tx<:Number,Ty<:Number} m, n = size(x) T = promote_type(Tx, Ty) - si = Vector{T}(undef, length(filter)-1) # tmp filter vector - tmpbuffer = Vector{T}(undef, max(n<<1, m)) # tmp storage vector + si = Vector{T}(undef, length(filter) - 1) # tmp filter vector + tmpbuffer = Vector{T}(undef, max(n << 1, m)) # tmp storage vector scfilter, dcfilter = WT.makereverseqmfpair(filter, fw, T) return _dwt!(y, x, filter, L, fw, dcfilter, scfilter, si, tmpbuffer) end function _dwt!(y::AbstractMatrix{<:Number}, x::AbstractMatrix{<:Number}, - filter::OrthoFilter, L::Integer, fw::Bool, - dcfilter::Vector{T}, scfilter::Vector{T}, - si::Vector{T}, tmpbuffer::Vector{T}) where T<:Number + filter::OrthoFilter, L::Integer, fw::Bool, + dcfilter::Vector{T}, scfilter::Vector{T}, + si::Vector{T}, tmpbuffer::Vector{T}) where {T<:Number} m, n = size(x) size(x) == size(y) || @@ -134,11 +134,11 @@ function _dwt!(y::AbstractMatrix{<:Number}, x::AbstractMatrix{<:Number}, throw(ArgumentError("size must have a sufficient power of 2 factor")) y === x && throw(ArgumentError("in array is out array")) - length(tmpbuffer) >= max(n<<1,m) || + length(tmpbuffer) >= max(n << 1, m) || throw(ArgumentError("length of tmpbuffer incorrect")) if L == 0 - return copyto!(y,x) + return copyto!(y, x) end row_stride = m #s = x @@ -149,9 +149,9 @@ function _dwt!(y::AbstractMatrix{<:Number}, x::AbstractMatrix{<:Number}, msub = m else lrange = L:-1:1 - nsub = div(n,2^(L-1)) - msub = div(m,2^(L-1)) - copyto!(y,x) + nsub = div(n, 2^(L - 1)) + msub = div(m, 2^(L - 1)) + copyto!(y, x) end inputArray = x @@ -159,9 +159,9 @@ function _dwt!(y::AbstractMatrix{<:Number}, x::AbstractMatrix{<:Number}, row_idx_func = i -> row_idx(i, m) col_idx_func = i -> col_idx(i, m) for l in lrange - tmpvec = unsafe_vectorslice(tmpbuffer, 1, msub) + tmpvec = unsafe_vectorslice(tmpbuffer, 1, msub) tmpvec2 = unsafe_vectorslice(tmpbuffer, 1, nsub) - tmpvec3 = unsafe_vectorslice(tmpbuffer, nsub+1, nsub) + tmpvec3 = unsafe_vectorslice(tmpbuffer, nsub + 1, nsub) if fw # rows dwt_transform_strided!(y, inputArray, msub, nsub, row_stride, row_idx_func, @@ -181,28 +181,28 @@ function _dwt!(y::AbstractMatrix{<:Number}, x::AbstractMatrix{<:Number}, dwt_transform_strided!(y, y, msub, nsub, row_stride, row_idx_func, tmpvec2, tmpvec3, filter, fw, dcfilter, scfilter, si) end - msub = (fw ? msub>>1 : msub<<1) - nsub = (fw ? nsub>>1 : nsub<<1) + msub = (fw ? msub >> 1 : msub << 1) + nsub = (fw ? nsub >> 1 : nsub << 1) end return y end # 3-D # writes to y -function _dwt!(y::AbstractArray{Ty, 3}, x::AbstractArray{Tx, 3}, - filter::OrthoFilter, L::Integer, fw::Bool) where {Tx<:Number, Ty<:Number} +function _dwt!(y::AbstractArray{Ty,3}, x::AbstractArray{Tx,3}, + filter::OrthoFilter, L::Integer, fw::Bool) where {Tx<:Number,Ty<:Number} m, n, d = size(x) T = promote_type(Tx, Ty) - si = Vector{T}(undef, length(filter)-1) # tmp filter vector - tmpbuffer = Vector{T}(undef, max(m, n<<1, d<<1)) # tmp storage vector + si = Vector{T}(undef, length(filter) - 1) # tmp filter vector + tmpbuffer = Vector{T}(undef, max(m, n << 1, d << 1)) # tmp storage vector scfilter, dcfilter = WT.makereverseqmfpair(filter, fw, T) return _dwt!(y, x, filter, L, fw, dcfilter, scfilter, si, tmpbuffer) end -function _dwt!(y::AbstractArray{<:Number, 3}, x::AbstractArray{<:Number, 3}, - filter::OrthoFilter, L::Integer, fw::Bool, - dcfilter::Vector{T}, scfilter::Vector{T}, - si::Vector{T}, tmpbuffer::Vector{T}) where T<:Number +function _dwt!(y::AbstractArray{<:Number,3}, x::AbstractArray{<:Number,3}, + filter::OrthoFilter, L::Integer, fw::Bool, + dcfilter::Vector{T}, scfilter::Vector{T}, + si::Vector{T}, tmpbuffer::Vector{T}) where {T<:Number} m, n, d = size(x) size(x) == size(y) || @@ -213,14 +213,14 @@ function _dwt!(y::AbstractArray{<:Number, 3}, x::AbstractArray{<:Number, 3}, throw(ArgumentError("size must have a sufficient power of 2 factor")) y === x && throw(ArgumentError("in array is out array")) - length(tmpbuffer) >= n<<1 || + length(tmpbuffer) >= n << 1 || throw(ArgumentError("length of tmpbuffer incorrect")) if L == 0 - return copyto!(y,x) + return copyto!(y, x) end row_stride = m - plane_stride = m*n + plane_stride = m * n if fw lrange = 1:L @@ -229,20 +229,20 @@ function _dwt!(y::AbstractArray{<:Number, 3}, x::AbstractArray{<:Number, 3}, dsub = d else lrange = L:-1:1 - msub = div(m,2^(L-1)) - nsub = div(n,2^(L-1)) - dsub = div(d,2^(L-1)) - copyto!(y,x) + msub = div(m, 2^(L - 1)) + nsub = div(n, 2^(L - 1)) + dsub = div(d, 2^(L - 1)) + copyto!(y, x) end inputArray = x for l in lrange - tmpcol = unsafe_vectorslice(tmpbuffer, 1, msub) - tmprow = unsafe_vectorslice(tmpbuffer, 1, nsub) - tmprow2 = unsafe_vectorslice(tmpbuffer, nsub+1, nsub) - tmphei = unsafe_vectorslice(tmpbuffer, 1, dsub) - tmphei2 = unsafe_vectorslice(tmpbuffer, dsub+1, dsub) + tmpcol = unsafe_vectorslice(tmpbuffer, 1, msub) + tmprow = unsafe_vectorslice(tmpbuffer, 1, nsub) + tmprow2 = unsafe_vectorslice(tmpbuffer, nsub + 1, nsub) + tmphei = unsafe_vectorslice(tmpbuffer, 1, dsub) + tmphei2 = unsafe_vectorslice(tmpbuffer, dsub + 1, dsub) if fw # planes for j in 1:nsub @@ -286,9 +286,9 @@ function _dwt!(y::AbstractArray{<:Number, 3}, x::AbstractArray{<:Number, 3}, tmphei, tmphei2, filter, fw, dcfilter, scfilter, si) end end - msub = (fw ? msub>>1 : msub<<1) - nsub = (fw ? nsub>>1 : nsub<<1) - dsub = (fw ? dsub>>1 : dsub<<1) + msub = (fw ? msub >> 1 : msub << 1) + nsub = (fw ? nsub >> 1 : nsub << 1) + dsub = (fw ? dsub >> 1 : dsub << 1) end return y end @@ -298,15 +298,15 @@ end # WPT # 1-D # writes to y -function _wpt!(y::AbstractVector{T}, x::AbstractVector{T}, filter::OrthoFilter, tree::BitVector, fw::Bool) where T<:Number - si = Vector{T}(undef, length(filter)-1) - ns = ifelse(fw, length(x)>>1, length(x)) +function _wpt!(y::AbstractVector{T}, x::AbstractVector{T}, filter::OrthoFilter, tree::BitVector, fw::Bool) where {T<:Number} + si = Vector{T}(undef, length(filter) - 1) + ns = ifelse(fw, length(x) >> 1, length(x)) snew = Vector{T}(undef, ns) scfilter, dcfilter = WT.makereverseqmfpair(filter, fw, T) return _wpt!(y, x, filter, tree, fw, dcfilter, scfilter, si, snew) end -function _wpt!(y::AbstractVector{T}, x::AbstractVector{T}, filter::OrthoFilter, tree::BitVector, fw::Bool, dcfilter::Vector{T}, scfilter::Vector{T}, si::Vector{T}, snew::Vector{T}) where T<:Number +function _wpt!(y::AbstractVector{T}, x::AbstractVector{T}, filter::OrthoFilter, tree::BitVector, fw::Bool, dcfilter::Vector{T}, scfilter::Vector{T}, si::Vector{T}, snew::Vector{T}) where {T<:Number} size(x) == size(y) || throw(DimensionMismatch("in and out array size must match")) @@ -314,12 +314,12 @@ function _wpt!(y::AbstractVector{T}, x::AbstractVector{T}, filter::OrthoFilter, throw(ArgumentError("in array is out array")) isvalidtree(y, tree) || throw(ArgumentError("invalid tree")) - if tree[1] && length(snew) < ifelse(fw, length(x)>>1, length(x)) + if tree[1] && length(snew) < ifelse(fw, length(x) >> 1, length(x)) throw(ArgumentError("length of snew incorrect")) end if !tree[1] - return copyto!(y,x) + return copyto!(y, x) end first = true @@ -329,9 +329,9 @@ function _wpt!(y::AbstractVector{T}, x::AbstractVector{T}, filter::OrthoFilter, while L > 0 ix = 1 k = 1 - Lfw = (fw ? Lmax-L : L-1) + Lfw = (fw ? Lmax - L : L - 1) nj = detailn(n, Lfw) - treeind = 2^(Lfw)-1 + treeind = 2^(Lfw) - 1 dx = first ? x : unsafe_vectorslice(snew, 1, nj) # dx will be overwritten if first while ix <= n @@ -361,15 +361,15 @@ end macro filtermainloop(si, silen, b, val) quote - @inbounds for j=2:$(esc(silen)) - $(esc(si))[j-1] = $(esc(si))[j] + $(esc(b))[j]*$(esc(val)) + @inbounds for j = 2:$(esc(silen)) + $(esc(si))[j-1] = $(esc(si))[j] + $(esc(b))[j] * $(esc(val)) end - @inbounds $(esc(si))[$(esc(silen))] = $(esc(b))[$(esc(silen))+1]*$(esc(val)) + @inbounds $(esc(si))[$(esc(silen))] = $(esc(b))[$(esc(silen))+1] * $(esc(val)) end end macro filtermainloopzero(si, silen) quote - @inbounds for j=2:$(esc(silen)) + @inbounds for j = 2:$(esc(silen)) $(esc(si))[j-1] = $(esc(si))[j] end @inbounds $(esc(si))[$(esc(silen))] = 0.0 @@ -385,45 +385,45 @@ end # ss : shift downsampling # based on Base.filt function filtdown!(f::AbstractVector{T}, si::AbstractVector{T}, - out::AbstractVector{<:Number}, iout::Integer, nout::Integer, - x::AbstractVector{<:Number}, ix::Integer, - shift::Integer=0, ss::Bool=false) where T<:Number - nx = nout<<1 + out::AbstractVector{<:Number}, iout::Integer, nout::Integer, + x::AbstractVector{<:Number}, ix::Integer, + shift::Integer=0, ss::Bool=false) where {T<:Number} + nx = nout << 1 silen = length(si) flen = length(f) - @assert length(x) >= ix+nx-1 - @assert length(out) >= iout+nout-1 - @assert (nx-1)>>1 + iout <= length(out) - @assert silen == flen-1 + @assert length(x) >= ix + nx - 1 + @assert length(out) >= iout + nout - 1 + @assert (nx - 1) >> 1 + iout <= length(out) + @assert silen == flen - 1 @assert shift <= 0 - fill!(si,0.0) + fill!(si, 0.0) istart = flen + Int(ss) - dsshift = (flen%2 + Int(ss))%2 # is flen odd, and shift downsampling + dsshift = (flen % 2 + Int(ss)) % 2 # is flen odd, and shift downsampling rout1, rin, rout2 = splitdownrangeper(istart, ix, nx, shift) @inbounds begin for i in rout1 # rout1 assumed to be in 1:istart-1 # periodic in the range [ix:ix+nx-1] - xatind = x[mod(i-1+shift, nx) + ix] + xatind = x[mod(i - 1 + shift, nx)+ix] @filtermainloop(si, silen, f, xatind) end # rin is inbounds in [ix:ix+nx-1] ixsh = -1 + shift + ix for i in rin - xatind = x[i + ixsh] + xatind = x[i+ixsh] # take every other value after istart - if (i+dsshift)%2 == 0 && i >= istart - out[(i-istart)>>1 + iout] = si[1] + f[1]*xatind + if (i + dsshift) % 2 == 0 && i >= istart + out[(i-istart)>>1+iout] = si[1] + f[1] * xatind end @filtermainloop(si, silen, f, xatind) end for i in rout2 # periodic in the range [ix:ix+nx-1] - xatind = x[mod(i-1+shift, nx) + ix] + xatind = x[mod(i - 1 + shift, nx)+ix] # take every other value after istart - if (i+dsshift)%2 == 0 && i >= istart - out[(i-istart)>>1 + iout] = si[1] + f[1]*xatind + if (i + dsshift) % 2 == 0 && i >= istart + out[(i-istart)>>1+iout] = si[1] + f[1] * xatind end @filtermainloop(si, silen, f, xatind) end @@ -438,14 +438,14 @@ function splitdownrangeper(istart, ix, nx, shift) ixsh = -1 + shift + ix if mod(shift, nx) + ix == 1 + ixsh # shift likely 0 inxi = 1 - iend = nx-1 + iend = nx - 1 while mod(iend - 1 + shift, nx) + ix != iend + ixsh iend -= 1 end return (0:-1, inxi:iend, (iend+1):(nx-1+istart)) elseif mod(istart - 1 + shift, nx) + ix == istart + ixsh inxi = istart - iend = nx-1+istart + iend = nx - 1 + istart while mod(iend - 1 + shift, nx) + ix != iend + ixsh iend -= 1 end @@ -465,71 +465,71 @@ end # ss : shift upsampling # based on Base.filt function filtup!(add2out::Bool, f::Vector{T}, si::Vector{T}, - out::AbstractVector{<:Number}, iout::Integer, nout::Integer, - x::AbstractVector{<:Number}, ix::Integer, - shift::Integer=0, ss::Bool=false) where T<:Number - nx = nout>>1 + out::AbstractVector{<:Number}, iout::Integer, nout::Integer, + x::AbstractVector{<:Number}, ix::Integer, + shift::Integer=0, ss::Bool=false) where {T<:Number} + nx = nout >> 1 silen = length(si) flen = length(f) - @assert length(x) >= ix+nx-1 # check array size - @assert length(out) >= iout+nout-1 # check array size - @assert (nx<<1-1)>>1 + iout <= length(out) # max for out index - @assert silen == flen-1 + @assert length(x) >= ix + nx - 1 # check array size + @assert length(out) >= iout + nout - 1 # check array size + @assert (nx << 1 - 1) >> 1 + iout <= length(out) # max for out index + @assert silen == flen - 1 @assert shift <= 0 - fill!(si,0.0) - istart = flen - shift%2 - dsshift = Int(ss)%2 # shift upsampling + fill!(si, 0.0) + istart = flen - shift % 2 + dsshift = Int(ss) % 2 # shift upsampling rout1, rin, rout2 = splituprangeper(istart, ix, nx, nout, shift) @inbounds begin xatind = 0.0 for i in rout1 # rout1 assumed to be in 1:istart-1 - if (i+dsshift)%2 == 0 + if (i + dsshift) % 2 == 0 @filtermainloopzero(si, silen) else # periodic in the range [ix:ix+nx-1] - xindex = mod((i-1)>>1+shift>>1,nx) + ix #(i-1)>>1 increm. every other + xindex = mod((i - 1) >> 1 + shift >> 1, nx) + ix #(i-1)>>1 increm. every other xatind = x[xindex] @filtermainloop(si, silen, f, xatind) end end - ixsh = shift>>1 + ix + ixsh = shift >> 1 + ix for i in rin - if (i+dsshift)%2 == 0 + if (i + dsshift) % 2 == 0 xatind = 0.0 else - xindex = (i-1)>>1 + ixsh + xindex = (i - 1) >> 1 + ixsh xatind = x[xindex] end if i >= istart if add2out - out[(i-istart) + iout] += si[1] + f[1]*xatind + out[(i-istart)+iout] += si[1] + f[1] * xatind else - out[(i-istart) + iout] = si[1] + f[1]*xatind + out[(i-istart)+iout] = si[1] + f[1] * xatind end end - if (i+dsshift)%2 == 0 + if (i + dsshift) % 2 == 0 @filtermainloopzero(si, silen) else @filtermainloop(si, silen, f, xatind) end end for i in rout2 - if (i+dsshift)%2 == 0 + if (i + dsshift) % 2 == 0 xatind = 0.0 else - xindex = mod((i-1)>>1+shift>>1,nx) + ix + xindex = mod((i - 1) >> 1 + shift >> 1, nx) + ix xatind = x[xindex] end if i >= istart if add2out - out[(i-istart) + iout] += si[1] + f[1]*xatind + out[(i-istart)+iout] += si[1] + f[1] * xatind else - out[(i-istart) + iout] = si[1] + f[1]*xatind + out[(i-istart)+iout] = si[1] + f[1] * xatind end end - if (i+dsshift)%2 == 0 + if (i + dsshift) % 2 == 0 @filtermainloopzero(si, silen) else @filtermainloop(si, silen, f, xatind) @@ -543,18 +543,18 @@ end # where mod((i-1)>>1+shift>>1, nx) + ix == (i-1)>>1 + ixsh function splituprangeper(istart, ix, nx, nout, shift) inxi = 0 - ixsh = shift>>1 + ix - if mod(shift>>1, nx) + ix == ixsh # shift likely 0 + ixsh = shift >> 1 + ix + if mod(shift >> 1, nx) + ix == ixsh # shift likely 0 inxi = 1 - iend = nout-1 - while mod((iend-1)>>1+shift>>1, nx) + ix != (iend-1)>>1 + ixsh + iend = nout - 1 + while mod((iend - 1) >> 1 + shift >> 1, nx) + ix != (iend - 1) >> 1 + ixsh iend -= 1 end return (0:-1, inxi:iend, (iend+1):(nout-1+istart)) - elseif mod((istart-1)>>1+shift>>1, nx) + ix == (istart-1)>>1 + ixsh + elseif mod((istart - 1) >> 1 + shift >> 1, nx) + ix == (istart - 1) >> 1 + ixsh inxi = istart - iend = nout-1+istart - while mod((iend-1)>>1+shift>>1, nx) + ix != (iend-1)>>1 + ixsh + iend = nout - 1 + istart + while mod((iend - 1) >> 1 + shift >> 1, nx) + ix != (iend - 1) >> 1 + ixsh iend -= 1 end return (1:inxi-1, inxi:iend, (iend+1):(nout-1+istart)) diff --git a/src/mod/transforms_lifting.jl b/src/Transforms/transforms_lifting.jl similarity index 69% rename from src/mod/transforms_lifting.jl rename to src/Transforms/transforms_lifting.jl index 5d4b56c..eeaa00d 100644 --- a/src/mod/transforms_lifting.jl +++ b/src/Transforms/transforms_lifting.jl @@ -7,20 +7,20 @@ ################################################################################## # for split and merge -reqtmplength(x::AbstractArray) = (size(x,1)>>2) + (size(x,1)>>1)%2 +reqtmplength(x::AbstractArray) = (size(x, 1) >> 2) + (size(x, 1) >> 1) % 2 # return scheme parameters adjusted for direction and type -function makescheme(::Type{T}, scheme::GLS, fw::Bool) where T<:Number +function makescheme(::Type{T}, scheme::GLS, fw::Bool) where {T<:Number} n = length(scheme.step) stepseq = Vector{WT.LSStep{T}}(undef, n) for i = 1:n - j = fw ? i : n+1-i + j = fw ? i : n + 1 - i stepseq[i] = WT.LSStep(scheme.step[j].steptype, - convert(Vector{T}, scheme.step[j].param.coef * (fw ? -1 : 1)), - scheme.step[j].param.shift) + convert(Vector{T}, scheme.step[j].param.coef * (fw ? -1 : 1)), + scheme.step[j].param.shift) end - norm1, norm2 = convert(T, fw ? scheme.norm1 : 1/scheme.norm1), - convert(T, fw ? scheme.norm2 : 1/scheme.norm2) + norm1, norm2 = convert(T, fw ? scheme.norm1 : 1 / scheme.norm1), + convert(T, fw ? scheme.norm2 : 1 / scheme.norm2) return stepseq, norm1, norm2 end @@ -28,14 +28,14 @@ end # 1-D # inplace transform of y, no vector allocation function _dwt!(y::AbstractVector{T}, scheme::GLS, L::Integer, fw::Bool, - tmp::Vector{T} = Vector{T}(undef, reqtmplength(y))) where T<:Number + tmp::Vector{T}=Vector{T}(undef, reqtmplength(y))) where {T<:Number} n = length(y) 0 <= L || throw(ArgumentError("L must be positive")) sufficientpoweroftwo(y, L) || throw(ArgumentError("size must have a sufficient power of 2 factor")) - length(tmp) >= n>>2 || + length(tmp) >= n >> 2 || throw(ArgumentError("length of tmp incorrect")) if L == 0 @@ -47,9 +47,9 @@ function _dwt!(y::AbstractVector{T}, scheme::GLS, L::Integer, fw::Bool, ns = n else lrange = L:-1:1 - ns = detailn(n, L-1) + ns = detailn(n, L - 1) end - half = ns>>1 + half = ns >> 1 s = y stepseq, norm1, norm2 = makescheme(T, scheme, fw) @@ -60,16 +60,16 @@ function _dwt!(y::AbstractVector{T}, scheme::GLS, L::Integer, fw::Bool, lift!(s, half, step.param, step.steptype) end normalize!(s, half, ns, norm1, norm2) - ns = ns>>1 - half = half>>1 + ns = ns >> 1 + half = half >> 1 else normalize!(s, half, ns, norm1, norm2) for step in stepseq lift!(s, half, step.param, step.steptype) end Util.merge!(s, ns, tmp) # inverse split - ns = ns<<1 - half = half<<1 + ns = ns << 1 + half = half << 1 end end return y @@ -81,12 +81,12 @@ end # oopv: the out of place location function unsafe_dwt1level!(y::AbstractArray{T}, iy::Integer, incy::Integer, oopc::Bool, oopv::FVector{T}, scheme::GLS, fw::Bool, stepseq::FVector, norm1::T, norm2::T, - tmp::FVector{T}) where T<:Number + tmp::FVector{T}) where {T<:Number} if !oopc oopv = y end ns = length(oopv) - half = ns>>1 + half = ns >> 1 if fw if oopc @@ -125,16 +125,16 @@ end # inplace transform of y, no vector allocation # tmp: size at least n>>2 # tmpvec: size at least n -function _dwt!(y::Matrix{T}, scheme::GLS, L::Integer, fw::Bool, tmp::Vector{T} = Vector{T}(undef, reqtmplength(y)), tmpvec::Vector{T} = Vector{T}(undef, size(y,1))) where T<:Number +function _dwt!(y::Matrix{T}, scheme::GLS, L::Integer, fw::Bool, tmp::Vector{T}=Vector{T}(undef, reqtmplength(y)), tmpvec::Vector{T}=Vector{T}(undef, size(y, 1))) where {T<:Number} - n = size(y,1) + n = size(y, 1) iscube(y) || throw(ArgumentError("array must be square/cube")) 0 <= L || throw(ArgumentError("L must be positive")) sufficientpoweroftwo(y, L) || throw(ArgumentError("size must have a sufficient power of 2 factor")) - (length(tmp) >= n>>2) || + (length(tmp) >= n >> 2) || throw(ArgumentError("length of tmp incorrect")) (length(tmpvec) >= n) || throw(ArgumentError("length of tmpvec incorrect")) @@ -143,14 +143,14 @@ function _dwt!(y::Matrix{T}, scheme::GLS, L::Integer, fw::Bool, tmp::Vector{T} = return y end row_stride = n - plane_stride = n*n + plane_stride = n * n if fw lrange = 1:L nsub = n else lrange = L:-1:1 - nsub = div(n,2^(L-1)) + nsub = div(n, 2^(L - 1)) end stepseq, norm1, norm2 = makescheme(T, scheme, fw) @@ -162,14 +162,14 @@ function _dwt!(y::Matrix{T}, scheme::GLS, L::Integer, fw::Bool, tmp::Vector{T} = for i in 1:nsub xi = row_idx(i, n) unsafe_dwt1level!(y, xi, row_stride, true, tmpsub, scheme, fw, - stepseq, norm1, norm2, tmp) + stepseq, norm1, norm2, tmp) end # columns for i in 1:nsub xi = col_idx(i, n) ya = unsafe_vectorslice(y, xi, nsub) unsafe_dwt1level!(ya, 1, 1, false, tmpsub, scheme, fw, - stepseq, norm1, norm2, tmp) + stepseq, norm1, norm2, tmp) end else # columns @@ -177,17 +177,17 @@ function _dwt!(y::Matrix{T}, scheme::GLS, L::Integer, fw::Bool, tmp::Vector{T} = xi = col_idx(i, n) ya = unsafe_vectorslice(y, xi, nsub) unsafe_dwt1level!(ya, 1, 1, false, tmpsub, scheme, fw, - stepseq, norm1, norm2, tmp) + stepseq, norm1, norm2, tmp) end # rows for i in 1:nsub xi = row_idx(i, n) unsafe_dwt1level!(y, xi, row_stride, true, tmpsub, scheme, fw, - stepseq, norm1, norm2, tmp) + stepseq, norm1, norm2, tmp) end end - nsub = (fw ? nsub>>1 : nsub<<1) + nsub = (fw ? nsub >> 1 : nsub << 1) end return y @@ -197,16 +197,16 @@ end # inplace transform of y, no vector allocation # tmp: size at least n>>2 # tmpvec: size at least n -function _dwt!(y::Array{T,3}, scheme::GLS, L::Integer, fw::Bool, tmp::Vector{T} = Vector{T}(undef, reqtmplength(y)), tmpvec::Vector{T} = Vector{T}(undef, size(y,1))) where T<:Number +function _dwt!(y::Array{T,3}, scheme::GLS, L::Integer, fw::Bool, tmp::Vector{T}=Vector{T}(undef, reqtmplength(y)), tmpvec::Vector{T}=Vector{T}(undef, size(y, 1))) where {T<:Number} - n = size(y,1) + n = size(y, 1) iscube(y) || throw(ArgumentError("array must be square/cube")) 0 <= L || throw(ArgumentError("L must be positive")) sufficientpoweroftwo(y, L) || throw(ArgumentError("size must have a sufficient power of 2 factor")) - (length(tmp) >= n>>2) || + (length(tmp) >= n >> 2) || throw(ArgumentError("length of tmp incorrect")) (length(tmpvec) >= n) || throw(ArgumentError("length of tmpvec incorrect")) @@ -215,14 +215,14 @@ function _dwt!(y::Array{T,3}, scheme::GLS, L::Integer, fw::Bool, tmp::Vector{T} return y end row_stride = n - plane_stride = n*n + plane_stride = n * n if fw lrange = 1:L nsub = n else lrange = L:-1:1 - nsub = div(n,2^(L-1)) + nsub = div(n, 2^(L - 1)) end stepseq, norm1, norm2 = makescheme(T, scheme, fw) @@ -234,20 +234,20 @@ function _dwt!(y::Array{T,3}, scheme::GLS, L::Integer, fw::Bool, tmp::Vector{T} for i in 1:nsub, j in 1:nsub xi = plane_idx(i, j, n) unsafe_dwt1level!(y, xi, plane_stride, true, tmpsub, scheme, fw, - stepseq, norm1, norm2, tmp) + stepseq, norm1, norm2, tmp) end # rows for i in 1:nsub, j in 1:nsub xi = row_idx(i, j, n) unsafe_dwt1level!(y, xi, row_stride, true, tmpsub, scheme, fw, - stepseq, norm1, norm2, tmp) + stepseq, norm1, norm2, tmp) end # columns for i in 1:nsub, j in 1:nsub xi = col_idx(i, j, n) ya = unsafe_vectorslice(y, xi, nsub) unsafe_dwt1level!(ya, 1, 1, false, tmpsub, scheme, fw, - stepseq, norm1, norm2, tmp) + stepseq, norm1, norm2, tmp) end else # columns @@ -255,23 +255,23 @@ function _dwt!(y::Array{T,3}, scheme::GLS, L::Integer, fw::Bool, tmp::Vector{T} xi = col_idx(i, j, n) ya = unsafe_vectorslice(y, xi, nsub) unsafe_dwt1level!(ya, 1, 1, false, tmpsub, scheme, fw, - stepseq, norm1, norm2, tmp) + stepseq, norm1, norm2, tmp) end # rows for i in 1:nsub, j in 1:nsub xi = row_idx(i, j, n) unsafe_dwt1level!(y, xi, row_stride, true, tmpsub, scheme, fw, - stepseq, norm1, norm2, tmp) + stepseq, norm1, norm2, tmp) end # planes for i in 1:nsub, j in 1:nsub xi = plane_idx(i, j, n) unsafe_dwt1level!(y, xi, plane_stride, true, tmpsub, scheme, fw, - stepseq, norm1, norm2, tmp) + stepseq, norm1, norm2, tmp) end end - nsub = (fw ? nsub>>1 : nsub<<1) + nsub = (fw ? nsub >> 1 : nsub << 1) end return y @@ -280,12 +280,12 @@ end # WPT # 1-D -function _wpt!(y::AbstractVector{T}, scheme::GLS, tree::BitVector, fw::Bool, tmp::Vector{T} = Vector{T}(undef, reqtmplength(y))) where T<:Number +function _wpt!(y::AbstractVector{T}, scheme::GLS, tree::BitVector, fw::Bool, tmp::Vector{T}=Vector{T}(undef, reqtmplength(y))) where {T<:Number} n = length(y) isvalidtree(y, tree) || throw(ArgumentError("invalid tree")) - length(tmp) >= n>>2 || + length(tmp) >= n >> 2 || throw(ArgumentError("length of tmp incorrect")) if !tree[1] @@ -299,10 +299,10 @@ function _wpt!(y::AbstractVector{T}, scheme::GLS, tree::BitVector, fw::Bool, tmp while L > 0 ix = 1 k = 1 - fw && (Lfw = Lmax-L) - !fw && (Lfw = L-1) + fw && (Lfw = Lmax - L) + !fw && (Lfw = L - 1) nj = detailn(n, Lfw) - treeind = 2^(Lfw)-1 + treeind = 2^(Lfw) - 1 while ix <= n if tree[treeind+k] @@ -320,7 +320,7 @@ end -function normalize!(x::AbstractVector{T}, half::Int, ns::Int, n1::T, n2::T) where T<:Number +function normalize!(x::AbstractVector{T}, half::Int, ns::Int, n1::T, n2::T) where {T<:Number} for i = 1:half @inbounds x[i] *= n1 end @@ -330,30 +330,30 @@ function normalize!(x::AbstractVector{T}, half::Int, ns::Int, n1::T, n2::T) wher return x end # out of place normalize from x to y -function normalize!(y::AbstractVector{T}, x::AbstractVector{T}, half::Int, ns::Int, n1::T, n2::T) where T<:Number +function normalize!(y::AbstractVector{T}, x::AbstractVector{T}, half::Int, ns::Int, n1::T, n2::T) where {T<:Number} for i = 1:half - @inbounds y[i] = n1*x[i] + @inbounds y[i] = n1 * x[i] end for i = half+1:ns - @inbounds y[i] = n2*x[i] + @inbounds y[i] = n2 * x[i] end return y end -function normalize!(y::AbstractArray{T}, iy::Int, incy::Int, x::AbstractVector{T}, half::Int, ns::Int, n1::T, n2::T) where T<:Number +function normalize!(y::AbstractArray{T}, iy::Int, incy::Int, x::AbstractVector{T}, half::Int, ns::Int, n1::T, n2::T) where {T<:Number} for i = 1:half - @inbounds y[iy + (i-1)*incy] = n1*x[i] + @inbounds y[iy+(i-1)*incy] = n1 * x[i] end for i = half+1:ns - @inbounds y[iy + (i-1)*incy] = n2*x[i] + @inbounds y[iy+(i-1)*incy] = n2 * x[i] end return y end -function normalize!(y::AbstractVector{T}, x::AbstractArray{T}, ix::Int, incx::Int, half::Int, ns::Int, n1::T, n2::T) where T<:Number +function normalize!(y::AbstractVector{T}, x::AbstractArray{T}, ix::Int, incx::Int, half::Int, ns::Int, n1::T, n2::T) where {T<:Number} for i = 1:half - @inbounds y[i] = n1*x[ix + (i-1)*incx] + @inbounds y[i] = n1 * x[ix+(i-1)*incx] end for i = half+1:ns - @inbounds y[i] = n2*x[ix + (i-1)*incx] + @inbounds y[i] = n2 * x[ix+(i-1)*incx] end return y end @@ -364,20 +364,20 @@ end # For predict: writes to range 1:half, reads from 1:2*half # For update : writes to range half+1:2*half, reads from 1:2*half for step_type in (WT.PredictStep, WT.UpdateStep) -@eval begin -function lift!(x::AbstractVector{T}, half::Int, - param::WT.LSStepParam{T}, steptype::$step_type) where T<:Number - lhsr, irange, rhsr, rhsis = getliftranges(half, length(param), param.shift, steptype) - coefs = param.coef - # left boundary - lift_perboundary!(x, half, coefs, lhsr, rhsis, steptype) - # main loop - lift_inbounds!(x, coefs, irange, rhsis) - # right boundary - lift_perboundary!(x, half, coefs, rhsr, rhsis, steptype) - return x -end -end # eval begin + @eval begin + function lift!(x::AbstractVector{T}, half::Int, + param::WT.LSStepParam{T}, steptype::$step_type) where {T<:Number} + lhsr, irange, rhsr, rhsis = getliftranges(half, length(param), param.shift, steptype) + coefs = param.coef + # left boundary + lift_perboundary!(x, half, coefs, lhsr, rhsis, steptype) + # main loop + lift_inbounds!(x, coefs, irange, rhsis) + # right boundary + lift_perboundary!(x, half, coefs, rhsr, rhsis, steptype) + return x + end + end # eval begin end # for function irlimits(half::Int, nc::Int, shift::Int) @@ -385,35 +385,35 @@ function irlimits(half::Int, nc::Int, shift::Int) # 1 <= i <= half # 1 <= i+1-1-shift <= half # 1 <= i+nc-1-shift <= half - return (max(shift+1, 1-nc+shift), - min(half+1+shift-nc, half+shift)) + return (max(shift + 1, 1 - nc + shift), + min(half + 1 + shift - nc, half + shift)) end function getliftranges(half::Int, nc::Int, shift::Int, steptype::WT.UpdateStep) # define index shift rhsis - rhsis = -shift-half + rhsis = -shift - half irmin, irmax = irlimits(half, nc, shift) if irmin > half || irmax < 1 irange = 1:0 # empty else irmin = max(irmin, 1) irmax = min(irmax, half) - irange = irmin + half:irmax + half + irange = irmin+half:irmax+half end # periodic boundary - if length(irange)==0 - lhsr = 1 + half:half + half - rhsr = 1 + half:0 + half + if length(irange) == 0 + lhsr = 1+half:half+half + rhsr = 1+half:0+half else - lhsr = 1 + half:irmin - 1 + half - rhsr = irmax + 1 + half:half + half + lhsr = 1+half:irmin-1+half + rhsr = irmax+1+half:half+half end return (lhsr, irange, rhsr, rhsis) end function getliftranges(half::Int, nc::Int, shift::Int, steptype::WT.PredictStep) # define index shift rhsis - rhsis = -shift+half + rhsis = -shift + half irmin, irmax = irlimits(half, nc, shift) if irmin > half || irmax < 1 irange = 1:0 # empty @@ -423,59 +423,59 @@ function getliftranges(half::Int, nc::Int, shift::Int, steptype::WT.PredictStep) irange = irmin:irmax end # periodic boundary - if length(irange)==0 + if length(irange) == 0 lhsr = 1:half rhsr = 1:0 else - lhsr = 1:irmin - 1 - rhsr = irmax + 1:half + lhsr = 1:irmin-1 + rhsr = irmax+1:half end return (lhsr, irange, rhsr, rhsis) end # periodic boundary -for (step_type, puxind) in ((WT.PredictStep, :(mod1(i+k-1+rhsis-half,half)+half)), - (WT.UpdateStep, :(mod1(i+k-1+rhsis,half))) ) -@eval begin -function lift_perboundary!(x::AbstractVector{T}, half::Int, - c::Vector{T}, irange::AbstractRange, rhsis::Int, ::$step_type) where T<:Number - nc = length(c) - for i in irange - for k in 1:nc - @inbounds x[i] += c[k]*x[$puxind] +for (step_type, puxind) in ((WT.PredictStep, :(mod1(i + k - 1 + rhsis - half, half) + half)), + (WT.UpdateStep, :(mod1(i + k - 1 + rhsis, half)))) + @eval begin + function lift_perboundary!(x::AbstractVector{T}, half::Int, + c::Vector{T}, irange::AbstractRange, rhsis::Int, ::$step_type) where {T<:Number} + nc = length(c) + for i in irange + for k in 1:nc + @inbounds x[i] += c[k] * x[$puxind] + end + end + return x end - end - return x -end -end # eval begin + end # eval begin end # for # main lift loop -function lift_inbounds!(x::AbstractVector{T}, c::Vector{T}, irange::AbstractRange, rhsis::Int) where T<:Number +function lift_inbounds!(x::AbstractVector{T}, c::Vector{T}, irange::AbstractRange, rhsis::Int) where {T<:Number} nc = length(c) if nc == 1 # hard code the most common cases (1, 2, 3) for speed c1 = c[1] for i in irange - @inbounds x[i] += c1*x[i+rhsis] + @inbounds x[i] += c1 * x[i+rhsis] end elseif nc == 2 - c1,c2 = c[1],c[2] - rhsisp1 = rhsis+1 + c1, c2 = c[1], c[2] + rhsisp1 = rhsis + 1 for i in irange - @inbounds x[i] += c1*x[i+rhsis] + c2*x[i+rhsisp1] + @inbounds x[i] += c1 * x[i+rhsis] + c2 * x[i+rhsisp1] end elseif nc == 3 - c1,c2,c3 = c[1],c[2],c[3] - rhsisp1 = rhsis+1 - rhsisp2 = rhsis+2 + c1, c2, c3 = c[1], c[2], c[3] + rhsisp1 = rhsis + 1 + rhsisp2 = rhsis + 2 for i = irange - @inbounds x[i] += c1*x[i+rhsis] + c2*x[i+rhsisp1] + c3*x[i+rhsisp2] + @inbounds x[i] += c1 * x[i+rhsis] + c2 * x[i+rhsisp1] + c3 * x[i+rhsisp2] end else for i in irange for k in 0:nc-1 - @inbounds x[i] += c[k]*x[i+k+rhsis] + @inbounds x[i] += c[k] * x[i+k+rhsis] end end end diff --git a/src/mod/transforms_maximal_overlap.jl b/src/Transforms/transforms_maximal_overlap.jl similarity index 90% rename from src/mod/transforms_maximal_overlap.jl rename to src/Transforms/transforms_maximal_overlap.jl index 2dc3643..9747ebb 100644 --- a/src/mod/transforms_maximal_overlap.jl +++ b/src/Transforms/transforms_maximal_overlap.jl @@ -8,7 +8,7 @@ scaling filters. Returns a tuple `(v, w)` of the scaling and detail coefficients at level `j+1`. """ function modwt_step(v::AbstractVector{T}, j::Integer, h::Array{S,1}, - g::Array{S,1}) where {T <: Number, S <: Number} + g::Array{S,1}) where {T<:Number,S<:Number} N = length(v) L = length(h) v1 = zeros(T, N) @@ -18,7 +18,7 @@ function modwt_step(v::AbstractVector{T}, j::Integer, h::Array{S,1}, w1[t] = h[1] * v[k] v1[t] = g[1] * v[k] for n in 2:L - k -= 2^(j-1) + k -= 2^(j - 1) if k <= 0 k = mod1(k, N) end @@ -44,7 +44,7 @@ coefficients for level j in column j. The scaling coefficients are in the last (L+1th) column. """ function modwt(x::AbstractVector{T}, wt::OrthoFilter, - L::Integer=maxmodwttransformlevels(x)) where T <: Number + L::Integer=maxmodwttransformlevels(x)) where {T<:Number} L <= maxmodwttransformlevels(x) || throw(ArgumentError("Too many transform levels (length(x) < 2^L)")) L >= 1 || throw(ArgumentError("L must be >= 1")) @@ -68,7 +68,7 @@ and returns a vector of the `j-1`th level scaling coefficients. The vectors `h` and `g` are the MODWT detail and scaling filters. """ function imodwt_step(v::AbstractVector{T}, w::AbstractVector{T}, j::Integer, - h::Array{S,1}, g::Array{S,1}) where {T<:Number, S<:Number} + h::Array{S,1}, g::Array{S,1}) where {T<:Number,S<:Number} length(v) == length(w) || throw(DimensionMismatch("Input array sizes must match")) length(h) == length(g) || @@ -80,7 +80,7 @@ function imodwt_step(v::AbstractVector{T}, w::AbstractVector{T}, j::Integer, k = t v0[t] = h[1] * w[k] + g[1] * v[k] for n in 2:L - k += 2^(j-1) + k += 2^(j - 1) if k > N k = mod1(k, N) end @@ -94,7 +94,7 @@ end Perform an inverse maximal overlap discrete wavelet transform (MODWT) of `xw`, the inverse of `modwt(x, wt, L)`. """ -function imodwt(xw::Array{T, 2}, wt::OrthoFilter) where T <: Number +function imodwt(xw::Array{T,2}, wt::OrthoFilter) where {T<:Number} g, h = WT.makereverseqmfpair(wt) g /= sqrt(2) h /= sqrt(2) diff --git a/src/Util/Util.jl b/src/Util/Util.jl new file mode 100644 index 0000000..881435c --- /dev/null +++ b/src/Util/Util.jl @@ -0,0 +1,44 @@ +module Util + +using LinearAlgebra: rmul!, norm +using DSP: conv + + + +include("non_dyadic.jl") +include("dyadic.jl") +include("util.jl") + + + +export + # dyadic + dyadicdetailindex, + dyadicdetailrange, + dyadicscalingrange, + dyadicdetailn, + ndyadicscales, + maxdyadiclevel, + tl2dyadiclevel, + dyadiclevel2tl, + # non-dyadic + detailindex, + detailrange, + detailn, + maxtransformlevels, + maxmodwttransformlevels, + # other util functions + mirror, + upsample, + downsample, + iscube, + isdyadic, + sufficientpoweroftwo, + wcount, + stridedcopy!, + isvalidtree, + maketree, + makewavelet, + testfunction + +end diff --git a/src/Util/dyadic.jl b/src/Util/dyadic.jl new file mode 100644 index 0000000..a04d93b --- /dev/null +++ b/src/Util/dyadic.jl @@ -0,0 +1,20 @@ +# Dyadic +# detail coef at level j location i (i starting at 1) -> vector index +dyadicdetailindex(j::Integer, i::Integer) = 2^j + i +# the range of detail coefs at level j +dyadicdetailrange(j::Integer) = (2^j+1):(2^(j+1)) +# the range of scaling coefs at level j +dyadicscalingrange(j::Integer) = 1:2^j +# number of detail coefs at level j +dyadicdetailn(j::Integer) = 2^j +# number of scales of dyadic length signal (n=2^J) +ndyadicscales(n::Integer) = round(Int, log2(n)) +ndyadicscales(x::AbstractArray) = ndyadicscales(size(x, 1)) +# the largest detail level +maxdyadiclevel(n::Integer) = ndyadicscales(n) - 1 +maxdyadiclevel(x::AbstractArray) = maxdyadiclevel(size(x, 1)) +# convert number of transformed levels L to minimum level j +tl2dyadiclevel(n::Integer, L::Integer) = ndyadicscales(n) - L +tl2dyadiclevel(x::AbstractVector, L::Integer) = tl2dyadiclevel(length(x), L) +# convert maximum level j to number of transformed levels L +dyadiclevel2tl(arg...) = tl2dyadiclevel(arg...) diff --git a/src/Util/non_dyadic.jl b/src/Util/non_dyadic.jl new file mode 100644 index 0000000..d6babdf --- /dev/null +++ b/src/Util/non_dyadic.jl @@ -0,0 +1,27 @@ +using ..Util + +# WAVELET INDEXING AND SIZES + +# Non-dyadic +# detail coef at level l location i (i starting at 1) -> vector index +detailindex(arraysize::Integer, l::Integer, i::Integer) = round(Int, arraysize / 2^l + i) +detailindex(x::AbstractArray, l::Integer, i::Integer) = detailindex(size(x, 1), l, i) +# the range of detail coefs at level l +detailrange(arraysize::Integer, l::Integer) = round(Int, (arraysize / 2^l + 1)):round(Int, arraysize / 2^(l - 1)) +detailrange(x::AbstractArray, l::Integer) = detailrange(size(x, 1), l) +# number of detail coefs at level l +detailn(arraysize::Integer, l::Integer) = round(Int, arraysize / 2^l) +detailn(x::AbstractArray, l::Integer) = detailn(size(x, 1), l) +# max levels to transform +maxtransformlevels(x::AbstractArray) = minimum(maxtransformlevels.(size(x))) +function maxtransformlevels(arraysize::Integer) + arraysize > 1 || return 0 + tl = 0 + while sufficientpoweroftwo(arraysize, tl) + tl += 1 + end + return tl - 1 +end + +maxmodwttransformlevels(x::AbstractArray) = floor(Int, log2(length(x))) +maxmodwttransformlevels(arraysize::Integer) = floor(Int, log2(arraysize)) diff --git a/src/mod/Util.jl b/src/Util/util.jl similarity index 53% rename from src/mod/Util.jl rename to src/Util/util.jl index 105de00..1e0c09a 100644 --- a/src/mod/Util.jl +++ b/src/Util/util.jl @@ -1,86 +1,4 @@ -module Util -export dyadicdetailindex, - dyadicdetailrange, - dyadicscalingrange, - dyadicdetailn, - ndyadicscales, - maxdyadiclevel, - tl2dyadiclevel, - dyadiclevel2tl, - # non-dyadic - detailindex, - detailrange, - detailn, - maxtransformlevels, - maxmodwttransformlevels, - # - mirror, - upsample, - downsample, - iscube, - isdyadic, - sufficientpoweroftwo, - wcount, - stridedcopy!, - isvalidtree, - maketree, - makewavelet, - testfunction - - -if VERSION >= v"0.7.0-DEV.986" - using DSP: conv -end -using LinearAlgebra - - -# WAVELET INDEXING AND SIZES - -# Non-dyadic -# detail coef at level l location i (i starting at 1) -> vector index -detailindex(arraysize::Integer, l::Integer, i::Integer) = round(Int, arraysize/2^l+i) -detailindex(x::AbstractArray, l::Integer, i::Integer) = detailindex(size(x,1), l ,i) -# the range of detail coefs at level l -detailrange(arraysize::Integer, l::Integer) = round(Int, (arraysize/2^l+1)) : round(Int, arraysize/2^(l-1)) -detailrange(x::AbstractArray, l::Integer) = detailrange(size(x,1), l) -# number of detail coefs at level l -detailn(arraysize::Integer, l::Integer) = round(Int, arraysize/2^l) -detailn(x::AbstractArray, l::Integer) = detailn(size(x,1), l) -# max levels to transform -maxtransformlevels(x::AbstractArray) = minimum(maxtransformlevels.(size(x))) -function maxtransformlevels(arraysize::Integer) - arraysize > 1 || return 0 - tl = 0 - while sufficientpoweroftwo(arraysize, tl) - tl += 1 - end - return tl - 1 -end - -maxmodwttransformlevels(x::AbstractArray) = floor(Int, log2(length(x))) -maxmodwttransformlevels(arraysize::Integer) = floor(Int, log2(arraysize)) - -# Dyadic -# detail coef at level j location i (i starting at 1) -> vector index -dyadicdetailindex(j::Integer,i::Integer) = 2^j+i -# the range of detail coefs at level j -dyadicdetailrange(j::Integer) = (2^j+1):(2^(j+1)) -# the range of scaling coefs at level j -dyadicscalingrange(j::Integer) = 1:2^j -# number of detail coefs at level j -dyadicdetailn(j::Integer) = 2^j -# number of scales of dyadic length signal (n=2^J) -ndyadicscales(n::Integer) = round(Int, log2(n)) -ndyadicscales(x::AbstractArray) = ndyadicscales(size(x,1)) -# the largest detail level -maxdyadiclevel(n::Integer) = ndyadicscales(n)-1 -maxdyadiclevel(x::AbstractArray) = maxdyadiclevel(size(x,1)) -# convert number of transformed levels L to minimum level j -tl2dyadiclevel(n::Integer, L::Integer) = ndyadicscales(n)-L -tl2dyadiclevel(x::AbstractVector, L::Integer) = tl2dyadiclevel(length(x), L) -# convert maximum level j to number of transformed levels L -dyadiclevel2tl(arg...) = tl2dyadiclevel(arg...) - +using ..Util # UTILITY FUNCTIONS @@ -88,14 +6,14 @@ dyadiclevel2tl(arg...) = tl2dyadiclevel(arg...) # are all dimensions equal length? function iscube(x::AbstractArray) for i = 1:ndims(x) - size(x,1) != size(x,i) && return false + size(x, 1) != size(x, i) && return false end return true end # are all dimensions dyadic length? function isdyadic(x::AbstractArray) for i = 1:ndims(x) - isdyadic(size(x,i)) || return false + isdyadic(size(x, i)) || return false end return true end @@ -105,36 +23,36 @@ isdyadic(n::Integer) = (n == 2^(ndyadicscales(n))) # must have 2^L as a factor. function sufficientpoweroftwo(x::AbstractArray, L::Integer) for i = 1:ndims(x) - sufficientpoweroftwo(size(x,i), L) || return false + sufficientpoweroftwo(size(x, i), L) || return false end return true end -sufficientpoweroftwo(n::Integer, L::Integer) = (n%(2^L) == 0) +sufficientpoweroftwo(n::Integer, L::Integer) = (n % (2^L) == 0) # mirror of filter -mirror(f::AbstractVector{<:Number}) = f .* (-1).^(0:length(f)-1) +mirror(f::AbstractVector{<:Number}) = f .* (-1) .^ (0:length(f)-1) # upsample function upsample(x::AbstractVector, sw::Int=0) - @assert sw==0 || sw==1 + @assert sw == 0 || sw == 1 n = length(x) - y = zeros(eltype(x), n<<1) + y = zeros(eltype(x), n << 1) sw -= 1 for i = 1:n - @inbounds y[i<<1 + sw] = x[i] + @inbounds y[i<<1+sw] = x[i] end return y end # downsample function downsample(x::AbstractVector, sw::Int=0) - @assert sw==0 || sw==1 + @assert sw == 0 || sw == 1 n = length(x) - @assert n%2 == 0 - y = zeros(eltype(x), n>>1) + @assert n % 2 == 0 + y = zeros(eltype(x), n >> 1) sw -= 1 for i in eachindex(y) - @inbounds y[i] = x[i<<1 + sw] + @inbounds y[i] = x[i<<1+sw] end return y end @@ -168,9 +86,9 @@ end function circshift!(a::AbstractVector, shift::Integer) atype = typeof(a) s = length(a) - shift = mod(shift,s) + shift = mod(shift, s) shift == 0 && return a - shift = ifelse(s>>1 < shift, shift-s, shift) # shift a smaller distance if possible + shift = ifelse(s >> 1 < shift, shift - s, shift) # shift a smaller distance if possible if shift < 0 tmp = a[1:-shift] for i = 1:s+shift @@ -191,14 +109,14 @@ function circshift!(b::AbstractVector, a::AbstractVector, shift::Integer) @assert length(a) == length(b) atype = typeof(a) s = length(a) - shift = mod(shift,s) - shift == 0 && return copyto!(b,a) - shift = ifelse(s>>1 < shift, shift-s, shift) # shift a smaller distance if possible + shift = mod(shift, s) + shift == 0 && return copyto!(b, a) + shift = ifelse(s >> 1 < shift, shift - s, shift) # shift a smaller distance if possible if shift < 0 for i = 1:s+shift @inbounds b[i] = a[i-shift] end - sh = s+shift + sh = s + shift for i = s+shift+1:s @inbounds b[i] = a[i-sh] end @@ -206,7 +124,7 @@ function circshift!(b::AbstractVector, a::AbstractVector, shift::Integer) for i = s:-1:1+shift @inbounds b[i] = a[i-shift] end - sh = -s+shift + sh = -s + shift for i = shift:-1:1 @inbounds b[i] = a[i-sh] end @@ -215,148 +133,148 @@ function circshift!(b::AbstractVector, a::AbstractVector, shift::Integer) end # put odd elements into first half, even into second half -function split!(a::AbstractVector{T}) where T<:Number +function split!(a::AbstractVector{T}) where {T<:Number} n = length(a) - nt = n>>2 + (n>>1)%2 + nt = n >> 2 + (n >> 1) % 2 tmp = Vector{T}(undef, nt) split!(a, n, tmp) return a end # split only the range 1:n -function split!(a::AbstractVector{T}, n::Integer, tmp::Vector{T}) where T<:Number +function split!(a::AbstractVector{T}, n::Integer, tmp::Vector{T}) where {T<:Number} @assert n <= length(a) - @assert n%2 == 0 + @assert n % 2 == 0 n == 2 && return a - nt = n>>2 + (n>>1)%2 + nt = n >> 2 + (n >> 1) % 2 @assert nt <= length(tmp) - for i=1:nt # store evens + for i = 1:nt # store evens @inbounds tmp[i] = a[i<<1] end - for i=1:n>>1 # odds to first part - @inbounds a[i] = a[(i-1)<<1 + 1] + for i = 1:n>>1 # odds to first part + @inbounds a[i] = a[(i-1)<<1+1] end - for i=0:nt - 1 # evens to end - @inbounds a[n-i] = a[n - i<<1] + for i = 0:nt-1 # evens to end + @inbounds a[n-i] = a[n-i<<1] end - copyto!(a,n>>1+1,tmp,1,nt) + copyto!(a, n >> 1 + 1, tmp, 1, nt) return a end # out of place split from a to b, only the range 1:n -function split!(b::AbstractVector{T}, a::AbstractVector{T}, n::Integer) where T<:Number +function split!(b::AbstractVector{T}, a::AbstractVector{T}, n::Integer) where {T<:Number} @assert n <= length(a) && n <= length(b) - @assert n%2 == 0 + @assert n % 2 == 0 if n == 2 b[1] = a[1] b[2] = a[2] return b end - h = n>>1 - for i=1:h # odds to b - @inbounds b[i] = a[(i-1)<<1 + 1] + h = n >> 1 + for i = 1:h # odds to b + @inbounds b[i] = a[(i-1)<<1+1] end - for i=h+1:n # evens to b - @inbounds b[i] = a[(i - h)<<1] + for i = h+1:n # evens to b + @inbounds b[i] = a[(i-h)<<1] end return b end # out of place split from a to b, only the range a[ia:inca:ia+(n-1)*inca] to b[1:n] -function split!(b::AbstractVector{T}, a::AbstractArray{T}, ia::Integer, inca::Integer, n::Integer) where T<:Number - @assert ia+(n-1)*inca <= length(a) && n <= length(b) - @assert n%2 == 0 +function split!(b::AbstractVector{T}, a::AbstractArray{T}, ia::Integer, inca::Integer, n::Integer) where {T<:Number} + @assert ia + (n - 1) * inca <= length(a) && n <= length(b) + @assert n % 2 == 0 if n == 2 b[1] = a[ia] - b[2] = a[ia + inca] + b[2] = a[ia+inca] return b end - h = n>>1 - inca2 = inca<<1 - for i=1:h # odds to b - @inbounds b[i] = a[ia + (i-1)*inca2] + h = n >> 1 + inca2 = inca << 1 + for i = 1:h # odds to b + @inbounds b[i] = a[ia+(i-1)*inca2] end iainca = ia + inca hp1 = h + 1 - for i=h+1:n # evens to b - @inbounds b[i] = a[iainca + (i - hp1)*inca2] + for i = h+1:n # evens to b + @inbounds b[i] = a[iainca+(i-hp1)*inca2] end return b end # inverse the operation of split! -function merge!(a::AbstractVector{T}) where T<:Number +function merge!(a::AbstractVector{T}) where {T<:Number} n = length(a) - nt = n>>2 + (n>>1)%2 + nt = n >> 2 + (n >> 1) % 2 tmp = Vector{T}(undef, nt) merge!(a, n, tmp) return a end # merge only the range 1:n -function merge!(a::AbstractVector{T}, n::Integer, tmp::Vector{T}) where T<:Number +function merge!(a::AbstractVector{T}, n::Integer, tmp::Vector{T}) where {T<:Number} @assert n <= length(a) - @assert n%2 == 0 + @assert n % 2 == 0 n == 2 && return a - nt = n>>2 + (n>>1)%2 + nt = n >> 2 + (n >> 1) % 2 @assert nt <= length(tmp) - copyto!(tmp,1,a,n>>1+1,nt) - for i=nt-1:-1:0 # evens from end - @inbounds a[n - 2*i] = a[n-i] + copyto!(tmp, 1, a, n >> 1 + 1, nt) + for i = nt-1:-1:0 # evens from end + @inbounds a[n-2*i] = a[n-i] end - for i=n>>1:-1:1 # odds from first part - @inbounds a[(i-1)<<1 + 1] = a[i] + for i = n>>1:-1:1 # odds from first part + @inbounds a[(i-1)<<1+1] = a[i] end - for i=nt:-1:1 # retrieve evens + for i = nt:-1:1 # retrieve evens @inbounds a[i<<1] = tmp[i] end return a end # out of place merge from a to b, only the range 1:n -function merge!(b::AbstractVector{T}, a::AbstractVector{T}, n::Integer) where T<:Number +function merge!(b::AbstractVector{T}, a::AbstractVector{T}, n::Integer) where {T<:Number} @assert n <= length(a) && n <= length(b) - @assert n%2 == 0 + @assert n % 2 == 0 if n == 2 b[1] = a[1] b[2] = a[2] return b end - h = n>>1 - for i=1:h # odds to b - @inbounds b[(i-1)<<1 + 1] = a[i] + h = n >> 1 + for i = 1:h # odds to b + @inbounds b[(i-1)<<1+1] = a[i] end - for i=h+1:n # evens to b - @inbounds b[(i - h)<<1] = a[i] + for i = h+1:n # evens to b + @inbounds b[(i-h)<<1] = a[i] end return b end # out of place merge from a to b, only the range a[1:n] to b[ib:incb:ib+(n-1)*incb] -function merge!(b::AbstractArray{T}, ib::Integer, incb::Integer, a::AbstractVector{T}, n::Integer) where T<:Number - @assert n <= length(a) && ib+(n-1)*incb <= length(b) - @assert n%2 == 0 +function merge!(b::AbstractArray{T}, ib::Integer, incb::Integer, a::AbstractVector{T}, n::Integer) where {T<:Number} + @assert n <= length(a) && ib + (n - 1) * incb <= length(b) + @assert n % 2 == 0 if n == 2 b[ib] = a[1] - b[ib + incb] = a[2] + b[ib+incb] = a[2] return b end - h = n>>1 - incb2 = incb<<1 - for i=1:h # odds to b - @inbounds b[ib + (i-1)*incb2] = a[i] + h = n >> 1 + incb2 = incb << 1 + for i = 1:h # odds to b + @inbounds b[ib+(i-1)*incb2] = a[i] end ibincb = ib + incb hp1 = h + 1 - for i=h+1:n # evens to b - @inbounds b[ibincb + (i - hp1)*incb2] = a[i] + for i = h+1:n # evens to b + @inbounds b[ibincb+(i-hp1)*incb2] = a[i] end return b @@ -364,18 +282,18 @@ end function stridedcopy!(b::AbstractVector{<:Number}, a::AbstractArray{<:Number}, ia::Integer, inca::Integer, n::Integer) - @assert ia+(n-1)*inca <= length(a) && n <= length(b) + @assert ia + (n - 1) * inca <= length(a) && n <= length(b) @inbounds for i = 1:n - b[i] = a[ia + (i-1)*inca] + b[i] = a[ia+(i-1)*inca] end return b end function stridedcopy!(b::AbstractArray{<:Number}, ib::Integer, incb::Integer, a::AbstractVector{<:Number}, n::Integer) - @assert ib+(n-1)*incb <= length(b) && n <= length(a) + @assert ib + (n - 1) * incb <= length(b) && n <= length(a) @inbounds for i = 1:n - b[ib + (i-1)*incb] = a[i] + b[ib+(i-1)*incb] = a[i] end return b end @@ -386,8 +304,8 @@ end function isvalidtree(x::AbstractVector, b::BitVector) ns = maxtransformlevels(x) nb = length(b) - nb == 2^(ns)-1 || return false - @assert (2^(ns-1)-1)<<1+1 <= nb + nb == 2^(ns) - 1 || return false + @assert (2^(ns - 1) - 1) << 1 + 1 <= nb for i in 1:2^(ns-1)-1 @inbounds if !b[i] && (b[i<<1] || b[i<<1+1]) @@ -406,9 +324,9 @@ s=:dwt, nodes corresponding to a dwt for first L levels equal 1, others 0 maketree(x::Vector, s::Symbol=:full) = maketree(length(x), maxtransformlevels(x), s) function maketree(n::Int, L::Int, s::Symbol=:full) ns = maxtransformlevels(n) - nb = 2^(ns)-1 + nb = 2^(ns) - 1 @assert 0 <= L <= ns - @assert (2^(ns-1)-1)<<1+1 <= nb + @assert (2^(ns - 1) - 1) << 1 + 1 <= nb b = BitArray(undef, nb) fill!(b, false) @@ -435,19 +353,19 @@ iterated with a cascade algorithm with N steps """ makewavelet(h, arg...) = makewavelet(h.qmf, arg...) function makewavelet(h::AbstractVector, N::Integer=8) - @assert N>=0 + @assert N >= 0 sc = norm(h) - h = h*sqrt(2)/sc + h = h * sqrt(2) / sc phi = copy(h) psi = mirror(reverse(h)) - for i=1:N + for i = 1:N phi = conv(upsample(phi), h) psi = conv(upsample(psi), h) end phi = phi[1:end-2^(N)+1] psi = psi[1:end-2^(N)+1] - return rmul!(phi,sc/sqrt(2)), rmul!(psi,sc/sqrt(2)), range(0, stop=length(h)-1, length=length(psi)) + return rmul!(phi, sc / sqrt(2)), rmul!(psi, sc / sqrt(2)), range(0, stop=length(h) - 1, length=length(psi)) end """ @@ -465,35 +383,35 @@ function testfunction(n::Int, ft::AbstractString) f = Vector{Float64}(undef, n) range = 0:1/n:1-eps() i = 1 - if ft=="Blocks" - tj = [0.1,0.13,0.15,0.23,0.25,0.4,0.44,0.65,0.76,0.78,0.81] - hj = [4,-5,3,-4,5,-4.2,2.1,4.3,-3.1,2.1,-4.2] + if ft == "Blocks" + tj = [0.1, 0.13, 0.15, 0.23, 0.25, 0.4, 0.44, 0.65, 0.76, 0.78, 0.81] + hj = [4, -5, 3, -4, 5, -4.2, 2.1, 4.3, -3.1, 2.1, -4.2] for t in range f[i] = 0 for j = 1:11 - f[i] += hj[j]*(1+sign(t-tj[j]))/2 + f[i] += hj[j] * (1 + sign(t - tj[j])) / 2 end i += 1 end - elseif ft=="Bumps" - tj = [0.1,0.13,0.15,0.23,0.25,0.4,0.44,0.65,0.76,0.78,0.81] - hj = [4,5,3,4,5,4.2,2.1,4.3,3.1,5.1,4.2] - wj = [0.005,0.005,0.006,0.01,0.01,0.03,0.01,0.01,0.005,0.008,0.005] + elseif ft == "Bumps" + tj = [0.1, 0.13, 0.15, 0.23, 0.25, 0.4, 0.44, 0.65, 0.76, 0.78, 0.81] + hj = [4, 5, 3, 4, 5, 4.2, 2.1, 4.3, 3.1, 5.1, 4.2] + wj = [0.005, 0.005, 0.006, 0.01, 0.01, 0.03, 0.01, 0.01, 0.005, 0.008, 0.005] for t in range f[i] = 0 for j = 1:11 - f[i] += hj[j]/(1+abs((t-tj[j])/wj[j]))^4 + f[i] += hj[j] / (1 + abs((t - tj[j]) / wj[j]))^4 end i += 1 end - elseif ft=="HeaviSine" + elseif ft == "HeaviSine" for t in range - f[i] = 4*sin(4*pi*t) - sign(t-0.3) - sign(0.72-t) + f[i] = 4 * sin(4 * pi * t) - sign(t - 0.3) - sign(0.72 - t) i += 1 end - elseif ft=="Doppler" + elseif ft == "Doppler" for t in range - f[i] = sqrt(t*(1-t))*sin(2*pi*1.05/(t+0.05)) + f[i] = sqrt(t * (1 - t)) * sin(2 * pi * 1.05 / (t + 0.05)) i += 1 end else @@ -501,5 +419,3 @@ function testfunction(n::Int, ft::AbstractString) end return f end - -end # module diff --git a/src/WT/WT.jl b/src/WT/WT.jl new file mode 100644 index 0000000..697be8f --- /dev/null +++ b/src/WT/WT.jl @@ -0,0 +1,13 @@ +module WT + +include("wt.jl") + +export + DiscreteWavelet, + FilterWavelet, + LSWavelet, + OrthoFilter, + GLS, + wavelet + +end diff --git a/src/WT/wt.jl b/src/WT/wt.jl new file mode 100644 index 0000000..547e8f8 --- /dev/null +++ b/src/WT/wt.jl @@ -0,0 +1,454 @@ +using ..WT +using ..Util + +import Base.length +using SpecialFunctions +using LinearAlgebra + + +# TYPE HIERARCHY + +abstract type DiscreteWavelet{T} end +abstract type ContinuousWavelet{Boundary,T} end +# discrete transforms via filtering +abstract type FilterWavelet{T} <: DiscreteWavelet{T} end +# discrete transforms via lifting +abstract type LSWavelet{T} <: DiscreteWavelet{T} end +# all wavelet types +#const WaveletTransformType = Union{DiscreteWavelet, ContinuousWavelet} + +"""Get wavelet type name.""" +function name(::DiscreteWavelet) end +"""Get wavelet filter length.""" +function length(::FilterWavelet) end + +struct FilterTransform end +struct LiftingTransform end +"""Transform by filtering.""" +const Filter = FilterTransform() +"""Transform by lifting.""" +const Lifting = LiftingTransform() + + +# BOUNDARY TYPES + +abstract type WaveletBoundary end +# periodic (default) +struct PerBoundary <: WaveletBoundary end +# zero padding +struct ZPBoundary <: WaveletBoundary end +# constant padding +#struct CPBoundary <: WaveletBoundary end +struct NullBoundary <: WaveletBoundary end +# symmetric boundary (as in the DCTII) +struct SymBoundary <: WaveletBoundary end +# and so on... + +const Periodic = PerBoundary() +const DEFAULT_BOUNDARY = PerBoundary() +const padded = ZPBoundary() +const NaivePer = NullBoundary() +const SymBound = SymBoundary() + +# WAVELET CLASSES + +""" +The `WaveletClass` type has subtypes `OrthoWaveletClass`, + `BiOrthoWaveletClass`, and `ContinuousWaveletClass` + +The `WT` module has for convenience constants defined named +as the class short name and optionally amended with a number +specifing the number of vanishing moments. + +A class can also be explicitly constructed as e.g. `Daubechies{4}()`. + +# Examples +`WT.db2`, `WT.haar`, `WT.cdf97` + +**See also:** `WT.class`, `WT.name`, `WT.vanishingmoments` +""" +abstract type WaveletClass end +abstract type OrthoWaveletClass <: WaveletClass end +abstract type BiOrthoWaveletClass <: WaveletClass end +abstract type ContinuousWaveletClass <: WaveletClass end + +# Single classes, orthogonal +for (TYPE, NAMEBASE, MOMENTS) in ( + (:Haar, "haar", 1), + (:Beylkin, "beyl", -1), # TODO moments + (:Vaidyanathan, "vaid", -1), # TODO moments +) + @eval begin + struct $TYPE <: OrthoWaveletClass end + class(::$TYPE) = $(string(TYPE)) + name(::$TYPE) = string($NAMEBASE) + vanishingmoments(::$TYPE) = $MOMENTS + end + CONSTNAME = Symbol(NAMEBASE) + @eval begin + const $CONSTNAME = $TYPE() # type shortcut + end +end + +# discrete ortho classes +for (TYPE, NAMEBASE, RANGE) in ( + (:Daubechies, "db", 1:10), + (:Coiflet, "coif", 2:2:8), + (:Symlet, "sym", 4:10), + (:Battle, "batt", 2:2:6), +) + @eval begin + struct $TYPE{N} <: OrthoWaveletClass end + class(::$TYPE) = $(string(TYPE)) + name(::$TYPE{N}) where {N} = string($NAMEBASE, N) + vanishingmoments(::$TYPE{N}) where {N} = N + end + for NUM in RANGE + CONSTNAME = Symbol(string(NAMEBASE, NUM)) + @eval begin + const $CONSTNAME = $TYPE{$NUM}() # type shortcut + end + end +end + +# Parameterized BiOrtho classes +for (TYPE, NAMEBASE, RANGE1, RANGE2) in ( + (:CDF, "cdf", [9], [7]), +) + @eval begin + struct $TYPE{N1,N2} <: BiOrthoWaveletClass end + class(::$TYPE) = $(string(TYPE)) + name(::$TYPE{N1,N2}) where {N1,N2} = string($NAMEBASE, N1, "/", N2) + vanishingmoments(::$TYPE{N1,N2}) where {N1,N2} = (N1, N2) + end + for i in length(RANGE1) + CONSTNAME = Symbol(string(NAMEBASE, RANGE1[i], RANGE2[i])) + @eval begin + const $CONSTNAME = $TYPE{$RANGE1[$i],$RANGE2[$i]}() # type shortcut + end + end +end + + + +# IMPLEMENTATIONS OF FilterWavelet + +""" +Wavelet type for discrete orthogonal transforms by filtering. + +**See also:** `GLS`, `wavelet` +""" +struct OrthoFilter{T<:WaveletBoundary} <: FilterWavelet{T} + qmf::Vector{Float64} # quadrature mirror filter + name::String # filter short name +end + +function OrthoFilter(w::WC, ::T=DEFAULT_BOUNDARY) where {WC<:OrthoWaveletClass,T<:WaveletBoundary} + name = WT.name(w) + if WC <: Daubechies + qmf = daubechies(vanishingmoments(w)) + else + qmf = get(FILTERS, name, nothing) + qmf == nothing && throw(ArgumentError("filter not found")) + end + # make sure it is normalized in l2-norm + return OrthoFilter{T}(qmf ./ norm(qmf), name) +end + +length(f::OrthoFilter) = length(f.qmf) +qmf(f::OrthoFilter) = f.qmf +name(f::OrthoFilter) = f.name + +"""Scale filter by scalar.""" +function scale(f::OrthoFilter{T}, a::Number) where {T<:WaveletBoundary} + return OrthoFilter{T}(f.qmf .* a, f.name) +end + +"""Quadrature mirror filter pair.""" +function makeqmfpair(f::OrthoFilter, fw::Bool=true, T::Type=eltype(qmf(f))) + scfilter, dcfilter = makereverseqmfpair(f, fw, T) + return reverse(scfilter), reverse(dcfilter) +end + +"""Reversed quadrature mirror filter pair.""" +function makereverseqmfpair(f::OrthoFilter, fw::Bool=true, T::Type=eltype(qmf(f))) + h = convert(Vector{T}, qmf(f)) + if fw + scfilter = reverse(h) + dcfilter = mirror(h) + else + scfilter = h + dcfilter = reverse(mirror(h)) + end + return scfilter, dcfilter +end + +#struct BiOrthoFilter{T<:WaveletBoundary} <: FilterWavelet{T} +# qmf1::Vector{Float64} # quadrature mirror filter 1 +# qmf2::Vector{Float64} # quadrature mirror filter 2 +# name::String # filter short name +# BiOrthoFilter(qmf1, qmf2, name) = new(qmf1, qmf2, name) +#end + + +# IMPLEMENTATIONS OF LSWavelet + +abstract type StepType end +struct PredictStep <: StepType end +struct UpdateStep <: StepType end +const Predict = PredictStep() +const Update = UpdateStep() + +struct LSStepParam{T<:Number} + coef::Vector{T} # lifting coefficients + shift::Int # + left shift, - right shift +end + +struct LSStep{T<:Number} + param::LSStepParam{T} + steptype::StepType +end + +function LSStep(st::StepType, coef::Vector{T}, shift::Int) where {T} + return LSStep{T}(LSStepParam{T}(coef, shift), st) +end + +length(s::LSStep) = length(s.param) +length(s::LSStepParam) = length(s.coef) + +""" +Wavelet type for discrete general (bi)orthogonal transforms +by using a lifting scheme. + +**See also:** `OrthoFilter`, `wavelet` +""" +struct GLS{T<:WaveletBoundary} <: LSWavelet{T} + step::Vector{LSStep{Float64}} # steps to be taken + norm1::Float64 # normalization of scaling coefs. + norm2::Float64 # normalization of detail coefs. + name::String # name of scheme +end + +function GLS(w::WC, ::T=DEFAULT_BOUNDARY) where {WC<:WaveletClass,T<:WaveletBoundary} + name = WT.name(w) + schemedef = get(SCHEMES, name, nothing) + schemedef == nothing && throw(ArgumentError("scheme not found")) + return GLS{T}(schemedef[1], schemedef[2], schemedef[3], name) +end + +name(s::GLS) = s.name + +# TRANSFORM TYPE CONSTRUCTORS + +""" + wavelet(c[, t=WT.Filter][, boundary=WT.Periodic][, s=Real]) + +Construct wavelet type where `c` is a wavelet class, +`t` is the transformation type (`WT.Filter` or `WT.Lifting`), +and `boundary` is the type of boundary treatment. In the continuous case, s>0 +is the number of wavelets between the octaves ``2^J`` and ``2^{J+1}`` +(defaults to 8, which is most appropriate for music and other audio). + +# Examples +```julia +wavelet(WT.coif6) +wavelet(WT.db1, WT.Lifting) +wavelet(WT.morl,4) +``` + +**See also:** `WT.WaveletClass` +""" +function wavelet end + +wavelet(c::K, boundary::WaveletBoundary=DEFAULT_BOUNDARY) where {K<:Union{BiOrthoWaveletClass,OrthoWaveletClass}} = wavelet(c, Filter, boundary) +wavelet(c::OrthoWaveletClass, t::FilterTransform, boundary::WaveletBoundary=DEFAULT_BOUNDARY) = OrthoFilter(c, boundary) +wavelet(c::WaveletClass, t::LiftingTransform, boundary::WaveletBoundary=DEFAULT_BOUNDARY) = GLS(c, boundary) + + +# ------------------------------------------------------------ + +# Compute filters from the Daubechies class +# N is the number of zeros at -1 +function daubechies(N::Int) + @assert N > 0 + # Create polynomial + C = Vector{Int}(undef, N) + @inbounds for n = 0:N-1 + C[N-n] = binomial(N - 1 + n, n) + end + + # Find roots in y domain (truncated binomial series; (1 - y)^{-N}) + Y = roots(C) + + # Find roots in z domain: + # z + z^{-1} = 2 - 4*y + # where y is a root from above + Z = zeros(ComplexF64, 2 * N - 2) + @inbounds for i = 1:N-1 + Yi = Y[i] + d = 2 * sqrt(Yi * Yi - Yi) + y2 = 1 - 2 * Yi + Z[i] = y2 + d + Z[i+N-1] = y2 - d + end + + # Retain roots inside unit circle + nr = 0 # count roots + @inbounds for i = eachindex(Z) + if abs(Z[i]) <= 1 + eps() + nr += 1 + end + end + + # Find coefficients of the polynomial + # (1 + z)^N * \prod_i (z - z_i) + R = Vector{ComplexF64}(undef, N + nr) + @inbounds for i = 1:N + R[i] = -1 + end + k = N + @inbounds for i = eachindex(Z) + if abs(Z[i]) <= 1 + eps() + k += 1 + R[k] = Z[i] + end + end + HH = vieta(R) + + # Normalize coefficients + rmul!(HH, 1 / norm(HH)) + return real(HH) +end + +# Compute roots of polynomial +# Input is a coefficient vector with highest powers first +function roots(C::AbstractVector) + A = compan(C) + return eigvals(A) +end + +# Create companion matrix for a polynomial +# Input is a coefficient vector with highest powers first +function compan(C::AbstractVector) + n = length(C) + A = zeros(n - 1, n - 1) + + if n > 1 + @inbounds A[1, :] .= -C[2:end] ./ C[1] + @inbounds A[2:n:end] .= 1 + end + return A +end + +# Vieta-like formula for computing polynomial coefficients from roots +# See +# http://www.mathworks.se/help/matlab/ref/poly.html +function vieta(R::AbstractVector) + n = length(R) + C = zeros(ComplexF64, n + 1) + C[1] = 1 + Ci::ComplexF64 = 0 + Cig::ComplexF64 = 0 + + @inbounds for k = 1:n + Ci = C[1] + for i = 1:k + Cig = C[i+1] + C[i+1] = Cig - R[k] * Ci + Ci = Cig + end + end + return C +end + + +# scaling filters h (low pass) +# the number at end of a filter name is the +# number of vanishing moments of the mother wavelet function +# sources: +# http://statweb.stanford.edu/~wavelab/ (Orthogonal/MakeONFilter.m) +# http://www.mathworks.com/matlabcentral/fileexchange/5502-filter-coefficients-to-popular-wavelets +### https://github.com/nigma/pywt/blob/master/src/wavelets_coeffs.template.h +# name => qmf +const FILTERS = Dict{String,Vector{Float64}}( + # Haar filter + "haar" => + [0.7071067811865475, 0.7071067811865475], + # Daubechies filters, see daubechies() + + # Coiflet filters + "coif2" => + [-0.072732619513, 0.337897662458, 0.852572020212, 0.384864846864, -0.072732619513, -0.015655728135], + "coif4" => + [0.0163873364635998, -0.0414649367819558, -0.0673725547222826, 0.3861100668229939, 0.8127236354493977, 0.4170051844236707, -0.0764885990786692, -0.0594344186467388, 0.0236801719464464, 0.0056114348194211, -0.0018232088707116, -0.0007205494453679], + "coif6" => + [-0.0037935128644910, 0.0077825964273254, 0.0234526961418362, -0.0657719112818552, -0.0611233900026726, 0.4051769024096150, 0.7937772226256169, 0.4284834763776168, -0.0717998216193117, -0.0823019271068856, 0.0345550275730615, 0.0158805448636158, -0.0090079761366615, -0.0025745176887502, 0.0011175187708906, 0.0004662169601129, -0.0000709833031381, -0.0000345997728362], + "coif8" => + [0.0008923136685824, -0.0016294920126020, -0.0073461663276432, 0.0160689439647787, 0.0266823001560570, -0.0812666996808907, -0.0560773133167630, 0.4153084070304910, 0.7822389309206135, 0.4343860564915321, -0.0666274742634348, -0.0962204420340021, 0.0393344271233433, 0.0250822618448678, -0.0152117315279485, -0.0056582866866115, 0.0037514361572790, 0.0012665619292991, -0.0005890207562444, -0.0002599745524878, 0.0000623390344610, 0.0000312298758654, -0.0000032596802369, -0.0000017849850031], + "coif10" => + [-0.0002120808398259, 0.0003585896879330, 0.0021782363583355, -0.0041593587818186, -0.0101311175209033, 0.0234081567882734, 0.0281680289738655, -0.0919200105692549, -0.0520431631816557, 0.4215662067346898, 0.7742896037334738, 0.4379916262173834, -0.0620359639693546, -0.1055742087143175, 0.0412892087544753, 0.0326835742705106, -0.0197617789446276, -0.0091642311634348, 0.0067641854487565, 0.0024333732129107, -0.0016628637021860, -0.0006381313431115, 0.0003022595818445, 0.0001405411497166, -0.0000413404322768, -0.0000213150268122, 0.0000037346551755, 0.0000020637618516, -0.0000001674428858, -0.0000000951765727], + # Symmlet filter + "sym4" => + [0.0455703458960000, -0.0178247014420000, -0.1403176241790000, 0.4212345342040000, 1.1366582434079999, 0.7037390686560000, -0.0419109651250000, -0.1071489014180000], + "sym5" => + [0.0276321529580000, -0.0298424998690000, -0.2479513626130000, 0.0234789231360000, 0.8965816483800000, 1.0230529668940000, 0.2819906968540000, -0.0553441861170000, 0.0417468644220000, 0.0386547959550000], + "sym6" => + [-0.0110318675090000, 0.0024999220930000, 0.0632505626600000, -0.0297837512990000, -0.1027249698620000, 0.4779043713330000, 1.1138927839260000, 0.6944579729580000, -0.0683231215870000, -0.1668632154120000, 0.0049366123720000, 0.0217847003270000], + "sym7" => + [0.0145213947620000, 0.0056713426860000, -0.1524638718960000, -0.1980567068070000, 0.4081839397250000, 1.0857827098140000, 0.7581626019640000, 0.0246656594890000, -0.0700782912220000, 0.0960147679360000, 0.0431554525820000, -0.0178704316510000, -0.0014812259150000, 0.0037926585340000], + "sym8" => + [-0.0047834585120000, -0.0007666908960000, 0.0448236230420000, 0.0107586117510000, -0.2026486552860000, -0.0866536154060000, 0.6807453471900000, 1.0991066305370001, 0.5153986703740000, -0.0734625087610000, -0.0384935212630000, 0.0694904659110000, 0.0053863887540000, -0.0211456865280000, -0.0004283943000000, 0.0026727933930000], + "sym9" => + [0.0019811937360000, 0.0008765025390000, -0.0187693968360000, -0.0163033512260000, 0.0427444336020000, 0.0008251409290000, -0.0771721610970000, 0.3376589236020000, 1.0152597908320000, 0.8730484073490000, 0.0498828309590000, -0.2708937835030000, -0.0257864459300000, 0.0877912515540000, 0.0125288962420000, -0.0145155785530000, -0.0006691415090000, 0.0015124873090000], + "sym10" => + [-0.0006495898960000, 0.0000806612040000, 0.0064957283750000, -0.0011375353140000, -0.0287862319260000, 0.0081528167990000, 0.0707035675500000, -0.0452407722180000, -0.0502565400920000, 0.5428130112130000, 1.0882515305000000, 0.6670713381540000, -0.1002402150310000, -0.2255589722340000, 0.0164188694260000, 0.0649509245790000, -0.0020723639230000, -0.0122206426300000, 0.0001352450200000, 0.0010891704470000], + # Battle-Lemarie filter + "batt2" => + [-0.0000867523000000, -0.0001586010000000, 0.0003617810000000, 0.0006529220000000, -0.0015570100000000, -0.0027458800000000, 0.0070644200000000, 0.0120030000000000, -0.0367309000000000, -0.0488618000000000, 0.2809310000000000, 0.5781630000000000, 0.2809310000000000, -0.0488618000000000, -0.0367309000000000, 0.0120030000000000, 0.0070644200000000, -0.0027458800000000, -0.0015570100000000, 0.0006529220000000, 0.0003617810000000, -0.0001586010000000, -0.0000867523000000], + "batt4" => + [0.0001033070000000, -0.0001642640000000, -0.0002018180000000, 0.0003267490000000, 0.0003959460000000, -0.0006556200000000, -0.0007804680000000, 0.0013308600000000, 0.0015462400000000, -0.0027452900000000, -0.0030786300000000, 0.0057993200000000, 0.0061414300000000, -0.0127154000000000, -0.0121455000000000, 0.0297468000000000, 0.0226846000000000, -0.0778079000000000, -0.0354980000000000, 0.3068300000000000, 0.5417360000000000, 0.3068300000000000, -0.0354980000000000, -0.0778079000000000, 0.0226846000000000, 0.0297468000000000, -0.0121455000000000, -0.0127154000000000, 0.0061414300000000, 0.0057993200000000, -0.0030786300000000, -0.0027452900000000, 0.0015462400000000, 0.0013308600000000, -0.0007804680000000, -0.0006556200000000, 0.0003959460000000, 0.0003267490000000, -0.0002018180000000, -0.0001642640000000, 0.0001033070000000], + "batt6" => + [0.0001011130000000, 0.0001107090000000, -0.0001591680000000, -0.0001726850000000, 0.0002514190000000, 0.0002698420000000, -0.0003987590000000, -0.0004224850000000, 0.0006355630000000, 0.0006628360000000, -0.0010191200000000, -0.0010420700000000, 0.0016465900000000, 0.0016413200000000, -0.0026864600000000, -0.0025881600000000, 0.0044400200000000, 0.0040788200000000, -0.0074684800000000, -0.0063988600000000, 0.0128754000000000, 0.0099063500000000, -0.0229951000000000, -0.0148537000000000, 0.0433544000000000, 0.0208414000000000, -0.0914068000000000, -0.0261771000000000, 0.3128690000000000, 0.5283740000000000, 0.3128690000000000, -0.0261771000000000, -0.0914068000000000, 0.0208414000000000, 0.0433544000000000, -0.0148537000000000, -0.0229951000000000, 0.0099063500000000, 0.0128754000000000, -0.0063988600000000, -0.0074684800000000, 0.0040788200000000, 0.0044400200000000, -0.0025881600000000, -0.0026864600000000, 0.0016413200000000, 0.0016465900000000, -0.0010420700000000, -0.0010191200000000, 0.0006628360000000, 0.0006355630000000, -0.0004224850000000, -0.0003987590000000, 0.0002698420000000, 0.0002514190000000, -0.0001726850000000, -0.0001591680000000, 0.0001107090000000, 0.0001011130000000], + # Beylkin filter + "beyl" => + [0.0993057653740000, 0.4242153608130000, 0.6998252140570000, 0.4497182511490000, -0.1109275983480000, -0.2644972314460000, 0.0269003088040000, 0.1555387318770000, -0.0175207462670000, -0.0885436306230000, 0.0196798660440000, 0.0429163872740000, -0.0174604086960000, -0.0143658079690000, 0.0100404118450000, 0.0014842347820000, -0.0027360316260000, 0.0006404853290000], + # Vaidyanathan filter + "vaid" => + [-0.0000629061180000, 0.0003436319050000, -0.0004539566200000, -0.0009448971360000, 0.0028438345470000, 0.0007081375040000, -0.0088391034090000, 0.0031538470560000, 0.0196872150100000, -0.0148534480050000, -0.0354703986070000, 0.0387426192930000, 0.0558925236910000, -0.0777097509020000, -0.0839288843660000, 0.1319716614170000, 0.1350842271290000, -0.1944504717660000, -0.2634948024880000, 0.2016121617750000, 0.6356010598720000, 0.5727977932110000, 0.2501841295050000, 0.0457993341110000]) + + +# biortho filters +# name => (qmf1,qmf2) +const BIFILTERS = Dict{String,NTuple{2,Vector{Float64}}}( + # test + "test" => + ([0.7071067811865475, 0.7071067811865475], + [0.7071067811865475, 0.7071067811865475])) + + +# in matlab (d,p) -> (predict, update) +const SCHEMES = Dict{String,NTuple{3,Any}}( + # cdf 5/3 -> bior 2.2, cdf 9/7 -> bior 4.4 + # Cohen-Daubechies-Feauveau [Do Quan & Yo-Sung Ho. Optimized median lifting scheme for lossy image compression.] + "cdf9/7" => ([LSStep(Update, [1.0, 1.0] * 1.5861343420604, 0), + LSStep(Predict, [1.0, 1.0] * 0.05298011857291494, 1), + LSStep(Update, [1.0, 1.0] * -0.882911075531393, 0), + LSStep(Predict, [1.0, 1.0] * -0.44350685204384654, 1)], + 1.1496043988603355, + 0.8698644516247099), + + # Haar, same as db1 + "haar" => ([LSStep(Predict, [-1.0], 0), + LSStep(Update, [0.5], 0)], + 0.7071067811865475, + 1.4142135623730951), + # Daubechies + "db1" => ([LSStep(Predict, [-1.0], 0), + LSStep(Update, [0.5], 0)], + 0.7071067811865475, + 1.4142135623730951), + "db2" => ([LSStep(Predict, [-1.7320508075688772], 0), + LSStep(Update, [-0.0669872981077807, 0.4330127018922193], 1), + LSStep(Predict, [1.0], -1)], + 0.5176380902050414, + 1.9318516525781364)) diff --git a/src/Wavelets.jl b/src/Wavelets.jl index 60eb86e..316f652 100644 --- a/src/Wavelets.jl +++ b/src/Wavelets.jl @@ -1,15 +1,66 @@ -__precompile__() - module Wavelets -using SpecialFunctions -include("mod/Util.jl") -include("mod/WT.jl") -include("mod/Transforms.jl") -include("mod/Threshold.jl") -include("mod/Plot.jl") -using Reexport -@reexport using .Util, .WT, .Transforms, .Threshold, .Plot +include("Util/Util.jl") +include("WT/WT.jl") +include("Transforms/Transforms.jl") +include("Threshold/Threshold.jl") +include("Plot/Plot.jl") + +using .Util +using .WT +using .Transforms +using .Threshold +using .Plot + + +export + # Submodules (except Util => would lead to conflict with DSP.jl (?) + WT, Transforms, Threshold, Plot, + + + # Util submodule : + # dyadic + dyadicdetailindex, dyadicdetailrange, dyadicscalingrange, + dyadicdetailn, ndyadicscales, maxdyadiclevel, + tl2dyadiclevel, dyadiclevel2tl, + + # non-dyadic + detailindex, detailrange, detailn, + maxtransformlevels, maxmodwttransformlevels, + + # other util functions + mirror, upsample, downsample, iscube, isdyadic, + sufficientpoweroftwo, wcount, stridedcopy!, + isvalidtree, maketree, makewavelet, testfunction, + + + # WT submodule : + DiscreteWavelet, FilterWavelet, LSWavelet, OrthoFilter, GLS, wavelet, + + + # Transforms submodule : + dwt, idwt, dwt!, idwt!, wpt, iwpt, wpt!, iwpt!, modwt, imodwt, + + + # Threshold submodule : + + # threshold + threshold!, threshold, HardTH, SoftTH, SemiSoftTH, SteinTH, + BiggestTH, PosTH, NegTH, + + # denoising + DNFT, VisuShrink, denoise, noisest, + + # basis functions + matchingpursuit, bestbasistree, + + # entropy + coefentropy, Entropy, ShannonEntropy, LogEnergyEntropy, + + + # Plot submodule : + wplotdots, wplotim + end diff --git a/src/mod/Threshold.jl b/src/mod/Threshold.jl deleted file mode 100644 index d985984..0000000 --- a/src/mod/Threshold.jl +++ /dev/null @@ -1,462 +0,0 @@ -module Threshold -export - # threshold - threshold!, - threshold, - HardTH, - SoftTH, - SemiSoftTH, - SteinTH, - BiggestTH, - PosTH, - NegTH, - # denoising - DNFT, - VisuShrink, - denoise, - noisest, - # basis functions - matchingpursuit, - bestbasistree, - # entropy - coefentropy, - Entropy, - ShannonEntropy, - LogEnergyEntropy -using ..Util, ..WT, ..Transforms - -using LinearAlgebra -using Statistics - -# THRESHOLD TYPES AND FUNCTIONS - -abstract type THType end -struct HardTH <: THType end -struct SoftTH <: THType end -struct SemiSoftTH <: THType end -struct SteinTH <: THType end -struct BiggestTH <: THType end -struct PosTH <: THType end -struct NegTH <: THType end - -const DEFAULT_TH = HardTH() - -# biggest m-term approximation (best m-term approximation for orthogonal transforms) -# result is m-sparse -function threshold!(x::AbstractArray{<:Number}, TH::BiggestTH, m::Int) - @assert m >= 0 - n = length(x) - m > n && (m = n) - ind = sortperm(x, alg=QuickSort, by=abs) - @inbounds begin - for i = 1:n-m - x[ind[i]] = 0 - end - end - return x -end - -# hard -function threshold!(x::AbstractArray{<:Number}, TH::HardTH, t::Real) - @assert t >= 0 - @inbounds begin - for i in eachindex(x) - if abs(x[i]) <= t - x[i] = 0 - end - end - end - return x -end - -# soft -function threshold!(x::AbstractArray{<:Number}, TH::SoftTH, t::Real) - @assert t >= 0 - @inbounds begin - for i in eachindex(x) - sh = abs(x[i]) - t - if sh < 0 - x[i] = 0 - else - x[i] = sign(x[i])*sh - end - end - end - return x -end - -# semisoft -function threshold!(x::AbstractArray{<:Number}, TH::SemiSoftTH, t::Real) - @assert t >= 0 - @inbounds begin - for i in eachindex(x) - if x[i] <= 2*t - sh = abs(x[i]) - t - if sh < 0 - x[i] = 0 - elseif sh - t < 0 - x[i] = sign(x[i])*sh*2 - end - end - end - end - return x -end - -# stein -function threshold!(x::AbstractArray{<:Number}, TH::SteinTH, t::Real) - @assert t >= 0 - @inbounds begin - for i in eachindex(x) - sh = 1 - t*t/(x[i]*x[i]) - if sh < 0 - x[i] = 0 - else - x[i] = x[i]*sh - end - end - end - return x -end - -# shrink negative elements to 0 -function threshold!(x::AbstractArray{<:Number}, TH::NegTH) - @inbounds begin - for i in eachindex(x) - if x[i] < 0 - x[i] = 0 - end - end - end - return x -end - -# shrink positive elements to 0 -function threshold!(x::AbstractArray{<:Number}, TH::PosTH) - @inbounds begin - for i in eachindex(x) - if x[i] > 0 - x[i] = 0 - end - end - end - return x -end - -# the non inplace functions -function threshold(x::AbstractArray{T}, TH::THType, t::Real) where T<:Number - y = Vector{T}(undef, size(x)) - return threshold!(copyto!(y,x), TH, t) -end -function threshold(x::AbstractArray{T}, TH::THType) where T<:Number - y = Vector{T}(undef, size(x)) - return threshold!(copyto!(y,x), TH) -end - - -# DENOISING - -abstract type DNFT end - -struct VisuShrink <: DNFT - th::THType # threshold type - t::Float64 # threshold for noise level sigma=1, use sigma*t in application - VisuShrink(th, t) = new(th, t) -end -# define type for signal length n -function VisuShrink(n::Int) - return VisuShrink(DEFAULT_TH, sqrt(2*log(n))) -end - -const DEFAULT_WAVELET = wavelet(WT.sym5, WT.Filter) # default wavelet type - -# denoise signal x by thresholding in wavelet space -# estnoise is (x::AbstractArray, wt::Union{DiscreteWavelet,Nothing}) -function denoise(x::AbstractArray, - wt::Union{DiscreteWavelet,Nothing}=DEFAULT_WAVELET; - L::Int=min(maxtransformlevels(x),6), - dnt::S=VisuShrink(size(x,1)), - estnoise::Function=noisest, - TI::Bool=false, - nspin::Union{Int,Tuple}=tuple([8 for i=1:ndims(x)]...) ) where S<:DNFT - iscube(x) || throw(ArgumentError("array must be square/cube")) - sigma = estnoise(x, wt) - - if TI - wt == nothing && error("TI not supported with wt=nothing") - y = zeros(eltype(x), size(x)) - xt = similar(x) - pns = prod(nspin) - - if ndims(x) == 1 - z = similar(x) - for i = 1:pns - shift = i - 1 - Util.circshift!(z, x, shift) - - Transforms.dwt_oop!(xt, z, wt, L) - threshold!(xt, dnt.th, sigma*dnt.t) - Transforms.idwt_oop!(z, xt, wt, L) - - Util.circshift!(xt, z, -shift) - arrayadd!(y, xt) - end - else # ndims > 1 - for i = 1:pns - shift = nspin2circ(nspin, i) - z = circshift(x, shift) - - Transforms.dwt_oop!(xt, z, wt, L) - threshold!(xt, dnt.th, sigma*dnt.t) - Transforms.idwt_oop!(z, xt, wt, L) - - z = circshift(z, -shift) - arrayadd!(y, z) - end - end - rmul!(y, 1/pns) - else # !TI - if wt == nothing - y = copy(x) - threshold!(y, dnt.th, sigma*dnt.t) - else - y = dwt(x, wt, L) - threshold!(y, dnt.th, sigma*dnt.t) - if isa(wt, GLS) - idwt!(y, wt, L) - else - y2 = idwt(y, wt, L) - y = y2 - end - end - end - - return y -end -# add z to y -function arrayadd!(y::AbstractArray, z::AbstractArray) - length(y) == length(z) || throw(DimensionMismatch("lengths must be equal")) - for i in eachindex(y) - @inbounds y[i] += z[i] - end - return y -end - - -# estimate the std. dev. of the signal noise, assuming Gaussian distribution -function noisest(x::AbstractArray, wt::Union{DiscreteWavelet,Nothing}=DEFAULT_WAVELET, L::Integer = 1) - if wt == nothing - y = x - else - y = dwt(x, wt, L) - end - dr = y[detailrange(y, L)] - return mad!(dr)/0.6745 -end -# Median absolute deviation -function mad!(y::AbstractArray) - m = median!(y) - for i in eachindex(y) - y[i] = abs(y[i]-m) - end - return median!(y) -end - -# convert index i to a circshift array starting at 0 shift -nspin2circ(nspin::Int, i::Int) = nspin2circ((nspin,), i) -function nspin2circ(nspin::Tuple, i::Int) - c1 = CartesianIndices(nspin)[i].I - c = Vector{Int}(undef, length(c1)) - for k in 1:length(c1) - c[k] = c1[k]-1 - end - return c -end - - -# BASIS FUNCTIONS - -# Matching Pursuit -# see: Mallat (2009) p.642 "A wavelet tour of signal processing" -# find sparse vector y such that ||x - f(y)|| < tol approximately -# f is the operation of a M by N (M= -1 - @assert tol > 0 - r = x - n = 1 - - if !oop - y = zeros(eltype(x), length(ft(x))) - else # out of place functions f and ft - y = zeros(eltype(x), N) - tmp = similar(x, N) - ftr = similar(x, N) - aphi = similar(x, length(x)) - end - spat = zeros(eltype(x), length(y)) # sparse for atom computation - nmax == -1 && (nmax = length(y)) - - while norm(r) > tol && n <= nmax - # find largest inner product - !oop && (ftr = ft(r)) - oop && ft(ftr, r, tmp) - i = findmaxabs(ftr) - - # project on i-th atom - spat[i] = ftr[i] - !oop && (aphi = f(spat)) - oop && f(aphi, spat, tmp) - spat[i] = 0 - - # update residual, r = r - aphi - broadcast!(-, r, r, aphi) - - y[i] += ftr[i] - n += 1 - end - return y -end -function findmaxabs(x::AbstractVector) - m = abs(x[1]) - k = 1 - @inbounds for i in eachindex(x) - if abs(x[i]) > m - k = i - m = abs(x[i]) - end - end - return k -end - - -# ENTROPY MEASURES - -abstract type Entropy end -struct ShannonEntropy <: Entropy end #Coifman-Wickerhauser -struct LogEnergyEntropy <: Entropy end - -# Entropy measures: Additive with coefentropy(0) = 0 -# all coefs assumed to be on [-1,1] after normalization with nrm -# given x and y, where x has "more concentrated energy" than y -# then coefentropy(x, et, norm) <= coefentropy(y, et, norm) should be satisfied. - -function coefentropy(x::T, et::ShannonEntropy, nrm::T) where T<:AbstractFloat - s = (x/nrm)^2 - if s == 0.0 - return -zero(T) - else - return -s*log(s) - end -end -function coefentropy(x::T, et::LogEnergyEntropy, nrm::T) where T<:AbstractFloat - s = (x/nrm)^2 - if s == 0.0 - return -zero(T) - else - return -log(s) - end -end -function coefentropy(x::AbstractArray{T}, et::Entropy, nrm::T=norm(x)) where T<:AbstractFloat - @assert nrm >= 0 - sum = zero(T) - nrm == sum && return sum - for i in eachindex(x) - @inbounds sum += coefentropy(x[i], et, nrm) - end - return sum -end - - -# find the best tree that is a subset of the input tree (use :full to find the best tree) -# for wpt -function bestbasistree(y::AbstractVector{T}, wt::DiscreteWavelet, L::Integer=maxtransformlevels(y), et::Entropy=ShannonEntropy()) where T<:AbstractFloat - bestbasistree(y, wt, maketree(length(y), L, :full), et) -end -function bestbasistree(y::AbstractVector{T}, wt::DiscreteWavelet, tree::BitVector, et::Entropy=ShannonEntropy()) where T<:AbstractFloat - - isvalidtree(y, tree) || throw(ArgumentError("invalid tree")) - - tree[1] || copy(tree) # do nothing - - x = copy(y) - n = length(y) - tmp = Vector{T}(undef, n) - ntree = length(tree) - entr_bf = Vector{T}(undef, ntree) - nrm = norm(y) - - Lmax = maxtransformlevels(n) - L = Lmax - k = 1 - while L > 0 - ix = 1 - Lfw = Lmax-L - nj = detailn(n, Lfw) - - @assert nj <= n - dtmp = Transforms.unsafe_vectorslice(tmp, 1, nj) - while ix <= n - @assert nj+ix-1 <= n - dx = Transforms.unsafe_vectorslice(x, ix, nj) - - entr_bf[k] = coefentropy(dx, et, nrm) - - dwt!(dtmp, dx, wt, 1) - copyto!(dx, dtmp) - - ix += nj - k += 1 - end - L -= 1 - end - - # entropy of fully transformed signal (end nodes) - n_af = 2^(Lmax-1) - entr_af = Vector{T}(undef, n_af) - n_coef_af = div(n, n_af) - for i in 1:n_af - range = (i-1)*n_coef_af+1 : i*n_coef_af - entr_af[i] = coefentropy(x[range], et, nrm) - end - - # make the best tree - besttree = copy(tree) - for i in 1:ntree - if (i>1 && !besttree[i>>1]) || !tree[i] # parent is 0 or input tree-node is 0 - besttree[i] = false - else - if entr_bf[i] <= bestsubtree_entropy(entr_bf, entr_af, i) - besttree[i] = false - else - besttree[i] = true - end - end - end - - @assert isvalidtree(y, besttree) - return besttree::BitVector -end - -# the entropy of best subtree -function bestsubtree_entropy(entr_bf::Array, entr_af::Array, i::Int) - n = length(entr_bf) - n_af = length(entr_af) - @assert isdyadic(n+1) - @assert isdyadic(n_af) - @assert n + 1 == n_af<<1 - - if n < (i<<1) # bottom of tree - sum = entr_af[i - n_af + 1] - else - sum = bestsubtree_entropy(entr_bf, entr_af, i<<1) - sum += bestsubtree_entropy(entr_bf, entr_af, i<<1+1) - end - - lowestentropy = min(entr_bf[i], sum) - return lowestentropy -end - - -end # module diff --git a/src/mod/Transforms.jl b/src/mod/Transforms.jl deleted file mode 100644 index 7cee871..0000000 --- a/src/mod/Transforms.jl +++ /dev/null @@ -1,239 +0,0 @@ -module Transforms -export dwt, idwt, dwt!, idwt!, - wpt, iwpt, wpt!, iwpt!, - modwt, imodwt -using ..Util, ..WT -using FFTW - -# TODO Use StridedArray instead of AbstractArray where writing to array. -# TODO change integer dependent wavelets to parametric types (see "Value types", https://docs.julialang.org/en/v1/manual/types/index.html#%22Value-types%22-1) -const DWTArray = AbstractArray -const WPTArray = AbstractVector -const ValueType = Union{AbstractFloat, Complex} - -const FVector = StridedVector # e.g., work space vectors - -# DWT - -""" - dwt(x, wt[, L=maxtransformlevels(x)]) - -Perform a discrete wavelet transform of the array `x`. -The wavelet type `wt` determines the transform type -(filter or lifting) and the wavelet class, see `wavelet`. - -The number of transformation levels `L` can be any non-negative -integer such that the size of `x` is divisible by `L`. -Performs the identity transformation if `L==0`. - -# Examples -```julia -dwt(x, wavelet(WT.coif6)) -``` - -**See also:** `idwt`, `dwt!`, `wavelet` -""" -function dwt end - -""" - idwt(x, wt[, L=maxtransformlevels(x)]) - -Perform an inverse discrete wavelet transform of the array `x`, -the inverse of `dwt(x, wt, L)`. - -**See also:** `dwt`, `idwt!` -""" -function idwt end - -""" - dwt!(y, x, wt::OrthoFilter[, L=maxtransformlevels(x)]) - - dwt!(y, wt::GLS[, L=maxtransformlevels(x)]) - -Same as `dwt` but without array allocation. -Perform "out of place" transform with a filter, or -a inplace transform with a lifting scheme. The difference -between the filter and lifting methods is due to the -structure of the transform algorithms. - -**See also:** `idwt!` -""" -function dwt! end - -""" - idwt!(y, x, wt::OrthoFilter[, L=maxtransformlevels(x)]) - - idwt!(y, wt::GLS[, L=maxtransformlevels(x)]) - -The inverse of `dwt!`. - -**See also:** `dwt!` -""" -function idwt! end - - -# WPT - -""" - wpt - -Perform a discrete wavelet packet transform of the array `x`. -**See also:** `dwt`, `wavelet` -""" -function wpt end - -""" - iwpt - -Perform an inverse discrete wavelet packet transform of the array `x`. -**See also:** `idwt`, `wavelet` -""" -function iwpt end - -""" - wpt! - -Same as `wpt` but without array allocation. -""" -function wpt! end - -""" - iwpt! - -Same as `iwpt` but without array allocation. -""" -function iwpt! end - - -# DWT (discrete wavelet transform) - -for (Xwt, Xwt!, _Xwt!, fw) in ((:dwt, :dwt!, :_dwt!, true), - (:idwt, :idwt!, :_dwt!, false)) -@eval begin - # filter - function ($Xwt)(x::DWTArray{T}, filter::OrthoFilter, - L::Integer=maxtransformlevels(x)) where T<:ValueType - y = Array{T}(undef, size(x)) - return ($_Xwt!)(y, x, filter, L, $fw) - end - function ($Xwt!)(y::DWTArray{<:ValueType}, x::DWTArray{<:ValueType}, filter::OrthoFilter, - L::Integer=maxtransformlevels(x)) - return ($_Xwt!)(y, x, filter, L, $fw) - end - # lifting - function ($Xwt)(x::DWTArray{T}, scheme::GLS, - L::Integer=maxtransformlevels(x)) where T<:ValueType - y = Array{T}(undef, size(x)) - copyto!(y, x) - return ($_Xwt!)(y, scheme, L, $fw) - end - function ($Xwt!)(y::DWTArray{T}, scheme::GLS, - L::Integer=maxtransformlevels(y)) where T<:ValueType - return ($_Xwt!)(y, scheme, L, $fw) - end -end # begin -end # for - -# WPT (wavelet packet transform) - -for (Xwt, Xwt!, _Xwt!, fw) in ((:wpt, :wpt!, :_wpt!, true), - (:iwpt, :iwpt!, :_wpt!, false)) -@eval begin - function ($Xwt)(x::WPTArray{T}, wt::DiscreteWavelet, - L::Integer=maxtransformlevels(x)) where T<:ValueType - return ($Xwt)(x, wt, maketree(length(x), L, :full)) - end - # filter - function ($Xwt)(x::WPTArray{T}, filter::OrthoFilter, - tree::BitVector=maketree(x, :full)) where T<:ValueType - y = Array{T}(undef, size(x)) - return ($_Xwt!)(y, x, filter, tree, $fw) - end - function ($Xwt!)(y::WPTArray{T}, x::WPTArray{T}, - filter::OrthoFilter) where T<:ValueType - return ($Xwt!)(y, x, filter, maketree(x, :full)) - end - function ($Xwt!)(y::WPTArray{T}, x::WPTArray{T}, - filter::OrthoFilter, L::Integer) where T<:ValueType - return ($Xwt!)(y, x, filter, maketree(length(x), L, :full)) - end - function ($Xwt!)(y::WPTArray{T}, x::WPTArray{T}, - filter::OrthoFilter, tree::BitVector) where T<:ValueType - return ($_Xwt!)(y, x, filter, tree, $fw) - end - # lifting - function ($Xwt)(x::WPTArray{T}, scheme::GLS, - tree::BitVector=maketree(x, :full)) where T<:ValueType - y = Array{T}(undef, size(x)) - copyto!(y, x) - return ($_Xwt!)(y, scheme, tree, $fw) - end - function ($Xwt!)(y::WPTArray{T}, scheme::GLS) where T<:ValueType - return ($Xwt!)(y, scheme, maketree(y, :full)) - end - function ($Xwt!)(y::WPTArray{T}, scheme::GLS, L::Integer) where T<:ValueType - return ($Xwt!)(y, scheme, maketree(length(x), L, :full)) - end - function ($Xwt!)(y::WPTArray{T}, scheme::GLS, tree::BitVector) where T<:ValueType - return ($_Xwt!)(y, scheme, tree, $fw) - end -end # begin -end # for - - -# DWTC (column-wise discrete wavelet transform) -#dwtc(::AbstractArray, ::DiscreteWavelet) -#idwtc(::AbstractArray, ::DiscreteWavelet) - -# SWT (stationary wavelet transform) -# swt(::AbstractVector, ::DiscreteWavelet) -# iswt(::AbstractVector, ::DiscreteWavelet) - -# Int -> Float -for Xwt in (:dwt, :idwt, :dwtc, :idwtc, :wpt, :iwpt) - @eval $Xwt(x::AbstractArray{<:Integer}, args...) = $Xwt(float(x), args...) -end - -# non-exported "out of place" functions -for (Xwt_oop!, Xwt!) in ((:dwt_oop!, :dwt!), (:idwt_oop!, :idwt!)) -@eval begin - # filter - function ($Xwt_oop!)(y::DWTArray{T}, x::DWTArray{T}, filter::OrthoFilter, - L::Integer=maxtransformlevels(x)) where T<:ValueType - return ($Xwt!)(y, x, filter, L) - end - # lifting - function ($Xwt_oop!)(y::DWTArray{T}, x::DWTArray{T}, scheme::GLS, - L::Integer=maxtransformlevels(x)) where T<:ValueType - copyto!(y, x) - return ($Xwt!)(y, scheme, L) - end -end # begin -end # for - -# Array with shared memory -function unsafe_vectorslice(A::Array{T}, i::Int, n::Int) where T - return unsafe_wrap(Array, pointer(A, i), n)::Vector{T} -end -function unsafe_vectorslice(A::StridedArray{T}, i::Int, n::Int) where T - return @view A[i:(i-1+n)] -end - -# linear indices of start of rows/cols/planes -# 2-D, size(A) = (m,n) -row_idx(i, m) = i -col_idx(i, m) = 1 + (i-1)*m -# 3-D, size(A) = (m,n,d) -row_idx(i, j, m, n=m) = row_idx(i, n) + (j-1)*n*m -col_idx(i, j, m, n=m) = col_idx(i, m) + (j-1)*n*m -plane_idx(i, j, m) = i + (j-1)*m - -# filter transforms -include("transforms_filter.jl") - -# lifting transforms -include("transforms_lifting.jl") - -include("transforms_maximal_overlap.jl") - -end # module diff --git a/src/mod/WT.jl b/src/mod/WT.jl deleted file mode 100644 index 5d10fec..0000000 --- a/src/mod/WT.jl +++ /dev/null @@ -1,493 +0,0 @@ -module WT -export DiscreteWavelet, - FilterWavelet, - LSWavelet, - OrthoFilter, - GLS, - wavelet - - -using ..Util -import Base.length -using SpecialFunctions -using LinearAlgebra - - -# TYPE HIERARCHY - -abstract type DiscreteWavelet{T} end -abstract type ContinuousWavelet{Boundary,T} end -# discrete transforms via filtering -abstract type FilterWavelet{T} <: DiscreteWavelet{T} end -# discrete transforms via lifting -abstract type LSWavelet{T} <: DiscreteWavelet{T} end -# all wavelet types -#const WaveletTransformType = Union{DiscreteWavelet, ContinuousWavelet} - -"""Get wavelet type name.""" -function name(::DiscreteWavelet) end -"""Get wavelet filter length.""" -function length(::FilterWavelet) end - -struct FilterTransform end -struct LiftingTransform end -"""Transform by filtering.""" -const Filter = FilterTransform() -"""Transform by lifting.""" -const Lifting = LiftingTransform() - - -# BOUNDARY TYPES - -abstract type WaveletBoundary end -# periodic (default) -struct PerBoundary <: WaveletBoundary end -# zero padding -struct ZPBoundary <: WaveletBoundary end -# constant padding -#struct CPBoundary <: WaveletBoundary end -struct NullBoundary <: WaveletBoundary end -# symmetric boundary (as in the DCTII) -struct SymBoundary <: WaveletBoundary end -# and so on... - -const Periodic = PerBoundary() -const DEFAULT_BOUNDARY = PerBoundary() -const padded = ZPBoundary() -const NaivePer = NullBoundary() -const SymBound = SymBoundary() - -# WAVELET CLASSES - -""" -The `WaveletClass` type has subtypes `OrthoWaveletClass`, - `BiOrthoWaveletClass`, and `ContinuousWaveletClass` - -The `WT` module has for convenience constants defined named -as the class short name and optionally amended with a number -specifing the number of vanishing moments. - -A class can also be explicitly constructed as e.g. `Daubechies{4}()`. - -# Examples -`WT.db2`, `WT.haar`, `WT.cdf97` - -**See also:** `WT.class`, `WT.name`, `WT.vanishingmoments` -""" -abstract type WaveletClass end -abstract type OrthoWaveletClass <: WaveletClass end -abstract type BiOrthoWaveletClass <: WaveletClass end -abstract type ContinuousWaveletClass <: WaveletClass end - -# Single classes, orthogonal -for (TYPE, NAMEBASE, MOMENTS) in ( - (:Haar, "haar", 1), - (:Beylkin, "beyl", -1), # TODO moments - (:Vaidyanathan, "vaid",-1), # TODO moments - ) - @eval begin - struct $TYPE <: OrthoWaveletClass end - class(::$TYPE) = $(string(TYPE)) - name(::$TYPE) = string($NAMEBASE) - vanishingmoments(::$TYPE) = $MOMENTS - end - CONSTNAME = Symbol(NAMEBASE) - @eval begin - const $CONSTNAME = $TYPE() # type shortcut - end - end - -# discrete ortho classes -for (TYPE, NAMEBASE, RANGE) in ( - (:Daubechies, "db", 1:10), - (:Coiflet, "coif", 2:2:8), - (:Symlet, "sym", 4:10), - (:Battle, "batt", 2:2:6), - ) - @eval begin - struct $TYPE{N} <: OrthoWaveletClass end - class(::$TYPE) = $(string(TYPE)) - name(::$TYPE{N}) where N = string($NAMEBASE,N) - vanishingmoments(::$TYPE{N}) where N = N - end - for NUM in RANGE - CONSTNAME = Symbol(string(NAMEBASE, NUM)) - @eval begin - const $CONSTNAME = $TYPE{$NUM}() # type shortcut - end - end -end - -# Parameterized BiOrtho classes -for (TYPE, NAMEBASE, RANGE1, RANGE2) in ( - (:CDF, "cdf", [9], [7]), - ) - @eval begin - struct $TYPE{N1, N2} <: BiOrthoWaveletClass end - class(::$TYPE) = $(string(TYPE)) - name(::$TYPE{N1, N2}) where {N1, N2} = string($NAMEBASE,N1,"/",N2) - vanishingmoments(::$TYPE{N1, N2}) where {N1, N2} = (N1, N2) - end - for i in length(RANGE1) - CONSTNAME = Symbol(string(NAMEBASE,RANGE1[i],RANGE2[i])) - @eval begin - const $CONSTNAME = $TYPE{$RANGE1[$i],$RANGE2[$i]}() # type shortcut - end - end -end - - - -# IMPLEMENTATIONS OF FilterWavelet - -""" -Wavelet type for discrete orthogonal transforms by filtering. - -**See also:** `GLS`, `wavelet` -""" -struct OrthoFilter{T<:WaveletBoundary} <: FilterWavelet{T} - qmf ::Vector{Float64} # quadrature mirror filter - name ::String # filter short name -end - -function OrthoFilter(w::WC, ::T=DEFAULT_BOUNDARY) where {WC<:OrthoWaveletClass, T<:WaveletBoundary} - name = WT.name(w) - if WC <: Daubechies - qmf = daubechies(vanishingmoments(w)) - else - qmf = get(FILTERS, name, nothing) - qmf == nothing && throw(ArgumentError("filter not found")) - end - # make sure it is normalized in l2-norm - return OrthoFilter{T}(qmf./norm(qmf), name) -end - -length(f::OrthoFilter) = length(f.qmf) -qmf(f::OrthoFilter) = f.qmf -name(f::OrthoFilter) = f.name - -"""Scale filter by scalar.""" -function scale(f::OrthoFilter{T}, a::Number) where T<:WaveletBoundary - return OrthoFilter{T}(f.qmf.*a, f.name) -end - -"""Quadrature mirror filter pair.""" -function makeqmfpair(f::OrthoFilter, fw::Bool=true, T::Type=eltype(qmf(f))) - scfilter, dcfilter = makereverseqmfpair(f, fw, T) - return reverse(scfilter), reverse(dcfilter) -end - -"""Reversed quadrature mirror filter pair.""" -function makereverseqmfpair(f::OrthoFilter, fw::Bool=true, T::Type=eltype(qmf(f))) - h = convert(Vector{T}, qmf(f)) - if fw - scfilter = reverse(h) - dcfilter = mirror(h) - else - scfilter = h - dcfilter = reverse(mirror(h)) - end - return scfilter, dcfilter -end - -#struct BiOrthoFilter{T<:WaveletBoundary} <: FilterWavelet{T} -# qmf1::Vector{Float64} # quadrature mirror filter 1 -# qmf2::Vector{Float64} # quadrature mirror filter 2 -# name::String # filter short name -# BiOrthoFilter(qmf1, qmf2, name) = new(qmf1, qmf2, name) -#end - - -# IMPLEMENTATIONS OF LSWavelet - -abstract type StepType end -struct PredictStep <: StepType end -struct UpdateStep <: StepType end -const Predict = PredictStep() -const Update = UpdateStep() - -struct LSStepParam{T<:Number} - coef ::Vector{T} # lifting coefficients - shift ::Int # + left shift, - right shift -end - -struct LSStep{T<:Number} - param::LSStepParam{T} - steptype::StepType -end - -function LSStep(st::StepType, coef::Vector{T}, shift::Int) where T - return LSStep{T}(LSStepParam{T}(coef, shift), st) -end - -length(s::LSStep) = length(s.param) -length(s::LSStepParam) = length(s.coef) - -""" -Wavelet type for discrete general (bi)orthogonal transforms -by using a lifting scheme. - -**See also:** `OrthoFilter`, `wavelet` -""" -struct GLS{T<:WaveletBoundary} <: LSWavelet{T} - step ::Vector{LSStep{Float64}} # steps to be taken - norm1 ::Float64 # normalization of scaling coefs. - norm2 ::Float64 # normalization of detail coefs. - name ::String # name of scheme -end - -function GLS(w::WC, ::T=DEFAULT_BOUNDARY) where {WC<:WaveletClass, T<:WaveletBoundary} - name = WT.name(w) - schemedef = get(SCHEMES, name, nothing) - schemedef == nothing && throw(ArgumentError("scheme not found")) - return GLS{T}(schemedef[1], schemedef[2], schemedef[3], name) -end - -name(s::GLS) = s.name - -# TRANSFORM TYPE CONSTRUCTORS - -""" - wavelet(c[, t=WT.Filter][, boundary=WT.Periodic][, s=Real]) - -Construct wavelet type where `c` is a wavelet class, -`t` is the transformation type (`WT.Filter` or `WT.Lifting`), -and `boundary` is the type of boundary treatment. In the continuous case, s>0 -is the number of wavelets between the octaves ``2^J`` and ``2^{J+1}`` -(defaults to 8, which is most appropriate for music and other audio). - -# Examples -```julia -wavelet(WT.coif6) -wavelet(WT.db1, WT.Lifting) -wavelet(WT.morl,4) -``` - -**See also:** `WT.WaveletClass` -""" -function wavelet end - -wavelet(c::K, boundary::WaveletBoundary=DEFAULT_BOUNDARY) where {K<:Union{BiOrthoWaveletClass, OrthoWaveletClass}}= wavelet(c, Filter, boundary) -wavelet(c::OrthoWaveletClass, t::FilterTransform, boundary::WaveletBoundary=DEFAULT_BOUNDARY) = OrthoFilter(c, boundary) -wavelet(c::WaveletClass, t::LiftingTransform, boundary::WaveletBoundary=DEFAULT_BOUNDARY) = GLS(c, boundary) - - -# ------------------------------------------------------------ - -# Compute filters from the Daubechies class -# N is the number of zeros at -1 -function daubechies(N::Int) - @assert N > 0 - # Create polynomial - C = Vector{Int}(undef, N) - @inbounds for n = 0:N-1 - C[N-n] = binomial(N-1+n, n) - end - - # Find roots in y domain (truncated binomial series; (1 - y)^{-N}) - Y = roots(C) - - # Find roots in z domain: - # z + z^{-1} = 2 - 4*y - # where y is a root from above - Z = zeros(ComplexF64, 2*N-2) - @inbounds for i = 1:N-1 - Yi = Y[i] - d = 2*sqrt( Yi*Yi - Yi ) - y2 = 1 - 2*Yi - Z[i] = y2 + d - Z[i+N-1] = y2 -d - end - - # Retain roots inside unit circle - nr = 0 # count roots - @inbounds for i = eachindex(Z) - if abs(Z[i]) <= 1 + eps() - nr += 1 - end - end - - # Find coefficients of the polynomial - # (1 + z)^N * \prod_i (z - z_i) - R = Vector{ComplexF64}(undef, N + nr) - @inbounds for i = 1:N - R[i] = -1 - end - k = N - @inbounds for i = eachindex(Z) - if abs(Z[i]) <= 1 + eps() - k += 1 - R[k] = Z[i] - end - end - HH = vieta( R ) - - # Normalize coefficients - rmul!(HH, 1/norm(HH)) - return real(HH) -end - -# Compute roots of polynomial -# Input is a coefficient vector with highest powers first -function roots(C::AbstractVector) - A = compan(C) - return eigvals(A) -end - -# Create companion matrix for a polynomial -# Input is a coefficient vector with highest powers first -function compan(C::AbstractVector) - n = length(C) - A = zeros(n-1, n-1) - - if n > 1 - @inbounds A[1,:] .= -C[2:end] ./ C[1] - @inbounds A[2:n:end] .= 1 - end - return A -end - -# Vieta-like formula for computing polynomial coefficients from roots -# See -# http://www.mathworks.se/help/matlab/ref/poly.html -function vieta(R::AbstractVector) - n = length( R ) - C = zeros(ComplexF64, n+1) - C[1] = 1 - Ci::ComplexF64 = 0 - Cig::ComplexF64 = 0 - - @inbounds for k = 1:n - Ci = C[1] - for i = 1:k - Cig = C[i+1] - C[i+1] = Cig - R[k] * Ci - Ci = Cig - end - end - return C -end - - -# scaling filters h (low pass) -# the number at end of a filter name is the -# number of vanishing moments of the mother wavelet function -# sources: -# http://statweb.stanford.edu/~wavelab/ (Orthogonal/MakeONFilter.m) -# http://www.mathworks.com/matlabcentral/fileexchange/5502-filter-coefficients-to-popular-wavelets -### https://github.com/nigma/pywt/blob/master/src/wavelets_coeffs.template.h -# name => qmf -const FILTERS = Dict{String, Vector{Float64}}( -# Haar filter -"haar" => -[0.7071067811865475,0.7071067811865475] -, -# Daubechies filters, see daubechies() - -# Coiflet filters -"coif2" => -[-0.072732619513,0.337897662458,0.852572020212,0.384864846864,-0.072732619513,-0.015655728135] -, -"coif4" => -[0.0163873364635998,-0.0414649367819558,-0.0673725547222826,0.3861100668229939,0.8127236354493977,0.4170051844236707,-0.0764885990786692,-0.0594344186467388,0.0236801719464464,0.0056114348194211,-0.0018232088707116,-0.0007205494453679] -, -"coif6" => -[-0.0037935128644910,0.0077825964273254,0.0234526961418362,-0.0657719112818552,-0.0611233900026726,0.4051769024096150,0.7937772226256169,0.4284834763776168,-0.0717998216193117,-0.0823019271068856,0.0345550275730615,0.0158805448636158,-0.0090079761366615,-0.0025745176887502,0.0011175187708906,0.0004662169601129,-0.0000709833031381,-0.0000345997728362] -, -"coif8" => -[0.0008923136685824,-0.0016294920126020,-0.0073461663276432,0.0160689439647787,0.0266823001560570,-0.0812666996808907,-0.0560773133167630,0.4153084070304910,0.7822389309206135,0.4343860564915321,-0.0666274742634348,-0.0962204420340021,0.0393344271233433,0.0250822618448678,-0.0152117315279485,-0.0056582866866115,0.0037514361572790,0.0012665619292991,-0.0005890207562444,-0.0002599745524878,0.0000623390344610,0.0000312298758654,-0.0000032596802369,-0.0000017849850031] -, -"coif10" => -[-0.0002120808398259,0.0003585896879330,0.0021782363583355,-0.0041593587818186,-0.0101311175209033,0.0234081567882734,0.0281680289738655,-0.0919200105692549,-0.0520431631816557,0.4215662067346898,0.7742896037334738,0.4379916262173834,-0.0620359639693546,-0.1055742087143175,0.0412892087544753,0.0326835742705106,-0.0197617789446276,-0.0091642311634348,0.0067641854487565,0.0024333732129107,-0.0016628637021860,-0.0006381313431115,0.0003022595818445,0.0001405411497166,-0.0000413404322768,-0.0000213150268122,0.0000037346551755,0.0000020637618516,-0.0000001674428858,-0.0000000951765727] -, -# Symmlet filter -"sym4" => -[0.0455703458960000,-0.0178247014420000,-0.1403176241790000,0.4212345342040000,1.1366582434079999,0.7037390686560000,-0.0419109651250000,-0.1071489014180000] -, -"sym5" => -[0.0276321529580000,-0.0298424998690000,-0.2479513626130000,0.0234789231360000,0.8965816483800000,1.0230529668940000,0.2819906968540000,-0.0553441861170000,0.0417468644220000,0.0386547959550000] -, -"sym6" => -[-0.0110318675090000,0.0024999220930000,0.0632505626600000,-0.0297837512990000,-0.1027249698620000,0.4779043713330000,1.1138927839260000,0.6944579729580000,-0.0683231215870000,-0.1668632154120000,0.0049366123720000,0.0217847003270000] -, -"sym7" => -[0.0145213947620000,0.0056713426860000,-0.1524638718960000,-0.1980567068070000,0.4081839397250000,1.0857827098140000,0.7581626019640000,0.0246656594890000,-0.0700782912220000,0.0960147679360000,0.0431554525820000,-0.0178704316510000,-0.0014812259150000,0.0037926585340000] -, -"sym8" => -[-0.0047834585120000,-0.0007666908960000,0.0448236230420000,0.0107586117510000,-0.2026486552860000,-0.0866536154060000,0.6807453471900000,1.0991066305370001,0.5153986703740000,-0.0734625087610000,-0.0384935212630000,0.0694904659110000,0.0053863887540000,-0.0211456865280000,-0.0004283943000000,0.0026727933930000] -, -"sym9" => -[0.0019811937360000,0.0008765025390000,-0.0187693968360000,-0.0163033512260000,0.0427444336020000,0.0008251409290000,-0.0771721610970000,0.3376589236020000,1.0152597908320000,0.8730484073490000,0.0498828309590000,-0.2708937835030000,-0.0257864459300000,0.0877912515540000,0.0125288962420000,-0.0145155785530000,-0.0006691415090000,0.0015124873090000] -, -"sym10" => -[-0.0006495898960000,0.0000806612040000,0.0064957283750000,-0.0011375353140000,-0.0287862319260000,0.0081528167990000,0.0707035675500000,-0.0452407722180000,-0.0502565400920000,0.5428130112130000,1.0882515305000000,0.6670713381540000,-0.1002402150310000,-0.2255589722340000,0.0164188694260000,0.0649509245790000,-0.0020723639230000,-0.0122206426300000,0.0001352450200000,0.0010891704470000] -, -# Battle-Lemarie filter -"batt2" => -[-0.0000867523000000,-0.0001586010000000,0.0003617810000000,0.0006529220000000,-0.0015570100000000,-0.0027458800000000,0.0070644200000000,0.0120030000000000,-0.0367309000000000,-0.0488618000000000,0.2809310000000000,0.5781630000000000,0.2809310000000000,-0.0488618000000000,-0.0367309000000000,0.0120030000000000,0.0070644200000000,-0.0027458800000000,-0.0015570100000000,0.0006529220000000,0.0003617810000000,-0.0001586010000000,-0.0000867523000000] -, -"batt4" => -[0.0001033070000000,-0.0001642640000000,-0.0002018180000000,0.0003267490000000,0.0003959460000000,-0.0006556200000000,-0.0007804680000000,0.0013308600000000,0.0015462400000000,-0.0027452900000000,-0.0030786300000000,0.0057993200000000,0.0061414300000000,-0.0127154000000000,-0.0121455000000000,0.0297468000000000,0.0226846000000000,-0.0778079000000000,-0.0354980000000000,0.3068300000000000,0.5417360000000000,0.3068300000000000,-0.0354980000000000,-0.0778079000000000,0.0226846000000000,0.0297468000000000,-0.0121455000000000,-0.0127154000000000,0.0061414300000000,0.0057993200000000,-0.0030786300000000,-0.0027452900000000,0.0015462400000000,0.0013308600000000,-0.0007804680000000,-0.0006556200000000,0.0003959460000000,0.0003267490000000,-0.0002018180000000,-0.0001642640000000,0.0001033070000000] -, -"batt6" => -[0.0001011130000000,0.0001107090000000,-0.0001591680000000,-0.0001726850000000,0.0002514190000000,0.0002698420000000,-0.0003987590000000,-0.0004224850000000,0.0006355630000000,0.0006628360000000,-0.0010191200000000,-0.0010420700000000,0.0016465900000000,0.0016413200000000,-0.0026864600000000,-0.0025881600000000,0.0044400200000000,0.0040788200000000,-0.0074684800000000,-0.0063988600000000,0.0128754000000000,0.0099063500000000,-0.0229951000000000,-0.0148537000000000,0.0433544000000000,0.0208414000000000,-0.0914068000000000,-0.0261771000000000,0.3128690000000000,0.5283740000000000,0.3128690000000000,-0.0261771000000000,-0.0914068000000000,0.0208414000000000,0.0433544000000000,-0.0148537000000000,-0.0229951000000000,0.0099063500000000,0.0128754000000000,-0.0063988600000000,-0.0074684800000000,0.0040788200000000,0.0044400200000000,-0.0025881600000000,-0.0026864600000000,0.0016413200000000,0.0016465900000000,-0.0010420700000000,-0.0010191200000000,0.0006628360000000,0.0006355630000000,-0.0004224850000000,-0.0003987590000000,0.0002698420000000,0.0002514190000000,-0.0001726850000000,-0.0001591680000000,0.0001107090000000,0.0001011130000000] -, -# Beylkin filter -"beyl" => -[0.0993057653740000,0.4242153608130000,0.6998252140570000,0.4497182511490000,-0.1109275983480000,-0.2644972314460000,0.0269003088040000,0.1555387318770000,-0.0175207462670000,-0.0885436306230000,0.0196798660440000,0.0429163872740000,-0.0174604086960000,-0.0143658079690000,0.0100404118450000,0.0014842347820000,-0.0027360316260000,0.0006404853290000] -, -# Vaidyanathan filter -"vaid" => -[-0.0000629061180000,0.0003436319050000,-0.0004539566200000,-0.0009448971360000,0.0028438345470000,0.0007081375040000,-0.0088391034090000,0.0031538470560000,0.0196872150100000,-0.0148534480050000,-0.0354703986070000,0.0387426192930000,0.0558925236910000,-0.0777097509020000,-0.0839288843660000,0.1319716614170000,0.1350842271290000,-0.1944504717660000,-0.2634948024880000,0.2016121617750000,0.6356010598720000,0.5727977932110000,0.2501841295050000,0.0457993341110000] - - -) - - -# biortho filters -# name => (qmf1,qmf2) -const BIFILTERS = Dict{String, NTuple{2, Vector{Float64}}}( -# test -"test" => -([0.7071067811865475,0.7071067811865475], - [0.7071067811865475,0.7071067811865475]) - -) - - -# in matlab (d,p) -> (predict, update) -const SCHEMES = Dict{String,NTuple{3, Any}}( -# cdf 5/3 -> bior 2.2, cdf 9/7 -> bior 4.4 -# Cohen-Daubechies-Feauveau [Do Quan & Yo-Sung Ho. Optimized median lifting scheme for lossy image compression.] -"cdf9/7" => ([ LSStep(Update, [1.0,1.0]*1.5861343420604, 0), - LSStep(Predict, [1.0,1.0]*0.05298011857291494, 1), - LSStep(Update, [1.0,1.0]*-0.882911075531393, 0), - LSStep(Predict, [1.0,1.0]*-0.44350685204384654, 1)], - 1.1496043988603355, - 0.8698644516247099) -, - -# Haar, same as db1 -"haar" => ([LSStep(Predict, [-1.0], 0), - LSStep(Update, [0.5], 0)], - 0.7071067811865475, - 1.4142135623730951) -, -# Daubechies -"db1" => ([ LSStep(Predict, [-1.0], 0), - LSStep(Update, [0.5], 0)], - 0.7071067811865475, - 1.4142135623730951) -, -"db2" => ([ LSStep(Predict, [-1.7320508075688772], 0), - LSStep(Update, [-0.0669872981077807,0.4330127018922193], 1), - LSStep(Predict, [1.0], -1)], - 0.5176380902050414, - 1.9318516525781364) - -) - - - - -end From 54f528186939b960a8288a8795223311f9af2490 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thimot=C3=A9=20Dupuch?= Date: Sat, 16 Aug 2025 17:29:40 +0200 Subject: [PATCH 2/3] added underscore to filenames --- src/Plot/Plot.jl | 2 +- src/Plot/{plot.jl => _plot.jl} | 0 src/Threshold/Threshold.jl | 2 +- src/Threshold/{threshold.jl => _threshold.jl} | 0 src/Transforms/Transforms.jl | 2 +- src/Transforms/{transforms.jl => _transforms.jl} | 0 src/Util/Util.jl | 2 +- src/Util/{util.jl => _util.jl} | 0 src/WT/WT.jl | 2 +- src/WT/{wt.jl => _wt.jl} | 0 10 files changed, 5 insertions(+), 5 deletions(-) rename src/Plot/{plot.jl => _plot.jl} (100%) rename src/Threshold/{threshold.jl => _threshold.jl} (100%) rename src/Transforms/{transforms.jl => _transforms.jl} (100%) rename src/Util/{util.jl => _util.jl} (100%) rename src/WT/{wt.jl => _wt.jl} (100%) diff --git a/src/Plot/Plot.jl b/src/Plot/Plot.jl index cac8d37..a8dad91 100644 --- a/src/Plot/Plot.jl +++ b/src/Plot/Plot.jl @@ -1,6 +1,6 @@ module Plot -include("plot.jl") +include("_plot.jl") export wplotdots, wplotim diff --git a/src/Plot/plot.jl b/src/Plot/_plot.jl similarity index 100% rename from src/Plot/plot.jl rename to src/Plot/_plot.jl diff --git a/src/Threshold/Threshold.jl b/src/Threshold/Threshold.jl index c96d3fe..470cc91 100644 --- a/src/Threshold/Threshold.jl +++ b/src/Threshold/Threshold.jl @@ -4,7 +4,7 @@ using LinearAlgebra: rmul!, norm using Statistics: median! -include("threshold.jl") +include("_threshold.jl") include("basis_functions.jl") include("denoising.jl") include("entropy.jl") diff --git a/src/Threshold/threshold.jl b/src/Threshold/_threshold.jl similarity index 100% rename from src/Threshold/threshold.jl rename to src/Threshold/_threshold.jl diff --git a/src/Transforms/Transforms.jl b/src/Transforms/Transforms.jl index e17ca42..8a08e82 100644 --- a/src/Transforms/Transforms.jl +++ b/src/Transforms/Transforms.jl @@ -1,6 +1,6 @@ module Transforms -include("transforms.jl") +include("_transforms.jl") include("transforms_filter.jl") include("transforms_lifting.jl") include("transforms_maximal_overlap.jl") diff --git a/src/Transforms/transforms.jl b/src/Transforms/_transforms.jl similarity index 100% rename from src/Transforms/transforms.jl rename to src/Transforms/_transforms.jl diff --git a/src/Util/Util.jl b/src/Util/Util.jl index 881435c..ec3ac30 100644 --- a/src/Util/Util.jl +++ b/src/Util/Util.jl @@ -7,7 +7,7 @@ using DSP: conv include("non_dyadic.jl") include("dyadic.jl") -include("util.jl") +include("_util.jl") diff --git a/src/Util/util.jl b/src/Util/_util.jl similarity index 100% rename from src/Util/util.jl rename to src/Util/_util.jl diff --git a/src/WT/WT.jl b/src/WT/WT.jl index 697be8f..2ecdb24 100644 --- a/src/WT/WT.jl +++ b/src/WT/WT.jl @@ -1,6 +1,6 @@ module WT -include("wt.jl") +include("_wt.jl") export DiscreteWavelet, diff --git a/src/WT/wt.jl b/src/WT/_wt.jl similarity index 100% rename from src/WT/wt.jl rename to src/WT/_wt.jl From a6dd8f62e8b0af8f10c0fe5ec593e7ad757416cc Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Sun, 17 Aug 2025 09:42:44 +0800 Subject: [PATCH 3/3] fix tests --- src/Plot/_plot.jl | 2 ++ test/util.jl | 1 + 2 files changed, 3 insertions(+) diff --git a/src/Plot/_plot.jl b/src/Plot/_plot.jl index 3046a74..1a8632d 100644 --- a/src/Plot/_plot.jl +++ b/src/Plot/_plot.jl @@ -2,6 +2,8 @@ using ..Util using ..WT using ..Transforms +using LinearAlgebra: norm + # PLOTTING UTILITIES # return levels and detail coefficient centers on the interval [0,r) above (>=) threshold t diff --git a/test/util.jl b/test/util.jl index f35fd52..99a987a 100644 --- a/test/util.jl +++ b/test/util.jl @@ -1,3 +1,4 @@ +using Wavelets.Util @testset "Indexing" begin @test dyadicdetailn(0) == 1