Skip to content

Commit

Permalink
add n-dimensional unwrap (#198)
Browse files Browse the repository at this point in the history
  • Loading branch information
platawiec authored and martinholters committed Jul 4, 2018
1 parent f92dcd7 commit 6e8901b
Show file tree
Hide file tree
Showing 8 changed files with 524 additions and 107 deletions.
3 changes: 2 additions & 1 deletion src/DSP.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
1 change: 1 addition & 0 deletions src/Filters/Filters.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
module Filters
using ..DSP: @importffts, mul!, rmul!
using ..Unwrap
using Polynomials, ..Util
import Base: *
using Compat: copyto!, undef
Expand Down
2 changes: 1 addition & 1 deletion src/Filters/response.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand Down
378 changes: 378 additions & 0 deletions src/unwrap.jl

Large diffs are not rendered by default.

65 changes: 1 addition & 64 deletions src/util.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,7 @@ using Compat: copyto!, ComplexF64, selectdim, undef
import Compat.LinearAlgebra.BLAS
@importffts

export unwrap!,
unwrap,
hilbert,
export hilbert,
Frequencies,
fftintype,
fftouttype,
Expand All @@ -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).
Expand Down
2 changes: 1 addition & 1 deletion test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down
140 changes: 140 additions & 0 deletions test/unwrap.jl
Original file line number Diff line number Diff line change
@@ -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
40 changes: 0 additions & 40 deletions test/util.jl
Original file line number Diff line number Diff line change
@@ -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)
Expand Down

0 comments on commit 6e8901b

Please sign in to comment.