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

Clemens Wegener clemens at chair.audio
Mon May 3 15:24:10 CEST 2021

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 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 <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, 

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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.puredata.info/pipermail/pd-list/attachments/20210503/aa56f9ee/attachment.htm>

More information about the Pd-list mailing list