A Novel Simulator of Nonstationary Random MIMO Channels in Rayleigh Fading Scenarios

For simulations of nonstationary multiple-input multiple-output (MIMO) Rayleigh fading channels in time-variant scattering environments, a novel channel simulator is proposed based on the superposition of chirp signals. This new method has the advantages of low complexity and implementation simplicity as the sum of sinusoids (SOS) method. In order to reproduce realistic time varying statistics for dynamic channels, an efficient parameter computation method is also proposed for updating the frequency parameters of employed chirp signals. Simulation results indicate that the proposed simulator is effective in generating nonstationary MIMO channels with close approximation of the time-variant statistical characteristics in accordance with the expected theoretical counterparts.


Introduction
Multiple-input multiple-output (MIMO) technologies are widely used to improve the channel capacity and frequency utilization ratio without increasing system bandwidth or transmitting power.The traditional MIMO channels are usually modeled as wide-sense stationary (WSS) Rayleigh fading processes [1,2].However, the WSS assumption is not always satisfied, because the time-variant scattering environments, such as vehicle-to-vehicle (V2V) channel [3] and high-speed train channel [4], will change the statistical properties of MIMO channel over time.This kind of non-WSS channel, also called no-stationary channel, plays an important role in designing, validating, and optimizing the performance of real world wireless communication systems [3][4][5][6].
Based on the modeling methodologies adopted, channel models can be classified into geometrically based stochastic models (GBSMs) and correlation-based stochastic models (CBSMs) [7].The nonstationary MIMO channels based on GBSMs have been investigated in [8][9][10].Xiao et al. [8] were concerned with the time-variant parameters such as anglesof-arrival (AoA) and angles-of-departure (AoD) irrespective of the scattering environment or clusters, while Wu et al. [9] focused on modeling the appearance and disappearance of clusters by Markov process or birth-death process.Moreover, Borhani and Pätzold [10] introduced a spatial Brownian path model to generate realistic moving trajectories of transceivers and employed it to model the time-variant channel by deriving the analytical expressions of AoA and angles-ofmotion (AOM).Researches have showed that GBSMs are suitable for modeling nonstationary MIMO channels; however, the complexity increases significantly when the number of clusters increases dramatically.
On the other hand, CBSMs are usually used to evaluate the performance of MIMO systems due to their low complexity and have been applied in MIMO channel models [11], as well as standard channel models such as IEEE 802.11TGn channel [12] and LTE-A channel [13].When adopting the conventional CBSMs in nonstationary cases, it is a prerequisite to know the variation of the statistical properties of channel with respect to the time and frequency [14,15], in order to construct a channel simulator generating nonstationary fading channel realizations.The simulation of CBSMs always requires generating multiple independent and identically distributed (i.i.d.) Rayleigh fading processes.In the past three decades, a large number of researches have showed that this simulation can be efficiently performed by a finite sum of sinusoids (SOS) or cisoids (SOC) [16][17][18][19][20][21][22][23][24][25].The mathematical reference model for this SOS method was initially proposed in [16] as Clarke's model, and the famous simulator, namely, Jakes' model, was developed in [17].Pop and Beaulieu [18] suggested an improved Jakes' model by applying random phase shifts to optimize the WSS property.Xiao and Zheng [19,20], as well as Zajić and Stüber [21], introduced some random simulation parameters to enhance the performance.
The aforementioned simulation models are all nonergodic models, which need the arithmetic average over a large amount of simulation trails (snapshots) in order to approximate the desired properties [22].To overcome this shortcoming, Pätzold and Wang [22][23][24][25] developed deterministic SOS simulators with low complexity.At the same time, the sumof-cisoids (SOC) simulators are found to be more flexible and suitable for channel simulation in the nonisotropic scattering environments [26][27][28][29].The accuracy of statistical channel properties obtained with both SOS and SOC simulators greatly depends on the simulation parameters computation method.The   -norm method [23,26] owns high accuracy of statistical properties but requires a time-consuming progress of parameter optimization.The method of equal area (MEA) and its derivatives [24,27,28] can reduce the computation complexity but demand plenty of cisoids to emulate the correlation properties precisely.Recently, the Riemann sum method (RSM) was proposed by Gutiérrez et al. in [29] to generate channels with accurate autocorrelation.However, this method fails to generate channels with desirable envelope distribution when the number of cisoids is small.
It should be noted that all above existing SOS or SOC simulators are only designed for stationary channels with constant statistical characteristics.For simulating nonstationary channels, the conventional Kronecker model is extended in this paper for reconstructing correlated nonstationary MIMO Rayleigh fading channel observed in time-variant scattering environments.In addition, we propose a new SOC simulator based on the sum of chirp signals to generate multiple independent identically distributed (i.i.d.) Rayleigh fading channels of time-variant statistical characteristics.By utilizing the linearly changing Doppler frequencies, the new method is applicable to reproducing the desirable time-variant Doppler power spectral density (DPSD), envelope probability distribution function (PDF), autocorrelation function (ACF), and cross-correlation function (CCF).
The remainder of this paper is organized as follows: in Section 2, we give a nonstationary MIMO channel model originating from CBSMs.The existing stationary Rayleigh fading channel models including the reference model and simulation model are reviewed in Section 3. The proposed nonstationary MIMO channel simulator with a new parameters computation method is presented in Section 4. Section 5 gives the performance evaluation results and analysis of the new simulator.Finally, conclusive remarks are made in Section 6.

Nonstationary MIMO Channel Model
A MIMO mobile to mobile (M2M) communication system is illustrated in Figure 1, where the solid line and dash line indicate the positions of transceiver antenna arrays at two different time instances.The different scattering environments make the MIMO channel's characters change over time.
Assuming a frequency nonselective MIMO channel system with  transmitting antennas and  receiving antennas, the MIMO channel matrix can be expressed as where ℎ , () denotes the impulse response of the subchannel between the uth ( = 1, 2, . . ., ) receiving antenna and the sth ( = 1, 2, . . ., ) transmitting antenna.Considering a Rayleigh fading channel, the envelope and phase distribution of ℎ , () obey, respectively, the Rayleigh and uniform distributions over (0, 2], while the statistical parameters such as variance may vary over time.The time-variant channel correlation matrix is defined as where vec(⋅) is the vector operator (stacking all elements of a matrix column-wise into a single vector), (⋅)  denotes the Hermitian transposition, and {⋅} represents the correlation coefficient.R H () is a  ×  matrix, and we can use it to remodel the nonstationary correlated MIMO channel as where (⋅) 1/2 denotes any matrix square root fulfilling and G is an  ×  random matrix with zero mean i.i.d.complex Gaussian process entries.Assuming the scattering environments around the receiver and the transmitter are mutually independent, the channel correlation matrix can be expressed by the Kronecker product of International Journal of Antennas and Propagation 3 the transmitter correlation matrix R Tx () and the receiver correlation matrix R Rx (), where ⊗ denotes the Kronecker product.

Multiple Stationary Rayleigh Fading Channels
The Kronecker model for a MIMO channel is related with the correlation matrix and complex Gaussian matrix.We can directly obtain the lower triangular matrix L by performing Cholesky decomposition on R H , L = R 1/2 H , and R H = LL  .If R H is nonpositive definite, the eigenvalue decomposition method [30] needs to be applied.In this section, a brief review of the SOC-based simulation method applicable to generating multiple i.i.d.stationary Rayleigh fading channels is provided, and its performance is evaluated.

The Reference Model.
Multiple Rayleigh fading processes are actually equivalent to a set of complex Gaussian random processes.Under this assumption, the kth Rayleigh fading process can be written as where  ≜ √ −1,  , (), and  , () are uncorrelated and zero mean real-valued Gaussian processes.Thus, the fading envelope   () ≜ |  ()| and phase   () ≜ arg {  ()} yield the following PDFs, respectively: where  2   denotes the variance of Rayleigh distribution.The autocorrelated Rayleigh fading process   () can be further characterized by its ACF      () ≜ { *  ()  ( + )}, while the cross-correlated Rayleigh fading processes are characterized by their CCF      () ≜ { *  ()  ( + )}, ∀ ̸ = .Here, (⋅) * is complex conjugation and {⋅} refers to the expectation of given argument.The DPSD      () describes the power distribution over Doppler frequency and can be calculated by performing the Fourier transform on the ACF; that is, Our aim is to generate  i.i.d.Rayleigh fading processes, which satisfy the following conditions: where  , () is the Doppler frequency PDF of the kth Rayleigh fading process.All elements of the antenna arrays share a common scattering environment where the smallscale characterization assumption applies.Thus, the  random processes have the same Doppler frequency PDF denoted as   ().
The Gaussian-shaped DPSD plays an important role for nonisotropic scattering environment as in the cases of aeronautical channels.Specifications of frequency-shifted Gaussian DPSD can also be found in the channel models for pan-European, terrestrial, and cellular GSM system [28].The theoretical expressions of Gaussian DPSD and ACF are given as follows: where   denotes the 3-dB-cut-off frequency and   ,   refer to the minimum and maximum frequency, respectively.For a symmetric Gaussian-shaped DPSD, the expression of ( 9) can be simplified as where   = − max and   =  max with  max =  0 ⋅V/ denoting the maximum Doppler frequency due to the movements of transceivers.Here  = 3 × 10 8 refers to the velocity of light, and  0 is the carrier frequency.

The SOC Simulation
Model.Based on the central limit theorem, a complex Gaussian random process can be modeled as the superposition of infinite properly weighted cisoids.
It is impossible to realize the simulation with infinite numbers of cisoids.Fortunately, most of statistical properties can be approximated accurately when the number reaches a certain threshold.The SOC simulator for the kth Rayleigh fading process is defined as [27] μ where  is the number of cisoids.The real-valued parameters, including gains ĉ, , phases θ, , and discrete Doppler frequencies, f, , remain constant during one simulation trial, which results in a deterministic channel in each snapshot.The gains ĉ, are given such that {ĉ 2 , } =  2 μ / for  = 1, 2, . . ., , where  2 μ is the average power of channel fading.The phases θ, are independent random variables uniformly distributed over [−, ).The discrete Doppler frequencies f, can be characterized by either      () or      ().How to generate these two functions properly is the key issue.
As a deterministic simulation model, the statistical properties of output fading processes should be calculated by using time average instead of statistical average.According to the definitions of ACF, CCF, and DPSD, we can calculate the time-averaged functions as follows: To satisfy with the requirements of ( 7), the following conditions for discrete Doppler frequency parameter should be fulfilled [29]:  2, which also shows the signal of each branch is a frequency modulated signal.
The envelope and phase PDFs are independent of Doppler frequencies according to (12)- (18).Thus, the simulated envelope and phase PDFs based on the proposed model are the same as those for the conventional SOC model.Because the conventional definition of DPSD is invalid for nonstationary processes, we redefine the DPSD as the short-time DPSD which can be calculated by the squared amplitude of signal's short-time Fourier transform (STFT) where (−) is an analysis time window, which is sufficiently short such that the process can be considered to be stationary.Now we can redefine the time-variant ACF as the inverse Fourier transform of  μ μ (, ) Similarly, the CCF of a nonstationary process can be defined as and only consider the cases with  = 0.

Parameters Computation.
The time-variant discrete frequency parameters will increase the complexity and uncertainty of our model, so it is very important to find an efficient  parameter computation method to update the value of frequency for each branch over time.
Assuming that the  Rayleigh fading channels have the same time-variant DPSD, we sample the continuous DPSD with a short-time interval   (namely, channel updating interval).The updating interval is short enough, so the DPSD is approximately constant provided    () is within the uth interval.In addition, we propose the discrete frequency parameters following the linear change, which is where   , means the initial value of the th Rayleigh fading channel and the th branch of signal and    indicates the slope of the th branch for  fading channels.
The slope    is calculated as where    is the computation result of MEA [28], which makes the area under the DPSD curve within the range For the special case of a symmetric DPSD, we have   0 = −  max and    =   max with   max being the maximum Doppler frequency.  refers to the approximate number of periods for the Doppler frequency parameters linearly varying between the upper and lower boundaries.In this paper, we propose   ≈ 10 to ensure the mutual independence among  fading channels.When the value of f , () exceeds the upper boundary or reduces below the lower boundary, the updating direction will be reversed.This step keeps the value of Doppler frequency confined within the upper and lower boundaries.We set the upper boundary    () as a straight line connecting    and  +1  ; that is, where   is the start time of th interval.Similarly, the lower boundary expression can be obtained by setting  =  − 1.
Finally, we set the initial value   , of the th interval equal to the end frequency in the last time interval For the case of  = 1,  1  , is specified as a random variable uniformly distributed over the range of [ 1  −1 ,  1  ).By applying the aforementioned computation method of (25)-( 29), we can observe that the discrete frequency parameters generated are continuous and meet the requirements of (20).Additionally, our computation method will make the signal of each branch in (21) as a chirp signal, which are widely adopted in radar systems [32].

Simulation and Validation
In this section, we investigate the performance of the proposed simulation model by generating a 2 × 2 MIMO channel with center frequency  0 = 5 GHz.Assuming that there are many scatters around the receiver but few scatters around the transmitter, the typical averaged correlation matrixes are [33] The time-variant DPSD showed in Figure 3 is assumed to be symmetric Gaussian DPSD with   =  max ln 2.
According to the proposed simulation method, we should generate  = 4 i.i.d.Rayleigh fading processes firstly.Then, ).These results are calculated from ( 25)-( 29) and (31) with the DPSD given in Figure 3 and they will be used in (21) to generate the nonstationary Rayleigh fading.
with the given Gaussian DPSD, we get the area under the th DPSD curve over [−∞,    ] as So the values of    can be obtained by applying a proper numerical root-finding technique on (31).Using  = 16,   = 10, and   = 10ms as an example, we can get the time-variant discrete frequency parameters from (25) to (29). Figure 4 illustrates the computation results of frequency parameters for 16 branches in solid lines, and the upper and lower boundaries of each branch are showed in dash lines.The four independent Rayleigh fading processes generated by using the proposed simulation model exhibit similar statistical properties for the envelopes PDF, ACF, and DPSD. Figure 5 illustrates the simulated envelope PDF of the first Rayleigh fading process.Here, we assume the Rayleigh fading channel with normalized power, and as a result, c, = 1/16 is used, and the time-averaged PDF for all time intervals is given by dotted line depicted in Figure 5.The figure demonstrates that the simulation result provides a good approximation to the theoretical Rayleigh PDF.We also use (22) and (23) to calculate the time-variant short-time DPSD and ACF, which are showed, respectively, in Figures 6 and 7.For comparison purposes, we also highlight the edge lines of desired DPSD from Figure 3 and plot the calculated theoretical ACF with dotted line.It is apparent that the simulated channel's DPSD and ACF change over time and fit well with the desired graphs.Figures 5-7 also showed that the statistical results Inserting the matrix G and (33) into (3), we can finally get the nonstationary mutual correlated MIMO channel.Figure 8 gives the varying fading envelopes of all subchannels within 200 ms.The simulation accuracy for the MIMO channel is determined by the statistical properties of its subchannels.
According to the theory of random process, the linear combination of four i.i.d.complex Gaussian processes will not change the characteristics of the envelopes PDF, ACF, and which also exhibits good approximation to the desired one as (32), and the mean of simulation errors is merely 0.46%.The maximum of error is about 1.13%, while the minimum of error is less than 0.01%.

Conclusion
In this paper, a new simulation model originating from SOC method has been proposed to generate multiple nonstationary Rayleigh fading processes, which are very useful for the simulation of MIMO channel in time-variant richscattering environments as in V-V communication scenarios.Meanwhile, we assume that the frequency parameters of our model change linearly over time, which makes the signal of each branch (or subchannel) in the form of chirp signals.Under this assumption, an efficient updating method for the frequency parameters is proposed, which ensures that the statistical properties of simulated channel fading coincide with the desired ones.Simulation results demonstrate that our proposed simulator is able to generate nonstationary MIMO channels with time-variant statistical characteristics, such as the envelopes PDF, ACF, DPSD, and CCF approximating with their theoretical counterparts accurately.

Figure 1 :
Figure 1: MIMO communication system with  transmitting antennas and  receiving antennas under mobile to mobile (M2M) scenario.The solid line and dash line indicate the positions of transceiver antenna arrays at two different time instances.

Figure 2 :
Figure 2: Simulation model for nonstationary Rayleigh fading by a finite sum of frequency modulated signals.|μ  ()| and arg {μ  ()} are the amplitude and phase of the simulated fading process, respectively.

Figure 3 :Figure 4 :
Figure 3: Time-variant Gaussian-shaped DPSD used to verify the performance of proposed simulator.This DPSD corresponds to the scenario of mobile receiver starting acceleration.

Figure 5 :Figure 6 :
Figure 5: Time-averaged envelope PDF of generated fading channel.The simulation time and channel updating interval are 200 ms and 10 ms, respectively.The Rayleigh fading channels during all intervals are assumed to have normalized power.The dotted line denotes the averaged envelope PDF for all time intervals.

Figure 7 :
Figure 7: Magnitude of time-variant ACF of generated fading channel.The dotted line denotes the theoretical ACF corresponding to the Gaussian-shaped DPSD in Figure 3.

Figure 8 :
Figure 8: Envelopes of four subchannels generated by using the proposed simulation model.

Figure 9 :
Figure 9: Magnitude of time-variant CCF between individual pairs out of four subchannels.DPSD.So we only need to evaluate the performance of crosscorrelation between individual pairs of subchannels.Figure9depicts the absolute value of the CCF over time.Furthermore, the average of varying CCF can be calculated to get the correlation matrix of the simulated MIMO channel as Figure 9   depicts the absolute value of the CCF over time.Furthermore, the average of varying CCF can be calculated to get the correlation matrix of the simulated MIMO channel as