Can someone please explain B-Splines in a simple fashion?

Does anyone have a *good* explaination of B-Splines?  Links to paper would
be great, but ideally I'm looking for someone to put B-Spline theory in
basic english!

I have looked at numerous papers and online documentation, but they
usually jump right to some advanced math and terminology.  My basic
understanding is that it is a piecewise polynomial and can be used for
smooth interpolation of signals.   In my previous job we used it in a
Sample Rate Converter, but looking at the code I never really understood
what the heck it was doing.

All I know was that it was a 5th order FIR filter where a set of equations
using alpha terms generated new coefficients.

Some simple things I'd like to know:  What is a knot?  How does a set of
piecewise polynomial coefficients get converted into FIR coefficients? 
What are good uses for B-splines?  etc.

Basically, I'd to know how a DSP Engineer could use this in the audio
realm.


Thanks,
Jeff

Did you check

http://en.wikipedia.org/wiki/B-spline

or

http://mathworld.wolfram.com/B-Spline.html

In article <8_ednZ2Udb8lmh7YnZ2dnUVZ_rGinZ2d@giganews.com>, "sparafucile17" 
<sparafucile17@hotmail.com> wrote:
>Does anyone have a *good* explaination of B-Splines?  Links to paper would
>be great, but ideally I'm looking for someone to put B-Spline theory in
>basic english!
>
>I have looked at numerous papers and online documentation, but they
>usually jump right to some advanced math and terminology.  My basic
>understanding is that it is a piecewise polynomial and can be used for
>smooth interpolation of signals.   In my previous job we used it in a
>Sample Rate Converter, but looking at the code I never really understood
>what the heck it was doing.
>
>All I know was that it was a 5th order FIR filter where a set of equations
>using alpha terms generated new coefficients.
>
>Some simple things I'd like to know:  What is a knot?  How does a set of
>piecewise polynomial coefficients get converted into FIR coefficients? 
>What are good uses for B-splines?  etc.
>
>Basically, I'd to know how a DSP Engineer could use this in the audio
>realm.
>
>
>Thanks,
>Jeff
>
>

>Did you check
>
>http://en.wikipedia.org/wiki/B-spline
>
>or
>
>http://mathworld.wolfram.com/B-Spline.html
>

I've seen the wikipedia site, the mathworld site is new.  This si pretty
much the same info I've seen elsewhere.  What I'd really like to know is
how to take these equations and convert this to a FIR filter.  I know it
has been done to make interpolation filters, but have no idea how the jump
was made from these equations.

Jeff

sparafucile17 wrote:
> >Did you check
> >
> >http://en.wikipedia.org/wiki/B-spline
> >
> >or
> >
> >http://mathworld.wolfram.com/B-Spline.html
> >
>
> I've seen the wikipedia site, the mathworld site is new.  This si pretty
> much the same info I've seen elsewhere.  What I'd really like to know is
> how to take these equations and convert this to a FIR filter.  I know it
> has been done to make interpolation filters, but have no idea how the jump
> was made from these equations.

there are a couple of online papers:

  http://yehar.com/dsp/deip.pdf

  http://www.mp3-tech.org/programmer/docs/resampler.pdf

and i can send you one that Duane Wise and i did nearly a decade ago
that might answer how polynomial interpolators (of which B-splines are
an example) can be viewed as polyphase FIR filters.

r b-j

>there are a couple of online papers:
>
>  http://yehar.com/dsp/deip.pdf
>
>  http://www.mp3-tech.org/programmer/docs/resampler.pdf
>
>and i can send you one that Duane Wise and i did nearly a decade ago
>that might answer how polynomial interpolators (of which B-splines are
>an example) can be viewed as polyphase FIR filters.
>
>r b-j
>
>

Robert, thanks I appreciate the extra links.  Thank you for re-phrasing my
question!  One of the things I was confused on was viewing polnomial
interpolators as polyphase FIR filters.  I'm not competent on all the
correct DSP terms.  If you send that paper (or post it on the web) that
would be great!

>>there are a couple of online papers:
>>
>>  http://yehar.com/dsp/deip.pdf
>>
>>  http://www.mp3-tech.org/programmer/docs/resampler.pdf
>>
>>and i can send you one that Duane Wise and i did nearly a decade ago
>>that might answer how polynomial interpolators (of which B-splines are
>>an example) can be viewed as polyphase FIR filters.
>>
>>r b-j
>>
>>

Robert,
I'm looking at the paper now and have a couple questions...  When you and
Duane say:

f-2(alpha)
f-1(alpha)
f0(alpha)
f1(alpha)

Do these correlate to FIR B-terms?  If so which ones are they for?  I also
notice that all interpolation types are in terms of alpha.  How do you know
what alpha is?

My previous understand was that alpha = (input Fs)/(output Fs).  So if I
was upsampling from 48k to 96k, then alpha would be 0.5.  Is this right?

Thanks,
Jeff

okay, just to get our terms down, i would recommend to look at this
periodic post i do occasionally about the Nyquist/Shannon Sampling and
Reconstruction Theorem and how that relates to the concept of
REsampling:

http://groups.google.com/group/comp.dsp/msg/e9b6488aef1e2580?fwc=1

some of the math in that has ended up in the Wikipedia article you've
cited.  the idea here is to think of resampling as hypothetically
reconstructing (using the latter "reconstruction half" of the Sampling
Theorem) and then sampling that reconstructed continuous-time signal
again, but this time at different sampling times.  this is described at
the latter "Interpolation, Sample Rate Conversion, etc" section of that
post.  now there you will see the reconstruction at some arbitrary time
t = (m+a)*T.  T is the original sampling period and we're busting up
t/T into integer and fractional parts: "m" and "a".  we have:

                +inf
    x((m+a)*T) = SUM{ x[k]*sinc((m+a)*T/T - k) }
                k=-inf

                +inf
               = SUM{ x[k]*sinc(a + m - k) }
                k=-inf

where t = (m + a)*T
  and m = floor(t/T)
  and a = fract(t/T) = t/T - floor(t/T)

and eventually:

                +inf
    x((m+a)*T) = SUM{ x[m+k] * sinc(a - k) }
                k=-inf

and, after approximating:

                 +16
    x((m+a)*T) = SUM{ x[m+k] * sinc(a-k) * w((a-k)/16) }
                k=-15

which is clearly an FIR filter where x[m+k] are the samples and

    h_a[k] = sinc(a+k) * w((a+k)/16)

are the coefficients (and the coefficient set depends on "a", the
amount of fractional delay) and discrete impulse response.

so "m" tells us which 32 adjacent samples to combine and "a" tells us
how to combine them to get x((m+a)*T) which is at a time (m+a)*T, which
is between time m*T and (m+1)*T because

    0 <= a < 1

(this "a" is the "alpha" in the Duane's and my paper i sent you.)  this
is what "polyphase" FIR interpolation is.  you actually have a
(possibly large) collection of FIR filters, which one depends on what
"a" (or "alpha") is.  so, in resampling, you're doing an FIR, but
you're really doing a bunch of simultaneous FIRs and picking your
output from one of these FIRs, depending on "a".  but you also
recognize that you really don't have to evaluate the FIRs that you
*don't* choose for the output, only the one you *do* choose (which is
determined by the fractional part of the instantaneous time, "a").

now, for simplification with no loss of generality, we can define our
unit of time to be T (normalizing it) and our counting or indexing of
these samples is sorta arbitary and we can start our counting wherever
we want, so we can define m to be zero and say we're interpolating x(a)
between x(0) and x(1).  if we're combining 32 samples, then we need
samples from x(-15), x(-14), ... to x(16).  since "a" is between 0 and
1, we need 16 samples on both sides of the x(a) we're trying to figure
out.  now this is for 32 samples.  in the paper, we reduced that, for
illustration, to 4 samples, two before x(a) and two after.  they are
x(-1), x(0), x(1), and x(2) and those "f_k(a)" factors should really be
called "h_a[k]", but there were two different authors, both with
different notational conventions (design by a committee of two).

i'm gonna leave that a little now and let it stew a bit, and we'll get
back to how to turn this polynomial interpolation into a polyphase FIR,
if you don't figure it out yourself.

sparafucile17 wrote:
> Robert,
> I'm looking at the paper now and have a couple questions...  When you and
> Duane say:
>
> f-2(alpha)
> f-1(alpha)
> f0(alpha)
> f1(alpha)
>
> Do these correlate to FIR B-terms?  If so which ones are they for?  I also
> notice that all interpolation types are in terms of alpha.  How do you know
> what alpha is?
>
> My previous understand was that alpha = (input Fs)/(output Fs).  So if I
> was upsampling from 48k to 96k, then alpha would be 0.5.  Is this right?

"alpha" can mean a lot of things. we throw the symbol around a lot.

r b-j

>
>i'm gonna leave that a little now and let it stew a bit, and we'll get
>back to how to turn this polynomial interpolation into a polyphase FIR,
>if you don't figure it out yourself.
>

Robert,
Yes, I'll need time to step through this and let it really sink in.  But
without doing some scratchpad math, I think my understanding was correct
on the surface.  This is one of those things where you see how someone
else has done it, see how it works, but don't know the thing you should
know the most... Why does it work?

I remember my former mentor explaining to me that B-splines have to
"center" on the point you're trying to interpolate.  And that alpha is a
value between +/- 1 that represents how far from your "center" is the
desired sample. (a point which is in between you're discrete samples)

In fact, in my previous job I saw where this was used for an Aync SRC and
alpha varied in order to accomodate two different clock domains.

Given that Christmas is near, I won't have time to go over this for a
couple weeks, but you should expect at least one follow-up on this...

Thanks for all your help,
Jeff

>now, for simplification with no loss of generality, we can define our
>unit of time to be T (normalizing it) and our counting or indexing of
>these samples is sorta arbitary and we can start our counting wherever
>we want, so we can define m to be zero and say we're interpolating x(a)
>between x(0) and x(1).  if we're combining 32 samples, then we need
>samples from x(-15), x(-14), ... to x(16).  since "a" is between 0 and
>1, we need 16 samples on both sides of the x(a) we're trying to figure
>out.  now this is for 32 samples.  in the paper, we reduced that, for
>illustration, to 4 samples, two before x(a) and two after.  they are
>x(-1), x(0), x(1), and x(2) and those "f_k(a)" factors should really be
>called "h_a[k]", but there were two different authors, both with
>different notational conventions (design by a committee of two).
>

Robert,

Tell me if I'm missing something here, but is it just that the transfer
function equals the f-terms in your paper?

    h_a[k] = f[k]?

Or is it just a replacement for the sinc function, indicating that these
f-terms still need to be windowed?  If I'm way off-base here please step
me through the last bit.  So far, I'm with you in the math.

Jeff

sparafucile17 wrote:
> >now, for simplification with no loss of generality, we can define our
> >unit of time to be T (normalizing it) and our counting or indexing of
> >these samples is sorta arbitary and we can start our counting wherever
> >we want, so we can define m to be zero and say we're interpolating x(a)
> >between x(0) and x(1).  if we're combining 32 samples, then we need
> >samples from x(-15), x(-14), ... to x(16).  since "a" is between 0 and
> >1, we need 16 samples on both sides of the x(a) we're trying to figure
> >out.  now this is for 32 samples.  in the paper, we reduced that, for
> >illustration, to 4 samples, two before x(a) and two after.  they are
> >x(-1), x(0), x(1), and x(2) and those "f_k(a)" factors should really be
> >called "h_a[k]", but there were two different authors, both with
> >different notational conventions (design by a committee of two).
> >
>
> Tell me if I'm missing something here, but is it just that the transfer
> function equals the f-terms in your paper?
>
>     h_a[k] = f[k]?
>
> Or is it just a replacement for the sinc function,

yes, we are replacing the windowed sinc function

> indicating that these
> f-terms still need to be windowed?

no we're replacing the windowed sinc.  no additional windowing is
happening.

>  If I'm way off-base here please step
> me through the last bit.  So far, I'm with you in the math.

you're not off-base.  i'm just having trouble explaining it.  these
polynomial interpolations *can* be thought of as an approximation to:

               +(N/2)
    x(m+a)  =    SUM{ x[m+n] * sinc(a-n) * w((a-n)/(N/2)) }
              n=-(N/2-1)

where the unit of time for "t" is chosen so that the sampling period
"T" is 1.  t = m+a where m = floor(t), the integer part of t. "a" is
the fractional part (0 <= a < 1). and N (which is even) is the number
of samples that we are combining to do our interpolation.  "m" tells us
which N adjacent samples to use and "a" tells us how we are combining
them.  the point of interpolation, t = m+a is in the middle of these N
adjacent samples so there are N/2 samples to the left of "t" and N/2
samples to the right.  we can set "m" to zero with no loss of
generality, so then we are evaluating:

             +(N/2)
    x(a)  =    SUM{ x[n] * sinc(a-n) * w((a-n)/(N/2)) }
            n=-(N/2-1)

or

             +(N/2)
    x(a)  =    SUM{ x[n] * h(a-n) }
            n=-(N/2-1)

where h(t) = sinc(t) * w(t/(N/2))  is an impulse response representing
the windowed sinc() or an impulse response of an actual LPF that we
want to be like the brick wall LPF (a perfect sinc) but is not a brick
wall LPF because it is truncated (which is what windows do).  so, if we
come up with another expression for h(t), we can compute the Fourier
Transform of it to see how good of a brick wall filter it is and then
see how much it attenuates the images (which, in the worst case, all
fold back into the baseband when this is resampled).

so how does this relate to the polynomial interpolation? dunno how to
explain this without just doing an example.  let's start with Lagrange
interpolation and N=4.

4 points: x[-1], x[0], x[1], x[2], and we're interpolating between x[0]
and x[1].  Lagrange interpolation creates an (N-1)th order polynomial
that goes through the N points or adjacent samples (and then we
evaluate it in the middle segment).

    x(a) =    x[-1] * {(a-2)*(a-1)*(a-0)}/{(-1-2)*(-1-1)*(-1-0)}
           +  x[0] * {(a-2)*(a-1)*(a-(-1))}/{(0-2)*(0-1)*(0-(-1))}
           +  x[1] * {(a-2)*(a-0)*(a-(-1))}/{(1-2)*(1-0)*(1-(-1))}
           +  x[2] * {(a-1)*(a-0)*(a-(-1))}/{(2-1)*(2-0)*(2-(-1))}

   for  0 <= a < 1

that's how Lagrange does things.  note that when a = 2, then x(a) =
x[2], when a = 1, then x(a) = x[1], when a = 0, x(a) = x[0], and when a
= -1, then x(a) = x[-1].  and in each term, there are three binomial
factors with a, so each term has a 3rd order polynomial of a and, if
x[-1], x[0], x[1], x[2] are given, then this whole thing adds up to a
3rd order polynomial of "a" where we evaluate "a" between 0 and 1 to
get x(a) between x[0] and x[1].  Are you with me so far?

now, to make the big comparison, compare the above expression for x(a)
to

             +(4/2)
    x(a)  =    SUM{ x[n] * h(a-n) }
            n=-(4/2-1)

or

    x(a) = x[-1]*h(a-(-1)) + x[0]*h(a-0) + x[1]*h(a-1) + x[2]*h(a-2)

    for  0 <= a < 1 .

so for 1 <= t <= 2
     then h(t) = h(a+1) = (a-2)*(a-1)*(a-0)/(-1-2)*(-1-1)*(-1-0)
     for  t = a+1

for 0 <= t <= 1
     then h(t) = h(a) = (a-2)*(a-1)*(a-(-1))/(0-2)*(0-1)*(0-(-1))
     for  t = a

for -1 <= t <= 0
     then h(t) = h(a-1) = (a-2)*(a-0)*(a-(-1))/(1-2)*(1-0)*(1-(-1))
     for  t = a-1

and for -2 <= t <= -1
     then h(t) = h(a-2) = (a-1)*(a-0)*(a-(-1))/(2-1)*(2-0)*(2-(-1))
     for  t = a-2


THAT defines h(t) for -2 <= t <= +2 and h(t) is zero for all other t
(no other samples contribute to x(a)).  then you can plot that h(t), it
looks like Figure 6 in Duane's and my paper and has a corresponding LPF
frequency response.

that frequency response puts deep notches in the gain at integer
multiples of the sampling frequency (where all of these images lie),
but B-splines put in deeper notches where those images are at the
expense of filtering your baseband signal more.  in fact, although
Lagrange and Hermite allow x(n) = x[n] if n is an integer, that is not
true for B-spline interpolation.

this is laborious and i can't do much more with this now.  let this
simmer a while.

r b-j