# Joint Delay Doppler Probability Density Functions for Air-to-Air Channels

Authors: Michael Walter; Dmitriy Shutin; Uwe-Carsten Fiebig
Journal: International Journal of Antennas and Propagation
 (2014)

---

## Abstract

Recent channel measurements indicate that the wide sense stationary uncorrelated scattering assumption is not valid for air-to-air channels. Therefore, purely stochastic channel models cannot be used. In order to cope with the nonstationarity a geometric component is included. In this paper we extend a previously presented two-dimensional geometric stochastic model originally developed for vehicle-to-vehicle communication to a three-dimensional air-to-air channel model. Novel joint time-variant delay Doppler probability density functions are presented. The probability density functions are derived by using vector calculus and parametric equations of the delay ellipses. This allows us to obtain closed form mathematical expressions for the probability density functions, which can then be calculated for any delay and Doppler frequency at arbitrary times numerically.


---

## Body

## 1. Introduction

Channel modeling plays a vital role in the development of new communication systems, since channel models describe the physical radio propagation environment between the transmitter and receiver and thus enable the development of simulation tools to facilitate the receiver design.Statistical channel modeling research goes back to the mid-sixties, when Bello in his seminal work [1] characterized the statistical properties of randomly time-variant linear channels. In particular, he introduced the idea of wide sense stationary uncorrelated scattering (WSSUS) channels and provided its characterization in terms of the scattering function and its autocorrelation amongst others; the former describes the statistical relationship between the delay and Doppler frequency of propagation paths between the transmitter and receiver. The latter describes the correlation of the channel in time. Such channel models were widely applied to describe the propagation environment for terrestrial wireless communication systems, see [2–4].Aeronautical communications did not undergo a similar explosive growth as terrestrial mobile communication systems in the last decades. A much slower innovation cycle for communication systems in aeronautics is due to increased safety requirements and the corresponding high costs for certification and world-wide deployment of a system. Therefore, the investigation of aeronautical channels has not been that extensive. However, the next generation air transportation systems [5] will rely heavily on modern digital communication systems, such as ADS-B, [6], AeroMACS [7], LDACS, [8], and VDL mode 2 [9]. It is foreseen that air-to-air communications will become an integral part of the future air transportation system. For this reason an accurate knowledge of the air-to-air channel becomes a crucial enabler for future aeronautical communication systems.According to Bello [10], a surface scatter channel for air-to-satellite communications includes three key contributions: a line of sight (LOS) component, a specular reflection component, and a diffuse scattering component; the appearance of these three components is also typical for air-to-air channels [11]. The LOS naturally possesses the strongest power when it is received. The specular reflection depends mainly on the reflection properties of the surface, where the reflection takes place. Over water, for example, the received power of the specular reflection comes close to the LOS component. Both LOS and specular reflection are shifted in Doppler frequency due to the movement of transmitter and receiver. The third component, however, introduces a spread of the Doppler frequency, since the scattering occurs over the whole ground. In this paper we concentrate on the diffuse scattering component of the channel, since it influences the propagation the most in terms of delay and Doppler spread.The modeling results presented by Bello in [10] as well as other works [12, 13] for aeronautical channels are strongly based on the WSSUS assumption. Yet, recent measurement results have indicated a strong nonstationarity of air-to-air channels [11]. Therefore new channel models are needed that do not rely on the WSSUS assumption.In order to address the nonstationarity, Bello proposed a quasi-wide sense stationary uncorrelated scattering (QWSSUS) channel [1]. Essentially, for QWSSUS channels the WSSUS assumption can be assumed to hold over a limited time and bandwidth. Later, the local scattering function (LSF) was introduced by Matz in [14] for non-WSSUS channels. There, the channel does not have to be uncorrelated or stationary at all. The LSF can be seen as a generalization of the classical scattering function by Bello; the LSF merges into the scattering function in case of WSSUS channels. Statistics of nonstationary channels can be obtained by using the theory of time-frequency analysis to derive the LSF [15].Alternatively, the non-WSSUS case can be addressed using a geometric stochastic channel model (GSCM) [16, 17]. Within the GSCM approach, we depart from a purely stochastic modeling and allow for some assumptions on the geometry of the propagation environment. Specifically, in the GSCM approach the locations of transmitter, receiver, and scatterers are incorporated in the model. Yet, the geometric information is used to derive statistical distributions of the delay and the Doppler frequency.GSCM modeling has proven to be very effective in vehicle-to-vehicle (V2V) environments [18] that exhibit non-WSSUS behavior as well. Liberti and Rappaport proposed a geometrically based model for single bounce scatterers for microcell and PCS systems in [19]. Since their model was intended for personal communication devices with low antenna elevation, it is only two-dimensional, while aeronautical channel models should be investigated in three dimensions. Some results for three-dimensional models have already been derived for mobile or vehicle-to-vehicle communications [18, 20–25]. However, to the best of our knowledge no investigations of an aeronautical three-dimensional GSCM have been done yet. The placement of a moving transmitter and moving receiver has to be carried out in a three-dimensional coordinate system with the scatterers being placed on the ground.In this work, we aim to demonstrate the advantages of a new geometric stochastic modeling approach for air-to-air channels. Furthermore, we think that the statistical theory of non-WSSUS channel modeling and the GSCM can be associated. The LSF and our joint delay Doppler probability density function (pdf) will be put in relation to each other.These joint delay Doppler pdfs allow to simulate the bit error rate of a time-variant air-to-air channel at arbitrary times. Since the channel strongly depends on the flight trajectories, real or simulated trajectories can be fed into our model to accommodate the geometry of the setup. The scatterers are placed according to an adequate statistical distribution. Since the scattering properties of the ground are dependent on its dielectric constant, conductivity, and surface roughness the joint pdf needs to be scaled accordingly, if the received power for a givenτ and f
d is of interest. Typical values of the dielectric constant and the conductivity can be found in Table 1 of [26] and typical values of the root mean square of the surface height slope, which is related to the surface roughness, are provided in [10].The paper is structured as follows. Mathematical expressions for the delay and Doppler frequency are derived in the first part of Section2. Based on these expressions closed form mathematical expressions to compute the joint delay Doppler pdf are presented in the second part of Section 2. In Section 3, we explain the different scenarios including the simulation parameters. The graphical presentation of the results is shown in Section 4. There, we display the delay-dependent Doppler pdfs and the joint delay Doppler pdfs for the scenarios described in Section 3. Section 5 concludes the paper.


## 2. 3D GSCM

We proceed with the derivation of the GSCM in two steps. Initially, we outline the geometry of a typical air-to-air propagation scenario. Then, we will give a stochastic interpretation of the resulting propagation environment.


### 2.1. Geometric Part

Consider two aircraft flying behind each other. Without loss of generality, we define a Cartesian coordinate system so that both aircraft are contained in thex
-
z plane. The origin of the coordinate system is selected to coincide with the specular reflection point; see Figure 1. Such a choice will merely simplify the notation complexity of the model later in the text.Figure 1
Aircraft positions in the Cartesian coordinate system.Assuming a constant velocity, the time-variant positionsx
t
(
t
) and x
r
(
t
) of both aircraft are given as
(1)
x
t
(
t
)
=
x
t
0
+
v
t
t
,
x
r
(
t
)
=
x
r
0
+
v
r
t
,
where x
t
0
=
[
d
t
x
,
0
,
h
t
]
T and x
r
0
=
[
d
r
x
,
0
,
h
r
]
T describe the positions of the transmitter and of the receiver at the beginning of the observation at t
=
0; compare Figure 1. The velocities of the aircraft are given as v
t
=
[
v
t
x
,
v
t
y
,
v
t
z
]
T for the transmitter and v
r
=
[
v
r
x
,
v
r
y
,
v
r
z
]
T for the receiver. The projections of the slant range of the transmitter to the specular reflection point and the receiver to the specular reflection point onto the ground can be expressed as
(2)
d
t
x
=
-
h
t
h
t
+
h
r
d
,
d
r
x
=
h
r
h
t
+
h
r
d
,
where d is the horizontal distance between the aircraft and h
t and h
r the altitudes of the transmitting and receiving aircraft, respectively; compare Figure 1.In order to determine the distance between the aircraft via the scatterers on the ground, one can apply simple geometrical considerations. Let us consider a single scatterer at some time-invariant positionx in space. The total distance between transmitter and receiver via this scatterer d
sc
(
x
,
t
) can be calculated as
(3)
d
sc
(
x
,
t
)
=
d
t
(
x
,
t
)
+
d
r
(
x
,
t
)
,
where d
t
(
x
,
t
)
=
∥
x
-
x
t
(
t
)
∥ is the distance from transmitter to the scatterer and d
r
(
x
,
t
)
=
∥
x
-
x
r
(
t
)
∥ is the distance from the receiver to the scatterer, respectively. The propagation delay τ
(
x
,
t
) of a signal traveling from one aircraft to the other via scatterers is obtained simply by dividing the distance d
sc
(
x
,
t
) by the speed of light c resulting in τ
(
x
,
t
)
=
d
sc
(
x
,
t
)
/
c.Due to the movement of both the transmitter and receiver, the received signal also experiences a Doppler shiftf
d. The Doppler frequency for two moving aircraft via a scatterer can be computed as
(4)
f
d
(
β
,
φ
)
=
(
∥
v
t
∥
cos
⁡
(
β
(
t
)
)
+
∥
v
r
∥
cos
⁡
(
φ
(
t
)
)
)
f
c
c
,
where β
(
t
) is the angle between the velocity vector of the transmitter and the direct line from the transmitter to the scatterer and φ
(
t
) is the angle between the velocity vector of the receiver and direct line from the receiver to the scatterer, respectively. As can be seen from (4), the Doppler frequency f
d depends on a superposition of two cosine functions with angles β
(
t
) and φ
(
t
). If we wish to derive the pdf of the Doppler frequencies for a particular delay p
(
f
d
∣
τ
), then, in general, the statistical distributions of these angles need to be taken into account. Yet, the nonlinear transformation of these angles through the cosine function makes such a computation in a closed form impossible. In order to simplify such a calculation, we use a vector calculus description of the Doppler frequency f
d. In [10], it is shown that the Doppler frequency can be computed as temporal derivative of the delay. This temporal derivative, however, can also be computed using the spatial derivative [27]. Specifically, f
d can be computed as
(5)
f
d
(
x
,
t
)
=
-
∂
τ
(
x
,
t
)
∂
t
f
c
=
(
v
t
T
∇
d
t
(
x
,
t
)
+
v
r
T
∇
d
r
(
x
,
t
)
)
f
c
c
.
The Doppler frequency f
d being evaluated in this way allows us to derive the probability density function of the Doppler frequency in closed form. The result of the calculation is shown for t
=
0 in Appendix A.Therefore, we need to consider all possible locations of scatterers that cause an equal delay of the received signal, that is, all those scatterers for whichd
sc
(
x
,
t
)
=
const in (3). Obviously, such scatterers will all lie on the surface of an ellipsoid, with the transmit and receive aircraft located in the foci of this ellipsoid. Without loss of generality, we determine the parameters of the ellipsoid at time t
=
0. For t
≠
0, the positions of transmitter and receiver change and the ellipsoid needs to be rotated correspondingly, yet the basic computational steps remain the same. It is therefore more convenient to use a moving coordinate system, where the origin stays at the specular reflection point, if other time instances t
≠
0 are considered.Let us now define the parameters of this ellipsoid. One of its parameters is the focus distancee
f, also known as linear eccentricity, which is the distance from the center of the ellipsoid to the focal points. We can describe the focus distance e
f mathematically as
(6)
e
f
=
∥
x
r
(
t
)
-
x
t
(
t
)
∥
2
|
t
=
0
.
The shape of the ellipsoid is determined by the three semiprincipal axes. The largest semiprincipal axis a
oid is also called the semimajor axis. The constant distance of scatterers to the two focus points depends on the semimajor axis, so that we obtain d
sc
(
x
,
t
=
0
)
=
2
a
oid with a
oid
∈
[
e
f
,
∞
). The length of the other two semiprincipal axes b
oid and c
oid must be the same due to the isotropy of speed of light c, meaning b
oid
=
c
oid. The parameter b
oid only depends on the focus distance e
f and the semimajor axis a
oid. Thus, b
oid can be computed as
(7)
b
oid
=
a
oid
2
-
e
f
2
with
b
oid
≤
a
oid
.
The quadratic surface equation of an ellipsoid with its center at the origin, the semimajor axis along the x-axis, and the semiminor axes along the y and z axes can then be described by
(8)
x
T
(
(
1
a
oid
)
2
0
0
0
(
1
b
oid
)
2
0
0
0
(
1
b
oid
)
2
)
x
=
1
.
When an air-to-air scenario is considered with aircraft flying at different altitudes, the resulting ellipsoid will be rotated and shifted compared to (8). In the equation it is accounted for by introducing a linear shift and rotation into the ellipsoid equation. Due to the chosen coordinate system origin, the rotation occurs around the y-axis and the shift along the x-axis and the z-axis. After the rotation and the shift, the equations of the ellipsoid in the new coordinates x
′ can be described as
(9)
x
′
=
(
x
-
x
s
)
T
(
cos
⁡
α
0
sin
α
0
1
0
-
sin
α
0
cos
⁡
α
)
.
Here, the angle α is the rotation angle, which is measured clockwise, between the semimajor axis a
oid and the x-axis. x
s is the offset of the center of the ellipsoid from the origin of the coordinate system. The angle α can be computed as
(10)
α
=
arccos
(
(
x
r
(
t
)
-
x
t
(
t
)
)
∥
x
r
(
t
)
-
x
t
(
t
)
∥
T
e
x
)
|
t
=
0
,
where e
x is the unit vector in x-direction. The offset x
s is calculated using
(11)
x
s
=
x
t
(
t
)
+
x
r
(
t
)
2
|
t
=
0
.
Given the parameters of the ellipsoid, we can determine the location of all those scatterers arriving with an equal delay at the receiver. In what follows, we derive the probability distribution of Doppler frequencies for a specific delay.


### 2.2. Delay-Dependent Probability Density Function of the Doppler Frequency

Let us begin with the mathematical description of the shifted and rotated ellipsoid that represents all scatterers with equal delay. In general, an ellipsoid arbitrarily placed in a chosen coordinate system, can be described by the following quadratic surface equation:(12)
x
′
T
(
a
f
g
f
b
h
g
h
c
)
x
′
+
2
(
p
q
r
)
x
′
+
d
=
0
with a, b, c, d, f, and g being the coefficients of the quadratic surface equation; these parameters determine hereby the shape of the geometric figure. For the air-to-air channel the scatterers are located on the ground, that is, on the x
-
y plane. The curves of constant delay on the x
-
y plane are ellipses formed by the intersection of the ellipsoid with the ground plane. The expressions for such ellipses can be determined by setting all parameters related to z in (12) to zero. Since in our case the aircraft are placed in the x
-
z plane the linear terms containing y in (12) are set to zero and the equation simplifies to
(13)
x
′
T
(
a
0
0
0
b
0
0
0
0
)
x
′
+
2
(
p
0
0
)
x
′
+
d
=
0
.
As can be seen from (13), the principal axes of the ellipse are aligned with the x and y coordinate axes, so that we can use the results from [27] in our further description of the Doppler frequency.The semimajor axisa
ell, the semiminor axis b
ell, and the offset from the origin [
x
c
,
y
c
] of the ellipse can be calculated as
(14)
a
ell
=
p
2
a
b
-
d
b
,
b
ell
=
(
p
b
)
2
-
d
a
,
x
c
=
-
p
a
,
y
c
=
0
.
Given these elliptic parameters, it is now possible to calculate the Doppler frequencies f
d for a given delay τ, since the delay determines the parameters of the ellipse; in other words, we evaluate the Doppler frequency on the ellipse of constant delay. It follows that the Doppler frequency becomes a function of the ellipsoidal parameter ϕ. The parametric walk through the Doppler frequency results in f
d
(
x
(
t
=
0
;
ϕ
∣
τ
)
) being a function of
(15)
x
(
0
;
ϕ
∣
τ
)
=
(
a
ell
cos
⁡
ϕ
b
ell
sin
ϕ
)
+
(
x
c
y
c
)
,
where x
(
0
;
ϕ
∣
τ
) describes the location of a scatterer on the delay ellipse at the parameter ϕ. The Doppler frequency function f
d
(
x
(
0
;
ϕ
)
∣
τ
) can be evaluated as shown in Appendix B for the case of the aircraft being aligned in the x
-
z plane with arbitrary velocity vectors v
t and v
r. Since the Doppler frequency is only dependent on ϕ for a constant delay τ, the distribution of the Doppler frequencies p
(
f
d
∣
τ
) can be derived. Yet, a nonlinear transformation of the angular distribution p
(
ϕ
∣
τ
) is required beforehand.For the derivation of the distributionp
(
ϕ
∣
τ
), we assume that the scatterers lying on an equal arc length of the ellipse cause equal power at the receiver. This means that the scatterers are uniformly distributed along the ellipse. Since the ellipse is not circular, the angular distribution is not uniform anymore. Nørklit showed in [22] how such a distribution on an ellipse is calculated. The differential arc length is described by
(16)
d
l
=
a
ell
1
-
ϵ
τ
2
cos
⁡
2
ϕ
d
ϕ
,
where ϵ
τ
=
a
ell
2
-
b
ell
2
/
a
ell is the numerical eccentricity of the ellipse, which depends on the ratio of the semimajor and the semiminor axis for a given delay. Under the assumption of a uniform scatterer distribution on the ellipse, we can claim that all arcs of the same length l are “equal” in a sense that the density of scatterers on these arcs is the same. Thus, the pdf p
(
l
∣
τ
) must be uniform; that is, p
(
l
∣
τ
)
=
L
τ
-
1, where L
τ is the circumference of the ellipse. The circumference of an ellipse is given by
(17)
L
τ
=
4
a
ell
∫
0
π
/
2
‍
1
-
ϵ
τ
2
cos
⁡
2
ϕ
d
ϕ
,
where the integral expression is known as a complete elliptic integral of the second kind [28].The probability density function ofϕ is obtained by the following probability density transformation:
(18)
p
(
ϕ
)
=
p
(
l
)
|
d
l
d
ϕ
|
.
By inserting (16) and (17) into (18), we obtain
(19)
p
(
ϕ
∣
τ
)
=
1
-
ϵ
τ
2
cos
⁡
2
ϕ
E
(
ϵ
τ
)
,
with E
(
ϵ
τ
)
=
∫
0
π
/
2
‍
1
-
ϵ
τ
2
cos
⁡
2
ζ
d
ζ being the complete elliptical integral of the second kind with a cosine instead of a sine inside the integral. Note that (19) is dependent on the size of the ellipse. The numerical eccentricity, meaning the ratio between a
ell and b
ell, has an influence on the distribution p
(
ϕ
∣
τ
). For large delays, ϵ
τ
→
0 and we obtain a uniform distribution.If the scattering is nonisotropic but has a directional characteristic due to the scatterers' distribution on the ground or due to the antenna diagram, directional probability distributions might fit better to characterize the parameterϕ. In order to take the possible impact of the directionality into account, we propose to use the von Mises distribution [24], which is also known as circular normal distribution, together with the angular distribution given in (19). After some calculations this leads to
(20)
p
(
ϕ
∣
τ
;
κ
,
μ
)
=
1
-
ϵ
τ
2
cos
⁡
2
ϕ
E
(
ϵ
τ
)
×
exp
⁡
(
κ
cos
⁡
(
(
π
/
E
(
ϵ
τ
)
)
(
1
-
ϵ
τ
2
cos
⁡
2
ϕ
-
μ
)
)
)
J
0
(
κ
)
,
with κ being the concentration parameter of the distribution and μ being the centrality parameter of the distribution. That means the place on the ellipse, where the dominant scatterers come from. The centrality parameter is given by μ
=
∫
0
ϕ
¯
‍
1
-
ϵ
τ
2
cos
⁡
2
ζ
d
ζ and J
0
(
κ
) is the zero order Bessel function of the first kind. The polar angle θ and the ellipse parameter ϕ are related by tan
θ
=
1
-
ϵ
τ
2
tan
ϕ. The new probability density function in (21) generalizes the circular normal distribution to a normal distribution on an ellipse. In general the parameters κ and μ have to be estimated based on measurement data, or otherwise specified for a particular propagation environment.After having derived the distributionp
(
ϕ
∣
τ
), the probability density function of the Doppler frequency can be calculated through a nonlinear transformation of the pdf. However, this transformation is slightly more elaborate as the relation ℱ
:
ϕ
↦
f
d
(
ϕ
∣
τ
) is a many-to-one correspondence. It means that several values of ϕ can result in the same Doppler frequency f
d. The transformation of the delay-dependent pdf p
(
ϕ
∣
τ
) into the delay-dependent Doppler pdf p
(
f
d
∣
τ
) has to take this into account. The Doppler pdf p
(
f
d
∣
τ
) is computed as
(21)
p
(
f
d
∣
τ
)
=
∑
ϕ
′
∈
ℱ
-
1
(
f
d
)
p
(
ϕ
′
∣
τ
)
|
d
f
d
/
d
ϕ
′
|
,
where ℱ
-
1 is the inverse relation f
d
(
ϕ
∣
τ
)
↦
ϕ. Unfortunately, the exact derivation of ℱ
-
1 is quite involved. Nonetheless, it is possible to claim that the inverse relation ℱ
-
1 exists and that the set ℱ
-
1
(
f
d
) is countable and finite. This claim involves the enumeration of the intersection points of the delay ellipse with the curves of constant Doppler frequency. This is followed by invoking Bézout’s Theorem [29].Despiteℱ
-
1 not being available in closed form, we can compute it numerically. For each Doppler frequency f
d
* we solve f
d
*
-
f
d
(
ϕ
∣
τ
)
=
0 for ϕ in the interval [
0
,
2
π
). In general, the number of solutions to f
d
*
-
f
d
(
ϕ
∣
τ
)
=
0 can be different for different Doppler frequencies f
d
*. The calculation for the Jacobian element d
f
d
/
d
ϕ is provided in Appendix C.Let us point out that this derivation is similar to the derivation of the Jakes Doppler spectrum in [30]. We use, however, vector calculus to obtain the Doppler frequency f
d; see (5). When the eccentricity ϵ
τ of the ellipse becomes zero, that is, in the case of a circle, which is obtained as τ
→
∞, the Doppler spectrum becomes independent of the delay τ, and the Jakes spectrum emerges. In [27] it is proofed that the delay-dependent Doppler frequency pdf p
(
f
d
∣
τ
) indeed converges to the Jakes spectrum for τ
→
∞. In the following section, we show that the delay-dependent Doppler pdfs can deviate considerably from the Gaussian and Jakes-shaped Doppler spectra typically assumed in the literature [10, 26].Once the distribution for one specific delay and one specific time instance has been calculated, it can be repeated for other delays and time instances. When a moving coordinate system is used, the position and velocity vectors have to be updated for each time instantt, but the calculations presented above remain the same. This enables us to determine time-variant, delay-dependent Doppler spectra p
(
t
;
f
d
∣
τ
). In order to obtain the joint delay Doppler probability density function a time-variant delay distribution p
(
t
;
τ
) is needed. The joint delay Doppler pdf results in p
(
t
;
τ
,
f
d
)
=
p
(
t
;
τ
)
p
(
t
;
f
d
∣
τ
). In the following we assume a uniform distribution of the delay.


### 2.3. Relationship between Joint PDF and Received Power

The computational rules considered above allow to compute the joint distribution of delay and Doppler frequencyp
(
t
;
τ
,
f
d
). Yet the computed pdf does not necessarily provide information about the power distribution of the scattered signals without further assumptions. In the WSSUS case it can be shown that the scattering function of the channel is proportional to the joint delay Doppler pdf [3]. Similarly, for QWSSUS channels, this proportionality applies for a limited time period and frequency range.If the WSSUS assumption is violated, the situation is more complicated. Due to the correlations between amplitudes of scattered signals, as well as their nonstationary behavior in time, the areas of high probability inp
(
t
;
τ
,
f
d
) might not necessarily coincide with the areas of higher received power for the same parts in the delay Doppler plane. The generalized local scattering function (GLSF) [14] that measures the average strength of scatterers for a specific delay and Doppler, which is the analogon to the scattering function for non-WSSUS channels, does not explicitly consider the joint delay Doppler pdf; the relationship between them is yet to be studied.However, for the important class of non-WSSUS channels satisfying the so-called doubly underspread (DU) property [14], we can claim that the joint pdf p
(
t
;
τ
,
f
d
) is at least approximately proportional to the GLSF. Indeed, as it has been shown in [14], a DU channel can be locally approximated by the QWSSUS model. The latter model would imply the exact proportionality of the pdf p
(
t
;
τ
,
f
d
) and power received from the corresponding parts in the delay Doppler plane. Thus, the pdf p
(
t
;
τ
,
f
d
) must be locally approximately proportional to the GLSF.


## 2.1. Geometric Part

Consider two aircraft flying behind each other. Without loss of generality, we define a Cartesian coordinate system so that both aircraft are contained in thex
-
z plane. The origin of the coordinate system is selected to coincide with the specular reflection point; see Figure 1. Such a choice will merely simplify the notation complexity of the model later in the text.Figure 1
Aircraft positions in the Cartesian coordinate system.Assuming a constant velocity, the time-variant positionsx
t
(
t
) and x
r
(
t
) of both aircraft are given as
(1)
x
t
(
t
)
=
x
t
0
+
v
t
t
,
x
r
(
t
)
=
x
r
0
+
v
r
t
,
where x
t
0
=
[
d
t
x
,
0
,
h
t
]
T and x
r
0
=
[
d
r
x
,
0
,
h
r
]
T describe the positions of the transmitter and of the receiver at the beginning of the observation at t
=
0; compare Figure 1. The velocities of the aircraft are given as v
t
=
[
v
t
x
,
v
t
y
,
v
t
z
]
T for the transmitter and v
r
=
[
v
r
x
,
v
r
y
,
v
r
z
]
T for the receiver. The projections of the slant range of the transmitter to the specular reflection point and the receiver to the specular reflection point onto the ground can be expressed as
(2)
d
t
x
=
-
h
t
h
t
+
h
r
d
,
d
r
x
=
h
r
h
t
+
h
r
d
,
where d is the horizontal distance between the aircraft and h
t and h
r the altitudes of the transmitting and receiving aircraft, respectively; compare Figure 1.In order to determine the distance between the aircraft via the scatterers on the ground, one can apply simple geometrical considerations. Let us consider a single scatterer at some time-invariant positionx in space. The total distance between transmitter and receiver via this scatterer d
sc
(
x
,
t
) can be calculated as
(3)
d
sc
(
x
,
t
)
=
d
t
(
x
,
t
)
+
d
r
(
x
,
t
)
,
where d
t
(
x
,
t
)
=
∥
x
-
x
t
(
t
)
∥ is the distance from transmitter to the scatterer and d
r
(
x
,
t
)
=
∥
x
-
x
r
(
t
)
∥ is the distance from the receiver to the scatterer, respectively. The propagation delay τ
(
x
,
t
) of a signal traveling from one aircraft to the other via scatterers is obtained simply by dividing the distance d
sc
(
x
,
t
) by the speed of light c resulting in τ
(
x
,
t
)
=
d
sc
(
x
,
t
)
/
c.Due to the movement of both the transmitter and receiver, the received signal also experiences a Doppler shiftf
d. The Doppler frequency for two moving aircraft via a scatterer can be computed as
(4)
f
d
(
β
,
φ
)
=
(
∥
v
t
∥
cos
⁡
(
β
(
t
)
)
+
∥
v
r
∥
cos
⁡
(
φ
(
t
)
)
)
f
c
c
,
where β
(
t
) is the angle between the velocity vector of the transmitter and the direct line from the transmitter to the scatterer and φ
(
t
) is the angle between the velocity vector of the receiver and direct line from the receiver to the scatterer, respectively. As can be seen from (4), the Doppler frequency f
d depends on a superposition of two cosine functions with angles β
(
t
) and φ
(
t
). If we wish to derive the pdf of the Doppler frequencies for a particular delay p
(
f
d
∣
τ
), then, in general, the statistical distributions of these angles need to be taken into account. Yet, the nonlinear transformation of these angles through the cosine function makes such a computation in a closed form impossible. In order to simplify such a calculation, we use a vector calculus description of the Doppler frequency f
d. In [10], it is shown that the Doppler frequency can be computed as temporal derivative of the delay. This temporal derivative, however, can also be computed using the spatial derivative [27]. Specifically, f
d can be computed as
(5)
f
d
(
x
,
t
)
=
-
∂
τ
(
x
,
t
)
∂
t
f
c
=
(
v
t
T
∇
d
t
(
x
,
t
)
+
v
r
T
∇
d
r
(
x
,
t
)
)
f
c
c
.
The Doppler frequency f
d being evaluated in this way allows us to derive the probability density function of the Doppler frequency in closed form. The result of the calculation is shown for t
=
0 in Appendix A.Therefore, we need to consider all possible locations of scatterers that cause an equal delay of the received signal, that is, all those scatterers for whichd
sc
(
x
,
t
)
=
const in (3). Obviously, such scatterers will all lie on the surface of an ellipsoid, with the transmit and receive aircraft located in the foci of this ellipsoid. Without loss of generality, we determine the parameters of the ellipsoid at time t
=
0. For t
≠
0, the positions of transmitter and receiver change and the ellipsoid needs to be rotated correspondingly, yet the basic computational steps remain the same. It is therefore more convenient to use a moving coordinate system, where the origin stays at the specular reflection point, if other time instances t
≠
0 are considered.Let us now define the parameters of this ellipsoid. One of its parameters is the focus distancee
f, also known as linear eccentricity, which is the distance from the center of the ellipsoid to the focal points. We can describe the focus distance e
f mathematically as
(6)
e
f
=
∥
x
r
(
t
)
-
x
t
(
t
)
∥
2
|
t
=
0
.
The shape of the ellipsoid is determined by the three semiprincipal axes. The largest semiprincipal axis a
oid is also called the semimajor axis. The constant distance of scatterers to the two focus points depends on the semimajor axis, so that we obtain d
sc
(
x
,
t
=
0
)
=
2
a
oid with a
oid
∈
[
e
f
,
∞
). The length of the other two semiprincipal axes b
oid and c
oid must be the same due to the isotropy of speed of light c, meaning b
oid
=
c
oid. The parameter b
oid only depends on the focus distance e
f and the semimajor axis a
oid. Thus, b
oid can be computed as
(7)
b
oid
=
a
oid
2
-
e
f
2
with
b
oid
≤
a
oid
.
The quadratic surface equation of an ellipsoid with its center at the origin, the semimajor axis along the x-axis, and the semiminor axes along the y and z axes can then be described by
(8)
x
T
(
(
1
a
oid
)
2
0
0
0
(
1
b
oid
)
2
0
0
0
(
1
b
oid
)
2
)
x
=
1
.
When an air-to-air scenario is considered with aircraft flying at different altitudes, the resulting ellipsoid will be rotated and shifted compared to (8). In the equation it is accounted for by introducing a linear shift and rotation into the ellipsoid equation. Due to the chosen coordinate system origin, the rotation occurs around the y-axis and the shift along the x-axis and the z-axis. After the rotation and the shift, the equations of the ellipsoid in the new coordinates x
′ can be described as
(9)
x
′
=
(
x
-
x
s
)
T
(
cos
⁡
α
0
sin
α
0
1
0
-
sin
α
0
cos
⁡
α
)
.
Here, the angle α is the rotation angle, which is measured clockwise, between the semimajor axis a
oid and the x-axis. x
s is the offset of the center of the ellipsoid from the origin of the coordinate system. The angle α can be computed as
(10)
α
=
arccos
(
(
x
r
(
t
)
-
x
t
(
t
)
)
∥
x
r
(
t
)
-
x
t
(
t
)
∥
T
e
x
)
|
t
=
0
,
where e
x is the unit vector in x-direction. The offset x
s is calculated using
(11)
x
s
=
x
t
(
t
)
+
x
r
(
t
)
2
|
t
=
0
.
Given the parameters of the ellipsoid, we can determine the location of all those scatterers arriving with an equal delay at the receiver. In what follows, we derive the probability distribution of Doppler frequencies for a specific delay.


## 2.2. Delay-Dependent Probability Density Function of the Doppler Frequency

Let us begin with the mathematical description of the shifted and rotated ellipsoid that represents all scatterers with equal delay. In general, an ellipsoid arbitrarily placed in a chosen coordinate system, can be described by the following quadratic surface equation:(12)
x
′
T
(
a
f
g
f
b
h
g
h
c
)
x
′
+
2
(
p
q
r
)
x
′
+
d
=
0
with a, b, c, d, f, and g being the coefficients of the quadratic surface equation; these parameters determine hereby the shape of the geometric figure. For the air-to-air channel the scatterers are located on the ground, that is, on the x
-
y plane. The curves of constant delay on the x
-
y plane are ellipses formed by the intersection of the ellipsoid with the ground plane. The expressions for such ellipses can be determined by setting all parameters related to z in (12) to zero. Since in our case the aircraft are placed in the x
-
z plane the linear terms containing y in (12) are set to zero and the equation simplifies to
(13)
x
′
T
(
a
0
0
0
b
0
0
0
0
)
x
′
+
2
(
p
0
0
)
x
′
+
d
=
0
.
As can be seen from (13), the principal axes of the ellipse are aligned with the x and y coordinate axes, so that we can use the results from [27] in our further description of the Doppler frequency.The semimajor axisa
ell, the semiminor axis b
ell, and the offset from the origin [
x
c
,
y
c
] of the ellipse can be calculated as
(14)
a
ell
=
p
2
a
b
-
d
b
,
b
ell
=
(
p
b
)
2
-
d
a
,
x
c
=
-
p
a
,
y
c
=
0
.
Given these elliptic parameters, it is now possible to calculate the Doppler frequencies f
d for a given delay τ, since the delay determines the parameters of the ellipse; in other words, we evaluate the Doppler frequency on the ellipse of constant delay. It follows that the Doppler frequency becomes a function of the ellipsoidal parameter ϕ. The parametric walk through the Doppler frequency results in f
d
(
x
(
t
=
0
;
ϕ
∣
τ
)
) being a function of
(15)
x
(
0
;
ϕ
∣
τ
)
=
(
a
ell
cos
⁡
ϕ
b
ell
sin
ϕ
)
+
(
x
c
y
c
)
,
where x
(
0
;
ϕ
∣
τ
) describes the location of a scatterer on the delay ellipse at the parameter ϕ. The Doppler frequency function f
d
(
x
(
0
;
ϕ
)
∣
τ
) can be evaluated as shown in Appendix B for the case of the aircraft being aligned in the x
-
z plane with arbitrary velocity vectors v
t and v
r. Since the Doppler frequency is only dependent on ϕ for a constant delay τ, the distribution of the Doppler frequencies p
(
f
d
∣
τ
) can be derived. Yet, a nonlinear transformation of the angular distribution p
(
ϕ
∣
τ
) is required beforehand.For the derivation of the distributionp
(
ϕ
∣
τ
), we assume that the scatterers lying on an equal arc length of the ellipse cause equal power at the receiver. This means that the scatterers are uniformly distributed along the ellipse. Since the ellipse is not circular, the angular distribution is not uniform anymore. Nørklit showed in [22] how such a distribution on an ellipse is calculated. The differential arc length is described by
(16)
d
l
=
a
ell
1
-
ϵ
τ
2
cos
⁡
2
ϕ
d
ϕ
,
where ϵ
τ
=
a
ell
2
-
b
ell
2
/
a
ell is the numerical eccentricity of the ellipse, which depends on the ratio of the semimajor and the semiminor axis for a given delay. Under the assumption of a uniform scatterer distribution on the ellipse, we can claim that all arcs of the same length l are “equal” in a sense that the density of scatterers on these arcs is the same. Thus, the pdf p
(
l
∣
τ
) must be uniform; that is, p
(
l
∣
τ
)
=
L
τ
-
1, where L
τ is the circumference of the ellipse. The circumference of an ellipse is given by
(17)
L
τ
=
4
a
ell
∫
0
π
/
2
‍
1
-
ϵ
τ
2
cos
⁡
2
ϕ
d
ϕ
,
where the integral expression is known as a complete elliptic integral of the second kind [28].The probability density function ofϕ is obtained by the following probability density transformation:
(18)
p
(
ϕ
)
=
p
(
l
)
|
d
l
d
ϕ
|
.
By inserting (16) and (17) into (18), we obtain
(19)
p
(
ϕ
∣
τ
)
=
1
-
ϵ
τ
2
cos
⁡
2
ϕ
E
(
ϵ
τ
)
,
with E
(
ϵ
τ
)
=
∫
0
π
/
2
‍
1
-
ϵ
τ
2
cos
⁡
2
ζ
d
ζ being the complete elliptical integral of the second kind with a cosine instead of a sine inside the integral. Note that (19) is dependent on the size of the ellipse. The numerical eccentricity, meaning the ratio between a
ell and b
ell, has an influence on the distribution p
(
ϕ
∣
τ
). For large delays, ϵ
τ
→
0 and we obtain a uniform distribution.If the scattering is nonisotropic but has a directional characteristic due to the scatterers' distribution on the ground or due to the antenna diagram, directional probability distributions might fit better to characterize the parameterϕ. In order to take the possible impact of the directionality into account, we propose to use the von Mises distribution [24], which is also known as circular normal distribution, together with the angular distribution given in (19). After some calculations this leads to
(20)
p
(
ϕ
∣
τ
;
κ
,
μ
)
=
1
-
ϵ
τ
2
cos
⁡
2
ϕ
E
(
ϵ
τ
)
×
exp
⁡
(
κ
cos
⁡
(
(
π
/
E
(
ϵ
τ
)
)
(
1
-
ϵ
τ
2
cos
⁡
2
ϕ
-
μ
)
)
)
J
0
(
κ
)
,
with κ being the concentration parameter of the distribution and μ being the centrality parameter of the distribution. That means the place on the ellipse, where the dominant scatterers come from. The centrality parameter is given by μ
=
∫
0
ϕ
¯
‍
1
-
ϵ
τ
2
cos
⁡
2
ζ
d
ζ and J
0
(
κ
) is the zero order Bessel function of the first kind. The polar angle θ and the ellipse parameter ϕ are related by tan
θ
=
1
-
ϵ
τ
2
tan
ϕ. The new probability density function in (21) generalizes the circular normal distribution to a normal distribution on an ellipse. In general the parameters κ and μ have to be estimated based on measurement data, or otherwise specified for a particular propagation environment.After having derived the distributionp
(
ϕ
∣
τ
), the probability density function of the Doppler frequency can be calculated through a nonlinear transformation of the pdf. However, this transformation is slightly more elaborate as the relation ℱ
:
ϕ
↦
f
d
(
ϕ
∣
τ
) is a many-to-one correspondence. It means that several values of ϕ can result in the same Doppler frequency f
d. The transformation of the delay-dependent pdf p
(
ϕ
∣
τ
) into the delay-dependent Doppler pdf p
(
f
d
∣
τ
) has to take this into account. The Doppler pdf p
(
f
d
∣
τ
) is computed as
(21)
p
(
f
d
∣
τ
)
=
∑
ϕ
′
∈
ℱ
-
1
(
f
d
)
p
(
ϕ
′
∣
τ
)
|
d
f
d
/
d
ϕ
′
|
,
where ℱ
-
1 is the inverse relation f
d
(
ϕ
∣
τ
)
↦
ϕ. Unfortunately, the exact derivation of ℱ
-
1 is quite involved. Nonetheless, it is possible to claim that the inverse relation ℱ
-
1 exists and that the set ℱ
-
1
(
f
d
) is countable and finite. This claim involves the enumeration of the intersection points of the delay ellipse with the curves of constant Doppler frequency. This is followed by invoking Bézout’s Theorem [29].Despiteℱ
-
1 not being available in closed form, we can compute it numerically. For each Doppler frequency f
d
* we solve f
d
*
-
f
d
(
ϕ
∣
τ
)
=
0 for ϕ in the interval [
0
,
2
π
). In general, the number of solutions to f
d
*
-
f
d
(
ϕ
∣
τ
)
=
0 can be different for different Doppler frequencies f
d
*. The calculation for the Jacobian element d
f
d
/
d
ϕ is provided in Appendix C.Let us point out that this derivation is similar to the derivation of the Jakes Doppler spectrum in [30]. We use, however, vector calculus to obtain the Doppler frequency f
d; see (5). When the eccentricity ϵ
τ of the ellipse becomes zero, that is, in the case of a circle, which is obtained as τ
→
∞, the Doppler spectrum becomes independent of the delay τ, and the Jakes spectrum emerges. In [27] it is proofed that the delay-dependent Doppler frequency pdf p
(
f
d
∣
τ
) indeed converges to the Jakes spectrum for τ
→
∞. In the following section, we show that the delay-dependent Doppler pdfs can deviate considerably from the Gaussian and Jakes-shaped Doppler spectra typically assumed in the literature [10, 26].Once the distribution for one specific delay and one specific time instance has been calculated, it can be repeated for other delays and time instances. When a moving coordinate system is used, the position and velocity vectors have to be updated for each time instantt, but the calculations presented above remain the same. This enables us to determine time-variant, delay-dependent Doppler spectra p
(
t
;
f
d
∣
τ
). In order to obtain the joint delay Doppler probability density function a time-variant delay distribution p
(
t
;
τ
) is needed. The joint delay Doppler pdf results in p
(
t
;
τ
,
f
d
)
=
p
(
t
;
τ
)
p
(
t
;
f
d
∣
τ
). In the following we assume a uniform distribution of the delay.


## 2.3. Relationship between Joint PDF and Received Power

The computational rules considered above allow to compute the joint distribution of delay and Doppler frequencyp
(
t
;
τ
,
f
d
). Yet the computed pdf does not necessarily provide information about the power distribution of the scattered signals without further assumptions. In the WSSUS case it can be shown that the scattering function of the channel is proportional to the joint delay Doppler pdf [3]. Similarly, for QWSSUS channels, this proportionality applies for a limited time period and frequency range.If the WSSUS assumption is violated, the situation is more complicated. Due to the correlations between amplitudes of scattered signals, as well as their nonstationary behavior in time, the areas of high probability inp
(
t
;
τ
,
f
d
) might not necessarily coincide with the areas of higher received power for the same parts in the delay Doppler plane. The generalized local scattering function (GLSF) [14] that measures the average strength of scatterers for a specific delay and Doppler, which is the analogon to the scattering function for non-WSSUS channels, does not explicitly consider the joint delay Doppler pdf; the relationship between them is yet to be studied.However, for the important class of non-WSSUS channels satisfying the so-called doubly underspread (DU) property [14], we can claim that the joint pdf p
(
t
;
τ
,
f
d
) is at least approximately proportional to the GLSF. Indeed, as it has been shown in [14], a DU channel can be locally approximated by the QWSSUS model. The latter model would imply the exact proportionality of the pdf p
(
t
;
τ
,
f
d
) and power received from the corresponding parts in the delay Doppler plane. Thus, the pdf p
(
t
;
τ
,
f
d
) must be locally approximately proportional to the GLSF.


## 3. Scenarios

In this section, we briefly describe the different scenarios that are used to evaluate the proposed algorithm. The carrier frequency is chosen to bef
c
=
250 MHz in order to compare our theoretical results with the measurement data provided in [11]. This frequency is also close to the aeronautical VHF band, where currently the aeronautical digital communication (VDL mode 2 [9]) as well as analog voice communication takes place in. Let us point out that this frequency shows good scattering properties from the ground.The antenna pattern is assumed to be uniform and therefore not considered in the model. The transmitter and receiver aircraft are described by single points in the Cartesian coordinate system. The parameters of the selected scenarios are summarized below.(i)
In Scenario I the two aircraft fly behind each other at the same altitude and with the same speed in the same direction. The transmitter is at a positionx
t
=
[
-
1175,0
,
600
]
T m. Its velocity is v
t
=
[
70,0
,
0
]
T m/s. The receiver is at a position x
r
=
[
1175,0
,
600
]
T m and flies with a velocity v
r
=
[
70,0
,
0
]
T m/s. The scattering is assumed to be homogeneous.
(ii)
In Scenario II the two aircraft fly behind each other at different altitudes but with the same speed in the same direction. The transmitter is at a positionx
t
=
[
-
1175,0
,
1600
]
T m and has a velocity v
t
=
[
70,0
,
0
]
T m/s. The receiver is at x
r
=
[
1175,0
,
600
]
T m and has a velocity v
r
=
[
70,0
,
0
]
T m/s. The scattering is assumed to be homogeneous.
(iii)
In Scenario III the two aircraft fly towards each other at the same altitude and with the same speed in opposite directions. The transmitter is at a positionx
t
=
[
-
1175,0
,
600
]
T m and has a velocity v
t
=
[
70,0
,
0
]
T m/s. The receiver is at x
r
=
[
1175,0
,
600
]
T m and flies with a velocity v
r
=
[
-
70,0
,
0
]
T m/s. The scattering is assumed to be homogeneous.
(iv)
In Scenario IV the two aircraft fly behind each other at the same altitude and with the same speed in the same direction. The transmitter is atx
t
=
[
-
1175,0
,
600
]
T m and has a velocity v
t
=
[
70,0
,
0
]
T m/s. The receiver is at a position x
r
=
[
1175,0
,
600
]
T m and has a velocity v
r
=
[
70,0
,
0
]
T m/s. The scattering is assumed to be heterogeneous. In order to model the heterogeneity of scatterers we propose to make use of the von-Mises distribution as in (21) with the concentration parameter κ
=
0.5 and mean centrality angle θ
¯
=
π
/
4.Scenarios I and II are selected in such a way as to allow for a direct comparison with the measurement results in [11]. Scenarios III and IV provide predictions for different aircraft setups. They are selected to show different aspects of the channel that are not visible in Scenarios I and II.


## 4. Results

In this section, we evaluate the joint delay Doppler pdfsp
(
t
;
τ
,
f
d
) for the four different scenarios outlined above. The evaluation shows that the pdfs look significantly different compared to each other depending on the geometric setup.The steps to compute the delay-dependent Doppler pdfp
(
t
;
τ
,
f
d
) for an arbitrary time instant can be roughly represented with the pseudocode given in Algorithm 1.Algorithm 1:Computation of p
(
t
;
τ
,
f
d
).
(1)for  τ
[
n
]
∈
[
τ
min
⁡
,
τ
max
⁡
]  do
(2)  Compute the parameters a
ell, b
ell and ϵ
τ of the ellipse corresponding to the delay τ
[
n
]
(3)  Use (19) or (21) to compute p
(
ϕ
∣
τ
[
n
]
) for the corresponding scenario
(4)  for  f
d
[
m
]
∈
[
f
d
min
⁡
,
f
d
max
⁡
]  do
(5)  Compute  p
(
t
;
f
d
[
m
]
∣
τ
[
n
]
)  by numerically solving (21) for each  f
d
[
m
].
(6)  end for
(7)end forFigure2(a) shows the joint delay Doppler pdf for Scenario I. Let us point out that this scenario can be compared to a V2V channel analyzed in [27, 31], there however evaluated for the two-dimensional case. Here, the Doppler frequency is only slowly increasing with the delay. The difference compared to the V2V scenario is the placement of the delay ellipse. In the vehicular case, the ellipse always comprises all possible Doppler frequencies. In the aeronautical scenario the ellipse results from the intersection of the ellipsoid with the ground plane and therefore starts with a single frequency for the specular reflection point. As the delay ellipse grows it includes larger Doppler frequencies. The maximum and minimum Doppler frequency, which is given by f
d
=
(
∥
v
t
∥
+
∥
v
r
∥
)
f
c
/
c, is only reached asymptotically. This is due to the fact that only scatterers directly in line with the velocity vectors of the aircraft cause such a Doppler frequency.Scenario I for aircraft positions ofx
t
0
=
[
-
1175,0
,
600
]
T m and x
r
0
=
[
1175,0
,
600
]
T m and velocity vectors of v
t
=
v
r
=
[
70,0
,
0
]
T m/s. The aircraft fly in the same direction; homogeneous scattering occurs all over the ground.
(a)
Joint delay Doppler pdfp
(
t
;
τ
,
f
d
)
(b)
Doppler pdfp
(
t
;
f
d
∣
τ
) at t
=
0 s for the different delays indicated in (a)In order to get a clearer understanding of how the Doppler pdfs develop for increasing delay, we show four selected pdfs in Figure2(b). The figure shows that the Doppler spectra get wider for increasing delay. Furthermore, close to the specular reflection point at τ
=
9 μs, the spectrum exhibits a distinctive W shape, which is similar to the V2V case. It deviates from the Jakes spectrum both in terms of maximum and minimum frequency and shape. For larger delays, for example, τ
=
16 μs, however, this difference becomes smaller. Asymptotically, for τ
→
∞ the Jakes spectrum is obtained as a special case.Results for Scenario II are presented in Figures3(a) and 3(b). Observe that the difference in aircraft altitude causes an asymmetry of the joint delay Doppler pdf. For increasing delay the width of the spectrum for positive Doppler frequencies grows faster than for negative Doppler frequencies. Furthermore, the probability of receiving negative Doppler frequencies is slightly higher than the ones for positive Doppler frequencies. Due to the higher altitude of the transmitter as compared to the Scenario I, the joint delay Doppler pdf starts at a larger delay.Scenario II for aircraft positions ofx
t
0
=
[
-
1175,0
,
1600
]
T m  and x
r
0
=
[
1175,0
,
600
]
T m and velocity vectors of v
t
=
v
r
=
[
70,0
,
0
]
T m/s. The aircraft fly in the same direction; homogeneous scattering occurs all over the ground.
(a)
Joint delay Doppler pdfp
(
t
;
τ
,
f
d
)
(b)
Doppler pdfp
(
t
;
f
d
∣
τ
) at t
=
0 s for the different delays indicated in (a)Figure3(b) shows that the Doppler pdfs exhibit a different shape compared to Scenario I. Since the altitudes of transmitter and receiver are different, the spectrum does not show the W shape as in Figure 2(b) but rather looks like a skewed Jakes spectrum for τ
=
11 μs. For larger delays, for example, for τ
=
18 μs, the similarity to the Jakes spectrum again becomes more apparent and the pdf becomes more symmetric. Finally, the Jakes spectrum is asymptotically reached for τ
→
∞ as in Scenario I.A different shape of the joint delay Doppler pdf is obtained for Scenario III. Figure4(a) shows the joint delay Doppler pdf for these conditions. There the highest Doppler frequency occurs for the specular reflection point at τ
=
8.8 μs. For increasing delay, the joint pdf gets wider, reaches its maximum width between τ
=
12 μs and τ
=
14 μs and then gets smaller again. Furthermore, the maximum and minimum Doppler frequency decrease for increasing delay. The Doppler pdfs in Figure 4(b) exhibit a similar shape as the ones from Figure 3(b) but mirrored with respect to the vertical axis. However, they exhibit only positive Doppler frequencies. The width and the position of the pdfs change for the different delays. The largest spread of Doppler frequencies for the shown Doppler pdfs occurs at τ
=
12 μs. The closest one to the origin can be seen for τ
=
16 μs. Asymptotically, the pdf reaches Jakes again for τ
→
∞, but the width of spectrum reaches zero due to the contrary orientation of the velocity vectors.Scenario III for aircraft positions ofx
t
0
=
[
-
1175,0
,
600
]
T m and x
r
0
=
[
1175,0
,
600
]
T m and velocity vectors of v
t
=
-
v
r
=
[
70,0
,
0
]
T m/s. The aircraft fly in opposite directions; homogeneous scattering occurs all over the ground.
(a)
Joint delay Doppler pdfp
(
t
;
τ
,
f
d
)
(b)
Doppler pdfp
(
t
;
f
d
∣
τ
) at t
=
0 s for the different delays indicated in (a)In Scenario IV, we evaluate the joint delay Doppler pdf for a nonuniform scatterer distribution. Due to the use of von-Mises distribution, the probability is concentrated at certain Doppler frequencies. The von-Mises distribution can be used to model the effect of different terrains, for example, water and meadow land, which exhibit different scattering properties. In Figure5(a), the magnitude of the joint delay Doppler pdf is shifted towards the region with better scattering properties, which in our case is the region with positive Doppler frequencies. The envelope of the joint delay Doppler pdf, however, stays the same as in Scenario I. Figure 5(b) shows this effect in more detail. The W shape of the spectrum is warped towards positive Doppler frequencies, which can be seen especially well for small delays like at τ
=
9 μs. This means that the terrain in front of the two aircraft possesses the better scattering properties. The warping effect becomes less prominent for increasing delay. At τ
=
16 μs, the pdf is almost symmetric again like the one in Figure 2(b). The width of the pdfs is exactly as in Figure 2(b).Scenario IV for aircraft positions ofx
t
0
=
[
-
1175,0
,
600
]
T m and x
r
0
=
[
1175,0
,
600
]
T m and velocity vectors of v
t
=
v
r
=
[
70
,
0,0
]
T m/s. The aircraft fly in the same direction; heterogeneous scattering occurs all over the ground.
(a)
Joint delay Doppler pdfp
(
t
;
τ
,
f
d
)
(b)
Doppler pdfp
(
t
;
f
d
∣
τ
) at t
=
0 s for the different delays indicated in (a)


## 5. Conclusion

In this paper, we derived new joint delay Doppler probability density functions to exactly describe the influence of scattering on aeronautical air-to-air channels. This was achieved by calculating time-variant delay-dependent Doppler pdfs taking the exact flight trajectories into account. Geometric and algebraic methods were applied to calculate the Doppler frequencies for a specific delay. Finally, stochastic methods were used to determine the delay-dependent Doppler probability density functions. We presented a new probability density function that generalizes the von-Mises distribution for an ellipse. This directional distribution of the scatterers influences the amplitudes of the joint delay Doppler pdf but not the shape.Furthermore, the results show that the specific flight trajectories have a direct influence on the shape of the joint delay Doppler pdf. Our theoretic approach shows that it is not adequate to use the Jakes and Gaussian Doppler spectra for nonstationary air-to-air channels.The presented nonstationary theoretical model can be used to simulate time-variant air-to-air channels at arbitrary times. Real flight trajectories can be fed into the model to account for the geometry and the joint pdfs can be scaled by ground dependent coefficients to determine the received power for a given delayτ and Doppler frequency f
d.



---

*Source: 814218-2014-04-03.xml*