Cell-Free Massive MIMO-OFDM Over Spatially Correlated Doubly Selective Channels

This paper investigates the impact of doubly selective fading channels on the design and performance analysis of cell-free massive multiple-input multiple-output (MIMO)-orthogonal frequency division multiplexing (OFDM) networks. Time and frequency selectivity, as well as the intercarrier interference (ICI) effects, are fully characterized, and the spatial channel cross-covariance matrices associated to each mobile station (MS)/access point (AP) pair and for every OFDM symbol/subcarrier conforming a particular resource block (RB) are derived. Assuming that these cross-covariance matrices are known at the APs, optimal minimum mean square error (MMSE)-based and suboptimal projected MMSE (PMMSE)-based double selectivity-aware channel estimators/predictors are designed. Spatial cross-covariance matrices and channel estimates/predictions are then used to design scalable double selectivity-aware partial-MMSE (P-MMSE), averaged P-MMSE (AP-MMSE) and statistically averaged P-MMSE (SAP-MMSE) combiners and precoders to be used during the uplink (UL) and downlink (DL) payload data transmission phases. Extensive numerical simulations are conducted to assess the performance of the proposed double selectivity-aware scalable cell-free massive MIMO systems under scenarios showing different degrees of time and frequency selectivity. Simulation results clearly show the performance advantages of double selectivity-aware channel estimation/prediction schemes and combining/precoding strategies when working in high-mobility multipath scenarios of practical interest. Moreover, it is shown that, even though suboptimal PMMSE/AP-MMSE and PMMSE/SAP-MMSE strategies perform relatively well under poor frequency and time selectivity conditions, they suffer from very significant performance losses with respect to the MMSE/P-MMSE scheme as frequency and/or time selectivity levels increase.


I. INTRODUCTION
Cellular massive multiple-input multiple-output (MIMO) systems [1], [2] have base stations (BSs) equipped with a very large number of antennas that simultaneously serve multiple mobile stations (MSs) on the same frequency-time resources. They are typically based on the availability of channel state information (CSI) at the BSs, which is acquired through time division duplex (TDD) operation and the transmission of The associate editor coordinating the review of this manuscript and approving it for publication was Wence Zhang . uplink (UL) pilot signals. Such a setting allows for very high spectral and energy efficiencies using simple linear signal processing at both the UL and downlink (DL) [2], [3]. Reminiscent of the network MIMO [4], [5], coordinated multipoint (CoMP) [6], distributed cooperation clustering (DCC) [7], [8] and, more recently, cloud radio access network (C-RAN) [9] architectures, cell-free massive MIMO [10], [11] has recently been proposed as a possible alternative to canonical cellular massive MIMO for beyond fifth generation (B5G) and sixth generation (6G) networks [12], [13], [14]. In a cell-free massive MIMO network, a very large number of geographically distributed access points (APs) cooperate to simultaneously serve a lower number of MSs on the same frequency-time resources. The APs are connected to central processing units (CPUs) via fronthaul links, which enable the coordination tasks. Moreover, making use of physical layer concepts drawn from the cellular massive MIMO concept, cell-free massive MIMO networks rely on a TDD protocol and exploit the transmission of UL pilot sequences to acquire CSI at the APs that is then used to design highly reliable UL combiners and DL precoders based on linear signal processing techniques [15].
Cellular and cell-free massive MIMO architectures have received a lot of research interest in the past decade and its theoretical foundations, as well as the key components of the channel estimation process, the UL and DL operation, and the spatial resource allocation algorithms, have already been settled (see, for instance, [3], [16], [17] for a detailed survey on these topics). Most of these results, however, have been substantiated assuming block-wise frequency-flat timeinvariant fading channels. That is, most of the aforementioned research works assume an absolutely constant channel fading affecting all the samples (or channel uses) making up a particular TDD resource block (RB). This is certainly a mathematically convenient assumption providing rather precise approximations in scenarios where the frequency-time resources comprising an RB are well within the limits of both the coherence bandwidth and the coherence time of the channel. In reality, however, wireless communication channels are frequency-selective and time-varying in nature due, on the one hand, to multipath propagation and, on the other hand, to the mobility of the transceivers and the surrounding objects. Moreover, supporting high-mobility multipath scenarios will be essential in future wireless networks aiming at providing ubiquitous connectivity. For instance, 5G new radio (5G NR) is touted as a pivotal enabling technology for the realization of connected and cooperative autonomous driving, high-speed train wireless communications support or air-ground network integration through the use of unmanned aerial vehicles (UAVs) [18], [19], [20]. In these scenarios, multipath propagation can lead to non-negligible levels of frequency selectivity. Moreover, the experienced Doppler frequencies, which may easily reach hundreds of Hertz (i.e., a substantial percentage of the subcarrier spacing in most of the operational orthogonal frequency division multiplexing (OFDM) standards) result in both intercarrier interference (ICI) and a significant degree of time selectivity [21], [22].
Among the very scarce previous art on the design and analysis of cellular massive MIMO networks on doubly-selective fading channels, we can only call attention to the works published by Yang et al. [23] and Gao et al. [24]. In these papers, the authors propose different channel estimation strategies exploiting either the learning capabilities of deep neural networks (DNNs) or the double-sparsity (i.e., beamspace and delay sparsity) characteristic of mmWave channels, respectively. The effects of channel aging (i.e., time selectivity) have also been investigated in the context of cellular massive MIMO in previous literature [25], [26], [27], [28]. For instance, in [25] Papazafeiropoulos investigated how the Doppler shift induced by the relative movement of MSs, as well as the phase noise due to non-ideal local oscillators, might impact the performance of a massive MIMO system due to channel aging, showing that the consequences of MS mobility are much more detrimental than those caused by oscillator phase noise even for moderate MS speeds. Moreover, reinforcing previous results presented with Kong et al. [26], he also showed that the transmit power of each MS can be scaled down at most by 1/ √ M (with M denoting the number of antennas at the BS), indicating that channel aging does not deteriorate the power scaling law. Kong et al. [26], as well as Truong and Heath in [27], also pointed out that the performance losses caused by channel aging can be compensated in part by applying CSI prediction. More recently, Yuan et al. [28] proposed leveraging the temporal channel correlation in massive MIMO networks with channel aging to design machine learning (ML)-based channel prediction schemes offering remarkable performance gains, over conventional systems, on the per-user achievable rates in both low and high mobility scenarios.
In the cell-free massive MIMO arena, the previous research work on this topic is even more sparse. As far as we know, there are no published papers on the performance analysis and design of cell-free massive MIMO networks over doubly selective fading channels. Moreover, the impact of channel aging and phase noise on the performance of these systems is only evaluated in [29] and [30]. In particular, Zheng et al. [29] derive novel closed-form UL and DL per-user achievable rate expressions that take aged CSI, as well as spatial correlation and pilot contamination, into account. In the UL, they consider large-scale fading decoding (LSFD) and matched filter (MF) receiver cooperation and small cell (SC) noncooperating systems. In the DL, the performance of both coherent and non-coherent transmission modes is analyzed. The authors prove that the cell-free massive MIMO network is more robust to channel aging than the SC system. In [30], Jiang et al. analyze the impact of channel aging and phase noise on the performance of zero-forcing precoding in cellfree massive MIMO scenarios. The authors claim that the exchange of CSI and precoded data via a fronthaul network introduces a significant delay that must be combated by using delay-aware strategies such as channel prediction.
Adding to previous research on cell-free massive MIMO, our contributions in this paper can be summarized as follows: • In contrast to most existing literature, which assumes an absolutely constant channel fading affecting all the samples (or channel uses) making up a particular TDD RB, this work addresses the impact of doubly selective fading channels on the design and performance analysis of cell-free massive MIMO-OFDM networks. The proposed formulation is a generalization of previous analytical approaches that can be particularized to any specific scenario of interest (e.g., time selective channels, frequency selective channels, non-selective channels). Note that when VOLUME 10, 2022 characterizing the time and frequency selectivity of the channel we are also implicitly characterizing the reciprocity error between the UL and DL channels.
• In classical cell-free massive MIMO networks, which typically assume non-selective channels, all the pilot symbols making up the training sequence allocated to a particular MS are assumed to be both intersymbol interference (ISI)-and ICI-free and affected by exactly the same channel response. The optimal channel estimation corresponding to a given MS is then performed at the APs by first projecting the received pilot samples onto the training sequence allocated to that particular MS and then implementing a minimum mean square error (MMSE) estimation process using the projected pilot signal samples. The transmission on a doubly selective fading channel, however, breaks the possible orthogonality among pilot sequences used by different MSs. Consequently, the projected pilot signal samples cannot be considered a set of sufficient statistics anymore and novel selectivity-aware channel estimation/prediction schemes must be devised. In this paper we characterize a double selectivity-aware MMSE channel estimator/ predictor that works on the non-projected received pilot samples. We also propose a low-complexity double selectivity-aware estimator/predictor that, inspired by the lower-dimensionality of the classical approach, processes the projected pilot received samples to obtain a suboptimal low-complexity projected minimum mean square error (PMMSE) estimation/prediction.
• Designed to be a scalable approximation of the optimal MMSE combining scheme, we characterize double selectivity-aware partial minimum mean square error (P-MMSE) combiners/precoders. The implementation of the selectivity-aware P-MMSE combiners/precoders, however, is based on computationally expensive mathematical calculations that have to be performed on a per frequency-time resource basis. Hence, in order to alleviate this computational burden, we also propose low-complexity statistically averaged suboptimal approaches that will be shown to perform remarkably well, specially under scenarios showing poor or even mild selectivity conditions.
• The performance of the proposed double selectivity-aware scalable cell-free massive MIMO systems is evaluated through extensive numerical Monte Carlo simulations under different scenarios. In particular, the normalized mean square error (NMSE) associated to the channel estimation/prediction process and the spectral efficiency provided by the system during both the UL and DL payload data transmission phases are fully assessed. More specifically, we consider the impact of time and frequency selectivity (i.e., Doppler shift associated to the mobility of the transceivers and the surrounding objects, and time dispersion associated to multipath propagation) on the implementation of channel estimation/prediction and combining/precoding strategies.
The rest of the paper is organized as follows. Section II describes the system and channel models. Section III is used to describe the UL training and payload data transmission phase, and the DL payload data transmission phase is described in IV. In Section V we present extensive numerical results, while our conclusions are recapped in Section VI.

II. SYSTEM AND CHANNEL MODELS
A. SYSTEM MODEL Let us consider a wideband cell-free massive MIMO-OFDM system where M geographically distributed APs, each equipped with an array of N antennas and connected to CPUs via error-free infinite capacity fronthaul links, cooperate to simultaneously serve K single-antenna MSs on the same frequency-time resources. In contrast to most seminal papers on cell-free massive MIMO, which assumed that all MSs were served by all APs in the network, a scalable user-centric DCC will be considered in this paper meaning that a particular MS k is only served by a subset M k = m k 1 , . . . , m k M k ⊆ {1, . . . , M } of M k ≤ M APs [17].
As shown in Fig. 1, the OFDM scheme is assumed to exploit N c useful subcarriers with a subcarrier spacing f . In words, transmissions between APs and MSs are organized in a TDD operation with frames, of duration T f , composed of time slots of duration T s . Each time slot is further arranged into two phases, namely, the UL phase (containing the UL training and the UL payload data transmission subphases), and the DL payload data transmission phase. During the UL and DL phases, as shown in Fig. 1, grids of N s × N u and N s ×N d frequency-time resources are organized in UL and DL RBs, respectively. The total number of RBs is equal to N RB = N c /N s and all the symbols located on the N c frequency resources (i.e., subcarriers) of a particular time resource on all RBs are simultaneously multiplexed on an OFDM symbol. Each OFDM symbol, with period T o , is formed by putting together a cyclic prefix (CP) of duration T CP , assumed to be longer than the maximum channel multipath delay to avoid inter-OFDM block interference (also known as ISI), and a data part of duration T = 1/ f . In order to clarify the mathematical notation that will be used in the following sections, a particular frequency-time resource will be indexed using the duple [q, n], where q ∈ {1, . . . , N c } and n ∈ {1, . . . , Nu + N d }. Note that the set of subcarriers making up the RB r can be indexed by the set Q r = {(r − 1)N s + 1, . . . , rN s }, where the RBs have been indexed by the set R = {1, . . . , N RB }, beginning with the RB containing the low frequency subcarriers and progressing up to the RB containing the high frequency subcarriers. Interestingly, in order to simplify notation, we can also use the index z ∈ Z = {1, . . . , N c (N u + N d )} to refer to a particular frequency-time resource of a time slot, where z = (q − 1)(N u +N d )+n and, equivalently, q = (z−1)/(N u +N d ) +1, and n = mod(z − 1, N u + N d ) + 1. Note, moreover, that can be used to denote the frequency of subcarrier q, where f 0 is the carrier frequency.

B. CHANNEL MODEL
The MSs are assumed to be able to move at high speeds and thus, the wideband massive MIMO channels among APs and MSs might experience spatially correlated doubly-selective fading. That is, they will be characterized by, on the one hand, fast fading parameters that, in contrast to the vast majority of previous works on this topic, cannot be assumed to be constant over neither the RB duration T s = T o (N u + N d ) nor the RB bandwidth B s = N s f and, on the other hand, largescale parameters that can be considered to be almost constant over time and frequency intervals much larger than T s and B s , respectively. The characteristics of wideband doubly-selective fading channels can be captured by adopting the typical tappeddelay-line representation used by most standardization bodies (see, for instance, [31] and references therein). Thus, the vector of signals describing the channel between MS k and AP m in the delay-time domain can be generically represented as where L denotes the number of propagation paths, the path gain vectors are modelled as h mkl ∼ CN (0, C mkl ), τ mkl and ν mkl represent, respectively, the delay and Doppler shift associated with the lth path, and δ(·) denotes the Dirac delta function. It is assumed in this paper that τ max = max τ mkl < T CP − βT o , where β is the roll-off parameter of the OFDM pulse shaping filter (note that this allows for the complete elimination of the ISI, also known as interblock interference (IBI)). The large-scale average gain of the lth propagation path between any of the antennas at the mth AP and MS k can be obtained as Correspondingly, the aggregate large-scale average channel gain between any of the antennas at the mth AP and MS k, including the pathloss and shadowing, is given by The dependency of the large-scale propagation gain on both the carrier frequency and the distance between APs and MSs will be fully specified in the numerical results section. Note that the path gain vectors characterizing the propagation paths between different AP/MS pairs are assumed to be uncorrelated, that is,

III. UL TRAINING AND PAYLOAD DATA TRANSMISSION PHASE
In the UL training and payload data transmission phase, a MS transmits the pilot symbols and UL payload data symbols over disjoint subsets of the frequency-time resources composing each of the RBs of a given time slot. Let us assume that N p out of the N u N s UL frequency-time resources making up any RB r ∈ R are used for pilot transmission. Frequency-time resources used for pilot transmission on this generic RB, which will be typically evenly spread over the corresponding frequency-time grid, will be indexed by the set Z p r = z p r1 , . . . , z p rN p , and the UL frequency-time resources used for data transmission will be indexed by the Denoting by s u k [q, n] the symbol (either a pilot symbol or a data symbol) transmitted by MS k on the frequency-time resource [q, n], the baseband OFDM signal transmitted by this particular MS can be generically expressed as where υ(t) is the OFDM pulse-shaping filter. Assuming the channel model described in (2), the baseband received signals VOLUME 10, 2022 at AP m can then be obtained as h mk l e j2π ν mk l t s u k (t − τ mk l ) + n um (t), (6) where n um (t) is a vector of independent identically distributed (iid) zero-mean complex valued low-pass equivalent additive white Gaussian noise (AWGN) stochastic processes. Let us now assume that N c samples of the received signals are taken in the time interval [nT o +t 0 , (n+1)T o ] with a sampling period 1/(N c f ) and T CP ≤ t 0 ≤ T o + 1/ f (this is tantamount to assume that the CP is removed and the received signals are sampled in the time interval where υ(t − nTo − τ mk l ) = 1 for all k ∈ {1, . . . , K } and l ∈ {1, . . . , L}). The µth sample vector in the nth OFDM symbol period, taken at time t = nT o + t 0 + µ/(N c f ), for any µ ∈ {0, N c − 1}, is given by where n um [µ, n] ∼ CN 0, σ 2 u I N . The discrete Fourier transform (DFT) of this set of N c vectors of samples provides the received signal vector on the frequency-time resource [q, n] as where η um [q, n] ∼ CN 0, σ 2 u I N . Now, using the identity and the received signal vector y um [q, n] can be rewritten in a more compact form as Equation (10)  For later mathematical convenience, the spatial crosscovariance matrix between two generic equivalent channel vectors h mk [q 1 , q 1 , n 1 ] and h mk [q 2 , q 2 , n 2 ] is defined at this point as Note, moreover, that E h mk [q 1 , q 1 , n 1 ]h H m k [q 2 , q 2 , n 2 ] = 0 for all (m, k) = (m , k ). The spatial cross-covariance matrices R mk [q 1 , q 1 , n 1 ; q 2 , q 2 , n 2 ] fall in the category of very large-scale parameters (i.e., vary at a much slower pace than the coherence interval) and thus, we assume that they can be estimated and made available wherever needed in the network [17]. In particular, these frequency-domain channel statistics can be obtained through the estimation of the power delay profile (PDP) using techniques such as those proposed, for instance, in [32], [33], and [34] and references therein.

A. UL TRAINING: CHANNEL ESTIMATION
The K MSs simultaneously transmit pilot sequences of N p symbols to the APs on each RB. Let us denote by ϕ k = ϕ k1 . . . ϕ kN p T ∈ C N p ×1 the training sequence allocated to MS k, where ϕ ki is the pilot symbol transmitted over the frequency-time resource z p ri ∈ Z p r , that is, s u k [q p ri , n p ri ] = P p ϕ k [q p ri , n p ri ] = P p ϕ ki , with ϕ k 2 = 1 and P p 1 Without loss of generality, from this point onwards we redefine the reference time so that t 0 = 0. denoting the power allocated by MSs to each pilot symbol. Due to simultaneous transmission of pilot sequences on each RB, there will be a set of subcarriers making up the OFDM symbol n p ri , which will be indexed using the set Q[n p ri ], that will convey pilot symbols assumed to be known at the APs. The other subcarriers of this particular OFDM symbol will convey data symbols, that is s u k [q, n p ri ] = √ p k d u k [q, n p ri ] for all q ∈ Q[n p ri ], where d u k [q, n p ri ] ∼ CN (0, 1) and p k is the transmit power allocated by MS k to each data symbol. Thus, the N × N p received UL pilot signal matrix at the mth AP on RB r is where the received pilot signal vector y um [q p ri , n p ri ] in (10) can be expressed as p k d uk q , n p ri h mk q p ri , q , n p ri The first term of y um [q p ri , n p ri ] is the desired pilot signal and the second and third terms are the ICI components caused, respectively, by pilot and data symbols transmitted on subcarriers other than q p ri .

1) MMSE CHANNEL ESTIMATION
Let us consider a Bayesian MMSE channel estimator, implemented under the assumption that the large-scale channel parameters (i.e., spatial covariance matrices) are perfectly known at the APs [35]. Let us first define the vector of received pilot signal samples as Using this definition, and considering the user-centric nature of the cell-free massive MIMO network under evaluation, the MMSE estimate of the channel vector describing the propagation between frequency-time resource [q , n] at MS k and frequency-time resource [q, n] of RB r at AP m ∈ M k can be derived as 2 [36], [37] where mkr [q, q , n] and mr are given in (17)-(20), respectively, as shown at the bottom of the next page, with δ ij denoting the Kronecker delta function.
is the channel estimation error, which can be shown to be statistically independent ofĥ mk [q, q , n] and dis- where the spatial cross-covariance matrix betweenh mk [q, q 1 , n] andh mk [q, q 2 , n] is defined, for later mathematical convenience, as

2) PMMSE CHANNEL ESTIMATION
It is well known that the optimal MMSE channel estimator for MS k under time and frequency flat fading channels can be implemented by first projecting the N × N p received UL pilot signal matrix Y p mr on the pilot sequence allocated to this particular MS. Subsequently, the MMSE estimator can be applied to the projected signal vector (see, for instance, [17] and references therein). Even though this low-complexity channel estimation approach is no longer optimal under either time of frequency selectivity, and owing to its much lower computational complexity, it might be worth exploring under which circumstances it would provide performance metrics similar to those obtained using the optimal MMSE approach. Thus, let us first obtain the projected pilot signal vector for MS k at a generic AP m ∈ M k (i.e., one of the APs serving MS k) aš 2 Note that for Doppler shifts of practical interest the coefficients ξ mk l [q − q] quickly decrease as the difference between subcarrier indexes (q − q) increases. Hence, even though it is possible to obtain estimates of the ICI channelsĥ mk [q, q , n] for all q = q, in order to keep the computational complexity at affordable levels, only a small subgroup of neighboring ICI subcarriers needs to be considered in a practical channel estimation process. For instance, assuming the use of a carrier frequency f 0 = 5 GHz, a subcarrier spacing f = 15 kHz, and vehicular speeds as high as v k = 300 km/h, the quotient between the maximum Doppler shift f d = v k f 0 c and the subcarrier spacing is less than 0.1 and, as it was shown by Cai and Giannakis in [22], when f d / f = 0.1, more than 98% of the energy allocated to a symbol transmitted on subcarrier q is distributed on this subcarrier and its two neighboring subcarriers. Thus, even though all analytical expressions used in this paper consider the whole set of neighboring subcarriers when calculating the ICI terms, in a practical system suffices to only consider a very small subset of neighboring subcarrier terms. In fact, all simulation results presented in this paper have been obtained by considering only two neighboring subcarriers when calculating the ICI.
Using this vector, the PMMSE estimate of the channel vector describing the propagation between MS k and AP m on the frequency-time resource [q, n] of RB r can be derived as 3 [36], [37] wherě 3 With a slight abuse of notation and in order to simplify the mathematical expressions used in the following subsections, the same variable is used to represent the channel vector estimates for both the MMSE and PMMSE channel estimators. andˇ mkr = E y p mkry × mr [q p ri , n p ri ; q p rj , n p rj ].
with mkr [q, q , n; q p ri ,  p k R mk [q p ri , q , n p ri ; q p ri , q , n p ri ] + σ 2 The received signal vector at the mth AP on a generic frequency-time resource [q, n] dedicated to UL payload data transmission can be expressed as in (10). The received signal vectors at the M APs in the network can be cooperatively exploited by using either centralized or distributed combining architectures. The maximal ratio (MR), the partial-regularized zero-forcing (P-RZF), the P-MMSE or the fully centralized MMSE are among the large catalogue of combiners that have been already proposed in the available literature on cell-free massive MIMO to jointly process the signal vectors received at the APs (for a thorough overview on this topic, see [17] and references therein). Even though the semi-analytical approach we are presenting in this paper could be extended to any of the previously mentioned combiners, only the centralized scalable P-MMSE combiner will be considered owing to its very good performance in terms of spectral efficiency. In order to implement a centralized combiner, APs are assumed to send the received signals (including both the pilot signals (or channel estimates) and the received UL payload data vectors) to the CPUs. Using where The achievable spectral efficiency of MS k can then be obtained as (measured in bps/Hz) [16] where the quotient ( accounts for the fraction of useful time dedicated to UL payload data transmission, and the UL instantaneous effective signal-to-interference-plus-noise ratio (SINR) for MS k on the frequency-time resource [q, n] can be written as with p k R M k k [q, q, n; q, q, n] +R M k k [q, q , n; q, q , n] The effective SINR in (30) where S k denotes the set of MSs that are served by partially the same APs as MS k, that is, S k = k : M k ∩ M k = ∅ . The fractional power control algorithm proposed and motivated by Nikbakht et al. [38] for a cell-free network using MRC can be also adapted to a scalable cell-free massive MIMO scenario using an arbitrary combining scheme [17]. In particular, the kth MS is allocated an UL transmit power where the exponent ν ∈ [−1, 1] can be used to modify the power control behavior [17]. In particular, if ν = 0 all the MSs transmit at the maximum power. If ν = −1, in contrast, this heuristic power allocation can be viewed as an approximation to the max-min spectral efficiency fairness power control strategy. Selecting a proper value of ν > 0, an approximation to the maximum sum spectral efficiency scheme can be obtained. The P-MMSE combining vector described in (33) requires the calculation of the inverse of a matrix of dimension M k N × M k N for each frequency-time resource. In order to alleviate this computational complexity, we propose a low-complexity approach where the suboptimal combining vector for MS k on the frequency-time resource z of RB r is obtained as where u k is calculated by averaging u k [q, n] over a predefined set of frequency-time resources Z a , that is, This will be termed as the averaged P-MMSE (AP-MMSE) combiner.
Even though the simplification introduced in the AP-MMSE combiner has reduced the computational complexity associated with the calculation of inverse matrices in (33), the multiplications of the channel estimation vectors appearing in (34) also have a very high computational complexity and, consequently, it may also be worth exploring an even simpler approach that avoids the need to calculate these vector products. Towards this end, we propose the statistically averaged P-MMSE (SAP-MMSE) combiner defined as where the productsĥ M k k [q, q , n]ĥ H M k k [q, q , n] are approximated by its statistical expectation (i.e., the crosscovariance matrix), that is, × R M k k [q, q , n; q, q , n] with Note, moreover, that we have used the fact that R M k k [q, q, n; q, q, n] is independent of [q, n], that is, R M k k [q, q, n; q, q, n] = R M k k . Remarkably, the fronthaul signaling associated to different combiner designs will be very similar and mostly affected by the particular orchestration of the functional split between the CPUs and APs, an issue that is beyond the scope of our research work.

C. NOTES ON THE COMPLEXITY OF THE CHANNEL ESTIMATION AND COMBINING/DECODING PROCESSES
The MMSE channel estimation process described in (16), which is executed at each of the APs in the network, requires of the inversion of an (NN p × NN p ) matrix that is then multiplied by an (NN p × 1) vector, for each of the RBs in the system. 4 The global complexity of these processes can then be roughly approximated as O MN RB N 3 N 3 p . Each of the resulting (NN p × 1) vectors is next multiplied by a matrix of dimension (N × NN p ) for each triplet [q, q , n] and for each MS and, hence, the corresponding computational complexity can be approximated as O MN RB K (N u + N d )N 2 c N 2 N p . The complexity of these signal processing steps can be alleviated by implementing the PMMSE channel estimator described in (23). The PMMSE estimator, for each MS and each allocated RB, calculates the inverse of an (N × N ) matrix that is then multiplied by an (N × 1) (34) for the set of frequency-time resources in Z a . Moreover, it only calculates the inverse of the corresponding average matrix with an approximate computational complexity O KM 3 k N 3 . Finally, the SAP-MMSE precoding scheme does not compute the outer product of vectors in (34), thus eliminating the computational complexity associated to the calculation of these matrices.

IV. DL PAYLOAD DATA TRANSMISSION PHASE
Under analogous assumptions as those adopted in the UL, the CPUs perform a centralized precoding strategy in the DL. That is, all APs in M k jointly implement the precoding vector resource can be written as where D m = k m 1 , . . . , k m D m is the set of D m MSs served by the mth AP. As in the UL, the data symbols are modelled as iid random variables distributed as s d k [q, n] ∼ CN (0, 1) and, moreover, the precoding vectors are designed to fulfil the large-scale power constraint where P d is the average available transmit power per-symbol at the APs. The received data signal at the kth MS on the frequency-time resource z ∈ Z d of RB r can then be expressed as with η d k [z] ∼ CN (0, σ 2 d ) accounting for the noise term at the receiver. Now, assuming that the MSs have only access to statistical CSI, the achievable spectral efficiency of MS k (measured in bps/Hz) can be computed by exploiting the widely used hardening bound [11], [15], [16], [17], [40] as with SINR d k given in (45), as shown at the bottom of the next page, where the expectations are computed with respect to the fast fading terms and will be obtained by using Monte-Carlo simulations as in, for instance, [15], [40], [41], [42]. Assuming, as in the UL, the implementation of centralized joint processing at the CPU, a nearly optimal normalized P-MMSE or low-complexity normalized AP-MMSE and SAP-MMSE precoding schemes can be obtained from (33), (36) and (38), respectively, as [17] where υ k is a non-negative real value that determines the total transmit power allocated to MS k from all the serving APs in M k . Aiming at fulfilling the power constraint (42), a scalable heuristic power control strategy can be implemented where VOLUME 10, 2022 the transmit power allocated to MS k is given by [17] the largest fraction of υ k that any of the APs can be allocated to transmit, the parameter ν ∈ [−1, 1] can be used, as in the UL, to modify the power control behavior, and the exponent κ ∈ [0, 1] is a parameter that reshapes the power allocation ratio among different MSs (see [17,Section 7.2.3] for further details).

V. NUMERICAL RESULTS
A cell-free massive MIMO simulation scenario is considered where APs and MSs are independently and uniformly distributed at random within a square coverage area of side D = 500 m. In order to avoid boundary effects and to approximate an infinitely large network, the nominal squared coverage area is wrapped-around by eight identical neighbor replicas. The association of MSs to APs as well as the pilot assignment to MSs are jointly performed based on the usercentric DCC formation algorithm described by Demir et al. [17,Section 4.4]. In implementing the UL and DL scalable heuristic fractional power control algorithms, the exponents ν and κ have been set to −0.5 and 0.5, respectively, to strive for similar spectral efficiencies for all MSs in the network. The large-scale channel propagation gain (measured in dB) between the mth AP and the kth MS is modelled as [17, Sect.

2.5.2]
β mk [dB] = −30.5 − 36.7 log 10 d mk + χ mk , where d mk (measured in meters) is the three-dimensional distance between these devices, taking a height difference equal to 10 meters into account, and χ mk ∼ N (0, σ 2 χ ) is the shadow fading component with standard deviation σ χ = 4 dB. The spatial correlation between the shadow fading samples χ mk and χ m k is modeled as where δ kk is the distance (measured in meters) between MS k and MS k , and δ d is the decorrelation distance, which has been set to 9 meters.
The multipath doubly-selective fading channels used in the simulations have been modeled based on three of the power delay profiles specified by the third generation partnership project (3GPP) [43], [44], namely, the extended pedestrian A (EPA), the extended vehicular A (EVA), and the extended typical urban (ETU) models, representing low, medium, and high delay spread environments, respectively. However, assuming the use of small-sized APs equipped with a uniform linear array (ULA) with N antenna elements spaced a half wavelength, 5 the (i, j)th element of the spatial correlation matrix C mkl is computed as [17, Sect. 2.5.3] where ϕ mkl and ϑ mkl denote the azimuth and elevation angles of the lth multipath component of the channel between MS k and AP m with respect to the broadside of the array, and f (ϕ mkl , ϑ mkl ) is the joint probability density function of ϕ mkl and ϑ mkl , which is assumed to be a jointly Gaussian distribution with ϕ mkl and ϑ mkl denoting the average azimuth and elevation angles, and σ ϕ ≥ 0 and σ ϑ ≥ 0 denoting the corresponding angular standard deviations (ASDs) (see [17,Sect. 2.5.3] for more details). The average azimuth and elevation angles ϕ mkl and ϑ mkl are assumed to be uniformly distributed in the ranges [0, 2π] and [0, π], respectively. The Doppler shift associated with the lth propagation path of the channel between the MS k and the AP m is generated as ν mkl = v k f 0 c cos(ϕ mkl ) sin(ϑ mkl ). Table 1 is used to summarize the default system parameters used to set-up the scenarios evaluated in the forthcoming subsections. Notably, most of these default simulation parameters are reminiscent of those considered in the latest releases of the 3rd Generation Partnership Project (3GPP) fifth Generation (5G) New Radio (NR) access technology [46]. In obtaining the numerical results, we have used Monte Carlo simulations with 100 different scenarios (random locations of APs and MSs) and 200 smallscale fading realizations per scenario. 5 As stated by Demir et al. [17], for small-sized APs, which are envisioned to be used in cell-free massive MIMO networks, even though they could be equipped with uniform planar arrays (UPAs), it is common practice to utilize ULAs. The analysis presented in this paper could be readily extended to scenarios considering the use of UPAs by using the spatial correlation matrices defined in [45,Section VIII].
where the expectations are obtained with respect to both the large scale parameters characterizing the propagation channel, and the fast scatter fading components experienced on either the UL or the DL sets of frequency-time resources. Results presented in Fig. 2 have been obtained assuming a cell-free massive MIMO network with M = 40 APs, all equipped with ULAs of N = 4 antenna elements, and providing service to K = 12 MSs that are allocated training sequences with N p = 12 pilot symbols. 6 It can be clearly observed in Fig. 2a that the NMSE associated to the channel estimation/prediction process increases with the time selectivity (i.e., with the Doppler shift induced by the speed of MSs). Notably, the NMSE of the channel estimates/predictions corresponding to the DL payload data transmission phase is considerably higher than that associated to the estimation/prediction of the channel experienced during the UL payload data transmission phase. Obviously, this is due to the fact that the pilot sequences and UL payload data are simultaneously transmitted on the same set of OFDM symbols, while the DL payload data are transmitted 6 In this case, pilot symbols have been located according to the set Z p r = {30, 32,34,36,85,87,89,91,142,144,146, 148}, that is, pilot symbols are distributed as evenly as possible over the frequency-time resources making up the UL phase of the RB. on different OFDM symbols. Thus, although the temporal correlation of the UL and DL channels can be considerably high, as the Doppler shift becomes increasingly significant, the covariance matrices available at the APs provide the estimators/predictors with decreasingly reliable information about the future behavior of the channel. 7 It is also interesting to note that, even though for very low Doppler shifts (i.e., in scenarios with low time selectivity) the suboptimal PMMSE estimator/predictor provides performance metrics very similar to those afforded by the MMSE strategy, as the relative speed of the MSs increases, the advantage of the MMSE channel estimator/predictors becomes increasingly apparent. Just to give a numerical example, in scenarios where the relative speed of the MSs is equal to 3 km/h, the MMSE and PMMSE schemes provide an NMSE of roughly 0.01 and 0.02, respectively, both in the UL and the DL. In scenarios where the relative speed of MSs is equal to 120 km/h, in contrast, these same metrics are approximately equal to 0.014 and 0.073 in the UL, and 0.039 and 0.34 in the DL. Figure 2b depicts results allowing the evaluation of the impact the frequency selectivity might have on the performance of the proposed channel estimators/predictors. As previously mentioned, the EPA, EVA, and ETU channel models represent low, medium, and high delay spread environments, respectively and, as can be observed in this figure, the NMSE associated to the channel estimation/prediction process clearly increases with the frequency selectivity. This is especially the case when estimating the DL channels and/or when using the PMMSE suboptimal strategy. As a matter of fact, when transmitting over the EPA channel with MSs moving at speeds of 3 km/h, that is, in scenarios where both frequency and time selectivities are very low, the MMSE and PMMSE channel estimators provide virtually identical performance metrics in both the UL and the DL. When transmitting over an EVA channel with MSs moving at very low speeds, that is, when only frequency selectivity can be considered relevant enough, the PMMSE estimation/prediction scheme starts to show considerably worse performance than the MMSE scheme. As the velocity of MSs in this same EVA environment also increases, the performance of both estimation/prediction schemes worsens, but the performance of the PMMSE scheme deteriorates in a much more evident way, especially in the DL. The performance degradation of the channel estimation metrics becomes even more apparent when working in high frequency selective scenarios as those characterized by the ETU channel model. For instance, the NMSE measured in the EPA scenario with MSs moving at relative speeds of 3 km/h is approximately equal to 0.009 in both 7 Note that the precoding matrices that will be used at the APs during the DL payload data transmission phase will be designed using the channel predictions obtained from the pilot symbol signals received during the UL training phase; thus, the use of a DL training phase would certainly improve the data detection process but it would not improve at all the quality of neither the channels estimated/predicted in the UL training phase nor the precoding process implemented in the DL payload data transmission phase. the UL and the DL and irrespective of whether the MMSE or the PMMSE channel estimators are used. In contrast, the NMSE experienced in the ETU scenario with the MSs moving at relative speeds of 50 km/h is equal to approximately 0.014 and 0.082 in the UL and 0.024 and 0.16 in the DL, when using the MMSE and PMMSE channel estimators/predictors, respectively.

B. UL AND DL SPECTRAL EFFICIENCIES
In Fig. 3 we compare the cumulative distribution functions (CDFs) of the UL and DL per-user achievable rates obtained when implementing the P-MMSE, AP-MMSE and SAP-MMSE combining/precoding strategies These performance metrics are obtained by using either MMSE or PMMSE channel estimators and transmitting on an EVA scenario with the MSs moving at relative speeds of 120 km/h. The first interesting conclusion that can be drawn from these results is that the spectral efficiencies achieved in the UL are higher than those attained in the DL. The causes of this behavior are diverse, and we can highlight three of them. First, the number of frequency-time resources devoted to the UL payload data transmission phase is greater than the number of resources devoted to the DL payload data transmission phase. Second, as previously stated, the quality of the UL channel estimates/predictions and associated combining processes is better than the quality of the DL channel estimates/predictions and associated precoding processes. Third, the scalable fractional power control strategies used in conjunction with the centralized combining/precoding schemes allow for a more efficient use of the power available at the MSs than that available at the APs. Another noteworthy conclusion that can be drawn from Fig. 3 is that, in scenarios with moderate to high channel selectivity, cell-free massive MIMO systems based on the PMMSE channel estimation/prediction provide significantly lower performance metrics than those based on the MMSE strategy (we will put more emphasis on this result when commenting the graphs presented in the next figures). Finally, the most notable result to be drawn from this figure is that, when compared to the P-MMSE combining/precoding strategy, the performance metrics attained by the AP-MMSE and SAP-MMSE combining/precoding schemes behave in a very different way. In particular, the AP-MMSE scheme reliably copes with double selectivity in the UL but suffers from considerable performance losses with respect to the P-MMSE scheme in the DL. The computationally simpler SAP-MMSE strategy, however, suffers from substantial performance losses due to double selectivity in both UL and DL in comparison to the designs based on AP-MMSE and SAP-MMSE. For example, the P-MMSE, AP-MMSE, and SAP-MMSE schemes provide 95%-likely per-user spectral efficiencies of approximately 0.9, 0.88, and 0.76 bit/s/Hz in the UL and 0.82, 0.75, and 0.71 in the DL when using MMSE estimation/prediction. However, when using a PMMSE scheme, these same metrics go down to 0.82, 0.81, and 0.73 bit/s/Hz in the UL and 0.66, 0.61, and 0.57 in the DL.
In Fig. 4 we evaluate the impact of double selectivity on the CDF of the per-user UL and DL achievable rates using either the MMSE or the PMMSE channel estimators combined with the P-MMSE precoding/combining strategy. As in the graphs shown in Fig. 3, it can be observed that the spectral efficiencies achieved in the UL are considerably higher than those attained in the DL regardless of the degree of frequency and/or time selectivity characterizing the propagation channels. Another interesting conclusion that can be deduced from results presented in these graphs is that the increased frequency selectivity that the cell-free massive MIMO system experiences when transmitting on EPA, EVA or even ETU scenarios, produces a much lower negative impact than that generated by the increase of time selectivity occurring when the velocity of the MSs shoots up from 3 to 50, 120, or even 300 km/h. The final result worth drawing attention to is,  as expected, that cell-free massive MIMO networks based on the PMMSE channel estimator/predictor offer much lower performance metrics than those based on the MMSE strategy. Performance losses are particularly dramatic in high-mobility VOLUME 10, 2022 scenarios showing a high degree of time selectivity. In fact, just to give a quantitative example, the 95%-likely per-user spectral efficiencies in the UL of an MMSE-based system are, approximately, 0.95, 0.94, 0.9, and 0.8 bit/s/Hz when the MSs move, at 3, 50, 120, and 300 km/h, respectively, while these same metrics go down to 0.92, 0.89, 0.83, and 0.59 bit/s/Hz, respectively, when the system is based on the use of PMMSE estimators. The performance losses are even more exaggerated when comparing the results obtained in the DL.
Although not shown due to lack of space, we would like to note that results have also been obtained by varying the number of antennas available at the APs and the length of the training sequences. Increasing the number of antennas is known to significantly improve spectral efficiency in both the UL and the DL, the improvements it offers in the channel estimation/prediction process, however, are almost negligible. Increasing the length of the training sequences slightly improves the performance of the channel estimators and the spectral efficiency in the DL (assuming that the value of N d is kept fixed), but these improvements come at the cost of a decrease of the spectral efficiency metrics observed in the UL.

VI. CONCLUSION
In this paper we have addressed the impact of doubly selective fading channels on the design and performance analysis of cell-free massive MIMO-OFDM networks. We have characterized the equivalent channel between a particular frequency-time resource at a MS and any other frequency-time resource at an AP. Then, spatial channel cross-covariance matrices associated to each MS/AP pair and for every OFDM symbol/subcarrier conforming a particular RB have been derived. These cross-covariance matrices have then been used to design optimal MMSE-based and suboptimal PMMSE-based double selectivity-aware UL channel estimators/predictors, as well as double selectivityaware P-MMSE, AP-MMSE and SAP-MMSE combiners and precoders used during the UL and DL payload data transmission phases. Extensive numerical simulations have been used to evaluate the performance of the proposed double selectivity-aware scalable cell-free massive MIMO systems under different scenarios. One important conclusion to be drawn from the obtained results is that the deleterious effects caused by doubly selective fading channels must be taken into account in order to design double selectivity-aware channel estimation/prediction schemes and combining/precoding strategies in high-mobility multipath scenarios of practical interest. In particular, and roughly speaking, MS velocities above 50 km/h are prone to cause a significant performance loss if the time and frequency selectivity are not rightly considered. Another interesting conclusion is that, even though suboptimal channel estimation/prediction strategies based on projected pilot signal samples perform relatively well under poor frequency and time selectivity conditions, they suffer from very significant performance losses as frequency and/or time selectivity levels increase. Moreover, in realistic systems implemented using design parameters compatible with 5G NR standards, the effects of time selectivity are much more apparent than those caused by frequency selectivity. As a result, those are the ones that should be given more relevance when designing the channel estimation/prediction schemes. Regarding the double selectivity-aware combining/precoding strategies proposed in this paper, we can conclude that the AP-MMSE and SAP-MMSE combining/precoding strategies, when combined with the MMSE channel estimation scheme, can reliably cope with poor to mild double selectivity, particularly in the UL, providing performance metrics not far from those attained by the much more computationally expensive P-MMSE scheme. However, they suffer from considerable performance degradation with respect to the P-MMSE scheme, especially in the DL, when using PMMSE channel estimation/prediction strategies and transmitting in high mobility scenarios.