Piecewise Iterative Extrapolation Method for Bandlimited Signal

: In some high spectrally efﬁcient communication structures, the Gerchberg-Papoulis extrapolation algorithm can be used to reconstruct nonorthogonal signals. However, the extrapolation efﬁciency and accuracy degrades when the energy of known signal is small or the iterative ﬁlter bandwidth is large. Focused on these drawbacks, a piecewise iterative extrapolation method is proposed to extrapolate bandlimited signals with higher extrapolation efﬁciency and accuracy. In this paper, the amplitude change of the unknown part over iterations is analyzed and the feasibility of the proposed extrapolation method is discussed. Numerical simulation shows the accuracy superiority of proposed method with the frequency offset compared to the GP algorithm with ﬁxed computational complexity. Besides, the feasibility of the proposed method is veriﬁed in fractional Fourier domain with performance advantages.


Introduction
The extrapolation for bandlimited signals is one of essential research objects in signal processing, wireless communication, and positioning scenarios where the transmitted signals are always bandlimited. For example, cosine signals and orthogonal frequency division multiplexing (OFDM) signals concentrate its information within a finite bandwidth in the frequency domain, as with the chirp-based signals in the fractional Fourier domain. Due to transmission interference or device failure, only partial data can be detected or used for subsequent processing [1][2][3]. Besides, the ever-evolving communications scenarios requires improved throughput to satisfy larger access demand, diversity, and mobility [4][5][6][7][8]. Signals are designed to be transmitted partially to obtain higher spectral efficiency according to the high-throughput requirements [8][9][10][11]. The reconstruction method in terms of the received partial signals is necessary to maintain the system workability.
The Gerchberg-Papoulis (GP) algorithm, as a classic iterative extrapolation method for frequency bandlimited signal, can recover the entire signal from a part of it in the time domain [12,13]. As proven in [12], the received signal can gradually approximate the original signal, as the iterations tend to be infinite. Due to its high practicability and reliability, the GP algorithm has been adopted in many applications such as spectral estimation, high spectral efficiency communication and signal reconstruction [13][14][15]. When extrapolating the discrete signals, the influence of analyticity on the convergence of the GP algorithm should be considered. The feasibility and convergence of a discrete GP algorithm is provided and the effective sampling period is analyzed in [16,17]. The core computation uses either Fourier transform (FT) for analytic signals or discrete Fourier transform (DFT) for discrete signals, and a congruent relationship is established with continuous or discrete prolate spheroidal functions. To enlarge the application scope,
As a generalized form of FT, p-order FrFT of a time domain signal f (t) is defined as where the kernel function K(p; t, u) is represented as where ϕ = pπ/2 is the rotation angle of the time-frequency plane, variant p ∈ R with the period of 4, and n ∈ Z. When p = 1 and ϕ = π/2, (1) reduces to FT, and when p = −1 and ϕ = −π/2, (1) reduces to inverse FT (IFT). If p = 0, F 0 (u) = f (t).
In terms of a σ-fractional bandlimited signal f (t), the values of its FrFT form should be zero outside the interval [−σ, σ] as According to Parseval principle, the energy should be finite where σ > 0 denotes the bandwidth of the original signal.
In the actual scenarios, some signals are not definitely bandlimited, but relatively. Take the OFDM signal as an example, the transmitted information is expected to be modulated into a series of orthogonal trigonometric basis functions. Due to the finite duration, energy dispersion exists in the actual OFDM signals, resulted in non-bandlimited characteristic (regardless of the lowpass filter processing). In that case, the OFDM signal can still be regarded as bandlimited since the dispersion is ignorable outside the original bandwidth. The compressed OFDM communication structure [9] is also proposed based on the relatively bandlimited characteristic.

Traditional GP Extrapolation Algorithm
GP extrapolation algorithm is proposed for analytic frequency bandlimited signals. To extrapolate an σ-bandlimited signal f (t), the observed g(t) is assumed to be where the operator U T {} remains the part of f (t) within [−T, T] and T is positive. The extrapolation is realized with iterations. For the first iteration i = 1, the initial signal is f 0 (t) = g(t). Defining f i−1 (t) as the output result of (i − 1)th iteration, and the ith iteration in GP algorithm consists of 6 main steps: Input f i−1 (t) as the initial signal of the ith iteration.
First, the signal f i−1 (t) is transformed into the frequency domain with FT as F i−1 (ω). Second, F i−1 (ω) is lowpass filtered within |ω| ≤ σ. Third, the signal is transformed to the time domain with IFT. Finally, the part within [−T, T] is replaced by the observed part g(t).
Output f i (t) as the result of the ith iteration.

Amplititude Variation during Extrapolation
The performance of the extrapolation method can be evaluated with mean-square error (MSE). In the case of Ω = σ for GP algorithm, the convergence can be explained with the relationship of f (t), f i−1 (t) and f i (t).
Define an MSE function as M{ f (t), The first greater-than relationship is attributed to that both f and g i is bandlimited within [−Ω, Ω], but f i is not. The lowpass filter maintains the bandlimited characteristics of g i in the 2nd step of GP algorithm in Section 2.2, thus the difference between f and f i is larger than that of f and g i . The second greater-than relationship is due to that the known signal within t ∈ [−T, T] replaces the extrapolated part in the 4th step. The known part is obviously more accurate than the extrapolated one hence reducing the MSE.
According to GP algorithm referred in Section 2.2, the iteration filter in step 2 is set to be with the same bandwidth as that of the original signal. However, the actual signal to be extrapolated is finite in the time domain and it is not strictly bandlimited. Since the signal bandwidth may shift, and the bandwidth of iteration filter may be chosen larger. In addition, the maximum iteration is always set considering the computation requirement, which is not a fixed number. Therefore, the influence of the iteration filter and maximum iteration (stopping criterions) on MSE and amplitude variation should be discussed, respectively.
As for a basic isometric time domain signal f (t) = cos(2πt), the variation in the signal waveform shown in Figure 1 can reflect the extrapolation process more intuitively. The original cosine signal is sampled uniformly by 2048 points within t ∈ (−4, 4), and the duration of signal to be extrapolated is determined by the truncated rate α. For example, duration of signal to be extrapolated is within t ∈ (−1, 1) for α = 0.25. Here the stopping criterions is set to be maximum iterations ITE, and the initial known signal g(t) = f (t), |t| < 1 . The bandwidth of iteration filter is chosen to be the center frequency of f 1 in the condition of Ω = σ. Two conditions are considered that are the waveforms under different maximum iterations (simplified as ITE in the simulation figures) and under different bandwidth values of iteration filter. The MSE of f and extrapolatedf = f ITE are shown in Figure 1, and waveforms of the extrapolated signal under different conditions are shown in Figure 2. The "r f /s " in the legend denotes the ratio of extrapolation filter bandwidth and signal bandwidth.   From the above analysis and simulation results, it can be seen that the GP algorithm obtains a convergent error through infinite iterations but the iterations are restrained in the actual extrapolation process. Therefore, two main problems should be focused. One is the slow convergence (MSE bottleneck) when the ratio of proportion of the known part is small or the iterations reach some threshold. Another is the inaccurate extrapolation (low MSE) when the iteration filter bandwidth is bigger than the original bandwidth. Considering the gradually decreasing amplitude of extrapolated signal, a piecewise iterative extrapolation method could start from the part closer to the known part and then extend the extrapolated range. The description, simulation results, and analysis of this method would be discussed in the following sections.

Explanation of Piecewise Extrapolation Method
To explain the extrapolation process more clearly, we notate the known part and the part to be extrapolated, respectively. As for a frequency-domain bandlimited signal, it can be divided into two parts: The operator (8) can also be expressed with convolution form as In Figure 1, all MSEs decrease with more iterations, reflecting a more accurate extrapolation. With the same r f /s , the MSE is smaller when the truncated rate is larger. However, the MSE bottleneck appears during each extrapolation, especially when r f /s = 1. For all truncated rate with r f /s = 1, the bottleneck appears within 400 iterations and appears earlier with larger truncated rate. According to the theory of GP algorithm, it is attributed to the very small eigenvalue extrapolated with large iterations cannot improve the accuracy obviously though the whole extrapolation is convergent [18]. With the same truncated rate α, the best extrapolation presents with r f /s = 1, and MSE performance deteriorates with larger ratio. Especially with small α (α = 0.25), the MSEs at r f /s = 1.1 and r f /s = 1.3 turn larger by 120% compared with r f /s = 1. Besides, the MSE improvement from r f /s = 1.3 to r f /s = 1.1 is slight, reflecting the extrapolation accuracy is difficult to improve by adjusting the extrapolation filter unless it is optimal (r f /s = 1). Considering the influence of finite duration on the signal bandwidth, it can be inferred that the optimal ratio is not 1 for some other signals. In that case, an extrapolation convergence bottleneck should be solved and the extrapolation at large iteration filter bandwidth needs to be improved. For a fair comparison, the MSE reflects the difference at the extrapolated location and is normalized, defined as In terms of the amplitude variation shown in Figure 2, the difference between the original signal and extrapolated signal reduces with more iterations. While in Figure 2a, the variation between ITE = 200 and 500 differs more than that between ITE = 500 and 1000, reflecting a gradually slower extrapolation convergence. Besides, the amplitude of the extrapolated part farther from the known part turns smaller and smaller. As for Figure 2b, the extrapolated waveform is closer the original one with less r f /s (where r f /s ≥ 1). The bandwidth of the extrapolation filter causes obvious influences on the amplitude, leading huge reduction of the amplitude of extrapolated signal when r f /s > 1 compared to that when r f /s = 1. Moreover, when the location of extrapolated waveform is closer to the known part, the amplitude is closer to the original waveform.
From the above analysis and simulation results, it can be seen that the GP algorithm obtains a convergent error through infinite iterations but the iterations are restrained in the actual extrapolation process. Therefore, two main problems should be focused. One is the slow convergence (MSE bottleneck) when the ratio of proportion of the known part is small or the iterations reach some threshold. Another is the inaccurate extrapolation (low MSE) when the iteration filter bandwidth is bigger than the original bandwidth. Considering the gradually decreasing amplitude of extrapolated signal, a piecewise iterative extrapolation method could start from the part closer to the known part and then extend the extrapolated range. The description, simulation results, and analysis of this method would be discussed in the following sections.

Explanation of Piecewise Extrapolation Method
To explain the extrapolation process more clearly, we notate the known part and the part to be extrapolated, respectively. As for a frequency-domain bandlimited signal, it can be divided into two parts: The operator R T {} remains the rest part. In addition, B Ω {} is an Ω-bandwidth lowpass filter operator denoted by (8) can also be expressed with convolution form as where i ∈ Z + , f 0 = g(t) = U T { f (t)}. The output signal isf (t) = f i (t) when meeting certain stopping criterion. Operating B Ω {} once will bring with one-time FT, one lowpass filtering, and another IFT. With the maximum iterations ITE and (2N + 1) sample points in the extrapolation, the complexity from complex multiplication computation is with the order O(2·ITE·(2N + 1) log(2N + 1)).
To M-piecewise extrapolate f (t) from g(t) = U T { f (t)}, the part to be extrapolated is firstly divided into M pieces. If it is uniformly divided into the M pieces, the whole extrapolation process is decomposed into M sub-processes as shown in Figure 3. In Figure 3, m−1 f is the initial signal of the mth-piece extrapolation, and m f is supposed to be obtained at the end of the iteration at current piece. It can be found the iterations happen within the current piece, so is the error. The length of each piece is ∆T at both sides. The initial known part of each piece m−1 f turns longer; if m−1 f is extrapolated well enough, the accuracy of the mth-piece extrapolation would be improved. It is resulted from that less error accumulates within the shorter extrapolated part and the proportion of known signal increases. In practice, the signal to be extrapolated is a sampled received signal. The detailed piecewise extrapolation processes are explained as Table 1. Step 1: For the mth piece, the m−1 f is zero-tapped on the extrapolated location to obtain the initial signal { m f 0 (n)} where |n| ≤ N m and N m = N + m∆N points. When m = 1, 0 f 0 (n) = g(n) |n| ≤ N.
Step 2: For the ith iteration, the initial signal f i−1 is processed with DFT operator, lowpass filter and an IDFT operator to obtain g i .
Step 3: Replace the initial signal of g i (n) with m f (n) to obtain f i (n) as Step 4: (Judgement 1) If the iteration reaches the stopping criterion, continue to step 5; if not, update i = i + 1 and return to step 2. In (10), the operator R N {} remains the rest points beyond (−N, N). The piecewise extrapolation degenerates into the discrete GP algorithm for M = 1.
The object of the above extrapolation method is the frequency domain bandlimited signal. Admittedly, the frequency bandlimited characteristic can be widely found in communication signals, but some signals show better energy accumulation characteristic in other domains. For instance, OCDM signals [25], with the orthogonal basis of chirp signals, is bandlimited in fractional Fourier domain.
As for a p-order FrFT bandlimited signal f (t), a finite signal g(t) = U T { f (t)} is observed. The piecewise extrapolation iteration could also be concluded with 4 steps, where the initial signal should transform to FrFT domain, then lowpass filtering, inversely transforming to time domain (IFrFT), and replacing by g(t) within [−T, T]. A new formulation to describe the ith iteration of mth piece is shown as where i ∈ Z + . The operator B  (12) where F p (u) is the p-order FrFT result of f (t) as calculated in (1), and the kernel function K(p; t, u) is presented in (2).
When g(t) is sampled, the piecewise extrapolation follows the discrete procedure of (11)- (12). Similar to Table 1, the transform in step 2 should be DFrFT and IDFrFT accordingly. The implementation of the discrete transform can be found in [31].

The Minimum Energy Least Squares Error Solution of the Extrapolation
After sampled with the period t s , the original signal f is a vector of N 0 × 1 where N 0 = T 0 /t s and g is a vector of N × 1 where N = T/t s . Define a truncator operator U with the size of N × N 0 , the (n, u n ) elements is 1 where n = 1, 2, . . . , N and the rest are zero. The known finite vector g can be expressed as g = U f . If f is well frequency bandlimited, it follows that where B represent a lowpass filter operator as B Ω {}. The minimum energy least squares error extrapolation resultf follows that argmin g − U Bf 2 [32]. As described in Table 1, the extrapolated result is 1 f actually for the first sub-process starting from The minimum energy least squares error extrapolation is considered as Since the discrete bandlimited signal f and the known part g can be expanded with discrete prolate spheroidal basis as [16], the result of (14) is expected to be . However, the matrix B N will become ill conditioned and hard to compute the inverse matrix when N increases [31]. The iterative process shown in Table 1 cannot directly obtain the minimum energy least square error extrapolation result, but its error of each piece tends to zero when iterations tend to infinity.

Stopping Criterion and Complexity
The stopping criterion of an iterative extrapolation can be set according to the maximum iterations or MSE threshold, and can be modified with system requirements. It is noted that the accuracy of the former piece greatly affects the latter piece. Considering frequency shift will influence the bandwidth σ of original signal f (t), the bandwidth of extrapolation filter is always set as Ω ≥ σ. As a result, Besides, the complexity of the extrapolation is restrained within a certain value in practice.
Suppose each piece iterates for ITE times, the complexity from complex multiplication When the extrapolation samples are divided into M pieces uniformly, N m increases uniformly from N 0 to N at the interval of ∆N = N−N 0 M , hence obtaining N m = N 0 + m∆N. More intuitively, the multiple of computation with different truncated rates α and pieces M according to (15) are listed in Table 2, respectively. The multiple in Table 2 is the ratio between the computation of M-piecewise extrapolation and 1-piecewise extrapolation (GP algorithm) with the same maximum iterations for each piece. It can be seen that the multiple value is lower than the number of pieces M. The smaller truncated rates α is, the lower the multiple value is.
In addition to the stopping criterion of maximum iterations, MSE threshold can be also considered. The error of each piece will tend to zero when iterations tend to infinity whereas the infinite extrapolation cannot be implemented actually, which leads residual error for each piece. It should be admitted that dividing the extrapolated location into pieces will accumulate residual errors when the pieces ahead are not extrapolated to some accuracy. However, the error could be corrected when reaching some MSE threshold, which is attributed to the residual error is correlated to the known part g [13].

Bandwidth of Iteration Filter
It is noted that the length of the signal and the lowpass filter should be calculated for different pieces.
To extrapolate the mth piece, the initial signal has (2N m−1 + 1) points and the result has (2N m + 1) points where N m−1 = N + (m − 1)∆N and N m = N + m∆N. Hence the initial signal is firstly zero-tapped to N m points at both sides, which is equivalent to N m t s , and the frequency space between two points in the frequency domain should be 1/(N m t s ). The length of the extrapolation filter with the bandwidth Ω should be corresponding to N Ω = ΩN m t s where is a ceil operator. Therefore, in the ith iteration of mth piece as shown in Table 1, the step 2 can also be expressed as where the f i−1 is a vector of (2N m + 1) × 1; the lowpass matrix B in (13) is equivalent to F −1 H N Ω F; F and F −1 are DFT matrix of (2N m + 1) × (2N m + 1) and its inverse matrix (IDFT matrix), respectively; and H N Ω denotes a lowpass filter matrix for spectral sample as Since the known part can be considered as g = U f , equivalent to a convolution in the frequency domain. The spectral leakage caused by convolution with discrete sinc signal will lead that the discrete bandwidth of original signal is longer than N σ = σN m t s . Moreover, the iteration filter bandwidth should consider the periodicity of the known part and the original signal.

Results and Discussion
Numerical examples are used to illustrate the performance of piecewise extrapolation method when the bandwidth of iteration filter is larger than σ and the truncated rate α is small. The original signals are selected from typical signal forms in communication field. Table 3 lists the explanation and the value of variables regarding GP algorithm or proposed method in the simulation. The number of points of truncated signal from original signal.
The number of points of the aim signal in mth piece extrapolation. N m M The number of pieces for piecewise extrapolation method. 1~6 r f /s The ratio between the bandwidth of iteration filter and original signal 1.0~1.5 ITE The maximum iterations of each piece for piecewise extrapolation method. 0~4000 Noted that the comparison of performance should be conducted under the same computational complexity for fairness. The relationship of the traditional GP algorithm and M-piecewise extrapolation method can be found in Table 2. For instance, the MSE value of GP at ITE = 1000 when α = 0.25 should be compared to that of 2-piecewise extrapolation method at ITE = 1000/1.59, 3-piecewise extrapolation method at ITE = 1000/2.18, and the rest can be carried out in the same manner.
A typical frequency bandlimited signal f (t) = sin c(t), |t| ≤ 4 is sampled for N = 1024 points. Figure 3 compares MSEs of different r f /s = 1 and different M at truncated rate α = 0.25, calculated by (6). The decreasing MSE with gradually slowing rate reflects the convergence of the extrapolation method in Figure 4. It reflects that piecewise extrapolation can improve the extrapolation accuracy of original GP algorithm at the same computation but that too many pieces may not lead further performance improvements.
In the cases of larger bandwidth of iteration filter as shown in Figure 4b,c, the accuracy of extrapolation turns worse but the improvements from M-piecewise method differ more apparently. Suppose the maximum frequency shift is 50%, it is reasonable to set the bandwidth of the iteration filter to 1.5 times. At r f /s = 1.5, the MSE at M = 2, 3, 4, 6 decreases by 5.3%, 11.1%, 17.8%, and 20.6%, respectively, compared to M = 1. Besides, the mild variation exists at different ITEs which are 2000 for r f /s = 1, 3000 for r f /s = 1.3 and more than 4000 for r f /s = 1.5. It illustrates that the M-piece extrapolation enables the MSE to descend longer and improve the accuracy of original GP algorithm.
Cosine signals, sine signals and their sum are common signal forms in wireless communication field. They are all periodic, which is different from the sinc signals above. In the compressive OFDM transmission system, the original OFDM signal can be truncated up to 0.5 to ignore the influence of self-interference on BER [15]. Therefore, the truncated rate α of those signals are chosen as 0.5 to simulate the MSE performance of piecewise extrapolation for compressed OFDM signals in Figure 5. The sum of cosine signals and cosine/sine signals belong to one-way mapping compressed OFDM signals and two-way mapping compressed OFDM signals [8], respectively, in Figure 5a,b. The original signal f is sampled for N = 1024 points. Note that the r f /s 1.1 and 1.3 represents that the maximum carrier frequency offset is 0.1 and 0.3 respectively, which is a little different from the explanation in Table 3.
The improvements of piecewise extrapolation for compressed OFDM signal on MSE is can be found in Figure 5, especially when the bandwidth of iteration filter is larger than the original bandwidth. For piecewise extrapolation with larger M, the mild-variation area of MSE appears later with lower MSE at the same r f /s and computational complexity. From Figure 5a,b, more complex mapping form leads higher MSE but not double at the same r f /s and M. Similarly, in Figure 4b, the MSE at 'ratio = 1.1, M = 4' is lower than that at 'r f /s = 1, M = 1', and so it is at 'ratio = 1.3, M = 4'. It can be seen that the truncation error caused by large iteration filter bandwidth can be suppressed by piecewise method.
To verify the feasibility of the proposed method for FrFT bandlimited signal, an LFM signal f = exp j2πω 0 t + jπkt 2 is selected since energy concentration characteristics of its FrFT spectrum at the fraction angle ϕ = −arccotk. The original signal is with ω 0 = 8, k = 1, |t| < 4 and is sampled for N = 1024 points. Considering that the discrete FrFT spectrum of LFM signal is enveloped in sinc shape, the original signal bandwidth is selected as its first zero-crossing point which is a little larger than the computational result from ϕ. Combined with the preliminary of FrFT and the piecewise extrapolation method, the MSEs at different r f /s and M's are illustrated in Figure 6. As shown in Figure 6, the piecewise extrapolation method enables to improve the accuracy of GP based extrapolation method (M = 1) for the OCDM signal, and the improvement is more obvious with larger r f /s . When comparing the MSE at the same r f /s but different M, the mild-variation area at different M appears at similar ITE but with different computational complexity and MSE values. For instance, the mild-variation area of MSE appears at ITE = 1000 for all M's at r f /s = 1.1, but the total iterations are different. Referred to Table 2, the total computation of M = 2 and M = 4 are 1.72 times and 3.17 times than that of M = 1, respectively. With a fixed computational complexity of ITE = 4000 and M = 1.0, the MSEs of 2-piece extrapolation and 4-piece extrapolation decrease by 60% and another 95% at r f /s = 1, 60% and another 88% at r f /s = 1.1, and 52% and another 46% at r f /s = 1.3. From the perspective of improvement degree, the proposed method behaves better at the r f /s closer to 1.

Conclusions
In this correspondence, an M-piecewise extrapolation method is proposed to improve the accuracy of a GP algorithm at the same computational complexity. The proposed method involves the extrapolation for Fourier domain and fractional Fourier domain bandlimited signals, of which the classical GP algorithm is a special case. Simulation results show the MSE superiority of the proposed method compared with GP algorithm or GPbased extrapolation algorithm for frequency bandlimited signals and FrFT bandlimited signals, respectively. The improvement from proposed method with 2~4 pieces is obvious in general with different bandwidth of iteration filter and truncated rates. When the GP algorithm reaches the MSE bottleneck, the MSE of proposed method is basically more than 5% lower than GP and can be improved further. When the proposed method reaches the MSE bottleneck, the MSE is more than 20% lower than GP at the same computational complexity. Moreover, the influence of the factors including number of pieces, the bandwidth of iteration filter and truncated rate on MSE and computational complexity is analyzed and quantified. When the iteration filter bandwidth is larger or the truncated rate is relatively small, the extrapolation accuracy can be lifted through increasing M under the same computational complexity (but too many pieces may not lead further performance improvement).
Author Contributions: Conceptualization, Y.Z. and X.S.; methodology, Y.Z.; writing-original draft preparation, Y.Z.; writing-review and editing, Y.Z., X.F., X.S. and X.L. All authors have read and agreed to the published version of the manuscript.