diff --git a/src/Filters/stream_filt.jl b/src/Filters/stream_filt.jl index 9dee1d20..17cf0eec 100644 --- a/src/Filters/stream_filt.jl +++ b/src/Filters/stream_filt.jl @@ -623,31 +623,38 @@ function filt!( return bufIdx end -function filt(self::FIRFilter{Tk}, x::AbstractVector{Tx}) where {Th,Tx,Tk<:FIRKernel{Th}} - bufLen = outputlength(self, length(x)) +function filt(self::FIRFilter, x::AbstractVector) + buffer, bufLen = allocate_output(self, x) + samplesWritten = filt!(buffer, self, x) + return checked_ff_output(buffer, self, bufLen, samplesWritten) +end + +function allocate_output(sf::FIRFilter{Tk}, x::AbstractVector{Tx}) where {Th,Tx,Tk<:FIRKernel{Th}} # In some cases when `filt(::FIRFilter{FIRArbitrary}, x)` is called # with certain values of `x`, `filt!(buffer, ::FIRFilter{FIRArbitrary}, x)` # tries to write one sample too many to the buffer and a `BoundsError` - # is thrown. Add one extra sample to catch these exceptional cases. + # is thrown. Add one extra sample to catch these exceptional cases. # # See https://github.com/JuliaDSP/DSP.jl/issues/317 # # FIXME: Remove this if and when the code in # `filt!(buffer, ::FIRFilter{FIRArbitrary}, x)` # is updated to properly account for pathological arbitrary rates. + outLen = outputlength(sf, length(x)) if Tk <: FIRArbitrary - bufLen += 1 + outLen += 1 end - buffer = Vector{promote_type(Th,Tx)}(undef, bufLen) - samplesWritten = filt!(buffer, self, x) + out = Vector{promote_type(Th, Tx)}(undef, outLen) + return out, outLen +end +function checked_ff_output(out, ::FIRFilter{Tk}, bufLen, samplesWritten) where Tk<:FIRKernel if Tk <: FIRArbitrary - samplesWritten == bufLen || resize!(buffer, samplesWritten) + samplesWritten == bufLen || resize!(out, samplesWritten) else - @assert samplesWritten == bufLen + samplesWritten == bufLen || throw(AssertionError("Length of resampled output different from expectation.")) end - - return buffer + return out end @@ -689,24 +696,35 @@ function resample(x::AbstractVector, rate::AbstractFloat, h::Vector, Nϕ::Intege _resample!(x, rate, FIRFilter(h, rate, Nϕ)) end -function _resample!(x::AbstractVector, rate::Real, self::FIRFilter) +function _resample!(x::AbstractVector{T}, rate::Real, sf::FIRFilter) where T + undelay!(sf) + outLen = ceil(Int, length(x) * rate) + xPadded = copyto!(allocate_inbuf(sf, length(x), T, outLen), x) + + y = filt(sf, xPadded) + return checked_resample_output(y, outLen) +end + +function undelay!(sf::FIRFilter) # Get delay, in # of samples at the output rate, caused by filtering processes - τ = timedelay(self) + τ = timedelay(sf) # Use setphase! to # a) adjust the input samples to skip over before producing and output (integer part of τ) # b) set the ϕ index of the PFB (fractional part of τ) - setphase!(self, τ) + setphase!(sf, τ) +end - # Calculate the number of 0's required - outLen = ceil(Int, length(x) * rate) - reqInlen = inputlength(self, outLen, RoundUp) - reqZerosLen = reqInlen - length(x) - xPadded = [x; zeros(eltype(x), reqZerosLen)] +function allocate_inbuf(sf::FIRFilter, xLen::Int, ::Type{T}, outLen::Int) where T<:Number + reqInlen = inputlength(sf, outLen, RoundUp) + xPadded = Vector{T}(undef, reqInlen) + xPadded[xLen+1:end] .= zero(T) + return xPadded +end - y = filt(self, xPadded) - @assert length(y) >= outLen - length(y) > outLen && resize!(y, outLen) +function checked_resample_output(y::AbstractVector, outLen::Int) + length(y) >= outLen || throw(AssertionError("Resample output shorter than expected.")) + resize!(y, outLen) return y end @@ -742,11 +760,21 @@ end resample(x::AbstractArray, rate::Real, args::Real...; dims) = _resample!(x, rate, FIRFilter(rate, args...); dims) -_resample!(x::AbstractArray, rate::Real, sf::FIRFilter; dims) = +function _resample!(x::AbstractArray{T}, rate::Real, sf::FIRFilter; dims::Int) where T + undelay!(sf) + size_v = size(x, dims) + outLen = ceil(Int, size_v * rate) + xPadded = allocate_inbuf(sf, size_v, T, outLen) + buffer, bufLen = allocate_output(sf, xPadded) + mapslices(x; dims) do v - reset!(sf) - _resample!(v, rate, sf) + undelay!(reset!(sf)) + copyto!(xPadded, v) + samplesWritten = filt!(buffer, sf, xPadded) + checked_ff_output(buffer, sf, bufLen, samplesWritten) + return checked_resample_output(buffer, outLen) end +end # # References diff --git a/test/resample.jl b/test/resample.jl index 4e06009d..4d3683a9 100644 --- a/test/resample.jl +++ b/test/resample.jl @@ -31,8 +31,8 @@ end @test resample(mat, rate; dims=1) == resample(mat, rate, h; dims=1) @test resample(mat, rate; dims=2) == resample(mat, rate, h; dims=2) - rate = 1//2 - x_ml = reference_vec("resample_x.txt") + rate = 1//2 + x_ml = reference_vec("resample_x.txt") h_ml = reference_vec("resample_taps_1_2.txt") y_ml = reference_vec("resample_y_1_2.txt") expected_result = [y_ml ℯ * y_ml] @@ -54,12 +54,12 @@ end expected_result_3d = permutedims(reshape(expected_result, (size(expected_result, 1), size(expected_result, 2), 1)), (3, 1, 2)) X_3d = permutedims(reshape(X, (size(X, 1), size(X, 2), 1)), (3, 1, 2)) - y_jl = resample(X_3d, rate, h_ml, dims=2) + y_jl = resample(X_3d, rate, h_ml; dims=2) @test y_jl ≈ expected_result_3d expected_result_3d = permutedims(expected_result_3d, (1, 3, 2)) X_3d = permutedims(X_3d, (1, 3, 2)) - y_jl = resample(X_3d, rate, h_ml, dims=3) + y_jl = resample(X_3d, rate, h_ml; dims=3) @test y_jl ≈ expected_result_3d end