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

[Charles Z Henry]

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.

For a variable delay line (like vd~), the paper's contents are more
relevant and maybe you should consider writing a vdsinc~ object next
(once you've optimized this one).  delread is more like a fixed,
stationary delay line, and I think it's better to default as a literal
interpolator (LP_SCALE = 1).

> Btw. I just implemented sharing of the interpolation table of the
> delreadsinc~ object according to your suggestions.
> It counts the number of references and frees the pointer when no object
> uses it anymore.

    I read through the changes.  For those reading along see the rest
at line 348 of d_delay.c at
https://github.com/chairaudio/pure-data/blob/feature/delreadsinc/src/d_delay.c
 (in progress, current code re: the tables pasted below).
    1. I wonder is this global table properly scoped for Pd.  I
seriously don't know.  But it could be used for multiple objects,
vdsinc~, tabreadsinc~, etc.... So you ought to think about how it
could be re-used between classes and choose a scalable approach now.
    2. I think the declared LP_SCALE variable is a bad approach as
currently implemented.  You can't change it at runtime.  It can be
used in the perform routine instead, with an additional argument to
set the cutoff frequency per object.  This is also good for your
testing, as you'll be able to put objects with different LP_SCALE
values side-by-side for comparison.
    3. I'm convinced this is the right approach for making an optimal
interpolator.  Adding more length pays off over upsampling, because
you evaluate the convolution at a fewer number of points.  Your
implementation needs a lot of tuning---you can squeeze out more
performance.  Focus on the inner most loops.  Also, you'll have to
compare some arm vs intel/amd platforms at some point.  I think you're
closer to the beginning than the end.  This is maybe not what you
wanted to research in the first place


...
typedef struct _sigvdsinc
{
t_object x_obj;
t_symbol *x_sym;
t_float x_sr; /* samples per msec */
int x_zerodel; /* 0 or vecsize depending on read/write order */
t_float x_f;
} t_sigvdsinc;

typedef struct _sigvdsinc_sharedtable
{
// sinc table is held in global variable
t_sample sinc_array[SINC_LEN];
// derivative of sinc funtion for interpolation of sinc function table
t_sample sinc_diff_array[SINC_LEN];
// reference count for shared pointer
int ref_count;
} t_sigvdsinc_sharedtable;

t_sigvdsinc_sharedtable *gSharedTable = NULL;

void sigvdsinc_initialize_sinc_table(t_sample* sinc_array, t_sample*
sinc_diff_array) {
sinc_array[0] = LP_SCALE;
sinc_diff_array[SINC_LEN] = 0;
float a0 = 0.35875;
float a1 = 0.48829;
float a2 = 0.14128;
float a3 = 0.01168;
for (int i=1; i<SINC_LEN; i++) {
float idx = 0.5*(M_PI*i/SINC_LEN - M_PI);
// four term blackmanharris after
https://www.mathworks.com/help/signal/ref/blackmanharris.html
float blackmanharris= a0 - a1*cos(2*idx) + a2*cos(4*idx) - a3*cos(6*idx);
// blackmanharris windowed sinc function
sinc_array[i] =
sin(LP_SCALE*M_PI*(float)i/STEPS_ZC)/(LP_SCALE*M_PI*(float)i/STEPS_ZC)
* blackmanharris * LP_SCALE;
// calculate derivative for table interpolation
sinc_diff_array[i-1] = sinc_array[i]-sinc_array[i-1];
}
}

...
void sigvdsinc_free(t_sigvdsinc *x)
{
// delete shared interpolation table if need be
if(gSharedTable != NULL){
if(gSharedTable->ref_count==1){
free((void*)(gSharedTable));
gSharedTable = NULL;;
}
else{
gSharedTable->ref_count--;
}
}
}