#  Spatial and Wavelength Division Joint Multiplexing System Design for
Visible Light Communications

Cheng Chen, Shenjie Huang, Iman Tavakkolnia, Majid Safari, and Harald Haas
Cheng Chen, Iman Tavakkolnia and Harald Haas are with the Department of
Electrical and Electrical Engineering, the University of Strathclyde, Glasgow,
Scotland, UK.Shenjie Huang and Majid Safari are with the institute for Digital
Communications, School of Engineering, University of Edinburgh, Edinburgh,
Scotland, UK.Harald Haas acknowledges the financial support from the Wolfson
Foundation, the Royal Society, and Engineering and Physical Sciences Research
Council (EPSRC) under Established Career Fellowship grant EP/R007101/1.

###### Abstract

The low-pass characteristics of front-end elements including light-emitting
diodes (LEDs) and photodiodes (PDs) limit the transmission data rate of
visible light communication (VLC) and Light Fidelity (LiFi) systems. Using
multiplexing transmission techniques, such as spatial multiplexing (SMX) and
wavelength division multiplexing (WDM), is a solution to overcome bandwidth
limitation.  However, spatial correlation in optical wireless channels and
optical filter bandpass shifts typically limit the achievable multiplexing
gain in SMX and WDM systems, respectively. In this paper, we consider a
multiple-input multiple output (MIMO) joint multiplexing VLC system that
exploits available degrees-of-freedom (DoFs) across space, wavelength and
frequency dimensions simultaneously. Instead of providing a new precoder/post-
detector design, we investigate the considered joint multiplexing system from
a system configuration perspective by tuning system parameters in both spatial
and wavelength domains, such as LED positions and optical filter passband. We
propose a novel spatial clustering with wavelength division (SCWD) strategy
which enhances the MIMO channel condition. We propose to use a state-of-the-
art black-box optimization tool: Bayesian adaptive direct search (BADS) to
determine the desired system parameters, which can significantly improve the
achievable rate. The extensive numerical results demonstrate the superiority
of the proposed method over conventional SMX and WDM VLC systems.

###### Index Terms:

Visible light communication, optical wireless communication, multiple-input
multiple-output, orthogonal frequency division multiplexing, spatial
multiplexing, wavelength division multiplexing.

##  I Introduction

With the development of information technology, an increasing number of
machine-type devices, wireless sensors and cloud services are deployed, which
further increases the demand on wireless network capability [1]. To meet the
need of future wireless services, various new technologies for high speed
wireless transmission have been proposed. visible light communication (VLC)
and light fidelity (LiFi) are among potential candidates [2]. Apart from the
use of the license-free optical spectrum and physical layer security feature,
VLC and LiFi can deliver multi-Gbps transmission data rate [3]. A major
challenge of developing high performance VLC and LiFi systems is the limited
modulation bandwidth of light-emitting diodes and photodiodes. One of the
solutions is to develop optical front-ends with a much wider bandwidth, such
as GaN-based micro-LEDs [4]. An alternative solution is to use spatial
multiplexing (SMX) or wavelength division division (WDM) techniques to
transmit data with multiple parallel channels, which can significantly boost
the aggregate data rate without expanding the modulation bandwidth.

Several early studies on SMX with/without precoding/post-detection have been
reported in [5, 6, 7]. To include the frequency-domain characteristics, VLC
systems with multiple-input multiple-output (MIMO)- orthogonal frequency
division multiplexing (OFDM) have been investigated in [8, 9]. In addition, a
coherent LiFi system with SMX has been considered in a recent study [10].
Regarding research on WDM-based VLC systems, many successful experimental
demonstrations were reported in the literature [11, 12]. An analytical work
shows that the achievable data rate is limited by inter-channel crosstalk when
a large number of wavelength divisions are used [13]. In order to mitigate the
crosstalk, the usage of signal processing in MIMO systems have been proposed
in a few studies [14, 15], where the WDM channels are treated as a colour MIMO
channel matrix and precoding/post-detection processes are used to diagonalise
the multiplexing channel matrix. A recent study shows that it is possible to
implement a WDM system without using optical filters [16].

###  I-A Related research

Regarding the combination of SMX and WDM VLC systems, MIMO-VLC systems with
multi-colour LEDs have been considered to use degree-of-freedoms in both
dimensions [17, 18]. These studies have proposed signal processing techniques
such as optimal precoder designs under lighting constraints [17] or a
chromaticity-adaptive generalised spatial modulation scheme [18]. In these
studies, advanced signal processing techniques are generally designed to
improve communication performance with a given set of system parameters.
Alternatively, the performance of VLC systems can also be improved by changing
system configurations, where key system parameters, such as LED position or
optical filter passband, are carefully selected so that the probability of
improved channel quality is increased. Several studies in multi-cell LiFi
systems investigated the optimal system configurations in terms of access
points spatial deployment and LED parameters so that the system reliability,
spectral efficiency or energy efficiency is maximised [19, 20, 21, 22].
Regarding the research on system configurations in VLC MIMO systems, angular
diversity receivers, mirror diversity receiver and irregular PD configurations
are investigated to improve MIMO channel condition [23, 24, 25]. System
configurations of wavelength domain parameters in WDM VLC systems have been
investigated in [13, 26]. Nevertheless, the above studies consider system
configurations with only a few parameters in either the spatial or wavelength
domain.

###  I-B Motivation and contributions

Due to significant spatial correlation, the number of parallel channels in SMX
VLC systems is limited [25]. On the other hand, it has been shown that the
passband of a thin film optical filter will shift to shorter wavelengths when
the light incident angle is greater than 0°, which causes a severe wavelength
mismatch in a WDM VLC system [26]. By considering a spatial and wavelength
division joint multiplexing VLC system, the MIMO channel with severe spatial
correlation can be decorrelated by the wavelength domain features. In
addition, the excessive inter-colour interference in wavelength domain due to
aforementioned passband shift issue can be mitigated by MIMO precoding and
post-processing blocks. Consequently, the resultant number of parallel
channels in joint multiplexing system can be increased and the corresponding
achievable rates can be improved. Despite the performance improvement from the
novel signal processing techniques, the combining features in spatial and
wavelength domains have not been comprehensively investigated in [17, 18]. The
research findings in [19, 20, 21, 22, 23, 24, 25, 26] also demonstrate the
importance of system configuration in VLC/LiFi systems, which has not been
explored in a spatial and wavelength domain joint multiplexing VLC system yet.
Furthermore, it is complicated to design a VLC multiplexing system using both
spatial and wavelength domain features efficiently. In this paper, a MIMO-OFDM
spatial and wavelength division joint multiplexing VLC system is thoroughly
studied from a system configuration perspective. To evaluate the impact of
various parameters, a detailed framework of a MIMO-OFDM joint multiplexing VLC
system is established considering the characteristics in the spatial,
wavelength and frequency domains, which is unavailable in the literature.
Based on the developed framework, the achievable rates with various system
configurations are evaluated and compared. In particular, a unique spatial
clustering with wavelength division (SCWD) configuration strategy is proposed
which can achieve higher achievable rates compare to the other benchmark
strategies. Furthermore, a Bayesian adaptive direct search (BADS) black-box
optimisation tool has been used to search for system parameters that offer
additional performance improvement. Compared to our previous study [27], a
more detailed system model and more scenarios with practical concerns such as
random user position/orientation are considered. The contributions of this
study are summarised as follows:

  * •

A detailed framework is established for characterising VLC MIMO-OFDM joint
multiplexing systems over space, wavelength and frequency domains. This
framework allows researchers to evaluate the performance of a joint
multiplexing system with a specific system configuration.

  * •

The system configurations of the considered VLC MIMO-OFDM joint multiplexing
system are thoroughly investigated with random user position and device
orientation. Both empirical parameter selections and parameter searching based
on BADS algorithm are considered.

  * •

Based on the idea of ‘division in either the spatial or wavelength domain’, a
SCWD strategy is proposed to efficiently use the DoF in both spatial and
wavelength domains. The performance of the joint multiplexing systems using
the SCWD strategy is compared with benchmark systems, which shows the
superiority of the joint multiplexing system over SMX and WDM techniques in
terms of achievable rate.

  * •

It has been found that systems with the SCWD strategy are superior to SMX and
WDM when two specific conditions are fulfilled: 1. When the achievable rates
of SMX and WDM are similar (with the same number of elements); 2. When the
multiplexing gain improvement from excessively increasing numbers of LEDs/PDs
is saturated for SMX and WDM techniques. Insights into system configuration
strategies and solutions are presented.

The remainder of this paper is arranged as follows. Section II presents the
considered MIMO-OFDM system model. The MIMO channel model considering
characteristics in space, wavelength and frequency domains are introduced in
Section III. The system configuration of the joint multiplexing system is
thoroughly investigated in Section IV. The conclusions are drawn in Section V.

##  II MIMO-OFDM System model

In this section, a conventional MIMO-OFDM VLC system model based on DC-biased
optical (DCO)-OFDM is considered [8]. A block diagram of the system is shown
in Fig. 1. Assuming there are $N_{\rm t}$ LEDs and $N_{\rm r}$ PDs, the
maximum supported number of data streams will be $I\leq\min(N_{\rm t},N_{\rm
r})$. Considering a $K$-point fast Fourier transform (FFT) operation, the
number of subcarriers carrying information bits is $\tilde{K}=K/2-1$. Firstly,
a modulation block maps information bits to $M$-ary quadrature amplitude
modulation (QAM) symbols. Then, after a power control module and a precoding
module, the signal vector on the $k$th subcarrier can be represented by:

$\displaystyle\mathbf{X}_{k}=\mathbf{F}_{k}\mathbf{Q}_{k}^{1/2}\mathbf{S}_{k},\leavevmode\nobreak\
\text{for}\leavevmode\nobreak\ k=1,2,\cdots,\tilde{K},$ (1)

where $\mathbf{S}_{k}\in\mathbb{C}^{I\times 1}$ is the modulated symbol vector
with unit variance; $\mathbf{Q}_{k}^{1/2}\in\mathbb{R}^{I\times I}$ is power
control diagonal matrix and $\mathbf{F}_{k}\in\mathbb{C}^{N_{\rm t}\times I}$
is a precoding matrix, which must be a unity matrix. Note that the QAM symbol
on the $k$th subcarrier of the $i$th data stream is defined as ${S}_{i}[k]$,
which is also the $i$th element of $\mathbf{S}_{k}$. Thus, the signal value on
the $k$th subcarrier for the $n_{\rm t}$th LED can be written as:

$\displaystyle{X}_{n_{\rm t}}[k]=\sum_{i=1}^{I}{F}_{n_{\rm
t},i}[k]\sqrt{q_{i}[k]}{S}_{i}[k],$ (2)

where $\sqrt{q_{i}[k]}$ is the $(i,i)$-th entry of $\mathbf{Q}_{k}^{1/2}$ and
${F}_{n_{\rm t},i}[k]$ is the $(n_{\rm t},i)$-th entry of $\mathbf{F}_{k}$. To
guarantee a real-value time-domain signal, the Hermitian symmetry condition
must be fulfilled, which requires: $X_{n_{\rm t}}[0]=X_{n_{\rm t}}[K/2]=0$ and
$X_{n_{\rm t}}[k]=X_{n_{\rm t}}^{*}[K-k]$ for $k=K/2+1,K/2+2,\cdots,K-1$.
Next, the frequency-domain signal is converted to time-domain for transmission
by using an inverse fast Fourier transform (IFFT) operation:

$\displaystyle x_{n_{\rm t}}[n]=\frac{1}{\sqrt{K}}\sum_{k=0}^{K-1}{X}_{n_{\rm
t}}[k]\mathrm{e}^{\frac{2\pi nkj}{K}},$ (3)

where $j=\sqrt{-1}$ is the imaginary number. Due to the limited dynamic range
of each LED, we normalise and constrain the signal variance not to be greater
than unity: $\mathbb{E}\left\\{x^{2}_{n_{\rm t}}[n]\right\\}\leq 1$, where
$\mathbb{E}\left\\{\cdot\right\\}$ refers to the expectation operator. Note
that the optical power constraint is considered in Section III-A. To avoid
inter-frame interference and inter-symbol interference (ISI), a cyclic prefix
(CP) is added to the beginning of each time-domain OFDM frame. After the
precoding and IFFT operations, a clipping operation is enforced to limit the
signal peak-to-average power ratio (PAPR):

Figure 1: Block diagram of a MIMO-OFDM VLC system based on DCO-OFDM.

$\displaystyle\hat{x}_{n_{\rm t}}[n]=\left\\{\begin{array}[]{lr}\kappa_{\rm
t}&:x_{n_{\rm t}}[n]\geq\kappa_{\rm t}\\\ x_{n_{\rm t}}[n]&:\kappa_{\rm
b}<x_{n_{\rm t}}[n]<\kappa_{\rm t}\\\ \kappa_{\rm b}&:x_{n_{\rm
t}}[n]\leq\kappa_{\rm b}\end{array}\right.,$ (7)

where $\kappa_{\rm t}$ and $\kappa_{\rm b}$ are the top and bottom clipping
levels. According to the Bussgang theorem, the non-linear clipping operation
can be approximated by:

$\displaystyle\hat{x}_{n_{\rm t}}[n]=\eta x_{n_{\rm
t}}[n]+\mathfrak{n}_{n_{\rm t}}^{\rm clip}[n],$ (8)

where $\eta$ is the clipping attenuation factor and $\mathfrak{n}_{n_{\rm
t}}^{\rm clip}[n]$ is the clipping noise at the $n_{\rm t}$th LED [28] which
follows a normal distribution with a zero mean and a variance of $\sigma_{\rm
clip}^{2}$. The value of $\mathfrak{n}_{n_{\rm t}}^{\rm clip}[n]$ and
$\sigma_{\rm clip}^{2}$ can be calculated analytically [29]. Then, the clipped
electric signal is converted to an optical signal and launched to the optical
wireless channel. At the receiver side, a fraction of the optical signals are
detected by the PDs. The detected signal by the $n_{\rm r}$th PD can be
calculated as:

$\displaystyle y_{n_{\rm r}}[n]=\mathfrak{n}^{\rm rx}_{n_{\rm
r}}[n]+\sum_{n_{\rm t}=1}^{N_{\rm t}}\hat{x}_{n_{\rm t}}[n]\otimes h_{n_{\rm
r},n_{\rm t}}[n],$ (9)

where $\mathfrak{n}^{\rm rx}_{n_{\rm r}}[n]$ is the receiver noise, $h_{n_{\rm
r},n_{\rm t}}[n]$ is the discrete channel impulse response between the $n_{\rm
r}$th PD and $n_{\rm t}$th LED and $\otimes$ refers to the convolution
operator. The receiver noise follows a normal distribution with a zero mean
and a variance of ${\sigma}_{{\rm rx},n_{\rm r}}^{2}$. The considered receiver
noise is composed of the shot noise and thermal noise. After the reception of
the signal, the CP is removed. The time-domain signal is converted back to the
frequency-domain using the FFT operation as:

$\displaystyle{Y}_{n_{\rm r}}[k]=\frac{1}{\sqrt{K}}\sum_{k=0}^{K-1}y_{n_{\rm
r}}[n]\mathrm{e}^{-\frac{2\pi nkj}{K}}.$ (10)

The addition of CPs leads to the circular convolution relationship between the
time-domain signal and the channel impulse response. This circular convolution
can be converted to a multiplication relationship in the frequency-domain.
Therefore, the conversion between ${Y}_{n_{\rm r}}[k]$ and ${X}_{n_{\rm
t}}[k]$ can also be evaluated in the frequency-domain directly. The received
signal vector on the $k$th subcarrier $\mathbf{Y}_{k}\in\mathbb{C}^{N_{\rm
r}\times 1}$ can be calculated by:

$\displaystyle\mathbf{Y}_{k}=\mathbf{H}_{k}\left(\eta\mathbf{X}_{k}+\mathbf{E}_{k}\right)+\mathbf{N}_{k},$
(11)

where $\mathbf{H}_{k}\in\mathbb{C}^{N_{\rm r}\times N_{\rm t}}$ is the
frequency-domain channel matrix, $\mathbf{E}_{k}\in\mathbb{C}^{N_{\rm t}\times
1}$ is the frequency-domain clipping noise vector and
$\mathbf{N}_{k}\in\mathbb{C}^{N_{\rm r}\times 1}$ is the frequency-domain
receiver noise vector. The $(n_{\rm r},n_{\rm t})$-th entry of
$\mathbf{H}_{k}$ is denoted as $H_{n_{\rm r},n_{\rm t}}[k]$, the $n_{\rm t}$th
entry of $\mathbf{E}_{k}$ is denoted as $N_{n_{\rm t}}^{\rm clip}[k]$ and the
$n_{\rm r}$th entry of $\mathbf{N}_{k}$ is denoted as $N_{n_{\rm r}}^{\rm
rx}[k]$. In addition, $H_{n_{\rm r},n_{\rm t}}[k]$, $N_{n_{\rm t}}^{\rm
clip}[k]$ and $N_{n_{\rm r}}^{\rm rx}[k]$ are FFT of $h_{n_{\rm r},n_{\rm
t}}[n]$, $\mathfrak{n}^{\rm clip}_{n_{\rm t}}[n]$ and $\mathfrak{n}^{\rm
rx}_{n_{\rm r}}[n]$, respectively. Finally, a MIMO post-detection matrix
$\mathbf{W}_{k}\in\mathbb{C}^{I\times N_{\rm r}}$ is used to retrieve the
transmitted data symbol vectors:

$\displaystyle\mathbf{Z}_{k}$
$\displaystyle=\mathbf{W}_{k}\mathbf{Y}_{k}=\eta\mathbf{W}_{k}\mathbf{H}_{k}\mathbf{F}_{k}\mathbf{Q}_{k}^{1/2}\mathbf{S}_{k}+\mathbf{W}_{k}\mathbf{H}_{k}\mathbf{E}_{k}+\mathbf{W}_{k}\mathbf{N}_{k}.$
(12)

Thus, the data symbol of the $i$th data stream on the $k$th subcarrier can be
written as:

$\displaystyle Z_{i}[k]=\eta\sum_{n_{\rm r}=1}^{N_{\rm r}}\sum_{n_{\rm
t}=1}^{N_{\rm t}}{W}_{i,n_{\rm r}}[k]H_{n_{\rm r},n_{\rm t}}[k]{F}_{n_{\rm
t},i}[k]\sqrt{q_{i}[k]}S_{i}[k]+\sum_{n_{\rm r}=1}^{N_{\rm r}}\sum_{n_{\rm
t}=1}^{N_{\rm t}}{W}_{i,n_{\rm r}}[k]H_{n_{\rm r},n_{\rm t}}[k]{N}^{\rm
clip}_{n_{\rm t}}[k]$ $\displaystyle+\eta\sum_{n_{\rm r}=1}^{N_{\rm
r}}\sum_{n_{\rm t}=1}^{N_{\rm t}}\sum_{\hat{i}=1,\hat{i}\neq
i}^{I}{W}_{i,n_{\rm r}}[k]H_{n_{\rm r},n_{\rm t}}[k]{F}_{n_{\rm
t},\hat{i}}[k]\sqrt{q_{\hat{i}}[k]}S_{\hat{i}}[k]+\sum_{n_{\rm r}=1}^{N_{\rm
r}}{W}_{i,n_{\rm r}}[k]{N}^{\rm rx}_{n_{\rm r}}[k],$ (13)

where ${W}_{i,n_{\rm r}}[k]$ is the $(i,n_{\rm r})$-th entry of
$\mathbf{W}_{k}$, the first term on the right-hand side of the equality is the
desired signal, the second term is the equivalent clipping noise, the third
term is the interference from other multiplexing channels and the last term
corresponds to the equivalent receiver noise. Therefore, the corresponding
signal-to-noise ratio (SNR) of the $i$th data stream on the $k$th subcarrier
can be calculated by:

$\displaystyle\gamma_{i}[k]$ $\displaystyle=\left(\eta^{2}\left|\sum_{n_{\rm
r}=1}^{N_{\rm r}}\sum_{n_{\rm t}=1}^{N_{\rm t}}{W}_{i,n_{\rm r}}[k]H_{n_{\rm
r},n_{\rm t}}[k]{F}_{n_{\rm
t},i}[k]\right|^{2}q_{i}[k]\right)\left(\sigma_{\rm
clip}^{2}\sum\limits_{n_{\rm r}=1}^{N_{\rm r}}\sum\limits_{n_{\rm
t}=1}^{N_{\rm t}}\left|{W}_{i,n_{\rm r}}[k]H_{n_{\rm r},n_{\rm
t}}[k]\right|^{2}\right.$
$\displaystyle\left.+\eta^{2}\sum\limits_{\hat{i}=1,\hat{i}\neq
i}^{I}\left|\sum\limits_{n_{\rm r}=1}^{N_{\rm r}}\sum\limits_{n_{\rm
t}=1}^{N_{\rm t}}{W}_{i,n_{\rm r}}[k]H_{n_{\rm r},n_{\rm t}}[k]{F}_{n_{\rm
t},\hat{i}}[k]\right|^{2}q_{\hat{i}}[k]+\sum\limits_{n_{\rm
r}=1}^{I}\left|{W}_{i,n_{\rm r}}[k]\right|^{2}{\sigma}_{{\rm rx},n_{\rm
r}}^{2}\right)^{-1},$ (14)

In this study, the well-known singular-value decomposition (SVD)-based
precoding and post-detection are used, which remove inter-channel interference
and convert MIMO channels to orthogonal parallel channels.

##  III Spatial and wavelength characteristics of MIMO-OFDM channel

In this section, we introduce the characteristics of the MIMO-OFDM channel
$\mathbf{H}_{k}$ introduced in Section II. The relationship between the
channel and parameters in the space, wavelength and frequency domains is
considered. Firstly, we consider the discrete samples forwarded to the $n_{\rm
t}$th LED and its driving circuit. The discrete samples $\hat{x}_{n_{\rm
t}}[n]$ are converted to a continuous analogue signal via a pulse shaping
process:

$\displaystyle\tilde{x}_{n_{\rm
t}}(t)=\sum_{n=-\infty}^{\infty}\hat{x}_{n_{\rm t}}[n]g(t-nT_{\rm s}),$ (15)

where $g(t)$ is the impulse response of the signal pulse and $T_{\rm s}$ is
the symbol period. Before feeding the analogue signal, $\tilde{x}_{n_{\rm
t}}(t)$ is amplified by a factor of $a_{n_{\rm t}}$ and a direct current
(DC)-bias of $b_{n_{\rm t}}$ is added in the driving circuit. The optical
signal of the $n_{\rm t}$th LED can be written as:

$\displaystyle p_{n_{\rm t}}(t,\lambda)=\mathcal{S}_{n_{\rm t}}^{\rm
led}(\lambda)\left(a_{n_{\rm t}}\tilde{x}_{n_{\rm t}}(t)+b_{n_{\rm
t}}\right)\otimes h^{\rm led}(t),$ (16)

where $\mathcal{S}_{n_{\rm t}}^{\rm led}(\lambda)$ is the normalized spectral
density of the $n_{\rm t}$th LED at wavelength $\lambda$ and $h^{\rm led}(t)$
is the impulse response of the LED. After the emission of the optical signal
to the wireless channel, a fraction of the signal is detected by the PDs on
the receiver side. The output photocurrent of the $n_{\rm r}$th PD can be
written as:

$\displaystyle i_{n_{\rm r}}(t)$ $\displaystyle=\int\limits_{\lambda_{\rm
min}}^{\lambda_{\rm max}}\mathcal{R}^{\rm pd}(\lambda)h^{\rm
pd}(t)\otimes{h}_{n_{\rm r},n_{\rm t}}^{\rm ow}(t,\lambda)\otimes p_{n_{\rm
t}}(t,\lambda)\mathrm{d}\lambda,$ (17)

where $\mathcal{R}^{\rm pd}(\lambda)$ is the spectral PD responsivity, $h^{\rm
pd}(t)$ is the low-pass impulse response of the PD and ${h}_{n_{\rm r},n_{\rm
t}}^{\rm ow}(t,\lambda)$ is the optical wireless channel impulse response.
Since the optical transmission operates in a wide spectrum region, the final
photocurrent is the result of an integration over the involved spectrum
region. Then the photocurrent is forwarded to a matched filter, where the
detected waveform is convolved with the signal pulse $g(t)$ and the discrete
signal is obtained by sampling at $nT_{\rm s}$. By applying a discrete-time
unit impulse function input $\hat{x}_{n_{\rm t}}[n]=\delta[n]$ to (15) and
inserting (15), (16) into (17), the continuous channel impulse response before
sampling can be calculated as:

$\displaystyle h_{n_{\rm r},n_{\rm t}}(t)=\int\limits_{\lambda_{\rm
min}}^{\lambda_{\rm max}}\mathcal{R}^{\rm pd}(\lambda)\mathcal{S}_{n_{\rm
t}}^{\rm led}(\lambda)g(t)\otimes h^{\rm pd}(t)\otimes{h}_{n_{\rm r},n_{\rm
t}}^{\rm ow}(t,\lambda)\otimes h^{\rm led}(t)\otimes\left(a_{n_{\rm
t}}g(t)+b_{n_{\rm t}}\right)\mathrm{d}\lambda.$ (18)

Noting that the convolution operation is with respect to $t$, but the integral
is with respect to $\lambda$. The Fourier transform of (18) can be calculated
as:

$\displaystyle H_{n_{\rm r},n_{\rm t}}(f)=\int\limits_{\lambda_{\rm
min}}^{\lambda_{\rm max}}\mathcal{R}^{\rm pd}(\lambda)\mathcal{S}_{n_{\rm
t}}^{\rm led}(\lambda)G(f)H^{\rm pd}(f){H}_{n_{\rm r},n_{\rm t}}^{\rm
ow}(f,\lambda)H^{\rm led}(f)\left(a_{n_{\rm t}}G(f)+b_{n_{\rm
t}}\delta(f)\right)\mathrm{d}\lambda,$ (19)

where $G(f)$, $H^{\rm pd}(f)$, ${H}_{n_{\rm r},n_{\rm t}}^{\rm ow}(f,\lambda)$
and $H^{\rm led}(f)$ are the FFT of $g(t)$, $h^{\rm pd}(t)$, ${h}_{n_{\rm
r},n_{\rm t}}^{\rm ow}(t,\lambda)$ and $h^{\rm led}(t)$, respectively. Based
on the definition of discrete-time Fourier transform (DTFT), the channel
transfer function after the matched filter sampling at $nT_{\rm s}$ is a
periodic summation of $H_{n_{\rm r},n_{\rm t}}(f)$ with a period of $1/T_{\rm
s}$ as: $H_{n_{\rm r},n_{\rm t}}^{1/T_{\rm
s}}(f)=\sum_{l=-\infty}^{\infty}H_{n_{\rm r},n_{\rm t}}(f-l/T_{\rm s})$. Thus
the channel transfer function between the $n_{\rm t}$th LED and the $n_{\rm
r}$th PD on the $k$th subcarrier can be calculated as:

$\displaystyle{H}_{n_{\rm r},n_{\rm t}}[k]$ $\displaystyle=H_{n_{\rm r},n_{\rm
t}}^{1/T_{\rm s}}\left(\frac{k}{KT_{\rm
s}}\right)=\sum_{l=-\infty}^{\infty}H_{n_{\rm r},n_{\rm
t}}\left(\frac{k}{KT_{\rm s}}-\frac{l}{T_{\rm s}}\right),$ (20)

for $k=0,1,\cdots,K-1$.  Although the limits of the summation in (20) are from
$-\infty$ to $\infty$, the band-limited signal pulse $G(f)$ makes most of the
terms in the summation equal zero. For example, the used root-square raised
cosine (RRC) pulse in this work makes the shifted channel frequency response
(19) equals zero for $f\in\left(-\infty,-\frac{\alpha+1}{2T_{\rm
s}}+\frac{l}{T_{\rm s}}\right)\cup\left(\frac{\alpha+1}{2T_{\rm
s}}+\frac{l}{T_{\rm s}},\infty\right)$. The subcarrier index of the transfer
function falls in the region of $[0,K-1]$, which is within the frequency range
of $[0,1/T_{\rm s})$. This implies that only the terms with $l=0$ and $1$ are
non-zero in (20).

(a)

(b)

Figure 2: (a) Electrical to optical signal conversion. (b) Spatial-domain
characteristics with LoS and NLoS channel components.

###  III-A LED electrical-to-optical conversion

In (16), the choice of $a_{n_{\rm t}}$ and $b_{n_{\rm t}}$ should map the peak
values of the input signal to the peak values of the optical output to
maximise the signal power, as illustrated in Fig. 2(a). Assuming the $n_{\rm
t}$th LED has a maximum and minimum optical levels of $p_{n_{\rm t}}^{\rm
max}$ and $p_{n_{\rm t}}^{\rm min}$, the following mapping equations can be
used: $a_{n_{\rm t}}\kappa_{\rm b}+b_{n_{\rm t}}=p_{n_{\rm t}}^{\rm
min},\leavevmode\nobreak\ a_{n_{\rm t}}\kappa_{\rm t}+b_{n_{\rm t}}=p_{n_{\rm
t}}^{\rm max}$, which lead to:

$\displaystyle a_{n_{\rm t}}=\frac{p_{n_{\rm t}}^{\rm max}-p_{n_{\rm t}}^{\rm
min}}{\kappa_{\rm t}-\kappa_{\rm b}},\leavevmode\nobreak\ b_{n_{\rm
t}}=\frac{p_{n_{\rm t}}^{\rm min}\kappa_{\rm t}-p_{n_{\rm t}}^{\rm
max}\kappa_{\rm b}}{\kappa_{\rm t}-\kappa_{\rm b}}.$ (21)

By using (16), the average optical power of the $n_{\rm t}$th LED can be
calculated by:

$\displaystyle\bar{p}_{n_{\rm
t}}=\lim_{T\to\infty}\int\limits_{-T}^{T}\int\limits_{\lambda_{\rm
min}}^{\lambda_{\rm max}}\frac{p_{n_{\rm
t}}(t,\lambda)}{2T}\mathrm{d}\lambda\mathrm{d}t=a_{n_{\rm
t}}\mathbb{E}\left\\{\hat{x}_{n_{\rm t}}[n]\right\\}+b_{n_{\rm t}},$ (22)

where $\mathbb{E}\left\\{\hat{x}_{n_{\rm t}}[n]\right\\}$ is the expectation
of the clipped signal, which can be evaluated analytically [29]. For
simplicity of the analysis, we consider a symmetric clipping ($\kappa_{\rm
t}=\kappa$ and $\kappa_{\rm b}=-\kappa$) and a zero minimum optical power
level ($p_{n_{\rm t}}^{\rm min}=0$). This leads to simplification to (21) and
(22) as $a_{n_{\rm t}}=p_{n_{\rm t}}^{\rm max}/2\kappa$, $b_{n_{\rm
t}}=p_{n_{\rm t}}^{\rm max}/2=\bar{p}_{n_{\rm t}}$ and $\bar{p}_{n_{\rm
t}}=b_{n_{\rm t}}$. If we consider $\bar{p}_{n_{\rm t}}$ as the given
parameter, the scaling factor and DC-bias can be calculated by $a_{n_{\rm
t}}=\bar{p}_{n_{\rm t}}/\kappa$ and $\bar{p}_{n_{\rm t}}$.

###  III-B Spatial-domain characteristics

The spatial-domain characteristics are primarily determined by the optical
wireless channel between the LEDs and the PDs, which can be decomposed into a
line-of-sight (LoS) component and non-line-of-sight (NLoS) components, as
shown in Fig. 2(b). The corresponding transfer function can be defined as:
${H}_{n_{\rm r},n_{\rm t}}^{\rm ow}(f,\lambda)={H}_{n_{\rm r},n_{\rm t}}^{\rm
LoS}(f,\lambda)+{H}_{n_{\rm r},n_{\rm t}}^{\rm NLoS}(f,\lambda)$. The LoS
corresponds to the signal propagation directly from the LEDs to the PDs, which
dominates the MIMO channel in most cases. The frequency response can be
calculated as [30]:

$\displaystyle{H}_{n_{\rm r},n_{\rm t}}^{\rm
LoS}(f,\lambda)=\frac{\left(m_{\rm led}+1\right)A_{\rm pd}\mathbf{1}_{\rm
v}}{2\pi D_{n_{\rm r},n_{\rm t}}^{2}}\exp(-j2\pi f\tau_{n_{\rm r},n_{\rm
t}})\mathcal{G}_{n_{\rm r}}^{\rm of}(\lambda,\psi_{n_{\rm r},n_{\rm
t}})\cos^{m_{\rm led}}\phi_{n_{\rm r},n_{\rm t}}\cos^{m_{\rm fov}}\psi_{n_{\rm
r},n_{\rm t}},$ (23)

where $A_{\rm pd}$ is the active area of the PD, $m_{\rm led}$ is the
Lambertian emission order of the LED, $m_{\rm fov}$ is the field-of-view (FoV)
coefficient [23], $\mathbf{1}_{\rm v}$ is a visibility function,
$\mathcal{G}_{n_{\rm r}}^{\rm of}(\lambda,\psi)$ is the transmittance of the
optical filter mounted on the $n_{\rm r}$th PD and $D_{n_{\rm r},n_{\rm t}}$,
$\phi_{n_{\rm r},n_{\rm t}}$, $\psi_{n_{\rm r},n_{\rm t}}$, $\tau_{n_{\rm
r},n_{\rm t}}$ are the Euclidean distance, radiant angle, incident angle, time
delay between the $n_{\rm t}$th LED and the $n_{\rm r}$th PD, respectively.
The Lambertian emission order $m_{\rm led}$ is related to the LED half-power
seminangle by $m_{\rm led}=-1/\log_{2}(\cos(\phi_{1/2}))$. The time delay can
be calculated by $\tau_{n_{\rm r},n_{\rm t}}=D_{n_{\rm r},n_{\rm t}}/c$, where
$c=3\times 10^{8}$ m/s is the speed of light. The visibility function is
defined as:

$\displaystyle\mathbf{1}_{\rm
v}=\left\\{\begin{array}[]{lr}1&:\phi<\pi/2\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ \psi<\pi/2\\\
0&:\text{otherwise}\end{array}\right.,$ (26)

which forces the channel to be zero when either $\phi$ or $\psi$ exceed
$pi/2$. The value of (23) is directly determined by the positions and
orientations of the LEDs and PDs. The trigonometric and distance terms in (23)
can be calculated by [31]:

$\displaystyle D_{n_{\rm r},n_{\rm t}}=\left|\left|\mathbf{p}_{n_{\rm t}}^{\rm
led}-\mathbf{p}_{n_{\rm r}}^{\rm pd}\right|\right|,\leavevmode\nobreak\
\cos\phi_{n_{\rm r},n_{\rm t}}=\frac{\mathbf{o}_{n_{\rm t}}^{\rm
led}}{D_{n_{\rm r},n_{\rm t}}}\cdot(\mathbf{p}_{n_{\rm r}}^{\rm
pd}-\mathbf{p}_{n_{\rm t}}^{\rm led}),\leavevmode\nobreak\ \cos\psi_{n_{\rm
r},n_{\rm t}}=\frac{\mathbf{o}_{n_{\rm r}}^{\rm pd}}{D_{n_{\rm r},n_{\rm
t}}}\cdot(\mathbf{p}_{n_{\rm t}}^{\rm led}-\mathbf{p}_{n_{\rm r}}^{\rm pd}),$
(27)

where $\mathbf{p}_{n_{\rm t}}^{\rm led}$ and $\mathbf{o}_{n_{\rm t}}^{\rm
led}$ are the position and orientation vectors of the $n_{\rm t}$th LED,
respectively; $\mathbf{p}_{n_{\rm r}}^{\rm pd}$ and $\mathbf{o}_{n_{\rm
r}}^{\rm pd}$ are the position and orientation vectors of the $n_{\rm r}$th
PD, respectively; $\\{\cdot\\}$ refers to the vector dot product and
$||\cdot||$ refers to the Euclidean norm. Note that the characteristics of
$\mathcal{G}_{n_{\rm r}}^{\rm of}(\lambda,\psi)$ are affected by features in
both the spatial and wavelength domains, which will be covered in Section
III-C.  The NLoS channel responses correspond to the signal propagation via
reflections by the room internal surfaces, which can be evaluated using an
efficient frequency-domain simulation method [30].

###  III-C Wavelength-domain characteristics

Regarding the wavelength-domain characteristics, analytical spectrum models
are used to improve the flexibility to configure the joint multiplexing
system. The LED normalised spectral intensity can be defined by [13]:

$\displaystyle\mathcal{S}_{n_{\rm t}}^{\rm
led}(\lambda)=\frac{\frac{2}{\sqrt{\pi}}\exp\left(-\frac{\left(\lambda-\lambda_{n_{\rm
t}}^{\rm
led,c}\right)^{2}}{\Delta\lambda_{0.5}^{2}}\right)+\frac{4}{\sqrt{\pi}}\exp\left(-\frac{5(\lambda-\lambda_{n_{\rm
t}}^{\rm
led,c})^{2}}{\Delta\lambda_{0.5}^{2}}\right)}{\Delta\lambda_{0.5}\left(\frac{2+\sqrt{5}}{\sqrt{5}}+\mathrm{erf}\left(\frac{\lambda_{n_{\rm
t}}^{\rm
led,c}}{\Delta\lambda_{0.5}}\right)+\frac{2}{\sqrt{5}}\mathrm{erf}\left(\frac{\sqrt{5}\lambda_{n_{\rm
t}}^{\rm led,c}}{\Delta\lambda_{0.5}}\right)\right)},$ (28)

where $\lambda_{n_{\rm t}}^{\rm led,c}$ is the central wavelength of the
$n_{\rm t}$th LED and $\Delta\lambda_{0.5}$ is a parameter determining the
spectrum shape of the LED, which is defined as:

$\displaystyle\Delta\lambda_{0.5}=\left\\{\begin{array}[]{lr}\frac{5.5\mathcal{K}_{\rm
B}T_{\rm j}}{\mathfrak{h}c}\left(\lambda_{n_{\rm t}}^{\rm
led,c}\right)^{2}&:\lambda_{n_{\rm t}}^{\rm led,c}\leq 560\leavevmode\nobreak\
\text{nm}\\\ \frac{2.5\mathcal{K}_{\rm B}T_{\rm
j}}{\mathfrak{h}c}\left(\lambda_{n_{\rm t}}^{\rm
led,c}\right)^{2}&:\lambda_{n_{\rm t}}^{\rm led,c}>560\leavevmode\nobreak\
\text{nm}\end{array}\right.,$ (31)

where $\mathcal{K}_{\rm B}=1.38\times 10^{-23}$ J/K is the Boltzmann’s
constant, $T_{\rm j}=300$ K is the active layer temperature and
$\mathfrak{h}=6.63\times 10^{-34}$ J/Hz is Planck’s constant. Note that
$\int_{\lambda_{\rm min}}^{\lambda_{\rm max}}\mathcal{S}^{\rm led}_{n_{\rm
t}}(\lambda)\mathrm{d}\lambda=1$. This model has been demonstrated to be
accurate compared to the off-the-shelf LED devices [13]. The PD spectral
responsivity can be defined by the following expression [32]:

$\displaystyle\mathcal{R}^{\rm pd}(\lambda)=\frac{\eta_{\rm
q}\mathfrak{q}\lambda}{\mathfrak{h}c},$ (32)

where $\eta_{\rm q}$ is the quantum efficiency and $\mathfrak{q}=1.6\times
10^{-19}$ C is the electric charge. Considering a thin-film optical bandpass
filter mounted on the $n_{\rm r}$th PD with a central passband wavelength of
$\lambda^{\rm of,c}_{n_{\rm r}}$ and a passband width of $\Delta\lambda^{\rm
of}$ at $0$°incident angle, the spectral transmittance of the filter can be
modelled by [26]:

$\displaystyle\mathcal{G}_{n_{\rm r}}^{\rm
of}(\lambda,\psi)=\left\\{\begin{array}[]{lr}\mathcal{G}_{\rm
T}&:\lambda_{n_{\rm r}}^{\rm of,l}(\psi)\leq\lambda\leq\lambda_{n_{\rm
r}}^{\rm of,r}(\psi)\\\ 0&:\text{otherwise}\end{array}\right.,$ (35)

where $\mathcal{G}_{\rm T}$ is the transmittance of the optical filter,
$\lambda_{n_{\rm r}}^{\rm of,l}(\psi)$ and $\lambda_{n_{\rm r}}^{\rm
of,r}(\psi)$ are left and right edges of the filter passband, which are
functions of the incident angle $\psi$:

$\displaystyle\lambda_{n_{\rm r}}^{\rm of,l}(\psi)=\left(\lambda^{\rm
of,c}_{n_{\rm r}}-\Delta\lambda^{\rm
of}/2\right)\sqrt{1-\sin^{2}\psi/\mathfrak{n}_{\rm e}^{2}},$ (36)
$\displaystyle\lambda_{n_{\rm r}}^{\rm of,r}(\psi)=\left(\lambda^{\rm
of,c}_{n_{\rm r}}+\Delta\lambda^{\rm
of}/2\right)\sqrt{1-\sin^{2}\psi/\mathfrak{n}_{\rm e}^{2}},$ (37)

and $\mathfrak{n}_{\rm e}$ is defined as the effective refraction index.
Several examples of (28), (32) and (35) have been depicted in Fig. 3(a). It is
worth noting that the passbands of the four plotted optical filters shift to
shorter wavelengths significantly when the light incident angle changes from
$\psi=0\degree$ to $\psi=60\degree$. In indoor VLC/LiFi applications, desired
detector alignment is unlikely. Therefore, it is important to consider this
bandpass shift characteristic.

(a)

(b)

Figure 3: (a) Illustration of wavelength-dependent quantities. 1) incident
angle $\psi=0\degree$ 2) incident angle $\psi=60\degree$. The central
wavelength of the LED spectra and optical filter passbands at 0°incident angle
are 437.5 nm, 512.5 nm, 587.5 nm and 662.5 nm. The optical filter passband
width is 75 nm. (b) 1) UE position and horizontal orientation geometry. 2) UE
rotation about $\mathtt{z}$-axis. 3) UE rotation about $\mathtt{x}$-axis. 4)
UE rotation about $\mathtt{y}$-axis.

###  III-D Receiver noise model

Assuming that the background light power is negligible compared to the optical
signal power, the dominant receiver noise components are the signal-dependent
shot noise and the thermal noise. Therefore, the receiver noise variance can
be calculated as:

$\displaystyle{\sigma}_{{\rm rx},n_{\rm r}}^{2}=2\mathfrak{q}\tilde{i}_{n_{\rm
r}}B_{\rm s}+\frac{4\mathcal{K}_{\rm B}\mathcal{T}_{\rm a}B_{\rm s}}{R_{\rm
L}},$ (38)

where $\tilde{i}_{n_{\rm r}}$ is the photocurrent of the $n_{\rm r}$th PD,
$B_{\rm s}=1/2T_{\rm s}$ is the signalling bandwidth and $\mathcal{T}_{\rm a}$
is the absolute temperature. For simplicity, the average photocurrent is used
for the value of $\tilde{i}_{n_{\rm r}}$.

##  IV  Performance evaluation with various system configurations

In Section III, it has been shown that the MIMO-OFDM channel is determined by
many different parameters in the spatial and wavelength domains. In this
section, we focus on exploring various system configurations. Specifically, we
evaluate the performance of joint multiplexing systems with different sets of
spatial and wavelength domain parameters, such as transmitter position
$\mathbf{p}_{n_{\rm t}}^{\rm led}$ or optical filter passband
$\Delta\lambda^{\rm of}$.  In a multiplexing system, the number of
multiplexing channels has a significant impact on the aggregate data rate,
which is determined by $\min(N_{\rm t},N_{\rm r})$. When evaluating the
performance of joint multiplexing systems scaling with the number of LEDs/PDs,
the data rate varies slightly if only one of the variables (number of LEDs or
PDs) is changed. To avoid generating trivial results that are less meaningful,
the numbers of LEDs and PDs are always identical ($N_{\rm t}=N_{\rm r}=I$) in
the considered MIMO-OFDM systems in the following sections. For the
convenience of description, we define a variable called the ‘number of
elements’ which is equivalent to the number of LEDs, PDs and data streams. In
addition, to reduce the complexity of the metric evaluation, a uniform power
allocation is used in this section: $q_{i}[k]=K/(K-2)$.

###  IV-A Performance metric: average achievable rate with random user
equipment (UE) position and orientation

Firstly, we define the average achievable rate as the performance metric.
Since the calculated SNR in (14) corresponds to complex bipolar QAM symbols
after a series of conversions, the achievable rate can be evaluated by [33]:

$\displaystyle C=\frac{1}{T_{\rm s}\left(K+N_{\rm
cp}\right)}\sum_{i=1}^{I}\sum_{k=1}^{\tilde{K}}\log_{2}\left(1+\frac{\gamma_{i}[k]}{\Gamma}\right),$
(39)

where $\gamma_{i}[k]$ is the SNR and $\Gamma$ is a gap factor between the
achievable rate and the Shannon capacity, which is empirically configured to
compensate the system impairment and imperfect constellation [33]. Most of the
parameters related to the SNR calculation are determined by the system
configuration only, such as spectrum intensity or radiation pattern of an LED.
However, the position and orientation of PDs are determined by the user
behaviour. Therefore, the randomness in UE position and orientation is
considered in the evaluation of the achievable rate.

The UE geometry is depicted in Fig. 3(b) 1). The three-dimensional (3D)
position vector of the UE is defined as $\mathbf{p}_{\rm rx}=[\mathtt{x}_{\rm
rx},\mathtt{y}_{\rm rx},\mathtt{z}_{\rm rx}]$, where the $\mathtt{x}_{\rm rx}$
and $\mathtt{y}_{\rm rx}$ are the two-dimensional (2D) coordinate of the UE
and $\mathtt{z}_{\rm rx}$ is the UE height from the floor level, which is
assumed to be deterministic in this study. It is also assumed that the size of
the receiver is small, relative to the link distance, so that the multiple PDs
mounted on the receiver are collocated. Thus, the position vector of the
$n_{\rm r}$th PD is the same as the UE position: $\mathbf{p}_{n_{\rm r}}^{\rm
pd}=\mathbf{p}_{\rm rx}$. In addition, an azimuth angle $\Omega$ is defined to
determine the horizontal direction of the UE.  Since there is no bias on the
user position and the azimuth angle, it is assumed that these quantities are
uniformly distributed, which have been widely accepted in the wireless
communication research community: $\mathtt{x}_{\rm
rx}\sim\mathcal{U}_{(0,W_{\rm room})}$, $\mathtt{y}_{\rm
rx}\sim\mathcal{U}_{(0,L_{\rm room})}$ and $\Omega\sim\mathcal{U}_{(0,2\pi]}$,
where $\mathcal{U}_{\mathcal{I}}$ refers to the continuous uniform
distribution with an interval $\mathcal{I}$. The notations $L_{\rm room}$ and
$W_{\rm room}$ correspond to the length and width of the room. In addition to
the random user position and horizontal direction, two user device orientation
scenarios are considered in this study. The random orientation of the UE is
defined by three rotation angles: $\beta_{\mathtt{z}}$, $\beta_{\mathtt{x}}$
and $\beta_{y}$ which correspond to the yaw, pitch and roll rotations about
the $\mathtt{z}$, $\mathtt{x}$ and $\mathtt{y}$-axes, as illustrated in Fig.
3(b) 2), 3) and 4). Then, the PD orientation vector $\mathbf{o}_{n_{\rm
r}}^{\rm pd}$ is determined by:

$\displaystyle{\left[\mathbf{o}_{n_{\rm r}}^{\rm
pd}\right]}^{\mathrm{T}}=\mathbf{R}_{\mathtt{z}}\mathbf{R}_{\mathtt{x}}\mathbf{R}_{\mathtt{y}}\left[\mathbf{o}_{n_{\rm
r}}^{\rm pd,\uparrow}\right]^{\mathrm{T}},$ (40)

where $\mathbf{o}_{n_{\rm r}}^{\rm pd,\uparrow}$ is the PD orientation vector
of the $n_{\rm r}$th PD when the UE has an upward facing orientation with the
normal vector of $\mathbf{o}_{\rm rx}=[0,0,1]$; $\mathbf{R}_{\mathtt{z}}$,
$\mathbf{R}_{\mathtt{x}}$ and $\mathbf{R}_{\mathtt{y}}$ are 3D rotation
matrices corresponding to $\beta_{\mathtt{z}}$, $\beta_{\mathtt{x}}$ and
$\beta_{y}$, respectively. They are defined by:

Table I: Upward and random orientations rotation angle characteristics | (a) upward | (b) random
---|---|---
| value |  mean $\mu$ |  standard deviation $\sqrt{2}b$
$\beta_{\mathtt{z}}$ | $\Omega-90\degree$ | $\Omega-90\degree$ | $3.67\degree$
$\beta_{\mathtt{x}}$ | $0\degree$ | $40.78\degree$ | $2.39\degree$
$\beta_{\mathtt{y}}$ | $0\degree$ | $-0.84\degree$ | $2.21\degree$

$\displaystyle\mathbf{R}_{\mathtt{z}}=\left[\begin{matrix}\cos\beta_{\mathtt{z}}&-\sin\beta_{\mathtt{z}}&0\\\
\sin\beta_{\mathtt{z}}&\cos\beta_{\mathtt{z}}&0\\\
0&0&1\end{matrix}\right],\leavevmode\nobreak\
\mathbf{R}_{\mathtt{x}}=\left[\begin{matrix}1&0&0\\\
0&\cos\beta_{\mathtt{x}}&-\sin\beta_{\mathtt{x}}\\\
0&\sin\beta_{\mathtt{x}}&\cos\beta_{\mathtt{x}}\end{matrix}\right],\leavevmode\nobreak\
\mathbf{R}_{\mathtt{y}}=\left[\begin{matrix}\cos\beta_{\mathtt{y}}&0&\sin\beta_{\mathtt{y}}\\\
0&1&0\\\ -\sin\beta_{\mathtt{y}}&0&\cos\beta_{\mathtt{y}}\end{matrix}\right].$
(41)

In the first scenario, an upward facing UE is considered, which corresponds to
large devices such as laptops. In this case, we have
$\beta_{\mathtt{z}}=\Omega-90\degree$ and
$\beta_{\mathtt{x}}=\beta_{\mathtt{y}}=0\degree$, as concluded in Table I (a).
In the second scenario, hand-held devices are considered with a random
orientation. The statistics of $\beta_{\mathtt{z}}$, $\beta_{\mathtt{x}}$ and
$\beta_{\mathtt{y}}$ have been obtained via experimental measurement and
modelled using Laplace and Gaussian distributions [34]. For sitting users, the
three rotation angles follow Laplace distributions:
$\beta\sim\mathcal{L}(\mu,b)$ with a mean value of $\mu$ and a standard
deviation of $\sqrt{2}b$. The means and standard deviations of different
rotation angles are concluded in Table I (b). By incorporating the UE
randomness, the average achievable rate can be evaluated by:

$\displaystyle\bar{C}=\mathbb{E}_{\mathtt{x}_{\rm rx},\mathtt{y}_{\rm
rx},\Omega}\left[C\right],\leavevmode\nobreak\
\bar{C}=\mathbb{E}_{\mathtt{x}_{\rm rx},\mathtt{y}_{\rm
rx},\Omega,\beta_{\mathtt{z}},\beta_{\mathtt{x}},\beta_{\mathtt{y}}}\left[C\right],$
(42)

for the cases of upward and random orientation scenarios, respectively. Due to
the complexity issues, these metrics have to be evaluated using a Monte Carlo
approach.

###  IV-B System configuration with BADS algorithm

In this paper, the average achievable rates (42) with several sets of
empirically selected parameters are evaluated. In addition, we aim to evaluate
the maximum potential of the proposed joint multiplexing system in terms of
achievable rate. With the number of parameters exceeding ten, it is very time-
consuming to find the best system configuration parameters with an exhaustive
search approach by testing all combinations of different parameter values.
Therefore, a suitable optimisation tool is used to efficiently search the
desired system configuration.

The optimisation objective function (42) is a complex function of the
considered system parameters by equations from (14) to (39). Additionally, the
numerical evaluation of (42) requires averaging over random user orientations
and positions with a Monte-Carlo approach. This leads to a noisy, high
dimensional and high complexity objective function, which cannot be solved by
a conventional gradient-based optimisation method, such as linear programming.
Instead, a state-of-the-art black-box optimisation tool is used to tackle this
problem. Black-box optimisation methods are used when the objective function
is unknown, non-smooth or complicated to evaluate, where the derivate of the
objective function is unavailable, unreliable or impractical to calculate
[35]. The used optimisation is known as BADS which combines the capability of
Bayesian optimisation in optimising expensive and noisy black-box functions
with the low computational cost of mesh adaptive direct search (MADS). It has
been demonstrated that the BADS algorithm outperforms a number of widely used
and state-of-the-art non-convex, derivative-free optimisation algorithms on
many practical problems [36]. This type of optimisation methods have wide
applications in the fields of computational neuroscience and machine learning
where models are evaluated via stochastic simulation or numerical
approximation [37].

Initialise: $n_{\rm iter}=0$, $\Delta^{\rm m}_{0}=2^{-10}$, $\Delta^{\rm
p}_{0}=1$, evaluate $\mathfrak{f}(\mathbf{x}_{0})$

while  _$n_{\rm iter}\leq N_{\rm iter,max}$ and $\Delta^{\rm
p}_{0}>\Delta^{\rm p}_{\rm th}$_ do

repeat

generate search input $\mathbf{x}_{\rm search}$ and evaluate
$\mathfrak{f}(\mathbf{x}_{\rm search})$

until  _$\mathfrak{f}(\mathbf{x}_{\rm search})$ provides suffient improvement
or reaches maximum search number $N_{\rm search}$_;

if  _Search stage is unsuccessful_ then

generate poll set $P_{n_{\rm iter}}$, sort it using the acquisition function
and evaluate $\mathfrak{f}(\mathbf{x})$ on $\mathbf{x}\in P_{n_{\rm iter}}$

end if

if  _$n_{\rm iter}$ th iteration is successful_ then

update incumbent $\mathbf{x}_{n_{\rm iter}+1}$

if  _Poll stage is successful_ then

$\Delta^{\rm m}_{n_{\rm iter}+1}=2\Delta^{\rm m}_{n_{\rm iter}}$, $\Delta^{\rm
p}_{n_{\rm iter}+1}=2\Delta^{\rm p}_{n_{\rm iter}}$

end if

else

$\Delta^{\rm m}_{n_{\rm iter}+1}=\frac{1}{2}\Delta^{\rm m}_{n_{\rm iter}}$,
$\Delta^{\rm p}_{n_{\rm iter}+1}=\frac{1}{2}\Delta^{\rm p}_{n_{\rm iter}}$

end if

$n_{\rm iter}=n_{\rm iter}+1$

end while

Algorithm 1 Bayesian Adaptive Direct Search [36]

The major operations in the BADS algorithm are introduced in Algorithm 1. The
BADS algorithm considers a possibly noisy objective function
$\mathfrak{f}(\mathbf{x})$ with an $N_{\rm d}$ dimension input variable vector
$\mathbf{x}\in\mathbb{R}^{N_{\rm d}}$. A lower bound vector $\mathbf{x}_{\rm
LB}$ and an upper bound vector $\mathbf{x}_{\rm UB}$ must be provided as the
constraints to the variable vector $\mathbf{x}$. An optional plausible bound
vectors $\mathbf{x}_{\rm PLB}$ and $\mathbf{x}_{\rm PUB}$ can also be provided
to indicate the regions where the optimal solutions are more likely. At the
beginning of the algorithm, an initial solution $\mathbf{x}_{0}$ will be
provided as a starting point. Then, the algorithm enters the search stage and
executes a series of efficient local Bayesian optimisations around the point
and try to find a solution which is better than the current one. Each new
evaluation of the cost function is used as a new sample of the data set to re-
train a Gaussian process model to be a more accurate surrogate. If the
algorithm fails to find a better solution in the search stage, it means the
current Gaussian process model is not useful in the optimisation. Then, the
algorithm can switch to the poll stage which uses model-free opportunistic
optimisation to compensate. The notations $\Delta^{\rm m}_{n_{\rm iter}}$ and
$\Delta^{\rm p}_{n_{\rm iter}}$ are the mesh and poll size of the $n_{\rm
iter}$th iteration. They function as the adaption step sizes, which vary in
the optimisation process. In unsuccessful iterations, $\Delta^{\rm m}_{n_{\rm
iter}}$ and $\Delta^{\rm p}_{n_{\rm iter}}$ are decreased to increase the
searching granularity. In successful iteration in the poll stage, they are
increased to speed up the convergence speed. Finally, the algorithm ends when
the poll size is smaller than a predefined threshold or the maximum number of
iterations is reached.

More details about the implementation, convergence analysis, optimality
condition and complexity of BADS can be found in [36, 38]. In this study, we
have used the BADS algorithm implementation in MATLAB provided by the authors
of [36]. It is worth noting that the BADS algorithm is a semi-local algorithm
which may generate a global or a local optimal solution. To increase the
chance of obtaining global optimal solutions, a multi-start strategy is used
and the number of iterations with different initial solutions for each
optimisation problem equals or greater than ten. In this study, the running
time of executing the BADS algorithm with multiple iterations varies from a
few seconds to tens of minutes with the increase of element number $I$.
However, the obtained system configuration parameters from searching the
problems with the BADS algorithm are static and used to implement the VLC
joint multiplexing system. These obtained parameters are fixed once the system
is implemented and will not change during the operation of the system.
Therefore, unlike many signal processing algorithms, there is little
computational complexity concerns about the use of BADS algorithm in this
work.

(a)

(b)

Figure 4: (a) Illustration of various LED layouts on the ceiling and PD
orientations on the mobile side with different configuration strategies: 1)
Division in space domain. 2) Division in wavelength domain. 3) Division in
both spatial and wavelength domains. 4) Division in either spatial or
wavelength domain. (b) Angular diversity receivers: 1), 2) and 3) show three
PR examples with $N_{\rm r}=4,8,16$. 4), 5) and 6) show three HR examples with
$N_{\rm r}=4,8,16$.

###  IV-C Overview of system configuration strategies

Now we consider four potential configuration strategies: 1) division in the
space domain, 2) division in the wavelength domain, 3) division in both the
spatial and wavelength domains, 4) division in either the spatial or
wavelength domain, as demonstrated in Fig. 4(a). The first two strategies use
the DoFs in the spatial or wavelength domain only, which are equivalent to SMX
and WDM transmission systems, respectively. These two strategies are referred
as space division (SD) and wavelength division (WD). They will be covered in
Section IV-D and IV-E, respectively. The third and fourth strategies use the
DoFs in the spatial and/or wavelength domain jointly. The third strategy
considers each optical element has a distinct spatial and wavelength domain
feature, as illustrated in Fig. 4(a) 3). This strategy leads to excessive
decorrelation of the channel matrix at the cost of inefficient use of spatial
and wavelength DoFs. In contrast, the fourth strategy considers to let each
optical element to have a unique feature in either the spatial or wavelength
domain, as illustrated in Fig. 4(a) 4). It is expected to provide enough
channel matrix decorrelation and also be able to support a greater number of
multiplexing channels compared to the SD and WD strategies. This joint
strategy will be covered in Section IV-F.

###  IV-D Spatial division strategy

Firstly, the configurable parameters considered in SD strategy is introduced.
Since only the DoF in spatial domain is utilised for multiplexing
transmission, wavelength domain parameters are empirically selected with
identical values for all LEDs and optical filters: $\Delta\lambda^{\rm
of}=300$ nm, $\lambda_{n_{\rm t}}^{\rm led,c}=550$ nm for $n_{\rm
t}=1,2,\cdots,N_{\rm t}$ and $\lambda_{n_{\rm r}}^{\rm of,c}=550$ nm for
$n_{\rm r}=1,2,\cdots,N_{\rm r}$. On the transmitter side, the position vector
of the $n_{\rm t}$th LED is defined as:

$\displaystyle\mathbf{p}_{n_{\rm t}}^{\rm led}=[\mathtt{x}_{n_{\rm
t}},\mathtt{y}_{n_{\rm t}},\mathtt{z}_{\rm tx}],$ (43)

where $\mathtt{z}_{\rm tx}$ is the height coordinate of the LEDs determined by
the height of the room; $\mathtt{x}_{n_{\rm t}}$ and $\mathtt{y}_{n_{\rm t}}$
are the 2D horizontal coordinates of the LED, which are configurable
parameters. Regarding the orientation of LEDs, a straight downward orientation
of $\mathbf{o}_{n_{\rm t}}^{\rm led}=[0,0,-1]$ is used to guarantee the
lighting performance and a wide coverage. Finally, the radiation pattern is
defined by the LED half-power semiangle $\phi_{1/2}$, which is another
important configurable parameter. On the receiver side, it is impractical to
control the exact position and orientation of each PD due to the random user
location and orientation, while the FoV coefficient $m_{\rm fov}$ can be
manipulated to control the light reception pattern, which is a configurable
parameter of interest. On the other hand, the small UE size assumption
prevents the effective spatial decorrelation using different PD positions. In
order to achieve low spatial correlation on the receiver side, two types of
angular diversity receivers are considered in this study: pyramid receiver
(PR) and hemispheric receiver (HR) [23]. In equation (40), the upward facing
PD orientation vector can be represented by $\mathbf{o}_{n_{\rm r}}^{\rm
pd,\uparrow}=\left[\cos\omega_{n_{\rm r}}\sin\theta_{n_{\rm
r}},\sin\omega_{n_{\rm r}}\sin\theta_{n_{\rm r}},\cos\theta_{n_{\rm
r}}\right]$, where $\omega_{n_{\rm r}}$ is a PD azimuth angle and
$\theta_{n_{\rm r}}$ is a PD elevation angle. In a PR, the values of
$\theta_{n_{\rm r}}$ is defined by a PR elevation angle $\theta_{\rm pd}$,
which is a suitable configurable parameter. The PD azimuth angles are
specified as:

$\displaystyle\omega_{n_{\rm r}}=\frac{2\pi}{N_{\rm r}}(n_{\rm r}-1).$ (44)

In a HR, the elevation angles are specified as:

$\displaystyle\theta_{n_{\rm r}}=\arccos\left(s_{n_{\rm r}}\right),$ (45)

with $s_{n_{\rm r}}=1-\frac{2(n_{\rm r}-1)}{2N_{\rm r}-1}$, and the azimuth
angles are specified as:

$\displaystyle\omega_{n_{\rm r}}=\left(\omega_{n_{\rm r}-1}+3.6\left(2N_{\rm
r}\left(1-s_{n_{\rm
r}}^{2}\right)\right)^{-1/2}\right)\mathrm{mod}\leavevmode\nobreak\ 2\pi,$
(46)

for $n_{\rm r}=2,\cdots,N_{\rm r}$ and $\omega_{1}=0$. Fig. 4(b) demonstrate
six examples of PRs and HRs with $N_{\rm r}=4,8,16$. In summary, the
considered configurable parameters in the SD strategy include:
$\mathtt{x}_{n_{\rm t}}$, $\mathtt{y}_{n_{\rm t}}$, $\phi_{1/2}$, $m_{\rm
fov}$, $(\theta_{\rm pd})$. Note that $\theta_{\rm pd}$ is only considered
when PRs are used.

Figure 5: Empirical configuration to 2D LED position layout.

####  IV-D1 Empirical configuration

In cases with the empirical configuration, the values of $\mathtt{x}_{n_{\rm
t}}$ and $\mathtt{y}_{n_{\rm t}}$ with different number of $N_{\rm t}$ are
specified according to the 2D layouts shown in Fig. 5. The intuition of this
LED spatial distribution is to distribute the LEDs with wider separations to
guarantee a low spatial correlation on the transmitter side. On the other
hand, the LEDs are distributed in a manner to let the system cover the entire
room. An LED halfpower semiangle of $\phi_{1/2}=60\degree$ is used to ensure a
wide coverage of each LED. An FoV coefficient of $m_{\rm fov}=1.4738$ is used
as this value has been found to achieve a good match to practical PD
characteristics [23]. In the case of PR, a PD elevation angle of $\theta_{\rm
pd}=40\degree$ is used as it has been demonstrated to achieve a high data rate
in a MIMO system [23].

####  IV-D2 Configuration with BADS algorithm

It is straightforward to use the BADS algorithm introduced in Section IV-B to
search system configuration parameters of interest with the following
procedures: 1. Formulate a general black-box optimisation problem with the
considered configurable parameters as input variables $\mathbf{x}$ and average
achievable data rate $\bar{C}$ as the objective function. 2. Define bound
vectors based on the practical constraints of the considered input variables.
3. Generate an initial solution vector $\mathbf{x}_{0}$, where each variable
is randomly selected from the values between the plausible bounds with equal
probability. 4. Execute the BADS algorithm multiple times with different
initial solutions. The black-box optimisation problem for system
configurations with the SD strategy is defined as:

$\displaystyle\underset{\mathcal{X}}{\text{maximise}}$
$\displaystyle\bar{C}^{\rm SD}(\mathcal{X})\leavevmode\nobreak\
\text{with}\leavevmode\nobreak\ \mathcal{X}=\left\\{\mathtt{x}_{n_{\rm
t}},\mathtt{y}_{n_{\rm t}},\phi_{1/2},m_{\rm fov},(\theta_{\rm pd})\right\\}$
(47) subject to $\displaystyle 0<\frac{\mathtt{x}_{n_{\rm t}}}{W_{\rm
room}}<1\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ n_{\rm
t}=1,\cdots,I,$ (48) $\displaystyle 0<\frac{\mathtt{y}_{n_{\rm t}}}{L_{\rm
room}}<1\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ n_{\rm
t}=1,\cdots,I,$ (49) $\displaystyle 0\degree<\phi_{1/2}\leq 60\degree,$ (50)
$\displaystyle m_{\rm fov}\geq 1,$ (51) $\displaystyle(0\degree\leq\theta_{\rm
pd}\leq 90\degree).$ (52)

Note that $\theta_{\rm pd}$ only exist in the optimisation with PR.

Table II:  Parameters | Symbol | Values | Parameters | Symbol | Values
---|---|---|---|---|---
Room width/length/height [m] | $W_{\rm room}/L_{\rm room}/H_{\rm room}$ | $\rm 5/5/3$ | DCO-OFDM clipping level | $\kappa$ | 3.2
Total optical power [W] | $\sum_{n_{\rm t}=1}^{N_{\rm t}}p_{n_{\rm t}}$ | 80 | PD physical area [$\rm cm^{2}$] | $A_{\rm pd}$ | $1$
Height of LED/PD [m] | $\mathtt{z}_{\rm tx}/\mathtt{z}_{\rm rx}$ | 3/0.75 | Number of subcarriers | $K$ | 256
Wavelength range [nm] | $\lambda_{\rm min}/\lambda_{\rm max}$ | 400/700 | Modulation bandwidth [MHz] | $1/2T_{\rm s}$ | 50
Optical filter transmittance | $\mathcal{G}_{\rm T}$ | 1 | LED/PD 3dB bandwidth [MHz] | $f_{\rm led,c}$ | $35/106$
PD quantum efficiency | $\eta_{\rm q}$ | 0.8 | Modulation gap factor | $\Gamma$ | 6.06 dB
Effective index | $\mathfrak{n}_{\rm e}$ | 2 | Cyclic-prefix | $N_{\rm cp}$ | 30

Fig. 6(a) shows the average achievable rate results against the number of
elements (up to 16) achieved by systems with SD strategy. The remaining system
parameters (irrelevant to DoF in the spatial domain) are listed in Table II.
Regarding the gap factor value, it has been shown in [39] that a gap factor of
6 dB is sufficient for $M$-QAM modulation to achieve a bit error rate (BER) at
$1\times 10^{-3}$ for $M=2$. With an increase in $M$, the required gap factor
decreases. This implies that by using a fixed gap factor of 6 dB is sufficient
for the calculated achievable rate with a BER of $1\times 10^{-3}$. Note that
as the number of elements changes, the total transmission optical power stays
the same and is equally distributed to each LED. Both cases with a PR and HR
are presented. In addition, the cases with upward facing/random orientation
scenarios and with empirical/optimised configurations are demonstrated. With
an increase in the number of elements, the average achievable rate also
increases due to more available multiplexing channels. However, a further
increase in the number of elements (beyond 5 to 7) leads to a less average
achievable rate improvement, especially for the cases with empirical
configurations using PR. This is because with a greater number of elements,
the channel spatial correlation also increases. This leads to a channel matrix
with worse channel conditions, which achieves fewer and weaker multiplexing
channels. Intuitively, the performance of systems with optimised
configurations are better than those with empirical configurations. With 16
elements, the average achievable rates by empirically configured systems are
lower than 1100 Mbps, while the systems with optimised configurations achieved
a sum rate in the range of 1100 to 1500 Mbps. On the other hand, systems with
random orientation scenarios received a performance penalty compared to cases
with upward facing receivers due to the more severe misalignment. Regarding
different types of angular diversity receivers, the PRs offer better
performance with a fewer number of elements, while the HR offers a slightly
better performance with a large numbers of elements.

###  IV-E Wavelength division strategy

In cases employing the WD strategy, the joint multiplexing system is
equivalent to a conventional WDM system with MIMO processing techniques, which
is similar to those in [14]. All spatial domain parameters are empirically
defined as follows: all LEDs are located in the centre of the room
$\mathbf{p}_{n_{\rm t}}^{\rm led}=[W_{\rm room}/2,L_{\rm
room}/2,\mathtt{z}_{\rm tx}]$ with an orientation of $\mathbf{o}_{n_{\rm
t}}^{\rm led}=[0,0,-1]$ and all PDs have the same orientation of
$\mathbf{o}_{n_{\rm r}}^{\rm pd,\uparrow}=[0,0,1]$. The wavelength domain
configurable parameters of interest in cases with the WD strategy include the
LED central wavelength $\lambda_{n_{\rm t}}^{\rm led,c}$, optical filter
passband centre $\lambda_{n_{\rm r}}^{\rm of,c}$ and optical filter passband
width $\Delta\lambda^{\rm of}$.

(a)

(b)

Figure 6: Average achievable rate against number of elements achieved by joint
multiplexing systems with (a) SD strategy and (b) WD strategy. ’upward’ and
’random’ refer to the upward facing receiver scenario and random orientation
scenario, respectively. ‘with proc’ and ‘w/o proc’ refer to with MIMO
processing and without MIMO processing, respectively.

####  IV-E1 With and without MIMO processing

With the WD strategy, the precoding and post-detection process transform the
‘colour’-MIMO channel matrix with inter-colour crosstalk to a diagonal matrix.
This operation is expected to mitigate the excessive crosstalk caused by the
bandpass shift phenomenon. To highlight this important feature, the results of
conventional WDM systems without MIMO processing are also included for
comparison. In the case of conventional WDM systems, the precoding and post-
detection matrices are set to identity matrices:
$\mathbf{F}_{k}=\mathbf{W}_{k}=\mathbf{I}$, where $\mathbf{I}$ is defined as
an identity matrix.

####  IV-E2 Empirical configuration

In the case of empirical configurations, the passband width of optical filters
are defined by:

$\displaystyle\Delta\lambda^{\rm of}$ $\displaystyle=(\lambda_{\rm
max}-\lambda_{\rm min})/I.$ (53)

The intuition of this setting is that the passband should be wide enough to
accept most of the signal power from one LED, but not too wide to receive too
much power from LEDs of other colours. In addition, the spectrum centre of the
$i$th LED and the passband centre of the $i$th optical filter are defined by:

$\displaystyle\lambda_{i}^{\rm led,c}$ $\displaystyle=\lambda_{i}^{\rm
of,c}=\lambda_{\rm min}+(i-0.5)\Delta\lambda^{\rm of}.$ (54)

These configurations aim to uniformly distribute the spectra of LEDs and the
passband of optical filters within the visible light spectrum.

####  IV-E3 Configuration with BADS algorithm

The parameter searching using the BADS algorithm with WD strategy is similar
to that in Section IV-D2 except that the black-box optimisation problem is
redefined as follows:

$\displaystyle\underset{\mathcal{X}}{\text{maximise}}$
$\displaystyle\bar{C}^{\rm WD}(\mathcal{X})\leavevmode\nobreak\
\text{with}\leavevmode\nobreak\ \mathcal{X}=\left\\{\lambda_{n_{\rm t}}^{\rm
led,c},\lambda_{n_{\rm r}}^{\rm of,c},\Delta\lambda^{\rm of}\right\\}$ (55)
subject to $\displaystyle\lambda_{\rm min}\leq\lambda_{n_{\rm t}}^{\rm
led,c}\leq\lambda_{\rm max}\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\
n_{\rm t}=1,\cdots,I,$ (56) $\displaystyle\lambda_{\rm min}\leq\lambda_{n_{\rm
r}}^{\rm of,c}\leq\lambda_{\rm max}\leavevmode\nobreak\
\text{for}\leavevmode\nobreak\ n_{\rm r}=1,\cdots,I,$ (57)
$\displaystyle\Delta\lambda^{\rm of}>0.$ (58)

Fig. 6(b) shows the average achievable rates achieved by joint multiplexing
systems with WD strategy against the number of elements. The remaining system
parameters are the same as those listed in Table II. It can be observed that
the highest average achievable rate by the systems without MIMO processing is
less than 400 Mbps at two elements. This demonstrates the severe performance
penalty caused by the bandpass shift issue and the excessive inter-colour
interference with a large number of WDs. Therefore, the MIMO processing is
important for joint multiplexing systems with WD strategy. In contrast, the
average rates achieved by systems with MIMO processing increase with the
number of elements consistently when the number of elements is small. However,
when the number of elements is between 7 to 9, a further increase in the
number of elements no longer provides an effective boost to the average
achievable rate. In actuality, the achievable rates drop with an increase of
the number of elements in the cases of empirical configurations. In the cases
of configuration with a BADS algorithm, the performance degradation with an
optical element increase can be avoided, but no further improvement can be
observed either, as shown in Fig. 6(b). This is because the use of optical
filters with a narrow passband actively block a significant amount of optical
power from LEDs. The narrower the optical filter passband, the more power loss
occurs, which leads to a more severe degradation in the overall SNRs and
achievable rate. In addition, the increased number of elements leads to a more
severe spectrum overlap between adjacent WD channels, which causes significant
correlations between WD channels. Consequently, the multiplexing channels with
smaller eigenvalues are severely degraded. Compared to cases with SD strategy,
all LEDs and PDs have the same orientation and positions. When there is a good
alignment, all multiplexing channels are achieving high SNR and capacity. When
there is a bad alignment, none of the channels show a good performance. The
issue is more severe in the case of random orientation scenarios. This
situation is quite different from the SD systems where the chance of at least
one or a few good channels is much higher due to the better spatial diversity.
Note that the performance gap between upward scenario and corresponding random
orientation scenario is smaller in SD systems.

###  IV-F Spatial clustering with wavelength division strategy

In Section IV-D and IV-E, the limitation of SD and WD strategies has been
demonstrated. In this section, the configurations based on the approach of
‘division in either the space or wavelength domain’ are considered. The
proposed strategy is named as SCWD. The basic idea of the SCWD scenario is to
divide the LEDs and PDs into multiple groups. The optical elements in each
group are clustered with the same positions and orientations, as illustrated
in Fig. 4(a) (d). Thus, the optical elements in different groups/clusters have
different spatial features to decorrelate the channel. On the other hand, LEDs
in the same cluster have different spectra, and optical filters with different
passbands are mounted on PDs in the same cluster. This wavelength division can
remove the channel correlation between optical elements within the same
group/cluster.

Next, the details about the proposed SCWD strategy will be introduced. Assume
that the MIMO-OFDM system has $I$ LEDs and PDs divided into $L$ clusters,
where $1\leq L\leq I$. The configuration of spatial variables is similar to
that introduced in Section IV-D except that the configuration is with respect
to each spatial cluster instead of each individual optical element. On the
transmitter side, the $l$th LED cluster has a position of:

$\displaystyle\mathbf{p}_{{\rm c},l}^{\rm
led}=[\mathtt{x}_{l},\mathtt{y}_{l},\mathtt{z}_{\rm tx}],$ (59)

for $l=1,2,\cdots,L$, where $\mathtt{x}_{l}$ and $\mathtt{y}_{l}$ are the
horizontal coordinates of the $l$th LED cluster, which are configurable
parameters. All LED clusters have the same orientation of $[0,0,-1]$. In the
case of an upward facing receiver, the $l$th PD cluster has an orientation
vector of:

$\displaystyle\mathbf{o}_{{\rm c},l}^{\rm
pd,\uparrow}=\left[\cos\omega_{l}\sin\theta_{l},\sin\omega_{l}\sin\theta_{l},\cos\theta_{l}\right],$
(60)

where $\omega_{l}$ and $\theta_{l}$ are the azimuth angle and elevation angle
of the $l$th PD cluster, which are defined based on the characteristics of PR
and HR using (44), (45) and (46). Similar to the SD strategy, $\phi_{1/2}$,
$\mathfrak{m}_{\rm fov}$ and $(\theta_{\rm pd})$ are also spatial domain
configurable parameters. The number of optical elements in the $l$th cluster
is defined as:

Initialise cluster index: $l=1$

for  _$i=1,2,\cdots,I$_ do

$M_{\rm sum}=\sum_{\hat{l}=1}^{l}M_{\hat{l}}$

if  _$i >M_{\rm sum}$_ then

Cluster index increment: $l=l+1$

end if

$m=i-M_{\rm sum}+M_{l}$

Variable mapping: $\mathtt{x}_{i}=\mathtt{x}_{l}$,
$\mathtt{y}_{i}=\mathtt{y}_{l}$, $\omega_{i}=\omega_{l}$,
$\theta_{i}=\theta_{l}$, $\lambda^{\rm led}_{i}=\lambda^{\rm led}_{{\rm
c},m}$, $\lambda^{\rm of}_{i}=\lambda^{\rm of}_{{\rm c},m}$

end for

Algorithm 2 Cluster-to-element variable mapping

$\displaystyle M_{l}$ $\displaystyle=\left\\{\begin{array}[]{lr}\lceil
I/L\rceil&:l\leq\mod(I,L)\\\ \lfloor
I/L\rfloor&:\text{otherwise}\end{array}\right.,$ (63)

which ensure the sizes of different clusters are similar to a maximum
difference of one element. Note that $\lceil\cdot\rceil$ and
$\lfloor\cdot\rfloor$ are defined as the ceiling and floor operators,
respectively. This also implies that the maximum required number of WD is
$\lceil I/L\rceil$. Within the $l$th cluster, the central wavelength of the
$m$th LED is defined as $\lambda^{\rm led}_{{\rm c},m}$ on the transmitter
side for $m=1,2,\cdots,M_{l}$. On the receiver side, the passband centre of
the $m$th optical filter is defined as $\lambda^{\rm of}_{{\rm c},m}$.
Opposite to the spatial parameter configurations, the wavelength domain
parameters are identical in different clusters. Similar to the WD strategy, in
addition to $\lambda^{\rm led}_{{\rm c},m}$ and $\lambda^{\rm of}_{{\rm
c},m}$, the optical filter passband with $\Delta\lambda^{\rm of}$ is also a
configurable parameter. By unifying some of the parameters, the system
configuration process can be significantly simplified. If the configuration is
with respect to each element, the number of configurable parameters will be
$4I+4$. With the cluster-based configuration, there are $2(L+\lceil
I/L\rceil)+4$ configurable parameters, which include: $\mathtt{x}_{l}$,
$\mathtt{y}_{l}$, $\lambda_{{\rm c},m}^{\rm led}$, $\lambda_{{\rm c},m}^{\rm
of}$, $\phi_{1/2}$, $\mathfrak{m}_{\rm fov}$, $(\theta_{\rm pd})$,
$\Delta\lambda^{\rm of}$.

With cluster-based parameters $\mathtt{x}_{l}$, $\mathtt{y}_{l}$,
$\omega_{l}$, $\theta_{l}$, $\lambda^{\rm led}_{{\rm c},m}$, $\lambda^{\rm
of}_{{\rm c},m}$, the corresponding parameters with respect to the $i$th
optical element $\mathtt{x}_{i}$, $\mathtt{y}_{i}$, $\omega_{i}$,
$\theta_{i}$, $\lambda^{\rm led}_{i}$, $\lambda^{\rm of}_{i}$ can be obtained
using Algorithm 2. With a given set of spatial and wavelength domain
parameters, the average achievable rate with a specific number of clusters $L$
can be calculated. However, it is not straightforward to select $L$ that
achieves the highest average achievable rate. Therefore, the average
achievable rates with all possible numbers of clusters $\bar{C}_{L}^{\rm
SCWD}$ are evaluated for $L=1,2,\cdots,I$. Finally, we find the result
achieves the highest average achievable rate among the evaluations with
different $L$ as:

$\displaystyle\bar{C}^{\rm SCWD}=\mathrm{max}\leavevmode\nobreak\
\left\\{\bar{C}_{L=1}^{\rm SCWD},\bar{C}_{L=2}^{\rm
SCWD},\cdots,\bar{C}_{L=I}^{\rm SCWD}\right\\}.$ (64)

####  IV-F1 Empirical configuration

Regarding the empirical configuration, some spatial domain parameters are the
same as those used in the SD strategy for a fairer comparison:
$\phi_{1/2}=60\degree$, $\mathfrak{m}_{\rm fov}=1.4738$, $\theta_{\rm
pd}=40\degree$. The values of $x_{l}$ and $y_{l}$ with different number of
clusters $L$ are selected based on Fig. 5. The configuration of $\lambda^{\rm
led}_{{\rm c},m}$, $\lambda^{\rm of}_{{\rm c},m}$ and $\Delta\lambda^{\rm of}$
are similar to (53) and (54) in the WD strategy empirical configuration:

$\displaystyle\Delta\lambda^{\rm of}$ $\displaystyle=(\lambda_{\rm
max}-\lambda_{\rm min})/\lceil I/L\rceil,$ (65) $\displaystyle\lambda_{{\rm
c},m}^{\rm led}$ $\displaystyle=\lambda_{{\rm c},m}^{\rm of}=\lambda_{\rm
min}+(m-0.5)\Delta\lambda^{\rm of}.$ (66)

After evaluation of $\bar{C}_{L}^{\rm SCWD}$ for all possible $L$, the highest
achievable rate $\bar{C}^{\rm SCWD}$ is selected using (64).

####  IV-F2 Configuration with BADS algorithm

The black-box optimisation problem for SCWD strategy is defined as:

$\displaystyle\underset{\mathcal{X}}{\text{maximise}}$
$\displaystyle\bar{C}_{\rm L}^{\rm SCWD}(\mathcal{X})\leavevmode\nobreak\
\text{with}\leavevmode\nobreak\
\mathcal{X}=\left\\{\mathtt{x}_{l},\mathtt{y}_{l},\lambda_{{\rm c},m}^{\rm
led},\lambda_{{\rm c},m}^{\rm of},\phi_{1/2},\mathfrak{m}_{\rm
fov},(\theta_{\rm pd}),\Delta\lambda^{\rm of}\right\\}$ (67) subject to
$\displaystyle 0<\frac{\mathtt{x}_{l}}{W_{\rm room}}<1\leavevmode\nobreak\
\text{for}\leavevmode\nobreak\ l=1,2,\cdots,L,$ (68) $\displaystyle
0<\frac{\mathtt{y}_{l}}{L_{\rm room}}<1\leavevmode\nobreak\
\text{for}\leavevmode\nobreak\ l=1,2,\cdots,L,$ (69)
$\displaystyle\lambda_{\rm min}\leq\lambda_{{\rm c},m}^{\rm
led}\leq\lambda_{\rm max}\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\
m=1,2,\cdots,M_{l},$ (70) $\displaystyle\lambda_{\rm min}\leq\lambda_{{\rm
c},m}^{\rm of}\leq\lambda_{\rm max}\leavevmode\nobreak\
\text{for}\leavevmode\nobreak\ m=1,2,\cdots,M_{l},$ (71) $\displaystyle
0\degree<\phi_{1/2}\leq 60\degree,$ (72) $\displaystyle m_{\rm fov}\geq 1,$
(73) $\displaystyle(0\degree\leq\theta_{\rm pd}\leq 90\degree)$ (74)
$\displaystyle\Delta\lambda^{\rm of}>0,$ (75)

Note that solving this problem only finds the highest average achievable rate
with $L$ clusters. Therefore, the BADS algorithm need to be executed for $I$
times for each initial solution so that the achievable rates for all possible
$L$ are obtained. It is challenging to include the parameter $L$ as one of the
optimisation variables, as the number of parameters in the system
configuration is a function of $L$.

(a)

(b)

(c)

(d)

Figure 7: Average achievable rate of joint multiplexing systems with SCWD
strategy against number of elements. (a) Empirical configurations with upward
receiver orientation. (b) Empirical configurations with random receiver
orientation. (c) BADS-based configurations with upward receiver orientation.
(d) BADS-based configurations with random receiver orientation.

Fig. 7 shows the average rates achieved by the joint multiplexing systems with
the SCWD strategy against different numbers of elements. In addition, the
results with SD and WD strategies are also included for comparison. The
remaining parameters are listed in Table II if they are not specified. The
four plots show cases with different receiver orientation scenarios and with
different types of configurations, respectively. It can be observed that with
an increase in the number of elements, the average achievable rate by the SCWD
strategy systems also increases consistently, even when the number of elements
is greater than 10. With 16 optical elements, the data rates achieved by SCWD
systems are in the range of 1300 to 2200 Mbps, while those achieved by either
SD or WD systems fall in the range of 500 to 1400 Mbps. The achievable rate
improvement is between 36% and 74% compared to the SD strategy, and the
improvement is between 47% and 135% compared to the WD strategy. In general,
the SCWD systems with PR shows a slightly better performance at a high number
of elements compared to those with HR.

###  IV-G Insights into the system parameter configurations and SCWD strategy

In this subsection, we will discuss the difference between the parameters used
in the empirical configurations and configurations with the BADS algorithm. We
will also explain why the SCWD strategy is able to provide extra multiplexing
gains when compared to SD and WD strategies. For simplicity, the cases with a
PR and the upward receiver orientation are demonstrated as an example.

(a)

(b)

Figure 8: (a) Comparison between empirical and BADS-based configuration
parameters used in 16-elements joint multiplexing systems with SCWD strategy.
The demonstrated systems use PR and upward receiver orientation scenario. 1)
LED cluster position layout. 2) Angular parameters. 3) Wavelength-dependent
parameters. (b) Achievable rate of systems using optimised configurations with
PR and upward receiver orientation scenario against number of elements in two
special conditions. Condition 1: larger room of size $\rm
10\leavevmode\nobreak\ m\times 10\leavevmode\nobreak\ m$. Condition 2:
extremely small space of size $\rm 0.05\leavevmode\nobreak\ m\times
0.05\leavevmode\nobreak\ m$. Default condition: room size of $\rm
5\leavevmode\nobreak\ m\times 5\leavevmode\nobreak\ m$.

The comparison between parameters used in empirical and BADS-based
configurations with 16 elements is shown in Fig. 8(a). In both configurations,
there are four spatial clusters, as demonstrated in Fig. 8(a) 1). Regarding
the geometric parameters, the spatial positions of LED clusters are closer to
the room centre in the case of BADS-based configuration compared to the
empirical configuration. In addition, a smaller half-power semiangle of
$\phi_{1/2}=35.2\degree$ is used on the transmitter side and a smaller FoV
coefficient of $m_{\rm fov}=1$ is used on the receiver side compared to the
empirical configuration, as shown in Fig. 8(a) 2). These parameter changes
lead to weaker spatial channel correlation and thereby achieving higher
overall multiplexing gain. Regarding the wavelength domain parameters, four
WDs are used within each spatial cluster, as shown in Fig. 8(a) 3). Compared
to the case of empirical configuration, all LED spectra shift to the longer
wavelengths with smaller separations in the case of BADS-based configuration.
This is because the PD has a higher responsivity at a longer wavelength and
the additional inter-colour interference introduced by smaller separation
causes negligible distortion. In addition, the passband of optical filters are
extended to longer wavelengths in the case of BADS-based configuration, which
compensate bandpass shift issue when light incident angle is large.

In order to demonstrate the characteristics of the SCWD strategy, we can
evaluate the performance of joint multiplexing systems using BADS-based
configurations with PR and upward receiver orientation scenario under several
conditions. In addition to cases with normal conditions (i.e. a room size of
$\rm 5\leavevmode\nobreak\ m\times 5\leavevmode\nobreak\ m\times
3\leavevmode\nobreak\ m$) shown in Fig. 7(c), we further evaluate two special
conditions. In the first condition, we show the transmissions in a large room
of size $\rm 10\leavevmode\nobreak\ m\times 10\leavevmode\nobreak\ m\times
3\leavevmode\nobreak\ m$. The remaining system parameters are identical to
those used in Sections IV-D, IV-E and IV-F. Due to the increase in coverage
area, the average link distance becomes greater and the misalignment issue due
to random user position and orientation is more severe. These factors lead to
a further decrease in the data rate achieved by joint multiplexing systems, as
shown in Fig. 8(b). The achievable rate decrease is more severe for the case
of WD strategy, where random users are more likely to receive a weaker signal
as all LEDs are located in the room centre. Therefore, the achievable rate of
the WD strategy is considerably lower than the case of SD strategy with the
same number of elements. Due to the inferior performance of system with WD
strategy, the SCWD strategy provides a decreased achievable rate improvement
against its counterparts compared to the case with $\rm 5\leavevmode\nobreak\
m\times 5\leavevmode\nobreak\ m\times 3\leavevmode\nobreak\ m$ room. In the
second special condition, we show that transmissions in an extremely small
indoor space of size $\rm 0.05\leavevmode\nobreak\ m\times
0.05\leavevmode\nobreak\ m\times 3\leavevmode\nobreak\ m$. With such a small
coverage area and the upward receiver scenario, all users experience a
minimised link distance with good alignment. Consequently, the passband shift
issue in systems with the WD strategy can be avoided. Note that a half-power
semiangle of $\phi_{1/2}=1.5\degree$ is also used in the case of WD strategy
to maximise the detected optical power. In the case of SD strategy,
performance degradation due to greater link distance and misalignment is
avoided. In addition, a BADS-based configuration solution with a very
collimated light beam and light reception pattern is obtained in case of the
SD strategy. Due to all the above factors, the joint multiplexing systems with
WD and SD strategies can operate in ideal conditions. Specifically, the
multiplexing gain scales almost linearly with the number of elements, as shown
in Fig. 8(b). In this condition, despite the minor performance gap between WD
and SD systems, the multiplexing gain of both systems are approaching the
theoretical limit of multiplexing systems. Consequently, the system with SCWD
strategy exhibits a similar ideal performance.

Despite the undesired performance exhibited in the two special conditions,
using a MIMO VLC system to cover an extremely large or extremely small area is
unlikely in practice. As long as the following two conditions are fulfilled,
SCWD strategy will provide a considerable improvement: 1) there is a minor
performance gap between joint multiplexing systems with SD and WD strategies,
2) multiplexing transmissions in individual spatial and wavelength domains are
inefficient. In general conditions of where a medium sized room is covered,
similar achievable rates are likely to be achieved by systems with SD and WD
strategies. In addition, the efficiency of VLC multiplexing systems is limited
by the random user position/orientation and channel correlation. Therefore,
the SCWD strategy is likely to provide a remarkable achievable rate
improvement, as shown in Fig. 7. The configuration problem in joint
multiplexing systems is similar to a power allocation problem in a multi-
carrier transmission system with two channels. If both channels have similar
qualities, loading them with power leads to a higher performance compared to
loading all the power to one channel. However, if one of the channels is
unreliable or the spectral efficiency scales linearly with the power increase,
loading two channels will no longer be beneficial.

###  IV-H Practical implementation and illumination

For implementations of the VLC MIMO joint multiplexing system with the
proposed SCWD strategy, there are several practical solutions. For example, it
is challenging to mount more than ten PDs on the same side of a portable
device, while it is easier to mount a cluster of a few PDs on each side of a
device with a different orientation [27]. For small devices, such as smart
phones, mounting one or two PDs on each side of the device is realistic. The
system with a greater number of PDs can be implemented on larger devices such
as tablets, laptop or Internet-of-things (IoT) device. On the transmitter
side, the SCWD strategy gather multiple LEDs with different wavelengths in the
same cluster. This feature provides an opportunity to generate white light
illumination by manipulating the intensity ratio of LEDs with different
colours [40]. In the cases where the LEDs in the cluster are insufficient to
provide the desired white light, additional LEDs with suitable wavelengths can
be installed in the LED cluster so that the combined light reaches the
illumination requirement.

##  V Conclusions

A framework for a VLC MIMO-OFDM joint multiplexing system with a channel
dependent on spatial, wavelength and frequency domain characteristics is
proposed. This study focused on investigating the suitable system
configurations for the VLC joint multiplexing system with different strategies
in order to achieve a higher multiplexing gain, thereby achieving a higher
aggregate transmission data rate. The proposed SCWD strategy shows a
considerable improvement in the average achievable rate compared to cases with
SD and WD strategies. The use of the BADS black-box optimisation algorithm
improves the achievable rates further. The simulation results demonstrate an
average achievable rate improvement of 36% to 135% with 16 LEDs/PDs by using
the proposed SCWD strategy compared to those achieved with SD and WD
strategies. This work also demonstrated the possibility of using a VLC MIMO-
OFDM system with limited modulation bandwidth to achieve multi-Gbps
transmission data rates, which shows the potential of VLC/LiFi systems in next
generation wireless networks in terms of achievable rate capability.  In
future studies, the scenarios with multiple access or multi-user MIMO will be
investigated.

## References

  * [1] W. Jiang, B. Han, M. A. Habibi, and H. D. Schotten, “The Road Towards 6G: A Comprehensive Survey,” _IEEE Open Journal of the Communications Society_ , vol. 2, pp. 334–366, 2021.
  * [2] H. Haas, L. Yin, Y. Wang, and C. Chen, “What is LiFi?” _J. Lightw. Technol._ , vol. 34, no. 6, pp. 1533–1544, Mar. 2016.
  * [3] J. Raring, C. Lee, M. S. Islim, A. Sparks, S. Videv, M. McLaurin, B. Shah, P. Rudy, and H. Haas, “>25 Gbit/s LiFi with Laser Based SMD White Light Source,” in _Optical Fiber Communications Conference and Exhibition_ , 2021, pp. 1–3.
  * [4] E. Xie, R. Bian, X. He, M. S. Islim, C. Chen, J. J. D. McKendry, E. Gu, H. Haas, and M. D. Dawson, “Over 10 gbps vlc for long-distance applications using a gan-based series-biased micro-led array,” _IEEE Photonics Technology Letters_ , vol. 32, no. 9, pp. 499–502, 2020.
  * [5] T. Fath and H. Haas, “Performance comparison of mimo techniques for optical wireless communications in indoor environments,” _IEEE Trans. Commun._ , vol. 61, no. 2, pp. 733–742, February 2013.
  * [6] K. Ying, H. Qian, R. J. Baxley, and S. Yao, “Joint optimization of precoder and equalizer in mimo vlc systems,” _IEEE J. Sel. Areas Commun._ , vol. 33, no. 9, pp. 1949–1958, Sep. 2015.
  * [7] P. M. Butala, H. Elgala, and T. D. C. Little, “SVD-VLC: A Novel Capacity Maximizing VLC MIMO System Architecture Under Illumination Constraints,” in _2013 IEEE Globecom Workshops (GC Wkshps)_ , 2013, pp. 1087–1092.
  * [8] Y. Hong, T. Wu, and L.-K. Chen, “On the Performance of Adaptive MIMO-OFDM Indoor Visible Light Communications,” _IEEE Photonics Technology Letters_ , vol. 28, no. 8, pp. 907–910, 2016.
  * [9] C. He, T. Q. Wang, and J. Armstrong, “Performance Comparison between Spatial Multiplexing and Spatial Modulation in Indoor MIMO Visible Light Communication Systems,” in _IEEE International Conf. on Commun._ , 2016, pp. 1–6.
  * [10] Y. Tan and H. Haas, “Coherent LiFi System With Spatial Multiplexing,” _IEEE Transactions on Communications_ , vol. 69, no. 7, pp. 4632–4643, 2021\.
  * [11] R. Bian, I. Tavakkolnia, and H. Haas, “15.73 Gb/s Visible Light Communication With Off-the-Shelf LEDs,” _Journal of Lightwave Technology_ , vol. 37, no. 10, pp. 2418–2424, 2019.
  * [12] G. Cossu, A. Khalid, P. Choudhury, R. Corsini, and E. Ciaramella, “3.4 Gbit/s Visible Optical Wireless Transmission based on RGB LED,” _Optics Express_ , vol. 20, no. 26, pp. B501–B506, 2012.
  * [13] L. Cui, Y. Tang, H. Jia, J. Luo, and B. Gnade, “Analysis of the Multichannel WDM-VLC Communication System,” _Journal of Lightwave Technology_ , vol. 34, no. 24, pp. 5627–5634, 2016.
  * [14] N. Omura, A. Higashi, J. Yabuuchi, T. Iwamatsu, and S. Oshiba, “Experimental Demonstration of OFDM Based WDM-MIMO Visible Light Communication System,” in _2018 Asia-Pacific Microwave Conference (APMC)_ , 2018, pp. 872–874.
  * [15] H. Lee, I. Lee, and S. H. Lee, “Deep Learning based Transceiver Design for Multi-Colored VLC Systems,” _Optics Express_ , vol. 26, no. 5, pp. 6222–6238, 2018.
  * [16] A. Burton, P. Chvojka, P. A. Haigh, Z. Ghassemlooy, and S. Zvanovec, “Optical Filter-Less WDM for Visible Light Communications Using Defocused MIMO,” _Electronics_ , vol. 10, no. 9, p. 1065, 2021.
  * [17] R. Wang, Q. Gao, J. You, E. Liu, P. Wang, Z. Xu, and Y. Hua, “Linear Transceiver Designs for MIMO Indoor Visible Light Communications Under Lighting Constraints,” _IEEE Trans. Commun._ , vol. 65, no. 6, pp. 2494–2508, 2017.
  * [18] Y. Xiao and Y.-J. Zhu, “Chromaticity-Adaptive Generalized Spatial Modulation for MIMO VLC With Multi-Color LEDs,” _IEEE Photonics Journal_ , vol. 11, no. 4, pp. 1–12, 2019.
  * [19] R. Sharma, A. C. Kumari, M. Aggarwal, and S. Ahuja, “Optimal LED Deployment for Mobile Indoor Visible Light Communication System: Performance Analysis,” _AEU-International Journal of Electronics and Communications_ , vol. 83, pp. 427–432, 2018.
  * [20] M. T. Niaz, F. Imdad, S. Kim, and H. S. Kim, “Deployment Methods of Visible Light Communication Lights for Energy Efficient Buildings,” _Optical Engineering_ , vol. 55, no. 10, p. 106113, 2016.
  * [21] M. S. Gismalla, M. F. L. Abdullah, M. S. Ahmed, W. A. Mabrouk, A.-F. Najib, E. Saeid, A. Supa’at, and B. Das, “Design and Analysis of Different Optical Attocells Deployment Models for Indoor Visible Light Communication System,” _International Journal of Integrated Engineering_ , vol. 13, no. 6, pp. 253–264, 2021.
  * [22] J. Xu, C. Gong, J. Luo, and Z. Xu, “LED Half-Power Angle Optimization for Ultra-Dense Indoor Visible Light Communication Network Deployment,” _IEEE Open Journal of the Communications Society_ , vol. 1, pp. 835–848, 2020\.
  * [23] A. Nuwanpriya, S. W. Ho, and C. S. Chen, “Indoor MIMO Visible Light Communications: Novel Angle Diversity Receivers for Mobile Users,” _IEEE J. Sel. Areas Commun._ , vol. 33, no. 9, pp. 1780–1792, Sep. 2015.
  * [24] K.-H. Park, W. G. Alheadary, and M.-S. Alouini, “A Novel Mirror Diversity Receiver for Indoor MIMO Visible Light Communication Systems,” in _2016 IEEE 27th Annual International Symposium on Personal, Indoor, and Mobile Radio Communications (PIMRC)_ , 2016, pp. 1–6.
  * [25] A. A. Purwita, A. Yesilkaya, I. Tavakkolnia, M. Safari, and H. Haas, “Effects of Irregular Photodiode Configurations for Indoor MIMO VLC with Mobile Users,” in _2019 IEEE 30th Annual International Symposium on Personal, Indoor and Mobile Radio Communications (PIMRC)_ , 2019, pp. 1–7.
  * [26] P. Ge, X. Liang, J. Wang, C. Zhao, X. Gao, and Z. Ding, “Optical Filter Designs for Multi-Color Visible Light Communication,” _IEEE Transactions on Communications_ , vol. 67, no. 3, pp. 2173–2187, 2019.
  * [27] C. Chen, I. Tavakkolnia, M. D. Soltani, M. Safari, and H. Haas, “Hybrid Multiplexing in OFDM-based VLC Systems,” in _2020 IEEE Wireless Communications and Networking Conference (WCNC)_ , 2020, pp. 1–6.
  * [28] D. Tsonev, S. Sinanovic, and H. Haas, “Complete Modeling of Nonlinear Distortion in OFDM-based Optical Wireless Communication,” _J. Lightw. Technol._ , vol. 31, no. 18, pp. 3064–3076, Sep. 15 2013.
  * [29] S. Dimitrov, S. Sinanovic, and H. Haas, “Clipping Noise in OFDM-Based Optical Wireless Communication Systems,” _IEEE Trans. Commun._ , vol. 60, no. 4, pp. 1072–1081, Apr. 2012.
  * [30] H. Schulze, “Frequency-Domain Simulation of the Indoor Wireless Optical Communication Channel,” _IEEE Trans. Commun._ , vol. 64, no. 6, pp. 2551–2562, June 2016.
  * [31] J. Barry, J. Kahn, W. Krause, E. Lee, and D. Messerschmitt, “Simulation of Multipath Impulse Response for Indoor Wireless Optical Channels,” _IEEE J. Sel. Areas Commun._ , vol. 11, no. 3, pp. 367–379, Apr. 1993.
  * [32] S. Hranilovic, _Wireless Optical Communication Systems_. Springer Science & Business Media, 2006.
  * [33] T. E. B. Cunha, J.-P. M. G. Linnartz, and X. Deng, “Achievable Rate of LED-based Distributed MIMO OWC Systems Under a Per-LED Power Constraint,” in _2021 17th International Symposium on Wireless Communication Systems (ISWCS)_ , 2021, pp. 1–6.
  * [34] M. D. Soltani, M. A. Arfaoui, I. Tavakkolnia, A. Ghrayeb, M. Safari, C. M. Assi, M. O. Hasna, and H. Haas, “Bidirectional Optical Spatial Modulation for Mobile Users: Toward a Practical Design for LiFi Systems,” _IEEE Journal on Selected Areas in Communications_ , vol. 37, no. 9, pp. 2069–2086, 2019\.
  * [35] C. Audet and W. Hare, _Derivative-free and Blackbox Optimization_. Springer, 2017, vol. 2.
  * [36] L. Acerbi and W. J. Ma, “Practical Bayesian Optimization for Model Fitting with Bayesian Adaptive Direct Search,” in _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , 2017, pp. 1834–1844.
  * [37] J.-A. Li, D. Dong, Z. Wei, Y. Liu, Y. Pan, F. Nori, and X. Zhang, “Quantum Reinforcement Learning During Human Decision-Making,” _Nature human behaviour_ , vol. 4, no. 3, pp. 294–307, 2020.
  * [38] C. Audet and J. E. Dennis Jr, “Mesh Adaptive Direct Search Algorithms for Constrained Optimization,” _SIAM Journal on optimization_ , vol. 17, no. 1, pp. 188–217, 2006.
  * [39] S. Mardanikorani, X. Deng, and J.-P. M. G. Linnartz, “Sub-Carrier Loading Strategies for DCO-OFDM LED Communication,” _IEEE Transactions on Communications_ , vol. 68, no. 2, pp. 1101–1117, 2020.
  * [40] C. Gong, S. Li, Q. Gao, and Z. Xu, “Power and Rate Optimization for Visible Light Communication System With Lighting Constraints,” _IEEE Transactions on Signal Processing_ , vol. 63, no. 16, pp. 4245–4256, 2015.

  *[VLC]: visible light communication
  *[LiFi]: light fidelity
  *[LED]: light-emitting diode
  *[SMX]: spatial multiplexing
  *[WDM]: wavelength division division
  *[MIMO]: multiple-input multiple-output
  *[OFDM]: orthogonal frequency division multiplexing
  *[PD]: photodiode
  *[SCWD]: spatial clustering with wavelength division
  *[BADS]: Bayesian adaptive direct search
  *[DoF]: degree-of-freedom
  *[DCO]: DC-biased optical
  *[FFT]: fast Fourier transform
  *[QAM]: quadrature amplitude modulation
  *[IFFT]: inverse fast Fourier transform
  *[ISI]: inter-symbol interference
  *[CP]: cyclic prefix
  *[PAPR]: peak-to-average power ratio
  *[SNR]: signal-to-noise ratio
  *[SVD]: singular-value decomposition
  *[DC]: direct current
  *[DTFT]: discrete-time Fourier transform
  *[RRC]: root-square raised cosine
  *[LoS]: line-of-sight
  *[NLoS]: non-line-of-sight
  *[FoV]: field-of-view
  *[UE]: user equipment
  *[3D]: three-dimensional
  *[2D]: two-dimensional
  *[MADS]: mesh adaptive direct search
  *[SD]: space division
  *[WD]: wavelength division
  *[PR]: pyramid receiver
  *[HR]: hemispheric receiver
  *[BER]: bit error rate
  *[IoT]: Internet-of-things