4 Recursive Linear Predictor
The main objective of this chapter is to design a linear predictor without having a priori knowledge about the correlation properties of the input signal. In the conventional linear predictor the known correlation sequence of the input data is used to determine the predictor coefficients. However, in this case the correlation sequence of the signal, which needs to be estimated, will be formed in a recursive way by using window-based techniques. Since the predictor updates its coefficients based on the recursively-tracked changes in the auto-correlation sequence, the predictor will be named a Recursive Linear Predictor (RLP). As the predictor coefficients are updated for every input sample, the recursive linear predictor is expected to perform better than the conventional linear predictor when the correlation properties of the input samples change with time. To make the RLP robust in the time-varying case, a forgetting factor - as in the WRLS [9] algorithm will be used while updating the correlation sequence. At the end of this chapter the designed predictor will be used to predict the envelope of a complex flat fading channel.
4.1 One Step Ahead Predictor
The block diagram of a one step ahead predictor is shown in Figure 4.1 where the output of the predictor, x( n) is the estimate for the present input sample x(n). The predictor ˆ output x( n) can be written in terms of the past values of the input sequence as shown below: ˆ
ˆ x (n) = −∑ a iM x(n − i) = −a x
i =1 T M
(4.1)
Here,
M − a = − a1
[
M − a2
M .... − a M
]
T
is the predictor coefficient vector and
x = [ x(n − 1) x(n − 2) L x(n − M )]T is the input vector. M is the order of the predictor.
51
The negative signs with all of the predictor coefficients are used for mathematical convenience and conform to current practice in the technical literature.
x(n)
x(n-1) z-1
− a1M
x(n-M) z-1
M − a2
z-1
M − a3
z-1
M − a M −1 M − aM
ˆ x( n)
Figure 4.1: Block Diagram of a One Step Ahead Predictor.
The general equations for a D step ahead linear predictor are given in Appendix A. The equations for the one step ahead predictor can be derived from these equations by simply setting D equal to 1 (one). The normal equation for the one step ahead predictor is:
M M 2 rxx (l ) + ∑ a k rxx (k − l ) = 0 , l = 1,2,…,L k =1
(4.2)
In matrix form, rxx (1) rxx (0) r * (1) rxx (0) xx M * M * rxx ( L − 1) rxx ( L − 2) or,
M L rxx ( M − 1) a1 rxx (1) M L rxx ( M − 2) a 2 = − rxx (2) M L M M r ( L ) M xx L rxx * ( L − M ) a M
(4.3)
Ra = b
(4.4)
In (4.2) and (4.3), M is the order of the predictor and L determines whether the predictor system normal equations are over-determined or not. If L=M, the system of equations represented by (4.3) is called the Yule-Walker equations [7]. In this case the correlation matrix R is not only Toeplitz but also Hermitian symmetric. If L>M, the system of normal equations becomes over-determined. Now the main objective of designing a linear
52
predictor is to determine a, the negative of the coefficient vector of the predictor, for a given auto-correlation function (ACF). If L = M, the Yule-Walker equations can be solved for a by inverting the correlation matrix R, and then multiplying it by b. The Yule-Walker equations can also be solved by using the Levinson algorithm, which exploits the symmetry properties of R, to solve (4.4) for a, and requires only 0.5M 2 multiplications and additions. So there is a significant reduction in the computational complexity from the usual matrix inversion using Gaussian elimination, which requires on the order of M3 multiplications. If the system is over-determined (L>M), instead of the usual matrix inverse, the pseudo-inverse of R is used to get a least squares solution to (4.4). The over-determined system generally yields better results than the symmetric system, when the input signal is noisy.
4.2 Recursive Method for Changing the Predictor Coefficients
In the absence of a priori knowledge, the auto-correlation function (ACF) of the input sequence to the predictor needs to be estimated to determine the coefficients of the linear predictor by using (4.3) or (4.4). In this section a recursive approach to estimating the ACF will be described where the auto-correlation matrix R and the auto-correlation vector b, needed for solving (4.4), will be updated for every new input sample. This recursive algorithm produces an estimate of the required a priori knowledge of the ACF. Again, for a signal with time varying statistical properties, the auto-correlation sequence needs to be estimated for the present condition. At the same time, the older past values of the estimated auto-correlation need to be discarded from the present statistics. For this purpose a windowbased technique will be used in the recursive method of estimation of the ACF. In the recursive method, the ACF at any instant is estimated by using the present data sequence and the past estimate for the ACF. Since we want the system to forget the effect of the past input samples, a window of length N is used to consider only the most recent N samples of input data. Let this sequence be xN(n). These N input samples are then used to estimate the present value of the auto-correlation sequence. An N length auto-correlation
53
(single-sided) sequence is then calculated by using xN(n). Let this auto-correlation sequence be ri(n). Θ(m) = DFT{x N (n)},
ri (n) = [ri (0) ri (1) L ri ( N − 1)]
ri ( p) = IDFT Θ(m)Θ* (m) N ,
{
}
m = 0,1, L, N − 1 p = 0,1, L , N − 1 (4.5)
The present value of the auto-correlation sequence is then estimated by using the previously estimated ACF and ri(n) in the following way. r ( n) = λ r (n − 1) + ri (n) , λ +1 0 ≤ λ ≤1
(4.6)
In (4.6) λ is the forgetting factor that has the same effect as the forgetting factor in the well-known exponentially weighted RLS algorithm. Since λ is less than or equal to 1, it will provide an exponentially decaying window effect, as shown in Figure 4.2. However the RLS kind of forgetting factor has a limitation working in the non-stationary environment. This is because it can only apply an exponential window to the time series and the exponential window not only has an infinite tail but also penalizes near present data too much (Figure 4.2).
1
λ
n -i
Number of iteration
n
Figure 4.2: RLS Window Function.
54
To improve the forgetting window function, a rectangular window of length NW will be used in addition to the exponential window to forget data in the far past completely. The algorithm equivalent to the RLS algorithm but using this kind of forgetting factor is called the Windowed Recursive Least Squares (WRLS) algorithm [9]. The window function of WRLS is shown in Figure 4.3. Due to the additional rectangular window, the WRLS algorithm performs better than RLS for non-stationary signal estimation [9]. Now if the WRLS based forgetting factor is used to estimate the present correlation sequence, update equation (4.6) needs to be changed in the following way:
r ( n) =
λ r(n − 1) + ri (n) − λ NW r (n − NW ) λ + 1 − λ NW
,
0 ≤ λ ≤1
(4.7)
In (4.7), NW is the length of the rectangular window as shown in Figure 4.3. Since the window used to collect the input data is of length N, at the beginning of the estimation process it is implicitly assumed that the correlation property of the input sequence will not be the same after N input samples. Therefore, while updating the correlation sequence r(n), the length of the rectangular window NW should reasonably be less than or equal to the length of the windowed input sequence xN(n), i.e., NW ≤ N . (4.8)
1
λ ≤ 1 n - N W+ 1 ≤ i ≤ n λ
n -i
NW
n– N W+ 1 N u m b e r of iteratio n
n
Figure 4.3: WRLS Window Function.
55
Since the rectangular window takes into account only the most recent NW samples and completely forgets the past information, a higher value of the forgetting factor λ can be used to give more weight to the immediate past samples in WRLS. The estimated auto-correlation sequence r(n) ((4.6) or (4.7)) is then used to form the correlation matrix R(n) and the correlation vector b(n) in (4.5) for the time instant n. The predictor coefficients for that time instant are then calculated by solving (4.5). The steps of the Recursive Linear Predictor algorithm, which updates the predictor coefficients for every time instant, can be summarized as follows:
Table 4.1: Algorithm for Recursive Linear Predictor.
Step 1
Initialization
x N (1) = [x(0) 0 L 0]T , N is the length of x N r(0) = L×M zero matrix
Step 2
Instantaneous Correlation Sequence
Θ(m) = DFT{x N (n)},
ri (n) = [ri (0) ri (1) L ri ( N − 1)] RLS Window: r (n) =
ri ( p) = IDFT Θ(m)Θ* (m) N ,
{
}
m = 0,1, L, N − 1 p = 0,1, L , N − 1
Step 3
Estimated Correlation Sequence
λ r (n − 1) + ri (n) , 1+ λ
0 ≤ λ ≤ 1 , or
λ r (n − 1) + ri (n) − λ NW r (n − N W ) WRLS Window: r (n) = 1+ λ rxx (1) rxx (0) r * (1) rxx (0) R = xx M M * * rxx ( L − 1) rxx ( L − 2) L rxx ( M − 1) rxx (1) L rxx ( M − 2) , b = − rxx (2) M L M r ( L) * xx L rxx ( L − M )
Step 4
Determination of the predictor coefficients
a = − R −1b, if L = M , or, a = − pseudoinverse(R ) × b, if L > M ˆ x(n + 1) = a * x N (n) Step 5 Step 6 Update Recursion x N (n + 1) = [x(n) first N − 1 elements of x N (n)] T Repeat from Step 2.
56
4.3 Prediction of the Fading Envelope
In this section, the fading envelope of a flat fading channel will be predicted by using the designed recursive linear predictor. Different ways of predicting the values of the fading envelope for a flat fading channel have been discussed [10]. If the fading envelope is sampled at a rate of more than two times fm, the maximum Doppler spread, it is possible to interpolate the value of the fading envelope at any point in between two samples [10]. This kind of prediction is possible because of the inherent correlation among the samples of the fading envelope. It has been shown [11] that the power spectrum of a flat fading envelope can be expressed by the following equation:
1 .5 f − fc S EC ( f ) = πf m 1 − f m 0 f − fc ≤ fm
2
,
(4.9)
, f − fc > fm
In (4.9), fc and fm are the carrier frequency and the maximum Doppler spread of the fading channel respectively. S EC ( f ) is the power spectrum of the fading envelope. The base-band equation can be derived from (4.9) by setting fc to zero, which yields:
1.5 , 2 S E B ( f ) = π f m − f 2 , 0
f ≤ fm
(4.10)
f > fm
Figure 4.4 shows the base-band power spectrum of a fading envelope for a maximum Doppler spread of 100 Hz. From (4.10) and from Figure 4.4, it is clear that the maximum frequency component in the fading envelope is fm. According to the Nyquist sampling theorem, the fading envelope can be represented completely by its samples as long as the sampling rate is higher than 2fm.
57
Figure 4.4: Base-Band Power Spectrum of the Fading Envelope with Doppler Spread of 100 Hz.
As the power spectral density (PSD) of the fading envelope is not white, there exists correlation between the time samples of the fading envelope. The correlation is stronger for a smaller fm, or higher fs. Therefore, if the fading envelope is sampled at a rate much higher than 2fm, the future value of the fading envelope can be predicted reasonably well from its present and past samples. A simple linear predictor (specifically a one step ahead predictor) can be used to predict the immediate next value of the fading envelope. To verify this statement to investigate the performance of the recursive linear predictor, we will generate the Rayleigh fading envelope for the non-stationary situation. The stationary fading envelope is a special case of the non-stationary envelope when all of the time-varying quantities are constants.
4.3.1 Non-Stationary Rayleigh Fading Envelope Simulator
The general block diagram of the Rayleigh fading simulator [12] is shown in Figure 4.5. While simulating this system, the main difficulty in the block diagram is the implementation of the base-band Doppler filter whose PSD is given by (4.10). The magnitude response of the filter can be derived by taking the square root of the right hand side of (4.10). The equation for the magnitude response of the Doppler filter is given below.
58
1.5 π 2 H D ( f ) = 4 fm − f 0
Baseband Gaussian Noise Source
Independent Independent Independent
2
, f ≤ fm , f > fm
(4.11)
Baseband Doppler Filter Absolute Value Operator
Fading Envelope
Baseband Gaussian Noise Source
Baseband Doppler Filter j
Figure 4.5: General Block Diagram of Rayleigh Fading Envelope Generator.
Since (4.11) results in infinity at f = fm and zero for f greater than fm, it is impossible to find a digital filter that provides a magnitude response equal to that given in (4.11). A frequency domain approach for implementing the fading generator (Clark’s model [12]) avoids the difficulties of finding the appropriate impulse response of the Doppler filter. Clark’s model also uses FFTs for efficiently calculating the fading envelope. The problem with Clark’s model is that this model is designed for processing a block of data with a constant Doppler spread fm. Therefore, Clark’s model cannot be used to simulate the fading envelope when the mobile velocity changes continuously with time. The relationship between the mobile velocity and the maximum Doppler spread of the fading envelope is given by the following equation. fm = fc v v = c λc (4.12)
In (4.12), v is the mobile velocity, c = 3 × 10 8 m/ sec is the velocity of light, and λc is the wavelength of the light for the carrier frequency fc. Any change in the mobile velocity during the observation time will result in a continuous change in the Doppler spread. Since the maximum Doppler spread determines the coefficients of the Doppler filters in Figure 4.5,
59
these filters will be time-varying if the mobile velocity changes with time. In the simulation, the Doppler filter was implemented by using a linear phase FIR filter designed using the frequency domain sampling method. The design procedure is described below. 1. Since the Doppler filter is being designed in the frequency domain, the total frequency band is defined from -fs/2 to fs/2 Hz, where fs is the sampling frequency of the filter. The magnitude response of the filter, defined over the total frequency band, is then sampled at K equally spaced frequencies. Therefore, the spacing between two adjacent frequency samples is: f ∆f = s K To make the calculation efficient, K is chosen so that log 2 K is an integer. 2. Since our objective is to design a real FIR filter, we first define K/2 samples for the positive frequencies only. The magnitude samples for the positive frequency band are defined as: H (k ) = H D (k∆f ), k = 0,1, L K −1 2 (4.14)
(4.13)
3. To make the filter realizable, the magnitude response of the filter should not be infinity at f = fm. The magnitude response at f = fm was determined by linearly interpolating the previous two frequency samples. Let, Nd be the frequency index corresponding to the Doppler frequency, i.e., N d = fm ± δ and δ < ∇f . In the ∆f
expression ± δ is used to compensate for the fact that fm may not correspond to any frequency domain sample point. H ( N d ) = 2 H ( N d − 1) − H ( N d − 2) (4.15)
4. Even though the magnitude response of the Doppler filter goes to zero immediately after f = fm, to reduce the ripples in the magnitude response of the designed filter, two
60
sample points is used to define the transition band of the filter. The values of the two transition band samples is defined as (experimentally determined for the lowest value of the pick ripple): H ( N d + 1) = 0.6 H ( N d ) H ( N d + 2) = 0.1H ( N d )
(4.16)
5. Again, to reduce the ripple in the pass-band of the filter [13], the sampled values of the magnitude response were changed in the following way: H alternate (k ) = (− 1)k H (k ), k = 0,1, L K −1 2 (4.17)
6. A linearly varying phase response was added to the magnitude samples. The phase response is zero at f = 0 (i.e., for k = 0) and is equal to -π at f = fs/2 (for, k = K/2). The frequency response of the filter for the first K/2 samples is given below:
−j 2π k K ,
H dopp (k ) = H alternate (k )e
k = 0,1, L
K −1 2
(4.18)
7. To make the phase response of the filter linear, the next half of the frequency sample array was filled up by taking the complex conjugate of the first half in reverse order, i.e.,
H dopp (
K K * + k ) = H dopp ( − k ), 2 2 K H dopp ( ) = 0, for this case. 2
k = 1,L
K −1 2
(4.19)
8. The impulse response of the desired filter, hdopp(n), is then calculated by taking the inverse DFT of H dopp (k ) . Ideally the impulse response of the filter is a real sequence. Due to the finite word length effect of the computing device, h(n) can be a complex sequence with a very small imaginary part. This imaginary part is then neglected.
61
For the non-stationary case, to get an output (fading envelope) with a smooth variation of its statistical properties a time-varying filter of length K is implemented. The coefficients of the time varying filter are the same as the coefficients of the Doppler filter hdopp(n) calculated for every instant of time with the instantaneous velocity. At instant m, let the filter coefficient vector be h(m) i.e., h(m) = [h(0) h(1) ……h(K-1)]T calculated at instant m. At that instant, let the filter state vector, or the input vector be w(m), where w (m) = [w(m) w(m − 1) L w(m − K + 1)] T (4.20)
In (4.20), {w(m)} is the complex i.i.d. zero mean Gaussian random sequence of variance 2 (Figure 4.5). The output of the time-varying filter at instant m is then calculated as follows: rcomp (m) = h T (m)w (m) (4.21)
For the required length of fading envelope, the complete sequence of {rcomp(m)} is first generated. The generated sequence is then normalized for an rms value of one. The absolute value of the resultant signal is the desired fading envelope. For the stationary case, the coefficients of the Doppler filter are calculated only once. All other procedures are similar to the non-stationary case.
4.3.2 Simulation Results
In this section simulation results will be presented to verify the performance of the recursive predictor under stationary and non-stationary conditions. Before going to those results, some intermediate results will be presented to show that the implemented Doppler filter is a good approximation for the desired non-realizable filter. The magnitude responses of the idealized (linear approximation near maximum Doppler frequency) Doppler filter and the implemented Doppler filter are shown in Figure 4.6. The parameters used to generate these two magnitude responses are: Sampling frequency of the filter, fs = 8000 Hz;
62
Mobile velocity, v = 50 mph; Carrier Frequency, fc = 900 MHz; Maximum Doppler spread,
fm = fc
v 50 × 1760 × 3 × 12 × 2.54 = 900 × 106 = 67.056Hz c 100 × 3600 × 3 × 108
Figure 4.6 shows that, except in the vicinity of fm, the maximum Doppler spread, the implemented filter has almost the same magnitude response as the idealized Doppler filter. The discrepancies near fm are due to changes (Section 4.3.1) made to make the filter realizable.
Figure 4.6: Magnitude Response of the Idealized and Simulated Doppler Filter.
4.3.2.1 Fading Envelope with Constant Velocity Since the maximum Doppler spread fm of the fading envelope is constant for a constant velocity, fm will be used as the independent variable for all of the results. A Rayleigh fading envelope with fm = 100 Hz and fs = 8 kHz was generated by using the procedure explained in Section 4.3.1. The fading envelope is shown in Figure 4.7. This generated envelope sequence is then used as the input to the recursive linear predictor. The output of the predictor when an RLS type window function is used to update the auto-correlation
63
sequence, along with the original predicted envelope, is shown in Figure 4.7. An overdetermined system was used to determine the predictor coefficients. The size of the R matrix used in the simulation was 20x7. The parameters used in the simulation are given below: Maximum Doppler spread, fm = 100 Hz; Sampling frequency, fs = 8 kHz; Length of the Predictor, M = 7; Forgetting factor, λ = 0.9.
Figure 4.7: Actual and Predicted Fading Envelope (RLS Type Window Used to Update the Correlation Sequence).
Figure 4.8 shows the actual and the predicted fading envelope when the WRLS type window function is used to update the auto-correlation sequence. The window size NW used in the simulation was 50. All other parameters used to generate Figure 4.8 were the same as those used to generate Figure 4.7. From Figure 4.7 and Figure 4.8, it is clear that both recursive linear predictor approaches can follow the input sequence. The performances derived from both approaches are compared on the basis of the prediction error, which is the difference between the original and the estimated envelope. To get a constant quantity for measuring the performance of the different approaches, the time
64
average of the squared prediction error was measured and the resultant quantity was named the Mean Squared Prediction Error (MSPE). Therefore, MSPE = 1 N
N −1 n =0
∑ e 2 ( n) =
1 N
N −1 n =0
ˆ ∑ (x(n) − x(n))2
(4.22)
Figure 4.8: Actual and Predicted Fading Envelope (WRLS Type Window Used to Update the Correlation Sequence).
The MSPE of the recursive linear predictor with two different updating windows are shown in Figure 4.92. Figure 4.9 shows that for the same Doppler spread, the MSPEs are different for different realizations of the fading envelope. Therefore, to get an estimate of MSPE, the simulation was run for 50 times with the same realizations for differently windowed estimate updates. Fifty different MSPEs were then used to calculate the ensemble averaged MSPE for both approaches. The resultant ensemble average will give us an estimate of the expected value of the MSPE and this result will be used to compare the performance of the different approaches. Figure 4.10 shows the estimated MSPE for both approaches for different values of maximum Doppler spread, fm.
2
Due to the transient effect, the first few samples of the predicted sequence are not the same as the original input samples. Therefore, the first 100 samples were not used in the calculation of the MSPE.
65
Figure 4.9: Mean Squared Prediction Error with RLS and WRLS Windows for Different Realizations (fm = 100 Hz).
From Figure 4.10 we observe that with the recursive linear predictor, the prediction error increases with an increase in fm, the maximum Doppler spread. This is because for higher fm, or in other words, for higher mobile velocity the randomness of the fading envelope increases. This increased randomness results in reduced correlation among the samples of the fading envelope and thus yields an increase in the prediction error. Figure 4.10 also shows that, for any fm, the recursive linear predictor performs better if the WRLS type window function is used for estimating the correlation sequence of the input.
Figure 4.10: Overall Mean Squared Prediction Error for Different fm.
66
4.3.2.2 Fading Envelope with Constant Acceleration The same block diagram as shown in Figure 4.5 was used to generate the Rayleigh fading envelope for constant acceleration. Since the mobile velocity is changing at every instant in time, the Doppler filter of Figure 4.5 was changed for every input sample. While generating the fading envelope for constant acceleration, the first 400 samples of the envelope were generated with constant velocity of 30 miles/hour. The acceleration was imposed after the first 400 samples. To reduce the number of samples, and consequently reduce the processing time, a lower sampling frequency (fs = 2000 Hz) was used in the simulation. The velocity profiles for different values of acceleration are shown in Figure 4.11. A portion of the fading envelope generated for an acceleration of 10 miles/sec2 is shown in Figure 4.12.
Figure 4.11: Velocity Profile of the Fading Envelope for Different Accelerations.
Figure 4.12 also shows the output of the recursive linear predictor. A WRLS type window function was used in the simulation to generate the predicted envelope shown in Figure 4.12. The simulation was performed with a window size NW = 15 and a forgetting factor λ = 0.95. The same window function was used to generate the mean squared prediction error (MSPE) for different values of acceleration. Similar to the case of constant velocity, the MSPE is different for different realizations of the fading envelope for the same velocity profile (i.e., for a constant acceleration). Therefore, ten different realizations were 67
used to estimate the MSPE value for a particular acceleration. The values of the MSPE for different accelerations are shown in Figure 4.13. Figure 4.13 also shows the MSPE vs. acceleration curve when the RLS type window is used to update the auto-correlation sequence of the input signal to the predictor. The value of the forgetting factor used in the simulation for the RLS type window was 0.9.
Figure 4.12: Rayleigh Fading Envelope for Constant Acceleration of 10 miles/sec2 and the Predicted Envelope. WRLS Type Window Used in the Predictor. Sampling Frequency, fs = 2000 Hz.
Figure 4.13: Mean Squared Prediction Error for Different Accelerations.
68
From Figure 4.13 it can be said that, unlike in the constant velocity case, the performance of the recursive linear predictor (measured by MSPE), either with the RLS window or the WRLS window, does not degrade with an increase in the acceleration. Figure 4.13 also shows that the predictor with the WRLS window performs approximately 1 dB better than the predictor with the RLS window when the mobile velocity changes with a constant acceleration. This improvement in the performance comes in a trade-off with the increased memory requirements of the WRLS window over the RLS window.
69