diff --git a/src/Filters/filt.jl b/src/Filters/filt.jl
index c8c96a876..9d8a40955 100644
--- a/src/Filters/filt.jl
+++ b/src/Filters/filt.jl
@@ -522,7 +522,7 @@ function fftfilt(b::AbstractVector{T}, x::AbstractArray{T},
     # FFT of filter
     filterft = similar(tmp2)
     copyto!(tmp1, b)
-    tmp1[nb+1:end] = zero(T)
+    tmp1[nb+1:end] .= zero(T)
     mul!(filterft, p1, tmp1)
 
     # FFT of chunks
@@ -533,8 +533,8 @@ function fftfilt(b::AbstractVector{T}, x::AbstractArray{T},
             xstart = off - nb + npadbefore + 1
             n = min(nfft - npadbefore, nx - xstart + 1)
 
-            tmp1[1:npadbefore] = zero(T)
-            tmp1[npadbefore+n+1:end] = zero(T)
+            tmp1[1:npadbefore] .= zero(T)
+            tmp1[npadbefore+n+1:end] .= zero(T)
 
             copyto!(tmp1, npadbefore+1, x, colstart+xstart, n)
             mul!(tmp2, p1, tmp1)
diff --git a/src/estimation.jl b/src/estimation.jl
index 2d585a8cd..ff8746422 100644
--- a/src/estimation.jl
+++ b/src/estimation.jl
@@ -1,5 +1,6 @@
 module Estimation
-import Base: *
+
+using Compat.LinearAlgebra: eig, svd
 
 export esprit
 
@@ -8,10 +9,10 @@ export esprit
 
 ESPRIT [^Roy1986] algorithm for frequency estimation.
 Estimation of Signal Parameters via Rotational Invariance Techniques
-    
+
 Given length N signal "x" that is the sum of p sinusoids of unknown frequencies,
 estimate and return an array of the p frequencies.
-    
+
 # Arguments
 - `x::AbstractArray`: complex length N signal array
 - `M::Integer`: size of correlation matrix, must be <= N.
@@ -21,7 +22,7 @@ estimate and return an array of the p frequencies.
       For faster execution (due to smaller SVD), use small M or small N-M
 - `p::Integer`: number of sinusoids to estimate.
 - `Fs::Float64`: sampling frequency, in Hz.
- 
+
 # Returns
 length p real array of frequencies in units of Hz.
 
diff --git a/src/unwrap.jl b/src/unwrap.jl
index c4aea218f..67ca19541 100644
--- a/src/unwrap.jl
+++ b/src/unwrap.jl
@@ -1,80 +1,80 @@
 module Unwrap
 
-using Compat: selectdim
+import Compat: selectdim
+import Compat
 
 export unwrap, unwrap!
 
 """
-    unwrap!(m; range=2pi, wrap_around=tuple(fill(false, N)...), seed=-1)
+    unwrap!(m; dims=nothing, range=2pi)
 
-In-place version of unwrap(m; range, wrap_around, seed)
+In-place version of [`unwrap`](@ref).
 """
-unwrap!
+function unwrap!(m::AbstractArray{T,N}; dims=nothing, range=2T(pi), kwargs...) where {T, 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, kwargs...)
+end
+
+"""
+    unwrap!(y, m; dims=nothing, range=2pi)
+
+Unwrap `m` storing the result in `y`, see [`unwrap`](@ref).
+"""
+function unwrap!(y::AbstractArray{T,N}, m::AbstractArray{T,N}; dims=nothing, range::Number=2T(pi), kwargs...) where {T, 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)
+    elseif dims == 1:N
+        unwrap_nd!(y, m; range=range, kwargs...)
+    else
+        throw(ArgumentError("`unwrap!`: Invalid dims specified: $dims"))
+    end
+    return y
+end
+
+unwrap_kernel(range=2π) = (x, y) -> y - round((y - x) / range) * range
 
 """
-    unwrap(m::AbstractArray{T, N}; range=2pi, wrap_around=tuple(fill(false, N)...), seed=-1)
+    unwrap(m; dims=nothing, range=2pi)
 
-Assumes `m` to be an array of values of arbitrary dimension that
-have been wrapped to be inside the given range (centered around
-zero), and undoes the wrapping by identifying discontinuities. To operate on a
-subset of dimensions of an array, simply pass in the corresponding `view`.
+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].
 
+If `dims` is a `Range`, then the data is unwrapped according to
+the number of dimensions passed to `dims`.
+
 # Arguments
-- `A::AbstractArray{T, N}`: Array to unwrap
-- `range::Number`: Range of wrapped array
-- `wrap_around::NTuple{N, Bool}`:  When an element of this tuple is `true`, the
+- `m::AbstractArray{T, N}`: Array to unwrap
+- `range=2pi`: Range of wrapped array
+- `wrap_around=`:  When an element of this tuple is `true`, the
     unwrapping process will consider the edges along the corresponding axis
     of the image to be connected
-- `seed::Int`: Unwrapping of arrays with dimension > 1 uses a random initialization. This
+- `seed=-1`: Unwrapping of arrays with dimension > 1 uses a random initialization. This
     sets the seed of the RNG.
 """
-unwrap
-
-function unwrap(m::AbstractArray{T, N}, args...; kwargs...) where {T<:AbstractFloat, N}
-    unwrap!(copy(m), args...; kwargs...)
-end
-
-### Deprecated unwrap function
-function unwrap!(m::Array{T}, dim::Integer; range::Number=2pi) where T<:AbstractFloat
-    Base.depwarn(string("`unwrap!(m::Array{T}, dim::Integer; range::Number=2pi) ",
-        " where T<:AbstractFloat` is deprecated, ",
-        "use `unwrap!(m; range)` instead if your data has more than one dimension. ",
-        "If your data is one-dimensional but embedded in a larger-dimensional array, ",
-        "pass instead a view of each individual vector you would like to unwrap."),
-         :unwrap!)
-    thresh = range / 2
-    if size(m, dim) < 2
-        return m
-    end
-    for i = 2:size(m, dim)
-        d = selectdim(m, dim, i:i) - selectdim(m, dim, i-1:i-1)
-        slice_tuple = ntuple(n->(n==dim ? (i:i) : (1:size(m,n))), ndims(m))
-        offset = floor.((d.+thresh) / (range)) * range
-        m[slice_tuple...] = m[slice_tuple...] - offset
+function unwrap(m::AbstractArray{T,N}; dims=nothing, range=2T(pi), kwargs...) where {T, 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
-    return m
+    unwrap!(similar(m), m; dims=dims, range=range, kwargs...)
 end
 
-function unwrap!(A::AbstractVector{T};
-                 range::Number=2*convert(T, pi),
-                 wrap_around::NTuple{1, Bool}=(false,),
-                 seed::Int=-1) where T
-    @inbounds previous_element = A[1]
-    difference = zero(previous_element)
-    periods = 0
-    @inbounds for i in 2:length(A)
-        difference = A[i] - previous_element
-        periods += ifelse(difference >  range/2, -1, 0)
-        periods += ifelse(difference < -range/2,  1, 0)
-        previous_element = A[i]
-        A[i] = previous_element + one(previous_element) * periods * range
-    end
-    return A
-end
+@deprecate(unwrap(m::AbstractArray, dim::Integer; range::Number=2pi),
+    unwrap(m, dims=dim, range=range))
 
 #= Algorithm based off of
  M. A. Herráez, D. R. Burton, M. J. Lalor, and M. A. Gdeisat,
@@ -119,10 +119,11 @@ function Edge{N}(pixel_image::AbstractArray, ind1::CartesianIndex{N}, ind2::Cart
 end
 @inline Base.isless(e1::Edge, e2::Edge) = isless(e1.reliability, e2.reliability)
 
-function unwrap!(wrapped_image::AbstractArray{T, N};
-                 range::Number=2*convert(T, pi),
-                 wrap_around::NTuple{N, Bool}=tuple(fill(false, N)...),
-                 seed::Int=-1) where {T, N}
+function unwrap_nd!(dest::AbstractArray{T, N},
+                    src::AbstractArray{T, N};
+                    range::Number=2*convert(T, pi),
+                    wrap_around::NTuple{N, Bool}=tuple(fill(false, N)...),
+                    seed::Int=-1) where {T, N}
 
     if seed != -1
         srand(seed)
@@ -130,10 +131,10 @@ function unwrap!(wrapped_image::AbstractArray{T, N};
 
     range_T = convert(T, range)
 
-    pixel_image = init_pixels(wrapped_image)
+    pixel_image = init_pixels(src)
     calculate_reliability(pixel_image, wrap_around, range_T)
     edges = Edge{N}[]
-    num_edges = _predict_num_edges(size(wrapped_image), wrap_around)
+    num_edges = _predict_num_edges(size(src), wrap_around)
     sizehint!(edges, num_edges)
     for idx_dim=1:N
         populate_edges!(edges, pixel_image, idx_dim, wrap_around[idx_dim], range_T)
@@ -141,9 +142,9 @@ function unwrap!(wrapped_image::AbstractArray{T, N};
 
     sort!(edges, alg=MergeSort)
     gather_pixels!(pixel_image, edges)
-    unwrap_image!(wrapped_image, pixel_image, range_T)
+    unwrap_image!(dest, pixel_image, range_T)
 
-    return wrapped_image
+    return dest
 end
 
 function _predict_num_edges(size_img, wrap_around)
@@ -171,9 +172,9 @@ function gather_pixels!(pixel_image, edges)
     end
 end
 
-function unwrap_image!(image, pixel_image, range)
-    @Threads.threads for i in eachindex(image)
-        @inbounds image[i] = range * pixel_image[i].periods + pixel_image[i].val
+function unwrap_image!(dest, pixel_image, range)
+    @Threads.threads for i in eachindex(dest)
+        @inbounds dest[i] = range * pixel_image[i].periods + pixel_image[i].val
     end
 end
 
diff --git a/src/util.jl b/src/util.jl
index 64401c7c4..c5a9d0208 100644
--- a/src/util.jl
+++ b/src/util.jl
@@ -1,6 +1,7 @@
 module Util
 using ..DSP: @importffts, mul!
 import Base: *
+import Compat
 using Compat: copyto!, ComplexF64, selectdim, undef
 import Compat.LinearAlgebra.BLAS
 @importffts
diff --git a/test/filter_conversion.jl b/test/filter_conversion.jl
index 217a34b0a..ec6121263 100644
--- a/test/filter_conversion.jl
+++ b/test/filter_conversion.jl
@@ -80,8 +80,8 @@ using DSP, Compat, Compat.Test, FilterTestHelpers, Polynomials
 
     # And with only poles (no zeros)
     m_sos_only_poles = copy(m_sos_full)
-    m_sos_only_poles[:, 1:2] = 0
-    m_sos_only_poles[:, 3] = 1
+    m_sos_only_poles[:, 1:2] .= 0
+    m_sos_only_poles[:, 3] .= 1
     @test m_sos_only_poles ≈ sosfilter_to_matrix(convert(SecondOrderSections, ZeroPoleGain(Float64[], p, k)))
 
     # Test that a filter with repeated zeros is handled properly
diff --git a/test/lpc.jl b/test/lpc.jl
index 8c44ed331..cf53e40ff 100644
--- a/test/lpc.jl
+++ b/test/lpc.jl
@@ -1,8 +1,8 @@
 # Filter some noise, try to learn the coefficients
-coeffs = [1, .5, .3, .2]
-x = filt([1], coeffs, randn(50000))
-
 @testset "$method" for method in (LPCBurg(), LPCLevinson())
+    coeffs = [1, .5, .3, .2]
+    x = filt([1], coeffs, randn(50000))
+
     # Analyze the filtered noise, attempt to reconstruct the coefficients above
     ar, e = lpc(x[1000:end], length(coeffs)-1, method)
 
diff --git a/test/unwrap.jl b/test/unwrap.jl
index 75e49f282..ec06f87a5 100644
--- a/test/unwrap.jl
+++ b/test/unwrap.jl
@@ -6,18 +6,37 @@ using DSP, Compat, Compat.Test
     @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)
@@ -35,6 +54,7 @@ using DSP, Compat, Compat.Test
     end
 end
 
+# tests for multi-dimensional unwrapping
 @testset "Unwrap 2D" begin
     types = (Float32, Float64, BigFloat)
     for T in types
@@ -43,10 +63,10 @@ end
         A_unwrapped = v_unwrapped .+ v_unwrapped'
         A_wrapped = A_unwrapped .% (2convert(T, π))
 
-        test_unwrapped = unwrap(A_wrapped)
+        test_unwrapped = unwrap(A_wrapped, dims=1:2)
         d = first(A_unwrapped) - first(test_unwrapped)
         @test (test_unwrapped + d) ≈ A_unwrapped
-        unwrap!(A_wrapped)
+        unwrap!(A_wrapped, dims=1:2)
         d = first(A_unwrapped) - first(A_wrapped)
         @test (A_wrapped + d) ≈ A_unwrapped
 
@@ -55,7 +75,7 @@ end
         test_range = convert(T, 2)
         A_wrapped_range = A_unwrapped_range .% test_range
 
-        test_unwrapped_range = unwrap(A_wrapped_range; 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
 
@@ -68,13 +88,13 @@ end
         # make periodic
         wa_uw[end, :] = wa_uw[1, :]
         wa_w = wa_uw .% (2π)
-        wa_test = unwrap(wa_w, wrap_around=wrap_around, seed=0)
+        wa_test = unwrap(wa_w, dims=1:2, wrap_around=wrap_around, seed=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 wrap_around does not recover the borders
-        wa_test_nowa = unwrap(wa_w)
+        wa_test_nowa = unwrap(wa_w, dims=1:2)
         @test !(wa_test_nowa[end, :] ≈ wa_test_nowa[1, :])
 
     end
@@ -86,35 +106,35 @@ end
     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 = linspace(zero(T), convert(T, 2π), 11)
+        grid = linspace(zero(T), 2convert(T, π), 11)
         f_uw = f.(grid, grid', reshape(grid, 1, 1, :))
-        f_wr = f_uw .% (2π)
-        uw_test = unwrap(f_wr)
+        f_wr = f_uw .% (2convert(T, π))
+        uw_test = unwrap(f_wr, dims=1:3)
         offset = first(f_uw) - first(uw_test)
-        @test isapprox(uw_test + offset, f_uw, atol=1e-8)
+        @test (uw_test+offset) ≈ f_uw rtol=eps(T) #oop, nowrap
         # test in-place version
-        unwrap!(f_wr)
+        unwrap!(f_wr, dims=1:3)
         offset = first(f_uw) - first(f_wr)
-        @test isapprox(f_wr + offset, f_uw, atol=1e-8)
+        @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 .% (2π)
-        uw_test = unwrap(f_wr, wrap_around=(true, true, false))
+        f_wr = f_uw .% (2convert(T, π))
+        uw_test = unwrap(f_wr, dims=1:3)
         offset = first(f_uw) - first(uw_test)
-        @test isapprox(uw_test + offset, f_uw, atol=1e-8)
+        @test (uw_test+offset) ≈ f_uw #oop, 2wrap
         # test in-place version
-        unwrap!(f_wr, wrap_around=(true, true, false))
+        unwrap!(f_wr, dims=1:3, wrap_around=(true, true, false))
         offset = first(f_uw) - first(f_wr)
-        @test isapprox(f_wr + offset, f_uw, atol=1e-8)
+        @test (f_wr+offset) ≈ f_uw #ip, 2wrap
 
         f_uw = f_wraparound3.(grid, grid', reshape(grid, 1, 1, :))
-        f_wr = f_uw .% (2π)
-        uw_test = unwrap(f_wr, wrap_around=(true, true, true))
+        f_wr = f_uw .% (2convert(T, π))
+        uw_test = unwrap(f_wr, dims=1:3, wrap_around=(true, true, true))
         offset = first(f_uw) - first(uw_test)
-        @test isapprox(uw_test + offset, f_uw, atol=1e-8)
+        @test (uw_test+offset) ≈ f_uw #oop, 3wrap
         # test in-place version
-        unwrap!(f_wr, wrap_around=(true, true, true))
+        unwrap!(f_wr, dims=1:3, wrap_around=(true, true, true))
         offset = first(f_uw) - first(f_wr)
-        @test isapprox(f_wr + offset, f_uw, atol=1e-8)
+        @test (f_wr+offset) ≈ f_uw #oop, 3wrap
     end
 end