Fix setphase! for pathological cases

[martinholters]

Hit this while trying to get a handle on #186.

For odd number of phases Nϕ and even length of the filter kernel, the old code would yield qidx above Nϕ, namely equal to Nϕ+1, which would result in a BoundsError during resample. Rewrite the code to reduce qidx to 1 in this case and increase inputDeficit by one instead. All other constellations should be unaffected.

For reviewing: I think what the code does is determining xThrowaway and kernel.ϕIdx such that they are both integer and xThrowaway + (kernel.ϕIdx-1) * kernel.Nϕ ≈ ϕ as closely as possible. This should have been the case before and after. But also, we require 1 ≤ kernel.ϕIdx ≤ kernel.Nϕ which the old code failed to ensure in pathological cases.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 97.90%. Comparing base (56773aa) to head (8e47ee5).

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #593      +/-   ##
==========================================
+ Coverage   97.81%   97.90%   +0.09%     
==========================================
  Files          19       19              
  Lines        3243     3246       +3     
==========================================
+ Hits         3172     3178       +6     
+ Misses         71       68       -3     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

---

🚨 Try these New Features:

• Flaky Tests Detection - Detect and resolve failed and flaky tests

[martinholters]

Of course, a similar issue can happen for FIRArbitrary, but only when manually calling setphase!, which probably doesn't happen too often, so this is less severe. I've included a fix nevertheless which is more direct but a bit less elegant.

[wheeheee]

is ϕ required to be positive throughout? if so, is it guaranteed to remain positive after this?

[martinholters]

Good question. First: No, ϕ can become negative here. The modf ensures 0 <= ϕ < 1, so after subtracting 1 it will become negative. In particular, we will have -1/(2Nϕ) <= ϕ < 1-1/(2Nϕ).

What has to be maintained throughout is that ϕIdx is within 1:Nϕ. Here it is initialized to round(ϕ*Nϕ + 1.0) which will fulfill the condition. (Without the subtraction, it could sometimes give Nϕ+1 which will now wrap around to 1.) Further note that ϕAccumulator is initialized to ϕ*Nϕ + 1.0, and we will obviously then fulfill 0.5 <= ϕAccumulator < Nϕ+0.5, where prior to this change, the bounds were 1 and Nϕ+1.

Now consider the update! function which is invoked for every output sample produced:

DSP.jl/src/Filters/stream_filt.jl

Lines 570 to 580 in 56773aa

	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
	end
	kernel.ϕIdx = floor(Int, kernel.ϕAccumulator)
	kernel.α = kernel.ϕAccumulator - kernel.ϕIdx
	end

Now for the worst-case, if ϕAccumulator == 0.5 and kernel.Δ < 0.5, the floor operation will give 0 and we run into an error. But can kernel.Δ < 0.5 happen? It is set to Nϕ/rate, so yes, if Nϕ < 0.5 * rate. And one can specify both upon construction, so it is possible to hit this case.

So this change may indeed lead to a new error - but only in cases that would have errored also without this change, just even earlier. I'm unsure we should land this as is or reduce this PR to only the first commit and worry about FIRArbitrary separately.

Finally, I think there is something fishy going on as the initialization does

DSP.jl/src/Filters/stream_filt.jl

Line 232 in 56773aa

kernel.ϕIdx = round(kernel.ϕAccumulator)

and the update does

DSP.jl/src/Filters/stream_filt.jl

Line 578 in 56773aa

kernel.ϕIdx = floor(Int, kernel.ϕAccumulator)

and it is this discrepancy that prevents this fix to work in all cases. I think they should either both be floor or both be round but so far, I lack understanding of the inner workings of FIRArbitrary to decide. But if it's floor, the fix proposed here would need adjustment, so maybe that's a reason to focus on Union{FIRInterpolator, FIRRational} here.

[martinholters]

As an example for a situation not fixed by this PR:

julia> let H = FIRFilter(122.2)
           setphase!(H, 0.99)
           @test length(filt(H, 1:2)) == 3
       end
Error During Test at REPL[2]:3
  Test threw exception
  Expression: length(filt(H, 1:2)) == 3
  BoundsError: attempt to access 38×32 Matrix{Float64} at index [1, 33] # master
  BoundsError: attempt to access 38×32 Matrix{Float64} at index [1, 0] # this PR

[martinholters]

Ok, thinking about this some more, I really think that setphase! should be consistent with update! here and use floor. Consider this:

julia> let H = FIRFilter([0, 0, 1], 1000.1)
           setphase!(H, 0.03)
           filt(H, [1.0])[1:10]
       end
10-element Vector{Float64}:
 0.96
 0.0
 0.02399360063993594
 0.05599040095990393
 0.08798720127987192
 0.1199840015998399
 0.1519808019198079
 0.18397760223977588
 0.21597440255974387
 0.24797120287971186

(The following samples ramp further up.) The initial 0.96 looks quite off. And indeed, replacing the round with floor, we get

julia> let H = FIRFilter([0, 0, 1], 1000.1)
           setphase!(H, 0.03)
           filt(H, [1.0])[1:10]
       end
10-element Vector{Float64}:
 0.0
 0.0
 0.02399360063993594
 0.05599040095990393
 0.08798720127987192
 0.1199840015998399
 0.1519808019198079
 0.18397760223977588
 0.21597440255974387
 0.24797120287971186

Ah, better. And considering that what setphase! stores in kernel.ϕIdx only affects the first sample, this is probably the way to go. I'll update the PR accordingly.

[martinholters]

Any objections against merging?

[wheeheee]

Not from me, although I only had time to look through it superficially.
Btw, this still doesn't fix the buffer off by one for FirArb right?

[martinholters]

Btw, this still doesn't fix the buffer off by one for FirArb right?

Nope. More on that on a PR I have incoming.