Performance Evaluation of CF-MMIMO Wireless Systems Using Dynamic Mode Decomposition

Abstract

Cell-Free Massive Multiple-Input–Multiple-Output (CF-MIMO) systems have transformed the landscape of wireless communication, offering unparalleled enhancements in Spectral Efficiency and interference mitigation. Nevertheless, the large-scale deployment of CF-MIMO presents significant challenges in processing signals in a scalable manner. This study introduces an innovative methodology that leverages the capabilities of Dynamic Mode Decomposition (DMD) to tackle the complexities of Channel Estimation in CF-MIMO wireless systems. By extracting dynamic modes from a vast array of received signal snapshots, DMD reveals the evolving characteristics of the wireless channel across both time and space, thereby promising substantial improvements in the accuracy and adaptability of channel state information (CSI). The efficacy of the proposed methodology is demonstrated through comprehensive simulations, which emphasize its superior performance in highly mobile environments. For performance evaluation, the most common techniques have been employed, comparing the proposed algorithms with traditional methods such as MMSE (Minimum Mean Squared Error), MRC (Maximum Ration Combining), and ZF (Zero Forcing). The evaluation metrics used are standard in the field, namely the Cumulative Distribution Function (CDF) and the average UL/DL Spectral Efficiency. Furthermore, the study investigates the impact of DMD-enabled Channel Estimation on system performance, including beamforming strategies, spatial multiplexing within realistic time- and delay-correlated channels, and overall system capacity. This work underscores the transformative potential of incorporating DMD into massive MIMO wireless systems, advancing communication reliability and capacity in increasingly dynamic and dense wireless environments.

1. Introduction

Cellular network deployments aim to handle the rapid increase in data traffic, with inter-cell interference being a significant challenge. Network multiple-input–multiple-output (MIMO) technology improves Spectral Efficiency by coordinating Access Points (APs). However, traditional methods require extensive sharing of channel state information (CSI) among APs, resulting in high signaling overhead and computational complexity. CF-MMIMO addresses these issues by distributing CSI tasks among APs, enabling coherent service to all User Equipment (UE) on a shared time–frequency resource. In this system, each AP is connected to a central processing unit (CPU) for signal coordination and processing. Unlike conventional distributed antenna systems, which operate within specific cells, cF-MMIMO eliminates the concept of defined cells, allowing all service antennas to coherently serve all User Equipment (UE).

The advent of Cell-free Multiple-Input–Multiple-Output (CF-MIMO) systems has been instrumental in redefining the capabilities of wireless networks [1,2], offering remarkable gains in data rates, coverage, and Spectral Efficiency. CF-MMIMO, an extension of MIMO and massive MIMO technologies [3], further amplifies these benefits by deploying an array of antennas at the base station, thereby enabling simultaneous communication with multiple users [4]. The synergy of CF-MMIMO and its dynamic communication environments hold the potential to revolutionize wireless systems, yet the complexity of accurate Channel Estimation in these scenarios remains a significant challenge.

In CF-MMIMO systems, the fundamental principle of exploiting spatial diversity and multiplexing gain is augmented by the presence of a large number of antennas at the base station [2,5]. This affords substantial improvements in the signal-to-interference-plus-noise ratio (SINR), enhancing both capacity and reliability. However, the practical implementation of CF-MMIMO introduces novel challenges, particularly in Channel Estimation, due to the sheer number of antennas and the dynamic nature of wireless channels [6].

The implementation of CF-MMassive MIMO is not a trivial problem; it has significant challenges. Coordinating and synchronizing a large number of geographically distributed antennas without a clear cell structure requires high precision in time and frequency synchronization to avoid interference and improve signal coherence. Implementing efficient and effective signal processing algorithms in a distributed manner across multiple antennas necessitates the development of robust algorithms capable of processing large volumes of data in real time. Managing the interference between users and antennas in a cell-free environment requires advanced beamforming techniques and interference cancellation methods. Accurately estimating communication channels in a dynamic and dispersed environment demands advanced Channel Estimation methods that can adapt to rapid variations in the channels. Providing a backhaul infrastructure capable of supporting the high volume of data exchanged between distributed antennas and the central processing unit needs high-capacity, low-latency backhaul links to ensure smooth and fast communication. Managing the energy consumption of a large number of distributed antennas requires energy-saving techniques and the optimization of resource usage. Developing and deploying the necessary hardware to support such a large and complex system needs advanced electronic components capable of operating in a distributed and large-scale environment. The costs associated with installing and maintaining such an extensive and complex infrastructure require significant investments and careful planning to ensure the economic viability of the system. Ensuring secure communications and maintaining user privacy in a widely dispersed environment necessitates the implementation of robust security protocols and efficient encryption mechanisms. These technical challenges underscore the complexity of implementing CF-MMIMO systems but also highlight the opportunities for technological innovations and improvements in wireless communications.

Traditional approaches to Channel Estimation often rely only on orthogonal pilot sequences to estimate the channel response [4]. However, these methods are susceptible to pilot contamination, leading to inaccurate channel estimates, especially in scenarios with a high density of antennas [7]. Some works suggest increasing the number of pilots to harness this issue, but the overhead incurred by pilot signaling escalates with the number of antennas, hampering the efficiency gains of massive MIMO.

In this paper, we propose a paradigm shift in Channel Estimation for CF-MMIMO systems by leveraging the power of DMD and two of its variants, in particular, Physical Informed Dynamic Mode Decomposition (PiDMD) and measure-preserving extended DMD (mpeDMD). DMD/PiDMD/mpeDMD are purely data-driven techniques that extract dynamic modes from time-series data, revealing the underlying dynamics of complex systems. Our approach emphasizes the ability of DMD/PiDMD/mpeDMD to capture the spatial–temporal evolution of the channel, offering a novel solution to the challenges posed by rapid variations in CF-MMIMO channels. Dynamic Mode Decomposition (DMD) is a data-driven technique utilized to decompose complex dynamical systems into spatial–temporal modes [8], effectively capturing the underlying patterns and dynamics from time-resolved data. PiDMD (Physics-informed DMD) [9] enhances this approach by incorporating prior knowledge of the governing physical laws, thereby improving the accuracy of the dynamic modes, particularly when data are limited or noisy. mpeDMD (Multi-parameter Extended DMD) further extends DMD by accommodating dynamical systems influenced by multiple parameters, allowing for a comprehensive analysis of how these parameters affect the observed dynamics [10].

This work also introduces a novel model for channel correlation [11], which is crucial in practical MIMO systems where correlation consistently exists. This model relies on delay–Doppler characteristics and draws inspiration from Orthogonal Time–Frequency Space (OTFS) [12]. To demonstrate the potential of DMD/PiDMD/mpeDMD in wireless networks, we will consider scenarios involving User Equipment (UE) with high mobility, which introduces the concept of channel aging [13].

In the remainder of this work, “DMD” will generally refer to DMD itself and its variants, including PiDMD and mpeDMD, unless otherwise specified.

For notational simplicity, we define a set of diagonal matrices $D_{kl}\in {\mathbb{C}}^{N\times N}$, where $k=1,\dots ,K$ and $l=1,\dots ,L$ represent the communication links between Access Points (APs) and User Equipment (UE). This set models link breakage as a uniformly distributed random process. Specifically, $D_{kl}$ is the identity matrix ${\mathbf{I}}_N$ if AP l is permitted to transmit to and receive signals from UE k and otherwise ${\mathbf{0}}_N$:

$$D_{kl}=\left\{\begin{array}{cc}{\mathbf{I}}_N\hfill & \mathrm{if}l\in {\mathcal{M}}_k,\hfill \\ {\mathbf{0}}_N\hfill & \mathrm{if}l\notin {\mathcal{M}}_k.\hfill \end{array}\right.$$

The elements of $\sqrt{{\eta }_k}{\varphi }_{t_k}$ are transmitted as the signal $s_k$ over ${\tau }_p$ consecutive samples. The received signal at AP l during the pilot transmission is denoted as $Y_l^{\mathrm{pilot}}\in {\mathbb{C}}^{N\times {\tau }_p}$ and is given by

$$Y_l^{\mathrm{pilot}}=\sum _{i=1}^K\sqrt{{\eta }_i}{h}_{il}{\varphi }_{t_i}^T+N_l$$

(1)

where $N_l\in {\mathbb{C}}^{N\times {\tau }_p}$ represents the receiver noise with i.i.d. elements distributed as ${\mathcal{N}}_{\mathbb{C}}(0,{\sigma }_{\mathrm{ul}}^2)$. The received uplink signal $Y_l^{\mathrm{pilot}}$ is used by AP l to estimate the channels $\{h_{kl}:l\in {\mathcal{M}}_k\}$ for all the UE it serves. This estimation can be performed either directly at AP l or at the CPU. In the latter case, AP l acts as a relay, forwarding the received pilot signal $Y_l^{\mathrm{pilot}}$ to the CPU via the fronthaul link. The CPU can then compute all the channel estimates $\{{\widehat{\mathbf{h}}}_{kl}:k=1,\dots ,K,l\in {\mathcal{M}}_k\}$ using the received pilot signals $\{{\mathbf{Y}}_l^{\mathrm{pilot}}:l\in {\mathcal{M}}_k\}$. However, since the channel vectors are independent, computing the channel estimates separately at each AP l in the set ${\mathcal{M}}_k$ does not affect optimality. The estimation results presented in this section apply to both approaches.

1.1. Related Work

The term “Cell-Free massive MIMO” began to gain prominence in the academic literature around 2017, as several articles were published exploring this concept and its applications in wireless communications [14]. Since then, there has been a growing interest and development in this research area, with numerous studies and advancements published in specialized journals and technical conferences. In the initial stages, Cell-Free MIMO emerged as a concept that seemed practically challenging to implement and theoretically difficult to resolve. A significant breakthrough was achieved in [5] by employing the Minimum Mean Square Error (MMSE) approach to address the massive MIMO problem. This foundational work set the stage for further explorations into the potential of cell-free architectures.

In [15], the modeling and analysis of cell-free massive MIMO (cF-MMIMO) systems are thoroughly examined. The study developed a mathematical framework to evaluate the performance and potential of these systems, focusing on the deployment of distributed antennas and emphasizing key metrics such as Spectral Efficiency, energy efficiency, and user fairness. Advanced signal processing techniques, such as MMSE and Zero Forcing (ZF), are incorporated to highlight the advantages over traditional cellular systems, including enhanced connectivity and reduced interference. Theoretical and simulation results underscore the practicality and benefits of cF-MMIMO in next-generation wireless networks. In [16], the application of deep learning techniques in cF-MMIMO systems is explored. This research investigated the unique challenges presented by these systems, such as efficient resource management, interference mitigation, and system performance optimization. By demonstrating how deep learning techniques can address design and optimization challenges—including resource allocation, spectrum management, and dynamic adaptation to changing communication environments—the study highlights the potential of deep learning to enhance efficiency, capacity, and flexibility in cF-MMIMO systems.

Further extending this line of inquiry, ref. [17] applied deep learning techniques for Channel Estimation in cF-MMIMO systems. This work investigated how deep neural networks can accurately estimate communication channels in cF-MMIMO environments by training models with channel data to predict conditions efficiently and accurately. The potential advantages of using deep learning techniques to enhance accuracy and efficiency in Channel Estimation are thoroughly explored, offering an innovative approach to optimizing performance and resource management in advanced wireless communication networks. Additionally, ref. [18] reviewed the application of DMD within the context of wireless communications, discussing how DMD can be used to analyze and model data generated by wireless systems, including channel characterization, signal prediction, resource optimization, and Spectral Efficiency enhancement.

In summary, the integration of Dynamic Mode Decomposition (DMD) into cell-free massive MIMO systems offers a transformative impact on wireless communications. The ability to predict and adapt to changing channel conditions opens up new avenues for optimizing transmission strategies, improving Spectral Efficiency, and ensuring reliable communication. Moreover, the use of advanced deep learning techniques further enhances the capacity and flexibility of these systems, paving the way for novel applications in wireless networks beyond 5G and 6G, where dynamic and agile communication environments are paramount. This comprehensive review of the literature underscores the significant progress made in this field and highlights the potential for future advancements and practical implementations in next-generation wireless networks.

1.2. Contributions

The contributions of this work are the following:

• Evaluating the application of DMD for Channel Estimation in CF-MMIMO wireless systems.
• Analyzing and valuating the behavior of DMD in a system that considers the correlation of the channel in terms of delay and Doppler.
• Evaluating the enhancement in overall system capacity achieved through DMD-based Channel Estimation.

1.3. Paper Structure

This paper is organized as follows.

Section 2, Materials and Methods, begins with an overview of the System Model, detailing the application of DMD for Channel Estimation in CF-MMIMO systems. It discusses the incorporation of channel correlation and introduces a novel proposed algorithm aimed at improving performance metrics.

Section 3, Channel Estimation, delves deeper into the methodologies employed for Channel Estimation using DMD techniques. It explores channel prediction and addresses the challenges posed by channel aging within CF-MMIMO environments.

Section 4, Uplink Transmission, focuses on the crucial components of uplink communication in CF-MMIMO systems. Topics covered include uplink training methods, Uplink Data Transmission strategies, and the evaluation of uplink Spectral Efficiency. Additionally, it examines key transmission techniques such as Maximum Ratio Combining (MRC), Minimum Mean Squared Error (MMSE) Channel Estimation, Full Zero Forcing (FZF), Partial Zero Forcing (PZF), and Optimal Least Squares Frequency Domain (LSFD) weights.

Section 5, Downlink Transmission, shifts attention to the complexities of Downlink Transmission. It contrasts Coherent and Non-Coherent Downlink transmission strategies, highlighting their respective advantages and challenges within the CF-MMIMO framework.

Section 6, Results, presents and analyzes empirical findings derived from the simulations and experiments conducted to validate the proposed methodologies and algorithms.

Section 7, Discussion, expands on the implications and insights drawn from the results, exploring broader implications and potential avenues for future research.

Section 8, Conclusion, offers a concise summary of the paper’s contributions, emphasizing how the adoption of DMD techniques in CF-MMIMO systems enhances Channel Estimation accuracy and transmission efficiency.

The CF-MMIMO architecture is illustrated in Figure 1. Several APs are randomly distributed within a circular area of radius D, accommodating N UEs. In principle, each UE communicates with every AP. A central CPU is capable of collecting data from each AP. This setup presents several challenges. For instance, the volume of data transmitted from the APs to the CPU grows substantially with the increasing number of UEs and APs. A centralized transmission strategy becomes impractical, and Channel Estimation becomes complex due to the large amount of data to be processed.

It is easy to realize that strategies that allow dealing with a huge amount of information while producing satisfactory results for good performance are highly desirable. This is precisely where DMD/PiDMD proves suitable for two main tasks: it optimally reduces the dimensions of the systems, and it leverages all available information, particularly the structure of the correlation matrix.

DMD methods are purely data-driven techniques. In general, data-driven methods in wireless networks lack accuracy because the number of samples is always small, especially for the estimation of the correlation matrix R, which is used for Channel Estimation in most methods. DMD is able, with a relatively small number of samples, to capture the spatial and temporal correlations of the wireless channel [8].

We assume that the base station ($B{S}_k$) and User Equipment ($U{E}_s$) are perfectly synchronized and operate using a pilot-based protocol. Each coherence block contains ${\tau }_p$ channel uses, with K designated for pilots. There are K orthogonal unit-norm pilot sequences, which are reused across cells. The pilot for UE k in cell j is represented by ${\varphi }_{jk}\in {\mathbb{C}}^K$ and satisfies ${∥{\varphi }_{jk}∥}^2=1$ [20].

2. Materials and Methods

2.1. Simulation Scenario

Our approach adheres to the standard massive MIMO protocol outlined in [21]. We adopt the notation from [5]; the system is illustrated in Figure 1, where a CF-MMIMO system operating under Time Division Duplexing (TDD) is used. All single-antenna User Equipment (UE) is simultaneously served by all Access Points (APs) (in this work Access Points and base stations (BSs) are exactly the same), L represents the number of APs, M is the number of antennas at each AP, and K is the number of UEs. The APs are linked to a central processing unit (CPU) through front haul connections. The communication channel undergoes continuous variations, segmented into ${\tau }_p$ channel uses for uplink pilots and ${\tau }_u$ for uplink data, with the total coherence interval given by ${\tau }_c={\tau }_p+{\tau }_u$ (Figure 2). We utilize a block-fading model, where $h_{kl}\in {\mathbb{C}}^{N\times 1}$ denotes the channel response between the k-th UE and the l-th AP, for $k=1,\dots ,K$ and $l=1,\dots ,L$. Each fading block is independently drawn from a Rayleigh distribution, as described in [19]. The number of time–frequency snapshots considered in this study is represented by ${\tau }_s$.

• ${\tau }_p$ symbols for uplink pilots;
• ${\tau }_u$ symbols for uplink data;
• ${\tau }_d$ symbols for downlink data.

The channel vector between $A{P}_l$ and $U{E}_k$ is denoted as $h_{l,k}\in {\mathbb{C}}^{M\times 1}$, incorporating both small-scale and large-scale fading characteristics. The following assumptions are made:

• The channel is reciprocal due to the precise calibration of the hardware chains, which can be achieved through standard methods [21].
• The channel is not static within a specific time–frequency interval, known as the coherence interval. It varies continually and independently between these intervals [3].
• The channel, $h_{l,k}\in {\mathbb{C}}^N$, follows a Rayleigh fading distribution, $\mathrm{CN}\sim (0,{\beta }_{l,k})$, where ${\beta }_{l,k}$ represents the constant large-scale fading coefficient (channel variance) between $A{P}_l$ and $U{E}_k$. This coefficient is consistent across antenna elements, meaning it does not depend on the antenna element index m, $m=1,\dots ,M$ [5].
• The large-scale fading coefficients are known in advance using DMD prediction at each AP. They change slowly and are influenced by UE mobility. Therefore, we assume that channel variances are estimated early in the process, and these estimates are then used to gauge the current channel response [11].
• Unlimited Front-Haul Network Capacity. The capability of CF-MMIMO with front-haul capacity limitations has been explored in [4,22].

The TDD frame consists of three phases [7]:

• Uplink (UL) pilot transmission, also known as UL training;
• Uplink Data Transmission; and
• Downlink (DL) data transmission.

We adopt the block-fading model, where $h_{kl}$∈${\mathbb{C}}^{N\times 1}$ represents the channel response between the k-th UE and the l-th AP, with $k=1,\dots ,K$ and $l=1,\dots ,L$. Within each block, an independent realization is drawn from an independent Rayleigh fading distribution:

$$h_{l,k}=N_c(0,{\beta }_{kl}{I}_N)$$

(2)

The large-scale fading coefficients ${\beta }_{l,k}$ incorporate pathloss and shadow fading as follows [7]:

$${\beta }_{l,k}=P{L}_{l,k}.{10}^{{\displaystyle \frac{{\sigma }_{sh}{z}_{l,k}}{10}}}$$

(3)

where $d_{l,k}$ is the distance between AP l and UE k, including the AP’s and UE’s heights. The shadow fading accounts for spatial correlation both between APs and between UEs.

We consider an area of radius R meters (1000 m). APs and UEs are uniformly distributed at random on the circle. For the simulation, the movement of the UEs is randomly positioned within the circle of radius R, and their motion is confined within it. A Brownian motion simulation is used for jumps to subsequent positions. The dwell time at each position depends on the velocity. The velocity of the UE is defined as a linear distribution within an interval $[v_{min},v_{max}]$ that is defined as desired. Figure 3 shows the random locations of APs and the paths followed by 40 UEs. All the results presented are formed over 100 network snapshots for different random realizations of the AP/UE locations.

The simulation scenario has been designed to identify the advantages and disadvantages of the proposed method compared to other methods commonly in wireless networks.

The main points to be evaluated in the proposed scenario are as follows.

• To evaluate the impact of using DMD on the accuracy and adaptability of CSI in CF-MMIMO systems, especially in highly dynamic scenarios.
• To evaluate the performance in terms of the accuracy and adaptability of DMD/PiDMD-based Channel Estimation compared to the most common methods used in wireless networks, such as MR, MMS, and other recently proposed methods, i.e., PZF, FZF, and PZFZF.
• To evaluate the overall system performance.

To achieve the stated objectives, the scenario used was formulated as follows.

• Randomly deploy the APs across the simulation area using an independent uniform distribution, employing a square grid within the circular simulation region.
• Sequentially introduce the UEs in a random fashion and simulate their movement using Brownian motion.
• Compute the distance between the chosen UE and each AP.
• Calculate the channel gain from the selected UE to each AP using Equation (3).
• Construct spatial correlation matrices $R_{kl}$ and estimation error correlation matrices $C_{kl}$.
• Generate the estimated channels $h_{kl}$, and utilize them to compute sample averages that approximate all the expectations in the SE expressions.
• Compute the CSI, uplink, and downlink capacity using the DMD algorithms.

2.2. DMD/PiDMD/mpEDMD for Channel Estimation in Wireless Systems

DMD/PiDMD/mpEDMD can be applied to channel estimates in wireless systems, including CF-MMIMO scenarios, to capture the dynamic behavior of the wireless channel, which highly improves the accuracy of CSI estimation. The following steps outline how DMD can be used for Channel Estimation following this procedure:

• Data Collection: Collect a series of snapshots of the received signals from the wireless channel over a period of time. These snapshots can be obtained through the transmission of known pilot symbols or training sequences.
• Data Matrices: Organize the received signal snapshots into two data matrices, often referred to as X and Y. These matrices represent the received signals at consecutive time steps.
• Singular Value Decomposition (SVD): Perform a singular value decomposition (SVD) on matrix X to factorize it into three matrices: U, $\Sigma$ (sigma), and $V^T$. The SVD captures the dominant modes of the received signals.
• DMD/PiDMD Modes and Eigenvalues: Calculate the DMD/PiDMD modes by using the relationship between matrices X and Y. These modes represent the coherent spatial patterns or oscillations in the channel behavior. The associated eigenvalues provide information about the growth rates and frequencies of these modes.
• Channel Estimation and Prediction: The DMD/PiDMD/mpEDMD modes and eigenvalues can be utilized for Channel Estimation and prediction. The DMD/PiDMD/mpEDMD modes help in identifying the spatial characteristics of the channel, while the eigenvalues offer insights into the temporal evolution of these characteristics.
• Adaptive Strategies: DMD-based Channel Estimation can be integrated into adaptive communication strategies. For instance, in MIMO systems, the identified modes can guide the selection of appropriate beamforming vectors, precoding matrices, or transmit power allocation, based on the changing channel conditions.
• Interpolation and Extrapolation: DMD/PiDMD/mpEDMD also enables interpolation and extrapolation of the channel behavior between and beyond the observed time instances. This can lead to more accurate CSI estimates during periods without pilot transmissions.
• Performance Evaluation: Validate the effectiveness of DMD-based Channel Estimation through simulations or real-world measurements. Compare the performance of DMD-based estimation with conventional methods to assess its benefits in terms of accuracy and adaptability.

It is important to note that while DMD-based methods have the potential to improve Channel Estimation in wireless systems, their success depends on factors such as the quality and quantity of data, the dynamic range of the channel variations, and the computational resources available for DMD calculations. In scenarios where the wireless channel experiences rapid variations, DMD and its variants can offer advantages by capturing the temporal evolution of the channel, leading to more accurate CSI estimates and enhanced system performance.

2.3. Impact of Channel Correlation

Small-scale fading arises from the presence of multiple propagation paths, resulting in signal replicas with varying phase shifts [23]. Within a local area, the variation in these phase shifts across different locations is influenced by the angle between the propagation paths and the line connecting those locations. This concept is illustrated in Figure 4, which depicts the line connecting two locations and a propagation path at an angle $\theta$ relative to it. The change in path length between the two locations is the projection of this line in the direction of the path [24]. Since this projection depends on $\varphi$, a key aspect of characterizing small-scale fading is a function that describes the angular distribution of paths, known as the Angular Spread Distribution (ASD). The ASD represents the received signal power as a function of the angle when acting as a receiver, or conversely, it represents the transmitted signal power reaching the receiver at the other end of the link as a function of the angle of departure. The ASD is a continuous function of angle, normalized to be interpreted as a probability density function (PDF), consistent with the unit-power normalization of small-scale fading. The standard deviation of the ASD, when viewed as a distribution, provides the root-mean-square (RMS) angle spread [24].

While uncorrelated Rayleigh fading is a common model due to its simplicity, it does not always accurately represent practical scenarios [24]. This paper proposes a Doppler–delay correlation model where $R_{\mathrm{delay}}$ is the delay correlation matrix and $R_{\mathrm{doppler}}$ is the Doppler correlation matrix. Incorporating delay and Doppler correlation into the model makes it more realistic, as these variables persist for longer periods, as demonstrated in Orthogonal Time–Frequency Space (OTFS) [25].

To begin, consider an uncorrelated scattering environment. The frequency–domain correlation $R_h(.)$ and the Power Delay Profile (PDP) can be related through a Fourier transform [24].

The distribution of received power as a function of delay, normalized to be interpreted as a PDF, is termed the Power Delay Profile (PDP) and is defined by $S_h\left(\tau \right)=E\left\{{\left|h\left(\tau \right)\right|}^2\right\}$ [26].

$$\begin{array}{cc}\hfill R_h\left({\Delta }_{\nu }\right)& =\mathbb{E}\left[h\left(\nu \right)h^{*}(\nu +{\Delta }_{\nu })\right]\hfill \end{array}$$

(4)

$$\begin{array}{cc}& =\mathbb{E}\left[{\int }_{-\infty }^{\infty }h\left({\tau }_0\right)e^{-j2\pi \nu {\tau }_0}d{\tau }_0{\int }_{-\infty }^{\infty }h\left({\tau }_1\right)e^{-j2\pi (\nu +{\Delta }_{\nu }){\tau }_1}d{\tau }_1\right]\hfill \end{array}$$

(5)

$$\begin{array}{cc}& ={\int }_{-\infty }^{\infty }{\mathbb{E}\left[\right|h\left(\tau \right)|}^2]e^{j2\pi {\Delta }_{\nu \tau }}d\tau .\hfill \end{array}$$

(6)

$$R_h\left({\Delta }_{\nu }\right)={\int }_{-\infty }^{\infty }\mathbb{E}\left[{\left|h\left(\tau \right)\right|}^2\right]e^{j2\pi {\Delta }_{\nu \tau }}d\tau ,$$

(7)

$$\begin{array}{cc}& ={\int }_{-\infty }^{\infty }{S}_h\left(\tau \right)e^{j2\pi {\Delta }_{\nu \tau }}d\tau .\hfill \end{array}$$

(8)

Assuming $\mathbb{E}\left[h\left({\tau }_0\right)h\left({\tau }_1\right)\right]=0$ for ${\tau }_0\ne {\tau }_1$, the condition outlined in (7) ensures that $\mathbb{E}\left[h\left({\tau }_0\right)h\left({\tau }_1\right)\right]=0$ for ${\tau }_0\ne {\tau }_1$. This condition, known as the uncorrelated scattering condition and specified in (8), is based on the idea that propagation mechanisms generating paths with different delays are inherently independent. When combined with the assumption of local stationarity, this condition forms the basis of traditional small-scale fading modeling within the wide-sense stationary uncorrelated scattering (WSSUS) framework [27]. Additionally, the uncorrelated scattering condition in the delay domain ensures that $R_h(.)$ depends only on $\Delta \nu$ rather than $\nu$, making the fading wide-sense stationary in the frequency domain as well. Due to the Fourier relationship between $R_h(.)$ and $S_h(.)$, it follows that $B_c\propto {\displaystyle \frac{1}{T_d}}$, with the exact proportionality determined by the characteristics of $R_h(.)$, $S_h(.)$, and the level of decorrelation represented by $B_c$ [24].

These conditions are often not met in practical scenarios [24]. For example, the assumption of independence between propagation mechanisms at ${\tau }_0$ and ${\tau }_1$ is frequently violated, despite the simplicity of this model, which has led to its widespread use. Introducing a more realistic model could significantly increase the complexity. This study investigates DMD-based approaches as a solution to this challenge and proposes a delay–Doppler correlation model.

The scattering function $S_h(\nu ,\tau )$ represents a combined Doppler–delay spectrum, which characterizes the distribution of power in the received signal concerning both Doppler frequency shift ($\nu$) and delay ($\tau$). This function is fundamental for understanding the signal’s spread across time and frequency resulting from multipath propagation. Mathematically, it is expressed as

$$S_h(\nu ,\tau )=\mathbb{E}\left[|\breve{h}{(\nu ,\tau )|}^2\right],$$

(9)

where $\breve{h}(\nu ,\tau )$ is the Doppler–delay spreading function representing the channel’s response at a given Doppler shift $\nu$ and delay $\tau$. The relationship to the time–frequency correlation function $R_h(\Delta t,\Delta f)$ is given by a two-dimensional Fourier transform [24]:

$$S_h(\nu ,\tau )={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{R}_h(\Delta t,\Delta f)e^{-j2\pi (\Delta t\nu -\Delta f\tau )}d\Delta td\Delta f.$$

(10)

This relationship allows the scattering function to encapsulate the statistical properties of the channel in both delay and Doppler domains. A more detailed proof is provided in Appendix A.

The Doppler correlation function, denoted as $R_{\nu }(\Delta \nu )$, quantifies the correlation between the channel’s response at different Doppler frequencies and is given by $R_{\nu }(\Delta \nu )=\mathbb{E}\left[h\left(\nu \right)h^{*}(\nu +\Delta \nu )\right]$, where $h\left(\nu \right)$ is the channel’s response at Doppler frequency $\nu$, and $h^{*}(\nu +\Delta \nu )$ is the complex conjugate of the channel’s response at $\nu +\Delta \nu$. The delay correlation function, often denoted as $R_h\left(\tau \right)$, quantifies the relationship between the channel’s impulse response at different delays and is defined as $R_h\left(\tau \right)=\mathbb{E}\left[h(t,\tau )h^{*}(t,\tau +\Delta \tau )\right]$, where $h(t,\tau )$ is the channel’s impulse response at time t with a delay $\tau$, and $h^{*}(t,\tau +\Delta \tau )$ is the complex conjugate of the channel’s response at a delay offset $\Delta \tau$.

2.4. Proposed Algorithm

• Uplink training phase: The received signal is mapped in the time–frequency plane; see Figure 5A. Then, the pilots are sampled. After taking k samples of pilots ${\varphi }_k$, R is computed as shown in (11), where each column of ${\Phi }_k$ matrix is given by ${\varphi }_n={\varphi }_1,n_1;{\varphi }_1,n_2;{\varphi }_1,n_3;\dots ;{\varphi }_1,n_k$, as shown in Figure 5B. Then, each column is used as shown in (12), and then DMD is applied to reduce the system dimension while preserving the main spatial and temporal characteristics. Subsequently, PiDMD and mpeDMD are applied. PiDMD is employed to ensure that the matrix takes the form of a Toeplitz matrix taking advantage of previous knowledge of the system. The mpeDMD technique is also capable of revealing the underlying Toeplitz structure, as will be demonstrated later.
  The corresponding approach to estimate the $M\times M$ correlation matrix R is to form the sample correlation matrix (11). As illustrated in Figure 5, the correlation matrix will be formed by the pilot’s samples (11):

  $${\widehat{R}}_{sample}={\displaystyle \frac{{\tau }_p}{N}}\sum _{s=1}^S{\varphi }_k{\varphi }_k^H.$$

  (11)
  The structure of $\varphi$ matrix is shown on (12). Each element of ${\widehat{R}}_{\mathrm{sample}}$ asymptotically converges to the corresponding element of R [11] by the law of large numbers. However, obtaining a sample correlation matrix with eigenvalues and eigenvectors that are well aligned with those of R is more challenging, as estimation errors in all $M^2$ elements of ${\widehat{R}}_{\mathrm{sample}}$ affect the eigenstructure. This is where PiDMD/mpeDMD comes into play. PiDMD/mpeDMD can accurately identify the eigenstructure of the data [9].
• Use PiDMD/mpeDMD: The purpose of this step is twofold; the first one is to keep only the most relevant components of the system and the second one is to achieve a better approximation of eigenvalues and eigenvectors of a reduced version of the system.

  $${\Phi }_k=\left[\begin{array}{cccc}|& |& |& |\\ {\varphi }_1& {\varphi }_2& \dots & {\varphi }_k\\ |& |& |& |\end{array}\right].$$

  (12)
• Channel Estimation: Once the correlation matrix is known in the whole network, it is possible to compute the Channel Estimation for every link.
• Design of Downlink precoders: With a good Channel Estimation, it is possible to use well-known downlink precoder design, locally at the Access Points or centrally at the CPU.
• Spectral Efficiency: Finally, the Spectral Efficiency can be computed and analyzed.

2.5. DMD/PiDMD/mpeDMD

2.5.1. Dynamic Mode Decomposition

DMD is a powerful data-driven technique that has found applications in various fields, ranging from fluid dynamics to neuroscience and beyond [28]. At its core, DMD seeks to capture the underlying dynamic behavior of a system from time-series data, offering insights into its modes of oscillation, growth rates, and spatial structures [8].

DMD operates on the premise that complex dynamics can be approximated by a set of coherent modes that evolve over time [8]. Given a sequence of data snapshots, DMD extracts these modes and associated frequencies to create a compact representation of the system’s behavior. This is particularly valuable in scenarios where the underlying dynamics are complex and difficult to model analytically.

DMD has demonstrated its efficacy in extracting relevant features from complex and high-dimensional data, making it particularly suitable for applications in which understanding system dynamics is crucial. In the context of CF-MMIMO wireless systems, DMD offers a novel approach to Channel Estimation by capturing the spatial–temporal evolution of the wireless channel’s behavior. By leveraging DMD, it becomes possible to adaptively track and predict the changing channel characteristics, thereby enhancing the efficiency and reliability of communication in dynamic environments [8]. The DMD algorithm is as follows.

$$y_{k+1}\approx A{y}_k,$$

(13)

$$\underset{A\in {\mathbb{R}}^{\mathbb{N}\times \mathbb{N}}}{min}\left|\right|y_{k+1}-A{y}_k{\left|\right|}_F.$$

(14)

which has the solution

$$A\approx y_{k+1}{y}_k^{†}=y_{k+1}W\Sigma V^{*}.$$

(15)

where † indicates the Moore–Penrose pseudoinverse. In (15),

$$W\in {\mathbb{R}}^{Nxr}V\in {\mathbb{R}}^{(m-1)xr},$$

(16)

$$y_k^{†}=V{\Sigma }^{-1}{W}^T.$$

(17)

$$\widehat{A}=y_{k+1}V{\Sigma }^{-1}{W}^T\in {\mathbb{R}}^{N\times N},$$

(18)

$$W^{T}AW=W{y}_{k+1}V{\Sigma }^{-1}{W}^{T}W\in {\mathbb{R}}^{N\times N},$$

(19)

$$\widehat{A}=U_r^{*}A{V}_r=W_r^T{y}_{k+1}{V}_r{\Sigma }_r^{-1}\in {\left(R\right)}^{r\times r}.$$

(20)

It is now feasible to compute the eigendecomposition of $\widehat{A}$. The eigenvectors of A can be approximated from the reduced eigenvectors $\Psi$ by (21):

$$\widehat{A}\Psi =\Psi \widehat{\Lambda }\Rightarrow \widehat{A}=X\Lambda X^{-1}.$$

(21)

If we start in state $u_0$, $u_k$ can be found using A, as shown next:

$$\begin{array}{cc}\hfill {\widehat{u}}_k& \approx \widehat{A}{\widehat{u}}_{k-1}.\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill {\widehat{u}}_k& \approx X\Lambda X^{-1}{u}_{k-1}.\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill X^{-1}{\widehat{u}}_k& \approx \Lambda X^{-1}{\widehat{u}}_{k-1}.\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill {\tilde{u}}_k& \approx \Lambda {\tilde{u}}_{k-1}.\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill & \approx \Lambda \Lambda {\tilde{u}}_{k-2}={\Lambda }^2{\tilde{u}}_{k-2}.\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill & \approx {\Lambda }^2\Lambda {\tilde{u}}_{k-3}={\Lambda }^3{\tilde{u}}_{k-3}.\hfill \\ \hfill & .\hfill \\ \hfill & .\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill & \approx {\Lambda }^k{\tilde{u}}_0.\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill {\tilde{u}}_k& \approx {\Lambda }^k{\tilde{u}}_0.\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill X^{-1}{\widehat{u}}_k& \approx {\Lambda }^k{X}^{-1}{\widehat{u}}_0.\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill {\widehat{u}}_k& \approx X{\Lambda }^k{X}^{-1}{\widehat{u}}_0.\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill W^T{\widehat{u}}_k& \approx X{\Lambda }^k{X}^{-1}{W}^T{\widehat{u}}_{0.}\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\widehat{u}}_k& \approx WX{\Lambda }^k{X}^{-1}{W}^T{\widehat{u}}_0.\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill u_k& =\Phi {\Lambda }^K{\Phi }^{†}{u}_0.\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill \varphi & =WX.\hfill \end{array}$$

(35)

Equation (37) is a method to compute b; alternative approaches to compute b can be found in [28,29,30].

$$\begin{array}{cc}\hfill b& ={\varphi }^{†}{u}_0.\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill u_k& =\varphi {\Lambda }^{K}b.\hfill \end{array}$$

(37)

where (37) predicts DMD at the $k_{th}$ time step [8].

This algorithm applied to correlation matrix R becomes

$${\Phi }_k=\left[\begin{array}{cccc}|& |& |& |\\ {{\varphi }_t}_1& {{\varphi }_t}_2& \dots & {{\varphi }_t}_{s-1}\\ |& |& |& |\end{array}\right].$$

(38)

$${\Phi }_{k+1}=\left[\begin{array}{cccc}|& |& |& |\\ {{\varphi }_t}_2& {{\varphi }_t}_3& \dots & {{\varphi }_t}_s\\ |& |& |& |\end{array}\right].$$

(39)

$$Y_k^{pilot}=H_k{\Phi }_k$$

(40)

$$R_k\approx E\left(Y_k^{pilot}{\left(Y_k^{pilot}\right)}^H\right).$$

(41)

$$R_{k+1}\approx E\left(Y_{k+1}^{pilot}{\left(Y_{k+1}^{pilot}\right)}^H\right).$$

(42)

$$R_{k+1}\approx A{R}_k.$$

(43)

$$A\approx R_{k+1}{R}_k^{†}=R_{k+1}V{\Sigma }^{†}{U}^{*}.$$

(44)

$$\widehat{A}=U_r^{*}A{V}_r=U_r^{*}{R}_{k+1}{V}_r{\Sigma }_r^{-1}.$$

(45)

$H_k$ values are the channel realizaton, ${\Phi }_k$ is the pilot matrix, and $R_k=U\Sigma V^{*}$ is the singular value decomposition. Once the singular value decomposition is computed, a reduced number of eigenvectors are picked to form $\widehat{A}$.

2.5.2. Physics-Informed Dynamic Mode Decomposition (piDMD)

Several variants of Dynamic Mode Decomposition (DMD) have been developed. When analyzing a system from a data-driven perspective, it is common to have partial knowledge of the system’s underlying physics. The initial step in physics-informed DMD (piDMD) requires the user to summarize the known or suspected physical properties of the system [9].

Once the physical principles to be enforced are determined, they must be translated into the matrix manifold to which the linear model should be constrained [9]. This step involves understanding the data collection procedure, such as the spatial grid or the energy inner product.

With a target matrix manifold defined, the relevant Procrustes problem [9] can be solved. Equation (46) is commonly known as the Procrustes problem in the literature. By appropriately modifying the optimization routine, it can be ensured that the resulting model adheres to the previously identified physical principles.

The final step of piDMD involves extracting physical information from the learned model [9]. This could include analyzing the spectrum or modes, computing the resolvent modes [31], making predictions, or exploring other diagnostics. When the matrix is Toeplitz, the algorithm takes the form of Equation (47). A Toeplitz matrix can be reduced to a circulant matrix [9].

$$\underset{R_{k+1}\in M}{arg\; min}\left|\right|R_{k+1}-A{R}_k{\left|\right|}_F,$$

(46)

where the A matrix is given by (44) or (45) if the reduced form is used. If matrix R is restricted to Toeplitz. then the problem can be written as (47). A key point to emphasize is that PiDMD leverages empirical data while simultaneously incorporating prior knowledge. Specifically, it enforces the condition that the matrix R exhibits a Toeplitz structure, thereby ensuring that this constraint is satisfied throughout the process.

$$\underset{R_{k+1} \mathrm{is} }{arg\; min}\left|\right|R_{k+L}-A{R}_k{\left|\right|}_F$$

(47)

$$R=\mathcal{F}diag\left(\widehat{r}\right){\mathcal{F}}^{-1}$$

(48)

Circulant matrices possess very useful properties, one of the most important being their ability to be diagonalized using the Fourier transform. In (48), ${\mathcal{F}}_{jk}=e^{2\pi i(j-2)(k-1)/n}/\sqrt{2}$ and ${\mathcal{F}}^{-1}={\mathcal{F}}^{*}$, with ${\widehat{r}}_j$ as the unknown eigenvalues and $R_{j,k}=r_{(j-k)modn}$. Once R is diagonalized, it can be written as

$$R=\left[\begin{array}{cccc}{r}_0& r_{-1}& \dots & r_{-(n-1)}\\ r_1& r_0& \ddots & ⋮\\ ⋮& \ddots & \ddots & r_{-1}\\ r_{n-1}& \dots & r_1& r_0\end{array}\right].$$

(49)

A Toeplitz matrix can be straightforwardly transformed into a circulant matrix form. Similarly to the circulant case, a Toeplitz matrix has elements that satisfy the relation $R_{i,j}=r_{i-j}$, with the distinction that in this case, the indices are not taken modulo n. As demonstrated in Equation (50), a Toeplitz matrix can be embedded within a circulant matrix. Once this transformation is achieved, the resulting circulant matrix can be diagonalized, as outlined above. Denoting the circulant matrix as C, this matrix is constructed using an auxiliary matrix, denoted $\tilde{R}$, whose entries are known a priori. The explicit forms of both matrices are illustrated in Equation (50).

$$\begin{array}{c}C=\left[\begin{array}{cc}R& \tilde{R}\\ \tilde{R}& R\end{array}\right]\\ \\ \\ \tilde{R}=\left[\begin{array}{cccc}0& r_{n-1}& \dots & r_1\\ r_{-(n-1)}& 0& \ddots & ⋮\\ ⋮& \ddots & \ddots & r_{-(n-1)}\\ r_{-1}& \dots & r_{-(n-1)}& 0\end{array}\right]\end{array}$$

(50)

Finally, if the matrix is circulant, the process for solving Equation (47) becomes significantly simplified, as detailed in [9]. The procedure is straightforward: first, a Toeplitz matrix is constructed, and then, for mathematical simplification, a circulant matrix is formed.

2.5.3. Measure-Preserving Extended Dynamic Mode Decomposition (mpeDMD)

This DMD variant is a powerful method that encodes geometric features, invariant measures, transient and long-time behavior, coherent structures, quasiperiodicity, etc. It is fully described in [10].

The basic requirement for the application of this method is having a measure-preserving system. An ergodic system falls into this category [10]. In this work, we are using a delay–Doppler channel model, which is in a relatively long period of time and is assumed to be ergodic, in which case this method is suitable for this scenario.

The Measure-Preserving Extended Dynamic Mode Decomposition (MPeDMD) is an extension of the DMD method used to analyze and extract dynamic modes from the sampled data of dynamical systems. mpeDMD is particularly useful for systems with ergodic characteristics and for preserving probability measures and ergodicity invariance properties. This method has the following steps:

• Data Collection: Dynamic data of a system are collected in the form of a snapshot matrix X, where each column represents a sample of the system state at a given time.
  Algorithm input: Snapshot data $\{x^{\left(m\right)},y^{\left(m\right)}\}=F\left(x^{\left(m\right)}\right)$, quadrature weights ${\left\{w_m\right\}}_{m=1}^M$, and a dictionary of functions ${\left\{\psi \right\}}_{j=1}^N$.

  $$\begin{array}{cc}\hfill {\Psi }_X& =\left(\begin{array}{c}\Phi \left(x^{\left(1\right)}\right)\\ ⋮\\ \Phi \left(x^{\left(M\right)}\right)\end{array}\right)\in {\mathbb{C}}^{M\times N},\hfill \end{array}$$

  (51)

  $$\begin{array}{cc}\hfill {\Psi }_Y& =\left(\begin{array}{c}\Phi \left(y^{\left(1\right)}\right)\\ ⋮\\ \Phi \left(y^{\left(M\right)}\right)\end{array}\right)\in {\mathbb{C}}^{M\times N}.\hfill \end{array}$$

  (52)
• Eigenvalue–Eigenvector Problem Formulation: The goal is to find dynamic modes $\varphi$ and eigenvalues $\lambda$ that satisfy the following equation:

  $$X^{\prime }\approx \Phi \Lambda .$$

  where $\Phi$ contains the dynamic modes and $\Lambda$ is a diagonal matrix of eigenvalues.
• Measure-Preserving Optimization: mpeDMD extends the standard DMD to preserve the system measure. This is achieved by optimizing the cost function to preserve the properties of the original system measure. The optimization problem is shown in (53)

  $${\mathbb{K}}_{EDMD}\in arg\underset{B\in {\mathbb{C}}^{N\times N}}{min}{\displaystyle \frac{1}{M}}\sum _{m=1}^M{∥\Psi \left(y^{\left(m\right)}\right)-\Psi \left(x^{\left(m\right)}\right)B∥}_2^2.$$

  (53)
  The Gram matrix is defined as $G={\Psi }_X^{*}W{\Psi }_X$ and the matrix $A={\Psi }_X^{*}W{\Psi }_Y$. Letting ^† denote the pseudoinverse, a solution to (53) is

  $${\mathrm{K}}_{EDMD}=G^{†}A={\left({\Psi }_X^{*}W{\Psi }_X\right)}^{†}\left({\Psi }_X^{*}W{\Psi }_Y\right).$$

  (54)

  $$={\left(\sqrt{W}{\Psi }_X\right)}^{†}\sqrt{W}{\Psi }_Y.$$

  (55)
• Selection of Relevant Modes: After calculating the dynamic modes and eigenvalues, the most relevant modes for system analysis are selected.
• Prediction and Analysis: Once the relevant modes are identified, they can be used to predict the future behavior of the system and to analyze its dynamic characteristics.

mpeDMD is especially useful for complex systems where ergodicity and measure-preservation properties are important, such as fluid systems, thermal systems, and quantum systems.

This method provides a powerful tool for the analysis of dynamical systems from sampled data and has found applications in a wide range of fields, including engineering, physics, biology, and economics.

The algorithm is as follows [10]:

• Compute $G={\Psi }^{*}XW\Psi X$ and $A={\Psi }^{*}XW\Psi Y$, where ${\Psi }_X,{\Psi }_Y$ are given in (51) and (52).
• Compute an SVD of $G^{-{\displaystyle \frac{1}{2}}}{A}^{*}{G}^{-{\displaystyle \frac{1}{2}}}=U_1\Sigma U_2^{*}$.
• Compute the eigendecomposition $U_2{U}_1^{*}=\widehat{V}\Lambda {\widehat{V}}^{*}$.
• Compute $K=G^{-{\displaystyle \frac{1}{2}}}{U}_2{U}_1^{*}{G}^{{\displaystyle \frac{1}{2}}}$ and $V=G^{-{\displaystyle \frac{1}{2}}}\widehat{V}$.

Output: Koopman matrix K, with eigenvectors V and eigenvalues $\Lambda$.

3. Channel Estimation

In this study, we analyze delay and Doppler-correlated Rayleigh fading channels, acknowledging that real-world channels inherently exhibit correlation [11]. The spatially correlated channel vector between the mth AP and the kth user, denoted as $g_{mk}\in {\mathbb{C}}^{N\times 1}$ (where $m=1,\dots ,M$ and $k=1,\dots ,K$), is described by (56), with ${\beta }_{mk}$ representing the large-scale fading coefficient.

$$g_{mk}={\beta }_{mk}^{1/2}{h}_{mk}.$$

(56)

CF-MMIMO systems exhibit significant macro-diversity and beneficial signal propagation [11], both of which are amplified by the deployment of a large number of Access Points (APs) [11]. Consequently, traditional Maximum Ratio (MR) processing [11] is commonly employed due to its optimal performance in environments with minimal inter-user interference. However, practical implementations often achieve substantial improvements in Spectral Efficiency (SE) through the use of Minimum Mean Square Error (MMSE) processing techniques, which effectively reduce inter-user interference [2,3]. The highest SE is attained through centralized processing, where the central processing unit (CPU) utilizes CSI to coordinate interference suppression among APs. Nevertheless, for practical reasons, the local execution of receiving, combining, and transmitting the precoding at each AP is preferred.

The following linear channel estimators are compared: Maximum Ratio (MR), Minimum Mean Squared Error (MMSE), Zero Forcing (ZF), and Partial Zero Forcing (PZF) channel estimators, in both traditional and DMD/PiDMD/mpeDMD fashions.

Figure 6 displays the norm of the eigenvalues of the Toeplitz structure of the matrix R considering the simulation scenario used in this work. In order to compute the Correlation matrix expression, (6) is used. This graph will be used in the next sections as a baseline for comparison.

Initially, it is crucial to determine whether there exists a relatively small group of eigenvalues that are significantly larger in magnitude than the others, which is essential for a proper DMD application. To demonstrate this empirically, several simulations are conducted using the scenario described in the Simulation Scenario section, employing a Monte Carlo algorithm with 500 scenarios. Figure 7 displays a typical graph illustrating the composition of normalized eigenvalues. Clearly visible in this graph is a small subset of eigenvalues that can characterize the system. This condition ensures the efficient use of the proposed algorithm.

3.1. Channel Estimation Using DMD

If the empirical correlation matrix formation is used [11], as specified in (11), a typical graph illustrating the norm of the magnitude of the matrix (with complex entries in the correlation matrix) can be observed in Figure 8. It is evident that the structure of the eigenvalues and eigenvectors lacks a defined pattern.

The application of DMD facilitates the extraction of the most representative spatio-temporal eigenvalues; however, starting from empirical data in the scenario proposed in this work, it does not yield a clearly defined matrix form. While the matrix R exhibits some improvement in structure using DMD, it remains a suboptimal approximation. The norm of the eigenvalues is illustrated in Figure 9. Although dimensionality has been reduced (arbitrarily set to 24 in this case), the larger eigenvalues are concentrated near the diagonal, and the overall configuration of matrix R remains significantly inadequate.

3.2. Channel Estimation Using PiDMD

If we know beforehand that the correlation structure follows a Toeplitz matrix form, we can employ PiDMD to generate a matrix with the desired structure using the empirical information available. Figure 10 illustrates the calculated correlation matrix under the specified simulation scenario.

The resulting Toeplitz matrix is not sufficiently satisfactory in terms of the eigenvalues obtained compared to the theoretical values. At this stage, following the proposed algorithm, DMD is initially used to reduce the order of the system, and subsequently, PiDMD values are applied to the resulting system. Figure 11 illustrates the outcome of the successive application of DMD and PiDMD.

3.3. Channel Estimation Using mpeDMD

As remarked in the previous section, the Toeplitz matrix obtained does not sufficiently match the theoretical eigenvalues from just using DMD. To address this, we first employ DMD to decrease the order of the system and then apply mpeDMD to the modified system. Figure 12 shows the results of sequentially using DMD and mpeDMD. mpeDMD is able to reuse the Toeplitz matrix shape of the empirical matrix R.

At this point, it is crucial to emphasize the importance and innovation of the proposed method: Figure 8 illustrates the empirical correlation matrix commonly suggested in the MIMO literature. The matrix $\mathbf{R}$ is fundamental for Channel Estimation; thus, starting from a vague approximation of $\mathbf{R}$ will inevitably lead to poor results. As demonstrated in this article, the proposed methods allow for the derivation of a much more accurate $\mathbf{R}$, leading to improved Channel Estimation and prediction. This results in a substantial enhancement in system performance through a highly efficient methodology capable of dealing with channels that closely resemble real-world conditions. Notably, this approach includes the treatment of correlation, which has been addressed in very few previous studies, and most importantly, it does so in a scalable and computationally feasible manner.

3.4. Channel Hardening and Favorable Propagation

The concepts of channel hardening and favorable propagation are crucial in MIMO systems [5]. We believe that a graphical representation of the behavior of the matrix $\mathbf{R}$ in the context of CF-MIMO is highly beneficial for a deeper understanding of these phenomena. In the Results section, the same experiments are qualitatively presented; however, this section provides a graphical intuition that illustrates these phenomena. Figure 13 and Figure 14 show the effect of the number of antennas and number of pilots, respectively, on $\mathbf{R}$. Channel hardening is a phenomenon in which the channel becomes more predictable as the number of antennas increases. On the other hand, favorable propagation leads to improved channel orthogonality, also due to the large number of antennas. Consequently, the number of pilots significantly impacts the formation of the $\mathbf{R}$ matrix.

3.5. Channel Prediction

DMD/PiDMD has demonstrated its efficacy in extracting intrinsic time and space modes from empirical measurements. In this section, we explore the effectiveness of DMD/PiDMD within the context of CF-MMIMO systems. Specifically, DMD methods prove to be an exceptionally powerful tool for channel prediction. This section will address the major challenges associated with channel prediction and discuss the utility of DMD/PiDMD/mpeDMD in overcoming these challenges.

3.6. Channel Aging

This section delves into the advanced methodologies for predicting channels in CF-MIMO systems. Conventional channel-prediction methods often rely on theoretical models, which might not fully encapsulate the intricacies of real-world scenarios. The relative movement between User Equipments (UEs) and Access Points (APs) leads to temporal fluctuations in the propagation environment, impacting channel coefficients within a single resource block. The channel realization $h_{kl}\left[n\right]$ is defined by its initial state $h_{kl}\left[0\right]$ and an innovation term [32]. Within this framework, we introduce and assess DMD/PiDMD-based predictors, benchmarking them against current techniques.

Historically, channel estimation and prediction in massive MIMO systems have employed methods such as Wiener or Kalman filtering, assuming model-based analytical channels [33,34,35]. For example, Zheng et al. [19] employed a simple first-order Gauss–Markov process for channel prediction, while more sophisticated autoregressive (AR) and autoregressive moving average (ARMA) models—linear stochastic models for correlated random processes—were explored in [35,36]. Despite their effectiveness, these approaches rely on simplified analytical models for long-term channel statistics, such as using a rectangular power spectrum to model temporal variations in channels, as discussed in [36], which may not fully encapsulate real-world conditions.

Previous prediction methods have also utilized autoregressive models derived from Yule–Walker equations [19], with additional research investigating Machine Learning (ML) approaches [37]. An ML-based predictor with noise pre-processing using linear Minimum Mean Square Error (LMMSE) has also been proposed.

Numerical analyses reveal that PiDMD/mpeDMD channel predictors offer a significant improvement in prediction accuracy and data rate compared to conventional channel models. While the ML-based predictor initially exhibits greater computational complexity than the DMD-based predictors, its operational complexity markedly decreases post-training, eventually surpassing the efficiency of the DMD/PiDMD-based methods.

The Kalman prediction model used for comparison is detailed in [35]. The correlation matrices employed for LMMSE Channel Estimation are empirically derived as follows:

$$\begin{array}{cc}\hfill C_{h_n{y}_n}& =\mathbb{E}\left[h_n{y}_n^H\right],\hfill \\ \hfill C_{y_n}& =\mathbb{E}\left[y_n{y}_n^H\right],\hfill \\ \hfill C_{h_n}& =\mathbb{E}\left[h_n{h}_n^H\right].\hfill \end{array}$$

The base station acquires the sampled auto-covariance $C_{y_n}$ through the following expression:

$$C_{y_n}={\displaystyle \frac{1}{N_s}}\sum _{i=1}^{N_s}{y}_i{y}_i^H$$

(57)

The LMMSE estimate $g_n$ of $y_n$ is computed as

$$\begin{array}{cc}\hfill g_n& =C_{h_n{y}_n}{C}_{y_n}^{-1}{y}_n\hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill & =C_{h_n}{\Psi }_n^H{({\Psi }_n{C}_{h_n}{\Psi }_n^H+I)}^{-1}{y}_n.\hfill \end{array}$$

(59)

Here, ${\Psi }_n$ is defined as

$${\Psi }_{li}^j={\left(\sum _{l,i\in \mathit{P}}{p}_{li}{\tau }_p{\mathbb{R}}_{li}^j+{\sigma }_{UL}^2{\mathbb{I}}_M^j\right)}^{-1}.$$

(60)

Following the algorithm proposed in this study, after obtaining empirical covariance matrices, PiDMD/mpeDMD is applied. The accuracy of the estimate diminishes as the time between pilot transmission and the considered channel realization increases. In DMD channel prediction, once matrix A is computed as $R_{k+1}\approx A{R}_k$, the prediction for $R_{k+2}$ becomes $R_{k+2}\approx A{R}_{k+1}$. Matrix R encapsulates the primary eigenstructure of the system and is reduced in order, thereby decreasing computational complexity. Numerical results highlight the superior performance of DMD/PiDMD compared to Kalman and Machine Learning (MLP) prediction models. The MLP model described can be found in [35].

4. Uplink Transmission

4.1. Uplink Training and Channel Estimation

The uplink training methodology employed in this study follows the approach detailed in [7], wherein each User Equipment (UE) simultaneously transmits its pilot sequence to the Access Points (APs) once per coherence block. It is assumed that there are ${\tau }_p\le K$ mutually orthogonal pilot signals, each of length ${\tau }_p$, allocated to the UEs. The pilot sequence ${\varphi }_{i_k}\in {\mathbb{C}}^{{\tau }_p\times 1}$ is transmitted by UE k, where $i_k$ denotes the index of the pilot used by UE k. The set $P_k\subseteq \{1,\dots ,K\}$ includes all UEs, including UE k, that utilize the same pilot as UE k, ensuring orthogonality among all pilot sequences.

The received signal $Z_l\in {\mathbb{C}}^{N\times {\tau }_p}$ at AP l is expressed as

$$Z_l=\sum _{k=1}^K\sqrt{p_k}{h}_k{\varphi }_{i_k}^H+N_l,$$

(61)

where $N_l\in {\mathbb{C}}^{N\times {\tau }_p}$ is a Gaussian noise matrix with i.i.d. $\mathcal{NC}(0,{\sigma }^2)$ elements, and $p_k$ represents the transmit power of UE k for uplink training. The MMSE estimate of $h_{kl}$ given $Z_l$ is provided in Equation (62) [2]:

$${\widehat{h}}_{kl}\triangleq {\displaystyle \frac{c_{kl}}{\sqrt{{\tau }_p}}}{Z}_l{\varphi }_{i_k},$$

(62)

where $c_{kl}$ is defined as

$$c_{kl}\triangleq {\displaystyle \frac{p_k{\tau }_p{\beta }_{kl}}{{\tau }_p{\sum }_{t\in P_k}{p}_t{\beta }_{kl}+{\sigma }^2}}.$$

(63)

It can be verified that the estimate ${\widehat{h}}_{kl}$ and the estimation error ${\tilde{h}}_{kl}=h_{kl}-{\widehat{h}}_{kl}$ are independent Gaussian vectors with distributions ${\widehat{h}}_{kl}\sim \mathcal{NC}(0,{\gamma }_{kl}{I}_N)$ and ${\tilde{h}}_{kl}\sim \mathcal{NC}(0,({\beta }_{kl}-{\gamma }_{kl})I_N)$, respectively, where

$${\gamma }_{kl}\triangleq \mathbb{E}\left[{∥{\widehat{h}}_{kl}∥}^2\right]={\displaystyle \frac{p_k{\tau }_p{\beta }_{kl}^2}{{\tau }_p{\sum }_{t\in P_k}{p}_t{\beta }_{tl}+{\sigma }^2}}.$$

(64)

The channel estimate ${\overline{H}}_l=Z_l\Phi \in {\mathbb{C}}^{N\times {\tau }_p}$ is utilized, where $\Phi =\left[{\varphi }_1,\dots ,{\varphi }_{{\tau }_p}\right]\in {\mathbb{C}}^{{\tau }_p\times {\tau }_p}$ as defined in Equation (38). By applying DMD, ${\overline{H}}_l$ is processed to estimate ${\widehat{h}}_{kl}$ as

$${\widehat{h}}_{kl}=c_{kl}{\overline{H}}_l{e}_{i_k},$$

(65)

where $e_{i_k}$ denotes the $i_k$-th column of the identity matrix $I_{{\tau }_p}$.

The second layer of centralized decoding involves processing the local estimates $\{{\widehat{s}}_{kl}:l=1,\dots ,L\}$ along with LSFD coefficients $\{a_{kl}:l=1,\dots ,L\}$ to compute ${\widehat{s}}_k={\sum }_{l=1}^L{a}_{kl}^{*}{\widehat{s}}_{kl}$. In this study, DMD is employed to enhance efficiency by capturing and utilizing the reduced eigenstructures of consecutive pilot snapshots.

4.2. Uplink Data Transmission

The uplink received signal at each Access Point (AP) is expressed as

$$\begin{array}{cc}\hfill s_k& =\sum _{l=1}^L{\widehat{s}}_{kl}=\sum _{l=1}^L{\mathbf{v}}_k^H{D}_k{\mathbf{y}}_{ul}^l.\hfill \end{array}$$

(66)

where $s_k\sim \mathcal{N}(0,p_k^{ul})$ denotes the information-bearing signal transmitted by User Equipment (UE) k with power $p_k^{ul}$, and ${\mathbf{N}}_l\sim \mathcal{N}(0,{\sigma }^2)$ represents the independent receiver noise. The matrix $D_k$ is employed to select the UEs served by each AP k.

Based on the signal in (66), the APs and the central processing unit (CPU) can decode the symbols using one of two methods: centralized or distributed. In the distributed approach, each AP independently decodes the information. Conversely, the centralized method gathers information from the APs, which is then processed at the CPU. This study employs the Least-Squares Frequency-Domain (LSFD) technique as the centralized method. The core idea of LSFD is that each AP computes local estimates of the desired data for all UEs in the first layer and forwards them to the CPU for final decoding in the second layer. Specifically, an estimate of the data symbol from UE k at AP l is determined by local linear combining using the vector ${\mathbf{v}}_{kl}\in {\mathbb{C}}^{N\times 1}$ in the first layer.

To detect the symbol transmitted by the k-th UE, the l-th AP calculates the product of the received signal ${\mathbf{y}}_l\left[n\right]$ and the conjugate of its locally estimated channel. The resulting value, ${\widehat{s}}_{kl}={\widehat{h}}_{kl}^H{\mathbf{y}}_l\left[n\right]$, is then sent to the CPU via the fronthaul. The CPU utilizes the coefficients $a_{kl}\left[n\right]$ to compute ${\widehat{s}}_k\left[n\right]$ [19] as follows:

$${\widehat{s}}_{kl}={\mathbf{v}}_{kl}^H{\mathbf{y}}_l={\mathbf{v}}_{kl}^H{\mathbf{h}}_{kl}{s}_k^H+{\mathbf{v}}_{kl}^H\sum _{t\in P_k}^K{\mathbf{h}}_{tl}{s}_t+{\mathbf{v}}_{kl}^H\sum _{t\notin P_k}^K{\mathbf{h}}_{tl}{s}_t+{\mathbf{v}}_{kl}^H{\mathbf{n}}_l.$$

(67)

The final estimate ${\widehat{s}}_k$ is computed as

$$\begin{array}{cc}\hfill {\widehat{s}}_k=& \sum _{l=1}^L{a}_{kl}^{*}{\mathbf{v}}_{kl}^H{\mathbf{h}}_{kl}{s}_k+\sum _{l=1}^L{a}_{kl}^{*}{\mathbf{v}}_{kl}^H\sum _{t\in P_k\setminus \left\{k\right\}}{\mathbf{h}}_{tl}{s}_t\hfill \\ \hfill & +\sum _{l=1}^L{a}_{kl}^{*}{\mathbf{v}}_{kl}^H\sum _{t\notin P_k}{\mathbf{h}}_{tl}{s}_t+\sum _{l=1}^L{a}_{kl}^{*}{\mathbf{v}}_{kl}^H{\mathbf{n}}_l.\hfill \end{array}$$

(68)

$$\begin{array}{c}\hfill {\widehat{s}}_k=\underset{\mathrm{A}}{\underbrace{{\mathbf{v}}_k^H{D}_k{\mathbf{b}}_{hk}{s}_k}}+\underset{\mathrm{B}}{\underbrace{{\mathbf{v}}_k^H{D}_k{\mathbf{e}}_{hk}{s}_k}}+\underset{\mathrm{I}}{\underbrace{\sum _{\begin{array}{c}i=1\\ i\ne k\end{array}}^K{\mathbf{v}}_k^H{D}_k{\mathbf{h}}_i{s}_i}}+\underset{\mathrm{N}}{\underbrace{{\mathbf{v}}_k^H{D}_k\mathbf{n}}},\end{array}$$

(69)

where A is the desired signal over the estimated channel, B is the desired signal over the unknown channel, I is interference, and N is noise.

The channel estimate ${\overline{H}}_l=Z_l\Phi \in {\mathbb{C}}^{N\times {\tau }_p}$ is then utilized, where $\Phi =\left[{\varphi }_1,\dots ,{\varphi }_{{\tau }_p}\right]\in {\mathbb{C}}^{{\tau }_p\times {\tau }_p}$ is defined in (38). Thus, the channel estimate ${\widehat{h}}_{kl}$ can be expressed as

$${\widehat{h}}_{kl}=c_{kl}{\overline{H}}_l{e}_{i_k},$$

(70)

where $e_{i_k}$ denotes the $i_k$-th column of the identity matrix $I_{{\tau }_p}$.

After estimating the local data, the second layer of centralized decoding is carried out using the local estimates $\{{\widehat{s}}_{kl}:l=1,\dots ,L\}$ and the LSFD coefficients $\{a_{kl}:l=1,\dots ,L\}$ to compute ${\widehat{s}}_k={\sum }_{l=1}^L{a}_{kl}^{*}{\widehat{s}}_{kl}$.

Dynamic Mode Decomposition (DMD) is employed to process S consecutive snapshots of the pilots, sequentially storing a reduced eigenstructure to estimate the channel, forming the matrix $\Phi$. This approach proves to be more efficient than the Least Mean Square Error (LMSE) method commonly used for Channel Estimation. The centralized algorithm sends the local channel estimates to the CPU to compute the optimal values of $a_{kl}$ for calculating ${\widehat{s}}_k\left[n\right]$, as shown in (71):

$$\begin{array}{c}\hfill {\widehat{s}}_{kl}=a_{kl}^{*}{\mathbf{h}}_{kl}^{DMD}{s}_k^H+\sum _{t\in P_k}^K{a}_{kl}^{*}{\mathbf{h}}_{tl}^{DMD}{s}_t+\sum _{t\notin P_k}^K{a}_{kl}^{*}{\mathbf{h}}_{tl}^{DMD}{s}_t+a_{kl}^{*}{\mathbf{n}}_l.\end{array}$$

(71)

4.3. Uplink Spectral Efficiency

According to the channel model in [23], the signal components are categorized as the desired signal (DS), beamforming uncertainty gain (BU), multiuser interference (UI), and noise interference (NI). These components are defined as follows:

$$\begin{array}{cc}\hfill D{S}_k& =\sqrt{\rho \eta }\mathbb{E}\left\{\sum _{m=1}^M{\widehat{g}}_{mk}^H{g}_{mk}\right\},\hfill \end{array}$$

(72)

$$\begin{array}{cc}\hfill B{U}_k& =\sqrt{\rho \eta }\left(\sum _{m=1}^M{\widehat{g}}_{mk}^H{g}_{mk}-\mathbb{E}\left\{\sum _{m=1}^M{\widehat{g}}_{mk}^H{g}_{mk}\right\}\right),\hfill \end{array}$$

(73)

$$\begin{array}{cc}\hfill U{I}_{k{k}^{\prime }}& =\sqrt{\rho }\mathbb{E}\left\{\sum _{m=1}^M\sqrt{{\eta }_k}{\widehat{g}}_{mk}^H{g}_{mk}\right\},\hfill \end{array}$$

(74)

$$\begin{array}{cc}\hfill N{I}_k& =\sum _{m=1}^M{\widehat{g}}_{mk}^H{w}_{u,m},\hfill \end{array}$$

(75)

where each expression corresponds to the terms in (77).

Invoking the arguments from [5], the achievable uplink Spectral Efficiency (SE) for $U{E}_k$ is given by (76), which can be written with the effective SINR, as shown in (78).

$$S{E}_k=\left(1-{\displaystyle \frac{{\tau }_p}{{\tau }_c}}\right){log}_2(1+\mathrm{SIN}\mathrm{R}),$$

(76)

$${\mathrm{SIN}\mathrm{R}}_k^{ul}={\displaystyle \frac{|D{S}_k{|}^2}{\mathbb{E}\left\{|B{U}_k{|}^2\right\}+\mathbb{E}\left\{|U{I}_{k{k}^{\prime }}{|}^2\right\}+\mathbb{E}\left\{|N{I}_k{|}^2\right\}}},$$

(77)

$$\begin{array}{c}\hfill {\mathrm{SIN}\mathrm{R}}_k^{ul}={\displaystyle \frac{p_k{\left|{\mathbf{a}}_k^H\mathbb{E}\left\{{\mathbf{g}}_{kk}\right\}\right|}^2}{{\mathbf{a}}_k^H\left({\sum }_{i=1}^K{p}_i\mathbb{E}\left\{{\mathbf{g}}_{ki}{\mathbf{g}}_{ki}^H\right\}-p_k\mathbb{E}\left\{{\mathbf{g}}_{kk}\right\}\mathbb{E}\left\{{\mathbf{g}}_{kk}^H\right\}+{\mathbf{F}}_k\right){\mathbf{a}}_k}}.\end{array}$$

(78)

where ${\mathbf{g}}_{ki}={[v_k^H{h}_{i1},\dots ,v_k^H{h}_{iL}]}^T$, $F_k=\mathrm{diag}\left(\right||v_{kL}{\left|\right|}^2,\dots ,\left|\right|v_{kL}{\left|\right|}^2)$, and ${\mathbf{a}}_k={[a_{k1},\dots ,a_{kL}]}^T$.

The effective SINR in (78) can be further maximized as shown in [7] by setting

$${\mathbf{a}}_k={\left(\sum _{t=1}^K{p}_{\mathrm{ul},t}{\mathbf{g}}_{kt}{\mathbf{g}}_{kt}^H+{\sigma }^2{\mathbf{F}}_k^{-1}\right)}^{-1}\mathbb{E}\left\{{\mathbf{g}}_{kk}\right\},$$

which leads to the maximum value [38]:

$${\mathrm{SIN}\mathrm{R}}_{\max ,k}={\displaystyle \frac{p_{\mathrm{ul},k}{\left|{\mathbf{g}}_{kk}^H\right|}^2}{{\sum }_{t=1}^K{p}_{\mathrm{ul},i}\mathbb{E}\left\{{\mathbf{g}}_{kt}{\mathbf{g}}_{kt}^H\right\}+{\sigma }^2{\mathbf{F}}_k-p_t^{\mathrm{ul}}\mathbb{E}\left\{{\mathbf{g}}_{kt}\right\}}}\mathbb{E}\left\{{\mathbf{g}}_{kk}\right\}.$$

(79)

4.4. Maximum Ratio Combining

The Maximum Ratio (MR) combining technique is the most straightforward method for combining received signals. It is defined as

$${\mathbf{v}}_{kl}^{\mathrm{MR}}={\mathbf{D}}_{kl}{\widehat{\mathbf{h}}}_{kl},$$

(80)

where ${\widehat{\mathbf{h}}}_{kl}=c_{kl}{\overline{\mathbf{H}}}_l{\mathbf{e}}_{i_k}$ represents the MMSE estimate, which is computationally efficient. This technique enhances the desired signal power in the numerator of Equation (78). The uplink Spectral Efficiency (SE) achievable with MR combining has been explicitly derived in [39].

In this research, MR combining [39] is utilized with ${\mathbf{v}}_{\mathrm{MR},kl}={\mathbf{D}}_{kl}{\mathbf{h}}_{kl}$. As a result, the expectations in Equation (78) are expressed as follows [4]:

$$\begin{array}{c}\hfill \mathbb{E}{\left\{g_{ki}\right\}}_l=\left\{\begin{array}{cc}\sqrt{{\eta }_k{\eta }_i{\tau }_p}{\mathbf{D}}_{kl}{\mathbf{R}}_{il}{\mathbf{R}}_{kl}{\Psi }_{t_{k}l}^{-1}{\mathbf{R}}_{kl},\hfill & \mathrm{if}i\in P_k\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.\end{array}$$

(81)

where ${\left[\mathbb{E}\left\{g_{ki}{g}_{ki}^H\right\}\right]}_{lr}={\left[\mathbb{E}\left\{g_{ki}\right\}\right]}_l{\left[\mathbb{E}\left\{g_{ki}\right\}\right]}_r^H$ for $r\ne l$, and

$${\mathbf{v}}_{\mathrm{MR},kl}={\mathbf{h}}_{kl}^{\mathrm{DMD}}.$$

(82)

Moreover,

$$\begin{array}{cc}\hfill \mathbb{E}\left\{g_{ki}{g}_{ki}^H\right\}& ={\eta }_k{\tau }_p\mathrm{tr}\left({\mathbf{D}}_{kl}{\mathbf{R}}_{il}{\mathbf{R}}_{kl}{\Psi }_{t_{k}l}^{-1}{\mathbf{R}}_{kl}\right)\hfill \\ \hfill & +\left\{\begin{array}{cc}{\eta }_k{\eta }_i{\tau }_p^2{\left|\mathrm{tr}\left({\mathbf{D}}_{kl}{\mathbf{R}}_{il}{\Psi }_{t_k}^{-1}{\mathbf{R}}_{kl}\right)\right|}^2,\hfill & \mathrm{if}i\in {\mathcal{P}}_k\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\hfill \end{array}$$

(83)

Additionally,

$${\left[F_k\right]}_{ll}={\sigma }_{ul}^2{\eta }_k{\tau }_p\mathrm{tr}\left({\mathbf{D}}_{kl}{\mathbf{R}}_{kl}{\Psi }_{t_{k}l}^{-1}{\mathbf{R}}_{kl}\right),$$

(84)

and

$${\Psi }_{li}^j={\left(\sum _{l,i\in \mathcal{P}}{p}_{li}{\tau }_p{\mathbb{R}}_{li}^j+{\sigma }_{UL}^2{\mathbb{I}}_M^j\right)}^{-1}.$$

(85)

When Dynamic Mode Decomposition (DMD) is applied, Equation (85) is adjusted to

$${{\Psi }_{li}^j}^{\mathrm{DMD}}={\left(\sum _{l,i\in \mathcal{P}}{p}_{li}{\tau }_p{\left({\mathbb{R}}_{li}^j\right)}^{\mathrm{DMD}}+{\sigma }_{UL}^2{\mathbb{I}}_M^j\right)}^{-1}.$$

(86)

4.5. Minimum Mean Squared Error Channel Estimator

According to the standard results in [3], the MMSE estimate of $h_{kl}$ based on the observed data is given by

$${\widehat{h}}_{lk}^j=\sqrt{p_{li}}{R}_{li}^j{\Psi }_{li}^j\left(y_{jli}^p\right).$$

(87)

The application of Dynamic Mode Decomposition (DMD) is illustrated in Equation (88).

$${\Psi }_{li}^j={\left(\sum _{l,i\in \mathcal{P}}{p}_{li}{\tau }_p{\left({\mathbb{R}}_{li}^j\right)}_L^{\mathrm{DMD}}+{\sigma }_{UL}^2{\left({{\mathbb{I}}_M}_j\right)}_L\right)}^{-1}.$$

(88)

Here, L represents the order of model reduction in the DMD model employed.

The estimation error [5], denoted as ${\tilde{h}}_{li}^j=h_{li}^j-\widehat{h_{li}^j}$, has a correlation matrix $\mathbb{C}=\mathbf{E}\{{\tilde{h}}_{li}^j{({\tilde{h}}_{li}^j)}^H\}$, which is expressed as follows [23]:

$${\mathbb{C}}_{li}^j={\mathbb{R}}_{li}^j-p_{li}{\tau }_p{\mathbb{R}}_{li}^j{\Psi }_{li}^j{\mathbb{R}}_{li}^j.$$

(89)

4.6. Full Zero Forcing

Several investigations have explored Zero Forcing (ZF) processing in Cloud-Fog Massive MIMO (CF-MMIMO) systems, as outlined in [40]. Centralized ZF methods, as described in [40,41], involve transmitting real-time Channel State Information (CSI) from all Access Points (APs) to a central processing unit (CPU) for designing ZF precoding or combining vectors. However, this centralized approach can result in substantial fronthaul traffic, making it less feasible as the number of User Equipments (UEs) increases. To mitigate this challenge, our study investigates distributed schemes where APs manage the processing independently.

Various distributed ZF precoding schemes, such as full-pilot ZF (FZF), partial FZF (PFZF), and protective PFZF (PPFZF), are presented in [7]. These approaches either fully distribute interference suppression or implement it in a coordinated and scalable manner, showing improved performance compared to the Maximum Ratio (MR) scheme. Additionally, local regularized ZF (LRZF) is evaluated as an upper benchmark in [7], despite its limitations in scalability.

The performance of FZF and an adapted LRZF precoder is further analyzed in [42] within NOMA-aided CF-MMIMO systems. However, the uplink scenario is not addressed in [42,43].

In the following sections, we focus on FZF and PFZF as distributed algorithms. For FZF, the local combining vector chosen by an AP for UE k, denoted as $v_{i,l}^{\mathrm{FZF}}\in {\mathbf{C}}^{N\times 1}$, is given by

$$v_{i,l}^{\mathrm{FZF}}=c_{i_k}{\theta }_{i_k}{\overline{H}}_l{\left(H_l^H{\overline{H}}_l\right)}^{-1},$$

(90)

where ${\theta }_{ikl}=E\left\{{\left|{\overline{H}}_{leik}\right|}^2\right\}={\displaystyle \frac{{\gamma }_{ikl}}{c_{ikl}^2}}$.

In our research, we apply Equation (91):

$$v_{i,l}^{\mathrm{FZF}\_\mathrm{DMD}}=c_{i_k}{\theta }_{i_k}{\overline{H}}^{\mathrm{DMD}}{\left(H_l^{\mathrm{DMD}}\right)}^H{\left({\overline{H}}_l^{\mathrm{DMD}}\right)}^{-1},$$

(91)

where $H^{\mathrm{DMD}}$ is determined using the proposed algorithms: DMD/PiDMD and DMD/ mpeDMD.

Matrix inversion is greatly simplified using Equation (92):

$$\begin{array}{cc}\hfill \tilde{A}& ={\tilde{U}}^{*}A\tilde{U}={\tilde{U}}^{*}{X}_0\tilde{V}{\tilde{\Sigma }}^{-1},\hfill \\ \hfill {\left(\widehat{A}\right)}^{-1}& ={\left({\widehat{U}}^{*}{X}_0\widehat{V}{\widehat{\Sigma }}^{-1}\widehat{U}\right)}^{-1}={\widehat{U}}^{*}{X}_0^{-1}{\widehat{V}}^{-1}\widehat{\Sigma }\widehat{U}.\hfill \end{array}$$

(92)

The inversion of matrix A is considerably more efficient with DMD due to the dimensionality reduction of matrix R and the use of Singular Value Decomposition (SVD), which simplifies the process. Channel prediction, taking these factors into account, is described by Equation (93):

$$H_{k+1}^{-1}=H_k^{-1}{A}^{-1}.$$

(93)

4.7. Partial Full Zero Forcing

The Full-Pilot Zero Forcing (FZF) combining method allocates ${\tau }_p$ degrees of freedom to address interference. However, since interference affecting User Equipment (UE) k primarily comes from a small subset of other UEs, we implement Partial-Pilot Zero Forcing (PFZF) combinations. This technique focuses on mitigating interference from UEs with strong channel gains while allowing for some interference from UEs with weaker gains. Consequently, Access Point (AP) l uses PFZF combining for strong UEs and Maximum Ratio (MR) combining for weaker ones.

The PFZF combining vector for UE $i_k$ at AP l, denoted as $v_{i_k,l}^{\mathrm{PFZF}}$, is given by

$$v_{i_k,l}^{\mathrm{PFZF}}=c_{i_k}{\theta }_{i_k}{\overline{H}}_l{E}_{S_l}{\left(E_{S_l}^H{\overline{H}}_l^H{\overline{H}}_l{E}_{S_l}\right)}^{-1},$$

(94)

where $E_{S_l}=(e_{r_{l,1}},\dots ,e_{r_{l,\tau S_l}})\in {\mathbb{C}}^{{\tau }_p\times \tau S_l}$ and $e_{r_{l,t}}$ denotes the t-th column of ${\mathbf{I}}_{\tau S_l}$.

Equation (94) is based on the derivation in [43]. In our study, we employ a similar formulation derived directly from the approach outlined in [43].

4.8. Spectral Efficiency and Optimal LSFD Weights

The received signal for User Equipment (UE) $i_k$ is described by [4]:

$${\widehat{s}}_{kl}={\mathbf{v}}_{kl}^H{\mathbf{D}}_k{\mathbf{y}}_{\mathrm{ul},l},$$

(95)

where ${\mathbf{D}}_k$ selects the signals meant for UE k.

By substituting (95) into (77) [4], the estimated signal ${\widehat{s}}_k$ for UE k in a distributed setting becomes

$${\widehat{s}}_k=\sum _{l=1}^L{a}_{kl}{\mathbf{v}}_{kl}^H{\mathbf{D}}_{kl}{\mathbf{h}}_{kl}{s}_k+\sum _{i=1}^K\sum _{\begin{array}{c}l=1\\ l\ne k\end{array}}^L{a}_{kl}{\mathbf{v}}_{kl}^H{\mathbf{D}}_{kl}{\mathbf{h}}_{il}{s}_i+n_{0k},$$

(96)

where $n_{0k}={\sum }_{l=1}^L{a}_{kl}{\mathbf{v}}_{kl}^H{\mathbf{D}}_{kl}{n}_l$ represents the noise term. The vector ${\mathbf{g}}_{ki}\in {\mathbb{C}}^L$ is defined as

$${\mathbf{g}}_{ki}=\left[\begin{array}{c}{\mathbf{v}}_{k1}^H{\mathbf{D}}_{k1}{\mathbf{h}}_{i1}\\ ⋮\\ {\mathbf{v}}_{kL}^H{\mathbf{D}}_{kL}{\mathbf{h}}_{iL}\end{array}\right].$$

This vector captures the combined receive channels between UE i and all APs serving UE k. In terms of this notation, (96) can be expressed as

$${\widehat{s}}_k=a_k^H{\mathbf{g}}_{kk}{s}_k+\sum _{i=1}^K\sum _{\begin{array}{c}l=1\\ l\ne k\end{array}}^L{a}_{kl}{\mathbf{g}}_{ki}{s}_i.$$

To compute the achievable Spectral Efficiency (SE) for UE k in a distributed context, we substitute (96) into the SE formula [4], leading to

$$\mathrm{SE}(u_{lk},d)={\displaystyle \frac{{\tau }_u}{{\tau }_c}}{log}_2\left(1+\mathrm{SIN}\mathrm{R}(u_{lk},d)\right)\mathrm{bit}/\mathrm{s}/\mathrm{Hz},$$

where the effective Signal-to-Interference-plus-Noise Ratio (SINR) is given by

$$\mathrm{SIN}\mathrm{R}(u_{lk},d)={\displaystyle \frac{p_k{\left|a_k^{H}E\left\{g_{kk}\right\}\right|}^2}{a_k^H\left({\sum }_{i=1}^K{p}_{i}E\left\{g_{ki}{g}_{ki}^H\right\}-p_{k}E\left\{g_{kk}\right\}E\left\{g_{kk}^H\right\}+F_k\right)a_k}},$$

(97)

with $F_k={\sigma }_{ul}^2\mathrm{diag}\left[E{\left\{n_k{D}_{k1}{v}_{k1}^k\right\}}^2,\dots ,E{\left\{n_k{D}_{kL}{v}_{kL}^k\right\}}^2\right]\in {\mathbb{R}}^{L\times L}$.

The combining vector ${\mathbf{v}}_{kl}$ that minimizes (97) is

$${\mathbf{v}}_{kl}^{\mathrm{L}\text{-}\mathrm{MMSE}}={\displaystyle \frac{p_k}{{\sum }_{i=1}^K{p}_i}}\left({\mathbf{b}}_{hl}{\mathbf{b}}_{il}^H+{\mathbf{C}}_{il}\right)+{\sigma }_{ul}^2{\mathbf{IN}}^{-1}{\mathbf{D}}_{kl}{\mathbf{b}}_{kl}.$$

The optimal effective SINR for UE k is achieved by

$$a_k^{\mathrm{opt}}={\displaystyle \frac{p_k}{{\sum }_{i=1}^K{p}_i}}E\left\{g_{ki}{g}_{ki}^H\right\}+F_k+{\sigma }_{Dk}^2{\mathbf{e}}^H{\left(E\left\{g_{kk}\right\}\right)}^{-1},$$

where ${\mathbf{e}}_{Dk}\in {\mathbb{R}}^{L\times L}$ is a diagonal matrix with the $(l,l)$-th element being one if $l\notin {\mathcal{M}}_k$ and zero otherwise.

5. Downlink Transmission

5.1. Coherent Downlink Transmission

Each downlink resource block allocates ${\tau }_c-{\tau }_p$ time slots for downlink data transmission. In the CF-MIMO system, coherent joint transmission is implemented, where every Access Point (AP) concurrently transmits the same data symbol to all User Equipments (UEs) utilizing Maximum Ratio precoding. The signal transmitted by AP l at time instant n is defined as [4]:

$$x_l\left[n\right]=\sqrt{{\displaystyle \frac{p_d}{K}}}\sum _{i=1}^K{\widehat{h}}_{il}\left[n\right]\sqrt{{\displaystyle \frac{{\mu }_{il}}{q_i\left[n\right]}}},$$

(98)

where $q_i\left[n\right]\sim \mathcal{CN}(0,1)$ represents the symbol transmitted to UE i, which remains consistent across all APs. Here, $p_d$ represents the maximum downlink transmission power for a single AP, and $0\le {\mu }_{il}\le 1$ values are power control coefficients selected to satisfy the downlink power constraint, ensuring $\mathbb{E}\left[|x_l{\left[n\right]|}^2\right]\le p_d$.

The received signal at the k-th UE at time n is described by Equation (99), where

$$\begin{array}{c}\hfill {\mathrm{SIN}\mathrm{R}}_{coh,k}\left[n\right]={\displaystyle \frac{{\rho }_k^2{\left|{\sum }_{l=1}^L\sqrt{{\mu }_{kl}\mathrm{tr}\left(Q_{kl}\right)}\right|}^2}{p_d{\sum }_{i=1}^L{\sum }_{j=1}^L{\mu }_{il}\mathrm{tr}\left(Q_{il}{R}_{kl}\right)+{\rho }_k^2{p}_d{\sum }_{i\in P_k\setminus \left\{k\right\}}{\left|{\sum }_{l=1}^L\sqrt{{\mu }_{il}}\mathrm{tr}\left({\overline{Q}}_{kil}\right)\right|}^2+{\sigma }_d^2}}.\end{array}$$

(99)

In this expression, ${\rho }_k$ denotes the effective SINR for UE k, and ${\sigma }_d^2$ represents the noise variance.

5.2. Non-Coherent Downlink Transmission

This section explores the concept of non-coherent joint transmission in downlink CF-MIMO systems. This method enables each Access Point (AP) to transmit distinct data symbols to each User Equipment (UE), thereby reducing the phase synchronization requirements compared to coherent transmission. Utilizing Maximum Ratio precoding, the signal transmitted by AP l is expressed as

$$x_l\left[n\right]=\sqrt{p_d}\sum _{i=1}^K{\widehat{h}}_{il}\left[\lambda \right]\sqrt{{\mu }_{il}}{q}_{li}\left[n\right],$$

(100)

The signal-to-interference-plus-noise ratio (SINR) received by UE k at time n is detailed in Equation (101), where

$$\begin{array}{c}\hfill {\mathrm{SIN}\mathrm{R}}_{nc,k}\left[n\right]={\displaystyle \frac{{\rho }_k^2{\sum }_{l=1}^L{\mu }_{kl}{\left|\mathrm{tr}\left(Q_{kl}\right)\right|}^2}{p_d{\sum }_{i=1}^L{\sum }_{j=1}^L{\mu }_{il}\mathrm{tr}\left(Q_{il}{R}_{kl}\right)+{\rho }_k^2{p}_d{\sum }_{i\in P_k\setminus \left\{k\right\}}{\sum }_{j=1}^L\sqrt{{\mu }_{il}}{\left|\mathrm{tr}\left({\overline{Q}}_{kil}\right)\right|}^2+{\sigma }_d^2}}.\end{array}$$

(101)

In this equation, ${\rho }_k$ represents the effective SINR for UE k, while ${\sigma }_d^2$ accounts for the noise variance.

5.3. Downlink Statistical Channel Cooperation Power Control

In [13], the power control policy is extended with a pre-determined function that includes global statistical channel information, leading to SCCPC coefficients as

$${\mu }_{kl}={\displaystyle \frac{{\beta }_k^{-1}}{{\sum }_{i=1}^K\mathrm{tr}\left(Q_{il}\right){\beta }_j^{-1}}},k=1,\dots ,K,l=1,\dots ,L,$$

(102)

where

$${\overline{\beta }}_k={\displaystyle \frac{{\sum }_{l=1}^L{\beta }_{kl}}{L}}.$$

(103)

In this work, ${\beta }_{kl}=\mathrm{Tr}\left(A\right)$, where A is given in (20).

6. Results

In this section, we evaluate the performance of MR, FZF, MMSE, and PZF in comparison to their corresponding DMD variants in the context of CF-MMIMO systems. We investigate how the distribution of antennas and the number of pilot sequences impact system performance, utilizing the closed-form SE expression with LSFD across various combining schemes. Our simulation setup mirrors that of a previous study [44], where L APs and K UEs are uniformly distributed within a circular simulation area with a radius of 1 km. The maximum transmit power for each UE is capped at $p_{\max }=100$ mW, with a channel bandwidth of $B=20$ MHz and a speed $v_{min}\le v\le v_{max}$. Each AP is equipped with N antennas. Pilot sequences, ${\tau }_p$, fewer than K, are randomly allocated to the UEs. Following the approach outlined in [43], the simulations use Montecarlo, with 256 scenarios, and each scenario is used 100 times. Based on the 3GPP model [45,46], the likelihood of a Line-of-Sight (LoS) path is influenced by the distance. The probability of having a LoS connection for the channel between UE $(l,i)$ and $BSj$ is

$$\begin{array}{c}\hfill Pr\left(LoS\right)=\left\{\begin{array}{cc}{\displaystyle \frac{300-d_{li}^j}{300}},\hfill & 0<d_{li}^j<300\mathrm{m}\hfill \\ 0,\hfill & d_{li}^j>300\mathrm{m}\hfill \end{array}\right.\end{array}$$

(104)

If the LoS path exists, then the corresponding large-scale fading coefficient is modeled (in dB) as

$${\beta }_{li}^j=-30.18-26{log}_{10}\left(d_{li}^j\right)+F_{li}^j,$$

(105)

We apply the heuristic UL power control policy from [3], where the transmit power UE k in the cell j is

$$\begin{array}{c}\hfill p_{jk}=\left\{\begin{array}{cc}{p}_{max}^{ul},\hfill & \Delta >{\displaystyle \frac{{\beta }_{j{k}^j}}{{\beta }_{j,min}^j}}\hfill \\ p_{max}^{ul}\Delta {\displaystyle \frac{{\beta }_{j,min}^j}{{\beta }_{jk}^j}},\hfill & \Delta \le {\displaystyle \frac{{\beta }_{jk}^j}{{\beta }_{j,min}^j}}.\hfill \end{array}\right.\end{array}$$

(106)

Various simulations are performed, taking into account the Doppler effect. Here, $T_s$ denotes the sampling time, while $f_{D,k}=\left({\displaystyle \frac{v_k{f}_c}{c}}\right)$ represents the Doppler shift for a UE moving at velocity $v_k$, where $f_c$ and c are the carrier frequency and the speed of light, respectively [19]. Our primary focus is on addressing the issue of channel aging and mitigating its effects through time-delay correlation. In all scenarios, the maximum normalized Doppler shift $f_D{T}_s$ is utilized [23].

6.1. Evaluation of the Applicability of DMD for Channel Estimation in CF-MMIMO Wireless Systems

In this section, we evaluate the applicability of Dynamic Mode Decomposition (DMD) for Channel Estimation in CF-MMIMO wireless systems through several experiments. Standard metrics, such as the Cumulative Distribution Function (CDF) of the per-UE uplink/downlink Spectral Efficiency (SE), are used for performance evaluation. These CDFs are applied to various methods, including MRC, MRC-DMD, FZF-DMD, MMSE, MMSE-DMD, PFZF, and PFZF-DMD. To assess adaptability, we employ a highly dynamic scenario as described in Section [referenced section]. Furthermore, PiDMD and mpeDMD are compared with traditional predictors, such as Machine Learning (ML) algorithms and the Kalman Filter (KF). The impact of the number of time–frequency snapshots, denoted as ${\tau }_p$, is also investigated.

In wireless networks, algorithm efficiency is commonly evaluated by calculating Spectral Efficiency. Figure 15 presents the CDF of the per-user downlink SE for different numbers of BS antennas using various combining schemes: MRC, MMSE, FZF, PFZF, and their DMD-enhanced versions—MRC-piDMD, MMSE-piDMD, FZF-piDMD, and PFZF-piDMD. The parameters used are $L=100$, $K=16$, $N=8$, ${\tau }_p=6$, ${\tau }_S=40$, $p_k=100$ mW, and $50\le v\le 150$ for each UE. There is a notable performance disparity between MR-based combinations and ZF-based schemes, particularly for UEs with high channel gains, largely due to inter-user interference. FZF and PiFZF schemes effectively mitigate this interference, with FZF combining leading to a 16% improvement in average SE. Moreover, FZF-PiDMD outperforms FZF, and PFZF-PiDMD achieves a 33%-likely SE improvement over PFZF, primarily by offering enhanced protection to weaker UEs with lower channel gains. The increase in the number of antennas per base station (BS) from 4 to 64 results in an average improvement in algorithm performance of 28%.

Figure 16 compares the UL CDF of the SEs with different estimators, MR, MR-DMD, ZF, ZF-DMD, MMSE, MMSE-DMD, and PZF, with PZF-DMD combining the scheme with L = 128, K = 128, N = 64, ${\tau }_p$ = 16 and $p_{max}$ = 100 mW with different number of UEs, K = [10, 40]. As the number of UEs decreases while maintaining the same pilot length, both interference and pilot contamination are reduced, resulting in higher SEs for each UE. In this case, the density of UEs increased raises the number of pilots. At first glance, it is surprising that MR-mpeDMD under these conditions presents very good results. While MR is optimal in the presence of noise, it is not as effective with both noise and interference. This reveals another important characteristic of MR-mpeDMD: by intrinsically separating stronger signals from weaker ones and treating the weaker signals as noise, it proves to be extremely efficient. MMSE-mpe, ZF-mpeDMD, and PZF-mpeDMD also benefit from the separation of strong and weak signals using mpeDMD.

Figure 17 depicts the average uplink Spectral Efficiency (SE) per User Equipment (UE) plotted as a function of the pilot sequence length, ${\tau }_p$. The simulation parameters include L = 100 Access Points (APs), N = 4 antennas per AP, K = 40 UEs, a pilot sequence length (${\tau }_p$) of 40, and spatially correlated Rayleigh fading with an Angular Spread Distribution (ASD) of ${\sigma }_{\varphi }={\sigma }_{\theta }={15}^{\circ }$ and $50<v<150$. FZF-mpeDMD shows great superiority over the other algorithms primarily due to the minimal number of pilots required for its operation, which allows for the optimization of data transmission. On the other hand, MMSE-mpeDMD presents almost constant behavior. The hypothesis is that MMSE is already optimal in environments with high interference, so adding DMD does not provide a significant improvement in SE. FZF-PiDMD also shows high performance around a length of 20 ${\tau }_p$, but then it decreases because the payload starts being affected by the increase in pilots.

Figure 18 shows the relative error of the average Spectral Efficiency (SE) derived from the asymptotic closed-form expression compared to the ergodic average SE, plotted against the number of User Equipments (UEs). The simulation parameters include the number of antennas (M) equal to the number of UEs (K), 64 Access Points (APs), a pilot sequence length (${\tau }_p$) equal to half the number of UEs (K/2), and a transmit power ($p_k$) of 100 mW for each UE. The results indicate superior performance of DMD methods. As seen in Figure 18, the approximated SE per UE becomes more accurate as N and K increase. MMSE is well known for its optimal performance in a noisy and interference-prone environment. In terms of error calculation, its performance is also very good; however, MMSE-PiDMD and MMSE-mpeDMD slightly surpass it. This is because, in addition to the optimal effect of MMSE, DMD is capable of discriminating those signals that are more relevant to the algorithm, ultimately reducing the error.

To assess adaptability, a highly dynamic scenario is utilized. PiDMD and mpeDMD are compared with predictors such as Machine Learning (ML) and Kalman Filter (KF) vs. UE speed. The results illustrated in Figure 19 indicate that piDMD demonstrates superior performance due to its computational speed and accuracy. While mpeDMD generally provides better accuracy for system identification compared to piDMD, the significant computational overhead required for its calculation greatly affects its overall performance. Nevertheless, DMD methods outperform traditional ones.

6.2. Performance Evaluation of DMD in a System Considering Delay–Doppler Correlated Channel

In this section, several experiments are conducted to demonstrate the impact of a system with delay–Doppler correlation and the performance of DMD-based algorithms in this specified scenario. Experiments aim to show Spectral Efficiency under conditions of a correlated channel, analyzing effects for both uplink and downlink transmissions. Correlation always exists in practice and is highly dependent on every communication scenario, primarily generated by ASD, which is quantified in this section. Different Doppler values in the channel are also studied. Some studies suggest mitigating the correlation by increasing the number of pilots ${\tau }_p$; therefore, the number of pilots is tested to determine this relationship among DMD algorithms. Additionally, the impact of channel correlation on coherent and non-coherent downlink transmissions is analyzed.

Channel correlation is largely influenced by ASD, making it insightful to investigate the relationship between Spectral Efficiency (SE) and ASD. Figure 20 depicts the average uplink SE per UE as a function of ASD for azimuth and elevation angles. In this graph, we use ${\sigma }_{\varphi }={\sigma }_{\theta }$. We consider different operation modes of CF-MMIMO with $L=100$, $N=4$, $K=40$, and ${\tau }_p=10$. Results for uncorrelated Rayleigh fading are shown as a reference. Surprisingly, mpeDMD can leverage correlated channels and outperform the case without correlation.

Figure 21 illustrates the 95%-likely per-user non-coherent downlink Spectral Efficiency (SE) against the value of the length of the frames ${\tau }_s$. The simulation parameters are set to $L=100$, $K=20$, $M=2$, and ASD = 30°. Simulations show that asymptotically, DMD-based methods converge faster as the delay increases. PiDMD and mpeDMD exhibit better performance than MR and FZF. MR and ZF have good performance when $f_D{T}_s=0$ but degrade rapidly as the Doppler effect increases. DMD models also degrade with a higher Doppler effect but much less than non-DMD schemes. The hypothesis for why this happens is that DMD methods separate strong signals from weaker ones. Weaker signals are treated as noise, and the stronger signals, which are just a few, are easily aligned, allowing good signal alignment at the receiver.

Figure 22 depicts the uplink (UL) transmission and its dependence on the number of pilot symbols (${\tau }_p$) with $L=100$, $K=20$, $N=2$, and ASD = 30°, considering a delay–Doppler correlated channel. The delay is uniformly distributed in $[0,0.01]$ ms, and the normalized Doppler is $f_D{T}_s=0.002$. The graph shows that piDMD and mpeDMD initially grow very rapidly, and then their growth becomes nearly negligible. MR, PZFZ, and ZF exhibit linear growth. This clearly indicates that DMD-based methods require fewer pilots; increasing the number of pilots beyond a certain limit does not substantially increase the SE.

Figure 23A illustrates the CDF of per-user uplink Spectral Efficiency (SE) with full power. Figure 23B shows the effect of normalized Doppler shifts of [$f_d{T}_s=0$, $f_d{T}_s=0.004$] on uplink SE. It is evident that DMD systems outperform non-DMD systems at both CDF and 95%-likely Uplink SE with Doppler effect. Increasing the normalized Doppler shift $f_d{T}_s$ from 0 to $0.004$ results in a 34% median SE loss for DMD, a 44% median SE loss for MRC, and a 60% median SE loss for PZFZ. This is because DMD leverages network-wide fading statistics to calculate coefficients, effectively minimizing interference and mitigating the effects of channel aging. It is noteworthy that the ZF systems, despite employing a stricter capacity bound through channel estimates in data detection, still perform poorly under channel aging. Another point to highlight with respect to the Doppler effect is that this effect is directly related to the speed of the UEs. Therefore, it is a direct indicator of the behavior of DMD in a static or dynamic system. As can be seen, the higher the speed, the better the performance of DMD compared to non-DMD schemes.

Figure 24 illustrates the CDF of the per-user downlink Spectral Efficiency (SE) for both coherent and non-coherent transmissions at full power, with normalized Doppler shifts of $f_d{T}_s=0$, $f_d{T}_s=0.0005$, $f_d{T}_s=0.001$, $f_d{T}_s=0.0015$, and $0.002$. The results demonstrate that coherent transmission achieves significantly higher SE compared to non-coherent transmission, irrespective of whether the UEs are stationary or in motion. An increase in the normalized Doppler shift $f_d{T}_s$ from 0 to $0.002$ leads to a 37% decrease in median SE for coherent transmission and a 42% decrease for non-coherent transmission. This disparity arises because, in non-coherent transmission, UEs detect signals from the APs sequentially, making them more susceptible to the impacts of channel aging.

6.3. Evaluation of the Performance in Overall System Capacity Achieved through DMD-Based Channel Estimation

In this section, we aim to determine the overall behavior of the network with the use of DMD. To achieve this, we compare the performance of DMD models with non-DMD models and appropriately adjust the most important parameters that affect network behavior. These parameters include the number of base stations (BSs) and the number of users. Additionally, in a highly concentrated scenario, we analyze the effect of the number of samples taken on the performance of the two major system groups we are using: DMD and non-DMD.

Figure 25 displays the Average Spectral Efficiency (SE) per User Equipment (UE) as a function of the number of base stations (BSs), denoted by L. The parameters used in the simulation include ${\tau }_p=7$, $K=10$, $N=8$, and $p_k=100$ mW for each UE. As the network becomes denser, local network traffic can increase, leading to an overall increase in global traffic. However, the rise in interference can cause a decrease in global traffic. Experimental results demonstrate that DMD optimizes global network traffic exceptionally well, achieving an almost linear increase with the density of both BSs and UEs, as shown in Figure 25 and Figure 26. This is attributed to the efficient handling and implicit power control in the DMD model, which optimally manages the separation of strong signals with appropriate handling and weak signals.

Figure 26 illustrates the CDF of the uplink (UL) sum Spectral Efficiency (SE) for a system with K = 50 users, plotted as a function of the number of base station (BS) antennas. The performance is evaluated for various channel estimators, providing insight into how the distribution of the SE of the system changes with different estimation techniques.

Figure 27 illustrates the average Spectral Efficiency (SE) per User Equipment (UE) against the pilot sequence length, ${\tau }_p$, with ${\tau }_S=20$, $K=10$, $N=8$, and $p_k=100$ mW for each UE. It shows the impact of Channel Estimation on the global network performance based on the number of pilots ${\tau }_p$. Generally, a higher number of pilots improves Channel Estimation; however, it also reduces the network’s effective traffic. The figure illustrates that non-DMD algorithms still increase global network traffic as the number of pilots increases. DMD schemes reach a point where the increase in pilots greatly decreases the global network traffic. Intuitively, this is because DMD requires significantly fewer pilots for accurate Channel Estimation.

6.4. Computational Intensity Analysis

In order to better understand the amount of processing required in each of the methods, three scenarios are used:

• Scenario 1: $L=40$, ${\tau }_p=16$, $K=100$, $N=8$, $p_k$ = 100 mW for each UE, $v=100$ k/h, delay = [0, 001] ms and Dopper $f_d{T}_s=0.002$.
• Scenario 2: $L=40$, ${\tau }_p=16$, $K=100$, $N=8$, $p_k$ = 100 mW for each UE, $v=100$ k/h, delay = 0 ms and Dopper $f_d{T}_s=0$.
• Scenario 3: $L=50$, ${\tau }_p=16$, $K=500$, $N=8$, $p_k$ = 100 mW for each UE, $v=150$ k/h, delay = [0, 001] ms and Dopper $f_d{T}_s=0.002$.

The results are obtained through a hundred simulations of a hundred scenarios [7].

The summary of this section’s simulations and analysis is shown in

Table 1. Computational Intensity. Part 1.

Modulation Scheme	Number of Operations	Processing Time *
MR	3N − 1	≈20.0 ms
		≈18.0 ms
		≈28.0 ms
FZF	$O(N{K}^2+K^3)$ **	≈26.5 ms
		≈31.6 ms
		≈41.2 ms
PFZF	$O\left(N^2\right)$	≈22.3 ms
		≈19.0 ms
		≈26.2 ms
MMSE	$O\left(K^3\right)+O\left(N{K}^2\right)$ ***	≈31.4 ms
		≈31.2 ms
		≈43.1 ms

* Time is shown sequentially according to the scenarios. ** In this case, the matrix is of size NxK. *** N is the number of antennas and K is the number of users.

,

Table 2. Computational Intensity 2. Part 2.

Modulation Scheme	Number of Operations	Processing Time *
MR-PiDMD	$O\left(n^3\right)+3N-1$	≈31.2 ms
		≈31.4 ms
		≈42.3 ms
FZF-PiDMD	$O\left(n^3\right)+3N-1+O(N{K}^2+K^3)$ **	≈28.9 ms
		≈27.8 ms
		≈38.2 ms
PFZF-PiDMD	$O\left(n^3\right)+O\left(n^2\right)$	≈35.2 ms
		≈36.2 ms
		≈45.2 ms
MMSE-PiDMD	$O\left(n^3\right)+3N-1O\left(K^3\right)+O\left(N{K}^2\right)$	≈31.2 ms
		≈35.9 ms
		≈46.7 ms

* Time is shown sequentially according to the scenarios. ** N is the number of antennas, and K is the number of users.

and

Table 3. Computational Intensity 3. Part 3.

Modulation Scheme	Number of Operations	Processing Time *
MR-mpeDMD	$O\left(n^3\right)+(3N-1)$	≈48.2 ms
		≈49.1 ms
		≈51.3 ms
FZF-mpeDMD	$O\left(n^3\right)+(N{K}^2+K^3)$	≈46.3 ms
		≈48.1 ms
		≈45.2 ms
PZFZ-mpeDMD	$O\left(n^3\right)+O\left(N^2\right)$	≈58.1 ms
		≈62.8 ms
		≈67.1ms
MMSE-mpeDMD	$O\left(n^3\right)+O\left(K^3\right)+O\left(N{K}^2\right)$	≈65.1 ms
		≈67.3 ms
		≈71.2 ms

* Time is shown sequentially according to the scenarios.

. The time presented is the time for each scenario.

Table 2 and Table 3 show the number of operations for matrices of the same order. However, the processing time is quite different. The reason is that the order of the operations needed for computation hides many minor operations, which are time-consuming.

7. Discussion

This study delves into the potential of DMD to revolutionize Channel Estimation within CF-MMIMO systems. Through extensive simulations and analysis, several critical insights have been gleaned.

Firstly, PiDMD/mpeDMD methods have consistently outperformed traditional techniques across various metrics. The core advantage of DMD lies in its ability to capture and utilize the dynamic modes from received signal snapshots, effectively tracking the evolving characteristics of wireless channels over time and space. This dynamic tracking capability has proven particularly beneficial in highly mobile scenarios, where traditional methods often struggle to maintain accuracy.

Our experiments illustrate that DMD-enabled Channel Estimation significantly enhances system performance. Notably, DMD methods demonstrate superior Spectral Efficiency and more effective interference mitigation compared to conventional techniques. This improvement is attributed to the robust extraction and application of CSI, which is crucial for optimal beamforming and spatial multiplexing strategies.

Additionally, the simulations have shown that DMD methods lead to better scalability in processing signals within large-scale CF-MMIMO deployments. This scalability is essential for practical implementations, as it ensures that the system can handle increased user densities and higher mobility without a corresponding increase in computational complexity.

The impact of DMD on system capacity is another highlight of this study. By accurately estimating and adapting to the channel conditions, DMD methods facilitate higher system capacity and improved communication reliability. This is particularly evident in dynamic and dense wireless environments, where the ability to swiftly adapt to changing conditions is paramount.

Moreover, the research underscores the importance of considering time- and delay-correlated channels in the evaluation of CF-MMIMO systems. The inclusion of realistic channel models in our simulations provides a more accurate depiction of system performance, ensuring that the results are applicable to real-world scenarios.

This study validates the transformative potential of incorporating DMD into CF-MMIMO systems. The demonstrated gains in Spectral Efficiency, interference mitigation, and system capacity underscore the significant advantages of DMD-based methods over traditional methods.

Future work may include topics such as Reconfigurable Intelligent Surface (RIS) [47,48], which is a very recent subject and will be explored in future works. In [49], massive synchrony in distributed Antenna Systems is explored, an issue that is a topic that could also be addressed using DMD methods, just to mention two topics among several others that could be studied using DMD and its variants. It is also worth noting that DMD-based methods themselves are a field of much recent research.

8. Conclusions

In this paper, we have presented a novel approach to address the challenges associated with Channel Estimation in CF-MMIMO wireless systems. Leveraging the power of DMD, we have introduced a data-driven methodology that captures the temporal dynamics of the wireless channel, offering a promising solution to the complexities posed by rapidly changing environments.

The integration of DMD/PiDMD/mpeDMD into the CF-MMIMO framework holds significant potential. By extracting dynamic modes from received signal data, we gain insights into the evolving behavior of the channel, enabling more accurate and adaptable CSI estimation. Through extensive simulations, we have demonstrated the effectiveness of the proposed DMD-based approach in enhancing system performance, particularly in scenarios with dynamic and dense wireless environments.

Our findings underscore the transformative impact of infusing Dynamic Mode Decomposition (DMD) into massive Multiple-Input–Multiple-Output (MIMO) wireless systems. The ability to predict and adapt to changing channel conditions opens up new avenues for optimizing transmission strategies, improving Spectral Efficiency, and ensuring reliable communication. Additionally, the integration of DMD may pave the way for novel applications in wireless systems beyond 5G and 6G, where dynamic and agile communication environments are paramount.

As the wireless landscape continues to evolve, our research indicates that integrating DMD with CF-MMIMO systems presents a promising avenue for achieving robust and efficient wireless communication. Nonetheless, further investigation is required to delve into the complexities of DMD-based Channel Estimation across a variety of real-world scenarios. Additionally, addressing challenges such as computational complexity and adaptation to different propagation environments remains imperative.

In closing, the fusion of advanced data-driven techniques such as DMD, PiDMD, and mpeDMD with cutting-edge wireless technologies holds the potential to reshape the future of wireless communications. This integration promises to usher in an era of adaptable, intelligent, and high-performance wireless networks.