Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix output length of resample #596

Merged
merged 6 commits into from
Dec 2, 2024
Merged

Fix output length of resample #596

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
baf7042
Fix `outputlength(::FIRRational, inputlength)`
martinholters Nov 20, 2024
1f6fd40
Improve `outputlength(::FIRArbitrary, inputlength)`
martinholters Nov 21, 2024
350ea0a
Support specifying rounding mode in `inputlength`
martinholters Nov 20, 2024
12eebb5
Let `resample` return expected length
martinholters Nov 20, 2024
9856def
simplify phase accumulator
wheeheee Nov 27, 2024
27898c7
cosmetics, `DomainError`
wheeheee Nov 27, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
96 changes: 49 additions & 47 deletions src/Filters/stream_filt.jl
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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)
Comment on lines +103 to +106
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These lengths now all match expectations:

julia> 35546 / 55.55
639.8919891989199

julia> 1822 * 0.9802414928649835
1786.0

julia> 16_367_000*2 * 10_000_000/16_367_000
2.0e7

julia> 1000 * 0.012
12.0

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
Loading