[PD] better tabread4~

[Charles Henry]

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