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

Authors: Qiuming Zhu; Xinglin Liu; Xuefeng Yin; Xiaomin Chen; Cuiwei Xue
Journal: International Journal of Antennas and Propagation
 (2016)

---

## Abstract

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.


---

## Body

## 1. 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–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–10]. Xiao et al. [8] were concerned with the time-variant parameters such as angles-of-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-of-motion (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.11 TGn 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–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–25] developed deterministic SOS simulators with low complexity. At the same time, the sum-of-cisoids (SOC) simulators are found to be more flexible and suitable for channel simulation in the nonisotropic scattering environments [26–29]. The accuracy of statistical channel properties obtained with both SOS and SOC simulators greatly depends on the simulation parameters computation method. The L
p-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 Section2, 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.


## 2. Nonstationary MIMO Channel Model

A MIMO mobile to mobile (M2M) communication system is illustrated in Figure1, 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.Figure 1
MIMO communication system withS transmitting antennas and U 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.Assuming a frequency nonselective MIMO channel system withS transmitting antennas and U receiving antennas, the MIMO channel matrix can be expressed as(1)
H
t
=
h
u
,
s
t
U
×
S
,where h
u
,
s
(
t
) denotes the impulse response of the subchannel between theuth (u
=
1,2
,
…
,
U) receiving antenna and thesth (s
=
1,2
,
…
,
S) transmitting antenna. Considering a Rayleigh fading channel, the envelope and phase distribution of h
u
,
s
(
t
) 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(2)
R
H
t
=
ρ
v
e
c
H
v
e
c
H
H
,where v
e
c
(
·
) is the vector operator (stacking all elements of a matrix column-wise into a single vector), (
·
)
H denotes the Hermitian transposition, and ρ
{
·
} represents the correlation coefficient. R
H
(
t
) is a U
S
×
U
S matrix, and we can use it to remodel the nonstationary correlated MIMO channel as (3)
v
e
c
H
=
R
H
1
/
2
v
e
c
G
,where (
·
)
1
/
2 denotes any matrix square root fulfilling R
H
1
/
2
(
R
H
1
/
2
)
H
=
R
H and G is an U
×
S 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 the transmitter correlation matrix R
T
x
(
t
) and the receiver correlation matrix R
R
x
(
t
), (4)
R
H
t
=
R
T
x
t
⊗
R
R
x
t
,where ⊗ denotes the Kronecker product.


## 3. 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 matrixL by performing Cholesky decomposition on R
H, L
=
R
H
1
/
2, and R
H
=
L
L
H. 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.


### 3.1. The Reference Model

Multiple Rayleigh fading processes are actually equivalent to a set of complex Gaussian random processes. Under this assumption, thekth Rayleigh fading process can be written as(5)
μ
k
t
=
μ
k
,
i
t
+
j
μ
k
,
q
t
,
k
=
1,2
,
…
,
K
,where j
≜
-
1, μ
k
,
i
(
t
), and μ
k
,
q
(
t
) are uncorrelated and zero mean real-valued Gaussian processes. Thus, the fading envelope ζ
k
(
t
)
≜
μ
k
t and phase ϕ
k
(
t
)
≜
a
r
g
{
μ
k
(
t
)
} yield the following PDFs, respectively: (6)
p
ζ
k
z
=
2
z
σ
μ
k
2
·
e
x
p
-
z
2
σ
μ
k
2
,
z
≥
0
,
p
ϕ
k
θ
=
1
2
π
,
θ
∈
-
π
,
π
,where σ
μ
k
2 denotes the variance of Rayleigh distribution. The autocorrelated Rayleigh fading process μ
k
(
t
) can be further characterized by its ACF r
μ
k
μ
k
(
τ
)
≜
E
{
μ
k
∗
(
t
)
μ
k
(
t
+
τ
)
}, while the cross-correlated Rayleigh fading processes are characterized by their CCF r
μ
k
μ
l
(
τ
)
≜
E
{
μ
k
∗
(
t
)
μ
l
(
t
+
τ
)
}
,
∀
k
≠
l. Here, (
·
)
∗ is complex conjugation and E
{
·
} refers to the expectation of given argument. The DPSD S
μ
k
μ
k
(
f
) describes the power distribution over Doppler frequency and can be calculated by performing the Fourier transform on the ACF; that is, S
μ
k
μ
k
(
f
)
=
∫
-
∞
∞
r
μ
k
μ
k
(
τ
)
e
x
p
{
-
j
2
π
f
τ
}
d
τ. Our aim is to generate K i.i.d. Rayleigh fading processes, which satisfy the following conditions: (7)
r
μ
k
μ
k
τ
=
σ
μ
k
2
∫
-
∞
∞
p
k
,
f
f
e
x
p
j
2
π
f
τ
d
f
,
S
μ
k
μ
k
f
=
σ
μ
k
2
p
k
,
f
f
,
r
μ
k
μ
l
τ
=
0
,
∀
k
≠
l
,where p
k
,
f
(
f
) is the Doppler frequency PDF of thekth Rayleigh fading process. All elements of the antenna arrays share a common scattering environment where the small-scale characterization assumption applies. Thus, the K random processes have the same Doppler frequency PDF denoted as p
f
(
f
).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:(8)
S
μ
μ
G
f
=
σ
μ
2
f
c
l
n
2
π
e
-
l
n
2
f
/
f
c
2
,
f
∈
f
l
,
f
u
,
(9)
r
μ
μ
G
τ
=
∫
f
l
f
u
S
μ
μ
G
f
e
x
p
j
2
π
f
τ
d
f
,where f
c denotes the 3-dB-cut-off frequency and f
l
,
f
u refer to the minimum and maximum frequency, respectively. For a symmetric Gaussian-shaped DPSD, the expression of(9) can be simplified as (10)
r
μ
μ
G
τ
=
σ
μ
2
e
-
π
f
c
/
l
n
2
τ
2
,where f
l
=
-
f
m
a
x and f
u
=
f
m
a
x with f
m
a
x
=
f
0
·
v
/
c denoting the maximum Doppler frequency due to the movements of transceivers. Here c
=
3
×
1
0
8 refers to the velocity of light, and f
0 is the carrier frequency.


### 3.2. 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 thekth Rayleigh fading process is defined as [27](11)
μ
^
k
t
=
∑
n
=
1
N
c
^
k
,
n
e
x
p
j
2
π
f
^
k
,
n
t
+
θ
^
k
,
n
,where N is the number of cisoids. The real-valued parameters, including gains c
^
k
,
n, phases θ
^
k
,
n, and discrete Doppler frequencies, f
^
k
,
n, remain constant during one simulation trial, which results in a deterministic channel in each snapshot. The gains c
^
k
,
n are given such that E
{
c
^
k
,
n
2
}
=
σ
μ
^
k
2
/
N for n
=
1,2
,
…
,
N, where σ
μ
^
k
2 is the average power of channel fading. The phases θ
^
k
,
n are independent random variables uniformly distributed over [
-
π
,
π
). The discrete Doppler frequencies f
^
k
,
n can be characterized by either S
μ
k
μ
k
(
f
) or r
μ
k
μ
k
(
f
). How to generate these two functions properly is the key issue.


### 3.3. Statistical Characteristics

From (11), we can derive the joint distribution p
μ
^
k
,
i
,
n
μ
^
k
,
q
,
n
(
z
1
,
z
2
) of two random variables μ
^
k
,
i
,
n
(
t
) and μ
^
k
,
q
,
n
(
t
) as [31] (12)
p
μ
^
k
,
i
,
n
μ
^
k
,
q
,
n
z
1
,
z
2
=
p
μ
^
k
,
i
,
n
z
1
·
δ
z
2
-
g
z
1
,where (13)
p
μ
k
,
i
/
q
,
n
z
1
/
2
=
1
π
c
k
,
n
2
-
z
1
/
2
2
,
g
z
1
=
c
k
,
n
2
-
z
1
2
,when z
1
/
2
<
c
k
,
n. The join characteristic function can be defined by the Fourier transform of joint probability density function as(14)
ψ
x
,
y
v
1
,
v
2
=
∬
p
x
,
y
z
1
,
z
2
e
j
2
π
v
1
z
1
+
v
2
z
2
d
z
1
d
z
2
.Combining with following equation(15)
ψ
x
,
y
v
1
,
v
2
=
∏
ψ
x
n
,
y
n
v
1
,
v
2
,where x
(
t
)
=
∑
n
=
1
N
x
n
(
t
), y
(
t
)
=
∑
n
=
1
N
y
n
(
t
), and x
n
(
t
) and y
n
(
t
) are i.i.d. random variables, we can get (16)
p
μ
^
k
,
i
μ
^
k
,
q
z
1
,
z
2
=
2
π
∫
0
∞
∏
n
=
1
N
J
0
2
π
c
k
,
n
v
J
0
2
π
v
z
1
2
+
z
2
2
v
d
v
.Finally, transforming the Cartesian coordination into the polar coordination with z
1
=
z
c
o
s
θ and z
2
=
z
s
i
n
θ, the PDFs of envelope and phase can be expressed, respectively, as(17)
p
ζ
^
k
z
=
z
2
π
2
∫
0
∞
∏
n
=
1
N
J
0
2
π
c
^
k
,
n
v
J
0
2
π
z
v
v
d
v
,
z
≥
0
,
(18)
p
ϕ
^
k
θ
=
1
2
π
,
θ
∈
-
π
,
πwhen c
^
k
,
n
=
σ
μ
^
k
/
N and N
→
∞ are considered; (17) is the same as the Rayleigh PDF.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:(19)
r
μ
^
k
μ
^
k
τ
=
∑
n
=
1
N
c
^
k
,
n
2
e
x
p
j
2
π
f
^
k
,
n
τ
,
S
μ
^
k
μ
^
k
f
=
∑
n
=
1
N
c
^
k
,
n
2
δ
f
-
f
^
k
,
n
,
r
μ
^
k
μ
^
l
τ
=
∑
n
=
1
N
∑
m
=
1
N
δ
f
^
k
,
n
-
f
^
l
,
m
e
x
p
-
j
2
π
f
^
l
,
m
τ
e
x
p
-
j
θ
^
k
,
n
-
θ
^
l
,
m
.To satisfy with the requirements of (7), the following conditions for discrete Doppler frequency parameter should be fulfilled [29]:(20)
f
k
,
n
≠
0
,
f
k
,
n
≠
f
k
,
m
,
f
k
,
n
≠
f
l
,
m
,
∀
k
,
n
,
l
,
m
,
k
≠
l
,
m
≠
n
.


## 3.1. The Reference Model

Multiple Rayleigh fading processes are actually equivalent to a set of complex Gaussian random processes. Under this assumption, thekth Rayleigh fading process can be written as(5)
μ
k
t
=
μ
k
,
i
t
+
j
μ
k
,
q
t
,
k
=
1,2
,
…
,
K
,where j
≜
-
1, μ
k
,
i
(
t
), and μ
k
,
q
(
t
) are uncorrelated and zero mean real-valued Gaussian processes. Thus, the fading envelope ζ
k
(
t
)
≜
μ
k
t and phase ϕ
k
(
t
)
≜
a
r
g
{
μ
k
(
t
)
} yield the following PDFs, respectively: (6)
p
ζ
k
z
=
2
z
σ
μ
k
2
·
e
x
p
-
z
2
σ
μ
k
2
,
z
≥
0
,
p
ϕ
k
θ
=
1
2
π
,
θ
∈
-
π
,
π
,where σ
μ
k
2 denotes the variance of Rayleigh distribution. The autocorrelated Rayleigh fading process μ
k
(
t
) can be further characterized by its ACF r
μ
k
μ
k
(
τ
)
≜
E
{
μ
k
∗
(
t
)
μ
k
(
t
+
τ
)
}, while the cross-correlated Rayleigh fading processes are characterized by their CCF r
μ
k
μ
l
(
τ
)
≜
E
{
μ
k
∗
(
t
)
μ
l
(
t
+
τ
)
}
,
∀
k
≠
l. Here, (
·
)
∗ is complex conjugation and E
{
·
} refers to the expectation of given argument. The DPSD S
μ
k
μ
k
(
f
) describes the power distribution over Doppler frequency and can be calculated by performing the Fourier transform on the ACF; that is, S
μ
k
μ
k
(
f
)
=
∫
-
∞
∞
r
μ
k
μ
k
(
τ
)
e
x
p
{
-
j
2
π
f
τ
}
d
τ. Our aim is to generate K i.i.d. Rayleigh fading processes, which satisfy the following conditions: (7)
r
μ
k
μ
k
τ
=
σ
μ
k
2
∫
-
∞
∞
p
k
,
f
f
e
x
p
j
2
π
f
τ
d
f
,
S
μ
k
μ
k
f
=
σ
μ
k
2
p
k
,
f
f
,
r
μ
k
μ
l
τ
=
0
,
∀
k
≠
l
,where p
k
,
f
(
f
) is the Doppler frequency PDF of thekth Rayleigh fading process. All elements of the antenna arrays share a common scattering environment where the small-scale characterization assumption applies. Thus, the K random processes have the same Doppler frequency PDF denoted as p
f
(
f
).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:(8)
S
μ
μ
G
f
=
σ
μ
2
f
c
l
n
2
π
e
-
l
n
2
f
/
f
c
2
,
f
∈
f
l
,
f
u
,
(9)
r
μ
μ
G
τ
=
∫
f
l
f
u
S
μ
μ
G
f
e
x
p
j
2
π
f
τ
d
f
,where f
c denotes the 3-dB-cut-off frequency and f
l
,
f
u refer to the minimum and maximum frequency, respectively. For a symmetric Gaussian-shaped DPSD, the expression of(9) can be simplified as (10)
r
μ
μ
G
τ
=
σ
μ
2
e
-
π
f
c
/
l
n
2
τ
2
,where f
l
=
-
f
m
a
x and f
u
=
f
m
a
x with f
m
a
x
=
f
0
·
v
/
c denoting the maximum Doppler frequency due to the movements of transceivers. Here c
=
3
×
1
0
8 refers to the velocity of light, and f
0 is the carrier frequency.


## 3.2. 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 thekth Rayleigh fading process is defined as [27](11)
μ
^
k
t
=
∑
n
=
1
N
c
^
k
,
n
e
x
p
j
2
π
f
^
k
,
n
t
+
θ
^
k
,
n
,where N is the number of cisoids. The real-valued parameters, including gains c
^
k
,
n, phases θ
^
k
,
n, and discrete Doppler frequencies, f
^
k
,
n, remain constant during one simulation trial, which results in a deterministic channel in each snapshot. The gains c
^
k
,
n are given such that E
{
c
^
k
,
n
2
}
=
σ
μ
^
k
2
/
N for n
=
1,2
,
…
,
N, where σ
μ
^
k
2 is the average power of channel fading. The phases θ
^
k
,
n are independent random variables uniformly distributed over [
-
π
,
π
). The discrete Doppler frequencies f
^
k
,
n can be characterized by either S
μ
k
μ
k
(
f
) or r
μ
k
μ
k
(
f
). How to generate these two functions properly is the key issue.


## 3.3. Statistical Characteristics

From (11), we can derive the joint distribution p
μ
^
k
,
i
,
n
μ
^
k
,
q
,
n
(
z
1
,
z
2
) of two random variables μ
^
k
,
i
,
n
(
t
) and μ
^
k
,
q
,
n
(
t
) as [31] (12)
p
μ
^
k
,
i
,
n
μ
^
k
,
q
,
n
z
1
,
z
2
=
p
μ
^
k
,
i
,
n
z
1
·
δ
z
2
-
g
z
1
,where (13)
p
μ
k
,
i
/
q
,
n
z
1
/
2
=
1
π
c
k
,
n
2
-
z
1
/
2
2
,
g
z
1
=
c
k
,
n
2
-
z
1
2
,when z
1
/
2
<
c
k
,
n. The join characteristic function can be defined by the Fourier transform of joint probability density function as(14)
ψ
x
,
y
v
1
,
v
2
=
∬
p
x
,
y
z
1
,
z
2
e
j
2
π
v
1
z
1
+
v
2
z
2
d
z
1
d
z
2
.Combining with following equation(15)
ψ
x
,
y
v
1
,
v
2
=
∏
ψ
x
n
,
y
n
v
1
,
v
2
,where x
(
t
)
=
∑
n
=
1
N
x
n
(
t
), y
(
t
)
=
∑
n
=
1
N
y
n
(
t
), and x
n
(
t
) and y
n
(
t
) are i.i.d. random variables, we can get (16)
p
μ
^
k
,
i
μ
^
k
,
q
z
1
,
z
2
=
2
π
∫
0
∞
∏
n
=
1
N
J
0
2
π
c
k
,
n
v
J
0
2
π
v
z
1
2
+
z
2
2
v
d
v
.Finally, transforming the Cartesian coordination into the polar coordination with z
1
=
z
c
o
s
θ and z
2
=
z
s
i
n
θ, the PDFs of envelope and phase can be expressed, respectively, as(17)
p
ζ
^
k
z
=
z
2
π
2
∫
0
∞
∏
n
=
1
N
J
0
2
π
c
^
k
,
n
v
J
0
2
π
z
v
v
d
v
,
z
≥
0
,
(18)
p
ϕ
^
k
θ
=
1
2
π
,
θ
∈
-
π
,
πwhen c
^
k
,
n
=
σ
μ
^
k
/
N and N
→
∞ are considered; (17) is the same as the Rayleigh PDF.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:(19)
r
μ
^
k
μ
^
k
τ
=
∑
n
=
1
N
c
^
k
,
n
2
e
x
p
j
2
π
f
^
k
,
n
τ
,
S
μ
^
k
μ
^
k
f
=
∑
n
=
1
N
c
^
k
,
n
2
δ
f
-
f
^
k
,
n
,
r
μ
^
k
μ
^
l
τ
=
∑
n
=
1
N
∑
m
=
1
N
δ
f
^
k
,
n
-
f
^
l
,
m
e
x
p
-
j
2
π
f
^
l
,
m
τ
e
x
p
-
j
θ
^
k
,
n
-
θ
^
l
,
m
.To satisfy with the requirements of (7), the following conditions for discrete Doppler frequency parameter should be fulfilled [29]:(20)
f
k
,
n
≠
0
,
f
k
,
n
≠
f
k
,
m
,
f
k
,
n
≠
f
l
,
m
,
∀
k
,
n
,
l
,
m
,
k
≠
l
,
m
≠
n
.


## 4. Nonstationary Channel Simulator

### 4.1. A New Simulation Model

The issue at hand is that the existing SOC simulator only applies to the stationary channel. To overcome this disadvantage, we propose a new simulator based on the sum of chirp signals as follows:(21)
μ
~
k
t
=
∑
n
=
1
N
c
~
k
,
n
e
x
p
j
2
π
∫
0
t
f
~
k
,
n
t
d
t
+
θ
~
k
,
n
,where the discrete frequency parameters are changed over time and characterized by the time-variant DPSD. It can be observed that our model will reduce to the original SOC model if f
~
k
,
n
(
t
) is time-invariant. The flow chart of the proposed simulation model is illustrated in Figure 2, which also shows the signal of each branch is a frequency modulated signal.Figure 2
Simulation model for nonstationary Rayleigh fading by a finite sum of frequency modulated signals.μ
~
k
t and a
r
g
{
μ
~
k
(
t
)
} are the amplitude and phase of the simulated fading process, respectively.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)(22)
T
μ
~
k
μ
~
k
t
,
f
≜
∫
-
∞
∞
μ
~
k
τ
e
-
j
2
π
f
τ
w
∗
t
-
τ
d
τ
,where w
(
t
-
τ
) 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 T
μ
~
k
μ
~
k
(
t
,
f
)
(23)
r
μ
~
k
μ
~
k
t
,
τ
=
∫
-
∞
∞
T
μ
~
k
μ
~
k
t
,
f
e
x
p
j
2
π
f
τ
d
f
.Similarly, the CCF of a nonstationary process can be defined as(24)
r
μ
~
k
μ
~
l
t
,
τ
≜
E
μ
~
k
∗
t
+
τ
2
μ
~
l
t
-
τ
2and only consider the cases with τ
=
0.


### 4.2. 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 theK Rayleigh fading channels have the same time-variant DPSD, we sample the continuous DPSD with a short-time interval T
u (namely, channel updating interval). The updating interval is short enough, so the DPSD is approximately constant provided S
μ
μ
u
(
f
) is within theuth interval. In addition, we propose the discrete frequency parameters following the linear change, which is(25)
f
~
k
,
n
u
t
=
b
k
,
n
u
+
s
n
u
t
,
t
∈
u
-
1
T
u
,
u
T
u
,
n
=
1,2
,
…
,
N
,
k
=
1,2
,
…
,
K
,
u
=
1,2
,
…
,where b
k
,
n
u means the initial value of the k
t
h Rayleigh fading channel and the n
t
h branch of signal and s
n
u indicates the slope of the n
t
h branch for K fading channels.The slopes
n
u is calculated as(26)
s
n
u
=
N
s
f
F
n
u
-
F
n
-
1
u
+
F
n
u
+
1
-
F
n
-
1
u
+
1
2
T
u
,where F
n
u is the computation result of MEA [28], which makes the area under the DPSD curve within the range [
F
n
-
1
u
,
F
n
u
) equal to σ
μ
~
2
/
N
(27)
∫
F
n
-
1
u
F
n
u
S
μ
~
μ
~
u
f
d
f
=
σ
μ
~
2
N
,
n
=
1,2
,
…
,
N
.For the special case of a symmetric DPSD, we haveF
0
u
=
-
f
m
a
x
u and F
N
u
=
f
m
a
x
u with f
m
a
x
u being the maximum Doppler frequency. N
s
f 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 N
s
f
≈
10 to ensure the mutual independence among K fading channels. When the value of f
~
k
,
n
u
(
t
) 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 B
n
u
(
t
) as a straight line connecting F
n
u and F
n
u
+
1; that is, (28)
B
n
u
t
=
F
n
u
+
1
-
F
n
u
T
u
t
-
t
u
+
F
n
u
,
n
=
1,2
,
…
,
N
,where t
u is the start time of u
t
h interval. Similarly, the lower boundary expression can be obtained by setting n
=
n
-
1.Finally, we set the initial valueb
k
,
n
u of the u
t
h interval equal to the end frequency in the last time interval (29)
b
k
,
n
u
=
U
F
n
-
1
u
,
F
n
u
u
=
1
f
~
k
,
n
u
-
1
u
=
2,3
,
…
,
∞
.For the case of u
=
1,  b
k
,
n
1 is specified as a random variable uniformly distributed over the range of [
F
n
-
1
1
,
F
n
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].


## 4.1. A New Simulation Model

The issue at hand is that the existing SOC simulator only applies to the stationary channel. To overcome this disadvantage, we propose a new simulator based on the sum of chirp signals as follows:(21)
μ
~
k
t
=
∑
n
=
1
N
c
~
k
,
n
e
x
p
j
2
π
∫
0
t
f
~
k
,
n
t
d
t
+
θ
~
k
,
n
,where the discrete frequency parameters are changed over time and characterized by the time-variant DPSD. It can be observed that our model will reduce to the original SOC model if f
~
k
,
n
(
t
) is time-invariant. The flow chart of the proposed simulation model is illustrated in Figure 2, which also shows the signal of each branch is a frequency modulated signal.Figure 2
Simulation model for nonstationary Rayleigh fading by a finite sum of frequency modulated signals.μ
~
k
t and a
r
g
{
μ
~
k
(
t
)
} are the amplitude and phase of the simulated fading process, respectively.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)(22)
T
μ
~
k
μ
~
k
t
,
f
≜
∫
-
∞
∞
μ
~
k
τ
e
-
j
2
π
f
τ
w
∗
t
-
τ
d
τ
,where w
(
t
-
τ
) 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 T
μ
~
k
μ
~
k
(
t
,
f
)
(23)
r
μ
~
k
μ
~
k
t
,
τ
=
∫
-
∞
∞
T
μ
~
k
μ
~
k
t
,
f
e
x
p
j
2
π
f
τ
d
f
.Similarly, the CCF of a nonstationary process can be defined as(24)
r
μ
~
k
μ
~
l
t
,
τ
≜
E
μ
~
k
∗
t
+
τ
2
μ
~
l
t
-
τ
2and only consider the cases with τ
=
0.


## 4.2. 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 theK Rayleigh fading channels have the same time-variant DPSD, we sample the continuous DPSD with a short-time interval T
u (namely, channel updating interval). The updating interval is short enough, so the DPSD is approximately constant provided S
μ
μ
u
(
f
) is within theuth interval. In addition, we propose the discrete frequency parameters following the linear change, which is(25)
f
~
k
,
n
u
t
=
b
k
,
n
u
+
s
n
u
t
,
t
∈
u
-
1
T
u
,
u
T
u
,
n
=
1,2
,
…
,
N
,
k
=
1,2
,
…
,
K
,
u
=
1,2
,
…
,where b
k
,
n
u means the initial value of the k
t
h Rayleigh fading channel and the n
t
h branch of signal and s
n
u indicates the slope of the n
t
h branch for K fading channels.The slopes
n
u is calculated as(26)
s
n
u
=
N
s
f
F
n
u
-
F
n
-
1
u
+
F
n
u
+
1
-
F
n
-
1
u
+
1
2
T
u
,where F
n
u is the computation result of MEA [28], which makes the area under the DPSD curve within the range [
F
n
-
1
u
,
F
n
u
) equal to σ
μ
~
2
/
N
(27)
∫
F
n
-
1
u
F
n
u
S
μ
~
μ
~
u
f
d
f
=
σ
μ
~
2
N
,
n
=
1,2
,
…
,
N
.For the special case of a symmetric DPSD, we haveF
0
u
=
-
f
m
a
x
u and F
N
u
=
f
m
a
x
u with f
m
a
x
u being the maximum Doppler frequency. N
s
f 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 N
s
f
≈
10 to ensure the mutual independence among K fading channels. When the value of f
~
k
,
n
u
(
t
) 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 B
n
u
(
t
) as a straight line connecting F
n
u and F
n
u
+
1; that is, (28)
B
n
u
t
=
F
n
u
+
1
-
F
n
u
T
u
t
-
t
u
+
F
n
u
,
n
=
1,2
,
…
,
N
,where t
u is the start time of u
t
h interval. Similarly, the lower boundary expression can be obtained by setting n
=
n
-
1.Finally, we set the initial valueb
k
,
n
u of the u
t
h interval equal to the end frequency in the last time interval (29)
b
k
,
n
u
=
U
F
n
-
1
u
,
F
n
u
u
=
1
f
~
k
,
n
u
-
1
u
=
2,3
,
…
,
∞
.For the case of u
=
1,  b
k
,
n
1 is specified as a random variable uniformly distributed over the range of [
F
n
-
1
1
,
F
n
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].


## 5. Simulation and Validation

In this section, we investigate the performance of the proposed simulation model by generating a2
×
2 MIMO channel with center frequency f
0
=
5
G
H
z. Assuming that there are many scatters around the receiver but few scatters around the transmitter, the typical averaged correlation matrixes are [33](30)
R
R
x
t
=
1
0.3
0.3
1
,
R
T
x
t
=
1
0.91
0.91
1
.The time-variant DPSD showed in Figure 3 is assumed to be symmetric Gaussian DPSD with f
c
=
f
m
a
x
l
n
2.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.According to the proposed simulation method, we should generateK
=
4 i.i.d. Rayleigh fading processes firstly. Then, with the given Gaussian DPSD, we get the area under the u
t
h DPSD curve over [
-
∞
,
F
n
u
] as(31)
G
F
n
u
=
∫
-
∞
F
n
u
S
μ
~
μ
~
u
f
d
f
=
σ
μ
~
2
2
1
+
e
r
f
F
n
u
f
c
l
n
2
,
n
=
1,2
,
…
,
N
.So the values of F
n
u can be obtained by applying a proper numerical root-finding technique on (31). Using N
=
16, N
s
f
=
10, and T
u
=
10
m
s 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.Figure 4
Calculated discrete frequency parameters for each branch (N
=
16). 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.The four independent Rayleigh fading processes generated by using the proposed simulation model exhibit similar statistical properties for the envelopes PDF, ACF, and DPSD. Figure5 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
~
k
,
n
=
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 obtained in each interval are very close to the theoretical value, which implies that our simulation model is applicable to reproducing the dynamic channels in the real world with tractable computational complexity.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 6
Time-variant DPSD of generated fading channel. The dotted line denotes the outline of desired Gaussian-shaped DPSD in Figure3.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 Figure3.Then, the four resultant i.i.d. Rayleigh fading processes are combined into the matrixG. The channel correlation matrix and the lower triangular matrix can be calculated by applying the Kronecker product and Cholesky decomposition as, respectively, (32)
R
H
=
1
0.3
0.91
0.273
0.3
1
0.273
0.91
0.91
0.273
1
0.3
0.273
0.91
0.3
1
,
(33)
L
=
1
0
0
0
0.3
0.954
0
0
0.91
0
0.415
0
0.273
0.868
0.124
0.396
.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 DPSD. So we only need to evaluate the performance of cross-correlation between individual pairs of subchannels. Figure 9 depicts the absolute value of the CCF over time.Figure 8
Envelopes of four subchannels generated by using the proposed simulation model.Figure 9
Magnitude of time-variant CCF between individual pairs out of four subchannels.Furthermore, the average of varying CCF can be calculated to get the correlation matrix of the simulated MIMO channel as(34)
R
~
H
=
1.0000
0.3003
0.9102
0.2753
0.3003
1.0000
0.2748
0.9099
0.9102
0.2748
1.0000
0.3034
0.2753
0.9099
0.3034
1.0000which 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%.


## 6. 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 rich-scattering 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.



---

*Source: 3492591-2016-10-26.xml*