[PD] better tabread4~

[Charles Henry]

More tables and math...

Let me go ahead and re-print the previous table of fourier transforms
and add some more functions to it, that will simplify things.

remember, these functions are non-zero on the interval from [-a,a] and
zero, elsewhere
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)

and a new set of functions that seem to be more useful, because terms
cancel.  I'll also include another column for the limit as w->0.

f(t)          | lim(w->0) F(w) |           F(w)
|t| - a      |        -a^2         |    2/w^2*(cos(aw)-1)
(|t| - a)^2 |      2/3*a^3      |    4a/w^2 - 4/w^3*sin(aw)
(|t| - a)^3 |     -1/2*a^4      |   -6a^2/w^2 - 12/w^4*(cos(aw)-1)
(|t| - a)^4 |      2/5*a^5      |    8a^3/w^2 - 48a/w^4 + 48/w^5*sin(aw)
(|t| - a)^5 |     -1/3*a^6      |   -10a^4/w^2 + 120a^2/w^4 +
240/w^6*(cos(aw)-1)
(|t| - a)^6 |      2/7*a^7      |    12a^5/w^2 - 240a^3/w^4 +
1440a/w^6 - 1440/w^7*sin(aw)
(|t| - a)^7 |     -1/4*a^8      |   -14a^6/w^2 + 420a^4/w^4 -
5040a^2/w^6 - 10080/w^8*(cos(aw)-1)

There is a general form for these functions, but I'm struggling to put
it in a good way.  I will show the last 4 functions F(w) to show the
pattern.  (if anyone wants to continue this work to bigger and bigger
polynomials, I hope this spares some headaches)

(|t| - a)^4 |    4*2*a^3/w^2  - 4!/(4-3)!*2*a/w^4 + 4!*2/w^5*sin(aw)
(|t| - a)^5 |   -5*2*a^4/w^2 + 5!/(5-3)!*2*a^2/w^4 + 5!*2/w^6*(cos(aw)-1)
(|t| - a)^6 |    6*2*a^5/w^2  - 6!/(6-3)!*2*a^3/w^4 +
6!/(6-5)!*2*a/w^6 - 6!*2/w^7*sin(aw)
(|t| - a)^7 |   -7*2*a^6/w^2 + 7!/(7-3)!*2*a^4/w^4  -
7!/(7-5)!*2*a^2/w^6 - 7!*2/w^8*(cos(aw)-1)

okay, to business.  We want to construct a 6-point polynomial, whose
spectrum falls off at a rate of 1/w^5.  I haven't exactly worked out
the rationale, but this math is tedious and slow and haven't gotten
much of a solid pattern established.

We will construct a fully constrained polynomial, using the functions
(|t|-a)^4 and (|t|-a)^5 using a=1, a=2, and a=3.

To try to keep the terms clear,
b4 is the coefficient of (|t| - 1)^4  on the interval [-1,1]
b5 is the coefficient of (|t| - 1)^5  on the interval [-1,1]
c4 is the coefficient of (|t| - 2)^4  on the interval [-2,2]
c5 is the coefficient of (|t| - 2)^5  on the interval [-2,2]
d4 is the coefficient of (|t| - 3)^4  on the interval [-3,3]
d5 is the coefficient of (|t| - 3)^5  on the interval [-3,3]

We need to set the following constraints:
1.  1/w^2 terms cancel
8*b4 - 10*b5 + 64*c4 - 160*c5 + 216*d4 - 810*d5 = 0

2.  1/w^4 terms cancel
-48*b4 + 120*b5 - 96*c4 + 480*c5 -144*d4 + 1080*d5 = 0

3.  f(0)=1
b4 - b5 + 16*c4 - 32*c5 + 81*d4 - 243*d5 = 1

4. f(1)=0
c4 - c5 + 16*d4 - 32*d5 = 0

5. f(2)=0
d4 - d5 = 0

6.  lim(w->0,  F(w)) = 1
2/5*b4 - 1/3*b5 + 64/5*c4 - 64/3*c5 + 486/5*d4 - 729/3*d5 = 1

Solve by linear algebra,
b4= -125/64
b5= -67/64
c4= 29/32
c5= 5/32
d4= 3/64
d5= 3/64

f(t)=  (|t| < 1) * [ -67/64*(|t| - 1)^5 - 125/64*(|t| - 1)^4 ]
     + (|t| < 2) * [ 5/32*(|t| - 2)^5 + 29/32*(|t| - 2)^4 ]
     + (|t| < 3) * [ 3/64*(|t| - 3)^5 + 3/64*(|t| - 3)^4 ]

F(w) =  1/w^5*( -375/4*sin(w) + 87/2*sin(2w) + 9/4*sin(3w))
         + 1/w^6*(405/2 - 1005/4*cos(w) + 75/2*cos(2w) + 45/4*cos(3w))

okay, so now, we've set the spectrum and impulse response 1st, and we
need to work backwards to find the polynomial interpolation.  We need
to re-write everything in terms of x on [0,1]

from the left:
----------------substitute t = -2 -x  (coefficient of g[-2])

3/64*(| -2 - x| - 3)^5 + 3/64*(| -2 -x| - 3)^4)
=3/64*(x-1)^5 + 3/64*(x-1)^4

----------------substitute t = -1 -x  (coefficient of g[-1])

3/64*(| -1 - x| - 3)^5 + 3/64*(| -1 -x| - 3)^4) + 5/32*(| -1 -x| -
2)^5 + 29/32*(|-1 - x| - 2)^4

=3/64*(x-2)^5 + 3/64*(x-2)^4 + 5/32*(x-1)^5 + 29/32*(x-1)^4

----------------substitute t = -x  (coefficient of g[0])

3/64*(| -x| - 3)^5 + 3/64*(| -x| - 3)^4) + 5/32*(| -x| - 2)^5 +
29/32*(| - x| - 2)^4 + 3/64*(| -x| - 3)^5 + 3/64*(| -x| - 3)^4 -
67/64*(| -x| - 1)^5 - 125/64*(| -x| - 1)^4

=3/64*(x-3)^5 + 3/64*(x-3)^4 + 5/32*(x-2)^5 + 29/32*(x-2)^4 -
67/64*(x-1)^5 - 125/64*(x-1)^4

---------------substitute t = 1 -x  (coefficient of g[1])

3/64*(|1 -x| - 3)^5 + 3/64*(|1 -x| - 3)^4) + 5/32*(|1 -x| - 2)^5 +
29/32*(|1 - x| - 2)^4 + 3/64*(|1 -x| - 3)^5 + 3/64*(|1 -x| - 3)^4 -
67/64*(|1 -x| - 1)^5 - 125/64*(|1 -x| - 1)^4
=3/64*(-x-2)^5 + 3/64*(-x-2)^4 + 5/32*(-x-1)^5 + 29/32*(-x-1)^4 -
67/64*(-x)^5 - 125/64*(-x)^4

= -3/64*(x+2)^5 + 3/64*(x+2)^4 - 5/32*(x+1)^5 + 29/32*(x+1)^4 +
67/64*x^5 - 125/64*x^4

---------------substitute t = 2 -x  (coefficient of g[2])

3/64*(| 2 - x| - 3)^5 + 3/64*(| 2 -x| - 3)^4) + 5/32*(| 2 -x| - 2)^5 +
29/32*(| 2 - x| - 2)^4
=3/64*(-x-1)^5 + 3/64*(-1x-1)^4 + 5/32*(-x)^5 + 29/32*(-x)^4

= -3/64*(x+1)^5 + 3/64*(x+1)^4 - 5/32*x^5 + 29/32*x^4

--------------substitute t = 3 -x  (coefficient of g[3])

3/64*(| 3 - x| - 3)^5 + 3/64*(| 3 -x| - 3)^4)
=3/64*(-x)^5 + 3/64*(-x)^4

= -3/64*x^5 + 3/64*x^4

I'm not sure yet how to condense it into code.  I'll come back to it
again on the weekend.

Chuck