Abstract
A channel simulator is an essential component in the development and accurate performance evaluation of wireless systems. A key technique for producing statistically accurate fading variates is to shape the flat spectrum of Gaussian variates using digital filters. This paper addresses various challenges when designing real and complex spectrum shaping filters with quantized coefficients for efficient realization of both isotropic and nonisotropic fading channels. An iterative algorithm for designing stable complex infinite impulse response (IIR) filters with fixedpoint coefficients is presented. The performance of the proposed filter design algorithm is verified with 16bit fixedpoint simulations of two example fading filters.
1. Introduction
Wireless communication systems must be designed to operate over radio channels in a wide variety of expected environments. Testing wireless transceivers is challenging due to unrepeatable and uncontrollable channel conditions. The initial performance verification of communication systems at the early stages of the design cycle is performed based on the channel characteristics defined by the underlying wireless communication standard. Therefore, accurate emulation of fading channels is a key step in the design and verification of wireless communication systems.
In the multipath propagation scenario, the received signal contains different faded copies of the transmitted signal. The effect of the multipath fading channel on the baseband signal can be modeled with a timevariant linear system with the following impulse response [1]: where is the number of independent paths, is the average attenuation of the th path, and and denote the complex gain and delay of the th path. Each path gain is commonly modeled as a complex Gaussian widesense stationary (WSS) process [2].
To simulate a fading channel, we need to generate a suitable sequence of complex path gains and then superimpose delayed replicas of the transmitted samples with the given delays and path attenuations . Two major techniques have been widely used for simulating fading channels. In the first approach, the socalled sumofsinusoids (SoS) model, the fading process is modeled as the superposition of a sufficiently large number of sinusoidal waves. This approach was originally proposed by Clarke [3] and later simplified by Jakes [4]. Over the past four decades several modified SoSbased models have been proposed (e.g., [5, 6]).
The second approach, which is used in this paper, is called the filterbased scheme. In this approach, to generate the inphase and quadrature components of complex fading process with a particular correlation between variates, a complex zeromean and unitvariance Gaussian random process with independent samples is passed through a spectrum shaping filter (SSF) [7] with an appropriate frequency response . A linear filtering operation on the complex Gaussian samples with flat power spectral density (PSD) yields samples that also have a Gaussian distribution, with spectrum , where is the PSD of the input samples and is the PSD of the output samples.
Compared to the SoSbased method, the filterbased fading simulation method is much trickier to design and implement. A filterbased simulator needs to be designed carefully because of possible instability and possible finite wordlength effects when implemented with fixedpoint arithmetic. On the other hand, filterbased fading simulation has several advantages over the SoSbased method. First, with the filterbased method, it is possible to simulate a wide range of power spectral densities. Second, the filters can be designed to provide a high level of statistical accuracy. Finally, the generated samples have very accurate Gaussian distribution.
In this paper, we present various considerations for designing stable fixedpoint spectrum shaping filters for modeling both isotropic and nonisotropic fading channels. We propose an iterative filter design technique. The leastsquares cost function is suggested in polar coordinates and augmented with a barrier function that keeps the poles (and potentially also the zeros) within the unit circle to enforce filter stability.
The rest of this paper is organized as follows. Section 2 discusses various considerations when designing isotropic fading channels. Section 3 presents filter design techniques for isotropic fading channels. Section 4 discusses the modeling of nonisotropic fading channels. In Section 5, we present our stable real and complex filter design techniques for stable fixedpoint implementations. Finally, Section 6 makes some concluding remarks.
2. Isotropic Fading Channels
Isotropic scattering refers to the case in which the distribution of the incident directions of the received multipath signals, or angle of arrival (AoA), are equally distributed. Assuming an isotropic scattering Rayleigh fading channel with an omnidirectional antenna at the receiver, the path gains are modeled using a unitvariance zeromean complex Gaussian process [2]. Both and are Gaussiandistributed independent stochastic processes that have the same autocorrelation function (ACF) , where is the maximum Doppler frequency and is the zerothorder Bessel function of the first kind [4]. The PSD associated with either the inphase or quadrature component of a complex fading signal has the well known Jakes' Ushaped bandlimited form with independent inphase and quadrature samples [8]: Fading samples can be generated by passing a stream of independent Gaussian samples through a SSF with the magnitude response equal to the square root of the magnitude of the PSD of the desired fading process (i.e., ).
The filtering process can be carried out in the time or frequency domains. In the frequency domain, the Gaussian samples are multiplied by . Then, an inverse fast Fourier transform (IFFT) can be applied to the resulting discrete spectrum to obtain time series fading samples [9]. The resulting series is still Gaussian by virtue of the linearity of the IFFT, and it has the desired Jakes' spectrum. One major disadvantage of the IFFT method is its blockoriented nature, which requires all fading coefficients to be generated and stored before the data is sent through the channel. This implies significant memory requirements and precludes unbounded continuous transmission, which is usually preferred in long running characterization applications such as hardware fading simulation.
In the time domain, the SSF can be implemented with either finiteimpulse response (FIR) filters [10, 11] or infinite impulse response (IIR) filters [12, 13]. When implementing fading channel simulators using FIR and IIR filters, it is important to note that the degree of the FIR filter is related to the time span of the truncated signal held in the filter and inversely proportional to the Doppler frequency. Specifically, implementation of an extremely narrowband digital filter with a sharp cutoff and a very large attenuation in the stopband requires a highorder FIR filter. Meeting the same specifications with an IIR filter typically requires fewer hardware resources than an FIR filter. In fact, utilizing both feedforward and feedback polynomials in an IIR filter permits steeper frequency rolloffs to be implemented for a given filter order than an FIR filter [14]. However, an FIR filter has no feedback and is thus inherently stable. As the coefficients are quantized in any fixedpoint implementation, the resulting numerical error is fed back in the IIR filter, possibly causing instability. Moreover, such effects can cause significant deviations from the expected response. Thus we must ensure that the designed SSF is stable.
The autoregressive (AR) modeling approach has also been proposed for generating fading processes by passing the noise through an allpole IIR filter [15]. To produce samples with accurate statistics, the AR model needs a large filter order, which greatly increases the number of required multiplications per output sample. Also, implementation of the AR fading simulator demands highly accurate floatingpoint variables, which makes it unappealing for compact fixedpoint implementations.
3. Spectrum Shaping Filter Design for Isotropic Fading Channels
In the design of the SSF, IIR filters are widely preferred to FIR filters due to their reduced computational complexity and their higher stopband attenuation. Moreover, since no constraints are imposed on the phase response, we choose to use IIR filters. We approximate the desired magnitude response of the SSF with an IIR filter of order . The magnitude response of the SSF is expressed as which is equivalent to the magnitude response of cascaded canonic secondorder sections (SOSs), also known as biquads. In (3), , , , and denote the coefficients and denotes the realvalued scaling factor of the th biquad.
For a typical wireless communication scenario, the Doppler frequency is significantly lower than the signal sample rate . Thus, SSF would have an extremely narrow bandwidth and a very sharp cutoff. We can reduce the complexity of SSF by designing it at a much lower sampling frequency, thereby reducing the required computations and also improving the accuracy of the designed filter. The resulting lowrate signal can then be interpolated to reach the target sampling frequency , where denotes the number of interpolation stages and is the interpolation factor at the th interpolation stage. Given the desired frequency response , we can find the IIR filter coefficients using the MATLAB function . This function uses doubleprecision floatingpoint variables and calculates the optimal values by minimizing the norm [16].
Due to the higher hardware cost and complexity of floatingpoint hardware, fixedpoint arithmetic is often preferred in verylargescale integration (VLSI) and fieldprogrammable gate array (FPGA) implementations. After the initial filter design, the filter coefficients and scaling factors are then quantized for fixedpoint implementation. For a compact hardware implementation, variables should be implemented with the minimum possible fixedpoint wordlength. However, reducing the wordlength impacts the response, and potentially the stability, of the designed IIR filter. When the filter coefficients are quantized from floatingpoint to fixedpoint, the poles and zeros of the system function typically shift to new positions in the plane. Unfortunately, this step can perturb the implemented frequency response from its intended response. If the designed IIR filter is extremely sensitive to coefficient changes, the filter response might not meet the target specifications or the filter might even become unstable [17]. At this step, the filter has to be tested with fixedpoint variables, and scaling factors are determined that sufficiently limit the magnitude of the generated samples to keep them within the representable range. To make sure that the filters are stable under quantization effects, we have designed the filters in fixedpoint format using Filter Design Toolbox [18], which offers bittrue implementations of SOSs with section scaling and reordering to obtain maximum accuracy.
An important observation is that if the stopband attenuation of the shaping filter is not sufficiently high, then the outofband noise that passes through the filter will degrade the accuracy of the statistics of the generated fading variates. Specifically, since designing a narrowband filter with a sharp cutoff and large attenuation invariably leads to a highorder filter, to obtain the closest approximation to the desired frequency response with a relatively small filter order, we only minimized the approximation error in the passband of the SSF. The lowpass filters utilized downstream in the interpolator stages can then be designed with extra attenuation over the transition region to ensure a sharp cutoff.
To generate the Rayleigh fading process, independent samples of a zeromean complex Gaussian process, generated by the Gaussian noise generator (GNG) block in Figure 1, pass through an SSF with a magnitude response equal to the square root of the magnitude of (2). Samples generated by the shaping filter at a low sampling frequency need to be oversampled and passed through lowpass filters in order to obtain the target sample rate . Let denote the integer upsampling factor. Interpolation of the signal is performed by inserting zeros between each pair of successive samples of and then passing the resulting stream through a lowpass filter with cutoff frequency radian/sec. In order to reduce the computational complexity, we perform the interpolation in multiple stages. Since only the amplitude response affects the correlation properties and no restrictions are imposed on the phase response, we use an elliptic IIR lowpass filter (EILPF) in the first stage. The lowpass filter has a symmetric frequency response and hence its poles and zeros appear in complex conjugate pairs, therefore, this filter can be realized using cascaded biquads. The SSF output is upsampled times and passed to the EILPF generating samples per second. The samples are further upsampled and interpolated using fadingspecific interpolation lowpass filters (SILPFs) to obtain the desired output sample rate . Since the maximum Doppler frequency is typically much smaller than the sampling frequency, we use a multiplicationfree SILPF with frequency response: where is the number of cascaded stages.
We propose to design the SSF at a sampling frequency , where . Choosing in this range satisfies the minimum Nyquist rate while keeping the computational complexity low. In addition, we have the opportunity to exploit powerof2 interpolation factors to further reduce the hardware complexity and simplify the filter design. The generated samples from the SSF are upsampled times and passed through the EILPF. Note that since the SILPF stages are designed to operate on narrowband signals, the first interpolation stage is positioned prior to the SILPF stages. Then, the samples are passed through successive SILPFs. The SILPF interpolates the signal times. Based on the processing architecture, the relation between and the target output sampling rate is . From here we have , where is an integer value in the range . Based on the maximum interpolation factor , where , each SILPF is assigned a specific interpolation factor. The minimum Doppler frequency that can be simulated by this system is
To ensure accuracy, we measured the ACF and the probability density function (PDF) of the generated fixedpoint fading samples against the ideal functions. Figure 2 shows the autocorrelation of the real part of the generated fading process. The Doppler frequency is set to Hz and the sample rate is million samples per second. This figure confirms a close match between the desired response and the generated results over up to 60 seconds. In another example, we measured the PDF for the amplitude of the generated samples with Hz and MHz. The plots in Figure 3 show that the measured PDF accurately matches the ideal Rayleigh PDF.
4. Modeling Nonisotropic Fading Channels
Isotropic scattering assumption has been challenged [19] due to the blockage of some propagation directions and antenna directivity, resulting in a nonuniform PDF for AoA at the receiver. Several nonuniform PDFs have been proposed in the literature to represent the AoA including the geometrically based PDFs [20], Gaussian PDF [21], quadratic PDF [22], Laplace PDF [23], cosine PDF [24], and von Mises PDF [25]. The von Mises PDF, which includes the uniform AoA distribution as a special case, is supported with empirical measurements of narrowband fading channels in [25]. Also it is argued that the von Mises PDF is attractive because it can approximate other nonuniform PDFs and can provide mathematical convenience for analysis [25].
When the scattering encountered in the propagation environment is nonisotropic [26, 27], the complex envelope of the fading process is where , are independent and identically distributed (i.i.d.) angles of arrival of the incoming wave at the receiver antenna with distribution , , are i.i.d. phases with uniform distribution over , and , are deterministic normalized complex constants that satisfy . The PSD function associated with is given by [25] where controls the beam width, denotes the average AoA of the scattered component, and is the th order modified Bessel function of the first kind. To obtain (7), it is assumed that the AoA of the scattered component is distributed with the von Mises/Tikhonov distribution [28] as follows: Note that when the beamwidth parameter is zero, the AoA has a uniform distribution over , and (7) reduces to Jakes’ “Ushaped” spectrum .
One of the important steps for generating the fading samples is designing the SSF. While most of the previously proposed filter design techniques might be appropriate for simulations of isotropic fading, they are not the best candidates for nonisotropic scattering scenarios, where AoA is not uniformly distributed or when the antennas are not omnidirectional. In contrast to isotropic fading, the PSD of fading samples in nonisotropic scattering (7) is not symmetric in general, and hence the filter coefficients are potentially complex [29] (e.g., the TGn fading channel models for simulating IEEE 802.11 radio propagation [30]). We will refer to filters with real coefficients as “real filters” and filters with complex coefficients “complex filters.”
Various methods such as nonlinear optimization [31], linear programming [32], and semidefinite programming [33] have been suggested for designing IIR digital filters in the complex Chebyshev sense. Leastsquares methods have also been applied extensively to design IIR filters [34]. For filter stabilization, several approaches have also been proposed. One method, proposed in [35], is to start the optimization from a stable point and then control the step size so that the solution never leaves the stable region. However, this method is computationally expensive and not easy to implement with traditional optimization procedures. In a second method, explicit constraints are imposed on the coefficients of the denominator of the transfer function [36]. This technique, however, has some limitations that affect the filter quality [37]. Finally, in a third method, the least squares cost function is modified so that the minimum always falls in the stable region [37]. This is ensured by adding a barrier function to the original cost function to avoid filter instability. To design the barrier function, an allpole proxy transfer function is formed consisting of all of the filter poles. The barrier function is basically the sum of the squared amplitude of a section of the impulse response of the proxy transfer function. If the filter is unstable, the trailing tail of the response will have (large) nonzero values.
For efficient filter implementation, the filter design method needs to produce filters that can be implemented with the minimum wordlength while providing the required accuracy. The bitprecision selected for implementing the SSF plays an important role in the accuracy, stability, and efficiency of the filter implementation. Various methods have been proposed for designing IIR filters with fixedpoint coefficients [38–41]. In [38], the optimal fixedpoint realvalued filter is found by reformulating a cost function to include the hardware complexity. The cost function is then minimized using the simulated annealing method [42]. In this approach, no constraints are imposed on the precision of the intermediate signals, which can result in inefficient implementations. The filter design method proposed in [39] is based on FIRtoIIR filter conversion. However, the resulting IIR filter has significantly higher complexity than its FIR prototype. In [40], IIR filter design is framed as a combinatorial optimization problem that is solved using a genetic algorithm. The design approach in [41] proposes bitflipping and the downhillsimplex method for fixedpoint IIR filter design.
5. Spectrum Shaping Filter Design for Nonisotropic Fading Channels
We represent the IIR filter as a product of firstorder sections (FOSs) [43]: where is a positive scaling factor, and are the th complex zero and pole, respectively, and is the number of FOSs, that is, the filter order. Here, we focus on designing IIR filters with a prescribed amplitude response. When the amplitude of the frequency response is symmetric, the poles and zeros of (9) appear as complex conjugate pairs and the IIR filter can be implemented using biquads.
We assume that the desired amplitude response is represented with samples , where and where is the desired response, is the normalized sampling frequency, and is the attenuation in the stopband. Similar to the work in [44], we express , where represents the product of FOSs in (9) and the column vector of length contains , , , and for . Next, to find the filter parameters we define the cost function as follows: where the weight vector allows us to emphasize the error minimization for certain frequency bands. Note that the sum of squared errors on a logarithmic scale is augmented by a parametric barrier function . Function is included to keep the poles (and zeros, if necessary) within the unit circle to enforce filter stability and is defined as where In (12), determines how fast the barrier function grows outside of the unit circle and the parameter determines an outer boundary for the poles and zeros. When , the barrier function tries to keep the poles within a circle of radius . Setting , on the other hand, forces both the poles and zeros into the same boundary. The barrier function (12) is especially useful when designing filters for fixedpoint implementation since it can be parameterized to keep the poles and zeros at any desired safe distance from the unit circle. Moreover, using this technique, the quantization noise can be reduced to acceptable levels. It can be shown that the variance of the quantization noise that originates in the th FOS when implemented in directformI (DFI) is where is the number of bits used to quantize the coefficients (and the intermediate variables) at the th stage. To derive (14), it is assumed that the quantization noise after each multiplier is uniformly distributed, widesense stationary white noise that is uncorrelated with both the input signal and the quantization noise in other stages. We also assume that the samples are truncated and represented in 2'scomplement. Controlling the maximum absolute value of the filter poles limits the quantization effects to an acceptable level.
The coefficients of the IIR filter are calculated using an iterative optimization algorithm. At each iteration, the optimum scaling factor is calculated as This expression for is found by differentiating (11) with respect to and setting the resulting expression to zero. To calculate the gradient vector the partial derivative of with respect to the th component of can be calculated as where , and each partial derivative of the barrier function is given by
Algorithm 1 summarizes the steps for our proposed iterative filter design. We utilized the ellipsoid algorithm [45] here for its simplicity; however, other optimization techniques could be used. This filter design algorithm can be parameterized to provide a close approximation of the desired response. The desired response , the fixedpoint format used for the filter coefficients, a weight vector , and the outer boundary for zeros and poles are passed to Algorithm 1. The filter design procedure can start with a reasonable order for the initial approximation. The filter order can be increased gradually if the desired filter characteristics are not met. The algorithm starts from an arbitrary point contained within the unit sphere and a relatively large ( times in this algorithm) initial ellipsoid matrix , where denotes the identity matrix. The algorithm then searches for the optimal solution within the present ellipsoid of feasible points. This algorithm then converges on the optimal solution by successively reducing the size of the ellipsoid by until it is small enough (i.e., the algorithm has converged) or when is within a chosen accuracy . The function represents the quantization effects that affect each element of in the Cartesian coordinate system (coefficients are transferred to Cartesian coordinates, quantized, and then transferred back to the polar coordinates). Note that stable real IIR filters can be designed with the above algorithm as well. To design such filters, the sample update is only performed for half of the poles and zeros, and the others are simply the complex conjugates of the updated samples.

An important point to note is that IIR filters are naturally susceptible to arithmetic overflow and instability due to the inherent feedback. Both the design and implementation of digital IIR filters must be carried out carefully to avoid such pitfalls. While scaling techniques are commonly used to keep the filter variables in range [46], a poor choice of scaling factor can result in a loss of signal precision and an increase in quantization noise. Another technique is to use more bits to represent intermediate signals. This method, however, is inefficient on DSPs with a fixed wordlength. Moreover, simply adding extra bits increases resource utilization in FPGA and ASIC implementations. We suggest to use (i) polezero ordering and/or (ii) augmenting auxiliary poles and zeros techniques for minimizing the range of intermediate signals that can be effectively used for reducing the signal range, overflow probability, and resource utilization in hardware implementations [47].
Considering a small section of an IIR filter, overflows are mainly caused by oscillations around resonant frequencies. Assuming a limited input signal range, we can reduce the signal range by reducing the maximum oscillation magnitude. The oscillation frequency is mainly determined by filter poles. Consider a single FOS of an IIR filter with only one pole at frequency . The output range of this section can be reduced significantly if the input signal to this section is attenuated around frequency . In polezero ordering technique, we need to implement the IIR filter with the DFI structure so that the input signal experiences a zero before it is affected by the pole. We thus sort the filter sections according to their pole magnitude in an ascending order and match the filter poles with the closest zeros. Further, signals are scaled in the different stages to assure filter stability.
The second method for reducing the signal range is to attenuate the input signal around the resonant frequencies of an IIR filter , that could potentially result in oscillation or overflow. This imposed distortion can be later compensated for with additional poles and/or zeros. This technique is applicable only if the input signal does not have a major frequency component around . One example is the implementation of narrowband lowpass filters with an approximate resonance frequency around DC. If the input signal does not have a DC component, it can be first passed through a highpass filter (difference) prior to being passed through the filter. The filter output can then be later compensated by passing the output signal through the integrator . The coefficient is intentionally added here since quantization noise and computational errors can render a perfect integrator (i.e., ) unstable.
When employed along with polezero ordering, the augmented poles and zeros are not included in the ordering process, and instead they keep their position in the DFI structure, that is, the augmented zero appears first and the pole appears last. This technique, when used in conjunction with polezero ordering and scaling, can provide efficient, accurate, and compact implementation of real and complex IIR filters.
To demonstrate the performance of our filter design procedure, we design two IIR filters with FOSs. Here, the SSF is approximated with where is the positive scaling factor, and are the complex zero and pole, respectively, and is the filter order. This filter can be realized as a cascade of FOSs. For isotropic scattering, however, poles and zeros of (20) appear in complex conjugate pairs and the shaping filter can be implemented using canonic secondorder sections. We set the parameters and , that is, the poles and zeros are bounded within a circle of radius . For all of the FOSs, the number of bits for representing each coefficient is set to . Figure 4 shows the frequency responses of the designed filters (with outofband attenuation ) as well as the desired responses. In this figure, filter (a) is a lowpass filter with normalized pass frequency , and filter (b) is a complex filter with frequency response for and for . As Figure 4 shows, the frequency responses of the designed fixedpoint filters closely match the desired responses. Figure 5 shows the position of the poles and zeros for the designed filters. Note that all of the poles and zeros lie within a circle of radius .
(a)
(b)
To illustrate the effectiveness of these range reduction techniques, we simulated an order Elliptic lowpass filter with sample rate Hz, pass frequency Hz, stop frequency Hz, passband ripple 1 dB, and stopband attenuation dB. We measured the maximum absolute range of variables by passing uncorrelated zeromean Gaussian samples through the designed filter. It is assumed that the input signal (white Gaussian noise) does not have a DC component. Table 1 shows the maximum absolute signal ranges for different filter implementations. The output of each FOS or SOS is scaled to lie within . As this table shows, the directformII (DFII) implementation requires the most number of bits (at least nine bits for the integer part). The DFI implementation of this filter with the polezero ordering (DFI + ord.) reduces the signal range significantly. Moreover, augmenting a zero at DC (DFI + ord. + aug.) can further reduce the signal range such that the minimum number of integer bits is reduced to five. In this example, the range reduction technique saves four bits in wordlength, which can significantly reduce the hardware complexity.
6. Conclusions
In this paper we discussed various considerations when designing accurate isotropic and nonisotropic fading channels. We investigated the filterbased approach, which can be customized to accurately provide the statistical properties required for simulating different fading scenarios. We proposed an iterative technique that can be used to design both complex and real IIR filters with fixedpoint coefficients. We augmented the leastsquares cost function with a specific barrier function to control the location of the poles (and potentially also the zeros) within the unit circle to ensure numerical stability. We also suggested two techniques for reducing the signal range. Signal range reduction is particularly important in hardware implementation since it can directly reduce hardware complexity. Simulation results showed that the proposed filter design technique can significantly reduce the required wordlength for fixedpoint IIR filter implementations.