[PD] 4-point interpolation changes timbre depending on sample rate

[Charles Z Henry]

> Am 2021-05-03 um 03:52 schrieb Charles Z Henry <czhenry at gmail.com>:
>
> On Sun, May 2, 2021 at 4:28 PM Clemens <clemens at chair.audio> wrote:
>
>
> When there is no aliasing to worry about, i might set the cutoff to pi
> again...
> On low sample rates (22kHz), the lower cutoff is definitely noticeable.
>
>
> I like this thesis you posted earlier
> https://www2.spsc.tugraz.at/www-archive/downloads/Mueller11_DopplerSRC_0.pdf
>
> f(at) <-> 1/|a|*F(w/a)
> This (and two formulas that follow) is listed as Smith's algorithm.  I
> actually got to speak with Julius Smith at the 2012 LAC about this
> formula.  I asked, can't we do any better in terms than O(a*n) number
> of computations for a>1?  He said nope!, but I still have some
> questions there.  Playback with speeds less than 1 always use O(n),
> rather than O(a*n).  I wrote an anti-aliasing external tabread4a~ that
> implements this formula per sample and works pretty well, except it
> becomes expensive when you transpose a few octaves up.
>
>
> I guess that’s the code you are referring to (?):
> https://lists.puredata.info/pipermail/pd-list/2007-03/048397.html
>
> Is it using a lookup of the sync table at certain fixed points?
> If so, would it be „compatible“ in the sense that these points are part
> of the lookup table in delreadsinc~?

It was more a test of this basic formula, rather than anything else.
It's the tabread4~ polynomial, not a sinc.  So I could do a direct
comparison to tabread4~, and show how smoothly it makes the transition
from low speeds to high speeds.  Any choice of the interp(x) function
yields the same basic trade-off: when you make your interpolating
function longer and flatter by factor a, proportional to playback
speeds greater than 1, the cutoff frequency goes to +pi/a rad/s.
Then, resampling at speed a is an evaluation of the convolution with
(1/a*interp(t/a)) with points spaced a apart and it stretches out that
spectrum back to +pi rad/s

The test patch was an obvious case, sweep playback of a 256 point
table with loud partials at f, 2f, 4f, 8f, and you'd see the partials
in an fft graph march right up to the Nyquist frequency and fall off
the edge.  At 8 octaves transposition, this object was doing a
1024-point interpolation (table was wrapped).  At that point, it was
completely silent and using 50% cpu doing very little else.
Meanwhile, you'd hear the aliasing effect from tabread4~ and all it's
partials still bouncing around in the passband.  So, the formula
scales well.  With a little bit of headroom from your LP_SCALE
parameter, there would be some improvement also.

float interp(float x)
{
    float absx=fabsf(x);
    return ((absx<2.0f)*((absx<1.0f)?(1-absx*(0.5f+absx*(1-0.5f*absx))):(1-absx*(1.833333f-absx*(1.0f-0.1666666f*absx)))));
}
...... and the innermost loop that uses interp(x) just uses an extra
term "diff" for difference between last two indexes
...
while((a++)<b)
sum_left+=(*(fp++))*interp((findex-a)/diff);
...
(there's also a sum_right and they are added together to form the
output per sample)
...
*out++=(sum_left+sum_right)/diff;


Your implementation is better as a starting point anyhow.  I've
suspected that tabulating the interpolation function would yield
better performance.