(based on a real story)
by Nikolai V. Shokhirev
An instrument can measure a spectrum at M channels:
In other words, there are the values F_{1}, F_{2}, . . , F_{M} at the points x_{1}, x_{2}, . . , x_{M}. It is necessary to get a spectrum at any arbitrary value x within the interval of measurements ( x_{1}< x_{ }< x_{M} ).
Let us present the spectrum as a polynomial
(1) |
The polynomial coefficients are the solution of the following system of linear equations:
(2) |
The polynomial (1) with coefficients (2) exactly fits the measurements.
Remark 1. The polynomial (1) is a linear combination of the power functions (1, x, x^{2}, . . ). If this basis set is not "natural" for the problem, then the interpolation between channels may be poor.
Remark 2. The matrix A_{mn} = x_{m}^{n} of the system of liner equations (2) is usually ill-conditioned. The reason for this is that in a relatively narrow interval of positive values of x all powers look alike and are almost linear dependent. This makes the solution unstable and inaccurate.
Trick: use y = x - x_{av} in (1) instead of x. Here x_{av} is the middle of the spectrum. This trick partly solve the problem for relatively small M.
Remark 3. The solution of (2) is a time consuming problem ( complexity of M ^{3 }). This solution should be done for each new spectrum.
1. Let us admit that there are strong physical reasons for a polynomial representation (1) for a spectrum.
2. Instead of using pure power basis function (1) we use the polynomials P_{k}(x)
(3) |
orthogonal on the set of measurements:
(4) |
3. The coefficients f_{n} can be easy found in one step:
(5) |
The solution (5) is fast, stable and accurate because of the orthogonality condition (5). For the solution it is actually necessary to know only the matrix P_{n}(x_{m}).
Initially we have a liner independent set of the power functions: 1, x, x^{2}, . . , x^{M-1} . We will construct an orthogonal polynomials recursively.
Suppose that we already constructed P_{0}, P_{1}, . ., P_{n} that are mutually orthogonal and normalized. The next polynomial is made of the new plover x^{n+1} and all previous polynomials:
(6) |
The coefficients q_{k} should be chosen to make (6) orthogonal to all previous polynomials:
(7) |
or
(8) |
The final stage is to normalize the new polynomial:
(9) |
The initial normalized polynomial is
(10) |
Remark 1. Eq. (8) is the subtraction of all previous P_{k}-components from x^{n+1} (actually from all functions of the old basis).
Remark 2. The above relatively complex procedure have to be performed only once for a given set of points x_{1}, x_{2}, . . , x_{M}.
Remark 3. This procedure can be applied to any set of linear independent basis functions.
1. Non-centered points
Let us consider four equidistant points: x_{1} = 1, x_{2} = 1.33333333333333, x_{3 }= 1.66666666666667, x_{4} = 2. The initial functions are the powers of x: 1, x, x^{2}, x^{3}. Below is the matrix of their scalar products <x^{k}|x^{n}>. The power functions are normalized: <x^{n}|x^{n}> = 1.
1 | x | x^{2} | x^{3} | |
1 | 1 | 0.970494958830946 | 0.904912290402514 | 0.833821472585545 |
x | 0.970494958830946 | 1 | 0.980330585324433 | 0.939926806394817 |
x^{2} | 0.904912290402514 | 0.980330585324433 | 1 | 0.988499732201293 |
x^{3} | 0.833821472585545 | 0.939926806394817 | 0.988499732201293 | 1 |
This matrix is related to x_{m}^{n} It is extremely ill-conditioned matrix. Its eigenvalues are: { } = {0.00000022343673427, 0.00172815339820463, 0.187758070383405, 3.81051154185105}. The condition number is defined as cond = Max()/Min(). For the above matrix cond = 1.7 E 7. The larger the condition number is the closer a matrix to a singular one. The condition number is a measure of noise and round off errors amplification. A matrix is considered as well-conditioned if cond << 10^{2}.
The coefficients of new orthogonal polynomials are:
C_{0} | C_{1} | C_{2} | C_{3} | |
P_{0} | 0.5 | 0 | 0 | 0 |
P_{1} | -2.01246117974981 | 1.34164078649987 | 0 | 0 |
P_{2} | 9.5 | -13.5 | 4.5 | 0 |
P_{3} | -61.0446557857441 | 131.257190279237 | -90.5607530887411 | 20.124611797498 |
Below is the matrix <P_{n}|P_{k}> :
P_{0} | P_{1} | P_{2} | P_{3} | |
P_{0} | 1.00000000000000 | 4.44 E-16 | -7.55 E-15 | 5.51 E-14 |
P_{1} | 4.44 E-16 | 1.00000000000000 | 6.11 E-16 | -1.66 E-14 |
P_{2} | -7.55 E-15 | 6.11 E-16 | 0.999999999999999 | -4.42 E-14 |
P_{3} | 5.51 E-14 | -1.66 E-14 | -4.42 E-14 | 0.999999999999935 |
This matrix is close to the identity matrix despite the large cond value.
2. Centered points (using the trick)
Now the points are: x_{1} = -0.5, x_{2} = -0.16666666666667, x_{3 }= 0.16666666666667, x_{4} = 0.5. In this case the matrix <x^{k}|x^{n}> is:
1 | x | x^{2} | x^{3} | |
1 | 1 | -1.12 E-16 | 0.78086880944303 | -1.18 E-16 |
x | -1.12 E-16 | 1 | -1.57 E-16 | 0.959737407008271 |
x^{2} | 0.78086880944303 | -1.57 E-16 | 1 | -2.76 E-16 |
x^{3} | -1.18 E-16 | 0.959737407008271 | -2.76 E-16 | 1 |
The even and odd powers are already orthogonal and the eigenvalues are: {0.0402625929917295, 0.21913119055697, 1.78086880944303, 1.95973740700827}. This gives cond = 48.7.
The constructed orthogonal polynomials include only the power of the same
parity:
C_{0} | C_{1} | C_{2} | C_{3} | |
P_{0} | 0.50000000000000 | 0 | 0 | 0 |
P_{1} | 5.58 E-17 | 1.34164078649987 | 0 | 0 |
P_{2} | -0.6250000000000 | 2.51 E-16 | 4.5000000000000 | 0 |
P_{3} | -1.67 E-16 | -4.58393935387457 | 1.33 E-15 | 20.1246117974981 |
The matrix <P_{n}|P_{k}> now is closer to the identity matrix:
P_{0} | P_{1} | P_{2} | P_{3} | |
P_{0} | 1 | 5.55 E-17 | -5.55 E-17 | 1.11 E-16 |
P_{1} | 5.55 E-17 | 1 | 0 | 6.38 E-16 |
P_{2} | -5.55 E-17 | 0 | 1 | 3.88 E-16 |
P_{3} | 1.11 E-16 | 6.38 E-16 | 3.88 E-16 | 1 |
The polynomial coefficients are calculated with higher accuracy.
©Nikolai V. Shokhirev, 2001-2004