That's a good reference. I was kind of interested in the hermite polynomials method, but I couldn't figure it out from tabread4h~.c.

On a tangential note, this entry seems interesting: http://en.wikipedia.org/wiki/Hermite_polynomials#Hermite_Functions_as_Eigenf... Probably ought to study that one... Those should have some useful spectral properties.

It's good to see that we haven't strayed too far from other people's work, yet. I had been looking into some orthogonal polynomials for interpolation, but I haven't made heads or tails out of that either.

Chuck

On Wed, Jun 25, 2008 at 12:48 PM, Matt Barber brbrofsvl@gmail.com wrote:

PS -- also for what it's worth, a formula which seems equivalent to the one used in [tabread4c~] is described here:

http://local.wasp.uwa.edu.au/~pbourke/other/interpolation/

it gives the following piece of code:

double CubicInterpolate( double y0,double y1, double y2,double y3, double mu) { double a0,a1,a2,a3,mu2;

mu2 = mu*mu; a0 = y3 - y2 - y0 + y1; a1 = y0 - y1 - a0; a2 = y2 - y0; a3 = y1;

return(a0*mu*mu2+a1*mu2+a2*mu+a3); }

Here's from tabread4c~, in case someone wants to see but couldn't find the file:

   a = wp[-1].w_float;
   b = wp[0].w_float;
   c = wp[1].w_float;
   d = wp[2].w_float;

           a0 = d - c - a + b;
           a1 = a - b - a0;
           a2 = c - a;
           *out++ = ((a0*frac+a1)*frac+a2)*frac+b;

Matt

On Wed, Jun 25, 2008 at 12:06 PM, Charles Henry czhenry@gmail.com wrote:

...

On Wed, Jun 25, 2008 at 10:32 AM, Matt Barber brbrofsvl@gmail.com wrote:

...

For what it's worth, here's supercollider's cubic interpolation function, which differs from csound's and Pd's, which I believe are equivalent:

static float cubicinterp(float x, float y0, float y1, float y2, float y3) { // 4-point, 3rd-order Hermite (x-form) float c0 = y1; float c1 = 0.5f * (y2 - y0); float c2 = y0 - 2.5f * y1 + 2.f * y2 - 0.5f * y3; float c3 = 0.5f * (y3 - y0) + 1.5f * (y1 - y2);

   return ((c3 * x + c2) * x + c1) * x + c0;

}

Matt

a0=(3b-a-3c+d)/2 (same as c3) a1=a-5b/2+2c-d/2 (same as c2) a2=(c-a)/2 (same as c1)

(b is c0)

f(x)=(((a0*x+a1)*x+a2)*x-+b

Actually that's just the same set of coefficients that I named in a previous post... just in a different form. So, that's just the sort of thing we could add to pd-extended.... It's good to know what other people are doing in their software too

Chuck

On 25/06/2008, at 22.33, cyrille henry wrote:

so finally, what should be the name of this object?

is it ok if i remove the other test i made and to use only this one?

I think you a free to name your code what you want. And also to
delete it. I however think people would find it interesting with a
"lib"* of tabread's with different interpolation schemes.

At least i think it would be cool with a collection where also the
interpolation methods of the "others" are included (fx, the SC one).

Thanks for publishing.

(* where ever its one-file-multi-class, one-class-per-file or multi- method-per-class.)

On Thu, Jun 26, 2008 at 12:36 PM, cyrille henry cyrille.henry@la-kitchen.fr wrote:

then, it's possibile to have interpolation using more points, and also interpolation that offer antialiasing.

BTW : i just commited a tabosc4aa~ (aa is for anti aliasing). this is also just a test. it's far from perfect, but lot's better than the original tabread4~ when dealing with high pitch transposition. for now, i did this with a classic 8x upsampling + biquad filter inside the tabread4aa~ (ok, this is more cpu intensive : about 1 to 1.5 % cpu) before i'm doing more work on this, did anyone have experience with this?

To me, the upsampling is not a good approach for anti-aliasing. This means, you have to interpolate, then filter, then sample again.

But you will wind up applying your filter at a lot of points, which will just be discarded.

In fewer steps, (I think) we can modify the interpolation function continuously as the speed changes. Rather than filtering on a long signal, just filter it in place. The part that gets tricky is that the DFT is a lot less relevant, and you have to rely upon continuous functions instead.

I didn't make any progress today... and friday is gone for sure :) I will calculate the spectrum for tabread4c~ on Saturday, I think.

I'm down with the project, but I can't really see the big picture for a tabread/vd/tabosc lib yet.

Chuck

Date: Fri, 27 Jun 2008 00:39:18 -0500 From: "Charles Henry" czhenry@gmail.com Subject: Re: [PD] better tabread4~ To: "Mathieu Bouchard" matju@artengine.ca, pd-list@iem.at Message-ID: 518fe7b20806262239y377186afka8b0f348ecca80d1@mail.gmail.com Content-Type: text/plain; charset=ISO-8859-1

On Thu, Jun 26, 2008 at 11:18 AM, Mathieu Bouchard matju@artengine.ca wrote:

...

I fucked up here. To get a C2 curve you may need to solve an equation system covering the whole table (!). Anyhow, a C1 system is fine enough for most uses, and it would be already much better than pd's.

The thing is that you can only match the 2nd derivatives if you let the 1st derivatives just match but freely float. Then there will be one curve going through all points of the whole table supposing that the 2nd derivative is zero at the beginning and end of the table. Clearly this is a wholly different game because you need to compute a 2nd table to remember what the 1st derivatives are supposed to be and then you can't change anything in the 1st table without recomputing the 2nd table from scratch, or something.

I get what you're saying now. I had to read it a couple times through to see :) You're referring to piecewise cubic polynomials, right? We would wind up with an overdetermined system of equations if we didn't float the 1st *and* 2nd derivatives, which would come out as a linear algebra problem of the size of the table.

but I think it gets even worse. There could be a non-zero null space to the problem. There are infinite solutions to interpolate a table full of zeros, with these conditions. What a mess :)

It's also bad because while a natural cubic spline is conceivable for a tabread (fixing the 2nd derivative to zero on both ends, reading in and keeping all the derived data in a buffer somewhere), you might need a different kind of spline (periodic?) for a tabosc~, and it shouldn't work at all for a vd~ since there is no codified beginning or end to the table (yes-no?).

...

...

or we could set x[0],x[1] and x'[0]=(x[1]-x[-1])/2 and x'[1]=(x[2]-x[0])/2 again, 4 constraints, cubic polynomial, etc...

Seems reasonable. What I want has to have constraints on x'[0] and x'[1]. Those would be a possibility. The problem is that it uses a gap of 2 samples instead of one, so it uses a "blurry" derivative, but the alternative is to have to pick between forward-difference and backwards-difference. The "blurry" derivative happens to be the average of the 1-sample forward-difference and backward-difference.

Which is equivalent to the slope between the 2-sample gap. This would have another advantage over forward- or backward-differences such that going through the table in reverse would produce a symmetric result. (or actually, would it matter after all, since the four points are in the same order whether you're going forward or backward through the table... ?)

...

...

and x''[0]=x[1]-2*x[0]+x[-1] and x''[1]=x[2]-2*x[1]+x[0] 6 constraints, 5th degree polynomial

I think that the replacement for tabread4~ should be another cubic, so that it takes almost the same time to compute it. What I said about C2 was based on a mistaken reading of webpages trying to refresh myself on splines. I should've been more careful.

Yeah, a cubic polynomial makes the most sense for small changes. I haven't ever heard of people interpolating 4 points with a 5th degree polynomial.... but I think I could make it work....

The following bit of code might work to that end as a test, borrowing Cyrille's general notation:

cminusb = c-b; aminusd = a-d;

a0 = aminusd + 3.0 * cminusb; a1 = -2.5f * aminusd - 7.5f * cminusb; a2 = 1.5f * aminusd + 4.5f * cminusb; a3 = 0.5f * (c + a) - b; a4 = 0.5f * (c - a); a5 = b;

*out++ = ((((a0*frac+a1)*frac+a2)*frac+a3)*frac+a4)*frac+a5;

The variables would have to be declared further up, obviously. Also, the compiler should optimize the a5 definition out and just use b (right?), so the above might be clearer and more explicit from a formal standpoint. I'm very prone to algebraic mistakes (especially on friday nights), so if someone else is interested you can check my work and see if you arrive at the same result (x''[n] = x[n-1] -2*x[n]

• x[n+1] , x'[n] = 0.5*(x[n+1] - x[n-1]) for both x[0] and x[1] )

-- the result may be able to be further optimized as well due to some redundancies in the coefficients. Following the naming scheme of the other tests, this might be [tabread4fi~] or some such (fi for "fifth"-order-polynomial).

As this line of experimentation proceeds, it might make sense to develop a set of benchmarks both for quality and performance. One place to start might be to test the residual error between all of the new and old [tabosc~] objects running through a cosine table and an [osc~] with the same frequency and phase, trying out different respective table sizes, and then further test with various cosinesum combinations.

Of course the "ear" test will probably determine things more, especially with sampled data, but it's still a little unclear to me exactly what these interpolations are "supposed" to do when the waveform has transients and discontinuities among the samples -- e.g. what bandwidth should result from moving through a table that's filled with white noise, or what should happen when moving slowly through a table that's filled with alternating 1 -1, or what should a snare-drum or bongo hit sound like at a fifth the speed? These seem to me to be more a matter of taste and interpretation than the cosine tests.

Thanks,

Matt

legends of the pd-list gag department:

Mattieu:

"We need a state-saving object class named [jesus]."

Iohannes:

"this might be the reason why i prefer [lop~] over [cool_filter~] and Pd over reactor."

Mathieu Bouchard wrote:

Contrary to myself in this case, Johannes is hardly joking.

but he is winking, which might qualify as well.

fgmasdr IOhannes

And now, for the sake of absurdity, here is the 7th-degree polynomial for four points which fixes x[-1], x[0], x[1], x[2], x'[0], x'[1], x''[0], and x''[1].

easy-to-read:

/*a7 = 1.25f*(b-c)+(5.f/12.f)*(d-a);*/ a7 = 1.25f*(b-c)+0.4166667f*(d-a);

a6 = 3.5f*a7; a5 = 2.1f*a7; a4 = a6; a3 = 3.1f*a7; a2 = 0.5f*(c+a)-b; a1 = 0.5f*(c-a); a0 = b;

*out++ = ((((((a7*frac-a6)*frac+a5)*frac+a4)*frac-a3)*frac+a2)*frac+a1)*frac+a0;

or somewhat optimized:

a7 = 1.25f*(b-c)+0.4166667f*(d-a); a4 = 3.5f*a7; a2 = 0.5f*(c+a)-b; a1 = 0.5f*(c-a);

*out++ = ((((((a7*frac-a4) * frac+2.1f*a7) * frac+a4) * frac-3.1*a7) * frac+a2) * frac+a1) * frac+b;

13 +'s and 14 *'s

Please test my algebra if you like, but I think it's right. It starts to get a little wobbly when there are more than two points in the same direction, which is to be expected (you can try with the 4-point table generated by the cosinesum method to an array, which makes 4 points and 3 guard points for a total of 7, to see the wobble). I worry that the oscillation will start to go nuts with the higher orders of polynomial -- I'd wonder if a six-point interpolation would gain much over a four-point, but we'd have to test... it might be worth investigating something other than the various kinds of polynomial interpolation, though.

If we find that some of these vastly outperform others for various tasks, it would be good to try to optimize the formulas as much as possible -- the formula in the current [tabread4~] in Pd is rather hard to read, but it's pretty efficient... we might be able to do something similar with these others.

Thanks,

Matt

PS -- Charles, you probably needn't spend your time analyzing this one, as it would most likely be a curiosity rather than anything usable.

On Sat, Jun 28, 2008 at 12:46 PM, Matt Barber brbrofsvl@gmail.com wrote:

On Sat, Jun 28, 2008 at 7:43 AM, cyrille henry cyrille.henry@la-kitchen.fr wrote:

...

Matt Barber a écrit : ...

...

The following bit of code might work to that end as a test, borrowing Cyrille's general notation:

cminusb = c-b; aminusd = a-d;

a0 = aminusd + 3.0 * cminusb; a1 = -2.5f * aminusd - 7.5f * cminusb; a2 = 1.5f * aminusd + 4.5f * cminusb; a3 = 0.5f * (c + a) - b; a4 = 0.5f * (c - a); a5 = b;

*out++ = ((((a0*frac+a1)*frac+a2)*frac+a3)*frac+a4)*frac+a5;

ok, i'll try that. but i don't think adjusting the 2nd derivative is the best thing to do. for me, having a 6 point interpolation would be more important.

well, we will see...

...

This would work, as well. Changing the coefficient order to match the previous code (sorry I bollixed that up before). a4 and a3 are simple ratios of a5, but keeping a2 as explicitly the fourth coefficient... dropping a0, aminusd, and cminusb:

a5 = a - d + 3.0f * (c - b); a2 = 0.5f * (c + a) - b; a1 = 0.5f * (c - a);

*out++ = ((((a5*frac-2.5f*a5)*frac+1.5*a5)*frac+a2)*frac+a1)*frac+b;

I count 11 adds and 10 multiplies, vs. 9 +'s and 7 *'s in the original Lagrange algorithm in Pd, and 10 +'s and 9 *'s in the "third order Hermite" example. I don't know if it's as simple as all that, though, but it would seem to be on at least a similar order of performance. BTW, the following for the C1 interpolation might be slightly more efficient (I left the original line commented out):

a1 = 0.5f * (c - a); /* a2 = a - 2.5 * b + 2.f * c - 0.5f * d; */ a3 = 0.5f * (d - a) + 1.5f * (b - c); a2 = a - b + a1 - a3;

*out++ = ((a3 * frac + a2) * frac + a1) * frac + b;

.

10 +'s 6 *'s

Matt

On Sat, Jun 28, 2008 at 7:43 AM, cyrille henry cyrille.henry@la-kitchen.fr wrote:

Matt Barber a écrit :

...

As this line of experimentation proceeds, it might make sense to develop a set of benchmarks both for quality and performance. One place to start might be to test the residual error between all of the new and old [tabosc~] objects running through a cosine table and an [osc~] with the same frequency and phase, trying out different respective table sizes, and then further test with various cosinesum combinations.

yes, a benchmarking tool would be good to make. but i would not use osc~ as a reference, as it's output come from a table interpolation. (if i did understand pd code, osc~ comes from a linear interpolation of a 512 points table).

(btw, if the interpolation lib is extended to a "better audio synthesis lib", then a better osc~ can be add there to)

the more i digg in pd audio code, the more i think it's important to make this kind of lib. But it would need lot's more work that i can do. and i also don't know much on this subject...

I did not know that about [osc~] -- I'm still going through the code, but [osc~] has often felt kind of "rough around the edges." I suppose that the error in various interpolation schemes from the tabosc objects would begin to converge with larger tables, so it's possible something like that could be used for the benchmark instead. With appropriate frequencies there could probably be some interesting fft tests as well.

I've always liked the sound of SC3's SinOsc ugen... it might be worth looking through that code to see how it differs from [osc~]'s.

Thanks,

Matt

hello Chuck,

i tested this. (and commited) i think tabread6c~ is a bit better than tabread4c~. but differences are more smaller

thx

Cyrille

Charles Henry a écrit :

On Sat, Jun 28, 2008 at 6:43 AM, cyrille henry cyrille.henry@la-kitchen.fr wrote:

...

ok, i'll try that. but i don't think adjusting the 2nd derivative is the best thing to do. for me, having a 6 point interpolation would be more important.

I put together a 6-point interpolation formula and analyzed it. For this I used a 5th degree polynomial, and 6 constraints:

(I want to change up the notation a bit, and not use the letters a, b, c, etc... when switching to 6-point. Y[-2],Y[-1],Y[0], Y[1], Y[2], Y[3] are the points from the table. a5 is the coefficient of x^5, a4 is the coeff. of x^4, ... a0 is a constant term. f(x) is the interpolation polynomial.)

f(0)=Y[0] f(1)=Y[1] f'(0)= 1/12*Y[-2] - 2/3*Y[-1] + 2/3*Y[1] - 1/12*Y[2] f'(1)= 1/12*Y[-1] - 2/3*Y[0] + 2/3*Y[2] - 1/12*Y[3] f''(0)= -1/12*Y[-2] + 4/3*Y[-1] - 5/2*Y[0] + 4/3*Y[1] - 1/12*Y[2] f''(1)= -1/12*Y[-1] + 4/3*Y[0] - 5/2*Y[1] + 4/3*Y[2] - 1/12*Y[3]

This uses improved approximations for the derivative. One advantage of going to 6-point interpolation is to get better numerical derivatives. These approximations of the 1st and 2nd derivatives are accurate up to a higher frequency than before. We can also continue to increase the number of points arbitrarily, without necessarily having to increase the degree of the polynomial. The degree of the polynomial is only determined by the number of constraints, not the number of points.

The coefficients used in this scheme are

a0= Y[0] a1= 1/12*Y[-2] - 2/3*Y[-1] + 2/3*Y[1] - 1/12*Y[2] a2= -1/24*Y[-2] + 2/3*Y[-1] - 5/4*Y[0] + 2/3*Y[1] - 1/24*Y[2] a3= -3/8*Y[-2] + 13/8*Y[-1] - 35/12*Y[0] + 11/4*Y[1] - 11/8*Y[2] + 7/24*Y[3] a4= 13/24*Y[-2] - 8/3*Y[-1] + 21/4*Y[0] - 31/6*Y[1] + 61/24*Y[2] - 1/2*Y[3] a5= -5/24*Y[-2] + 25/24*y[-1] - 25/12*Y[0] + 25/12*Y[1] - 25/24*Y[2] + 5/24*Y[3]

After that, I continued with the impulse response calculations and spectral response calculations, which are a bit disappointing. I'll spare you the equations (for now) and post the graphs. The new traces for the 6-point interpolator are shown in green. It's a little bit hard to see, but the things to look for are the rate at which the graph falls off and the locations of the peaks. The 6-point function has a flatter spectrum, which comes up closer to the Nyquist frequency, and falls off faster. These are the key characteristics of the spectrum we want. The green trace falls off according to 1/w^4, compared to 1/w^3 for tabread4c~ and 1/w^2 for tabread4~

You can see the impulse response in the first graph along with the spectrum. The log vs. dB scale is used same as before, and secondly, I've posted a linear graph, so you can see the difference between functions near the Nyquist frequency (x=pi).

It gives me some ideas for another 6-point scheme, more like tabread4c~, which will fall off at a rate of 1/w^5 and have more notches in the frequency response. I'll work on it a bit, and see how it goes.

Chuck

---

---

---

---

Pd-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list

ok, cool

now, it would also be nice to have a good band limited table reader...

cyrille

Charles Henry a écrit :

Hey, Cyrille,

I kind of thought so... we are quickly running into the law of diminishing returns. I was up late, last night, working on the analysis some more. I think I can have another 6-point version with better characteristics tonight.

Chuck

On Tue, Jul 8, 2008 at 11:35 AM, cyrille henry cyrille.henry@la-kitchen.fr wrote:

...

hello Chuck,

i tested this. (and commited) i think tabread6c~ is a bit better than tabread4c~. but differences are more smaller

thx

Cyrille

Charles Henry a écrit :

...

On Sat, Jun 28, 2008 at 6:43 AM, cyrille henry cyrille.henry@la-kitchen.fr wrote:

...

ok, i'll try that. but i don't think adjusting the 2nd derivative is the best thing to do. for me, having a 6 point interpolation would be more important.

I put together a 6-point interpolation formula and analyzed it. For this I used a 5th degree polynomial, and 6 constraints:

(I want to change up the notation a bit, and not use the letters a, b, c, etc... when switching to 6-point. Y[-2],Y[-1],Y[0], Y[1], Y[2], Y[3] are the points from the table. a5 is the coefficient of x^5, a4 is the coeff. of x^4, ... a0 is a constant term. f(x) is the interpolation polynomial.)

f(0)=Y[0] f(1)=Y[1] f'(0)= 1/12*Y[-2] - 2/3*Y[-1] + 2/3*Y[1] - 1/12*Y[2] f'(1)= 1/12*Y[-1] - 2/3*Y[0] + 2/3*Y[2] - 1/12*Y[3] f''(0)= -1/12*Y[-2] + 4/3*Y[-1] - 5/2*Y[0] + 4/3*Y[1] - 1/12*Y[2] f''(1)= -1/12*Y[-1] + 4/3*Y[0] - 5/2*Y[1] + 4/3*Y[2] - 1/12*Y[3]

This uses improved approximations for the derivative. One advantage of going to 6-point interpolation is to get better numerical derivatives. These approximations of the 1st and 2nd derivatives are accurate up to a higher frequency than before. We can also continue to increase the number of points arbitrarily, without necessarily having to increase the degree of the polynomial. The degree of the polynomial is only determined by the number of constraints, not the number of points.

The coefficients used in this scheme are

a0= Y[0] a1= 1/12*Y[-2] - 2/3*Y[-1] + 2/3*Y[1] - 1/12*Y[2] a2= -1/24*Y[-2] + 2/3*Y[-1] - 5/4*Y[0] + 2/3*Y[1] - 1/24*Y[2] a3= -3/8*Y[-2] + 13/8*Y[-1] - 35/12*Y[0] + 11/4*Y[1] - 11/8*Y[2] + 7/24*Y[3] a4= 13/24*Y[-2] - 8/3*Y[-1] + 21/4*Y[0] - 31/6*Y[1] + 61/24*Y[2] - 1/2*Y[3] a5= -5/24*Y[-2] + 25/24*y[-1] - 25/12*Y[0] + 25/12*Y[1] - 25/24*Y[2] + 5/24*Y[3]

After that, I continued with the impulse response calculations and spectral response calculations, which are a bit disappointing. I'll spare you the equations (for now) and post the graphs. The new traces for the 6-point interpolator are shown in green. It's a little bit hard to see, but the things to look for are the rate at which the graph falls off and the locations of the peaks. The 6-point function has a flatter spectrum, which comes up closer to the Nyquist frequency, and falls off faster. These are the key characteristics of the spectrum we want. The green trace falls off according to 1/w^4, compared to 1/w^3 for tabread4c~ and 1/w^2 for tabread4~

You can see the impulse response in the first graph along with the spectrum. The log vs. dB scale is used same as before, and secondly, I've posted a linear graph, so you can see the difference between functions near the Nyquist frequency (x=pi).

It gives me some ideas for another 6-point scheme, more like tabread4c~, which will fall off at a rate of 1/w^5 and have more notches in the frequency response. I'll work on it a bit, and see how it goes.

Chuck

---

---

---

---

Pd-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list

Sorry it took me so long to make something usable out of that mess. I played around with factoring but it seems like it got me nowhere, so I finally just multiplied out all the polynomials to get the usual form.

(given input points g[-2] g[-1] g[0] g[1] g[2] g[3]

a5 = 3/64*g[-2] + 13/64*g[-1] - 27/32*g[0] + 27/32*g[1] - 13/64*g[2]

• 3/64*g[3]

a4 = -3/16*g[-2] - 19/64*g[-1] + 63/32*g[0] - 9/4*g[1] + 23/32*g[2] + 3/64*g[3] a3 = 9/32*g[-2] - 9/16*g[-1] + 9/16*g[1] - 9/32*g[2] a2 = -3/16*g[-2] + 5/4*g[-1] - 17/8*g[0] + 5/4*g[1] - 3/16*g[2] a1 = 3/64*g[-2] - 19/32*g[-1] + 19/32*g[1] - 3/64*g[2] a0=g[0]

output[x]=((((a5*x+a4)*x+a3)*x+a2*x)+a1)*x+a0

and I did some analysis of the function:

This function is continuous up to the 3rd derivative with derivative approximations:

g'(0)=3/64*g[-2] - 19/32*g[-1] + 19/32*g[1] - 3/64*g[2] g''(0)=-3/8*g[-2] + 5/2*g[-1] - 17/4*g[0] + 5/2*g[1] - 3/8*g[2] g'''(0)=27/16*g[-2] - 27/8*g[-1] + 27/8*g[1] - 27/16*g[2]

but here's the rub. These approximations of the derivatives are horrible. They have terrible spectral response and are not very good for higher frequencies. I'm not sure what this all means in terms of how they sound, but I've got a solid grasp on how this problem works.

1st off: the number of computations is roughly linearly proportional to the number of points, and the degree of the polynomial. 2nd: High frequency response can be obtained by increasing the number of points, beyond the number of points required to constrain the problem for a given degree of polynomial. 3rd: The set of functions specifying the impulse response as sums of (|t|-a)^n*(|t| < a) should be used to construct interpolating polynomials for two clear reasons. First, the lowest degree of polynomial, n, that is used determines the number of continuous derivatives (for n=2, there is 1 continuous derivative, for n=4, there are 3 continuous derivatives). Second, n determines the fastest possible rate of attenuation in the stop-band (for n=2, 1/w^3, for n=4, 1/w^5, etc...)

In the accompanying graphs, the newest spectrum has been added in magenta.

And the question is, where do we go from here.... are there any remaining problems with tabread's? Is the high-frequency response good enough? Do we need faster attenuation?

I think there is little point in trying to increase the rate of attenuation. 1/w^3 is good for a fast interpolator..... 1/w^5 should be good enough for a high-accuracy interpolator (in my opinion).... so if this were carried out to 8-point, 10-point and so on.... we could get better high-frequency response..... ahhhh..... I don't know!

Chuck

On Sat, Jul 19, 2008 at 3:39 PM, Matt Barber brbrofsvl@gmail.com wrote:

Chuck,

Thanks again for this. Quick question: out of curiosity, how much would this differ from the one which has the standard derivative approximations?

Hey, Matt

It is a trade-off. In exchange for getting more roll-off above the Nyquist frequency, we get poorer performance in the high frequencies.

Also, if one wanted to put together the one with the standard approximations, would you use the best approximations available for each derivative, or would you use the ones which come from the same "series" of approximations? I don't know how to call them, but one series of approximation derivations need a 3 points for 1st and 2nd derivatives, and 5 points for 3rd and 4th -- while the next series up needs 5 points for 1st and 2nd and 7 for 3rd and 4th -- can you mix these freely in a 6-point interpolation using the 5-point approximations for everything?

To use better approximations for the derivative in the 6-point setting, we need to go back to using a 3rd degree polynomial. This leads to additional free constraints that can be used to set the 1st derivative with better high-frequency response. We lose the continuity of the second & third derivatives in the process.

We can mix our approximations somewhat freely. I haven't come up with a good rationale for choosing the derivative coefficients, yet....

I guess one important next direction is to work on the anti-aliasing problem -- you mentioned modulating the interpolation coefficients depending on the speed through the table -- would this be a continuous thing, or would there be a pre-defined set of ideal functions among which to choose? Or would this be a matter of figuring out the linear combination of the appropriate anti-aliasing filter (which might need to change with each sample?) and a standard interpolation function? (or am I totally misunderstanding?)

Thanks again,

Matt

The interpolation function is a filter. There would be no need to have an anti-aliasing filter and and interpolation function--there's just the one function. We use the fast interpolating function at speeds <= 1. But we need a general interpolation function as a function of speed that converges to the original function as the speed decreases to 1. This would provide the needed generality and flexibility, while having the same general characteristics of the fast interpolating function on which it is based. I'm open to any ideas on this thing... I think I need to take my eyes off of interpolation for a while, and stop beating up the pd list with tables :)

I've got two basic ideas that I'm playing with. The first is to modify the interpolation function continuously adding a series of "bumps" that are spaced exponentially outward from the original function. If there's some good spectral properties, there could be a way to make a smooth transition and hold the number of calculations to O(log(speed)) instead of O(speed)

My second idea is to replace the points and their derivatives, with filters (low-pass filters for the points and band-pass filters for the derivatives). Then, fit a polynomial as before and interpolate. Like existing schemes, this could be turned into continuous functions for impulse response, which vary as functions of speed.

Any ideas?

Chuck

On Sat, Jul 19, 2008 at 10:33 PM, Matt Barber brbrofsvl@gmail.com wrote:

Right -- wouldn't this be equivalent to doing the (defined) interpolation and the anti-aliasing as a filter in one step? You're modulating the interpolating function to include the effects of the appropriate anti-aliasing filter -- like a one-step sample rate converter. Except, the ratio between the source and target rates is variable. Is this an inappropriate way to be thinking about it?

Yes, that's it.

I guess one problem is how "speed" is measured -- do you just use absolute value of the index delta from one sample to the next (what happens when the index is not a linear function of time)? Or could you fill something like a delay line with past index positions and then use those to find speed as a three- or five-point approximation of the first derivative -- this would add a few samples of delay but might give a better estimate of "speed." Sorry to be dense with the questions, but I want to keep up the best I can. =o)

Yes, that's exactly right, also. The input is a sequence of table indexes, so the speed is the first derivative of the input. There is something problematic, when the user wants to jump between different positions in the table. Those instances shouldn't be treated as playing at high speed--it would just be glitchy. I would suggest a look-ahead method to figure out the speed, when there are rapid changes, to avoid having an error.

...

I've got two basic ideas that I'm playing with. The first is to modify the interpolation function continuously adding a series of "bumps" that are spaced exponentially outward from the original function. If there's some good spectral properties, there could be a way to make a smooth transition and hold the number of calculations to O(log(speed)) instead of O(speed)

My second idea is to replace the points and their derivatives, with filters (low-pass filters for the points and band-pass filters for the derivatives). Then, fit a polynomial as before and interpolate. Like existing schemes, this could be turned into continuous functions for impulse response, which vary as functions of speed.

Any ideas?

Can you give a quick example of the form of each idea?

Not particularly... I'm only working with intuition for the idea so far. It's going to take some insight or inspiration. Do you have any ideas?

In the first, are you adding "bumps" to the interpolator's impulse response? In the second are you saying you would replace a point with the impulse response of a low-pass filter (e.g. in a 5th-degree polynomial with coefficients a0 a1 a2 a3 a4 and a5, instead of matching a0+a1+a2+a3+a4+a5 with y[1] you'd match it with an impulse response centered on y[1])? Would the algebra still be such that you could keep the form for derivatives of the polynomials (in the last example, 2*a2+6*a3+12*a4+20*a5 as the 2nd derivative at y[1]) even though you're matching them with something other than an approximation?

Feeling my way through,

Matt

On Sat, 28 Jun 2008, Matt Barber wrote:

It's also bad because while a natural cubic spline is conceivable for a tabread (fixing the 2nd derivative to zero on both ends, reading in and keeping all the derived data in a buffer somewhere), you might need a different kind of spline (periodic?) for a tabosc~

Yes, the periodic version of the natural cubic spline doesn't have 2nd derivative constraints on endpoints, and instead matches the 1st derivatives of both ends together. It's a small change in one way, but it's not in another, because you can't use the shortcuts associated with tridiagonal matrices.

and it shouldn't work at all for a vd~ since there is no codified beginning or end to the table (yes-no?).

Right. But actually, it's not really useful to go on about this, it was just a mistake of mine because I misread a page.

Which is equivalent to the slope between the 2-sample gap. This would have another advantage over forward- or backward-differences such that going through the table in reverse would produce a symmetric result.

Yes.

(or actually, would it matter after all, since the four points are in the same order whether you're going forward or backward through the table... ?)

This is not the reason, it's because people have more-or-less-defined expectations about an interpolator, and if one assumes symmetricity because of a sample that has been flipped backwards and it introduces a 1-delay sample... one got to know... maybe.

_ _ __ ___ _____ ________ _____________ _____________________ ... | Mathieu Bouchard - tél:+1.514.383.3801, Montréal, Québec

An idea on naming:

It should be possible to exhaust all the different schemes for interpolations involving a given number of points, and it should be possible to do so hierarchically. For four points, e.g., we have:

3rd degree polynomial, setting:

A) x[-1], x[0], x[1], x[2] (Pd's Lagrange default - match point values). B) x[0], x[1], x'[0], x'[1] ("third order Hermite" example in SC3, current [tabosc4c~] - match inner point values and first derivatives)

5th degree polynomial, setting:

C) x[-1], x[0], x[1], x[2], x'[0], x'[1] (match all points, first derivatives at inner points) D) x[0], x[1], x'[0], x'[1], x''[0], x''[1] (match inner points, first derivatives, second derivatives)

7th degree polynomial, setting:

E) x[-1], x[0], x[1], x[2], x'[0], x'[1], x''[0], x''[1] (match all points, first and second derivatives at inner points)

(missing any?)

For polynomial interpolation using four points, if the above is right there are 5 ways to do it, and they are ordered first by degree of polynomial, then from fewest to greatest number of derivatives matched. I tend to agree with other posters who suggested that this kind of organization might best lend itself to one object with interpolation type specified by an argument or message.

The argument could follow one of a couple standards: 1) a numeral or letter in order, as above, or 2) a descriptive argument. The descriptive argument could work something like this:

[tabread4~ array] -- default to cubic Lagrange [tabread4~ array 3] -- explicit cubic Lagrange [tabread4~ array 3d] -- cubic, C1 [tabread4~ array 5dd] -- fifth-degree, C2

etc.

[tabread4~ array 3dd] -- cubic, C2 -- doesn't exist, object doesn't create -- might post a list of available options.

The d or dd could be replaced with whatever, as long as it was descriptive of the level of continuity.

Cons first:

1. No matter what, the argument will be somewhat opaque, and some

people (myself included) would probably have a hard time easily remembering what the arguments mean, which ones are available, and why.

2. List of possible arguments grows a great deal with the number of

points in the interpolation. Also, the first and second derivatives have more numerical approximations with a larger set of points -- one would assume you'd go with the most accurate. Do the third, fourth, etc. derivatives ever get used to derive the polynomial?

3. Having them organized hierarchically by type might imply to some

users that they are also ordered by quality (Lagrange worst, 7th-degree best), while this might not be the case at all.

4. We haven't proven that an exhaustive set is even desirable. Also,

it's unclear which should be the default from a theoretical standpoint (though obviously the current standard should be default for historical purposes).

5. This would force the user to apply an array name if they wanted to

set the interpolation. Same with the class method "set." Also, having a numeral or a symbol to call the interpolation (3 vs. 3d) might make things tough to implement -- this could be changed.

Some Pros:

1. This scheme could be imported directly into Pd (vanilla) with very

little cost, since the default would just be the normal Pd interpolation; this, provided Miller approved such an expansion of something relatively low-level. It would necessitate no change in the implementation of array methods (cosinesum, etc.) since all of these would be 4-point. It would also move the choice of "best" from programmers to users; the benchmarks could be referenced in the documentation and suggestions could be made for different applications (giving both a "what's it do" and a "what's it for" kind of documentation), and give users of vanilla access to a wider selection of interpolation.

2. The external library of other interpolations could be organized

around how many points used ([tabread2~] for linear, [tabread6~] for 6-point). Naming becomes easier and more explicit. Objects which mimic array methods could be made for any of the various interpolations -- [arraycosinesum array size interpolpoints <list of harmonic amps>]

3. The formulas themselves could be tucked into functions for easy

use in other external libraries, but you might not want to introduce a function call for every interpolation into the main tabread, tabosc, and vd object classes.

I'd welcome any comments.

Thanks,

Matt

Frank Barknecht wrote:

An alternative might be "keyword arguments" as used in [declare -lib x -path y] which is the same as [declare -path y -lib x] regardless of order.

hmm, probably a bad example. even though i am not so intimate with [declare], i know for sure that starting Pd with "-lib x -path y" does something different than starting it with "-path y -lib x".

i think this is intentional and definitely has it's uses (another reason why i think that the current preferences system is less than optimal)

fgamdsr IOhannes

On Mon, 30 Jun 2008, IOhannes m zmoelnig wrote:

i know for sure that starting Pd with "-lib x -path y" does something different than starting it with "-path y -lib x".

Does it? then why does -lib store its data in sys_externlist and that pd's main function looks up the libs only after the pdsettings, pdrc and commandline arguments are all loaded?

_ _ __ ___ _____ ________ _____________ _____________________ ... | Mathieu Bouchard - tél:+1.514.383.3801, Montréal, Québec

For polynomial interpolation using four points, if the above is right there are 5 ways to do it, and they are ordered first by degree of polynomial, then from fewest to greatest number of derivatives matched. I tend to agree with other posters who suggested that this kind of organization might best lend itself to one object with interpolation type specified by an argument or message.

The argument could follow one of a couple standards: 1) a numeral or letter in order, as above, or 2) a descriptive argument. The descriptive argument could work something like this:

[tabread4~ array] -- default to cubic Lagrange [tabread4~ array 3] -- explicit cubic Lagrange [tabread4~ array 3d] -- cubic, C1 [tabread4~ array 5dd] -- fifth-degree, C2

etc.

[tabread4~ array 3dd] -- cubic, C2 -- doesn't exist, object doesn't create -- might post a list of available options.

The d or dd could be replaced with whatever, as long as it was descriptive of the level of continuity.

After further thought, there are some more problems with this kind of naming scheme. While the above would be exhaustive for the four-point interpolations, it would not be for a 6-point one. With 6-point interpolations you'd have the opportunity to match more than two first and second derivatives, so there are two ways, for example, to make a 9th-degree C2 curve -- one that fixes all six points, two first derivatives and two second derivatives, and one that fixes four points, four first derivatives, and two second derivatives. An alternate notation would be something like:

[tabread4~ arrayname 4] -- fix all four points (default). [tabread4~ arrayname 2 2] -- fix points at x=0 and 1, and first derivatives at x=0 and 1 [tabread4~ arrayname 4 2 2] fix points at -1 0 1 and 2, f' at 0 and 1, and f'' at 0 and 1

[tabread6~ arrayname 6 4 2 2] fix all 6 points, all 4 first derivatives, and two second and third derivatives (13th-degree polynomial).

Obviously this starts to get ridiculous. In fact if we wanted an exhaustive set of 6-point interpolations, I count 13 of them, assuming you don't ever want to constrain more first derivatives than points, or more second derivatives than first, etc., and assuming you only constrain up to the 2nd derivative. If you want to start constraining 3rd and 4th derivatives (possible with 6 points but not with 4), then it goes up to 25 different interpolations, and up to a 17th-degree polynomial (fixing 6 points, 4 each of the 1st and 2nd derivatives, and 2 each of the 3rd and 4th... yikes!). That would allow for the following: 1 C0 option, 5 C1 options, 7 C2 options, 6 C3 options, and 6 C4 options.

It also starts to put the onus on the user to figure out what arguments to use, and requires that they know a little about the math behind it, which probably shouldn't be necessary. Maybe several classes could be made for the library -- for 4-point you could retain [tabread4~] as the C0 option, then have a [tabread4c1~] option for the continuous first derivative with an argument for the two different ways of doing it (there's a cubic one and a 5th-degree one). Then [tabread4c2~] for the 2 C2 options. This might help keep the classes from being too bloated with options, and would keep things organized by differing levels of continuity, but it means that more classes would have to be maintained...

Not knowing enough of the math to test out the responses, I just don't know which direction to pursue, or whether an exhaustive set of options would make this more of an engineering curiosity than a library for high-quality audio. My intuition is that there are diminishing returns once you start getting higher than 5th- or 7th-degree polynomials, but I also generally hate arbitrary constraints.

Any ideas??

Thanks,

Matt

a good place to find different kind of audio interpolator : http://www.student.oulu.fi/~oniemita/dsp/deip.pdf

cyrille

Matt Barber a écrit :

...

For polynomial interpolation using four points, if the above is right there are 5 ways to do it, and they are ordered first by degree of polynomial, then from fewest to greatest number of derivatives matched. I tend to agree with other posters who suggested that this kind of organization might best lend itself to one object with interpolation type specified by an argument or message.

The argument could follow one of a couple standards: 1) a numeral or letter in order, as above, or 2) a descriptive argument. The descriptive argument could work something like this:

[tabread4~ array] -- default to cubic Lagrange [tabread4~ array 3] -- explicit cubic Lagrange [tabread4~ array 3d] -- cubic, C1 [tabread4~ array 5dd] -- fifth-degree, C2

etc.

[tabread4~ array 3dd] -- cubic, C2 -- doesn't exist, object doesn't create -- might post a list of available options.

The d or dd could be replaced with whatever, as long as it was descriptive of the level of continuity.

After further thought, there are some more problems with this kind of naming scheme. While the above would be exhaustive for the four-point interpolations, it would not be for a 6-point one. With 6-point interpolations you'd have the opportunity to match more than two first and second derivatives, so there are two ways, for example, to make a 9th-degree C2 curve -- one that fixes all six points, two first derivatives and two second derivatives, and one that fixes four points, four first derivatives, and two second derivatives. An alternate notation would be something like:

[tabread4~ arrayname 4] -- fix all four points (default). [tabread4~ arrayname 2 2] -- fix points at x=0 and 1, and first derivatives at x=0 and 1 [tabread4~ arrayname 4 2 2] fix points at -1 0 1 and 2, f' at 0 and 1, and f'' at 0 and 1

[tabread6~ arrayname 6 4 2 2] fix all 6 points, all 4 first derivatives, and two second and third derivatives (13th-degree polynomial).

Obviously this starts to get ridiculous. In fact if we wanted an exhaustive set of 6-point interpolations, I count 13 of them, assuming you don't ever want to constrain more first derivatives than points, or more second derivatives than first, etc., and assuming you only constrain up to the 2nd derivative. If you want to start constraining 3rd and 4th derivatives (possible with 6 points but not with 4), then it goes up to 25 different interpolations, and up to a 17th-degree polynomial (fixing 6 points, 4 each of the 1st and 2nd derivatives, and 2 each of the 3rd and 4th... yikes!). That would allow for the following: 1 C0 option, 5 C1 options, 7 C2 options, 6 C3 options, and 6 C4 options.

It also starts to put the onus on the user to figure out what arguments to use, and requires that they know a little about the math behind it, which probably shouldn't be necessary. Maybe several classes could be made for the library -- for 4-point you could retain [tabread4~] as the C0 option, then have a [tabread4c1~] option for the continuous first derivative with an argument for the two different ways of doing it (there's a cubic one and a 5th-degree one). Then [tabread4c2~] for the 2 C2 options. This might help keep the classes from being too bloated with options, and would keep things organized by differing levels of continuity, but it means that more classes would have to be maintained...

Not knowing enough of the math to test out the responses, I just don't know which direction to pursue, or whether an exhaustive set of options would make this more of an engineering curiosity than a library for high-quality audio. My intuition is that there are diminishing returns once you start getting higher than 5th- or 7th-degree polynomials, but I also generally hate arbitrary constraints.

Any ideas??

Thanks,

Matt

---

Pd-list@iem.at mailing list UNSUBSCRIBE and account-management -> http://lists.puredata.info/listinfo/pd-list

On Wed, Jul 2, 2008 at 1:26 PM, Mathieu Bouchard matju@artengine.ca wrote:

On Tue, 1 Jul 2008, Matt Barber wrote:

...

Any ideas??

Just drop the idea of matching more than two sample points. It's what makes [tabread4~] miss the opportunity to be C1, but it's also in exchange for pretty much nothing. Well, maybe it's not nothing, but I still have no clue about what's the point of matching x[t-1] and x[t+2] for a curve that will only be used from x[t+0] to x[t+1]. When you get to x[t+2], one almost-arbitrarily different cubic has just passed, and you're entering another almost-arbitrarily different cubic, so the cubic used between x[t+0] and x[t+1] seems irrelevant at x[t+2].

Perhaps you have a totally different way of explaining it that would show that two adjacent cubic's continuations into each other are closely related and have special meaning, but if you do, please speak up cause I don't see any of it.

Nope, I don't have any explanation; I'm just a kid! =o)

Seriously though, I tend to agree with you -- this should explain my unease about searching for every polynomial possibility with a certain number of points. I want to help out as much as I can, but I just don't want to be the one to close a door on an option. I am only qualified to deliver some of the formulae and maybe do some of the programming, but I don't pack the mathematical guns to do the kinds of analytical work Chuck has been doing.

On the other hand, doesn't [tabread4~]'s Lagrange interpolator have a continuous 2nd derivative while the [tabread4c~] Hermite one does not? I don't know what that would mean spectrally, if anything. It's the "almost" in the "almost-arbitrary" curves you mention that I don't know how to gauge -- intuitively one could imagine that the more pieces of its surrounding environment are matched, the better the interpolation, but I certainly wouldn't put money on that argument. The other side of the coin is, is it optimal that any points at all should be matched? I could imagine the existence of interpolations which deliver all kinds of spectral benefits but which aren't constrained to pass through any particular sample value...

Also fairly intuitively one could imagine that some of this is signal-dependent, which would seem to stress the need for careful analysis and benchmarking. Curiosity would seem to demand testing a few options to feel out a direction before getting rid of any, testing some absurd options to confirm absurdity. Of course, there's more than a strong possibility that this wheel has already been invented several times over and that the answers have been thoroughly established somewhere. The thing that keeps bothering me is that the smart people who designed Pd and Csound both arrived at the Lagrange interpolation, and those who designed SC3 arrived at the Hermite, and I am certainly not qualified to second-guess any of them. I can do the algebra, though. =oD

Matt

...

On the other hand, doesn't [tabread4~]'s Lagrange interpolator have a continuous 2nd derivative while the [tabread4c~] Hermite one does not?

No. A Lagrange interpolator on N points is a polynomial of degree N-1, and so its Nth derivative is a flat zero function without holes, and so it is infinitely differentiable. However, those are pieced together as a disparate mosaïc in a way that is not even C1 (continuous 1st derivative), which is what prompted Cyrille to work on a replacement in the first place. Note that a discontinuous 1st derivative implies that all other orders of derivatives are discontinuous.

I'm with you on the general piecewise Lagrange not being C1, but I don't think it follows that all other orders are discontinuous -- can't they alternate? At any rate, check out the 2nd derivatives of the piecewise cubic Lagrange. I believe that at x=0 it will be y[-1]

• 2*y[0] + y[1], while at x=1 it will be y[0] - 2*y[1] + y[2].

Therefore, since the terms match at the points on adjacent pieces, the 2nd derivative is continuous even though the first isn't. I'd imagine you could run into this kind of phenomenon especially with piecewise functions. Not sure what it means for the spectral response of the interpolator, though.

Matt

On Thu, 2008-07-03 at 04:01 -0400, Matt Barber wrote:

...

...

On the other hand, doesn't [tabread4~]'s Lagrange interpolator have a continuous 2nd derivative while the [tabread4c~] Hermite one does not?

No. A Lagrange interpolator on N points is a polynomial of degree N-1, and so its Nth derivative is a flat zero function without holes, and so it is infinitely differentiable. However, those are pieced together as a disparate mosaïc in a way that is not even C1 (continuous 1st derivative), which is what prompted Cyrille to work on a replacement in the first place. Note that a discontinuous 1st derivative implies that all other orders of derivatives are discontinuous.

I'm with you on the general piecewise Lagrange not being C1, but I don't think it follows that all other orders are discontinuous -- can't they alternate? At any rate, check out the 2nd derivatives of the piecewise cubic Lagrange. I believe that at x=0 it will be y[-1]

• 2*y[0] + y[1], while at x=1 it will be y[0] - 2*y[1] + y[2].

Therefore, since the terms match at the points on adjacent pieces, the 2nd derivative is continuous even though the first isn't. I'd imagine you could run into this kind of phenomenon especially with piecewise functions. Not sure what it means for the spectral response of the interpolator, though.

yo, i am not too much a math guy, so correct me, if i am talking non-sense, but doesn't the a derivative describe the slope of of the original function at any point? if so, a function with one ore more discontinuities cannot have continuous derivative, because a jump at a certain point would result in a infinitely high value at this point of the derivative. one could argue, that in an analogue continuous world - if the jump is short enough - the peak would be too short to be noticed, but this certainly wouldn't be true in a digital, time discrete domain.

after all, i still don't get, how it could be figured out in the digital domain, whether a curve is continuous or not.

roman

___________________________________________________________ Der frühe Vogel fängt den Wurm. Hier gelangen Sie zum neuen Yahoo! Mail: http://mail.yahoo.de

The reason why the spectrum is necessary for the interpolator is not exactly obvious or simple. It has to do with manipulations in the time domain and what they do, in the frequency domain.

The basis of sampling starts with the following function, a dirac comb

g(t)= sum(k=-inf:inf, dirac-delta( (t-k/fs)) )

It is a series of pulses spaced 1/fs apart. It's Fourier transform

G(f)=integral( g(t)*e^(-i*2*pi*f*t) *dt) G(f)= integral( sum(k=-inf:inf, dirac-delta( (t-k/fs)) )*e^(-i*2*pi*f*t)*dt)

Interchange the order of the integral and summation: G(f)= sum(k=-inf:inf, integral(dirac-delta( (t-k/fs)) )*e^(-i*2*pi*f*t)*dt) G(f)= sum(k=-inf:inf, e^(-i*2*pi*f*k/fs) ) G(f)= 1+ sum(k=1:inf, 2*cos(2*pi*f*k/fs) ) Here we have a periodic function, with period fs. We can re-write this function as a Fourier series of a dirac-delta, in the frequency domain.

integral(-fs/2:fs/2 , dirac-delta(f)*1 df) / integral(-fs/2:fs/2, 1*1 df) = 1/fs so, 1=fs*integral(-fs/2:fs/2 , dirac-delta(f)*1 df) / integral(-fs/2:fs/2, 1*1 df) * 1

likewise, 2*cos(2*pi*f*k/fs)=fs*integral(-fs/2*fs/2: dirac-delta(f)*cos(2*pi*f*k/fs) df)/integral(-fs/2:fs/2: cos(2*pi*f*k/fs)*cos(2*pi*f*k/fs) df) * cos(2*pi*f*k/fs)

This shows that G(f)=fs*dirac-delta(f) on the interval [-fs/2,fs/2] and by periodic extension, we get G(f)=fs*sum(k=-inf:inf, dirac-delta(f-fs*k)

Suppose we have a band-limited function on f=[-fs/2:fs/2] in the time domain, h(t). we can get a sampled version of h(t) by multiplying by g(t) s(t)=h(t)*g(t) = sum(k=-inf:inf, dirac-delta(t-k/fs)*h(t-k/fs)) And multiplication in the time domain is convolution in the frequency domain. S(f)=conv(H(f), G(f)) S(f)=fs*sum(k=-inf:inf, H(f-k*fs))

The result is a periodic spectrum which repeats itself every fs Hz. We still have all the information from h(t), except we no longer know which band of frequencies that information comes from. In real valued frequencies, it could be [0, fs/2], or [fs/2, fs] or [fs, 3*fs/2], etc...

We can only reconstruct a continuous signal choosing one band of frequencies. Our ideal interpolator is:

i(t)=1/fs*sinc(fs*t) where sinc(x)=sin(pi*x)/(pi*x) I(f)={ 1/fs, -fs/2 < f < fs/2 and 0, elsewhere

Convolution in the time domain is multiplication in the frequency domain, so h(t)=conv(i(t), s(t)) means that we multiply our spectrum from -fs/2:fs/2 by factor fs, and everwhere else 0.

Because h(t) was bandlimited on -fs/2:fs/2 , we now have back our original signal up to a difference of a finite number of points.

|h(t) - conv(i(t), s(t))|^2 = 0 (This is by the way, how we can resolve the problem of holes and discontinuities, that vary between functions at only a finite number of points--no jump discontinuities, those have infinite spectral content and no asymptotes, either. We can consider two functions f(x), g(x) to be congruent up to a finite number of points if |f(x)-g(x)|^2=0.)

Okay, so back to the subject at hand. The situation in variable speed playback is this:

We have a continuous function (a sound) in the time domain. We sample it by multiplying by our dirac comb and convolve by an interpolation function. The time-domain function and its spectrum has an impact on the quality of the reconstructed function.

When we play our sound back at a variable speed, we are interpolating (evaluating our convolution) at a series of specified points. prior to playback at specified points, we have a continuous function with spectrum defined on (-inf:inf). When we turn that continuous function into a discrete series again, we stretch or compress this spectrum and resample, producing aliasing and other artifacts of the spectrum.

If we playback at a speed A:

m(t)=conv(h(t)*g(t), i(t) ) n(t)=m(A*t) Then, N(f)=M(f/A)

When we playback at slow speeds, A<1, we map the spectrum of our reconstructed signal onto a smaller interval. e.g. [-fs/2,fs/2] |-> [-fs/2A, fs/2A] This is not much of a problem, except that we now have some high frequencies from above fs/2 now that appear in our desired range of frequencies.

At speeds A>1, we get aliasing. My big idea for anti-aliasing tabread's is to modify the interpolation function continuously with speed changes, so that even as the spectrum is compressed, the cutoff frequency of the interpolation function stays exactly the same. The simplest way to do this is to stretch the function in time.

1/A*i(t/A) has the same spectrum when played back at speed A as the original function i(t)

But as I've said, this tends to be an expensive method to use, and I'm still looking for alternatives.

If you want to compute the spectrum of interpolating polynomials, I have come up with the following method (shown by example).

starting with the tabread4c~ polynomial: g(x)=(-1/2*f[-1] + 3/2*f[0] - 3/2*f[1] + 1/2*d)*x^3 + (f[-1] - 5/2*f[0] + 2*f[1] - 1/2*f[2])*x^2 + (-1/2*f[-1] + 1/2*f[1]) + f[0]

Re-write terms as products with f[-1],f[0],f[1], and f[2] g(x)=(-1/2*x^3 + x^2 - 1/2*x)*f[-1] + (3/2*x^3 - 5/2*x^2 + 1)*f[0] + (-3/2*x^3 + 2*x^2 + 1/2*x)*f[1] + (1/2*x^3 - 1/2*x^2)*f[2]

Now we have something that looks like a convolution. Except... the functions multiplied by each sample are all in terms of a fraction, x. We need to change coordinates to get a centered impulse response of the interpolation function. So, for the first term, we have t= -1 - x (the new coordinate minus the old coordinate)

We have that x belongs to [0,1] t= -1 -x, for t belonging to [-2,-1] t= -x, for t belonging to [-1,0] t= 1 -x, for t belonging to [0,1] t= 2 -x, for t belonging to [1,2]

We make those substitutions for each term in the previous form to find our interpolation function.

i(t)= { 1/2*t^3 + 5/2*t^2 + 4*t + 2, -2<t<-1 -3/2*t^3 - 5/2*t^2 + 1 , -1<t<0 3/2*t^3 - 5/2*t^2 + 1 , 0<t<1 -1/2*t^3 + 5/2*t^2 - 4*t + 2 , 1<t<2 0, elsewhere

This would be a difficult function to take the fourier transform of. Fortunately, it is symmetric and can be condensed.

i(t)={ 3/2*|t|^3 - 5/2*t^2 + 1 , |t|<1 and -1/2*|t|^3 + 5/2*t^2 - 4*|t| + 2 , 1<|t|<2

further, to make it easier to take the transform, we write everything as terms that can be tabulated.

i(t)= (|t| < 1) * (2*|t|^3 - 5*t^2 + 4*|t| - 1) + (|t| < 2) * (-1/2*|t|^3 + 5/2*t^2 - 4*|t| + 2)

One of the terms (on the outermost interval) is exactly the same. The other one is the difference of the two functions.

Next thing is, we write a table of Fourier transforms on the interval [-a,a]. We have to go through each term of the previous function and write down its Fourier transform and add them all up.

f(t) | F(w) 1 | 2/w*sin(aw) |t| | 2a/w*sin(aw) + 2/w^2*(cos(aw)-1) t^2 | 2a^2/w*sin(aw) + 4a/w^2*cos(aw) - 4/w^3*sin(aw) |t|^3 | 2a^3/w*sin(aw) + 6a^2/w^2*cos(aw) - 12a/w^3*sin(aw) - 12/w^4*(cos(aw)-1)

(I have a much longer table, but this is all we will need for this example)

so we make a list of terms and transforms:

on [-1,1]

2|t^3| | 4/w*sin(w) + 12/w^2*cos(w) - 24/w^3*sin(w) - 24/w^4*(cos(w)-1) -5t^2 | -10/w*sin(w) - 20/w^2*cos(w) + 20/w^3*sin(w) 4|t| | 8/w*sin(w) + 8/w^2*(cos(w)-1) -1 | -2/w*sin(w)

on [-2,2]

-1/2*|t|^3 | -8/w*sin(2w) - 12/w^2*cos(2w) + 12/w^3*sin(2w) + 6/w^4*(cos(2w)-1) 5/2*t^2 | 20/w*sin(2w) + 20/w^2*cos(2w) - 10/w^3*sin(2w) -4|t| | -16/w*sin(2w) - 8/w^2*(cos(2w)-1) 2 | 4/w*sin(2w)

Adding up all these terms gives:

I(w) = 1/w^3*(2*sin(2w) - 4*sin(w)) + 1/w^4*(18 - 24*cos(w) + 6*cos(2w))

You can see that a lot of terms cancel. So, we can also change the type of problem and work backwards, setting the spectrum and working backwards to find the corresponding impulse response, and interpolation polynomial.

That's all the major parts of the interpolation theory that I'm working with. I hope that it helps you to see how this process works and how to find the analytical, exact frequency response and impulse response functions of polynomial interpolators.

Later, Chuck

On Sat, Jul 5, 2008 at 1:42 PM, Matt Barber brbrofsvl@gmail.com wrote:

Chuck, thanks for the explanation! One question about this -- would working back from the spectrum possibly mean that you wouldn't necessarily match any of the points (would this mean you could say the impulse response would not have a value of 1 at t=0)? This would be an interesting proposition, and maybe it might help temper the tendency of polynomial interpolation to deliver values >1 or <-1 when the curve has a maximum or minimum between samples? I'm not sure if that would be the result, but I like the sentiment. =o)

Yes, there are some situations where we could construct something that looks good on the face of things, but actually behaves poorly. However, there are a number of free variables involved. I will try to work through an example (so that I can figure it out).

Anyway, thanks again.

Matt