add n-dimensional unwrap

src/DSP.jl

@@ -46,14 +46,15 @@ if VERSION >= v"0.7.0-DEV.602"
 end
 
 include("util.jl")
+include("unwrap.jl")
 include("windows.jl")
 include("periodograms.jl")
 include("Filters/Filters.jl")
 include("lpc.jl")
 include("estimation.jl")
 
 using Reexport
-@reexport using .Util, .Windows, .Periodograms, .Filters, .LPC, .Estimation
+@reexport using .Util, .Windows, .Periodograms, .Filters, .LPC, .Unwrap, .Estimation
 
 include("deprecated.jl")
 end

src/Filters/Filters.jl

@@ -1,5 +1,6 @@
 module Filters
 using ..DSP: @importffts, mul!, rmul!
+using ..Unwrap
 using Polynomials, ..Util
 import Base: *
 using Compat: copyto!, undef

src/Filters/response.jl

@@ -54,7 +54,7 @@ or frequencies `w` in radians/sample.
 """
 function phasez(filter::FilterCoefficients, w = Compat.range(0, stop=π, length=250))
     h = freqz(filter, w)
-    unwrap(-atan2.(imag(h), real(h)))
+    unwrap(-atan2.(imag(h), real(h)); dims=ndims(h))
 end
 
 

src/util.jl

@@ -6,9 +6,7 @@ using Compat: copyto!, ComplexF64, selectdim, undef
 import Compat.LinearAlgebra.BLAS
 @importffts
 
-export  unwrap!,
-        unwrap,
-        hilbert,
+export  hilbert,
         Frequencies,
         fftintype,
         fftouttype,
@@ -29,67 +27,6 @@ export  unwrap!,
         polyfit,
         shiftin!
 
-"""
-    unwrap!(m; dims=nothing, range=2pi)
-
-In-place version of [`unwrap`](@ref).
-"""
-function unwrap!(m::AbstractArray{<:Any,N}; dims=nothing, range=2pi) where N
-    if dims === nothing && N != 1
-        Base.depwarn("`unwrap!(m::AbstractArray)` is deprecated, use `unwrap!(m, dims=ndims(m))` instead", :unwrap!)
-        dims = N
-    end
-    unwrap!(m, m; dims=dims, range=range)
-end
-
-"""
-    unwrap!(y, m; dims=nothing, range=2pi)
-
-Unwrap `m` storing the result in `y`, see [`unwrap`](@ref).
-"""
-function unwrap!(y::AbstractArray{<:Any,N}, m::AbstractArray{<:Any,N}; dims=nothing, range::Number=2pi) where {N}
-    if dims === nothing
-        if N != 1
-            throw(ArgumentError("`unwrap!`: required keyword parameter dims missing"))
-        end
-        dims = 1
-    end
-    if dims isa Integer
-        Compat.accumulate!(unwrap_kernel(range), y, m, dims=dims)
-    else
-        # this is where #198 can hook in
-        throw(ArgumentError("`unwrap!`: Invalid dims specified: $dims"))
-    end
-    return y
-end
-
-unwrap_kernel(range=2π) = (x, y) -> y - round((y - x) / range) * range
-
-@deprecate(unwrap!(m::AbstractArray, dim::Integer; range::Number=2pi),
-    unwrap!(m, dims=dim, range=range))
-
-
-"""
-    unwrap(m; dims=nothing, range=2pi)
-
-Assumes m (along the single dimension dims) to be sequences of values that
-have been wrapped to be inside the given range (centered around
-zero), and undoes the wrapping by identifying discontinuities.
-
-A common usage is for a phase measurement over time, such as when
-comparing successive frames of a short-time-fourier-transform, as
-each frame is wrapped to stay within (-pi, pi].
-"""
-function unwrap(m::AbstractArray{<:Any,N}; dims=nothing, range=2pi) where {N}
-    if dims === nothing && N != 1
-        Base.depwarn("`unwrap(m::AbstractArray)` is deprecated, use `unwrap(m, dims=ndims(m))` instead", :unwrap)
-        dims = ndims(m)
-    end
-    unwrap!(similar(m), m; dims=dims, range=range)
-end
-
-@deprecate(unwrap(m::AbstractArray, dim::Integer; range::Number=2pi),
-    unwrap(m, dims=dim, range=range))
 
 function hilbert(x::StridedVector{T}) where T<:FFTW.fftwReal
 # Return the Hilbert transform of x (a real signal).

test/runtests.jl

@@ -2,7 +2,7 @@ using Compat, DSP, AbstractFFTs, FFTW, Compat.Test
 
 testfiles = [ "dsp.jl", "util.jl", "windows.jl", "filter_conversion.jl",
     "filter_design.jl", "filter_response.jl", "filt.jl", "filt_stream.jl",
-    "periodograms.jl", "resample.jl", "lpc.jl", "estimation.jl"]
+    "periodograms.jl", "resample.jl", "lpc.jl", "estimation.jl", "unwrap.jl"]
 
 for testfile in testfiles
     eval(@__MODULE__, :(@testset $testfile begin include($testfile) end))

test/unwrap.jl

@@ -0,0 +1,140 @@
+using DSP, Compat, Compat.Test
+
+@testset "Unwrap 1D" begin
+    @test unwrap([0.1, 0.2, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4]
+    @test unwrap([0.1, 0.2 + 2pi, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4]
+    @test unwrap([0.1, 0.2 - 2pi, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4]
+    @test unwrap([0.1, 0.2 - 2pi, 0.3 - 2pi, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4]
+    @test unwrap([0.1 + 2pi, 0.2, 0.3, 0.4]) ≈ [0.1 + 2pi, 0.2 + 2pi, 0.3 + 2pi, 0.4 + 2pi]
+    @test unwrap([0.1, 0.2 + 6pi, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4]
+
+    test_v = [0.1, 0.2, 0.3 + 2pi, 0.4]
+    res_v = unwrap(test_v)
+    @test test_v ≈ [0.1, 0.2, 0.3 + 2pi, 0.4]
+    res_v .= 0
+    unwrap!(res_v, test_v)
+    @test res_v ≈ [0.1, 0.2, 0.3, 0.4]
+    @test test_v ≈ [0.1, 0.2, 0.3 + 2pi, 0.4]
+    unwrap!(test_v)
+    @test test_v ≈ [0.1, 0.2, 0.3, 0.4]
+
+    # test unwrapping within multi-dimensional array
+    wrapped = [0.1, 0.2 + 2pi, 0.3, 0.4]
+    unwrapped = [0.1, 0.2, 0.3, 0.4]
+    wrapped = hcat(wrapped, wrapped)
+    unwrapped = hcat(unwrapped, unwrapped)
+    @test unwrap(wrapped, dims=2) ≈ wrapped
+    @test unwrap(wrapped, dims=1) ≈ unwrapped
+    @test unwrap!(copy(wrapped), dims=2) ≈ wrapped
+    @test unwrap!(copy(wrapped), dims=1) ≈ unwrapped
+
+    # this should eventually default to the multi-dimensional case
+    @test_throws ArgumentError unwrap!(similar(wrapped), wrapped)
+
+    # test unwrapping with other ranges
+    unwrapped = [1.0:100;]
+    wrapped = Float64[i % 10 for i in unwrapped]
+    @test unwrap(wrapped, range=10) ≈ unwrapped
+
+    # test generically typed unwrapping
+    types = (Float32, Float64, BigFloat)
+    for T in types
+        srand(1234)
+        A_unwrapped = collect(Compat.range(0, stop=4convert(T, π), length=10))
+        A_wrapped = A_unwrapped .% (2convert(T, π))
+
+        @test unwrap(A_wrapped) ≈ A_unwrapped
+        unwrap!(A_wrapped)
+        @test A_wrapped ≈ A_unwrapped
+
+        A_unwrapped_range = collect(Compat.range(0, stop=4, length=10))
+        test_range = convert(T, 2)
+        A_wrapped_range = A_unwrapped_range .% test_range
+        @test unwrap(A_wrapped_range; range=test_range) ≈ A_unwrapped_range
+    end
+end
+
+# tests for multi-dimensional unwrapping
+@testset "Unwrap 2D" begin
+    types = (Float32, Float64, BigFloat)
+    for T in types
+        srand(1234)
+        v_unwrapped = collect(Compat.range(0, stop=4convert(T, π), length=7))
+        A_unwrapped = v_unwrapped .+ v_unwrapped'
+        A_wrapped = A_unwrapped .% (2convert(T, π))
+
+        test_unwrapped = unwrap(A_wrapped, dims=1:2)
+        d = first(A_unwrapped) - first(test_unwrapped)
+        @test (test_unwrapped .+ d) ≈ A_unwrapped
+        unwrap!(A_wrapped, dims=1:2)
+        d = first(A_unwrapped) - first(A_wrapped)
+        @test (A_wrapped .+ d) ≈ A_unwrapped
+
+        v_unwrapped_range = collect(Compat.range(0, stop=4, length=7))
+        A_unwrapped_range = v_unwrapped_range .+ v_unwrapped_range'
+        test_range = convert(T, 2)
+        A_wrapped_range = A_unwrapped_range .% test_range
+
+        test_unwrapped_range = unwrap(A_wrapped_range, dims=1:2; range=test_range)
+        d = first(A_unwrapped_range) - first(test_unwrapped_range)
+        @test (test_unwrapped_range .+ d) ≈ A_unwrapped_range
+
+        # Test circular_dims
+        # after unwrapping, pixels at borders should be equal to corresponding pixels
+        # on other side
+        circular_dims = (true, true)
+        wa_vec = Compat.range(0, stop=4convert(T, π), length=10)
+        wa_uw = wa_vec .+ zeros(10)'
+        # make periodic
+        wa_uw[end, :] = wa_uw[1, :]
+        wa_w = wa_uw .% (2π)
+        wa_test = unwrap(wa_w, dims=1:2, circular_dims=circular_dims, rng=MersenneTwister(0))
+        # with wrap-around, the borders should be equal, but for this problem the
+        # image may not be recovered exactly
+        @test wa_test[:, 1] ≈ wa_test[:, end]
+        @test wa_test[end, :] ≈ wa_test[1, :]
+        # In this case, calling unwrap w/o circular_dims does not recover the borders
+        wa_test_nowa = unwrap(wa_w, dims=1:2)
+        @test !(wa_test_nowa[end, :] ≈ wa_test_nowa[1, :])
+
+    end
+end
+
+@testset "Unwrap 3D" begin
+    types = (Float32, Float64, BigFloat)
+    f(x, y, z) = 0.1x^2 - 2y + 2z
+    f_wraparound2(x, y, z) = 5*sin(x) + 2*cos(y) + z
+    f_wraparound3(x, y, z) = 5*sin(x) + 2*cos(y) - 4*cos(z)
+    for T in types
+        grid = Compat.range(zero(T), stop=2convert(T, π), length=11)
+        f_uw = f.(grid, grid', reshape(grid, 1, 1, :))
+        f_wr = f_uw .% (2convert(T, π))
+        uw_test = unwrap(f_wr, dims=1:3)
+        offset = first(f_uw) - first(uw_test)
+        @test (uw_test.+offset) ≈ f_uw rtol=eps(T) #oop, nowrap
+        # test in-place version
+        unwrap!(f_wr, dims=1:3)
+        offset = first(f_uw) - first(f_wr)
+        @test (f_wr.+offset) ≈ f_uw rtol=eps(T) #ip, nowrap
+
+        f_uw = f_wraparound2.(grid, grid', reshape(grid, 1, 1, :))
+        f_wr = f_uw .% (2convert(T, π))
+        uw_test = unwrap(f_wr, dims=1:3)
+        offset = first(f_uw) - first(uw_test)
+        @test (uw_test.+offset) ≈ f_uw #oop, 2wrap
+        # test in-place version
+        unwrap!(f_wr, dims=1:3, circular_dims=(true, true, false))
+        offset = first(f_uw) - first(f_wr)
+        @test (f_wr.+offset) ≈ f_uw #ip, 2wrap
+
+        f_uw = f_wraparound3.(grid, grid', reshape(grid, 1, 1, :))
+        f_wr = f_uw .% (2convert(T, π))
+        uw_test = unwrap(f_wr, dims=1:3, circular_dims=(true, true, true))
+        offset = first(f_uw) - first(uw_test)
+        @test (uw_test.+offset) ≈ f_uw #oop, 3wrap
+        # test in-place version
+        unwrap!(f_wr, dims=1:3, circular_dims=(true, true, true))
+        offset = first(f_uw) - first(f_wr)
+        @test (f_wr.+offset) ≈ f_uw #oop, 3wrap
+    end
+end

test/util.jl

@@ -1,45 +1,5 @@
 using DSP, Compat, Compat.Test
 
-@testset "unwrap" begin
-    @test unwrap([0.1, 0.2, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4]
-    @test unwrap([0.1, 0.2 + 2pi, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4]
-    @test unwrap([0.1, 0.2 - 2pi, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4]
-    @test unwrap([0.1, 0.2 - 2pi, 0.3 - 2pi, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4]
-    @test unwrap([0.1 + 2pi, 0.2, 0.3, 0.4]) ≈ [0.1 + 2pi, 0.2 + 2pi, 0.3 + 2pi, 0.4 + 2pi]
-    @test unwrap([0.1, 0.2 + 6pi, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4]
-
-    test_v = [0.1, 0.2, 0.3 + 2pi, 0.4]
-    res_v = unwrap(test_v)
-    @test test_v ≈ [0.1, 0.2, 0.3 + 2pi, 0.4]
-    res_v .= 0
-    unwrap!(res_v, test_v)
-    @test res_v ≈ [0.1, 0.2, 0.3, 0.4]
-    @test test_v ≈ [0.1, 0.2, 0.3 + 2pi, 0.4]
-    unwrap!(test_v)
-    @test test_v ≈ [0.1, 0.2, 0.3, 0.4]
-
-    # test multi-dimensional unwrapping
-    wrapped = [0.1, 0.2 + 2pi, 0.3, 0.4]
-    unwrapped = [0.1, 0.2, 0.3, 0.4]
-    wrapped = hcat(wrapped, wrapped)
-    unwrapped = hcat(unwrapped, unwrapped)
-    @test unwrap(wrapped, dims=2) ≈ wrapped
-    @test unwrap(wrapped, dims=1) ≈ unwrapped
-    @test unwrap!(copy(wrapped), dims=2) ≈ wrapped
-    @test unwrap!(copy(wrapped), dims=1) ≈ unwrapped
-
-    # these should be implemented as part of #198
-    @test_throws ArgumentError unwrap(wrapped, dims=1:2)
-    @test_throws ArgumentError unwrap!(similar(wrapped), wrapped, dims=1:2)
-    # this should eventually default to the above
-    @test_throws ArgumentError unwrap!(similar(wrapped), wrapped)
-
-    # test unwrapping with other ranges
-    unwrapped = [1.0:100;]
-    wrapped = Float64[i % 10 for i in unwrapped]
-    @test unwrap(wrapped, range=10) ≈ unwrapped
-end
-
 @testset "hilbert" begin
     # Testing hilbert transform
     t = (0:1/256:2-1/256)