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

[Martin Peach]

I think there's a niche for some kind of sinc convolver to be used in
bandlimiting arbitrary waveforms.

Martin

On Mon, May 3, 2021 at 9:26 AM Clemens Wegener <clemens at chair.audio> wrote:
>
> 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
>
> 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--;
> }
> }
> }
>
>
> _______________________________________________
> Pd-list at lists.iem.at mailing list
> UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list