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