# Synchronization Scheme based on Pilot Sharing in Cell-Free Massive MIMO
Systems

Qihao Peng, Hong Ren, Zhendong Peng, Cunhua Pan, ,
Maged Elkashlan, , Dongming Wang, , Jiangzhou Wang, , and Xiaohu You Q. Peng
and M. Elkashlan are with the School of Electronic Engineering and Computer
Science at Queen Mary University of London, U.K. (e-mail: {q.peng,
maged.elkashlan}@qmul.ac.uk). Zhendong Peng is with the Department of
Electrical and Computer Engineering, The University of British Columbia,
Vancouver, BC V6T 1Z4, Canada (e-mail: zhendongpeng@ece.ubc.ca). Jiangzhou
Wang is with the School of Engineering, University of Kent, CT2 7NT
Canterbury, U.K. (e-mail: j.z.wang@kent.ac.uk). H. Ren, C. Pan, Dongming Wang,
and Xiaohu You are with National Mobile Communications Research Laboratory,
Southeast University, Nanjing, China.
(e-mail:{hren,cpan,wangdm,xhyu}@seu.edu.cn). (Corresponding author: Hong Ren
and Cunhua Pan)

###### Abstract

This paper analyzes the impact of pilot-sharing scheme on synchronization
performance in a scenario where several slave access points (APs) with
uncertain carrier frequency offsets (CFOs) and timing offsets (TOs) share a
common pilot sequence. First, the Cramer-Rao lower bound (CRB) with pilot
contamination is derived for pilot-pairing estimation. Furthermore, a maximum
likelihood algorithm is presented to estimate the CFO and TO among the pairing
APs. Then, to minimize the sum of CRBs, we devise a synchronization strategy
based on a pilot-sharing scheme by jointly optimizing the cluster
classification, synchronization overhead, and pilot-sharing scheme, while
simultaneously considering the overhead and each AP’s synchronization
requirements. To solve this NP-hard problem, we simplify it into two sub-
problems, namely cluster classification problem and the pilot sharing problem.
To strike a balance between synchronization performance and overhead, we first
classify the clusters by using the K-means algorithm, and propose a criteria
to find a good set of master APs. Then, the pilot-sharing scheme is obtained
by using the swap-matching operations. Simulation results validate the
accuracy of our derivations and demonstrate the effectiveness of the proposed
scheme over the benchmark schemes.

###### Index Terms:

Cell-free massive MIMO, synchronization, pilot assignment, graph theory,
coordinated integrated sensing and communications.

##  I Introduction

Recently, cell-free massive multiple-input multiple-output (mMIMO) systems
have emerged as a promising technology for the sixth generation (6G) systems
to simultaneously meet the ever-growing demand on data rate, latency, and
decoding error probability [1, 2, 3, 4]. By geographically deploying a
multitude of distributed access points (APs) and eliminating the cell
boundary, the distributed APs can seamlessly cooperate to mitigate inter-cell
interference [5, 6, 7], and then jointly perform the coordinated multi-point
transmission/reception (CoMP) by exploiting the multiple enhanced path gains,
which consequently improves the spectral efficiency and user experience.
However, as each AP is based on its own local oscillator (LO), the carrier
frequency offset (CFO) and timing offset (TO) among distributed APs leads to
asynchronous reception, which severely deteriorates the cooperative
transmission and reception performance [8, 9]. Therefore, to reap the
promising gains provided by CoMP, it is crucial to address the fundamental
synchronization challenges in distributed architectures.

One method is to connect all distributed APs by physical wires, i.e., coaxial
cables or optical fibers. Though this method seems logically straightforward,
it has economic and technical challenges for practical implementation.
Specifically, the deployment of extensive wired infrastructure incurs
substantial costs in terms of installation and maintenance, especially in
scenarios where APs are distributed over large geographical areas or in
challenging terrains. Furthermore, the utilization of wired connections not
only constrains the system’s scalability and flexibility but also renders it
highly vulnerable to environmental damage. Another solution is to deploy a
global positioning system (GPS) disciplined oscillator at each AP as an
accurate time source for individual devices [10], which is also typically
expensive and impractical.

To address the above issues, extensive contributions have been devoted to the
realization of wireless synchronization [11, 12, 13, 14, 15, 16, 17, 18, 19].
Generally, the widely used methods can be divided into two categories, namely
reference signal-based synchronization [11, 12, 13, 15, 14, 16] and pilot-
based synchronization [17, 19]. Specifically, a master AP transmits an out-of-
band reference signal, while slave APs can synchronize the frequency and phase
by tracking the reference signal [11]. This approach still causes various
errors due to different LOs. To mitigate the effect of LO, the authors of [12]
proposed an AirShare scheme where the distributed receivers adjust LO to
generate a frequency equal to the difference of two tones. Then, this strategy
was further investigated in [13, 14] and also adopted in distributed microwave
power transfer systems [15]. However, these synchronization schemes based on
reference signals utilize extra bandwidth and circuits, which causes lower
spectral efficiency and higher power consumption. To tackle this issue, the
authors of [17] proposed a pilot-pair synchronization strategy where the slave
APs estimate the CFO and TO by receiving the burst pilot from the master APs.
Then, the design of optimal pairs was studied in [18]. Inspired by the pilot-
pair scheme, a similar strategy was applied to the distributed RadioWeaves
infrastructure [19]. However, these synchronization schemes based on pairing
pilot [17, 18, 19] did not analyze the impact of pilot allocation on the
estimation performance.

Currently, various pilot allocation schemes have been adopted in cell-free
mMIMO [20, 21, 22, 23, 24, 25]. Particularly, it was demonstrated that pilot
contamination is related to the distance among users that reuse the pilots
[20]. Then, to minimize pilot contamination, the system performance can be
maximized by using the Tabu-search-based algorithm [21] and the Hungarian
algorithm [22]. Furthermore, to further reduce the pilot overhead, the pilot
allocation relying on graphic theory was proposed in [23]. To meet the diverse
requirements of users, pilot sharing schemes based on graphic theory were
investigated in [24, 25]. Although pilot allocation schemes have been
extensively studied in cell-free mMIMO systems and a similar pilot sharing
scheme can be applied for pilot allocation, the relationship between
synchronization performance and pilot allocation is still unknown.
Furthermore, it is still challenging to design the optimal synchronization
scheme with limited overhead.

To fill this gap, we first derive the Cramer-Rao bound (CRB) for the CFO and
TO synchronization and then reveal the relationship between the estimation
error and synchronous pilot allocation. Finally, a master-slave AP
synchronization scheme based on a novel pilot-sharing strategy is proposed to
enhance the synchronization performance with limited overhead. Our
contributions are summarized as follows:

  1. 1.

By considering the impact of pilot-sharing schemes on estimation errors, we
analyze a more practical scenario where several slave APs with uncertain CFOs
and TOs share a common pilot sequence. Then, a new CRB for pilot-pair
estimation is derived while considering the pilot contamination. Finally, a
maximum likelihood (ML) algorithm is provided to estimate the CFO and TO among
two pairing APs.

  2. 2.

As for the pilot-pairing synchronization scheme, our objective is to minimize
the sum of CRBs by jointly optimizing clustering and the pilot allocation
scheme while simultaneously considering all APs’ synchronization requirements
and limited overhead. To solve this NP-hard problem, we simplify it into two
sub-problems, including the cluster classification and the pilot sharing
scheme. For the first sub-problem, the requirements on CFO and TO are
transformed into the desired signal-to-interference -and-noise ratio (SINR),
based on which we derive the maximum distance in the cluster. Then, the
clusters are classified by iteratively performing the K-means algorithm, and a
good set of master APs is found by using the proposed criteria. To solve the
second sub-problem, a pilot sharing scheme based on the swap-matching
operations is proposed to maximize the synchronization performance with
limited overhead.

  3. 3.

Our simulation results validate the accuracy of our derivations but also
demonstrate that increasing the synchronous pilot length can decrease the
estimation error. Furthermore, classifying clusters can shorten the distance
between master APs and slave APs, which can enhance the path gain and
significantly improve the estimation performance. Finally, our simulation
results also confirm its superiority to other pilot-pair algorithms.

The remainder of this paper is organized as follows. In Section II, the system
model is provided, and the CRB and ML estimations are given, respectively. By
optimizing the cluster classification, pilot length and pilot sharing scheme,
the sum of CRBs is minimized in Section III. Our simulation results are
presented in Section IV. Finally, our conclusions are drawn in Section V.

##  II System Model

Figure 1: Illustration of pilot-pair synchronization in cell-free mMIMO
systems.

Since each AP has a unique oscillator and the synchronization is non-ideal,
the very small residual CFO and TO will cause slow relative drift among
different APs, which has a severe impact on cooperative transmissions. To
prevent these errors from accumulating over time, the APs need to exchange the
information and perform synchronization frequently. However, owing to the
excessive number of APs, it is impractical for all APs to exchange the
synchronization information with the central processing unit (CPU) directly by
wireless channels, which leads to high overhead and low efficiency.

As illustrated in Fig. 1, to solve the scalability of cell-free mMIMO, $K$ APs
can be divided into $K_{C}$ non-overlapping clusters, denoted as
${\mathcal{C}}_{k}$, $\forall k\in\\{1,2,\cdots,K_{C}\\}$, and ${\cal
C}_{k}\cap{\cal C}_{k^{\prime}}=\emptyset$, if $k\neq k^{\prime}$.
Furthermore, we assume that each cluster has only one mater AP and several
slave APs. We denote the master AP in the $k$-th cluster as $j_{k}$ and the
set of slave APs in the $k$-th cluster as ${\mathcal{S}}_{k}$ with ${\cal
S}_{k}\cup j_{k}={\cal C}_{k}$ and ${\cal S}_{k}\cap j_{k}=\emptyset$. To
synchronize all APs with low overhead, slave APs of ${\cal S}_{k}$ first
transmit the burst pilot to master APs, while the master APs estimate the CFO
and TO based on the received signals and then exchange information with the
CPU. Then, the CPU performs the synchronization operation by transmitting the
commands to the master APs.

We consider MIMO-OFDM APs, wherein the antennas at each AP operate according
to a common LO. Owing to the slow relative drift among different APs and
imperfect synchronization, we denote the CFO and TO between the $i$-th AP and
the $j_{k}$-th AP as $\Delta f_{i,j_{k}}$ and $\Delta\tau_{i,j_{k}}$,
respectively. Furthermore, since the antennas at each AP are based on a common
local oscillator, the CFO and TO between the transmit antennas of the $i$-th
AP and the $j_{k}$-th AP’s receive antennas are the same. Therefore, for ease
of derivation, we apply the received signal based on the single-antenna model
for analysis.

We consider joint TO and CFO pilot-aided estimation for a transmitter-receiver
pair, where the channel is multi-path time-invariant with additive white
Gaussian noise (AWGN), with unknown path coefficients and multiple path
delays. Generally, the number of the channel taps $L_{i,j_{k}}$ between the
$i$-th slave AP and the $j_{k}$-th mater AP does not exceed the length of the
cyclic prefix (CP). The multi-path channel impulse response from the $i$-th
slave AP to the $j_{k}$-th master AP can be represented by

${h_{i,j_{k}}}\left(t\right)=\sum\limits_{l=0}^{{L_{i,j_{k}}}-1}{\alpha_{i,j_{k}}^{l}\sqrt{{\beta_{i,j_{k}}}}\psi^{l}_{i,j_{k}}\delta\left({t-\varsigma_{i,j_{k}}^{l}}\right)},$
(1)

where $\alpha_{i,j_{k}}^{l}$ and $\varsigma_{i,j_{k}}^{l}$, $\forall
l\in\\{0,1,2,\cdots,L_{i,j_{k}}-1\\}$, are the path gain coefficient and path
delay of the $l$-th channel tap, respectively. $\beta_{i,j_{k}}$ is the
distance-dependent path loss, $\psi^{l}_{i,j_{k}}\in\mathcal{CN}(0,1)$ denotes
the small-scale fading factor, and $\delta(t)$ is the impulse response.

In the phase of network setup, the slave APs transmit the pilot sequence for
synchronization based on a pilot burst transmission. Then, the $j_{k}$-th
master AP receives time-domain baseband signal from several slave APs that
share the common pilot sequence, which is given by

$\begin{split}{y_{j_{k}}}\left(t\right)&=\sum\limits_{i\in{{\cal
P}_{p}}}h_{i,j_{k}}(t)\otimes
s_{p}(t){e^{j2\pi\Delta{f_{i,j_{k}}}t}}+{w_{i,j_{k}}}\left(t\right)\\\
&=\sum\limits_{i\in{{\cal
P}_{p}}}{\sum\limits_{l=0}^{{L}_{i,j_{k}}-1}{h_{i,j_{k}}^{l}{s_{p}}\left({t-{\Delta\tau_{i,j_{k}}}-\varsigma_{i,j_{k}}^{l}}\right)}{e^{j2\pi\Delta{f_{i,j_{k}}}t}}}+{w_{i,j_{k}}}\left(t\right),\end{split}$
(2)

where ${\cal P}_{p}$ is the set of slave APs that share the $p$-th pilot
sequence,
$h^{l}_{i,j_{k}}=\alpha_{i,j_{k}}^{l}\sqrt{{\beta_{i,j_{k}}}}\psi^{l}_{i,j_{k}}$
is the $l$-th path coefficient, $\otimes$ is the convolution operator, and
$w_{i,j_{k}}\in\mathcal{CN}(0,\sigma^{2})$ represents the additive noise.
$s_{p}(t)=\frac{P_{t}}{{\sqrt{{T_{s}}}}}\sum\limits_{n=0}^{{N}-1}{s_{p}^{n}\prod\left({\frac{t}{{{T_{s}}}}-n}\right)}$
denotes the $p$-th pilot sequence, $s_{p}^{n}$ is the time-domain chip symbol
with ${\bf s}_{p}=\\{s_{p}^{0},s_{p}^{1},...,s_{p}^{N-1}\\}$, $T_{s}$ is the
chip interval, $P_{t}$ is the transmission power, and $N$ is the length of
pilot sequence. $\prod(t)$ denotes the rectangular pulse, which can be
expressed as

$\prod\left(t\right)=\left\\{{\begin{array}[]{*{20}{c}}{1,}&{0\leq t\leq 1}\\\
{0,}&{{\rm{otherwise}}}.\end{array}}\right.$ (3)

Then, the receiver performs chip-matched filtering, and sampling per one chip
by substituting $t=mT_{s}$, $m\in\\{0,\cdots,M-1\\}$. For ease of expression,
we denote $\nabla(t)=\prod(t)\otimes\prod(t)$ as the convolution of two
rectangular pulses. The received discrete-time signal is given by [17]

$\begin{split}{y_{j_{k}}}\left[m\right]=\sum\limits_{i\in{{\cal
P}_{p}}}{\left({\sum\limits_{l=0}^{{L_{i,j_{k}}}-1}{h_{i,j_{k}}^{l}\sum\limits_{n=0}^{N-1}{s_{p}^{n}\nabla\left({m-n-\frac{{\Delta{\tau_{i,j_{k}}}+\varsigma_{i,j_{k}}^{l}}}{{{T_{s}}}}}\right)}}}\right){e^{j2\pi\Delta{f_{i,j_{k}}}{T_{s}}m}}}+{w_{i,j_{k}}}\left[m\right].\end{split}$
(4)

By using a coarse frame synchronization protocol, the receiver has a coarse
knowledge of the frame timing such that it expects to find the pilot burst in
a given interval including the guard time. The master AP collects $M$ samples
such that the sampled interval of duration $MT_{s}$ contains the pilot
sequence.

Next, we collect the $M$ received samples into an $M\times 1$ vector. The
$j_{k}$-th master AP’s received signal is defined as
${\bf{y}}_{j_{k}}={\left[{y_{j_{k}}\left[0\right],y_{j_{k}}\left[1\right],\ldots,y_{j_{k}}\left[{{M}-1}\right]}\right]^{T}}$.
Then, the CFO matrix is defined as

${\bf{F}}\left({{\Delta
f_{i,{j_{k}}}}}\right)={\rm{diag}}\left({1,{e^{j2\pi{\Delta
f_{i,{j_{k}}}}{T_{s}}}},\ldots,{e^{j2\pi\left({M-1}\right){\Delta
f_{i,{j_{k}}}}{T_{s}}}}}\right).$ (5)

The $M\times\left({M-N+1}\right)$ matrix of the $p$-th pilot sequence is given
by

${{\bf{S}}_{p}}=\left[{\begin{array}[]{*{20}{c}}{s_{p}^{0}}&0&\ldots&0\\\
{s_{p}^{1}}&{s_{p}^{0}}&{}\hfil&\vdots\\\
\vdots&{s_{p}^{1}}&{}\hfil&{s_{p}^{0}}\\\
{s_{p}^{{N}-1}}&\vdots&\ddots&{s_{p}^{1}}\\\ 0&{s_{p}^{N-1}}&{}\hfil&\vdots\\\
\vdots&\vdots&{}\hfil&{}\hfil\\\
0&0&\cdots&{s_{p}^{{N}-1}}\end{array}}\right].$ (6)

The $[n,l]$-th element of $(M-N+1)\times L_{i,{j_{k}}}$ convolution matrix
${\bf U}(\Delta\tau_{i,{j_{k}}})$ is defined as

$[{\bf
U}(\Delta\tau_{i,{j_{k}}})]_{n,l}=\nabla\Big{(}{n-1-\frac{{{\Delta\tau_{i,{j_{k}}}}}}{{{T_{s}}}}-\frac{{\varsigma_{i,{j_{k}}}^{l}}}{{{T_{s}}}}}\Big{)}.$
(7)

Based on the above definitions, we arrive at the vector observation model as
follows

${{\bf{y}}_{j_{k}}}=\sum\limits_{i\in{{\cal P}_{p}}}{{\bf{F}}\left({{\Delta
f_{i,{j_{k}}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({{\Delta\tau_{i,{j_{k}}}}}\right){{\bf{h}}_{i,{j_{k}}}}}+{{\bf{w}}_{i,{j_{k}}}},$
(8)

where
${{\bf{w}}_{i,{j_{k}}}}={\left[{w_{i,{j_{k}}}\left[0\right],w_{i,{j_{k}}}\left[1\right],\ldots,w_{i,{j_{k}}}\left[{{M}-1}\right]}\right]^{T}}$
is the vector of noise, and ${{\bf{h}}_{i,{j_{k}}}}$ is a vector that collects
the coefficients of ${L_{i,{j_{k}}}}$ paths, denoted as
${{\bf{h}}_{i,{j_{k}}}}={\left[{h_{i,{j_{k}}}^{0},h_{i,{j_{k}}}^{1},\cdots,h_{i,{j_{k}}}^{{L_{i,{j_{k}}}}-1}}\right]^{T}}$.
The derivations of (8) can be found in Appendix A.

##  III Cramer-Rao bound with Pilot Contamination

Based on the given model, we jointly estimate the synchronization parameters
and derive the CRB for the CFO and TO. However, it is challenging to obtain
CFO and TO from all slave APs, and thus the master APs can only estimate the
synchronization parameters of the slave AP in the common cluster by treating
the pilot contamination from other clusters as noise. The received signal from
the $j_{k}$-th master AP can be rewritten as

$\begin{split}{{\bf{y}}_{j_{k}}}&={\bf{F}}\left({{\Delta
f_{i_{j_{k}}^{p},j_{k}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({{\Delta\tau_{i_{j_{k}}^{p},j_{k}}}}\right){{\bf{h}}_{i_{j_{k}}^{p},j_{k}}}+\underbrace{\sum\limits_{i\in{\left\\{{{{\cal
P}_{p}}\backslash i_{j_{k}}^{p}}\right\\}}}{{\bf{F}}\left({{\Delta
f_{i,j_{k}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right){{\bf{h}}_{i,j_{k}}}}+{{\bf{w}}_{i,j_{k}}}}_{{\bm{\bar{w}}}_{i,j_{k}}},\end{split}$
(9)

where $i_{j_{k}}^{p}$ is the slave AP using the $p$-th pilot sequence in
cluster ${\cal C}_{k}$ 111Similar to the frequency multiplexing, it is
impossible to reuse the common pilot sequence in the same cluster, which is
discussed in Section IV..

###  III-A Cramer-Rao Bound

In this subsection, we derive the CRB for the synchronization parameters with
pilot contamination. However, owing to the unknown TO and CFO from the slave
APs that share the common pilot sequence, it is challenging to derive the
exact distribution. By assuming that the pilot contamination follows the
complex Gaussian distribution 222Owing to no cell boundary and a large number
of APs, a pilot sequence needs to be shared lots of times, and thus the pilot
contamination can be approximated to follow a multivariate Gaussian
distribution, by using the central limit theorem., we adopt the well-known
conclusion for the Fisher matrix of the Gaussian probability density function
to derive CRB.

In our model, we denote $\bm{\Xi}=\text{var}\\{{\bm{\bar{w}}}_{i,{j_{k}}}\\}$
as the variance of ${\bm{\bar{w}}}_{i,{j_{k}}}$ and define
$\bm{m}(\bm{\theta})={\bf{F}}\left({{\Delta
f_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({{\Delta\tau_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{h}}_{i_{j_{k}}^{p},{j_{k}}}}$,
where the elements of $\bm{\theta}$ are given by

$\theta_{1}=\frac{{{\Delta\tau_{i_{j_{k}}^{p},{j_{k}}}}}}{{{T_{s}}}},\theta_{2}={\Delta
f_{i_{j_{k}}^{p},{j_{k}}}}T_{s},\theta_{2l+3}=\text{Re}\\{h_{i_{j_{k}}^{p},{j_{k}}}^{l}\\},\theta_{2l+3}=\text{Im}\\{h_{i_{j_{k}}^{p},{j_{k}}}^{l}\\},\forall
l\in\\{0,...,{{L}_{i_{j_{k}}^{p},{j_{k}}}-1}\\}.$ (10)

Then, the $[n,m]$-th element of the Fisher matrix can be expressed as

$[{\bf
J}]_{n,m}(\bm{\theta})=2\text{Re}{\Big{\\{}(\frac{\partial\bm{m}(\bm{\theta})}{\partial\theta_{n}})^{H}\bm{\Xi}^{-1}(\frac{\partial\bm{m}(\bm{\theta})}{\partial\theta_{m}})\Big{\\}}}.$
(11)

Based on the above discussion, we then derive the variance of pilot
contamination by assuming that the TO and CFO follow the uniform distribution
of $\frac{{{\Delta\tau_{i,{j_{k}}}}}}{{{T_{s}}}}\in[0,\eta]$,
$\eta\in{\mathbb{N}}^{+}$, 333Since the coarse synchronization is performed
periodically by synchronization protocol, the residual TO drift within several
chips. and ${\Delta f_{i,{j_{k}}}}\in[-f_{\max},f_{\max}]$,
$\forall{i\in{\left\\{{{{\cal P}_{p}}\backslash i_{j_{k}}^{p}}\right\\}}}$,
respectively.

Obviously, the mean of ${{\bm{\bar{w}}}_{i,{j_{k}}}}$ is $\bm{0}$ as all
variables are independent distribution with zero mean. The variance of
${{\bf{F}}\left({{\Delta
f_{i,j_{k}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right){{\bf{h}}_{i,j_{k}}}}$,
$\forall{i\in{\left\\{{{{\cal P}_{p}}\backslash i_{j_{k}}^{p}}\right\\}}}$,
can be expressed as

$\begin{split}&\text{var}\\{{{\bf{F}}\left({{\Delta
f_{i,j_{k}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right){{\bf{h}}_{i,j_{k}}}}\\}\\\
=&\mathbb{E}\\{{{\bf{F}}\left({{\Delta
f_{i,j_{k}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right){{\bf{h}}_{i,j_{k}}}}({{\bf{F}}\left({{\Delta
f_{i,j_{k}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right){{\bf{h}}_{i,j_{k}}}})^{H}\\}\\\
=&\mathbb{E}\\{{{\bf{F}}\left({{\Delta
f_{i,j_{k}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right){{\bf{h}}_{i,{j_{k}}}}}{{\bf{h}}^{H}_{i,{j_{k}}}}{\bf{U}}^{H}\left({{\Delta\tau_{i,j_{k}}}}\right){{\bf{S}}_{p}^{H}}{\bf{F}}^{H}\left({{\Delta
f_{i,j_{k}}}}\right)\\}.\end{split}$ (12)

In the following, we aim to derive (12). To alleviate the pilot contamination,
the distance among the slave APs that share the common pilot sequence should
be as far as possible, resulting in severe path loss and lower received signal
power. Furthermore, since the APs are deployed in relatively high positions,
the APs can transmit signal without blockage issues. Based on these reasons,
the pilot contamination from other clusters predominantly depends on the main
path with path delay $\varsigma_{i,{j_{k}}}^{0}=\frac{{{d_{i,{j_{k}}}}}}{c}$,
$\forall{i\in{\left\\{{{{\cal P}_{p}}\backslash i_{j_{k}}^{p}}\right\\}}}$,
where $d_{i,{j_{k}}}$ is the distance between the $i$-th slave AP and the
${j_{k}}$-th master AP and $c$ is the speed of light. For this special case,
${\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right)\in{\mathbb{R}}^{(M-N+1)\times
1}$ is a vector and ${\bf h}_{i,j_{k}}\in{\mathbb{C}}^{1\times 1}$ is
$\alpha_{i,j_{k}}^{0}\sqrt{{\beta_{i,j_{k}}}}\psi^{0}_{i,j_{k}}$. Then, owing
to the independent variables of ${\bf{F}}\left({{\Delta f_{i,j_{k}}}}\right)$
and ${\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right)$, we have

$\begin{split}&\text{var}\\{{{\bf{F}}\left({{\Delta
f_{i,j_{k}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right){{\bf{h}}_{i,{j_{k}}}}}\\}\\\
=&\alpha_{i,{j_{k}}}^{0}{{\beta_{i,{j_{k}}}}}\mathbb{E}\Big{\\{}{{\bf{F}}\left({{\Delta
f_{i,j_{k}}}}\right){{\bf{S}}_{p}}\mathbb{E}\\{{\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right)}{\bf{U}}^{H}\left({{\Delta\tau_{i,j_{k}}}}\right)\\}{{\bf{S}}_{p}^{H}}{\bf{F}}^{H}\left({{\Delta
f_{i,j_{k}}}}\right)\Big{\\}}.\end{split}$ (13)

Next, we derive the result of
$\mathbb{E}\\{{\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right){\bf{U}}^{H}\left({{\Delta\tau_{i,j_{k}}}}\right)\\}$,
denoted as $\bf A$. Owing to the TO between the $i$-th slave AP and the
${j_{k}}$-th master AP, we have the following three cases.

Case 1:
${0\leq{{\Delta\tau_{i,j_{k}}}}<1-\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}}}$,
we define
${\bf\bar{A}}(0)={\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right){\bf{U}}^{H}\left({{\Delta\tau_{i,j_{k}}}}\right)$
and have

$\begin{split}{\bf\bar{A}}(0)&={\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right){\bf{U}}^{H}\left({{\Delta\tau_{i,j_{k}}}}\right)={[{0,1-x_{i,j_{k}}^{0},x_{i,j_{k}}^{0},0,\ldots,0}]^{T}}[{0,1-x_{i,j_{k}}^{0},x_{i,j_{k}}^{0},0,\ldots,0}],\end{split}$
(14)

where $x_{i,{j_{k}}}^{0}$ is
$x_{i,{j_{k}}}^{0}={{\Delta\tau_{i,j_{k}}}}+\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}}$.
As can be seen from (14), there are only four non-zero elements, which are
denoted as

$\begin{split}[{\bf\bar{A}}(0)]_{2,2}&=(1-x_{i,{j_{k}}}^{0})^{2},[{\bf\bar{A}}(0)]_{3,3}=(x_{i,{j_{k}}}^{0})^{2},[{\bf\bar{A}}(0)]_{2,3}=[{\bf\bar{A}}(0)]_{3,2}=(1-x_{i,{j_{k}}}^{0})(x_{i,{j_{k}}}^{0}).\end{split}$
(15)

Case 2:
${\upsilon-1-\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}}\leq{{\Delta\tau_{i,j_{k}}}}<\upsilon-\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}}}$,
$\forall\upsilon\in\\{2,\cdots,{\eta}\\}$, we denote
${\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right){\bf{U}}^{H}\left({{\Delta\tau_{i,j_{k}}}}\right)$
as ${\bf\bar{A}}(\upsilon-1)$ and have

$\begin{split}{\bf\bar{A}}(\upsilon-1)&={\bf{U}}\left({{\Delta\tau_{i,j_{l}}}}\right){\bf{U}}^{H}\left({{\Delta\tau_{i,j_{l}}}}\right)\\\
&={[{\underbrace{0,\cdots,0}_{\upsilon},\upsilon-
x_{i,j_{k}}^{0},x_{i,j_{k}}^{0}+1-\upsilon,\ldots}]^{T}}[{\underbrace{0,\cdots,0}_{\upsilon},\upsilon-
x_{i,j_{k}}^{0},x_{i,j_{k}}^{0}+1-\upsilon,\ldots}].\end{split}$ (16)

Similar to the result given in (15), there are only four non-zero elements,
which are denoted as

$\begin{split}[{\bf\bar{A}}(\upsilon-1)]_{\upsilon+1,\upsilon+1}&=(\upsilon-
x_{i,{j_{k}}}^{0})^{2},[{\bf\bar{A}}(\upsilon-1)]_{\upsilon+2,\upsilon+2}=(x_{i,{j_{k}}}^{0}+1-\upsilon)^{2},\\\
[{\bf\bar{A}}(\upsilon-1)]_{\upsilon+1,\upsilon+2}&=[{\bf\bar{A}}(\upsilon-1)]_{\upsilon+2,\upsilon+1}=(\upsilon-
x_{i,{j_{k}}}^{0})(x_{i,{j_{k}}}^{0}+1-\upsilon).\end{split}$ (17)

Case 3:
${\eta-\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}}\leq{{\Delta\tau_{i,j_{k}}}}<\eta}$,
we define
${\bf\bar{A}}(\eta-\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}})={\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right){\bf{U}}^{H}\left({{\Delta\tau_{i,j_{k}}}}\right)$
and have

$\begin{split}{\bf\bar{A}}(\eta-\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}})=&{\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right){\bf{U}}^{H}\left({{\Delta\tau_{i,j_{k}}}}\right)\\\
=&{[{\underbrace{0,\cdots}_{\eta+1},{\eta+1}-x_{i,j_{k}}^{0},x_{i,j_{k}}^{0}-\eta,\ldots}]^{T}}[{\underbrace{0,\cdots}_{\eta+1},{\eta+1}-x_{i,j_{k}}^{0},x_{i,j_{k}}^{0}-\eta,\ldots}]\end{split}.$
(18)

The non-zero elements of
${\bf\bar{A}}(\eta-\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}})$ are given
by

$\begin{split}[{\bf\bar{A}}(\eta-\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}})]_{\eta+2,\eta+2}&=(\eta+1-x_{i,{j_{k}}}^{0})^{2},[{\bf\bar{A}}(\eta-\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}})]_{\eta+3,\eta+3}=(x_{i,{j_{k}}}^{0}-\eta)^{2},\\\
[{\bf\bar{A}}(\eta-\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}})]_{\eta+2,\eta+3}&=[{\bf\bar{A}}(\eta-\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}})]_{\eta+3,\eta+2}=(\eta+1-x_{i,{j_{k}}}^{0})(x_{i,{j_{k}}}^{0}-\eta).\end{split}$
(19)

By using the integration and the above results,
$\mathbb{E}\\{{\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right){\bf{U}}^{H}\left({{\Delta\tau_{i,j_{k}}}}\right)\\}$
is given by

$\begin{split}{\bf
A}=&\mathbb{E}\\{{\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right){\bf{U}}^{H}\left({{\Delta\tau_{i,j_{k}}}}\right)\\}=\mathbb{E}\\{{\bf\bar{A}}(0)+\sum_{v=2}^{\eta}{\bf\bar{A}}(\upsilon-1)+{\bf\bar{A}}(\eta-\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}})\\}.\end{split}$
(20)

Based on these discussions, if $\eta>1$, the elements of ${\bf A}$ are given
by

$\begin{split}[{\bf
A}]_{2,2}={\frac{(1-{\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}}})^{3}}{3\eta}},[{\bf
A}]_{2,3}=[{\bf
A}]_{3,2}={\frac{{1-3{{\left({\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}}}\right)}^{2}}+2{{\left({\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}}}\right)}^{3}}}}{6\eta}},\end{split}$
(21) $\begin{split}[{\bf
A}]_{3,3}={\frac{{2-{{\left({\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}}}\right)}^{3}}}}{3\eta}},[{\bf
A}]_{\upsilon+1,\upsilon+2}=[{\bf
A}]_{\upsilon+2,\upsilon+1}=\frac{1}{6\eta},\forall\upsilon\in\\{2,\cdots,\eta\\},\end{split}$
(22) $\begin{split}[{\bf
A}]_{\upsilon+2,\upsilon+2}=\frac{2}{3\eta},\forall\upsilon\in\\{2,\cdots,\eta-1\\},\end{split}$
(23) $\begin{split}[{\bf
A}]_{\eta+2,\eta+2}=\frac{2-(1-{\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}}})^{3}}{3\eta},[{\bf
A}]_{\eta+2,\eta+3}=[{\bf
A}]_{\eta+3,\eta+2}={\frac{{3{{\left({\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}}}\right)}^{2}}-2{{\left({\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}}}\right)}^{3}}}}{6\eta}},\end{split}$
(24) $\begin{split}[{\bf
A}]_{\eta+3,\eta+3}={\frac{{{{\left({\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}}}\right)}^{3}}}}{3\eta}}.\end{split}$
(25)

For the special case of $\eta=1$, the $[3,3]$-th element of ${\bf A}$ is
different from the result given in (22), which is given by

$[{\bf
A}]_{3,3}={\frac{{1-3{{\left({\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}}}\right)}^{2}}+3{\frac{{\varsigma_{i,{j_{k}}}^{0}}}{{{T_{s}}}}}}}{3}}.$
(26)

Then, since the $[n,m]$-th element of ${\bf{F}}\left({{\Delta
f_{i,j_{k}}}}\right){{\bf{S}}_{p}}$ is $e^{j2\pi(n-1){\Delta
f_{i,{j_{k}}}}{T_{s}}}{s}_{p}^{n-m}$, we have

$\begin{split}&[{{\bf{F}}\left({{\Delta
f_{i,j_{k}}}}\right){{\bf{S}}_{p}}\mathbb{E}\\{{\bf{U}}\left({{\Delta\tau_{i,{j_{k}}}}}\right)}{\bf{U}}^{H}\left({{\Delta\tau_{i,j_{k}}}}\right)\\}]_{n,m}=\sum_{k^{\prime}=1}^{M-N+1}{[{\bf{F}}\left({{\Delta
f_{i,j_{k}}}}\right){{\bf{S}}_{p}}]_{n,k^{\prime}}[{\bf A}]_{k^{\prime},m}}\\\
=&\sum_{k^{\prime}=m-1}^{m+1}{[{\bf A}]_{k^{\prime},m}e^{j2\pi({n-1}){\Delta
f_{i,{j_{k}}}}{T_{s}}}{s}_{p}^{n-k^{\prime}}},\end{split}$ (27)

where $s^{n-k^{\prime}}_{p}$ is $0$ when $n-k^{\prime}<0$ and
$n-k^{\prime}>N-1$.

By substituting the result of (27) into (13), we have

$\begin{split}&[{{\bf{F}}\left({{\Delta
f_{i,j_{k}}}}\right){{\bf{S}}_{p}}\mathbb{E}\\{{\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right)}{\bf{U}}^{H}\left({{\Delta\tau_{i,j_{k}}}}\right)\\}{{\bf{S}}_{p}^{H}}{\bf{F}}^{H}\left({{\Delta
f_{i,j_{k}}}}\right)]_{n,m}\\\
=&\sum\limits_{k^{\prime}=1}^{M-N+1}{[{{\bf{F}}\left({{\Delta
f_{i,j_{k}}}}\right){{\bf{S}}_{p}}\mathbb{E}\\{{\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right)}{\bf{U}}^{H}\left({{\Delta\tau_{i,j_{k}}}}\right)\\}]_{n,k^{\prime}}[{{\bf{S}}_{p}^{H}}{\bf{F}}^{H}\left({{\Delta
f_{i,j_{k}}}}\right)]_{k^{\prime},m}}\\\
=&\sum\limits_{k^{\prime}=1}^{M-N+1}{\sum_{k_{1}=k^{\prime}-1}^{k^{\prime}+1}{[{\bf
A}]_{k_{1},k^{\prime}}e^{j2\pi({n-1}){\Delta
f_{i,{j_{k}}}}{T_{s}}}{s}_{p}^{n-k_{1}}}e^{-j2\pi(m-1){\Delta
f_{i,{j_{k}}}}{T_{s}}}{s}_{p}^{m-k^{\prime}}}\\\
=&\sum\limits_{k^{\prime}=2}^{M-N+1}{\sum_{k_{1}=k^{\prime}-1}^{k^{\prime}+1}{[{\bf
A}]_{k_{1},k^{\prime}}e^{j2\pi({n-m}){\Delta
f_{i,{j_{k}}}}{T_{s}}}{s}_{p}^{n-k_{1}}}{s}_{p}^{m-k^{\prime}}}.\end{split}$
(28)

Then, based on the integration of CFO $\Delta f_{i,{j_{k}}}$, we have

$\begin{split}\int_{-f_{\max}}^{f_{\max}}{\frac{e^{j2\pi(n-m)\Delta
f_{i,{j_{k}}}T_{s}}}{2f_{\max}}d_{\Delta
f_{i,{j_{k}}}}=\left\\{{\begin{array}[]{*{20}{c}}{1,}&{m=n}\\\ {0,}&{m\neq
n}.\end{array}}\right.}\end{split}$ (29)

As a result, $\text{var}\\{{{\bf{F}}\left({{\Delta
f_{i,{j_{k}}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({{\Delta\tau_{i,{j_{k}}}}}\right){{\bf{h}}_{i,{j_{k}}}}}\\}$
is a diagonal matrix, denoted as ${\bf D}_{i,{j_{k}}}$, whose $[n,n]$-th
element can be expressed as

$\begin{split}[{{\bf
D}_{i,{j_{k}}}}]_{n,n}&\buildrel\Delta\over{=}[\text{var}\\{{{\bf{F}}\left({{\Delta
f_{i,j_{k}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({{\Delta\tau_{i,j_{k}}}}\right){{\bf{h}}_{i,{j_{k}}}}}\\}]_{n,n}\\\
&=\alpha_{i,{j_{k}}}^{0}{{\beta_{i,{j_{k}}}}}\sum\limits_{k^{\prime}=2}^{M-N+1}{\sum_{k_{1}=k^{\prime}-1}^{k^{\prime}+1}{[{\bf
A}]_{k_{1},k^{\prime}}{s}_{p}^{n-k_{1}}}{s}_{p}^{n-k^{\prime}}}.\end{split}$
(30)

By substituting (30) into (13), we obtain $\bm{\Xi}$, which can be written as

$\begin{split}\bm{\Xi}=\sum\limits_{i\in{\left\\{{{{\cal P}_{p}}\backslash
i_{j_{k}}^{p}}\right\\}}}{{{\bf{D}}_{i,{j_{k}}}}}+{\sigma^{2}}{\bf{I}}_{M}.\end{split}$
(31)

In the following, we derive the partial derivatives of $\bm{\theta}$. Define
${\bf M}={\text{diag}}(0,1,\cdots,M)$ and let ${\bf
e}_{l}^{i^{p}_{j_{k}},{j_{k}}}$ denote the $l$-th column of unit matrix ${\bf
I}_{L_{i^{p}_{j_{k}},{j_{k}}}}$. Denote
${\bf{U}}^{\prime}\left({{\Delta\tau_{i^{p}_{j_{k}},{j_{k}}}}}\right)$ as the
partial derivation with respect to ${{\Delta\tau_{i^{p}_{j_{k}},{j_{k}}}}}$,
whose $l$-th column has two non-zeros elements equal to $-1$ at component
$\left\lceil{1+\frac{{{\Delta\tau_{i^{p}_{j_{k}},{j_{k}}}}}}{{{T_{s}}}}+\frac{{\varsigma_{i^{p}_{j_{k}},{j_{k}}}^{l}}}{{{T_{s}}}}}\right\rceil$
and equal to $1$ at component
$1+\left\lceil{1+\frac{{{\Delta\tau_{i^{p}_{j_{k}},{j_{k}}}}}}{{{T_{s}}}}+\frac{{\varsigma_{i^{p}_{j_{k}},{j_{k}}}^{l}}}{{{T_{s}}}}}\right\rceil$.
Based on these definitions, we have

$\frac{{\partial{\bf{m}}\left(\bm{\theta}\right)}}{{\partial\theta_{1}}}={\bf{F}}\left({\Delta{f_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{S}}_{p}}{\bf{U}}^{\prime}\left({\Delta{\tau_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{h}}_{i_{j_{k}}^{p},{j_{k}}}},$
(32)
$\frac{{\partial{\bf{m}}\left(\bm{\theta}\right)}}{{\partial\theta_{2}}}=j2\pi{\bf{MF}}\left({\Delta{f_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({\Delta{\tau_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{h}}_{i_{j_{k}}^{p},{j_{k}}}},$
(33)
$\frac{{\partial{\bf{m}}\left(\bm{\theta}\right)}}{{\partial\theta_{2l+3}}}={\bf{F}}\left({\Delta{f_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({\Delta{\tau_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{h}}_{i_{j_{k}}^{p},j_{k}}}{\bf
e}_{l}^{i_{j_{k}}^{p},{j_{k}}},\forall
l\in\\{0,1,\cdots,L_{i_{j_{k}}^{p},{j_{k}}}\\},$ (34)
$\frac{{\partial{\bf{m}}\left(\bm{\theta}\right)}}{{\partial\theta_{2l+4}}}=j{\bf{F}}\left({\Delta{f_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({\Delta{\tau_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{h}}_{i_{j_{k}}^{p},{j_{k}}}}{\bf
e}_{l}^{i_{j_{k}}^{p},{j_{k}}},\forall
l\in\\{0,1,\cdots,L_{i_{j_{k}}^{p},{j_{k}}}\\}.$ (35)

Finally, by defining the $M\times 2(L_{i_{j_{k}}^{p},{j_{k}}}+1)$ matrix
$\bm{\Omega}=[\frac{{\partial{\bf{m}}\left(\bm{\theta}\right)}}{{\partial\theta_{1}}},\frac{{\partial{\bf{m}}\left(\bm{\theta}\right)}}{{\partial\theta_{2}}},\cdots,\frac{{\partial{\bf{m}}\left(\bm{\theta}\right)}}{{\partial\theta_{2L_{i_{j_{k}}^{p},{j_{k}}}+2}}}]$,
the Fisher matrix can be written as

${\bf J}(\bm{\theta})=2\text{Re}\\{{\bm{\Omega}}^{H}\bm{\Xi}{\bm{\Omega}}\\}.$
(36)

Finally, the CRBs of CFO and TO are obtained by the first two diagonal
elements of ${\bf J}^{-1}(\bm{\theta})$. These CRBs can be regarded as the
theoretical lower bound of mean square errors (MSEs) between the real
parameters and estimated results. Furthermore, as given the result in (36), it
is worth emphasizing that the estimation accuracy is strongly related to the
pilot contamination, especially for the low noise power regime.

###  III-B Maximum Likelihood Estimation

In this subsection, we jointly estimate CFO and TO between the $i_{j_{k}}^{p}$
slave AP and the ${j_{k}}$-th master AP, assuming that the noise power
$\sigma^{2}$ and the path delays $\varsigma_{i_{j_{k}}^{p},{j_{k}}}^{l}$,
$\forall l\in\\{0,1,2,...,L_{i_{j_{k}}^{p},{j_{k}}}-1\\}$ are known 444Due to
the fact that the $i_{j_{k}}^{p}$-th slave AP and the ${j_{k}}$-th mater AP
are in the same cluster, the distances from the $i_{j_{k}}^{p}$-th slave AP to
the ${j_{k}}$-th mater AP is closer than that of the slave APs in other
clusters and the variation of multi-path delay is relatively slow compared to
inter-cluster ones. Furthermore, since the synchronization is performed
periodically in a short time interval and all APs are deployed in high
locations with relatively stable propagation paths, these factors lead to the
known path delays.. Then, the ML estimation can be obtain by minimizing the
square distance, which is given by

$\left\|{{{\bf{y}}_{j_{k}}}-{\bf{F}}\left({\Delta{f_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({\Delta{\tau_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{h}}_{i_{j_{k}}^{p},{j_{k}}}}}\right\|^{2}.$
(37)

For given $\frac{{\Delta{\tau_{i_{j_{k}}^{p},{j_{k}}}}}}{{{T_{s}}}}$ and
$\Delta{f_{i_{j_{k}}^{p},{j_{k}}}}{T_{s}}$, the estimation of ${\bf
h}_{i_{j_{k}}^{p},{j_{k}}}$ can be obtained by

$\begin{split}{\bf\hat{h}}_{i_{j_{k}}^{p},{j_{k}}}&=\Big{(}{\bf{U}}^{H}\left({\Delta{\tau_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{S}}_{p}^{H}}{\bf{F}}^{H}\left({\Delta{f_{i_{j_{k}}^{p},{j_{k}}}}}\right){\bf{F}}\left({\Delta{f_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({\Delta{\tau_{i_{j_{k}}^{p},{j_{k}}}}}\right)\Big{)}^{-1}\\\
&\quad\times{\bf{U}}^{H}\left({\Delta{\tau_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{S}}_{p}^{H}}{\bf{F}}^{H}\left({\Delta{f_{i_{j_{k}}^{p},{j_{k}}}}}\right){\bf
y}_{j_{k}}.\end{split}$ (38)

Then, by substituting the result of (38) into (37), we can estimate $\Delta
f_{i_{j_{k}}^{p},{j_{k}}}T_{s}$ and
$\frac{\Delta\tau_{i_{j_{k}}^{p},{j_{k}}}}{T_{s}}$ by maximizing the following
quadratic form

$\begin{split}&{\bf
y}_{j_{k}}{\bf{F}}\left({\Delta{f_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({\Delta{\tau_{i_{j_{k}}^{p},{j_{k}}}}}\right)\Big{(}{\bf{U}}^{H}\left({\Delta{\tau_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{S}}_{p}^{H}}{\bf{F}}^{H}\left({\Delta{f_{i_{j_{k}}^{p},{j_{k}}}}}\right){\bf{F}}\left({\Delta{f_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({\Delta{\tau_{i_{j_{k}}^{p},{j_{k}}}}}\right)\Big{)}^{-1}\\\
\times&{\bf{U}}^{H}\left({\Delta{\tau_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{S}}_{p}^{H}}{\bf{F}}^{H}\left({\Delta{f_{i_{j_{k}}^{p},{j_{k}}}}}\right){\bf
y}_{j_{k}}.\end{split}$ (39)

This can be regarded as maximizing the energy of the projection ${\bf
y}_{j_{k}}$ onto the column space of the signal matrix
${\bf{F}}\left({\Delta{f_{i_{j_{k}}^{p},{j_{k}}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({\Delta{\tau_{i_{j_{k}}^{p},{j_{k}}}}}\right)$.
As a result, the ML estimator is obtained by searching over a two-dimensional
sufficiently fine grid of points with respect to the variables $\Delta
f_{i_{j_{k}}^{p},{j_{k}}}T_{s}$ and
$\frac{\Delta\tau_{i_{j_{k}}^{p},{j_{k}}}}{T_{s}}$.

##  IV Synchronization based on Pilot Assignment

Based on the result given in (36), we devise a synchronization scheme based on
pilot sharing to minimize the sum of CRBs while simultaneously satisfying the
requirement for all APs with limited overhead.

###  IV-A Problem Formulation

In this subsection, we jointly optimize the cluster classification, pilot
overhead, and pilot sharing scheme to minimize the synchronization errors. By
denoting the synchronization requirement of the $i$-th slave AP and the
$j_{k}$-th master AP as ${\rm CRB}^{\rm Req}_{i,j_{k}}$,
$j_{k}\in{\mathcal{C}}_{k}$, $\forall i\in\\{{\mathcal{C}}_{k}\backslash
j_{k}\\}$, $k=\\{1,\cdots,K_{C}\\}$, the minimization of sum CRBs can be
formulated as

$\displaystyle\mathop{\min}\limits_{K_{C},{\mathcal{C}_{k},k\in\left\\{1,\cdots,K_{C}\right\\}},{\tau},{{{\cal
P}_{p}},p\in\left\\{{1,\cdots,{\tau}}\right\\}}}$
$\displaystyle\sum\limits_{k=1}^{K_{C}}{\sum\limits_{i\in\\{{\cal
C}_{k}\backslash j_{k}\\}}{{{\rm{CRB}}_{i,j_{k}}(\Delta
f_{i,j_{k}},\Delta\tau_{i,j_{k}})}}}$ s. t. $\displaystyle\qquad{\rm
CRB}_{i,j_{k}}(\Delta f_{i,j_{k}},\Delta\tau_{i,j_{k}})\leq{\rm
CRB}_{i,j_{k}}^{\rm Req},$ $\displaystyle\qquad\qquad
j_{k}\in{\mathcal{C}}_{k},\forall i\in\\{{\mathcal{C}}_{k}\backslash
j_{k}\\},k\in\\{1,\cdots,K_{C}\\},$ (40a) $\displaystyle\qquad\qquad
2K_{C}+\tau\leq L_{S},$ (40b)

where ${{{\rm{CRB}}_{i,j_{k}}(\Delta f_{i,j_{k}},\Delta\tau_{i,j_{k}})}}$ is
the theoretical lower bound for estimation errors of CFO and TO between the
$i$-th slave AP and the $j_{k}$ master AP, $\tau$ is the pilot overhead, and
$L_{S}$ is the total overhead for synchronization. Constraint (40a) means that
each AP should satisfy the minimal requirements for synchronization, and
constraint (40b) represents that the total overhead, including the pilot
overhead and information exchange, should be less than $L_{S}$.

As can be seen from (40), the synchronization performance will be better if
there are more master APs and only one slave AP. However, it is unrealistic to
deploy more master APs and only one slave AP in the cell-free mMIMO systems,
which leads to significant overhead for exchanging information and impractical
implementation. Furthermore, considering the limited channel resources and
high requirements of synchronization, it is challenging to strike a balance
between the synchronization performance and overhead, i.e., information
exchange and pilot overhead. To address this issue, we aim to devise a
synchronization scheme with low overhead while simultaneously satisfying the
requirements of each AP.

To address this NP-hard problem, we first need to classify the clusters
according to APs’ synchronization requirements, and then allocate the pilot
sequences to all slave APs. Therefore, Problem (40) can be simplified into two
sub-problems, including adaptive cluster classification problem and pilot
sharing scheme, which is depicted in Fig. 2.

Figure 2: Illustration for solving the sum of CRBs.

###  IV-B Adaptive Cluster Classification

In this subsection, we find a good way to classify the clusters with a good
set of master APs while considering the overhead of information exchange. By
ignoring the pilot allocation strategy, Problem (40) can be simplified into

$\displaystyle\mathop{\min}\limits_{K_{C},{\mathcal{C}_{k},k\in\left\\{1,\cdots,K_{C}\right\\}}}$
$\displaystyle\sum\limits_{k=1}^{K_{C}}{\sum\limits_{i\in\\{{\cal
C}_{k}\backslash j_{k}\\}}{{{\rm{CRB}}_{i,j_{k}}(\Delta
f_{i,j_{k}},\Delta\tau_{i,j_{k}})}}}$ s. t. $\displaystyle\qquad{\rm
CRB}_{i,j_{k}}(\Delta f_{i,j_{k}},\Delta\tau_{i,j_{k}})\leq{\rm
CRB}_{i,j_{k}}^{\rm Req},$ $\displaystyle\qquad\qquad
j_{k}\in{\mathcal{C}}_{k},\forall i\in\\{{\mathcal{C}}_{k}\backslash
j_{k}\\},l\in\\{1,\cdots,K_{C}\\},$ (41a) $\displaystyle\qquad\qquad
2K_{C}<L_{S}.$ (41b)

Obviously, minimizing the sum of CRBs can be obtained by simply having one
slave AP with the largest path gain, but this violates the constraint (41b)
owing to the heavy burden on information exchange. Conversely, deploying an
insufficient number of clusters may increase the average distance between
master AP and slave APs in each cluster, which leads to a significant decrease
in synchronization performance. Furthermore, even though the clusters are
classified, it is still worth exploring which AP can be chosen as the master
AP to minimize the CRB. To address these issues, an adaptive cluster
classification based on the synchronization requirements is proposed, and a
criterion that finds a good set of master APs is to minimize the sum of CRBs.

In the following, our first step is to investigate how to find a reasonable
number of clusters according to the stringent synchronization performance.
Since estimation errors decrease with the signal-to-noise ratio (SNR) (dB)
[17], the synchronization performance can be transformed into the requirement
of signal-to-interference-plus-noise ratio (SINR) (dB) by treating the pilot
contamination as a part of the noise, which can be expressed as

$\begin{split}{\rm CRB}_{i,j_{k}}(\Delta
f_{i,j_{k}},\Delta\tau_{i,j_{k}})\leq{\rm CRB}_{i,j_{k}}^{\rm
Req}\Rightarrow{\rm SINR}_{i,j_{k}}\geq{\rm SINR}_{i,j_{k}}^{\rm
Req}.\end{split}$ (42)

Then, based on the required SINR, we derive the size of each cluster, which
can be used to determine the number of clusters. Specifically, the maximum
intra-cluster distance can be found by relying on the relationship between the
distance and path loss. For the ideal case, i.e., no multiple paths and no
pilot contamination, the expectation of ideal SINR can be written as

${\rm SINR}^{\rm ideal}=10\log_{10}\Big{(}\frac{{\rm Tr}\\{{\bf D}^{\rm
ideal}\\}}{{M\sigma^{2}}}\Big{)},$ (43)

where ${\bf D}^{\rm ideal}$ is $\text{var}\Big{\\{}{{\bf{F}}\left({{\Delta
f_{i^{p}_{j_{k}},j_{k}}}}\right){{\bf{S}}_{p}}{\bf{U}}\left({{\Delta\tau_{i^{p}_{j_{k}},j_{k}}}}\right){{\bf{h}}_{i^{p}_{j_{k}},{j_{k}}}}}\Big{\\}}$,
whose proof is omitted owing to the same derivation of (30). Then, the maximum
intra-cluster distance is determined by the minimum ${{\rm MSE}}^{\rm
Req}_{i,j_{k}}$, $\forall i\in\\{{\mathcal{C}}_{k}\backslash j_{k}\\}$,
$\forall k\in\\{1,\cdots,K_{C}\\}$. For ease of expression, we denote the
${\rm SINR}^{\min}$ as the equivalent condition of the minimum ${{\rm
MSE}}^{\rm Req}_{i,j_{k}}$. By assuming that
$\alpha_{i^{p}_{j_{k}},{j_{k}}}^{0}$ is 1 and substituting (43) into (42), we
have

$\begin{split}&{\rm SINR}^{\rm ideal}\geq{\rm
SINR}^{\min}\Rightarrow\beta^{\min}\geq\frac{M\sigma^{2}10^{\frac{{\rm
SINR}^{\min}}{10}}}{{\rm Tr}\Big{\\{}{\mathbb{E}}\big{\\{}{\bf F}(\Delta
f_{i^{p}_{j_{k}},j_{k}}){\bf S}_{p}{\bf A}{\bf S}_{p}^{H}{\bf F}^{H}(\Delta
f_{i^{p}_{j_{k}},j_{k}})\big{\\}}\Big{\\}}}.\end{split}$ (44)

By using the three-slope channel model given in [2], the maximum intra-cluster
distance can be readily obtained by

$\begin{split}{\rm
dis}^{\max}={\rm{1}}{{\rm{0}}^{-\frac{2}{7}\log\frac{{M\sigma^{2}10^{\frac{{\rm
SINR}^{\min}}{10}}}}{{\rm Tr}\Big{\\{}{\mathbb{E}}\\{{\bf F}(\Delta
f_{i^{p}_{j_{k}},j_{k}}){\bf S}_{p}{\bf A}{\bf S}_{p}^{H}{\bf F}^{H}(\Delta
f_{i^{p}_{j_{k}},j_{k}})\\}\Big{\\}}}-\frac{L^{\rm loss}}{{35}}}},\end{split}$
(45)

where $L^{\rm loss}$ is the constant factor. Although we cannot directly
determine the number of clusters, we can use this maximum intra-cluster
distance to determine whether the clusters satisfy the synchronization
requirements.

Based on the above discussions, for a fixed number of cluster $K_{C}$, the
minimum sum of CRBs can be obtained by maximizing the sum of SINRs, which can
be expressed as

$\displaystyle\mathop{\max}\limits_{{\mathcal{C}_{k},k\in\left\\{1,\cdots,K_{C}\right\\}}}$
$\displaystyle\sum\limits_{k=1}^{K_{C}}{\sum\limits_{i\in\\{{\cal
C}_{k}\backslash j_{k}\\}}{{{\rm{SINR}}_{i,j_{k}}}}}$ s. t.
$\displaystyle\qquad{\rm SINR}_{i,j_{k}}\geq{\rm SINR}_{i,j_{k}}^{\rm
Req},j_{k}\in{\mathcal{C}}_{k},\forall i\in\\{{\mathcal{C}}_{k}\backslash
j_{k}\\},k\in\\{1,\cdots,K_{C}\\},$ (46a) $\displaystyle\qquad\qquad
2K_{C}<L_{S}.$ (46b)

Although it is impossible to derive a closed-form expression of ${\rm
SINR}_{i,j_{k}}$, the result of the pilot contamination reveals that the power
primarily depends on distance-based large-scale fading factors. Consequently,
the maximum sum of SINRs can be obtained by shortening the distance between
the slave AP and master AP. Then, the cluster can be classified by minimizing
the sum of intra-cluster distances while ensuring that the AP-to-centroid
distance in the cluster is no larger than the maximum intra-cluster distance
${\rm dis}^{\max}$, which can be written as

$\displaystyle\min\limits_{K_{C},{\cal C}_{1},\cdots,{\cal C}_{K_{C}}}$
$\displaystyle\sum_{k=1}^{K_{C}}{\sum_{i\in{\cal C}_{k}}}|(x_{k},y_{k})-{\rm
Loc}_{i}|,$ s. t. $\displaystyle\qquad|(x_{k},y_{k})-{\rm Loc}_{i}|\leq{\rm
dis}^{\max},\forall i\in{\mathcal{C}}_{k},k\in\\{1,\cdots,K_{C}\\},$ (47a)

where $(x_{k},y_{k})$ means the optimal 2-dimensional location of the $k$-th
cluster’s centroid that has the minimal sum of distances to all APs in the
$k$-th cluster, and ${\rm Loc}_{i}$ is the 2-dimensional location of the
$i$-th AP. However, for the given number of clusters $K_{C}$, it is still
challenging to obtain the optimal solution of Problem (47) since it is an the
NP-hard problem. To address this issue, we classify the clusters by using the
K-means algorithm, whose details are given in [26]. By combining the K-means
algorithm with the maximum intra-cluster distance, we can gradually classify
the clusters by gradually increasing the number of clusters $K_{C}$ until the
maximum intra-cluster distances of all clusters are no more than the value
given in (45).

Based on the given clusters, a good set of master APs should be found to
minimize the sum of CRBs. Similar to the clustering approach, the master AP
and the centroid of the cluster should have the same characteristics, thereby
maximizing the sum of path gains in the cluster and improving the
synchronization performance. Mathematically, a criteria that finds a good set
of master APs can be written as

${j_{k}}=\mathop{\arg\min}\limits_{{j_{k}}\in{{\cal
C}_{k}}}\left\\{{\sum\limits_{i\in\left\\{{{{\cal
C}_{k}}\backslash{j_{k}}}\right\\}}{\left\|{{{\rm{Loc}}_{{j_{k}}}}-{{\rm{Loc}}_{i}}}\right\|}}\right\\},\forall
k\in\left\\{{1,\cdots,{K_{C}}}\right\\}.$ (48)

By using this criterion of (48), the master AP has the minimum sum distances
to all slave APs, which leads to larger path gain and better synchronization
performance.

Based on the above discussions, we propose an adaptive cluster classification
based on the K-means algorithm, which is detailed in Algorithm 1. Initially,
we divide all APs in this area into $K_{C}$ clusters and find the maximum
intra-cluster distance between the master AP and the slave AP of each cluster.
Then, if the maximum inter-cluster distance of all clusters is less than the
defined value ${\rm dis}^{\max}$, the adaptive cluster classification based on
the K-means algorithm terminates. Otherwise, the number of clusters will
increase by 1, and then the clusters will be classified again until the
maximum intra-cluster distances of all clusters are less than the specific
value ${\rm dis}^{\max}$.

Algorithm 1 Adaptive Cluster Classification based on K-means Algorithm

1: Initialize the iterative number $n=1$ and the number of clusters
$K_{C}^{n}$ = 2;

2: Initialize the cluster based on the K-means algorithm and find the set of
master APs.

3: Calculate the maximum intra-cluster distance between the master AP and the
slave AP of each cluster, denoted as ${\rm dis}_{k,n}^{\max}$, $\forall
k\in\\{1,\cdots,K_{C}^{n}\\}$;

4: while ($\min\Big{\\{}{\rm dis}_{1,n}^{\max},\cdots,{\rm
dis}_{K_{C}^{n},n}^{\max}\Big{\\}}>{\rm dis}^{\max}$) do

5: Update $n=n+1$, $K_{C}^{n}=K_{C}^{n-1}+1$;

6: Update the cluster classification
${\mathcal{C}}^{n}=\Big{\\{}{\mathcal{C}}_{1}^{n}\cup\cdots{\mathcal{C}}_{K_{C}^{n}}^{n}\Big{\\}}$
with $K_{C}^{n}$, by using the K-means algorithm;

7: Find the set of master APs by using (48);

8: Calculate the maximum intra-cluster distance between the master AP and the
slave AP of each cluster, denoted as ${\rm dis}_{k,n}^{\max}$, $\forall
k\in\\{1,\cdots,K_{C}^{n}\\}$;

9: end while

###  IV-C Pilot Sharing Algorithm

Based on the given cluster ${\mathcal{C}}_{k}$, $\forall
k\in\\{1,\cdots,K_{C}\\}$, we jointly optimize the pilot overhead $\tau$ and
the pilot sharing scheme, ${\cal P}_{p}$, $\forall p\in\\{1,\cdots,\tau\\}$,
to minimize the estimation errors. Mathematically, minimizing the sum of CRBs
can be expressed as

$\displaystyle\mathop{\min}\limits_{{\tau},{{{\cal
P}_{p}},p\in\left\\{{1,\cdots,{\tau}}\right\\}}}$
$\displaystyle\sum\limits_{k\in\left\\{{1,2,\cdots,K_{C}}\right\\}}{\sum\limits_{i\in\\{{\cal
C}_{k}\backslash j_{k}\\}}{{{\rm{CRB}}_{i,j_{k}}(\Delta
f_{i,j_{k}},\Delta\tau_{i,j_{k}})}}}$ s. t. $\displaystyle\qquad{\rm
CRB}_{i,j_{k}}(\Delta f_{i,j_{k}},\Delta\tau_{i,j_{k}})\leq{\rm
CRB}_{i,j_{k}}^{\rm Req},$ $\displaystyle\qquad\qquad
j_{k}\in{\mathcal{C}}_{k},\forall i\in\\{{\mathcal{C}}_{k}\backslash
j_{k}\\},k\in\\{1,\cdots,K_{C}\\},$ (49a) $\displaystyle\qquad\qquad\tau\leq
L_{S}-2K_{C}.$ (49b)

As can be seen from (49), the optimal solution can be obtained if all slave
APs are allocated unique pilot sequences. However, this orthogonal pilot
strategy is impractical owing to significant pilot overhead. To address this
issue, we devise a pilot-sharing scheme by searching for the optimal solution
with the given clusters.

In the following, we first simplify the pilot-sharing scheme by transforming
it into a graph coloring problem. Specifically, two rules are defined to share
the pilot sequence. For the first rule, to reduce the pilot contamination, the
slave APs in the common cluster cannot share the same pilot sequence, which
can be formulated as

${\bf s}_{p}^{H}{\bf s}_{p^{\prime}}=0,\quad\forall i^{p}\quad{\rm
and}\quad\forall i^{p^{\prime}}\in{\mathcal{C}}_{k},\quad\forall
k\in\\{1,\cdots,K_{C}\\},\quad{\rm and}\quad i^{p}\neq i^{p^{\prime}},$ (50)

where ${\bf s}_{p}\in{\mathbb{C}}^{N\times 1}$ is the orthogonal pilot
sequence allocated to the $i^{p}$-th slave AP. By denoting ${\cal S}={\cal
S}_{1}\cup{\cal S}_{2}\cdots{\cal S}_{K_{C}}$, we define a binary matrix $\bf
B\in{\mathbb{N}}^{|\cal S|\times|\cal S|}$, where $|\cal S|$ denotes the
number of element in ${\cal S}$, to indicate whether the $i$-th device is a
potential candidate for sharing the common pilot sequence with the $i$-th
device. The $i$-th row and $i^{\prime}$-th column element of the matrix $\bf
B$ is given by

${b_{i,i^{\prime}}}=\left\\{{\begin{array}[]{*{20}{c}}{1,}&i\cup
i^{\prime}\notin{{\cal C}_{k,}}\forall
k\in\left\\{{1,\cdots,{K_{C}}}\right\\},{\text{and}}\quad i\neq i^{\prime}\\\
{0,}&{{\text{otherwise,}}}\end{array}}\right.$ (51)

$b_{i,i^{\prime}}=1$ means that the $i$-th slave AP and the $i^{\prime}$-th
slave AP are in a common cluster, and thus these two APs should be allocated
unique orthogonal pilot sequences. By contrast, $b_{i,i^{\prime}}=0$ denotes
that the $i^{\prime}$-th slave and $i$-th slave AP belong to different
clusters, and it is possible to share the common pilot sequence between the
$i$-th slave AP and the $i^{\prime}$-th slave AP.

The second rule is to allocate the pilot sequence in a fair manner.
Specifically, the $p$-th pilot sequence can only be reused no more than
$N^{\max}$ times, which can be written as

$|{\mathcal{P}}_{p}|\leq N^{\max},\forall p\in\\{1,\cdots,\tau\\}.$ (52)

In this way, the worst-case scenario where a single pilot sequence is reused
by too many devices can be avoided.

(a) Undirected graph based on given clusters, , where S-AP and M-AP represent
the slave AP and master AP, respecively.

(b) Pilot assignment based on the Dsatur algorithm, where the same color
indicates that the same pilot sequence is used.

(c) Pilot allocation after combining the Dsatur algorithm with the swap-
matching operation.

Figure 3: Undirected graph and pilot allocation schemes.

By using the above two rules, allocating the pilot sequence can be transformed
into the graph coloring problem, which can be effectively solved by using the
Dsatur algorithm [27]. As depicted in Fig. 3, it is worth noting that a pilot
allocation strategy based on the Dsatur algorithm can follow the above two
rules by using the lowest pilot overhead. However, this algorithm did not
consider the estimation errors. For example, based on the result given in
(30), the synchronization performance based on the pilot allocation strategy
given in Fig. 3 (b) is not optimal, as slave AP 5 and slave AP 2 that share
the same pilot sequence are in close proximity, which leads to severe pilot
contamination and poor synchronization performance. Therefore, the objective
function given in (49) cannot be minimized by simply using the Dsatur
algorithm. To address this issue, we introduce the switch matching and two-
sided exchange into our pilot sharing scheme.

In the following, we introduce the concept of the swap-matching operation.

Definition 1: Given the pilot sharing scheme ${\cal
P}=\\{{\mathcal{P}}_{1},\cdots,{\mathcal{P}}_{\tau}\\}$ with
$i\in{\mathcal{P}}_{n}$, $i^{\prime}\in{\mathcal{P}}_{m}$, and
$i\notin{\mathcal{P}}_{m}$, $i^{\prime}\notin{\mathcal{P}}_{n}$, a swap-
matching $\mu_{i,i^{\prime}}$ is defined by the function
${\mathcal{P}}_{n}(\mu_{i,i^{\prime}})={\mathcal{P}}_{n}\backslash i\cup
i^{\prime}$ and
${\mathcal{P}}_{m}(\mu_{i,i^{\prime}})={\mathcal{P}}_{m}\backslash
i^{\prime}\cup i$, and the exchanged pilot sharing scheme is defined as ${\cal
P}_{\mu_{i,i^{\prime}}}=\\{{\mathcal{P}}_{1},\cdots,{\mathcal{P}}_{n}(\mu_{i,i^{\prime}}),{\mathcal{P}}_{m}(\mu_{i,i^{\prime}}),{\mathcal{P}}_{\tau}\\}$.

To be more specific, the swap-matching means a matching generated by swap
operation 555The swap operation is a two-sided version of the “exchange”
considered in [28]., in which two slave APs using different pilot sequences
switch their pilot sequences while keeping all other pilot allocation
assignments the same.

However, considering the previously defined two rules, the swap operation may
not be allowed by all slave APs. Specifically, after this swap-matching
operation, if the common pilot sequence is shared by the slave APs in the same
cluster, then this operation is not permitted, which is detailed in the
following.

Definition 2: Given the pilot sharing scheme
$\\{{\mathcal{P}}_{1},\cdots,{\mathcal{P}}_{\tau}\\}$ with
$i\in{\mathcal{P}}_{n}$, $i^{\prime}\in{\mathcal{P}}_{m}$, and
$i\notin{\mathcal{P}}_{m}$, $i^{\prime}\notin{\mathcal{P}}_{n}$, a swap-
blocking pair $(i,i^{\prime})$ satisfies that:

(1): $\exists a,b\in{\mathcal{P}}_{m}(\mu_{i,i^{\prime}})$ or $\exists
a,b\in{\mathcal{P}}_{n}(\mu_{i,i^{\prime}})$, $a\cup b\in{\cal C}_{k}$,
$\forall k\in\\{1,\cdots,K_{C}\\}$.

Compared to the conventional swap-blocking pair, we do not care whether this
swap-matching operation, in which all slave APs are involved, is two-sided
exchange-stable. This is due to the fact that the two-sided exchange-stable
match can only converge to the locally optimal pilot sharing scheme 666Similar
to the extremum of a multivariate function, it is impractical to obtain the
globally optimal solution along the direction of one variable’s partial
derivative while holding the other variables constant. Therefore, given the
other pilot sharing scheme, the swap-matching operation can only maximize its
own benefits while ignoring the other slave APs’ synchronization performance..
To search the globally optimal matching with the given pilot overhead $\tau$,
a swap-matching operation based on a simulated annealing method is proposed.
Specifically, we will determine whether to update the pilot sharing scheme
with a probability $P_{a}$, which is given by

$P_{a}=\frac{1}{1+e^{-T\times[{\rm Obj}({\cal P}_{\mu_{i,i^{\prime}}})-{\rm
Obj}({\cal P}_{\mu})]}},$ (53)

where $T$ is a probability parameter. ${\rm Obj}({\cal
P}_{\mu_{i,i^{\prime}}})$ is denoted as the sum of SINR based on the given
pilot sharing scheme ${\cal P}_{\mu_{i,i^{\prime}}}$, which can be written as

${\rm{Obj}}({\cal
P}{{}_{{\mu_{i,i^{\prime}}}}})=\left\\{{\begin{array}[]{*{20}{c}}{\sum\limits_{k\in\left\\{{1,\cdots,{K_{C}}}\right\\}}{\sum\limits_{i_{1}\in\left\\{{{{\cal
C}_{k}}\backslash{j_{k}}}\right\\}}{{\rm{SIN}}{{\rm{R}}_{i_{1},{j_{k}}}},}}}&{{\rm{SIN}}{{\rm{R}}_{i_{1},{j_{k}}}}\geq{\rm{SINR}}_{i_{1},{j_{k}}}^{{\rm{Req}}}}\\\
{0,}&{{\rm{SIN}}{{\rm{R}}_{i_{1},{j_{k}}}}<{\rm{SINR}}_{i_{1},{j_{k}}}^{{\rm{Req}}}}.\end{array}}\right.$
(54)

As can be seen from the probability in (53), if the swapped pilot-sharing
scheme is better than the previous one, then there is a high probability that
the exchanged scheme will be accepted. Conversely, if any AP cannot satisfy
the synchronization requirements, then the pilot allocation scheme may hold
the same with high probability.

Algorithm 2 Pilot-Sharing based on Swap-Matching Algorithm

1: Initialize the first iteration number $n_{1}=1$, the maximum objective
function ${\rm Obj}_{\max}=0$, ${\rm Flag}=1$, and the probability parameter
$T$ is 0.01;

2: Initialize the maximum pilot overhead $\tau^{(n_{1})}=\max\\{{|{\cal
S}_{1}|},\cdots,{|{\cal S}_{K_{C}}|}\\}$, and the maximum number for reusing
pilot sequences $N^{(n_{1})}_{\max}=\left\lceil{\frac{{\left|{{{\cal
S}_{1}}\cup{{\cal S}_{1}}\cdots{{\cal
S}_{{K_{{C}}}}}}\right|}}{{\tau^{(n_{1})}}}}\right\rceil$;

3: Use the Dastur algorithm to allocate the pilot sequences, denoted as ${\cal
P}^{(n_{1})}=\\{{\cal P}_{1},\cdots,{\cal P}_{\tau^{(n_{1})}}\\}$;

4: while (${\rm Flag}==1$) do

5: if ($\tau^{(n_{1})}>L_{S}-2K_{C}$) then

6: ${\rm Flag}=0$;

7: else

8: Set the second iteration number $n_{2}$ to 1 and initialize the pilot
sharing scheme as ${\cal P}^{(n_{1},n_{2})}={\cal P}^{(n_{1})}$;

9: while $n_{2}\leq N_{2}^{\max}$ do

10: Update $n_{2}=n_{2}+1$, randomly choose the proper slave APs and pilot
sequence based on the two rules given in (50) and (52), and denote the
switched pilot sharing scheme as ${\cal P}^{{\rm sw}}$;

11: Calculate the objective functions for both two pilot sharing schemes,
${\cal P}^{(n_{1},n_{2}-1)}$ and ${\cal P}^{{\rm sw}}$, which are denoted as
${\rm Obj}^{(n_{1},n_{2}-1)}$ and ${\rm Obj}^{{\rm sw}}$, and update pilot
sharing scheme ${\cal P}^{(n_{1},n_{2})}$ by executing the swap-matching
operation with probability $P_{a}$;

12: if (${\rm Obj}^{{\rm sw}}>{\rm Obj}^{(n_{1},n_{2}-1)}$) then

13: Update the pilot sharing scheme ${\cal P}={\cal P}^{{\rm sw}}$ and ${\rm
Obj}_{\max}={\rm Obj}^{{\rm sw}}$;

14: else

15: Update the pilot sharing scheme ${\cal P}={\rm Obj}^{(n_{1},n_{2}-1)}$ and
${\rm Obj}_{\max}={\rm Obj}^{(n_{1},n_{2}-1)}$;

16: end if

17: end while

18: Update $n_{1}=n_{1}+1$, $\tau_{p}^{(n_{1})}=\tau^{(n_{1}-1)}+1$,
$N^{(n_{1})}_{\max}=\left\lceil{\frac{{\left|{{{\cal S}_{1}}\cup{{\cal
S}_{1}}\cdots{{\cal S}_{{K_{{C}}}}}}\right|}}{{\tau^{(n_{1})}}}}\right\rceil$;

19: Randomly select ($N^{(n_{1})}_{\max}-1$) slave APs with low SINRs, replace
their pilots with the $\tau^{(n_{1})}$-th pilot sequence, and denote the new
pilot sharing scheme as ${\cal P}^{(n_{1})}$;

20: end if

21: end while

In this way, we can keep tracking the optimal matching found in an iterative
manner in a feasible region, i.e., $\tau\leq L_{S}-2K_{C}$, even if the
utility of the current matching is not locally optimal, which is detailed in
Algorithm 2. Specifically, based on the coloring theory, in the $n_{1}$-th
iteration, the pilot overhead $\tau^{(n_{1})}$ is initialized as the maximum
number of slave APs in one cluster, and the maximum number for reusing one
pilot sequence is determined by the number of slave APs and the pilot
overhead, which is defined as

$N^{(n_{1})}_{\max}=\left\lceil{\frac{{\left|{{{\cal S}_{1}}\cup{{\cal
S}_{1}}\cdots{{\cal S}_{{K_{{C}}}}}}\right|}}{{\tau^{(n_{1})}}}}\right\rceil,$
(55)

where $\left\lceil\cdot\right\rceil$ means the ceiling function operation. If
the total overhead is less than $L_{S}$, we begin to search the optimal
solution. With the given pilot overhead $\tau^{(n_{1})}$, the swap-matching
operation is repeated until the number of iterations $n_{2}$ exceeds a
specified number $N_{2}^{\max}$. Furthermore, after each swap-matching
operation, the pilot sharing scheme with better performance is updated. In
this way, the optimal pilot sharing scheme with pilot overhead
$\tau^{(n_{1})}$ can be found based on the synchronization requirements. Then,
to keep tracking the optimal solution, a new pilot sequence will be introduced
and allocated to the slave APs with low synchronization performance.
Particularly, we update the number of iterations by $n_{1}=n_{1}+1$, the pilot
overhead by $\tau^{(n_{1})}=\tau^{(n_{1}-1)}+1$, and the number for reusing
pilot sequence by $N^{(n_{1})}_{\max}=\left\lceil{\frac{{\left|{{{\cal
S}_{1}}\cup{{\cal S}_{1}}\cdots{{\cal
S}_{{K_{{C}}}}}}\right|}}{{\tau^{(n_{1})}}}}\right\rceil$. Then, we randomly
select ($N_{\max}^{(n_{1})}-1$) slave APs with low SINRs, denoting their
indices as $\\{i^{p_{1}},\cdots,i^{p_{N_{\max}^{(n_{1})}-1}}\\}$ and their
pilot sequences as $\\{{p_{1}},\cdots,{p_{N_{\max}^{(n_{1})}-1}}\\}$,
${p_{m}}\neq{p_{m^{\prime}}}$ if $m\neq m^{\prime}$. Furthermore, to follow
the first rule, for given any two selected APs, $i^{p_{m}}$ and
$i^{p_{m^{\prime}}}$, $\forall
m,m^{\prime}\in\\{1,\cdots,{N_{\max}^{(n_{1})}-1}\\}$, $m\neq m^{\prime}$, two
APs are not in the common cluster, which is given by

$i^{p_{m}}\cup i^{p_{m}^{\prime}}\notin{\cal C}_{k},\forall
k\in\\{1,\cdots,K_{C}\\}.$ (56)

After that, the new pilot sequence is assigned to the selected slave APs, and
the pilot sharing scheme in the $n_{1}$-th iteration is updated, which is
given by

$\begin{split}{{\cal P}^{(n_{1})}_{{p_{m}}}}&=\left\\{{{{\cal
P}^{(n_{1}-1)}_{{p_{m}}}}\backslash i^{p_{m}}}\right\\},\forall
m\in\left\\{{1,\cdots,{N^{({n_{1}})}_{{{\max}}}}-1}\right\\},\\\ {{\cal
P}^{(n_{1})}_{\tau^{(n_{1})}}}&=\left\\{i^{p_{1}},\cdots,i^{p_{{N^{(n_{1})}_{\max}-1}}}\right\\}.\end{split}$
(57)

Finally, based on the given clusters $\\{{\cal C}_{1},\cdots,{\cal
C}_{K_{C}}\\}$, the optimal pilot sharing scheme is obtained by searching the
feasible solutions.

To find the optimal solution to Problem (40), the above two algorithms are
executed iteratively and the optimal solution is reserved, which is detailed
in Algorithm 3.

Algorithm 3 Resource Allocation for Problem (40)

1: Initialize iterative number $n=0$ and the number of clusters $K_{C}^{n}=2$
;

2: Initialize the $\tau^{n}=0$ and sum of CRBs as ${\rm CRBs}^{(n)}=+\infty$;

3: while ($\tau^{n}+2K^{n}_{C}\leq L_{S}$) do

4: Update $n=n+1$, execute Algorithm 1 with given $K_{C}^{n}$, and obtain the
clusters $\\{{\cal C}_{1},\cdots,{\cal C}_{K_{C}^{n}}\\}$;

5: Execute Algorithm 2 to obtain the pilot allocation scheme in the $n$-th
iteration, denoted as $\\{{\cal P}_{1},\cdots,{\cal P}_{\tau^{n}}\\}$;

6: Calculate the sum of CRBs, denoted as ${\rm CRBs}^{(n)}$;

7: Update $K_{C}^{n+1}=K_{C}^{n}+1$;

8: if (${\rm CRBs}^{(n)}<{\rm CRBs}^{(n-1)}$) then

9: Update cluster classification ${\cal C}=\\{{\cal C}_{1},\cdots,{\cal
C}_{K_{C}^{n}}\\}$ and pilot allocation ${\cal P}=\\{{\cal P}_{1},\cdots,{\cal
P}_{\tau^{n}}\\}$;

10: end if

11: end while

###  IV-D Algorithm Analysis

In this subsection, the complexity of our proposed algorithm is analyzed. For
the cluster classification algorithm, the average complexity is $O(L_{C}\times
K\times I_{c})$, where $I_{c}$ is the number of iterations. For the pilot-
sharing scheme, we need to search the optimal solution subject to the pilot
overhead budget. Furthermore, for a given pilot overhead, the swap-matching
operation needs to be executed $N_{2}^{\max}$ times. The complexity of the
pilot-sharing scheme is $O(N_{2}^{\max}\times I_{\tau})$, where $I_{\tau}$ is
the number of total iterations. Therefore, the complexity of our proposed
algorithm lies in the number of both iterations of cluster classification and
swap-matching operations.

##  V Simulation Results

In this section, the performance of our proposed pilot allocation strategy is
numerically evaluated and discussed.

###  V-A Simulation Parameters

There are $K$ APs that are uniformly positioned constellation points in a
smart factory of size $0.15\times 0.15$ ${\text{km}}^{2}$. The large-scale
fading factors are based on the Hata-COST231 propagation model of [2], where
the height of each AP is 15 m. Specifically, the three-slope path loss (dB)
can be expressed as

${\rm{P}}{{\rm{L}}_{m,k}}=\left\\{{\begin{array}[]{*{20}{l}}{\begin{array}[]{*{20}{l}}{L_{\rm{loss}}+35{{\log}_{10}}\left({{d_{m,k}}}\right),}\\\
{L_{\rm{loss}}+15{{\log}_{10}}\left({{d_{1}}}\right)+20{{\log}_{10}}\left({{d_{m,k}}}\right),}\\\
{L_{\rm{loss}}+15{{\log}_{10}}\left({{d_{1}}}\right)+20{{\log}_{10}}\left({{d_{0}}}\right),}\end{array}}&{\begin{array}[]{*{20}{c}}{{d_{m,k}}>{d_{1}}}\\\
{{d_{0}}<{d_{m,k}}\leq{d_{1}}},\\\
{{d_{m,k}}\leq{d_{0}}}\end{array}}\end{array}}\right.$ (58)

where $d_{m,k}$ is the distance between the $m$-th AP and the $k$-th AP,
$L_{\rm{loss}}$ is a constant of 140.7 (dB), while $d_{0}$ and $d_{1}$ are
0.01 km and 0.05 km, respectively. Furthermore, each AP’s transmission power
is 1 W.

###  V-B Derivation and Monte-Carlo Simulations

Figure 4: $[{{{\bf{D}}_{i,{j_{k}}}}}]_{n,n}$ based on Monte-Carlo simulation
and derived result. Figure 5: Power of Pilot Contamination based on Monte-
Carlo simulation and derived result.

To check whether our derivations can provide a more convenient expression for
the pilot sharing scheme, we first check the $[n,n]$-th element of ${\bf
D}_{i,j_{k}}$ based on Monte-Carlo simulations and derived result. Fig. 4
depicts the randomly chosen elements of ${\bf D}_{i,j_{k}}$ with different
CFOs, TOs, and pilot lengths. Monte-Carlo simulation is consistent with our
derived closed-form expression, which demonstrates the accuracy of our
derivations.

To investigate the relationship between the distance-based large-scale fading
factor and the power of pilot contamination, we depict the power of derived
pilot contamination and that based on Monte-Carlo simulations. Fig. 5 shows
the power with a pilot sequence $\bf s$ formed by various lengths $N$ in the
binary phase shift keying (BPSK) constellation. Obviously, our derived results
are consistent with Monte-Carlo simulations, demonstrating that our derived
result can provide an explicit expression for the latter pilot-sharing scheme.
Furthermore, we observe that the power of pilot contamination increases with
the pilot length $N$. This is due to the fact that the power is equal to the
trace of the matrix, which increases with the pilot length $N$. More
importantly, it is worth noting that the power of pilot contamination has no
significant relationship between the CFO and TO but significantly decreases
with the increased distance, which demonstrates the effectiveness of our
derivations.

###  V-C Cramer-Rao Bound and ML Estimation

Figure 6: CRB and ML estimation versus SNRs with path delay [0, 0.7049$T_{s}$,
2.1230$T_{s}$, 2.7063$T_{s}$], path gain [1, 0.8443,0.4119,0.3223], and large-
scale fading factor 5.16 $\times 10^{-10}$.

To verify our derivations, we first check the theoretical lower bound, i.e.,
CRB, and the MSE based on ML estimation. Fig. 6 depicts the synchronization
performance with various SNRs. Furthermore, to obtain the inverse matrix, the
large-scale fading factor of each path is normalized by dividing the maximum
large-scale fading factor from the slave AP to master AP. Apparently, the CRBs
and MSEs of ML estimation decrease with the increased SNR, especially for the
results without pilot contamination. This phenomenon arises due to the
diminished interference between the desired signal and noise when the SNR is
high, resulting in better synchronization performance. However, owing to the
pilot sharing, the lower bound based on the pilot contamination no longer
linearly decreases with SNR, but instead experiences a flat trend. This is due
to the fact that as the noise power decreases, the main factor affecting the
estimation performance is no longer the noise, but the pilot contamination.

###  V-D Performance of the Proposed Algorithm

Figure 7: Sum of CRBs versus the number of APs with $N=256$, $\eta=3$,
$f_{\max}=3$, ${\rm SINR}^{\min}=15$ dB.

To investigate the impact of clustering on the synchronization performance, we
assume that the $K$ APs are randomly distributed in this area and all slave
APs are allocated to a unique pilot sequence. Furthermore, to illustrate the
performance of our adaptive cluster classification, we also compare our
approach with the following two benchmarks:

Benchmark 1: All slave APs transmit pilot sequences directly to the CPU, whose
2-dimensional location is [0,0].

Benchmark 2: The clusters are classified by the K-means algorithm while the
master APs are randomly chosen in each cluster.

By averaging over 100 simulation results where the path delays and path gain
coefficients are randomly generated, Fig. 7 illustrates the sum of CRBs of all
slave APs. Obviously, the synchronization performance slightly increases with
the number of APs. Furthermore, it is worth noting that our approach is
superior to all benchmarks. This is due to the fact that the cluster
classification could shorten the distance between master and slave APs,
thereby improving the path gain and enhancing the estimation performance.
Additionally, in each cluster, the synchronization performance is further
improved by reducing the distance between master APs and slave APs.

(a) Synchronization performance.

(b) Synhchronization overhead.

Figure 8: Sum of CRBs versus the number of APs with $N=256$, $\eta=3$,
$f_{\max}=3$, ${\rm SINR}^{\min}=15$ dB, $L_{S}=[25,28,31,34,37,40]$.

To further investigate the performance of our pilot sharing algorithm, we also
compare our algorithm with the following pilot allocation schemes.

Scheme 1: All salve APs transmit orthogonal pilot sequences directly to the
CPU, whose 2-dimensional location is [0,0].

Scheme 2: The clusters are classified by the K-means algorithm while the
master APs are randomly chosen in each cluster, and the pilot sharing scheme
is obtained by using the Dastur algorithm.

Scheme 3: The clusters are classified by using our adaptive cluster algorithm,
and the pilot sharing scheme is obtained by using the Dastur algorithm.

Fig. 8(a) illustrates the averaged sum of CRBs under the case of various
numbers of APs, while Fig. 8(b) shows the averaging overhead of each scheme.
As can be seen from Fig. 8(a), the synchronization performance slightly
increases with the number of APs. For Scheme 1, owing to the lower path gain,
the estimation errors are unable to meet the synchronization requirements.
Even though sharing pilots causes pilot contamination between APs, the
distance among APs that use the common pilot sequence is maximized to meet the
synchronization requirements. Furthermore, benefiting from our proposed
algorithm, the estimation errors are comparatively stable in the demand range.

##  VI Conclusions

To investigate the impact of synchronous pilot sharing on the estimation
performance, we derived the CRB while considering the pilot contamination, and
the results revealed that the estimation error is strongly related to the
power of pilot contamination, especially for the high SNR regime. Furthermore,
a maximum likelihood algorithm was presented to estimate the CFO and TO among
the pairing APs. Then, to minimize the sum of CRBs, a synchronization strategy
based on the pilot-sharing scheme was devised by jointly optimizing the
cluster classification, synchronization overhead, and pilot-sharing scheme.
This NP-hard problem was simplified into two sub-problems, including cluster
classification and the pilot sharing scheme. Specifically, the clusters were
classified by combining the K-means algorithm with the minimum distance sum
criterion, while the pilot sharing scheme was obtained using swap-matching
operations. Simulation results validated the accuracy of our derivations and
demonstrated the effectiveness of the proposed scheme over other pilot-pair
algorithms.

##  Appendix A Proof of (8)

The $[k^{\prime},l]$-th element of
${{{\bf{S}}_{p}}{\bf{U}}(\Delta{\tau_{i,{j_{k}}}})}$ is given by

$\begin{split}{\left[{{{\bf{S}}_{p}}{\bf{U}}(\Delta{\tau_{i,{j_{k}}}})}\right]_{k^{\prime},l}}&=\sum\limits_{n=1}^{M-N+1}{{{\left[{{{\bf{S}}_{p}}}\right]}_{k^{\prime},n}}{{\left[{{\bf{U}}(\Delta{\tau_{i,{j_{k}}}})}\right]}_{n,l}}}\\\
&=\sum\limits_{n=1}^{M-N+1}{{{\left[{{{\bf{S}}_{p}}}\right]}_{k^{\prime},n}}\nabla(n-1-\frac{{\Delta{\tau_{i,{j_{k}}}}+\varsigma_{i,{j_{k}}}^{l}}}{{{T_{s}}}})}\\\
&=\sum\limits_{n=1}^{M-N+1}{s_{p}^{k^{\prime}-n}\nabla(n-1-\frac{{\Delta{\tau_{i,{j_{k}}}}+\varsigma_{i,{j_{k}}}^{l}}}{{{T_{s}}}})},\end{split}$
(59)

where $s^{k^{\prime}-n}_{p}$ is $0$ when $k^{\prime}-n<0$ and
$k^{\prime}-n>N-1$.

Based on the definition of $\nabla\left(t\right)$, we have

$\nabla\left(t\right)=\left\\{{\begin{array}[]{*{20}{c}}{t,}&{0<t\leq 1}\\\
{2-t,}&{1<t\leq 2}\\\ {0,}&{\rm otherwise}.\end{array}}\right.$ (60)

By substituting (60) into (59), we have

${\left[{{{\bf{S}}_{p}}{\bf{U}}(\Delta{\tau_{i,{j_{k}}}})}\right]_{k^{\prime},l}}=s_{p}^{k^{\prime}-{x_{1}}}\nabla({x_{1}}-1-\frac{{\Delta{\tau_{i,{j_{k}}}}+\varsigma_{i,{j_{k}}}^{l}}}{{{T_{s}}}})+s_{p}^{k^{\prime}-{x_{1}}-1}\nabla({x_{1}}-\frac{{\Delta{\tau_{i,{j_{k}}}}+\varsigma_{i,{j_{k}}}^{l}}}{{{T_{s}}}}),$
(61)

where ${x_{1}}$ is
$\left\lceil{1+\frac{{\Delta{\tau_{i,{j_{k}}}}+\varsigma_{i,{j_{k}}}^{l}}}{{{T_{s}}}}}\right\rceil$.

Then, letting $k^{\prime}-x_{1}=n$, $k^{\prime}\in\\{1,\cdots,M\\}$, we have

$\begin{split}{\left[{{{\bf{S}}_{p}}{\bf{U}}(\Delta{\tau_{i,{j_{k}}}})}\right]_{k^{\prime},l}}&=s_{p}^{n}\nabla(k^{\prime}-1-n-\frac{{\Delta{\tau_{i,{j_{k}}}}+\varsigma_{i,{j_{k}}}^{l}}}{{{T_{s}}}})+s_{p}^{n-1}\nabla(k^{\prime}-n-\frac{{\Delta{\tau_{i,{j_{k}}}}+\varsigma_{i,{j_{k}}}^{l}}}{{{T_{s}}}})\\\
&=\sum\limits_{n=0}^{N-1}{s_{p}^{n}\nabla(k^{\prime}-1-n-\frac{{\Delta{\tau_{i,{j_{k}}}}+\varsigma_{i,{j_{k}}}^{l}}}{{{T_{s}}}})}.\end{split}$
(62)

Finally, we substitute $m=k^{\prime}-1$ into (62) and complete the proof of
(8) by combining the results of ${\bf F}(\Delta f_{i,j_{k}})$ and ${\bf
h}_{i,j_{k}}$.

## References

  * [1] G. Interdonato, E. Björnson, H. Quoc Ngo, P. Frenger, and E. G. Larsson, “Ubiquitous cell-free massive MIMO communications,” _EURASIP J. Wireless Commun. and Netw._ , vol. 2019, no. 1, pp. 1–13, Dec. 2019.
  * [2] H. Q. Ngo, A. Ashikhmin, H. Yang, E. G. Larsson, and T. L. Marzetta, “Cell-free massive MIMO versus small cells,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 3, pp. 1834–1850, 2017.
  * [3] Q. Peng, H. Ren, C. Pan, N. Liu, and M. Elkashlan, “Resource allocation for uplink cell-free massive MIMO enabled URLLC in a smart factory,” _IEEE Trans. Commun._ , vol. 71, no. 1, pp. 553–568, 2022.
  * [4] C.-X. Wang, X. You, X. Gao, X. Zhu, Z. Li, C. Zhang, H. Wang, Y. Huang, Y. Chen, H. Haas _et al._ , “On the road to 6G: Visions, requirements, key technologies and testbeds,” _IEEE Commun. Sur. Tuts._ , vol. 25, no. 2, pp. 905–974, 2nd Quar. 2023.
  * [5] H. Q. Ngo, L.-N. Tran, T. Q. Duong, M. Matthaiou, and E. G. Larsson, “On the total energy efficiency of cell-free massive MIMO,” _IEEE Trans. Green Commun. Netw._ , vol. 2, no. 1, pp. 25–39, Mar. 2018.
  * [6] E. Björnson and L. Sanguinetti, “Making cell-free massive MIMO competitive with MMSE processing and centralized implementation,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 1, pp. 77–90, 2019.
  * [7] ——, “Scalable cell-free massive MIMO systems,” _IEEE Trans. Commun._ , vol. 68, no. 7, pp. 4247–4261, Jul. 2020.
  * [8] H. Yan and I.-T. Lu, “Asynchronous reception effects on distributed massive MIMO-OFDM system,” _IEEE Trans. Commun._ , vol. 67, no. 7, pp. 4782–4794, 2019.
  * [9] J. Li, M. Liu, P. Zhu, D. Wang, and X. You, “Impacts of asynchronous reception on cell-free distributed massive MIMO systems,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 10, pp. 11 106–11 110, 2021.
  * [10] H. Sallouha, A. Chiumento, S. Rajendran, and S. Pollin, “Localization in ultra narrow band IoT networks: Design guidelines and tradeoffs,” _IEEE Internet Things J._ , vol. 6, no. 6, pp. 9375–9385, 2019.
  * [11] H. V. Balan, R. Rogalin, A. Michaloliakos, K. Psounis, and G. Caire, “Airsync: Enabling distributed multiuser MIMO with full spatial multiplexing,” _IEEE/ACM Trans. Netw._ , vol. 21, no. 6, pp. 1681–1695, 2013.
  * [12] O. Abari, H. Rahul, D. Katabi, and M. Pant, “Airshare: Distributed coherent transmission made seamless,” in _Proc. IEEE Conf. Comput. Commun. (INFOCOM)_. IEEE, 2015, pp. 1742–1750.
  * [13] K. Alemdar, D. Varshney, S. Mohanti, U. Muncuk, and K. Chowdhury, “Rfclock: Timing, phase and frequency synchronization for distributed wireless networks,” in _Proc. Int. Conf. Mobile Computing and Netw._ , 2021, pp. 15–27.
  * [14] M. Rashid and J. A. Nanzer, “Frequency and phase synchronization in distributed antenna arrays based on consensus averaging and kalman filtering,” _IEEE Trans. Wireless Commun._ , vol. 22, no. 4, pp. 2789–2803, 2022.
  * [15] K. Matsuura, K. Shin, D. Kobuchi, Y. Narusue, and H. Morikawa, “Synchronization strategy for distributed wireless power transfer with periodic frequency and phase synchronization,” _IEEE Commun. Lett._ , vol. 27, no. 1, pp. 391–395, 2022.
  * [16] S. R. Mghabghab, A. Schlegel, and J. A. Nanzer, “Adaptive distributed transceiver synchronization over a 90 m microwave wireless link,” _IEEE Trans. Antenn. Propag._ , vol. 70, no. 5, pp. 3688–3699, 2022.
  * [17] R. Rogalin, O. Y. Bursalioglu, H. Papadopoulos, G. Caire, A. F. Molisch, A. Michaloliakos, V. Balan, and K. Psounis, “Scalable synchronization and reciprocity calibration for distributed multiuser MIMO,” _IEEE Trans. Wireless Commun._ , vol. 13, no. 4, pp. 1815–1831, 2014.
  * [18] Y. Feng, W. Zhang, Y. Ge, and H. Lin, “Frequency synchronization in distributed antenna systems: Pairing-based multi-CFO estimation, theoretical analysis, and optimal pairing scheme,” _IEEE Trans. Commun._ , vol. 67, no. 4, pp. 2924–2938, 2018.
  * [19] U. K. Ganesan, R. Sarvendranath, and E. G. Larsson, “Beamsync: Over-the-air synchronization for distributed massive MIMO systems,” _IEEE Trans. Wireless Commun._ , 2023.
  * [20] M. Attarifar, A. Abbasfar, and A. Lozano, “Random vs structured pilot assignment in cell-free massive MIMO wireless networks,” in _Proc. IEEE Int. Conf. Commun. Workshops (ICCW)_ , Jul. 2018, pp. 1–6.
  * [21] H. Liu, J. Zhang, X. Zhang, A. Kurniawan, T. Juhana, and B. Ai, “Tabu-search-based pilot assignment for cell-free massive MIMO systems,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 2, pp. 2286–2290, Feb. 2020.
  * [22] S. Buzzi, C. D’Andrea, M. Fresia, Y.-P. Zhang, and S. Feng, “Pilot assignment in cell-free massive MIMO based on the hungarian algorithm,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 1, pp. 34–37, Jan. 2021.
  * [23] H. Liu, J. Zhang, S. Jin, and B. Ai, “Graph coloring based pilot assignment for cell-free massive MIMO systems,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 8, pp. 9180–9184, Aug. 2020.
  * [24] W. Zeng, Y. He, B. Li, and S. Wang, “Pilot assignment for cell free massive MIMO systems using a weighted graphic framework,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 6, pp. 6190–6194, Jun. 2021.
  * [25] Q. Peng, H. Ren, M. Dong, M. Elkashlan, K.-K. Wong, and L. Hanzo, “Resource allocation for cell-free massive MIMO-aided URLLC systems relying on pilot sharing,” _IEEE J. Sel. Areas Commun._ , vol. 41, no. 7, pp. 2193–2207, Jul. 2023.
  * [26] T. Kanungo, D. M. Mount, N. S. Netanyahu, C. D. Piatko, R. Silverman, and A. Y. Wu, “An efficient k-means clustering algorithm: Analysis and implementation,” _IEEE Trans. Pattern Anal. Machine Intelligence_ , vol. 24, no. 7, pp. 881–892, 2002.
  * [27] C. Pan, H. Mehrpouyan, Y. Liu, M. Elkashlan, and N. Arumugam, “Joint pilot allocation and robust transmission design for ultra-dense user-centric TDD C-RAN with imperfect CSI,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 3, pp. 2038–2053, 2018.
  * [28] E. Bodine-Baron, C. Lee, A. Chong, B. Hassibi, and A. Wierman, “Peer effects and stability in matching markets,” in _International Symposium on Algorithmic Game Theory_. Springer, 2011, pp. 117–129.