In the early 60’s David Slepian and his colleagues discovered the bandlimited function that is maximally concentrated, in the mean-square sense, within a given time interval; this function is the prolate spheroidal wave function (PSWF) of zero-order.
The linear prolate functions (LPFs) are the one-dimensional version of the prolate spheroidal functions and they form sets of bandlimited functions which are orthogonal and complete over a finite interval. Moreover, unlike other functions, they are also complete and orthogonal over the infinite interval. An additional property is that the finite Fourier transform (FT) of a linear prolate function is proportional to the same prolate function. Although there are other functions which are their own infinite Fourier transform, only the prolate functions enjoy the property for the finite transform: this property uniquely defines the prolate functions  . Associated with each function, there is an eigenvalue and a free parameter which is a useful descriptor of sys- tem performance  . Some of the mentioned mathematical properties make the prolate functions easily applicable to optics  . In particular, we are interested in the problem of determining a bandlimited function from the knowledge of a finite segment of the function, since it is relevant in many practical situations from application to filters in communication systems  to optical systems when, for example, due to intrinsic instrumental limits, only limited observation data are available.
Specifically, in the research area of bandlimited signal extrapolation, there have been contributions with iterative and non-iterative algorithms for extrapolation of signals in the LCT (linear canonical transform) domain that is a generalization of the Fourier transform. The challenges of convergence of algorithms based on the Gerchberg-Papoulis (GP) algorithm  and an applica- tion to high frequencies have been extensively investigated  . However, appro- aches based on the use of the prolate spheroidal wave functions  need to provide efficient ways to compute the prolate functions.
In this paper, we benefit from a proprietary algorithm developed theoretically and implemented numerically by Cada  , for accurate generation of linear prolate functions with desired high precision to use LPFs for signal extrapolation. In what follows, we introduces the basics of signal expansion using the linear prolate functions in Section 2; Section 3 presents our approach to signal extrapolation based on LPFs. In Section 4 and Section 5 simulation results, error analysis and numerical examples are presented and discussed. Finally conclu- sions are drawn in Section 6.
2. Signal Expansion
As sets of bandlimited functions, orthogonal on the finite interval and orthonormal on the infinite interval, the linear prolate functions can be successfully used for the expansion of a generally complex, bandlimited function :
the representation is valid for all , the bandwidth parameter is where represents the finite bandwidth or a cutoff frequency, and is the time interval. The function is supposed to be -bandlimited. Adopting the criterion of a minimized mean-square error, the expansion coefficients in (1) are given by:
There is an alternative way to derive the coefficients using only the values of within the finite observation interval and introducing the eigenvalues . After multiplying (1) by , integrating as reported in  and from the orthogonality properties of the LPFs which are valid on both the finite and infinite interval, one can obtain:
The latest expression for together with (1) states that when the bandlimited function is known over a finite interval of extend then is theoretically known everywhere if one can accurately calculate the coefficients for , the functions and the eigenvalues . can be regarded as the index of energy concentration of each function in . Therefore:
Accurate computing of , and for (where is known as the critical value) turns the orthogonal expansion expression presented in (1) into a signal extrapolation problem. Indeed, for any LPFs set with a fixed , the energy concentration of the functions within decreases as the order increases and for , the signal’s maximum concentration reaches the boundary of the observation interval. Hence the summation of is mostly determined by the first terms whose
individual value is very close to 1, and the series converges to a finite value , as extensively analyzed in  . Accurate estimation of the overlap integral for high orders of becomes then a
challenging problem of high-precision numerical integration with an absolute necessity of having with a high precision as well  .
3. Signal Extrapolation
Our main objective in signal extrapolation using linear prolate functions aims to take advantage of a generalized expression stated in  , never exploited so far, for the coefficients in (3) which enables the finite interval to not be nece- ssarily symmetric with respect to the origin. Hence, for a general interval , relation (4) becomes:
Substituting , the following expression for is obtained:
Thanks to the significant generalization for the calculation of coefficients , one can use (6) to perform the signal extrapolation on beyond the observation interval using an iterative approach.
Specifically, we start from the assumption that the function is perfectly known in the interval . We use (6) with to extrapolate the signal by an interval using already proposed algorithms to obtain accurate calculations of till . Instead of pursuing the more challenging computing for , Formula (6) is re-applied for to extrapolate the signal by an additional interval via accurate calcula- tions of to . The procedure is then repeated for the -th iteration and up to the number of iterations that has been set. Specifically, at each iteration , we form the function which becomes the input for iteration :
with being defined in (6). For the sake of simplicity is chosen to be the same at every iteration. Also,
4. Numerical Results
A LPFs set with bandwidth parameter and is used as the orthogonal basis for the proposed extrapolation method. The functions are discretized in time at a sampling rate of 0.001 for numerical implementation and each discrete sample has a high numerical precision greater than 100 digits. The software Mathematica characterized by high precision computing has been used for the simulations. Extrapolation is carried out on the -bandlimited test function shown in Figure 1.
4.1. Perfect Knowledge of in the Integration Interval
In order to test the proposed approach for signal extrapolation as described in Section 3, we consider the ideal case first. This assumption means that at each iteration of the extrapolation, the function in the integral in (6) is known in Mathematica user-defined precision. Figure 2 shows reconstruction/ extrapolation after the first iteration when the time shift is , the truncation value is and the Mathematica built-in interpolation and integration functions are used to calculate the overlap integral. In Figure 3, we show results for signal reconstruction/extrapolation after the second iteration with an actual application of our time shift approach; specifically and . In both cases, the extrapolation succeeds to extrapolate the signal over an addi-
Figure 1. in the interval .
Figure 2. Extrapolated , no time shift.
Figure 3. Extrapolated , shift .
tional interval which is up to the of half of the time range where the function is known. The presence of the truncation error is discussed in Section 5.
4.2. Estimate of in the Integration Interval
We consider a more realistic case when the function in the integral in (6) is known in Mathematica user-defined precision only for . Figure 4(b) shows signal reconstruction/extrapolation after 16 iterations with a total time shift and .
Figure 4(c) shows results after 36 iterations with a total time shift and . Each progressive shift is equal to . At every iteration , we use the function as the new input to (6) to make the integral calculation successful. The piecewise polynomial interpolation method present- ed in  has been applied for the accurate computing of the overlap integral and the LPFs set with bandwidth parameter has been used. Indeed, in terms of the normalized mean-squared error (NMSE), the method in  performs superiorly when compared to the iterative approach proposed in  and the generalized PSWFs (prolate spheroidal wave functions) expansion method proposed in  . Specifically, for comparison purposes, Figure 4(a) is
(a) (b)(c) (d)
Figure 4. Extrapolation outputs. (a) known in , ; (b) known in , ; (c) known in , ; (d) Reconstruction cases in (b) and (c) (solid line).
obtained by setting and hence , and reproduces results presented in  .
In Figure 5, details are shown for the extrapolation of the portion of the function in the time interval . Despite the effect of accumulated errors, it verifies that given the same truncation value , the shift-approach ( ) outperforms the reference approach ( ) when the extrapolation capability of the reference approach vanishes.
The difference quotients in Table 1 calculated between time instant 1.73 and time instant 1.78 are indices of the curves slope and show that the shift-approach follows better the slope of the exact function. In Figure 6, extrapolation details are shown for the portion of the function in the time interval . Despite the effect of accumulated errors, the shift-approach for outper- forms the shift-approach for . The difference quotients in Table 2 calcu- lated between time instants 1.90 and 2.00 show that the shift-approach for follows better the slope of the exact function.
Table 1. Difference quotient.
Figure 5. Detailed curves.
Figure 6. Detailed curves.
Table 2. Difference quotient.
5. Error Analysis
The proposed method is subject to an inherent series truncation error. Its mean squared error expression is the following, after an extrapolation interval :
The first term in the summation in (9) represents the error in the fit of (defined in (6)) to within the interval . Specifically, as reported in  , the calculation of into the sum form in (10) follows from the fact that are orthogonal on the interval . The term is the truncation error in
the extrapolation interval and then depends on the quantity . It is reasonable to consider the first term not critical for truncation values above , when the energy factor rapidly approaches zero (an example is shown in Figure 7 for ). The last statement has to satisfy the condition that for a
Figure 7. Eigenvalues vs. order for and .
sufficiently large , the integral tends to zero faster than
the corresponding at the denominator of the products . This consideration motivated our work and indirectly highlights again that calculating accurate coefficients for large is critical since both overlap integrals and eigenvalues become very small quantities. This has been a known problem since the 60’s of the last century. The critical aspect is an accurate calculation of , which is now possible  . Under the assumption that the function is initially known in , the sum of truncation errors, computer roundoff and analog to digital conversion errors makes the practical implementation of the proposed iterative approach subject to a total error that has been measured as normalized mean-square error (NMSE) between the original function and the extrapolated signal as:
It is clear from Table 3 that the zero-shift still gives the best results up to the point it is capable to extrapolate correctly (e.g., ); increasing will extend and it will still be the best extrapolation up to that point. However, if is a limit, as in practice always is, our shifting method performs better. This can be interpreted as a consequence of the multiplication by the quantity in (6) which means moving the energy of the LPFs accordingly to the shift of the function. Hence results confirm that our approach works in principle and, as pointed out in  , the chosen algorithm for computing the overlap integral shows its sensitivity to reduced accuracy and noise when applied to an already extrapolated portion of the function. It is also important observing that the reported NMSE values have been calculated for following
Table 3. Normalized mean-square error (NMSE).
notation in (8). The extrapolation can be optimized by reconstructing the extra- polated function as concatenation of the known function with segments of optimum estimates. However, we still observe numerical inaccuracies oc- curring at the points of concatenations, which is presently under investigation.
In this paper, we have proposed and implemented a low complexity iterative algorithm for bandlimited signal extrapolation based on orthogonal projections over real-valued eigenvectors: the linear prolate functions. The method is valid for an arbitrary large range of frequencies with immediate applications in signal processing. The main contribution of our work is a theoretical derivation such that given a truncated series (up to a selectable value) of prolate functions, it is possible to extrapolate the bandlimited function (initially known in a limited time interval) elsewhere if each extrapolated portion of the function is subject only to moderate series truncation errors. These errors are controllable by the depth of extrapolation at each iteration. By doing so and with the aim of finding an alternative solution to the initial problem of implementing an accurate summation of infinite terms, we have investigated a property of the signal de- composition formula based on LPFs according to which the integration interval does not need to be symmetric with respect to the origin while time-shifted prolate functions are used in the summation. Also, we have investigated the effects of different sources of errors by implementing and analyzing the iterative algorithm as a generalization of the special case presented in  . Our method has shown to outperform concurrent approaches in terms of the normalized mean-square error of the extrapolated signal.