I think now would be a good place to pause the implementation and discuss 
if there is a need for this kind of algorithm in the PD community. 

We saw some use cases for the Whittaker-Shannon interpolation where
we gain in quality and/or speed. Namely waveguides and pitch shifters. 

Is there anything else, where we would like to use this interpolation kernel?
Like in general resampling?
In the tabread~ for sample playback?
Are there really quality or speed issues that we could solve there?

For our use case the code I submitted is good enough. But I would be happy 
to spend more time optimizing if there is a need / a broad use of that algorithm.
In that case we need another restructuring and some help from somebody
who is very proficient in writing pd source code. :)

Chuck, I commented your last message below!

Am 2021-05-03 um 03:52 schrieb Charles Z Henry <czhenry@gmail.com>:

On Sun, May 2, 2021 at 4:28 PM Clemens <clemens@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~?


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).

Ok. That’s what I thought as well. Still I would need to read a bit more  
about the subject to understand how the variable delay case relates to 
down-/upsampling. In tabread4a~ you derive the sampling factor from 
the difference in delay time, right?


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.

Good question. :) I’m not sure if Miller has suggestions or examples how 
global variables are handled in PD.

   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.

That’s right, LP_SCALE should be tunable in real-time. 
I think that would need only little change to the code. 

   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

When it gets down to processor specific optimization, I’m not sure how 
much the compiler does already and if it’s really worth the effort. 
I would investigate this if we find that we want to put this code into pd.

Best wishes, 
Clemens



...
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--;
}
}
}