Linear Prolate Functions for Signal Extrapolation with Time Shift

We propose a low complexity iterative algorithm for band limited signal extrapolation. The extrapolation method is based on the decomposition of finite segments of the signal via truncated series of real-valued linear prolate functions. Our theoretical derivation shows that given a truncated series (up to a selectable value) of prolate functions, it is possible to extrapolate the band limited function elsewhere if each extrapolated portion of the function is subject only to moderate truncation errors that we quantify in this paper. The effects of different sources of errors have been analyzed via extensive simulations. We have investigated a property of the signal decomposition formula based on linear prolate functions whereby the integration interval does not need to be symmetric with respect to the origin while time-shifted prolate functions are used in the series.


Introduction
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 [1].Associated with each function, there is an eigenvalue ( ) n c λ and a free parameter c which is a useful descriptor of system performance [2].Some of the mentioned mathematical properties make the prolate functions easily applicable to optics [3].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 [4] 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 [5] and an application to high frequencies have been extensively investigated [6].However, approaches based on the use of the prolate spheroidal wave functions [7] 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 [8], 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 conclusions are drawn in Section 6.

Signal Expansion
As sets of bandlimited functions, orthogonal on the finite interval and orthonormal on the infinite interval, the linear prolate functions ( ) , n c t ψ can be successfully used for the expansion of a generally complex, bandlimited function ( ) the representation is valid for all t , the bandwidth parameter is where 0 Ω represents the finite bandwidth or a cutoff frequency, and 0 t is the time interval.The function ( ) f t is supposed to be 0 Ω -bandlimited.Adopting the criterion of a minimized mean-square error, the expansion coefficients n a in (1) are given by: There is an alternative way to derive the coefficients n a using only the values of ( )

Signal Extrapolation
Our main objective in signal extrapolation using linear prolate functions aims to take advantage of a generalized expression stated in [3], never exploited so far, for the coefficients in (3) which enables the finite interval 0 2t to not be nece- ssarily symmetric with respect to the origin.Hence, for a general interval Substituting y t T ′ = − , the following expression for ( ) Thanks to the significant generalization for the calculation of coefficients . The procedure is then repeated for the i -th iteration and up to the number of iterations that has been set.Specifically, at each iteration i , we form the function which becomes the input for iteration with (6).For the sake of simplicity t ∆ is chosen to be the same at every iteration.Also, ( ) ( )

Numerical Results
A LPFs set with bandwidth parameter ( )

Perfect Knowledge of ( ) f t 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 ( ) known in Mathematica user-defined precision.tional interval which is up to the 60% of half of the time range where the function is known.The presence of the truncation error is discussed in Section 5.

Estimate of ( ) f t in the Integration Interval
We consider a more realistic case when the function ( ) 6) is known in Mathematica user-defined precision only for 0 T = .Figure 4(b) shows signal reconstruction/extrapolation after 16 iterations with a total time shift ( ) as the new input to (6) to make the integral calculation successful.The piecewise polynomial interpolation method presented in [10] has been applied for the accurate computing of the overlap integral and the LPFs set with bandwidth parameter 20π c = has been used.Indeed, in terms of the normalized mean-squared error (NMSE), the method in [10] performs superiorly when compared to the iterative approach proposed in [6] and the generalized PSWFs (prolate spheroidal wave functions) expansion method proposed in [7].Specifically, for comparison purposes, Figure 4 ) outperforms the reference approach ( 1 i = ) 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 36 i = outperforms the shift-approach for 16 i = .The difference quotients in Table 2 calculated between time instants 1.90 and 2.00 show that the shift-approach for 36 i = follows better the slope of the exact function.

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 e T : The first term in the summation in (9) represents the error in the fit of ( ) . Specifically, as reported in [1], the calculation of ).The last statement has to satisfy the condition that for a 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., 0.65 e T = ); increasing N will extend e T and it will still be the best extrapolation up to that point.However, if N 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 ( ) 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 [10], 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 ( ) ( ) ( ) i f t as concatenation of the function with segments of optimum estimates.However, we still observe numerical inaccuracies occurring at the points of concatenations, which is presently under investigation.

Conclusion
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 decomposition 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 [10].Our method has shown to outperform concurrent approaches in terms of the normalized mean-square error of the extrapolated signal.
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 c , the energy concentration of the functions within [ ] of n becomes then a challenging problem of high-precision numerical integration with an absolute precision as well[8].

20π c = and 0 1 t
= 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 0 Ω -bandlimited test function shown in Figure1.

Figure 4 ( 1 i
c) shows results after 36 iterations with a total time shift 0+ , we use the function
Figure 7 for 20π c =