EP4680994A1 - Method of joint communication and environment perception, and system implementing the method - Google Patents

Abstract

A method of processing wireless communication signals for use in joint communication and environment perception in a region of interest (ROI) is presented. The method alternatingly uses a first and a second iterative process for determining an initial representation of a voxelated environment of the ROI and for determining the transmitted signals. The first and second iterative processes may implement a Gaussian approximation belief propagation process.

Description

202401191 1 METHOD OF JOINT COMMUNICATION AND ENVIRONMENT PERCEPTION, AND SYSTEM IMPLEMENTING THE METHOD FIELD OF THE INVENTION The invention relates to the field of environment sensing or environment mapping, in particular to such sensing using wireless communication signals. More specifically, the invention relates to a method of processing communication signals for environment perception, to a computer program product implementing the method, to a computer-readable storage medium storing the computer program product, and to a system configured to execute the method. Throughout this specification the term environment sensing will be used for the various expressions widely used for capturing information about an environment for creating a three-dimensional representation thereof. NOTATIONS Scalar values are denoted herein by lowercase letters in italics, as in x, while complex vectors and matrices are denoted by boldface lowercase and uppercase letters, as in x and X, respectively. (·)^T and (·)^* denote the transposition and complex conjugation operators respectively, and diag(·) denotes the diagonalization operator. | · | denotes the absolute value operator whereas || ||ℓ denotes the ℓ-th norm. ^x(x) and Varx(x) respectively denote the expectation and variance operator of x with respect to the distribution of x given by ℙx(x). ℝ and ℂ denote the real and complex number fields respectively, and ^(μ, ν) and ^ ^(μ, ν) denotes the real and complex Gaussian distributions with mean μ and variance ν. BACKGROUND Joint Communication and Sensing (JCAS) is a technique in wireless communications with the objective of retrieving information about the environment from the signal scattering which is present in the effective channel state information (CSI), e.g., due to objects in the environment, blockage, user activity, etc., while simultaneously achieving data communication. Most known JCAS methods exploit radar technology to infer information about the environment. This is also known as joint radar and communication (JRC). 202401191 2 Various methods are known in JRC, including alternating or sharing spectrum between radar and communication signals, using standard radar signals to embed information, extracting radar parameters from standard communication signals, or even designing new waveforms suited for both tasks. These techniques are highly based on conventional radar signal processing (e.g., ambiguity function estimation) and dependent on the radar frequency-delay properties and prone to similar challenges. Spatial representations of environments are widely used in the field of robotic vision, localizing and mapping, as well as in computer graphics and medical imaging. One known technique involves discretely approximating a true environment by a voxelated occupancy grid. Typically, voxels are regular cubes arranged in a 3D space covering the environment, which are assigned information about the space they occupy. The resulting discreteness of the voxelated space is well suited for 3D modelling where, depending on the size of the unit voxel, objects in the region of interest may be represented by rough estimates or more accurate shapes. Figure 1 a) shows an exemplary environment that is represented, in voxelated representations having different resolutions, through use of voxels of different sizes in figures 1 b) and 1 c). In figures 1 b) and c) 3D voxelated occupancy grids are shown, where the total region of interest (ROI) is defined as a cuboidal space of dimensions L_x ^ L_y ^ L_z, each denoting the dimensional lengths of the x, y, z-axes in meters, respectively. The entire ROI is subdivided into a grid consisting of NV ≜ Nx · Ny · Nz voxels, where N_x ≜ ^{^^} ^_^, N_y ≜ denote the number of voxels, or partitions, per x, y, z -axes, respectively, and L_V is the edge length of a voxel cube in meters, corresponding to the image resolution. If represented as a tensor of three dimensions (Nx ^ Ny ^ Nz), the voxelated occupancy grid directly represents a discretized model of the ROI, as shown in Fig.1 b) and c), where the size of the voxels corresponds to the image resolution. The elements of the three-dimensional tensor indicate the occupancy of the voxels, and thus, whether that portion of the space is empty or filled with a given material. 202401191 3 In addition to the estimation accuracy flexibility provided by voxelated grids, one of the biggest advantages of the grid-based models is the simplicity of the data representation, especially as opposed to 3D vertex-based or point cloud-based methods, which makes voxelated grids even more favourable from a machine- learning point of view. Recent developments in communication technology have identified modelling the space between a transmitter and a receiver by a voxelated environment as beneficial for JCAS, sometimes also referred to as integrated sensing and communication, or ISAC. Throughout this specification, the terms JCAS and ISAC are used interchangeably. To this end, the electromagnetic scattering behaviour of objects in the true environment may be added to the 3D geometric information as provided by the classic voxelated occupancy grid, for tailoring the communication channel modelling method that is utilized in a wireless communication scenario. For considering the scattering behaviour, each voxel is assigned a voxel occupancy coefficient _vk ∈ {0, 1} with k ∈ {1, …, NV }, where v^k = 0 indicates that the k-th voxel is empty, i.e., the corresponding environment is free-space, and v^k = 1 indicates that the k -th voxel is occupied by a scatterer object, e.g., the table, chair or the object on the wall, as illustrated in Fig.1. Note that the binary voxel occupancy coefficients may be extended to complex voxel scattering coefficients, i.e., v^k ≜ βk·(^)*+, ∈ ℂ, to also capture the effect incurred to the reflected electromagnetic waves by the occupied voxels. The constants depend not only on the material itself, but also on the frequency and the angle of incidence of propagating signals, and capture the effect of the material occupying a given voxel onto the electromagnetic wave reflected or refracted by it. In other words, the values of the voxel scattering coefficients are expected to be highly dependent on the electromagnetic characteristics of the scatterer object and the impinging wave, which may be empirically measured and modelled as a function of the properties such as frequency and material. 202401191 4 An exemplary regular voxelated 3D space with some scatterer objects accommodating a single access point (AP), a single reconfigurable intelligent surface (RIS), and multiple single-antenna user equipment (UEs) is discussed by X. Tong, Z. Zhang, J. Wang, C. Huang and M. Debbah in “Joint Multi-User Communication and Sensing Exploiting Both Signal and Environment Sparsity,” IEEE Journal of Selected Topics in Signal Processing, vol.15, no.6, pp.1409-1422, Nov.2021. The multiple UEs are communicating to the AP via sparse code multiple access (SCMA) over multiple frequency subcarriers and over multiple transmission instances, via line-of-sight (LOS) paths and non-line-of-sight (NLOS) paths from the UEs, to the scatterers, to the RIS, then finally to the AP. Figure 2 a) and b) shows a general concept of LOS paths and NLOS paths in a voxelated space, where all paths are available and the reflection angles are within a range that actually reflects impinging electromagnetic waves. In figure 2 a) the LOS path is the direct path between the user equipment (UE) and the access point (AP), indicated by the solid line, while the two occupied voxels in the ROI reflect signals emitted by the UE towards the AP. The dashed lines represent the NLOS UE-to- voxel path, and the dotted lines represent the NLOS voxel-to-AP path. In figure 2 b) one of the voxels that was filled is now empty, which will obviously have an effect on the scattering coefficients of that voxel. The dashed line from the empty voxel to the AP represents the path that was present in figure 2 a) and is not to be taken as an existing path. Figure 3 shows a schematic representation of the 3D space considered in the known exemplary system and method, including the RIS. Here, the signal reflected off the RIS towards the AP is shown in a dash-dotted line, to highlight its specific origin. According to the known exemplary system model and assumptions, the signals received at the AP would not only carry the UEs’ payload, or data, but also contain the scattered path information which can be used to retrieve, or perceive, the environment. To aid the estimation, UEs are assumed to transmit a known length of pilot symbols at the beginning of the transmission interval, i.e., a known preamble. The known method aims to return the estimates of the transmitted SCMA codes of the UEs, i.e., user detection, or identification, and a binary voxel occupancy grid 202401191 5 which corresponds to the estimation of the environment, e.g., 0 for a voxel representing empty space and 1 for a voxel representing an occupied space. The known system may be represented mathematically by the following model Y = (H + PRQVB) [XpX] + W where Y is the received signal matrix over all receive antennas and APs across all symbol instances, H is the UE-to-AP LOS channel matrix, P is the RIS-to-AP NLOS channel matrix, R is the RIS reflection coefficient matrix, Q is the voxel-to-RIS NLOS channel matrix, B is the UE-to-voxel NLOS channel matrix, V is the voxel occupancy matrix, Xp is the pilot symbol matrix, X is the data symbol matrix, and W is the AWGN noise matrix. Since the matrices P, R, Q, and B are known, the model can be simplified to Y = (H + AVB) [XpX] + W where A = PRQ is the effective voxel-to-AP NLOS channel matrix. From this, the main joint communication and sensing problem is formulated, i.e., to estimate the environment matrix V and data symbol matrix X, with known channel matrices H, A, B, and known pilot matrix Xp. To do so the prior art method uses three modules: ^ an environment estimation module, using a linear generalized approximate message passing (GAMP) algorithm which uses the known channels, with given data symbols X, either pilots or estimated symbols from a previous iteration, for estimating the environment V, on the model Y = (H + AVB) [X] + W ^ an effective channel reconstruction module which, by integrating known channels and the estimated environment and using V estimated by the GAMP algorithm, combines the information with known H, A, B to calculate the effective total channel 202401191 6 .^/ = H + AVB ^ a data signal estimation module which, using a linear SCMA message passing algorithm (MPA) and the estimated effective channel, estimates the unknown data symbols X, 0 = .^/1 + 3 These three modules are iterated in a sliding window fashion on the time domain as illustrated in Figure 4, such that first, the environment estimation is carried out using the known pilot symbols, but for the data signal estimation, the window slides for estimating the unknown data symbols, and so on. However, like other known method of JCAS in voxelated environments, the prior art method discussed further above relies on the idealised assumptions that all paths between all transmit and receive antennas are fully available, i.e., no paths are blocked, that the channel gains for UEs-to-AP LOS paths are known, that the channel gains for the UEs-to-voxels, voxels-to-RIS, RIS-to-AP NLOS paths are known, that the RIS reflection coefficients are known, that the voxelated environment model is binary, i.e., only discrete occupancy values 0 or 1 are possible, that all scattered paths have realistic reflection angles, that an SCMA communication scheme is used, and that only a single AP is present. These assumptions pose severe limitations to applying the known method to actual 3D real-world environments, giving rise to the need for an improved method of processing wireless communication signals for use in environment sensing, a receiver configured to execute the improved method, and a communication system comprising one or more such receivers. SUMMARY OF THE INVENTION This need is addressed by the method of claim 1, the system of claim 7, and the computer program product of claim 9. A corresponding computer-readable storage medium and a vehicle comprising an improved receiver in accordance with the invention are presented in claim 10. Advantageous embodiments and developments are described in respective dependent claims. 202401191 7 Prior to describing the various aspects of the invention, the concept of unavailable paths and unrealistic paths in the voxelated environment will be elucidated. Figure 5 a) and b) show exemplary voxelated spaces where some paths may be rendered infeasible or unavailable due to various physical phenomena. In both, figure 5 a) and b), the direct solid line from user equipment UE to the access point AP represents a LOS path, which is available. It is obvious that, if an occupied voxel lies directly in line with the path between the E and the AP (not shown in the figure), the corresponding path will be unavailable. If the scattering angle between the incident and the reflected path is too large and exceeds a critical angle, e.g., an angle close to 180°, the corresponding NLOS path will be unrealistic and, in effect, be unavailable, as shown in figure 5 a). The dashed line from voxel to the AP represents the reflected voxel-to- AP NLOS path that has an unrealistic reflection angle, i.e., is not available. The determination of the critical angle is dependent on the electromagnetic properties of the RF wave and the environment, e.g., on the operating frequency, and hence the present invention proposes a simplified model to approximately incorporate such phenomena into the channel matrices of the voxelated grid environment model. Unavailability of a path may also be caused by a skewed surface, as shown in figure 5 b), where the signal is reflected away from the AP. Like in figure 5 a), the theoretically possible NLOS path reflected from the empty voxel in figure 2 b) will be unavailable due to the object in the voxel reflecting the incident wave away. Further prior to describing the various aspects of the invention, an underlying channel model will be described. It is assumed that the ROI comprises NU single-antenna UEs, and NA multi-antenna APs equipped with NR receive antennas each. As illustrated in Fig.5, the effective channels between the UEs and APs contain two distinct types of components, the LOS components, which is the direct path between the UEs and APs, and the NLOS components, which encompass the paths reflected at occupied voxels corresponding to parts of objects in the environment, i.e., paths scattered by voxels representing scatterer objects, as described further above. It is also assumed that the frequency band is sufficiently high that the power of paths reflected more than once is negligible such that NLOS components may be decomposed into two sub-paths, the UE-to-voxel sub-path and the voxel-to-AP sub- 202401191 8 path, which together with the voxel scattering coefficient comprises the aggregate NLOS channel. In light of the above, the effective channel between the NU single-antenna UEs and the ensemble of NANR receive antennas of all APs (i.e., for all NR antennas per each of the NA APs), is given by (1) where G ∈ ℂ^7879:7; is the effective channel matrix, H ∈ ℂ^7879:7;, A ∈ ℂ^{7879:7^}, and B ∈ ℂ^{7^:7;} are the constituting channel matrices for the UE-to-AP LOS path, voxel- to-AP NLOS sub-path, and UE-to-voxel NLOS sub-path, respectively, and v ∈ ℂ^{7^:<} is the vector containing all scattering coefficients of the voxelated grid. The elements of the channel matrices H, A, and B are assumed to follow a zero-mean complex Normal distribution with variances =_> ^?, and =_A ^?, respectively. It is noted that the channel model presented above can be straightforwardly extended to a multi-carrier system, yielding G(f) = H(f) + A(f) · diag(v(f) · B(f) ∀f ∈ℱ (2) with f denoting a specific frequency in the set ℱ of all subcarrier frequencies. It is important to notice that the channel model in the equation shown above assumes that all paths between UEs, APs, and voxels are fully available, like in the known system discussed above in the background section. However, in reality many paths may be rendered unavailable due to various physical phenomena, e.g., blockage by air-borne particles, absorption, path loss, or by the finite resolution of the voxelated model itself. Also, as described above, if the angle between the incident and reflected path is too large and exceeds the critical angle, the corresponding NLOS path will not be available, as shown in Fig.5 a). 202401191 9 Likewise, if an occupied voxel is directly in line with the LOS path, the corresponding path will not be available. The determination of the critical angle is dependent on the electromagnetic properties of the RF wave and the environment, e.g., on the operating frequency, and hence a simplified model is proposed to approximately incorporate such phenomena into the channel matrices of the voxelated grid environment model. As realistically modelling all paths may be computationally prohibitive, embodiments of the present invention may consider the statistical feasibility of paths for improving the voxelated grid model. To this end, the physical phenomena occurring at the reflection of propagating waves are considered, in particular, the fact that for any given frequency: i) a critical angle θ* exists such that, as illustrated in Fig.5 a), if the incidence angle θ > θ*, the wave is absorbed rather than reflected, and consequently the corresponding voxel-to-AP NLOS sub-path is not available, and ii) the curvature of the surface exposed to the impinging wave may be such that no signal is reflected towards an AP, as illustrated in Fig.5 b). As mentioned before, the complexity of modelling such phenomena at each voxel may be far too complex to carry out, especially if the resolution of the voxelated ROI is large. Employing a statistical approach for considering the angle between the impinging and reflected waves at each voxel, hereafter referred to as the scattering angle, and consequently, the availability of each voxel-to-AP NLOS sub-path, will significantly reduce the computational complexity. In light of the above, embodiments of the invention may use the following stochastic- geometric model for integrating the aforementioned phenomena into the channel matrices of the voxelated grid environment model. First, the positions of the UEs and the APs are discretized into the 3D grid of the voxelated environment model, such that their positions may be equivalently described by the voxel coordinates. Note that it is also assumed that the multiple 202401191 10 antennas of the respective APs are fully located within a single voxel, such that their angles-of-arrival (AoA) are assumed to be identical, while having different channel path coefficients. Given the 3D coordinates of the UE, AP, and the voxel as cU = [xU, yU, zU]^T ∈ ℝ³, cA = [xA, yA, zA]^T ∈ ℝ³, and cV = [xV , yV , zV ]^T ∈ ℝ³, respectively, the scattering angle θ of the path at the voxel may be computed as where the arccos(·) operator denotes the inverse cosine trigonometric function. The empirical probability distribution functions (PDFs) of the scattering angles θ can be obtained by evaluating the foregoing equation for all possible combinations of admissible locations of UE, AP and voxel within the ROI, respectively given by cU, cA and cV, and examples of the latter for an environment voxelated at various resolutions are shown in figure 6. It is visible in figure 6 that for a sufficiently large NV the distribution of scattering angles θ can be well modelled by a mixture of two beta distributions, namely, fx(x) = γ·β(a1, b1) + (1-γ)·(a2, b2), with support x ∈ [0; 180°] and where γ is a weighing factor and the quantities a₁, a₂, b₁, are shape parameters optimised to match the empirical data obtained by evaluating the foregoing equation, with cU, cA and cV taken randomly within the voxelated grid. Utilizing this empirical stochastic-geometric approach, the voxelated environment models used so far can be improved by the incorporation of random blockages of NLOS sub-paths, in proportion to the complement cumulative distribution of the approximated beta mixture PDF, and in accordance to the scattering angles at each voxel, as a function of a selected critical angle. As mentioned before, the scattering angle θ at all voxels is calculated for all pairs of UEs and APs. By introducing an arbitrary critical angle θ_crit ∈ [0°, 180°], a scattered path is determined to be unavailable if θ > θ_crit, and the corresponding paths are removed from the NLOS channel matrices. 202401191 11 The effect of the critical angle θcrit on the severity of puncturing on the channel matrices is illustrated in figure 7, where the result is obtained by numerical evaluations with random numbers and positions of UEs and APs for various voxelated grid resolutions, with the average channel puncturing rate normalized by the total number of path vertices (NA + NR). As expected, the puncturing rate shows a smooth increase between θcrit = 180° with no puncturing, and θcrit = 0° with full puncturing. This inherently includes the blocked LOS case, since the blockages may be regarded as a case with θ = 180°. As mentioned above, the true critical angle is dependent on the elaborate electromagnetic properties of the environment, but its determination is out of scope of this invention. However, a very interesting behaviour is observed where the severity of puncturing becomes effectively invariant to the resolution of the voxelated environment model, i.e., the size of the voxels, converging to the same relationship at a sufficiently high resolution. As illustrated in figure 7, the convergent curve may be very closely approximated by a scaled Gaussian curve, where a heuristic search yields the optimal parameterization of ^(−9.8, 54²) with a scaling factor . Drawing from the above, an efficient model of the puncturing behaviour is proposed by introducing the feasibility coefficient ξ ∈ [0, 1] which follows a Bernoulli distribution with probability pξ obtained by evaluating the scaled Gaussian distribution Then, independent and identically distributed (i.i.d.) feasibility coefficients are multiplied to each element of the channel matrices, and capture the behaviour of the infeasible paths being unavailable and punctured. In the following discussion an uplink communication scenario between the group of NU UEs and a total of NA APs, under the models described above, is considered, with the NA APs connected to a central processing unit (CPU) via error-free fronthaul links of unlimited throughput, such that the received signals at all NANR receive antennas are aggregated without loss of information or delay. Then, the aggregated received signal matrix Y, over NT discrete transmission instances (symbol slots) is given by 0 = R1 + 3 ∈ ℂ^7879:7S (4) 202401191 12 where R ∈ ℂ^7879:7; effective channel matrix as described above, 1 ∈ ℂ^7;:7S is the transmit signal matrix collecting the symbols from all NU UEs, each drawn from the constellation ^ of cardinality N ^, and 3 ∈ ℂ^7879:7S is the receiver- side additive white Gaussian noise (AWGN) matrix with independent and identically distributed (i.i.d.) elements drawn from ^ ^(0, N₀), where N₀ is the noise variance. The transmit signal X comprises of a pilot block 1_V ∈ ℂ^7;:7W, such that where NP and ND denote the number of symbol slots allocated to the pilot and data sequences, respectively, with NT = NP + ND, and where the pilot symbol matrix XP is assumed to be perfectly known at the CPU. In view of the equations provided above, the goal of the ISAC method presented hereafter can be concisely stated. The communication objective of the CPU is to estimate the unknown data symbol matrix XD, under the knowledge of only the pilot symbols in XP, after the estimation of the channel matrix G. In turn, the sensing objective is to extract the voxelated model of the environment as the vector of occupancy coefficients v, from said channel matrix G. By combining the previously presented channel decomposition model for G, the received signal model for Y and the transmit signal model for X, the overall system model becomes (6) where the unknown variables of interest are the environment (voxel coefficients) vector v and the data symbol matrix XD. In the following description it is assumed that the LOS channel H, and the UE-to- voxel and voxel-to-AP sub-path of components A and B in the foregoing equation for Y are known, which still leaves an atypical relationship between the two variables diag(v) and XD. In particular, the latter unknowns are related, under the foregoing 202401191 13 equation for Y, by an asymmetric bilinear system, requiring sophisticated computation methods to be either decoupled or jointly estimated. In light of the above, hereafter an ISAC method for the joint estimation problem of the asymmetric bilinear system expressed by the foregoing equation for Y is proposed, which leverages the well-known Gaussian belief propagation (GaBP) message passing (MP) framework. The proposed method, dubbed Alternating Linear ISAC (AL-ISAC), incorporates two separate linear estimation modules for each of the unknown variables v and XD, estimating these unknown variables in an alternate fashion via a feedback chain between the two modules. Both estimation modules are based on equation (6) and can respectively described as: - an iterative linear GaBP MP method for estimating the environment vector v, given the transmit signal matrix X, and conversely, and - an iterative linear GaBP MP method for estimating the transmit signal matrix X, given the environment vector v. The two linear GaBP methods and the constituting MP rules are derived hereinafter, before the construction of the full AL-ISAC method encompassing the two derived methods is described. The first iterative linear GaBP MP method for the environment vector v operates on only one unknown variable, so that in order to estimate _v, the entire transmit signal matrix X must be assumed known, in addition to the known channel matrices H, A, and B. Assuming knowledge of X, the system described by the foregoing equation for Y may be reformulated as (7) where, since the channel matrix \ ∈ ℂ^{7879:7^} and matrix products .1 ∈ ℂ^7879:7S and ]1 ∈ ℂ^{7^:7S} are known, the described system is linear on v, to which a corresponding factor graph as shown in figure 8 may be obtained. 202401191 14 Each element ym,t of the receive signal of Y, with m ^ {1, …, NANR} and t ^ {1, …, NT} corresponds to the factor nodes, represented by the squares in the figure, and each element vk of the unknown environment variable v, with k ^ {1, …, NV} corresponds to the variable nodes, represented by the circles in the figure. In turn, each (m,t)-th factor node on the factor graph has a corresponding soft-replica of each variable node element vk, denoted by with the corresponding mean-squared-error (MSE) given by Utilizing soft-replicas and their MSEs, the factor nodes perform soft interference cancellation (IC) on the received signal ym,t for each variable vk, yielding the IC symbol where xn,t and wm,t are the (n,t)-th and (m,t)-th elements of X and W, respectively, with n ^ {1, …, NU}, and where the auxiliary variable ck,t ≜ represents the aggregated incident signal at the k-th voxel from all UEs. Next, by leveraging the central limit theorem (CLT), the sum of difference and noise terms are approximated as a complex Gaussian scalar, such that the PDF of the interference-cancelled symbols qr_b ^g :_d,e can be modelled as with the corresponding variance s_b ^t :_d,e given by 202401191 15 The conditional variances for all vk are computed by all factor nodes, and the message is sent to the corresponding variable nodes. Consequently, the k-th variable node obtains the NANRNT conditional variances from all factor nodes, from which the extrinsic belief is computed. In GaBP, self-interference is suppressed by excluding the conditional PDF at the k-th variable node to yield with the PDF where the extrinsic mean and variance is respectively given by Finally, by following the Bayes rule, the updated posterior may be obtained by combining the PDF of the extrinsic belief and the prior distribution of vk, from which the updated soft-replica is obtained as where the normalizing factor in the denominator is the integrated updated posterior over the complex field. Similarly, the updated error variance of the soft-replica is obtained by evaluating Given the information of the voxel coefficient distributions, i.e., binary coefficients with a discrete prior given by a Bernoulli distribution with occupancy probability 202401191 16 |_g ≜ ℙ_g, ⁽`<sub>b</sub> = 1<sup>)</sup>, the soft-replica and its MSE can be efficiently obtained in closed form, respectively, given by The updated soft-replica and the MSE of each variable node are then transmitted back to all factor nodes for the next iteration of the GaBP MP computation method. After a given number of GaBP iterations to refine the soft-estimates, a belief consensus is taken at each variable node across the soft-replicas to obtain a single estimate ṽ by with consensus mean {}<sub>b</sub> <sup>g</sup> and variance expressed as which is consequently used to yield the final estimate by Equations (8) to (20) fully describe the linear GaBP MP method for estimating the voxel environment v, i.e. the voxel coefficients, which is summarized as follows, with reference to the steps of the method 200 shown in figure 9. Inputs to the method are the received signal matrix Y, the channel matrices H, A, and B, the transmit signal matrix X, the noise variance N0, and the prior distribution of the environment voxels ℙ<sub>^,</sub> <sup>(</sup>`_b ⁾. The iterative process begins with an initialisation, since the soft-replicas of a`<sub>b:d,e</sub> and the corresponding MSEs from earlier iterations are not yet available. After computing, in step 202, a first signal matrix C as per C ≜ BX, which describes the effective signals as modified by the channel between the UE and a voxel, in step 202401191 17 204, the soft-replica at all variables is initialised as = ^<sub>g,</sub> <sup>[</sup>`_b ^]∀^, ^, ^, and in step 206 the MSEs at all variable nodes are initialized ∀^, ^, ^ as per equation (8). The initialisation is executed ∀^. The following steps are iterated until a termination criterion is met. In step 210 the soft-IC is performed on the received signal ym,t using equation (9), yielding qr_b ^g :_d,e . In step 212 the corresponding conditional variance s_b ^t :_d,e is determined as per equation (11). In step 214 the extrinsic mean {_b ^g :_d,e and variance is determined as per equation (13). In step 216 the new soft-replica a`<sub>b:d,e</sub> and the MSE f<sub>b</sub> <sup>g</sup> :<sub>d,e</sub> are computed as per equations (14) and (15), respectively, and in step 218 the soft- replica a`_b:d,e and the MSE f_b ^g :_d,e are updated via damping, e.g., as discussed by P. Som, T. Datta, A. Chockalingam and B. S. Rajan, in “Improved large-MIMO detection based on damped belief propagation,” IEEE Inf. Theory Workshop on Inf. Theory, 2010, pp.1–5. The damping factor may be selected from an interval η ^ [0, 1] and serves for preventing early convergence to a local optimum. It is to be noted that each of steps 210 through 218 is carried out ∀^, ^, ^. In step 220 a check is performed to find out if a termination criterion is met. In the negative case, “no”- branch of step 220, steps 210 through 218 are repeated, using the respective most recent values of the soft-replicas and the MSEs as input. In the positive case, “yes”- branch of step 220, the consensus mean {}_b ^g and variance ^^/ _b ^g are computed ∀^ as per equation (19) in step 222, and in step 224 the final soft-estimate ṽ is computed as per equation (20). The termination criterion for the iteration loop may comprise a predetermined maximum number of iterations, or a predetermined convergence threshold which may be evaluated in terms of the soft-replica MSE. Note that the steps following meeting the termination criterion are carried out ∀^. It is important to note that the signal matrix X is assumed given – i.e., X is not estimated by the method – and therefore is kept constant throughout the iterations, as seen by the precomputation of the effective signals ck,t. Next, the complementary second linear GaBP MP method dedicated to the estimation of the transmit signal matrix X for a given v will be described. 202401191 18 In order to derive the linear GaBP MP method for estimating the signal matrix X for a given v, the overall system model derived further above and represented by equation (6) is first reduced to the aggregated received signal matrix Y initially presented in equation (4), with the effective channel G ≜ H+A diag(v)B. The corresponding factor graph is illustrated in Fig.10, where each element xn,t of the unknown signal matrix X with n ∈{1, …, N_U }, is a variable node, represented by the circles. Note that in this case, since the variable X is two-dimensional, the linear system results in a factor graph that is separated into “pages”, such that the variable nodes and factor nodes having different time indices t ∈{1, …, N_T } are independent and that the messages are only exchanged by nodes with the same time index t. Other than this separation of factor graphs, the derivation of the MP rules is similar to that of the linear GaBP MP method for the environment vector v presented further above. The soft-replica of the transmit signal matrix element xn,t to the (m,t)-th factor node is denoted by with the corresponding MSE given by The soft-replica and the MSE are used in the soft-IC of the received signals for xn,t, following and the conditional variance is given by The conditional PDFs are combined with self-interference cancellation at the variable nodes to yield the extrinsic beliefs ℓ_^ ^{^} _,e:d,e following 202401191 19 with the extrinsic mean In turn, since the symbols have a uniformly discrete prior from the symbol constellation χ , with the symbol probability ℙx(x) = 1/| χ |, their soft-replicas and MSEs are obtained by For the particular case of M-ary quadrature amplitude modulation (M -QAM) with M = 4, the soft-replica and MSE computations reduce to efficient closed-form expressions given by where E ^ ≜ ^x [|x|²] denotes the average symbol power of the constellation ^, and tanh(·) denotes the trigonometric hyperbolic tangent function. Finally, the consensus PDF, which is taken after the iterations is given by (31) with the consensus mean {}_^ ^{^} ,_e and variance expressed as 202401191 20 yielding the soft estimate Equations (21) to (33) fully describe the linear GaBP MP method for estimating the transmit signal matrix X given the environment vector v, which is summarized below, with reference to the steps of the method 300 shown in figure 11. Inputs to the method are the received signal matrix Y, the channel matrices H, A, and B, the environment vector v, the noise variance N0, and the prior distribution of the transmit symbols Like computation method 1, computation method 2 begins with an initialisation, in which the effective channel matrix G≜ H + A diag(v) B is computed in step 302, followed by initialising the soft-replica at all variable nodes in step 304 as = ^_^^,y^o_^,e^, and by initialising, in step 306, the MSE at all variable nodes as per equation (21). The initialisation is executed ∀^, ^, ^. The following steps are iterated until a termination criterion is met. In step 310 the received signal after soft-IC is computed as per equation (22), and the corresponding conditional variance is computed in step 312, using equation (24). Next, in step 314 the extrinsic mean and variance are computed as per equation (26), which permits computing, in step 316, the new soft-replica as per equations (27) and (28), respectively. Now the soft-replica oa_^,e:d,e and the can be updated via damping in step 318, similar to the updating process in computation method 1 discussed further above. Here too, damping serves for preventing early convergence to a local optimum. It is to be noted that each of steps 310 through 318 is carried out ∀^, ^, ^. In step 320 a check is performed to find out if a termination criterion is met. In the negative case, “no”-branch of step 320, steps 310 through 318 are repeated, using the respective most recent values of the soft-replicas and the MSEs as input. In the positive case, “yes”-branch of step 320, the consensus mean {}_b ^g and variance ^^/ _b ^g are computed ∀^ as per equation (32) in step 322, and in step 324 the final soft- 202401191 21 estimate o}_^,e is computed as per equation (33). The final soft-estimate o}_^,e is then projected, in step 326, to the symbol constellation ^, and the projected o}_^,e is output, in step 328, as final hard estimate. The termination criterion for the iteration loop may comprise a predetermined maximum number of iterations, or a predetermined convergence threshold which may be evaluated in terms of the soft-replica MSE. c The previously described first and second methods now need to be combined for completing the proposed full AL-ISAC method, which estimates both unknown variables. As mentioned before, each of the first and second estimation methods presented before is adapted to estimate only one of the two variables v or X, assuming in each case availability of full information on the respective other variable. However, the very inherent problem of ISAC in the overall system model derived further above is that neither of the variables are fully known such that the linear methods may not be directly applied for estimation. To address the problem, the proposed AL-ISAC method successively invokes the two linear GaBP MP methods for estimating the two sets of variables. This requires separating the received signal is into blocks corresponding to the pilot phase and the data phase, as 0_T = (. + \diag 0_V = (. + \diag where Y ≜ [YPYD] and W ≜ [WPWD], as defined with X ≜ [XPXD]. First, using only the pilot phase of the system as represented by equation (34a), the first linear GaBP MP process for estimating the voxel environment v is invoked for estimating the initial environment vector ^}_^^^e with the pilot block _XP as known input signal matrix. This includes separating the pilot signals Y_P and data signals Y_D contained in the received signals Y. Next, using only the data phase of the system as represented by equation (34b), i.e., using only the received data signals YD, the second linear GaBP MP process for estimating the transmit signal matrix X is 202401191 22 invoked for estimating the unknown data block X^ D using the initial environment estimate ^}_^^^e as the known input environment vector. Finally, the environment vector is obtained by invoking the first linear GaBP MP process once again, for estimating the voxel environment v, but with the initial environment estimate ^}_^^^e as the initialization value of the soft-replicas at all factor nodes, and [XPX^ D], i.e., the known pilot signals XP and the previously estimated data signals 1^/ _V, as the known input signal matrix. The estimated environment vector ṽ and the estimated transmit signal matrix 1^/, comprising the known pilot signals XP and the previously estimated data signals 1^/ _V, are output. A schematic block chart of the proposed AL-ISAC method is illustrated in figure 12. Despite having a potential complexity advantage, especially for scenarios with large numbers of UEs, as will be shown further down, the alternating approach of the AL- ISAC method has the drawback of causing a heavy dependence on the length of the pilot sequence XP, which, as will likewise be shown further below, affects the performance of both environment and data signal estimation, in addition to the obvious trade-off concerning the total communication throughput. In Table I, the complexity order of the proposed method and its constituent modules is given, in terms of the system size parameters NU, NV, NA, NR, NP, NT with ρ ≜ NP/NT denoting the pilot length ratio, and λ denoting the number of iterations. For comparison purposes, the complexity of a Bi-Linear ISAC (Bi-ISAC) as described in German patent publication DE 102022212615 A1, the content of which is hereby incorporated by reference, is likewise provided in the table. Linear GaBP on v (first process), O [λNU(NVNANRNP)²] only pilot phase Linear GaBP on X (second process), O [λ(NUNANRND)²] only data phase 202401191 23 Linear GaBP on v (first process), O [λNU(NVNANRNT)²] pilot and data phase Full AL-ISAC method ^{^ ^} O [λ(NUNVNANR)²H^{7U^7S ?} 7_; + u_VM] = O [λ(NUNVNANRNT)²H(1 j ^)^? + ^{^^^<} 7_; M] Full BI-ISAC method O [λ(NUNVNANRNT)²] (DE 102022212615 A1) Table I: Complexity orders of the proposed estimation modules and methods The comparison of the full order of complexity shows that the Bi-ISAC has a second order complexity dependent on all system size parameters and a linear complexity on the number of iterations. Interestingly, it is found that the AL-ISAC requires the same order of complexity as the Bi-ISAC, except for the scaling factor (1 j ^)^? + ; , which is dependent on NU and ρ. The latter is therefore a measure of the relative complexity of the proposed method and the reference method, which is plotted in figure 13, for various values of NU as a function of ρ. Considering first the case of NU = 1, it is seen that the relative complexity is larger than 1 for all values of pilot length ratio ρ ∈ [0, 1], indicating that in a single-UE case, the AL-ISAC algorithm requires a higher complexity over the Bi-ISAC algorithm, regardless of the amount of pilot symbols. It is also seen, however, that in the multi- UE scenario of NU ≥ 2 both methods have the same complexity order if the pilot 7_;)¤¥7_; ^{^} )7_;)<¤ length ratio ^_¢£ = _7;^< ; and further, that when ^ > ^_¢£, the relative complexity drops below 1, indicating that the AL-ISAC achieves a lower complexity than the Bi-ISAC with increasing NU and sufficient ρ. Figure 13 also indicates that the relative complexity follows a truncated quadratic behaviour with the minimum value occurring at the specific pilot ratio ^_d^^ = which converges to 1 for N → ∞. The result implies that a sufficiently high pilot length ratio is required for the AL-ISAC method to achieve the optimal complexity, especially with increasing NU. In conclusion, the Bi-ISAC method enjoys a complexity that is 202401191 24 robust to the number of UEs and the length of pilots (throughput), whereas the AL- ISAC method has a complexity that decreases with the number of UEs, but at the cost of throughput. In the following section the proposed method, AL-ISAC, and the BI-ISAC method introduced for comparison purpose in the preceding section, are compared against the method presented by X. Tong, Z. Zhang, J. Wang, C. Huang, and M. Debbah, in “Joint multi-user communication and sensing exploiting both signal and environment sparsity,” IEEE J. Sel. Topics Signal Process., vol.15, no.6, pp.1409–1422, 2021, which also utilizes a voxelated grid to perform ISAC via leveraging multiple linear MP algorithms. It is noted, however, that this known method utilizes the support of a fully known intelligent reflective surface (IRS) in the ROI, and strongly relies on sparsity in the received signal, offered by means of an SCMA interface between the UEs and the AP. Contrary to that, the method proposed herein as well as the BI-ISAC method are not limited to such specific conditions. Rather, the proposed method as well as the BI-ISAC method operate over fully dense received signals and do not need support of IRSs. Due to this distinction in system set-up, an adequate system parametrization must be considered for a fair comparison of the method proposed herein against the known method. Specifically, in the SCMA scheme employed in the known method, each single-antenna UE transmits an M-bit code by utilizing df out of R orthogonal frequency bands, which is received by a single AP with NR receive antennas. Since the considered system for the methods proposed herein is a single-frequency model as represented in equation (4), the diversity gain between each UE and the CPU is mimicked by setting the number of APs as NA = df , such that NANR = df ·NR, i.e., the number of nodes and edges of the final factor graph is the same in both systems. Figure 14 illustrates the sensing and communication performances of the ISAC method proposed herein in comparison to the BI-ISAC method as well as the known ISAC algorithm in terms of the minimum square error (MSE) and symbol error rate (SER), respectively. First, in Fig.14a, the sensing performance, i.e., the MSE of the environment estimation in the form of voxel occupancy coefficients, is evaluated. Under equivalent system parameters, in particular with a pilot ratio of ρ = 0.5 (a value 202401191 25 taken from the paper discussing the known ISAC method for enabling direct comparison), the MSE of the Bi-ISAC method is found to significantly outperform that of the known method at all signal-to-noise ratio (SNR) values, while the AL-ISAC is found to be slightly outperformed by the latter at the high SNRs regime. In Fig.14b, the communication performances of the ISAC systems are evaluated in terms of the SER of the estimated symbols. It can be seen that the method proposed herein as well as the BI-ISAC method exhibit a superior symbol estimation performance compared to the known method at ρ = 0.5, with the Bi-ISAC exhibiting the additionally desirable feature that no error floors are observed even at higher SNRs. All in all, the results corroborate the initial claim that the method proposed herein generally outperform the known reference method in both, sensing and communication functionalities. In view of the excellent performance of the ISAC method proposed herein in comparison with the reference methods, in the following section the robustness of the proposed method is analysed. While the MSE of the estimated voxel coefficients and the SER of the estimated communication symbols were used for evaluating the performance for sensing and communication functions and for comparing with the known methods, in the following section different performance metrics are used. For the sensing function in particular, it is noted that metrics used for radar-based ISAC cannot be used directly due to the unique voxelated occupancy grid-based approach exploited herein. It is, therefore, useful to instead introduce a new metric, referred to as the voxel-occupancy-error-rate (VOER), which measures the rate of false-positive (FP) and false-negative (FN) elements, defined as the incorrect estimation of an occupied voxel element in presence of an empty ground-truth, and the incorrect estimation of an empty voxel element in the presence of an occupied ground-truth, respectively. Mathematically, the VOER is defined as VOER ≜ ^^[‖^ j ^}^‖ _v ^] ∕ u_®, where v is the ground truth and ^} is the voxel coefficient estimate vector, respectively, while ‖·‖_v denotes the ℓ0-norm of a vector. 202401191 26 Note that for the trivial all-empty (or “blind”) estimator, which returns ^} = 0_7^:<, the VOER reduces to |_g ≜ ^^[‖^^‖ _v ^]⁄ u_® = VOER_¢d°e±, which is the average sparsity of the environment. This figure can, therefore, be used as an absolute reference of performance, in the sense ²³|´ ≪ VOER_¢d°e± indicating a good sensing performance for a considered ISAC method. Finally, instead of the SER often used in related literature, here the communication performance of the proposed ISAC method in terms of the more descriptive bit error rate (BER) is evaluated, which is defined where Be denotes the number of erroneously detected data bits of XD, and B is the total number of bits conveyed in _XD. In the first evaluation, the effect of varying the amount of pilot symbols - as captured by the parameter ρ - onto the performance of the method proposed herein is studied. The results on the sensing performance, shown in Fig.15a, indicate that, with the same values of ρ, the Bi-ISAC method outperforms the AL-ISAC method over the entire SNR range. In other words, for a given SNR the Bi-ISAC method requires a much lower number of pilots to achieve the same VOER performance as the AL-ISAC method. The effect of the pilot ratio on the communications performance in terms of BER is evaluated in figure 15b. It can be seen that tor small pilot ratios, i.e., ρ < 0.2, the AL- ISAC method proposed herein achieves a similar performance as the reference Bi- ISAC method, but that for larger pilot ratios, i.e., ρ > 0.3, the AL-ISAC method outperforms the Bi-ISAC method in moderate SNR cases, i.e., 5dB. This result, which may appear to be counterintuitive, is actually expected and explained by the fact that linear MP methods of the AL-ISAC method are constructed on the assumption of perfect symbol knowledge, i.e., 0 uncertainty for the symbol estimates, which assumption is increasingly met with large pilot ratios. Ultimately, however, at significantly high SNRs, e.g., SNR ≥ 15dB, the Bi-ISAC method is again shown to outperform the AL-ISAC method, which is a direct consequence of the error-floor behaviour exhibited by the AL-ISAC method, as already observed in figure 14b. 202401191 27 Finally, figure 16 elucidates the effect of random channel blockages. In particular, the figure compares the VOER and BER performances of the proposed ISAC method and the reference Bi-ISAC method with respect to the critical angle θ *, which determines the channel blockage rate following the stochastic-geometric empirical model derived further above (see figure 6). The environment sensing performance illustrated in figure 16a exhibits a similar behaviour with respect to the effect of ρ, with the Bi-ISAC method achieving a superior performance for all cases. In addition, the curves of the Bi-ISAC method show a slower increase in gradient compared to those of the AL-ISAC method, which indicates the higher robustness of the Bi-ISAC method to path blockages. The superior robustness of the Bi-ISAC method at short pilot lengths and against random channel blockages can be accredited to the increased number of edges in the full factor graph arising from the bilinear representation of the system. Each factor node of the bilinear factor graph is connected to a significantly larger number of variable nodes compared to the factor nodes of the linear factor graphs, which implies more remaining edges for stable message passing even when a large number of edges are removed. Message passing over the pruned graph is only feasible when there is sufficient pilot data still connected to the main graph, therefore making the Bi-ISAC scheme more dependent on the pilot ratio for stability. As for the communications performances, compared in Fig.16b, it is found that the behaviour of both methods differs with the pilot ratio. For low pilot ratios, the AL-ISAC method is shown to slightly outperform the Bi-ISAC method, with similar robustness to the path blockages, whereas for high pilot ratios, the Bi-ISAC method exhibits a significantly superior performance over the AL-ISAC method. In accordance with a first aspect of the present invention a method of processing wireless communication signals for use in joint communication and environment perception in a region of interest (ROI) is proposed. The ROI, which is represented by a voxelated environment (VE) comprising voxels arranged in a three-dimensional grid, comprises N_A ≥ 1 access points and N_U ≥ 1 UEs. All APs serving the ROI are 202401191 28 communicatively connected to a central processing unit (CPU), and each of the NA access points has NR ≥ 1 antennas. The method comprises receiving, at the CPU, a signal matrix Y representing NT ≥ 1 transmission instances received at NR antennas respectively associated with the NA APs via LOS paths and NLOS paths of the VE. The N_T transmission instances carry a plurality of transmit symbols x_n,t comprising pilot signals and data signals, sent by all of the N_U UEs. The method further comprises receiving, prior channel matrices for the LOS paths and the NLOS paths, a prior distribution ℙ_^, ⁽`<sub>b</sub> <sup>)</sup> of voxels of the VE, a noise variance N0, a matrix XP representing the pilot signals used in the transmissions and a prior distribution transmit symbols. The method yet further comprises separating pilot signals YP and data signals YD comprised in the received signals Y, and subjecting only the received pilot signals YP to a first iterative process of estimating an environment vector v representing the VE of the ROI from known transmit symbols X<sub>P</sub>, for obtaining an initial estimated environment vector ṽinit. Only the received data signals Y<sub>D</sub> are then subjected to a second iterative process of estimating transmit data signals 1<sup>/</sup> <sub>V</sub> in a known environment, using the initial vector ṽinit previously estimated as known environment input. Next, the pilot signals X<sub>P</sub> and the previously estimated data signals 1<sup>/</sup> <sub>V</sub> are subjected again to the first iterative process of estimating a vector v representing the VE of the ROI, using the initial environment vector ṽinit previously determined as initialization vector. Finally, the estimated environment vector ṽ and the estimated transmit signal matrix 1<sup>/</sup> are output. In one or more embodiments of the method the prior distribution ℙ<sub>^,</sub>(`_b) of voxels in the VE is provided from a prior execution of the method, from an empirical stochastic process based on the locations of APs and UEs in the VE and the geometric LOS and NLOS paths between them, or is a random distribution. The empirical stochastic process may comprise a process as described in International patent application no. PCT/EP2024/056582, the content of which is hereby incorporated by reference. The prior distribution ℙ_^, ⁽`<sub>b</sub> <sup>)</sup> provides some information about the scattering or blocking properties of the respective voxel and determines, inter alia, the channel matrices for the NLOS paths. 202401191 29 In one or more embodiments of the method the prior channel matrices for the LOS and NLOS paths are provided from a prior execution of the method or through a geometrical analysis of the VE based on the locations of APs and UEs in the VE, in the latter alternative at least the NLOS paths optionally being punctured in accordance with the prior distribution ℙ<sub>^,</sub> <sup>(</sup>`_b ⁾ of voxels in the VE or with a random distribution. In one or more embodiments of the method the prior distribution of the transmit symbols is provided from a statistical analysis of prior executions of the method or is a random or pseudo-random distribution of the transmit symbols comprised in the constellation of possible transmit symbols. In one or more embodiments of the method the first iterative process of estimating an environment vector v representing the VE of the ROI from known transmit symbols comprises a first linear GaBP process. The first linear GaBP process comprises, after an initialisation, computing a received signal after soft interference cancellation, computing a conditional variance for the received signal, computing an extrinsic mean and variance for each voxel based on the conditional variance, computing a new soft replica and corresponding error for each voxel based on the extrinsic mean and variance and a distribution for the respective voxel from a prior iteration, and updating the soft replica and the error via damping. The computations and updating are repeated until a termination criterion is met. After the termination criterion is met, a consensus mean and variance for each voxel are computed, and ultimately the final soft estimate for each voxel is computed. In one or more embodiments of the method the second iterative process of estimating transmit data signals 1^/ _V in a known environment a second linear GaBP process. The second linear GaBP process comprises, after an initialisation computing a received signal after soft interference cancellation, computing a conditional variance for the received signal, computing an extrinsic mean and variance for the signal based on the conditional variance, computing a new soft replica and corresponding error for the signal based on the extrinsic mean and variance and a distribution for the respective voxel from a prior iteration, and updating the soft replica and the error via damping. The computations and updating are 202401191 30 repeated until a termination criterion is met. After the termination criterion is met, a consensus mean and variance for the signal and the final soft estimate for the signal are computed. The soft estimate is projected to the symbol constellation, and the projected symbol is output as the hard estimate. The termination criteria of the first and/or the second iterative process can include, for example, a predetermined numerical iteration limit, or a convergence of the respective estimate within a predetermined range or below a predetermined value. Such convergence criterion can be fulfilled, e.g., when the average change between consecutive post-iteration estimates is below the predetermined value. In accordance with a second aspect of the present invention a system for wireless communication is presented. The system comprises a CPU and one or more APs communicatively connected to the CPU. The CPU comprises a microprocessor, volatile and non-volatile memory and a communication interface for communicating with one or more APs. The elements or components of the CPU are communicatively connected via one or more data or signal lines or buses. Each of the APs has at least one antenna, circuitry for processing radio frequency signals, a microprocessor, volatile and non-volatile memory and a communication interface for communicating with the CPU. The elements or components of the respective AP are connected via one or more data and/or signal lines or buses. The non-volatile memory of the each AP stores computer program instructions which, when executed by the microprocessor, causes the respective AP to transmit received transmission instances to the CPU. The non-volatile memory of the CPU stores computer program instructions which, when executed by the microprocessor, configure components of the CPU to implement or carry out a method in accordance with the first aspect of the invention. In one or more embodiments of the system the circuitry for processing radio frequency signals of the APs comprises a low noise amplifier and/or a mixer configured for providing a representation of a received signal at an intermediate frequency. 202401191 31 The methods described hereinbefore may be represented by computer program instructions. Thus, in accordance with a third aspect of the present invention a computer program product comprises computer program instructions which, when executed by a microprocessor of or functionally coupled with a CPU of a system in accordance with the second aspect of the invention, cause the processor and/or the CPU to carry out a method in accordance with the first aspect of the invention and/or, when executed by a microprocessor of or functionally coupled with an AP of the system in accordance with the second aspect of the invention, cause the processor and/or the AP to transmit received transmission instances to the CPU. The computer program instructions may be retrievably stored or transmitted on a computer-readable medium or data carrier. The medium or the data carrier may by physically embodied, e.g., in the form of a hard disk, solid state disk, flash memory device or the like. However, the medium or the data carrier may also comprise a modulated electro-magnetic, electrical, or optical signal that is received by the computer by means of a corresponding receiver, and that is transferred to and stored in a memory of the computer. The method presented herein represents a new JCAS or ISAC method in which a voxelated 3D representation of a ROI is extracted from the scattering features present in the effective CSI, utilizing the same physical layer communications air interface of an uplink connection between multiple UEs, acting as transmitters, and one or more APs, acting as receivers. The ISAC method, dubbed AL-ISAC, relies on a modular feedback structure in which the transmit data and the environment are estimated alternatingly. Computer simulations of the method show to outperform known methods in accurately recovering the transmitted data, as well as in obtaining a voxelated 3D image of the environment. An analysis of the computational complexities of the proposed method reveals distinct advantages of the scheme, namely, that AL-ISAC offers lower complexity especially in scenarios with large numbers of UEs. The method presented herein may be used for providing joint communication and environment perception in scenarios with multiple independent users and multiple cooperating receivers, e.g., fully connected and automated factories, warehouses, 202401191 32 etc., with centralized processing, such as industrial edge computing. In particular, the present invention can be used in scenarios with stationary APs, communicating and detecting an environment by mobile vehicles communicating to roadside units (RSUs) while achieving vehicular/pedestrian detection, by multiple vehicles cooperatively sensing an environment and road conditions without RSUs, by multiple connected UEs (Bluetooth, Wi-Fi, IoT, etc.) for passively sensing an environment (i.e., without the use of sensing specific signals), and the like. The use-cases and benefits of the methods of communication and environment sensing functionalities presented herein involve and can be achieved within, respectively, the same wireless interface. The present invention advantageously removes the limitation of the communication and access scheme to specific transmission schemes found in known methods, such as the SCMA transmission scheme in the known method discussed in the background section, thereby inter alia dispensing with the requirement of multiple frequency subcarriers in deployment and, thus, permitting the robust application of JCAS in various situations. Further, the present invention lifts the confinement imposed on the detection sliding window length found in prior art methods, which is determined by the pilot length. Yet further, the present invention provides a system model that is no longer limited to a single AP, and single antenna UEs, permitting exploitation of larger receive and transmit diversity as well as multiple-input and multiple-output (MIMO) techniques, which refers to a practical technique for sending and receiving more than one data signal simultaneously over the same radio channel by exploiting multipath propagation. In addition, the present invention removes the dependency on a single RIS, which dependency restricts some known methods to specific scenarios in which such single RIS is available. BRIEF DESCRIPTION OF THE DRAWING 202401191 33 The figures in the attached drawing are used for detailing aspects of the present invention. In the figures Fig.1 shows an exemplary environment with objects in a region of interest, and voxelated representations thereof, Fig.2 shows a representation of the general concept of LOS and NLOS paths in a voxelated space, Fig.3 shows a schematic representation of the 3D space considered in a known system and method, Fig.4 shows process modules of a prior art JCAS method that are iterated in a sliding window fashion, Fig.5 shows a schematic representation of unavailable and unrealistic communication paths in the voxelated space, Fig.6 shows empirical beta mixture modelling of scattering angle distributions, Fig.7 shows the effect of the critical angle θcrit on the severity of puncturing on the channel matrices for different voxelated grid resolutions Fig.8 shows the factor graph of the linear system formulated for the estimation of the voxel coefficients v, Fig.9 shows an exemplary flow diagram of method steps invoked for the estimation of the voxel coefficients v, Fig.10 shows the factor graph of the linear system formulated for the estimation of the signal matrix X, Fig.11 shows an exemplary flow diagram of method steps invoked for the estimation of the signal matrix X, Fig.12 shows a schematic flow diagram of the proposed AL-ISAC method, Fig.13 shows the relative order of complexity of the proposed method for varying values of ρ and NU, Fig.14 shows the MSE and SER performances of the known method and the and Bi-ISAC method against the proposed AL-ISAC method, Fig.15 shows the performance of the proposed AL-ISAC method in comparison to the Bi-ISAC method in a system with parameters: N_U = 4, N_RN_R = 12, N_T = 100, NV = 512, Ev = 1,5%, at varying SNR values, as a function of the pilot length ratio ρ, 202401191 34 Fig.16 shows the performance of the proposed AL-ISAC method in comparison to the Bi-ISAC method in a system with parameters: _NU = 4, _NT = 100, _ρ = 0.5, NV = 512; Ev = 1.5%; SNR = 15dB, with random channel path blockage as a function of the critical angle θ*, Fig.17 shows an exemplary block diagram of a system for wireless communication in accordance with the second aspect of the invention, and Fig.18 shows an exemplary flow diagram of the full AL-ISAC method in accordance with the first aspect of the invention. DETAILED DESCRIPTION OF EMBODIMENTS Figures 1 to 16 have been described further above and will not be discussed again. Figure 17 shows an exemplary block diagram of a system for wireless communication in accordance with the second aspect of the invention. The exemplary wireless communication system comprises a central processing unit (CPU) 400 and two access points (AP) 500 communicatively connected to the CPU 400, indicated by the arrows having dash-dotted lines. The CPU 400 comprises a microprocessor 402, volatile 404 and non-volatile memory 406 and a communication interface 408 for communicating with the two APs 500. The elements or components of the CPU 400 are communicatively connected via one or more data or signal lines or buses 410. Each of the APs 500 has at least one antenna 502, circuitry 504 for processing radio frequency signals, a microprocessor 506, volatile 508 and non- volatile memory 510 and a communication interface 512 for communicating with the CPU 400. The elements or components of the respective AP 500 are connected via one or more data and/or signal lines or buses 514. The non-volatile memory 510 of each AP 500 stores computer program instructions which, when executed by the microprocessor 506, cause the APs 500 to transmit transmission instances received via the one or more antennas 502 to the CPU 400. The non-volatile memory 406 of the CPU 400 stores computer program instructions which, when executed by the microprocessor 402, configure components of the CPU 400 to implement or carry out the method in accordance with the first aspect of the invention. Figure 18 shows an exemplary flow diagram of the full AL-ISAC method 100 in accordance with the first aspect of the invention. In step 110 a signal matrix Y 202401191 35 representing NT ≥ 1 transmission instances received at NR antennas respectively associated with the NA APs 500 via LOS paths and NLOS paths of the VE, the NT transmission instances carrying a plurality of transmit symbols xn,t comprising pilot signals and data signals, sent by all of the N_U UEs, prior channel matrices for the LOS paths and the NLOS paths, a prior distribution ℙ_^,(`_b) of voxels of the VE, a noise variance N0, a matrix XP representing the pilot signals used in the transmissions and a prior distribution ℙ_^^,y^o_^,e^ of the transmit symbols are received at the CPU 400. In step 120 pilot signals Y_P and data signals Y_D comprised in the received signals Y are separated. In step 130 only the received pilot signals Y_P are subjected to the first iterative process 200 of estimating an environment vector v representing the VE of the ROI from known transmit symbols XP, for obtaining an initial estimated environment vector ṽinit. In step 140 only the received data signals YD are subjected to the second iterative process 300 of estimating transmit data signals 1^/ _V in a known environment, using the initial vector ṽinit previously estimated as known environment input. In step 150 the pilot signals X_P and the previously estimated data signals 1^/ _V are submitted to the first iterative process 200 of estimating a vector v representing the VE of the ROI, using the initial environment vector ṽinit previously determined as initialization vector. Finally, in step 160 the estimated environment vector ṽ and the estimated transmit signal matrix 1^/ are output.

202401191 36 LIST OF REFERENCE NUMERALS (PART OF THE DESCRIPTION 100 method 400 CPU 110 receive signal matrix 402 microprocessor 120 separate pilot and data signals 404 volatile memory 130 subject pilot signals to first 406 non-volatile memory method 408 comm. I/F 140 subject pilot signals to second 410 signal/data line/bus method 150 subject pilot and data signals to 500 AP first method 502 antenna 160 output estimated environment 504 RF circuitry vector and signal matrix 506 microprocessor 200 first iterative method 508 volatile memory 202 – 510 non-volatile memory 224 method steps 512 comm. I/F 300 second iterative method 514 signal/data line/bus 302 – 328 method steps

Claims

202401191 37 CLAIMS 1. A method (100) of processing wireless communication signals for use in joint communication and environment perception in a region of interest (ROI) comprising NA ≥ 1 access points (AP) (500) and NU ≥ 1 UEs (UE), each of the NA APs (500) having NR ≥ 1 antennas and all APs (500) serving the ROI being communicatively connected to a central processing unit (CPU) (400), the region of interest being represented by a voxelated environment (VE) comprising voxels arranged in a three-dimensional grid, the method comprising: - receiving (110), at the CPU (400), a signal matrix (Y) representing N_T ≥ 1 transmission instances received at N_R antennas respectively associated with the N_A APs (500) via line-of-sight (LOS) paths and non-line- of-sight (NLOS) paths of the VE, the N_T transmission instances carrying a plurality of transmit symbols (xn,t) comprising pilot signals and data signals, sent by all of the NU UEs, prior channel matrices for the LOS paths and the NLOS paths, a prior distribution (ℙ_^, ⁽`<sub>b</sub> <sup>)</sup>) of voxels of the VE, a noise variance (N<sub>0</sub>), a matrix (X<sub>P</sub>) representing the pilot signals used in the transmissions and a prior distribution of the transmit symbols, - separating (120) pilot signals (YP) and data signals (YD) comprised in the received signals (Y), - subjecting (130) only the received pilot signals (Y<sub>P</sub>) to a first iterative process (200) of estimating an environment vector (v) representing the VE of the ROI from known transmit symbols (XP), for obtaining an initial estimated environment vector (ṽinit), - subjecting (140) only the received data signals (YD) to a second iterative process (300) of estimating transmit data signals (1<sup>/</sup> <sub>V</sub>) in a known environment, using the initial vector (ṽinit) previously estimated as known environment input, - subjecting (150) the pilot signals (X<sub>P</sub>) and the previously estimated data signals (1<sup>/</sup> <sub>V</sub>) to the first iterative process (200) of estimating a vector (v) representing the VE of the ROI, using the initial environment vector (ṽinit) previously determined as initialization vector, and 202401191 38 - output (160) the estimated environment vector (ṽ) and the estimated transmit signal matrix (1<sup>/</sup>). 2. The method (100) of claim 1, wherein the prior distribution (ℙ<sub>^,</sub> <sup>(</sup>`_b ⁾) of voxels in the VE is provided from a prior execution of the method (100), from an empirical stochastic process based on the locations of APs (500) and UEs in the VE and the geometric LOS and NLOS paths between them, or is a random distribution. 3. The method (100) of claim 1 or 2, wherein the prior channel matrices for the LOS and NLOS paths are provided from a prior execution of the method (100) or through a geometrical analysis of the VE based on the locations of APs (500) and UEs in the VE, in the latter alternative the prior distribution (ℙ_^, ⁽`<sub>b</sub> <sup>)</sup>) of voxels in the VE optionally being punctured in accordance with the prior distribution (ℙ<sub>^,</sub>(`_b)) of voxels in the VE or with a random distribution. 4. The method (100) of any one or more of claims 1 to 3, wherein the prior distribution transmit symbols is provided from a statistical analysis of prior executions of the method (100) or is a random or pseudo- random distribution of the transmit symbols comprised in the constellation of possible transmit symbols. 5. The method (100) of any one or more of claims 1 to 4, wherein the first iterative process (200) of estimating an environment vector (v) representing the VE of the ROI from known transmit symbols comprises a first linear Gaussian approximation belief propagation (GaBP) process, the first linear GaBP process comprising, after an initialisation: - computing a received signal after soft interference cancellation, - computing a conditional variance for the received signal, - computing an extrinsic mean and variance for each voxel based on the conditional variance, - computing a new soft replica and corresponding error for each voxel based on the extrinsic mean and variance and a distribution for the respective 202401191 39 voxel from a prior iteration, - updating the soft replica and the error via damping, repeating the computations and updating until a termination criterion is met, - computing a consensus mean and variance for each voxel, and - computing the final soft estimate for each voxel. 6. The method (100) of claim 1, wherein the second iterative process (300) of estimating transmit data signals (1^/ _V) in a known environment a second linear Gaussian approximation belief propagation (GaBP) process, the second linear GaBP process comprising, after an initialisation: - computing a received signal after soft interference cancellation, - computing a conditional variance for the received signal, - computing an extrinsic mean and variance for the signal based on the conditional variance, - computing a new soft replica and corresponding error for the signal based on the extrinsic mean and variance and a distribution for the respective voxel from a prior iteration, - updating the soft replica and the error via damping, repeating the computations and updating until a termination criterion is met, - computing a consensus mean and variance for the signal, - computing the final soft estimate for the signal, - projecting the soft estimate to the symbol constellation, and - outputting the projected symbol as the hard estimate. 7. A system for wireless communication comprising a central processing unit (CPU) (400) and one or more access points (AP) (500) communicatively connected to the CPU (400), wherein the CPU (400) comprises a microprocessor (402), volatile (404) and non-volatile memory (406) and a communication interface (408) for communicating with one or more APs (500), the elements or components being communicatively connected via one or more data or signal lines or buses (410), wherein each of the APs (500) has at least one antenna (502), circuitry (504) for processing radio frequency signals, a microprocessor (506), volatile (508) and non-volatile memory (510) and a communication interface (512) for communicating with 202401191 40 the CPU (400), the elements or components of the respective AP (500) being connected via one or more data and/or signal lines or buses (514), wherein the non-volatile memory (510) of each AP (500) stores computer program instructions which, when executed by the microprocessor (506), cause the APs (500) to transmit received transmission instances to the CPU (400), wherein the non-volatile memory (406) of the CPU (400) stores computer program instructions which, when executed by the microprocessor (402), configure components of the CPU (400) to implement or carry out a method of any one of the preceding claims 1 to 6. 8. The system of claim 7, wherein the circuitry (504) for processing radio frequency signals of the APs (500) comprises a low noise amplifier and/or a mixer configured for providing a representation of a received signal at an intermediate frequency. 9. A computer program product comprising computer program instructions, which, when executed by a microprocessor (402) of or functionally coupled with a central processing unit (CPU) (400) of a system in accordance with claim 7 or 8, cause the processor and/or the CPU (400) to carry out a method of any one of claims 1 to 6 and/or, when executed by a microprocessor (506) of or functionally coupled with an access point (AP) (500) of the system in accordance with claim 7 or 8, cause the processor and/or the AP (500) to transmit received transmission instances to the CPU (400). 10. Computer readable medium or data carrier retrievably transmitting or storing the computer program product of claim 9.