Skip to content

Commit

Permalink
Fix output length of resample (#596)
Browse files Browse the repository at this point in the history 
* Fix `outputlength(::FIRRational, inputlength)`

* Improve `outputlength(::FIRArbitrary, inputlength)`

This may still give results that are inconsistent with how many samples `filt!`
produces as the latter accumulates time-deltas (modulo Nϕ) which may be
numerically different. However, `outputlength(H, length(x)) == filt!(y, H, x)`
should hold more often this way.

Also adjust `inputlength` accordingly.

* Support specifying rounding mode in `inputlength`

* Let `resample` return expected length

Fixes #186.

* simplify phase accumulator

modf

* cosmetics, `DomainError`

---------

Co-authored-by: wheeheee <[email protected]>
  • Loading branch information
@martinholters @wheeheee
martinholters and wheeheee authored Dec 2, 2024
1 parent ed6c5f6 commit 1edcf64
Show file tree
Hide file tree
Showing 3 changed files with 87 additions and 54 deletions.
96 changes: 49 additions & 47 deletions src/Filters/stream_filt.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,7 +111,7 @@ function FIRArbitrary(h::Vector, rate_in::Real, Nϕ_in::Integer)
pfb = taps2pfb(h, Nϕ)
dpfb = taps2pfb(dh, Nϕ)
tapsPerϕ = size(pfb, 1)
ϕAccumulator = 1.0
ϕAccumulator = 0.0
ϕIdx = 1
α = 0.0
Δ = Nϕ/rate
Expand Down Expand Up @@ -210,26 +210,26 @@ end
# setphase! set's filter kernel phase index
#
function setphase!(kernel::FIRDecimator, ϕ::Real)
ϕ >= zero(ϕ) || throw(ArgumentError("ϕ must be >= 0"))
ϕ >= zero(ϕ) || throw(DomainError(ϕ, "ϕ must be >= 0"))
xThrowaway = round(Int, ϕ)
kernel.inputDeficit += xThrowaway
nothing
end

function setphase!(kernel::Union{FIRInterpolator, FIRRational}, ϕ::Real)
ϕ >= zero(ϕ) || throw(ArgumentError("ϕ must be >= 0"))
ϕ >= zero(ϕ) || throw(DomainError(ϕ, "ϕ must be >= 0"))
xThrowaway, ϕIdx = divrem(round(Int, ϕ * kernel.Nϕ), kernel.Nϕ)
kernel.inputDeficit += xThrowaway
kernel.ϕIdx = ϕIdx + 1
nothing
end

function setphase!(kernel::FIRArbitrary, ϕ::Real)
ϕ >= zero(ϕ) || throw(ArgumentError("ϕ must be >= 0"))
ϕ >= zero(ϕ) || throw(DomainError(ϕ, "ϕ must be >= 0"))
(ϕ, xThrowaway) = modf(ϕ)
kernel.inputDeficit += round(Int, xThrowaway)
kernel.ϕAccumulator = ϕ*(kernel.Nϕ) + 1.0
kernel.ϕIdx = floor(kernel.ϕAccumulator)
kernel.ϕAccumulator = ϕ * kernel.Nϕ
kernel.ϕIdx = 1 + floor(Int, kernel.ϕAccumulator)
kernel.α = modf(kernel.ϕAccumulator)[1]
nothing
end
Expand Down Expand Up @@ -257,7 +257,7 @@ function reset!(kernel::FIRDecimator)
end

function reset!(kernel::FIRArbitrary)
kernel.ϕAccumulator = 1.0
kernel.ϕAccumulator = 0.0
kernel.ϕIdx = 1
kernel.α = 0.0
kernel.inputDeficit = 1
Expand Down Expand Up @@ -332,11 +332,11 @@ function outputlength(kernel::FIRDecimator, inputlength::Integer)
end

function outputlength(kernel::FIRRational, inputlength::Integer)
outputlength(inputlength-kernel.inputDeficit+1, kernel.ratio, 1)
outputlength(inputlength-kernel.inputDeficit+1, kernel.ratio, kernel.ϕIdx)
end

function outputlength(kernel::FIRArbitrary, inputlength::Integer)
ceil(Int, (inputlength-kernel.inputDeficit+1) * kernel.rate)
ceil(Int, (inputlength-kernel.inputDeficit+1) * kernel.rate - kernel.ϕAccumulator / kernel.Δ)
end

function outputlength(self::FIRFilter, inputlength::Integer)
Expand All @@ -347,63 +347,58 @@ end
#
# Calculates the input length of a multirate filtering operation,
# given the output length
# With RoundDown, inputlength returns the largest input length such that the
# actual output length will be at most the given one.
# With RoundUp, inputlength returns the shortest input length such that the
# actual output length will be at least the given one.
#

function inputlength(outputlength::Int, ratio::Union{Integer,Rational}, initialϕ::Integer)
function inputlength(outputlength::Int, ratio::Union{Integer,Rational}, initialϕ::Integer, r::RoundingMode=RoundDown)
interpolation = numerator(ratio)
decimation = denominator(ratio)
inLen = (outputlength * decimation + initialϕ - 1) / interpolation
floor(Int, inLen)
d = r == RoundUp || r == RoundFromZero ? decimation : 1
inLen = (outputlength * decimation + initialϕ - d) / interpolation
round(Int, inLen, r)
end

function inputlength(::FIRStandard, outputlength::Integer)
function inputlength(::FIRStandard, outputlength::Integer, ::RoundingMode=RoundDown)
outputlength
end

function inputlength(kernel::FIRInterpolator, outputlength::Integer)
inLen = inputlength(outputlength, kernel.interpolation, kernel.ϕIdx)
function inputlength(kernel::FIRInterpolator, outputlength::Integer, r::RoundingMode=RoundDown)
inLen = inputlength(outputlength, kernel.interpolation, kernel.ϕIdx, r)
inLen += kernel.inputDeficit - 1
end

function inputlength(kernel::FIRDecimator, outputlength::Integer)
inLen = inputlength(outputlength, 1//kernel.decimation, 1)
function inputlength(kernel::FIRDecimator, outputlength::Integer, r::RoundingMode=RoundDown)
inLen = inputlength(outputlength, 1//kernel.decimation, 1, r)
inLen += kernel.inputDeficit - 1
end

function inputlength(kernel::FIRRational, outputlength::Integer)
inLen = inputlength(outputlength, kernel.ratio, kernel.ϕIdx)
function inputlength(kernel::FIRRational, outputlength::Integer, r::RoundingMode=RoundDown)
inLen = inputlength(outputlength, kernel.ratio, kernel.ϕIdx, r)
inLen += kernel.inputDeficit - 1

end

# TODO: figure out why this fails. Might be fine, but the filter operation might not being stepping through the phases correcty.
function inputlength(kernel::FIRArbitrary, outputlength::Integer)
inLen = floor(Int, outputlength/kernel.rate)
function inputlength(kernel::FIRArbitrary, outputlength::Integer, r::RoundingMode=RoundDown)
d = r == RoundUp || r == RoundFromZero ? 1 : 0
inLen = floor(Int, (outputlength - d + kernel.ϕAccumulator / kernel.Δ)/kernel.rate) + d
inLen += kernel.inputDeficit - 1
end

function inputlength(self::FIRFilter, outputlength::Integer)
inputlength(self.kernel, outputlength)
function inputlength(self::FIRFilter, outputlength::Integer, r::RoundingMode=RoundDown)
inputlength(self.kernel, outputlength, r)
end


#
# Calculates the delay caused by the FIR filter in # samples, at the input sample rate, caused by the filter process
#

function timedelay(kernel::Union{FIRRational, FIRInterpolator, FIRArbitrary})
(kernel.hLen - 1)/(2.0*kernel.Nϕ)
end

function timedelay(kernel::Union{FIRStandard, FIRDecimator})
(kernel.hLen - 1)/2
end


function timedelay(self::FIRFilter)
timedelay(self.kernel)
end

timedelay(kernel::Union{FIRRational,FIRInterpolator,FIRArbitrary}) =
(kernel.hLen - 1) / (2 * kernel.Nϕ)
timedelay(kernel::Union{FIRStandard,FIRDecimator}) = (kernel.hLen - 1) / 2
timedelay(self::FIRFilter) = timedelay(self.kernel)

#
# Single rate filtering
Expand Down Expand Up @@ -565,18 +560,18 @@ end
# Arbitrary resampling
#
# Updates FIRArbitrary state. See Section 7.5.1 in [1].
# [1] uses a phase accumilator that increments by Δ (Nϕ/rate)
# [1] uses a phase accumulator that increments by Δ (Nϕ/rate)

function update!(kernel::FIRArbitrary)
kernel.ϕAccumulator += kernel.Δ

if kernel.ϕAccumulator > kernel.Nϕ
kernel.xIdx += div(kernel.ϕAccumulator-1, kernel.Nϕ)
kernel.ϕAccumulator = mod(kernel.ϕAccumulator-1, kernel.Nϕ) + 1
if kernel.ϕAccumulator >= kernel.Nϕ
Δx, kernel.ϕAccumulator = divrem(kernel.ϕAccumulator, kernel.Nϕ)
kernel.xIdx += Int(Δx)
end

kernel.ϕIdx = floor(Int, kernel.ϕAccumulator)
kernel.α = kernel.ϕAccumulator - kernel.ϕIdx
kernel.α, foffset = modf(kernel.ϕAccumulator)
kernel.ϕIdx = 1 + Int(foffset)
end

function filt!(
Expand Down Expand Up @@ -645,7 +640,11 @@ function filt(self::FIRFilter{Tk}, x::AbstractVector{Tx}) where {Th,Tx,Tk<:FIRKe
buffer = Vector{promote_type(Th,Tx)}(undef, bufLen)
samplesWritten = filt!(buffer, self, x)

samplesWritten == bufLen || resize!(buffer, samplesWritten)
if Tk <: FIRArbitrary
samplesWritten == bufLen || resize!(buffer, samplesWritten)
else
@assert samplesWritten == bufLen
end

return buffer
end
Expand Down Expand Up @@ -694,11 +693,14 @@ function resample(x::AbstractVector, rate::Real, h::Vector)

# Calculate the number of 0's required
outLen = ceil(Int, length(x)*rate)
reqInlen = inputlength(self, outLen)
reqInlen = inputlength(self, outLen, RoundUp)
reqZerosLen = reqInlen - length(x)
xPadded = [x; zeros(eltype(x), reqZerosLen)]

filt(self, xPadded)
y = filt(self, xPadded)
@assert length(y) >= outLen
length(y) > outLen && resize!(y, outLen)
return y
end

function resample(x::AbstractVector, rate::Real)
Expand Down
5 changes: 2 additions & 3 deletions test/filt_stream.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -362,9 +362,8 @@ end
end
end

# check that these don't throw; the output should actually probably be longer
@test resample(1:2, 3, [zeros(2); 1; zeros(3)]) == [1, 0, 0, 2] # [1, 0, 0, 2, 0, 0]
@test resample(1:2, 3//2, [zeros(2); 1; zeros(3)]) == [1, 0] # [1, 0, 0]
@test resample(1:2, 3, [zeros(2); 1; zeros(3)]) == [1, 0, 0, 2, 0, 0]
@test resample(1:2, 3//2, [zeros(2); 1; zeros(3)]) == [1, 0, 0]
let H = FIRFilter(2.22)
setphase!(H, 0.99)
@test length(filt(H, 1:2)) == 3
Expand Down
40 changes: 36 additions & 4 deletions test/resample.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,10 +100,10 @@ end
@testset "arbitrary ratio" begin
# https://github.com/JuliaDSP/DSP.jl/issues/317
@testset "Buffer length calculation" begin
@test length(resample(sin.(1:1:35546), 1/55.55)) == 641
@test length(resample(randn(1822), 0.9802414928649835)) == 1787
@test length(resample(1:16_367_000*2, 10_000_000/16_367_000)) == 20_000_001
@test resample(zeros(1000), 0.012) == zeros(13)
@test length(resample(sin.(1:1:35546), 1/55.55)) == 640
@test length(resample(randn(1822), 0.9802414928649835)) == 1786
@test length(resample(1:16_367_000*2, 10_000_000/16_367_000)) == 20_000_000
@test resample(zeros(1000), 0.012) == zeros(12)
end
end

Expand Down Expand Up @@ -151,3 +151,35 @@ end
@test isapprox(rc, Nϕ/2, rtol=0.001)
end
end

@testset "inputlength" begin
# FIRDecimator, FIRInterpolator, FIRRational
for _ in 1:1000
M=rand(1:10)//rand(1:10)
H=FIRFilter(zeros(rand(1:100)), M)
if M != 1
setphase!(H, 10*rand())
end
yL = rand(1:100)
@test outputlength(H, inputlength(H, yL)) <= yL < outputlength(H, inputlength(H, yL)+1)
@test outputlength(H, inputlength(H, yL, RoundUp)-1) < yL <= outputlength(H, inputlength(H, yL, RoundUp))
end

# FIRArbitrary
for _ in 1:1000
M = 10*rand()
H = FIRFilter(zeros(rand(1:100)), M)
setphase!(H, 10*rand())
yL = rand(1:100)
@test outputlength(H, inputlength(H, yL)) <= yL < outputlength(H, inputlength(H, yL)+1)
@test outputlength(H, inputlength(H, yL, RoundUp)-1) < yL <= outputlength(H, inputlength(H, yL, RoundUp))
end
let
M = 2.0
H = FIRFilter(resample_filter(M), M)
setphase!(H, timedelay(H))
yL = 200
@test outputlength(H, inputlength(H, yL)) <= yL < outputlength(H, inputlength(H, yL)+1)
@test outputlength(H, inputlength(H, yL, RoundUp)-1) < yL <= outputlength(H, inputlength(H, yL, RoundUp))
end
end

0 comments on commit 1edcf64

Please sign in to comment.