Click here for the enhanced PDF version of this article
For a pdf version of this article, including diagrams.

In wireless transmission scenarios where a receiver is in motion relative to a transmitter with no line-of-sight path between their antennas, Rayleigh fading is often a good approximation of realistic channel conditions. The term Rayleigh fading channel refers to a multiplicative distortion h(t) of the transmitted signal s(t), as in y(t) = h(t) { s(t) + n(t), where y(t) is the received waveform and n(t) is the noise. The design of a simulator presented here mimics a sampled version of the channel waveform h(t) in a statistically accurate and computationally efficient fashion. This simulator exists as one of the channel models in the Visual System Simulator (VSS).

The channel waveform h(t) is modeled as a wide-sense-stationary complex Gaussian process with zero mean, which makes the marginal distributions of the phase and amplitude at any given time uniform and Rayleigh, respectively, hence the term Rayleigh fading. The autocorrelation properties of the random process h(t) are governed by the Doppler frequency ƒ_D, as in [10]:

where J_o(·) is the zero-order Bessel function of the first kind. This gives rise to the well-known, non-rational power spectral density (PSD) of the channel process

Even when the transmission rate is high enough to make the channel frequency-selective (i.e., more than one time-varying channel taps, unlike above), a good channel model is the WSSUS model of Bello [8]. According to this model, each channel tap is a zero-mean, complex Gaussian random process like h(t) described above, uncorrelated with — and thus also independent from — any other tap process, but having time-auto-correlation as described by Equation 1. It is important to note that even if there is a significant line-of-sight component in the channel, the only change in the above model of multiplicative distortion is that each channel coefficient (tap) now has non-zero mean. However, its correlation properties are still described by Equation 1, and its PSD is given by Equation 2 — to within an additive constant.

The problem of efficiently generating one — or in the case of frequency-selective channels, more than one — random complex Gaussian process with the above statistics for simulating a wireless channel has been approached in three main ways. The first approach, known as the sum of sinusoids was proposed by Jakes in [10] and amounts to the superposition of a number of sinusoids having equal amplitudes and random uniformly distributed phases to generate h(t). This approach was refined in [9] and extended somewhat to generate multiple, uncorrelated fading processes needed in frequency-selective channels. It was recently improved with additional randomization in [6] and [11]. This method, however, is computationally cumbersome due to the large number of expensive sin(.) function calls needed in the simulation.

The second approach to generating correlated complex Gaussians with the PSD in Equation 2 was presented in [3]. The idea was to multiply a series of independent complex Gaussian variables by a frequency mask equal to the square root of the spectral shape in Equation 2. Then the resulting sequence is zero-padded, and an inverse FFT is applied. The resulting series of variables is still Gaussian by virtue of the linearity of the IFFT, and has the desired spectrum (Equation 2), and hence also the autocorrelation of Equation 1. Computationally this approach is efficient because the heaviest effort is required by the IFFT, which only costs O(N log_εN) operations, where N is the number of time-domain sampled Rayleigh channel coefficients.

One disadvantage of the IFFT method is its block-oriented nature, requiring all channel coefficients to be generated and stored before the data are sent through the channel, which increases memory requirements and precludes continuous transmission.

The third approach, also followed here, adopts the filtering of independent Gaussian variables, easily generated in a pseudorandom Gaussian generator, by a filter with frequency response equal to the square root of Equation 2. Because the filter performs a linear operation, the resulting sequence remains Gaussian, with a spectrum S_out(ƒ) = S_in(ƒ)|H(ƒ)|², where |H(ƒ)|² is the squared magnitude response of the filter, chosen equal to Equation 2. If the input sequence is independent (flat PSD), the spectral shape of the output Gaussian sequence will follow Equation 2 as desired. Two main remaining design challenges remain.

First, note that the sampled channel waveform is band-limited to a discrete frequency of ƒ_DT, where 1/T is the sampling rate. For many common applications, this discrete frequency is very small: for instance, for a system at 900 MHz, with a vehicular speed of 50 mph, ƒ_D = 67 Hz, and at a sampling rate of, say, 1 Msample/sec, the discrete Doppler rate becomes ƒ_DT = 0.000067. This means that the filter used for spectral shaping of the sequence would have to be narrowband, with a very sharp cut-off and infinite attenuation in the stopband. While these requirements would lead to an impractically long FIR filter, they can be satisfied a lot more easily with an IIR filter combined with a nearly ideal interpolator, similar to [1]. The combination of IIR filter and polyphase interpolator employed in VSS is more accurate and computationally faster than an FIR filter.

The second main source of problems is that the square root of the spectrum in Equation 2 — the target filter magnitude response — is irrational and does not lend itself to any of the straightforward filter design methods. All-pole or autoregressive (AR) filter design [4] requires high orders to approximate the autocorrelation at large lags. Fortunately, a technique developed in a 1970 paper by Steiglitz [5] allows the design of an IIR filter with arbitrary magnitude response. This technique is more general than is needed here, and we present it briefly, adjusted to the needs of this particular problem.

Simulator description

The Rayleigh fading simulator is composed of two basic parts and a scaling factor, as shown in the block diagram of Figure 1. Sampled complex white Gaussian noise w[n] is generated and fed to a fixed IIR filter F(z). This filter is low-order and is designed such that its squared magnitude response approximates the spectrum in Equation 2 for a given discrete maximum Doppler rate of ρ = 0.2. Before the polyphase interpolator, a scaling factor A normalizes the power of the final resulting channel waveform h[k] = h(kT) to 1.0.

Following that, a polyphase interpolator ideally resamples the process such that any discrete Doppler rate desired by the user can be approximated with minimum computational effort. Clearly, the sampling rate at the output is intended to be several times larger than that of the input. The two are related by the interpolation factor I. This structure can be replicated for as many paths (i.e., complex channel taps) as the user desires, each with possibly distinct Doppler rates, to simulate a time-varying frequency-selective channel.

Strictly speaking, the above structure is not flexible enough to produce a Rayleigh fading waveform with arbitrary discrete Doppler rate ƒ_DT, since only rational fractions of the chosen fixed rate (ƒ_DT)^o = ρ = 0.2 can result from this filtering/interpolation combination. For example, while ƒ_DT = 0.05 and ƒ_DT = 0.001 are simulated exactly by interpolating I = 4 and I = 200 times the output process of the IIR filter, it is cumbersome to simulate exactly ƒ_DT = 0.0773 (interpolation ratio of 773/2000) and outright impossible for ƒ_DT = π/100, using the above technique. We argue, however, that this poses no serious limitation to a systems designer simulating the effects of fast fading on a mobile receiver, for two main reasons.

First, accuracy in the Doppler rate is not of primary importance in fading simulations testing the reliability of a wireless link. If a design works with ƒ_DT = 0.02 and ƒ_DT = 0.04, both of which can be easily and accurately simulated (with I = 10 and I = 5), then simulating ƒ_DT = π/100 becomes a mere esoteric exercise. Secondly, although the model provided in VSS implements only integer interpolating factors of ρ = (ƒ_DT)^o = 0.2 (e.g., using the FLATRAYLEIGH model in VSS2004, the simulation for ƒ_DT = 0.0773 would be indistinguishable from that for ƒ_DT = 0.06667 = 0.2/3), the routines provided for the design of the IIR filter in the next section can implement any fixed Doppler rate, and then the implementation routines can be modified accordingly. For example, the irrational ƒ_DT = π/100 can be implemented by designing an IIR filter with bandwidth π/10 = 0.31415 and then interpolating by a factor of I = 10. The software routines available at the Web site shown in [2] implement only integer interpolation factors I, taken to be

However, extension to any rational I is straightforward. In fact, VSS2004 software includes blocks that perform interpolation of the data stream by any rational number. In the FLATRAYLEIGH model, any Doppler rate resulting in ƒ_DT ≥ 0.2 is simulated as if ƒ_DT 5 0.2 (without any interpolation after the basic IIR filter), and any Doppler resulting in ƒ_DT ≤ 0.00001 is simulated as if ƒ_DT = 0.00001 (strongly correlated Rayleigh, almost constant with time).

The interpolator is implemented as a polyphase filter with a windowed sinc(.) function impulse response. We have found that about G = 7 one-sided periods of the impulse response are enough to yield good results with small computational complexity. In this polyphase implementation, and counting only real multiplications, the cost per complex output sample at the desired Doppler rate ƒ_DT = ρ/I is

where K is the number of biquads in the IIR filter. A small constant overhead must be added to the cost of Equation 3 because of a short initialization phase, where the simulator runs idle for a few samples to eliminate transient behavior in the IIR filter and the interpolator. The memory required is insignificant. For a typical wireless channel case of, say, ƒ_DT = 0.01, i.e., I = 20, and for the K = 7 biquads we implemented, the cost evaluates to 30.8 real multiplications per sample. This compares favorably against both FIR filtering (an order of 31 would be too small for any ƒ_DT) and the “sum-of-sinusoids” method in [10]. For perspective, with a block-IFFT solution, where for a block of, say, N = 10,000 samples, the total cost is 2 · (Nlog₂N) real multiplications, amounting to 26.5 real multiplications per complex output sample — even ignoring the cost of the initial Doppler mask multiplications as well as the large memory requirement. Therefore, the complexity of the proposed method is linear and actually slightly improves for decreasing ƒ_DT, as can be seen in Equation 3.

IIR filter design

To design an IIR filter approximating the square root of the magnitude frequency response of Equation 2, or, more accurately, the square root of the spectrum of Equation 2 as it is translated into discrete frequency, we follow an approach proposed by K. Steiglitz [5] in 1970. An IIR filter of order 2K, synthesized as a cascade of K second-order canonic sections (biquads), is designed. We use an iterative optimization procedure known as the ellipsoid algorithm [2] to search for the optimum real coefficients a_k, b_k, c_k, d_k, k = 1,…,K and the scaling factor A, such that the magnitude response of the filter

for z = є^јω approaches the desired magnitude response Y ^d(ω). For completeness, we summarize the algorithm in [5] for the IIR filter design below.

where we define M + 1 frequency points in the Nyquist interval Wi ε [0, 1], i = 0, 1, …,M, with W_i = i/M, and their discrete frequency counterparts ω_i = πW_i, ω_i ε [0,π]. Also define z_i = є^јωi = є^јπWi. Then, for the chosen discrete cuttoff frequency ρ = (ƒ_DT)^o = 0.2 define

as the largest frequency index not exceeding the chosen discrete Doppler rate π. According to [3], the desired response to be approximated by our filter is:

where the response for i = L was determined from the requirement that the area under the sampled spectrum be equal to the theoretical case given in Equation 2, as shown in [3].

Next, define the vector x of length 4K, containing the filter coefficients a_k, b_k, c_k, d_k, k = 1,…,K. Express H(z) = AF(z, x), where F(z) is the product of biquadratic transfer functions in Equation 4, excluding the scale factor A. The filter design problem amounts to minimizing the squared error:

Since the optimum scaling factor A^o is positive, we can differentiate Q(A; x) in Equation 5 with respect to |A| and set to zero, to obtain:

Now we have to optimize R(x) = Q(A^o, x) in 4K dimensions. If we compute the gradient each partial derivative is given by:

To evaluate Equation 7 for the frequency index i = 0,…,M and for the biquad index k = 1,…,K, we have:

At this point, the basic quantities for iterative optimization exist, with the vector x as a variable. Obviously, the filter design has been done once and off-line. An important detail is that the algorithm may converge to unstable or non-minimum phase filters, which is undesirable. If the optimization leads to such biquads, we invert the misbehaving poles or zeros to bring them back into the unit circle and restart the optimization from that point, as recommended in [5]. In this fashion the filter magnitude response remains unchanged and presumably close to the desired, and the final resulting filter is a stable and minimum-phase design. Figure 2 shows the resulting magnitude response for a filter with K = 7 biquads and specified ellipsoid accuracy of ε = 0.01. The number of frequency points was M + 1 = 501. To make discrepancies clear we plot the response in dB.

Performance evaluation

The simulator was implemented in C++ as a block inside VSS2004 in its form discussed above. It operates in block mode (with user-specified block length), with the behavior between blocks also specified by the user. In continuous mode, waveform continuity is maintained across different blocks. In independent mode, the fading waveform begins independently for each block, thereby simulating potential frequency hopping, which gives rise to independent channel realizations for every data block. To test the block we used it first in continuous mode, to get correlation properties and the level crossing rate (LCR) and check the results against theoretical predictions. Next, we used the block in block-independent mode and performed symbol-error-rate (SER) simulations with QPSK and 16-QAM, checking the simulated SER against known results.

Correlation and LCR

Some of the comparisons that can be made are summarized in Figure 3, where we plot the real and imaginary parts of the autocorrelation of the simulated complex channel waveform. The real part agrees well with the theoretical autocorrelation given by Equation 1. The imaginary part of the autocorrelation should ideally be zero everywhere, as implied by Equation 1, and indeed we observe it remains low for the simulated waveform. The autocorrelations in Figure 3, both theoretical and simulated, were computed for ƒ_DT = 0.05, but similar, good agreement is true for any other discrete Doppler rate. Hence, we conclude that the approach taken is rather successful in that with reasonable computational complexity, the simulated fading waveform captures the statistical behavior expected by theory.

Another statistical characteristic often measured in Rayleigh fading waveforms is their LCR, N_R. This is defined in [10] as the expected rate at which the magnitude of the fading waveform crosses a specified signal level R in the positive direction. Its formal definition for a continuous time fading waveform is:

is the joint density function of the amplitude of the fading and its time-derivative, also given in [10]. From the definition and for the spectrum in Equation 2, the LCR is:

where λ is the signal level being crossed normalized with the rms of the fading waveform, λ = R/R_rms. For three discrete Doppler rates ƒ_DT = 0.05, 0.01 and 0.002 (interpolation factors for our simulator structure I = 4, 20, 100, respectively), we obtain simulated results for the LCR NR. These appear in Figure 4 against the theoretical values of Equation 13.

From Figure 4 we conclude that the transfer of the level crossing notion from the continuous to the discrete (sampled in time) world is not without consequences. For the finely sampled fading waveform with the low Doppler rate ƒ_DT = 0.002, the agreement is excellent for any crossing level because given the slow variation of the fading and the density of the sampling, the representation of the continuous-time properties in the discrete-time world is accurate. For the middle Doppler rate ƒ_DT = 0.01, the result deviates somewhat from theory for low λ but is almost identical to a plot in [4] using a large AR model and IFFT to generate the Rayleigh fading waveform. For fast Doppler ƒ_DT = 0.05, however, and for the low crossing levels (more than 20 dB below the rms amplitude) the more sparse sampling of the fading waveform causes underestimation of the LCR. This is not particular to our simulator, but rather a general problem, arising due to sampling. The reason it surfaces at high Doppler rate and low λ is that the sparsely sampled waveform “misses” some of the low (and also short in duration) minima of the implied continuous waveform. A portion of the low extremes is unavoidably overlooked by the sampled fading because of their small duration in continuous time, leading to the underestimation of the LCR.

SER simulations

We used the system depicted in Figure 5 in VSS2004 to run SER simulations using the VSS block FLATRAYLEIGH as well as the VSS block AWGN, in addition to QPSK and 16-QAM transmitters and receivers, also implemented as VSS blocks. The receiver is assumed to have perfect knowledge of the channel coefficient h(t), which is used to invert the received noisy samples. This inversion is often referred to as a single-tap equalizer for flat fading.

For uncoded M-QAM in i.i.d. Rayleigh fading, the SER is given by (8.106) in [7], for perfect channel knowledge. We plot the theoretical curves for M = 4 (QPSK) and M = 16 (16-QAM) against results from simulation using VSS. It can be seen that the agreement is good, with only small discrepancies arising because of the non-independent nature of the fading, which can lead to bursts of errors for periods of contiguous low channel values.

Conclusions

A simulation structure can be used for correlated complex baseband Rayleigh fading waveforms. A design method for arbitrary IIR filters and the ellipsoid algorithm were applied to obtain a fixed IIR filter (in biquad-form for robustness). This filter is used to generate a fading waveform at a fixed Doppler rate. Subsequently, fast polyphase interpolation accommodates the variety of Doppler rates that may be desired. This fixed filter-variable interpolator system ensures adaptability, good statistical properties, and computational efficiency. Simulated results compare well with theoretically expected respective ones, which confirms the validity of the proposed simulation technique.

References

1. A. Anastasopoulos and K. M. Chugg. “An Efficient Method for Simulation of Frequency-selective Isotropic Rayleigh Fading,” Proc. of Veh. Tech. Conf., Phoenix, Ariz., May 1997, p. 539-543.

2. C. Komninakis. “A Fast and Accurate Rayleigh Fading Simulator.” IEEE Global Communications Conference, GLOBECOM 2003, GC01-8, San Francisco, Dec. 1-5, 2003, Vol. 6, p. 3306-3310.

3. D. J. Young and N. C. Beaulieu. “On the Generation of Correlated Rayleigh Random Variates by Inverse Discrete Fourier Transform,” Proc. of ICUPC - 5th International Conference on Universal Personal Communications, vol. 1, Cambridge, Mass., Sept. 29.-Oct. 2. 1996, p. 231-235.

4. K. E. Baddour and N. C. Beaulieu. “Autoregressive Models for Fading Channel Simulation,” IEEE Global Telecommunications Conference, Globecom 2001, 2:1187-1192.

5. K. Steiglitz. “Computer-aided Design of Recursive Digital Filters,” IEEE Trans. on Audio and Electroacoustics, AU- 18:123-129, June 1970. Reprinted in Digital Signal Processing, L. R. Rabiner and C. R. Rader, editors, IEEE Press, 1972, p. 143-149.

6. M. F. Pop and N. C. Beaulieu. “Design of Wide-sense Stationary Sum-of-sinusoids Fading Channel Simulators,” in Proc. Of IEEE ICC, Vol. 2, April 2002, p. 709-716.

7. M. K. Simon and M. S. Alouini. Digital Communication over Fading Channels - A Unified Approach to Performance Analysis, John Wiley and Sons, Inc., 2000.

8. P. A. Bello. “Characterization of Randomly Time-variant Linear Channels,” Trans. on Communication Systems, CS (11):360-393, December 1963.

9. P. Dent, G. E. Bottomley and T. E. Croft. “Jakes Fading Model Revisited,” Electronics Letters, 19 (13):1162-1163, June 1993.

10. W. C. Jakes, Jr. Microwave Mobile Communications, John Wiley & Sons, N.Y., 1974.

11. Y. R. Zheng and C. Xiao. “Simulation Models with Correct Statistical Properties for Rayleigh Fading Channels,” IEEE Trans. on Comm., 51(6):920-8, June 2003.