From baf7042a13a84291f59474f3340988d33910bf39 Mon Sep 17 00:00:00 2001 From: Martin Holters Date: Wed, 20 Nov 2024 16:10:48 +0100 Subject: [PATCH 1/6] Fix `outputlength(::FIRRational, inputlength)` --- src/Filters/stream_filt.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/Filters/stream_filt.jl b/src/Filters/stream_filt.jl index 10d424ce6..98152229e 100644 --- a/src/Filters/stream_filt.jl +++ b/src/Filters/stream_filt.jl @@ -332,7 +332,7 @@ 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) @@ -645,7 +645,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 From 1f6fd401a8efc62904e7ae5e6df64e5ead39b0c2 Mon Sep 17 00:00:00 2001 From: Martin Holters Date: Thu, 21 Nov 2024 12:43:22 +0100 Subject: [PATCH 2/6] Improve `outputlength(::FIRArbitrary, inputlength)` MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- src/Filters/stream_filt.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Filters/stream_filt.jl b/src/Filters/stream_filt.jl index 98152229e..6e9ad8fc4 100644 --- a/src/Filters/stream_filt.jl +++ b/src/Filters/stream_filt.jl @@ -336,7 +336,7 @@ function outputlength(kernel::FIRRational, inputlength::Integer) 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 - 1) / kernel.Δ) end function outputlength(self::FIRFilter, inputlength::Integer) @@ -378,7 +378,7 @@ 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) + inLen = floor(Int, (outputlength + (1 - kernel.ϕAccumulator) / kernel.Δ)/kernel.rate) inLen += kernel.inputDeficit - 1 end From 350ea0a0ef471fe65c9a27613cea27925f2ade26 Mon Sep 17 00:00:00 2001 From: Martin Holters Date: Wed, 20 Nov 2024 16:18:57 +0100 Subject: [PATCH 3/6] Support specifying rounding mode in `inputlength` --- src/Filters/stream_filt.jl | 36 ++++++++++++++++++++---------------- test/resample.jl | 32 ++++++++++++++++++++++++++++++++ 2 files changed, 52 insertions(+), 16 deletions(-) diff --git a/src/Filters/stream_filt.jl b/src/Filters/stream_filt.jl index 6e9ad8fc4..19299eafb 100644 --- a/src/Filters/stream_filt.jl +++ b/src/Filters/stream_filt.jl @@ -347,43 +347,47 @@ 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 + (1 - kernel.ϕAccumulator) / kernel.Δ)/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 - 1) / 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 diff --git a/test/resample.jl b/test/resample.jl index 2c8cbb312..ecf7df331 100644 --- a/test/resample.jl +++ b/test/resample.jl @@ -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 From 12eebb55247458d916d3a4b8aa84bdbe265621b2 Mon Sep 17 00:00:00 2001 From: Martin Holters Date: Wed, 20 Nov 2024 16:19:53 +0100 Subject: [PATCH 4/6] Let `resample` return expected length Fixes #186. --- src/Filters/stream_filt.jl | 7 +++++-- test/filt_stream.jl | 5 ++--- test/resample.jl | 8 ++++---- 3 files changed, 11 insertions(+), 9 deletions(-) diff --git a/src/Filters/stream_filt.jl b/src/Filters/stream_filt.jl index 19299eafb..958a37d5a 100644 --- a/src/Filters/stream_filt.jl +++ b/src/Filters/stream_filt.jl @@ -702,11 +702,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) diff --git a/test/filt_stream.jl b/test/filt_stream.jl index 1c17e1c03..8aa0d17b3 100644 --- a/test/filt_stream.jl +++ b/test/filt_stream.jl @@ -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 diff --git a/test/resample.jl b/test/resample.jl index ecf7df331..202496815 100644 --- a/test/resample.jl +++ b/test/resample.jl @@ -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 From 9856defdb9346467601b849b830a2ec3e4513ac4 Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Wed, 27 Nov 2024 15:46:40 +0800 Subject: [PATCH 5/6] simplify phase accumulator modf --- src/Filters/stream_filt.jl | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/Filters/stream_filt.jl b/src/Filters/stream_filt.jl index 958a37d5a..c37ddb3c2 100644 --- a/src/Filters/stream_filt.jl +++ b/src/Filters/stream_filt.jl @@ -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 @@ -228,8 +228,8 @@ function setphase!(kernel::FIRArbitrary, ϕ::Real) ϕ >= zero(ϕ) || throw(ArgumentError("ϕ 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 @@ -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 @@ -336,7 +336,7 @@ function outputlength(kernel::FIRRational, inputlength::Integer) end function outputlength(kernel::FIRArbitrary, inputlength::Integer) - ceil(Int, (inputlength-kernel.inputDeficit+1) * kernel.rate - (kernel.ϕAccumulator - 1) / kernel.Δ) + ceil(Int, (inputlength-kernel.inputDeficit+1) * kernel.rate - kernel.ϕAccumulator / kernel.Δ) end function outputlength(self::FIRFilter, inputlength::Integer) @@ -382,7 +382,7 @@ end function inputlength(kernel::FIRArbitrary, outputlength::Integer, r::RoundingMode=RoundDown) d = r == RoundUp || r == RoundFromZero ? 1 : 0 - inLen = floor(Int, (outputlength - d + (kernel.ϕAccumulator - 1) / kernel.Δ)/kernel.rate) + d + inLen = floor(Int, (outputlength - d + kernel.ϕAccumulator / kernel.Δ)/kernel.rate) + d inLen += kernel.inputDeficit - 1 end @@ -569,18 +569,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!( From 27898c73f7692610485ff1efcc526f89910e61cb Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Wed, 27 Nov 2024 16:13:41 +0800 Subject: [PATCH 6/6] cosmetics, `DomainError` --- src/Filters/stream_filt.jl | 23 +++++++---------------- 1 file changed, 7 insertions(+), 16 deletions(-) diff --git a/src/Filters/stream_filt.jl b/src/Filters/stream_filt.jl index c37ddb3c2..7fb0fecf0 100644 --- a/src/Filters/stream_filt.jl +++ b/src/Filters/stream_filt.jl @@ -210,14 +210,14 @@ 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 @@ -225,7 +225,7 @@ function setphase!(kernel::Union{FIRInterpolator, FIRRational}, ϕ::Real) 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ϕ @@ -395,19 +395,10 @@ 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