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.

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

WIP/RFC: Fix output length of resample #596

Open
wants to merge 6 commits into
base: master
Choose a base branch
from

Conversation

martinholters
Copy link
Member

This comprises a few changes that can be considered in on their own, but later ones depend on earlier ones for proper results. In commit order:

  • outputlength(::FIRRational, inputlength) is fixed to correctly predict the number of valid samples produced by filt!. (For context: filt! expects a sufficiently large output buffer and returns the number of samples actually written.)
  • outputlength(::FIRArbitrary, inputlength) is improved to do this correctly more often. Unfortunately, it can still be off-by-one sometimes, which seems almost inevitable. The reason is -- somewhat simplified -- that filt! increases the time incrementally, which may be numerically different from considering the whole time-span at once. More complicated in reality, but conceptually: Starting at sample 1 and going in steps of 0.1, how many steps do you need to reach sample 2? Ten, because 10 * 0.1 == 1? No, one more, because 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 < 1. Now, to be sure to have outputlength give the correct result, I think one would have to replicate the same iterative index/phase computations, which seems excessive.
  • Add rounding mode support to inputlength, specifically for RoundUp and RoundDown (where the latter is the old behavior and still the default). With RoundDown, inputlength returns the largest input length such that the actual output length will be at most the given one, i.e. outputlength(H, inputlength(H, yL)) <= yL < outputlength(H, inputlength(H, yL)+1). OTOH, with RoundUp, inputlength returns the shortest input length such that the actual output length will be at least the given one, i.e. outputlength(H, inputlength(H, yL, RoundUp)-1) < yL <= outputlength(H, inputlength(H, yL, RoundUp)).
  • In resample, use inputlength in RoundUp mode to pad the input to the minimum length to get at least ceil(Int, length(x)*rate) output samples, then trim the output to the required size as necessary.

Some points I'd especially appreciate feedback on:

  • Are the current tests sufficient or should I add more? If so, which?
  • Are we confident with the @asserts or should that be more defensive?
  • Is the RoundingMode API for inputlength appropriate? Only RoundUp and RoundDown really make sense there, I guess, so maybe not? (Or should we just drop the RoundDown behavior which might not be needed at all?)

Fixes #186.

test/resample.jl Outdated Show resolved Hide resolved
Comment on lines +103 to +106
@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)
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

Copy link

codecov bot commented Nov 25, 2024

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 97.90%. Comparing base (ed6c5f6) to head (3573b6e).

Additional details and impacted files
@@           Coverage Diff           @@
##           master     #596   +/-   ##
=======================================
  Coverage   97.90%   97.90%           
=======================================
  Files          19       19           
  Lines        3248     3252    +4     
=======================================
+ Hits         3180     3184    +4     
  Misses         68       68           

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

test/filt_stream.jl Outdated Show resolved Hide resolved
@wheeheee
Copy link
Contributor

After some simplification, I think there is a way to make sure that the output length matches. Is there an obstruction to generating ϕAccumulator with a pre-calculated range? Don't have time to try it out now, however.

@martinholters
Copy link
Member Author

martinholters commented Nov 27, 2024

After some simplification, I think there is a way to make sure that the output length matches.

For FIRArbitrary you mean, I suppose? That would be great. I'm curious to see it!

Is there an obstruction to generating ϕAccumulator with a pre-calculated range?

Not sure what you're trying to say here. I guess we could consider changing the state from ϕAccumulator to a counter and then calculate ϕAccumulator from something like mod(counter*Δ, Nϕ). Then it would be much easier to determine the final value for that counter, I think. EDIT: But of course, in the long run, with large values of the counter, the accuracy of counter*Δ may become arbitrarily bad.

I tried to avoid putting too many changes in this PR for sake of reviewability, and rewriting the inner workings of FIRArbitrary is somewhat on the border. It's not immediately required for fixing #186, but surely closely related, so might make sense.

@wheeheee
Copy link
Contributor

wheeheee commented Nov 28, 2024

Yeah, that was what I meant. Didn't consider the numerical error...if only there was a better way to calculate mod(init_acc + counter*Δ, Nϕ).
About the points you raised, I actually think the current implementation makes sense, but for more documentation.
EDIT: oh, it's still WIP, never mind.

@martinholters
Copy link
Member Author

One more thing I'm not 100% sure about: Currently, the target output length is ceil(length(input), ratio). Why ceil? Why not round?

kernel.ϕAccumulator = mod(kernel.ϕAccumulator-1, kernel.Nϕ) + 1
if kernel.ϕAccumulator >= kernel.Nϕ
Δx, kernel.ϕAccumulator = divrem(kernel.ϕAccumulator, kernel.Nϕ)
kernel.xIdx += round(Int, Δx)
Copy link
Member Author

Choose a reason for hiding this comment

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

Suggested change
kernel.xIdx += round(Int, Δx)
kernel.xIdx += Int(Δx)

just to make it explicit that this should not actually round, but just be a type conversion?

kernel.ϕIdx = floor(Int, kernel.ϕAccumulator)
kernel.α = kernel.ϕAccumulator - kernel.ϕIdx
kernel.α, foffset = modf(kernel.ϕAccumulator)
kernel.ϕIdx = 1 + round(Int, foffset)
Copy link
Member Author

Choose a reason for hiding this comment

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

Suggested change
kernel.ϕIdx = 1 + round(Int, foffset)
kernel.ϕIdx = 1 + Int(foffset)

likewise

martinholters and others added 6 commits November 28, 2024 13:58
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.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Resample doesn't resample at given rate
2 participants