Merge a14484f into 2cce19f

JuliaDSP · Nov 10, 2014 · 4e4e289 · 4e4e289
2 parents 2cce19f + a14484f
commit 4e4e289
Show file tree

Hide file tree

Showing 5 changed files with 129 additions and 9 deletions.
diff --git a/README.md b/README.md
@@ -94,13 +94,17 @@ wavelet(name::String; transform::String="filter", boundary::String="per")
 waveletfilter(name::String; boundary::String="per")
 waveletls(name::String; boundary::String="per")
 # DWT (discrete wavelet transform)
-dwt(x::AbstractArray, wt::DiscreteWavelet, L::Integer=nscales(size(x,1)))
-idwt(x::AbstractArray, wt::DiscreteWavelet, L::Integer=nscales(size(x,1)))
+dwt(x::AbstractArray, wt::DiscreteWavelet, L::Integer=nscales(x))
+idwt(x::AbstractArray, wt::DiscreteWavelet, L::Integer=nscales(x))
 dwt!(y::AbstractArray, x::AbstractArray, filter::OrthoFilter, L::Integer, fw::Bool)
 dwt!(y::AbstractArray, scheme::GLS, L::Integer, fw::Bool)
 # DWTC (column-wise discrete wavelet transform)
-dwtc(x::AbstractArray, wt::DiscreteWavelet, L::Integer=nscales(size(x,1)))
-idwtc(x::AbstractArray, wt::DiscreteWavelet, L::Integer=nscales(size(x,1)))
+dwtc(x::AbstractArray, wt::DiscreteWavelet, L::Integer=nscales(x))
+idwtc(x::AbstractArray, wt::DiscreteWavelet, L::Integer=nscales(x))
+# WPT (wavelet packet transform)
+wpt(x::AbstractArray, wt::DiscreteWavelet, L::Integer=nscales(x))
+iwpt(x::AbstractArray, wt::DiscreteWavelet, L::Integer=nscales(x))
+wpt!(y::AbstractArray, x::AbstractArray, filter::OrthoFilter, L::Integer, fw::Bool)
 ```
 
 #### Wavelet class information

diff --git a/src/transforms.jl b/src/transforms.jl
@@ -1,6 +1,8 @@
 module Transforms
 using ..Util, ..WaveletTypes
-export dwt, idwt, dwt!, dwtc, idwtc
+export  dwt, idwt, dwt!, 
+        dwtc, idwtc, wpt, 
+        iwpt, wpt!
 
 # TODO Use StridedArray instead of AbstractArray where writing to array.
 
@@ -34,7 +36,7 @@ export dwt, idwt, dwt!, dwtc, idwtc
 ## GENERAL TRANSFORM FUNCTIONS
 
 # general functions for all types and dimensions
-for Xwt in (:dwt, :idwt, :dwtc, :idwtc)
+for Xwt in (:dwt, :idwt, :dwtc, :idwtc, :wpt, :iwpt)
 @eval begin
     # assume full transform
     $Xwt(x::AbstractArray, wt::DiscreteWavelet) = $Xwt(x, wt, nscales(size(x,1)))
@@ -44,7 +46,7 @@ end
 end
 
 # 1-D, 2-D MRA
-# methods with Array allocation
+# DWT methods with Array allocation
 for (Xwt, fw) in ((:dwt, true), (:idwt, false))
 @eval begin
     function $Xwt{T<:FloatingPoint}(x::AbstractArray{T}, wt::FilterWavelet, L::Integer)
@@ -59,6 +61,17 @@ for (Xwt, fw) in ((:dwt, true), (:idwt, false))
     end
 end
 end
+# 1-D
+# WPT methods with Array allocation
+for (Xwt, fw) in ((:wpt, true), (:iwpt, false))
+@eval begin
+    function $Xwt{T<:FloatingPoint}(x::AbstractArray{T}, wt::FilterWavelet, L::Integer)
+        y = Array(T, size(x))
+        wpt!(y, x, wt, L, $fw)
+        return y
+    end
+end
+end
 
 # column-wise transforms or color images, transform each x[:,...,:,i] separately
 for (Xwtc, Xwt) in ((:dwtc, :dwt), (:idwtc, :idwt))

diff --git a/src/transforms_filter.jl b/src/transforms_filter.jl
@@ -22,7 +22,7 @@ function makereverseqmf(h::AbstractVector, fw::Bool, T::Type=eltype(h))
     return scfilter, dcfilter
 end
 
-# 1-d
+# 1-D
 # writes to y
 function dwt!{T<:FloatingPoint}(y::AbstractVector{T}, x::AbstractVector{T}, filter::OrthoFilter, L::Integer, fw::Bool)
     si = Array(T, length(filter)-1)       # tmp filter vector
@@ -94,7 +94,7 @@ function unsafe_dwt1level!{T<:FloatingPoint}(y::AbstractVector{T}, x::AbstractVe
     return y
 end
 
-# 2-d
+# 2-D
 # writes to y
 function dwt!{T<:FloatingPoint}(y::Matrix{T}, x::AbstractMatrix{T}, filter::OrthoFilter, L::Integer, fw::Bool)
     n = size(x,1)
@@ -181,6 +181,80 @@ function dwt!{T<:FloatingPoint}(y::Matrix{T}, x::AbstractMatrix{T}, filter::Orth
     return y
 end
 
+# WPT
+# 1-D
+# writes to y
+function wpt!{T<:FloatingPoint}(y::AbstractVector{T}, x::AbstractVector{T}, filter::OrthoFilter, L::Integer, fw::Bool)
+    si = Array(T, length(filter)-1)       # tmp filter vector
+    if L > 1
+        if fw
+            ns = length(x)>>1
+        else
+            ns = length(x)
+        end
+    else
+        ns = 0
+    end
+    snew = Array(T, ns)
+    scfilter, dcfilter = makereverseqmf(filter.qmf, fw, T)
+
+    wpt!(y, x, filter, L, fw, dcfilter, scfilter, si, snew)
+    return y
+end
+function wpt!{T<:FloatingPoint}(y::AbstractVector{T}, x::AbstractVector{T}, filter::OrthoFilter, L::Integer, fw::Bool, dcfilter::Vector{T}, scfilter::Vector{T}, si::Vector{T}, snew::Vector{T})
+    n = length(x)
+    J = nscales(n)
+    @assert size(x) == size(y)
+    @assert isdyadic(y)
+    @assert 0 <= L <= J
+    is(y,x) && error("input vector is output vector")
+    L == 0 && return copy!(y,x) 
+
+    if fw
+        unsafe_dwt1level!(y, x, filter, fw, dcfilter, scfilter, si)
+        if L > 1  # recursion
+            @assert length(snew) >= n>>1
+            nj = detailn(J-1)
+            dx = unsafe_vectorslice(snew, 1, nj)
+            dy = unsafe_vectorslice(y, detailindex(J-1,1), nj)
+            copy!(dx, dy)
+            # detail
+            wpt!(dy, dx, filter, L-1, fw, dcfilter, scfilter, si, snew)
+            dy = unsafe_vectorslice(y, 1, nj)
+            copy!(dx, dy)
+            # scaling
+            wpt!(dy, dx, filter, L-1, fw, dcfilter, scfilter, si, snew)
+        end
+    else
+        if L == 1
+            unsafe_dwt1level!(y, x, filter, fw, dcfilter, scfilter, si)
+        else
+            @assert length(snew) >= n
+            first = true
+            while L > 0
+                ix = 1
+                nj = detailn(tl2level(n, L-1))
+                dx = unsafe_vectorslice(snew, 1, nj)
+                while ix <= n
+                    dy = unsafe_vectorslice(y, ix, nj)
+                    if first
+                        dx = unsafe_vectorslice(x, ix, nj)
+                    else
+                        copy!(dx, dy)
+                    end
+                    unsafe_dwt1level!(dy, dx, filter, fw, dcfilter, scfilter, si)
+                    ix += nj
+                end
+                L -= 1
+                first = false
+            end
+        end
+    end
+
+    return y
+end
+
+
 macro filtermainloop(si, silen, b, val)
     quote
         @inbounds for j=2:$(esc(silen))

diff --git a/src/util.jl b/src/util.jl
@@ -17,6 +17,7 @@ detailn(j::Integer) = 2^j
 # number of scales of dyadic length signal (n=2^J)
 nscales(n::Integer) = int(log2(n))
 nscales(x::Vector) = nscales(length(x))
+nscales(x::Array) = nscales(size(x,1))
 # the largest detail level
 maxlevel(n::Integer) = nscales(n)-1
 maxlevel(x::Vector) = maxlevel(length(x))

diff --git a/test/transforms.jl b/test/transforms.jl
@@ -168,6 +168,34 @@ EE = Exception
 @test_throws EE dwt(randn(4,4),uwt)
 
 
+# ============= WPT ================
+
+wf = wavelet("db2")
+x = randn(16)
+
+L = 1
+wp = wpt(x,wf,L)
+dw = dwt(x,wf,L)
+@test_approx_eq wp dw
+@test_approx_eq iwpt(wp,wf,L) x
+
+L = 2
+wp = wpt(x,wf,L)
+dw = dwt(x,wf,L)
+dw2 = copy(dw)
+dw2[9:end] = dwt(dw[9:end],wf,1)
+@test_approx_eq dw[1:8] wp[1:8]
+@test_approx_eq dw2 wp
+@test_approx_eq iwpt(wp,wf,L) x
+
+L = 3
+wp = wpt(x,wf,L)
+dw = dwt(x,wf,L)
+@test_approx_eq dw[1:4] wp[1:4]
+@test_approx_eq dwt(dw2[5:8],wf,1) wp[5:8]
+@test_approx_eq dwt(dw2[9:12],wf,1) wp[9:12]
+@test_approx_eq dwt(dw2[13:16],wf,1) wp[13:16]
+@test_approx_eq iwpt(wp,wf,L) x
 
 # ============= tranform low level functions ================