Integrated Sensing and Communications for 3D Object Imaging via Bilinear Inference

We consider an uplink integrated sensing and communications (ISAC) scenario where the detection of data symbols from multiple user equipment (UEs) occurs simultaneously with a three-dimensional (3D) estimation of the environment, extracted from the scattering features present in the channel state information (CSI) and utilizing the same physical layer communications air interface, as opposed to radar technologies. By exploiting a discrete (voxelated) representation of the environment, two novel ISAC schemes are derived with purpose-built message passing (MP) rules for the joint estimation of data symbols and status (filled/empty) of the discretized environment. The first relies on a modular feedback structure in which the data symbols and the environment are estimated alternately, whereas the second leverages a bilinear inference framework to estimate both variables concurrently. Both contributed methods are shown via simulations to outperform the state-of-the-art (SotA) in accurately recovering the transmitted data as well as the 3D image of the environment. An analysis of the computational complexities of the proposed methods reveals distinct advantages of each scheme, namely, that the bilinear solution exhibits a superior robustness to short pilots and channel blockages, while the alternating solution offers lower complexity with large number of UEs and superior performance in ideal conditions.

A part of this article has been presented at the 2022 56th Asilomar Conference on Signals, Systems, and Computers, [1].
Within that context, a new research field called integrated sensing and communications (ISAC) [17]- [19], also known as joint communication and sensing (JCAS) [20]- [24], has recently gained significant attention as a promising technology to fulfill such requirements and enable new applications for B5G and 6G systems.In particular, ISAC technology seeks to enhance B5G and 6G systems by enabling both communication and environment sensing functionalities under the same wireless interface, thus realizing the concepts of ambient-sensing and environment-aware radio [25], which are crucial to emerging applications such as autonomous driving (AD) [26] and drone networking (DN) [27], besides offering new means to optimize system performance.
For instance, in B5G and 6G systems operating at high-frequency channels, which are sensitive to path-dependent scattering [28], [29], environment parameters of interest include not only the "usual" channel state information (CSI), but also the positions of users and obstacles that may lead to path blockages.In such systems, ISAC is an alternative to image-based path-blockage prediction approaches, crucial to mitigate the deleterious effects of blockages [30]- [32].
The prominent challenge of ISAC arises from the fact that two independently-developed wireless technologies, namely, wireless communications and radar systems [33], [34], are both fundamentally based on distinct air-interfaces, such that a concurrent deployment of existing waveforms is bound to suffer from performance degradation of both functions due to interference.
Aiming to frontally address this issue, the earliest family of ISAC approaches known as the radar and communication coexistence (RCC) [35], [36], consequently focused on minimizing the interference and maximizing the cooperation between the independently-operated communications and radar subsystems sharing the same frequency spectrum.While the RCC strategy addresses succeeds in managing interference and harmoniously allocating radio resources to operate both subsystems, the approach achieves relatively low spectral efficiencies and does not alleviate hardware costs, whose components remain separate for both subsystems [23], [35], [36].
In light of the fundamental drawbacks of RCC, techniques integrating the communication and radar functions into a single wireless interface have been recently proposed to truly realize the ISAC goal of jointly offering communication and sensing capabilities in a single system.
Literature [22]- [24] classifies such techniques into three types: a) Radar-centric ISAC (RC-ISAC) schemes, which realizes an additional communication functionality over typical radar waveforms, for example by utilizing index modulation (IM) to encode information into multipleinput multiple-output (MIMO) radar signals employing, e.g., the carrier agile phased array radar (CAESAR) waveform [37], or the frequency modulated continuous waveforms (FMCW) [38]; b) Communication-centric ISAC (CC-ISAC) schemes, which realizes an additional environment sensing functionality over standard communication waveforms, typically by extracting the radar parameters such as Doppler-shift and delay from waveforms designed fundamentally for communications functions, such as the IEEE 802.11ad waveform [39], the orthogonal frequency-division multiplex (OFDM) waveform [40], or the orthogonal time frequency space (OTFS) waveform [41]; and c) Dual-function radar communications (DFRC) schemes, which while not having exclusive boundaries with aforementioned RC-ISAC and CC-ISAC categories, are based on waveforms that can be adaptively or jointly optimized between the two functionalities [42], e.g., via waveforms designed based on mutual information [43] or the multi-beam approach in the mmWave bands [44].But still, all these approaches are related by the fact that the target sensing is enabled by the fundamental radar relationship between measurable physical quantities and information on the target [33], [34], that is, target range can be extracted from the delay in the signal, velocity from the Doppler frequency, and bearing from the angle-of-arrival (AoA).
Concomitant with the aforementioned methods, contributions have also been recently made to realize environment sensing capabilities not via target detection with radar, but rather by new standards of environment sensing information such as ambient human activity [45], [46], and three-dimensional (3D) environment image [47]- [49].In particular, the latter family of works exploit the voxelated occupancy grid [50]- [53] to discretize the environment into 3D cubic units of space representing its state (solid or void).Such methods exhibit a unique advantage in that the extracted information not only describes the location of objects, but also their 3D shape and orientation, in any desired resolution according to the voxelated model, enabling useful applications such as 3D environment mapping and ray-tracing propagation modeling [54], [55].
However, wireless voxelated imaging technology is still a very new notion in context of ISAC, due to the inherently convoluted channel paths arising from the discretization of the environment scatterers, and the fundamental challenge that the unknown information symbols must be simultaneously recovered on top of the very large number of environment voxels.
In light of the above challenges, we offer in this article the following contributions: • An extension of the discrete voxelated environment model utilized in [47]- [49] is introduced, in which an empirical stochastic-geometric approach is incorporated to capture the viability of non-line-of-sight (NLOS) paths, in addition to an extension where the occupancy coefficients are not limited to 0 or 1, but instead can take on non-binary complex values, enabling reflection losses and phase shifts of reflected waves to be modeled 1 .
• A novel, scalable CC-ISAC scheme is proposed, in which the 3D voxelated environment image and the transmit symbols are estimated alternately via two distinct linear modules; • A novel, high-performance CC-ISAC scheme is proposed, in which the the 3D voxelated image of the environment and the transmit data symbols are concurrently estimated via a single message passing (MP) module which leverages a bilinear inference framework; • Insights on the advantages of the two proposed CC-ISAC algorithms are provided via performance assessment and comparison against the state-of-the-art (SotA), which highlights the robustness of the bilinear method against short pilot lengths and path blockages, and the accuracy of the alternating solution in systems with many user equipment (UEs).
Notation: Scalar values are denoted by slanted lowercase letters, as in x, while complex vectors and matrices are denoted by boldface lowercase and uppercase letters, as in x and X, respectively.The transposition, complex conjugation, diagonalization, absolute value, and ℓ-th norm operators are denoted by (•) T , (•) * , diag(•), | • |, and || • || ℓ respectively, while E x (x) and Var x (x) respectively denote the expectation and variance of x with respect to its distribution P x (x).The sets of real and complex numbers are denoted by R and C, respectively, and N (µ, ν) and CN (µ, ν) denote the real and complex Normal distributions with mean µ and variance ν.

II. SYSTEM MODEL
The system model considered throughout the article, consists of three parts: a) the environment model, where a voxelated occupancy grid is used to discretely approximate the true environment and its scattering properties, b) the channel model, which defines the unique channel paths arising from the voxelated environment model, and c) the signal model, where the uplink communication scheme between the UEs and the access points (APs) is described.

A. Environment Voxelation Model
The 3D image of an environment can be represented via a number of techniques, including the well-known point-cloud and the ray-tracing methods, which are often utilized in robotics, machine vision and computer graphics [54], [56].Another well-known method, however, is 1 We clarify that although the message-passing rules derived in this article are for this extended paradigm, such that the proposed algorithms are fully generalized, binary-valued occupancy coefficients are considered in Section IV for the purpose of evaluation performance, in order to enable direct comparison with state-of-the-art (SotA) methods.
the voxelated occupancy grid [50]- [53], where the region of interest (ROI) is represented as a cuboidal space of dimensions L x × L y × L z , each denoting the lengths of the x, y, z-axes in meters, respectively, as depicted in Figure 1.In this model, the ROI is subdivided into a regular grid consisting of N V ≜ N x • N y • N z voxels, where N x ≜ Lx L V , N y ≜ Ly L V , and N z ≜ Lz L V denote the number of partitions along the x, y and z-axes, respectively, and L V is the edge length of a voxel in meters, which therefore corresponds to the image resolution.
The environment is then represented by a 3D tensor of dimensions (N x × N y × N z ), whose elements indicate the occupancy of the voxels, and thus whether that portion of the space is empty or filled with a given material.One may therefore consider, in general, each k-th voxel to be represented by an occupancy (a.k.a.scattering) coefficient v k , with k ∈ {1, • • • , N V }, where v k = 0 indicates that the k-th voxel is empty, while an occupied voxel is indicated by a complex number i.e., v k ≜ β k • e −jω k ∈ C. In such a model, the constants β k and ω k , which dependent not only on the material itself, but also on the frequency and the angle of incidence of propagating signals [29], capture the effect of the material occupying a given voxel onto the electromagnetic wave reflected or refracted by it.For the sake of reducing the complexity of the ISAC algorithm to be later introduced, however, we will in this article consider a simplified model whereby the occupancy coefficients v k take on discrete real values in the interval ∈ [0, 1].
Since phase rotations due to reflections are captured by channel estimation, this simplification is equivalent to the assumptions that the ISAC waveform is narrowband, so that frequencydependence can be ignored [29].The incorporation of the geometry of the interaction between propagating waves and occupied voxels will be described in Subsection II-C.

B. Channel Model
Consider a scenario in which the ROI contains N U single-antenna UEs, and N A multi-antenna APs equipped with N R receive antennas each.As illustrated in Fig. 2, the effective channels between the UEs and APs contain two distinct types of components, namely, line-of-sight (LOS) components which are direct paths between the UEs and APs, and NLOS components which encompass paths reflected at occupied voxels corresponding to parts of the environment, as described in Subsection II-A.Assuming that the operating frequency band is sufficiently high that the power of paths reflected more than once is negligible [28], [57], NLOS components may be decomposed into two subpaths, the UE-to-voxel subpath and the voxel-to-AP subpath, which together with the voxel scattering coefficient comprises the aggregate NLOS channel.
In light of the above, the effective channel between the N U single-antenna UEs and the ensemble of N A N R receive antennas of all APs can be compactly described by where and UE-to-voxel NLOS subpath, respectively, whose elements are assumed to follow zero-mean complex Normal distributions with variances σ 2 H , σ 2 A , and σ 2 B , respectively; while v ∈ C N V ×1 is the vector containing all scattering coefficients of the voxelated grid.
Although not further exploited in this article, we emphasize that the channel model in equation ( 1) can be straightforwardly extended to a multi-carrier scenario by simply introducing frequencyselectivity such that the environment variables are functions of the carrier frequency, i.e., with f denoting a specific frequency in the set F of all subcarrier frequencies.
We leave details for follow-up work, but it shall become evident that under such an extended model, the sensing component of the ISAC algorithm to be introduced in Section III can also be extended for even better performance, for instance by exploiting knowledge of channel correlation across carriers [58], and to incorporate additional features, such as the estimation of material types based on the frequency-dependence observed on estimated scattering coefficients [29].Since the above also requires incorporating message-passing rules that exploit cross-carrier correlation to the algorithms, such extension will be pursued in a future work, and only the frequencyindependent (i.e., single-carrier) model of equation ( 1) will be assumed in this article.
subpath Fig. 2: Illustration of the LOS and NLOS channel components and their subpaths.

C. Stochastic Geometric Environment Model
Notice that the channel model summarized by equation (1) implies that all paths between UEs, APs, and voxels are available.Although such an assumption is common in related literature (see e.g.[47]- [53]), in practice, many subpaths may not be available due to either physical phenomena (e.g.blockage by air-borne particles, absorption, or path loss) or the finite resolution of the voxelated model itself.In order to capture such realistic behavior, the work in [49] considers the occlusion effect of waves reaching the voxels, where the UE-to-voxel subpaths are assumed to be unavailable if a voxel is present nearby the path.This perturbation effect was approximated by applying a 3D Gaussian kernel convolution to the channel matrix, but it was acknowledged [49] that this approach is a highly simplified model of the true physical phenomena.Building on the latter, we therefore seek to contribute to improving the voxelated grid model by considering the statistical feasibility of paths.
To this end, we refer to the physical phenomena occurring at the reflection of propagating waves, in particular, the fact that for any given frequency: a) a critical angle θ * exists such that, as illustrated in Fig. 3a, if the incidence angle θ > θ * , the wave is absorbed rather than reflected, and consequently the corresponding voxel-to-AP NLOS subpath is not available [29]; and b) the curvature of the surface exposed to the impinging wave may be such that no signal is reflected towards an AP 2 [50], as illustrated in Fig. 3b.But since the complexity of modeling such phenomena at each voxel is far too complex to carry out, especially if the resolution of the voxelated ROI is large, we instead employ a statistical approach whereby the angle between the 2 Although the phenomenon in Fig. 3b would reduce to the phenomenon in Fig. 3a for infinitely small voxels, such extreme resolution leads to prohibitive complexity of the algorithms, so that modeling both phenomena distinctly is preferred in practice.impinging and reflected waves at each voxel, hereafter referred to as the scattering angle, and consequently, the availability of each voxel-to-AP NLOS subpath, are considered.
In light of the above, the following stochastic-geometric model is proposed to integrate 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 ROI, such that their positions may be described by voxel coordinates 3 Denoting the 3D coordinates of an UE, an AP, and an environment voxel, respectively by the scattering angle θ of the NLOS path reflected at the voxel is given by where arccos(•) denotes the inverse cosine trigonometric function.
The empirical 4 probability distribution functions (PDFs) of the scattering angles θ can be obtained by evaluating equation ( 3) for all possible combinations of admissible locations of UE, AP and voxel within the ROI, respectively given by c U , c A and c V , and examples of the latter for an environment voxelated at various resolutions are shown in Fig. 4. 3 It is also assumed that the multiple antennas of the APs are placed within a single voxel, such that their AoA are assumed to be identical, albeit each with a different channel path coefficient. 4In principle, the analytical distribution of scattering angles in eq. ( 3) can be derived, either by studying the N V 3 constituent angles within the highly subdivided geometry of the grid, or via a stochastic geometry-based grid analysis [59].To the best of the authors' knowledge, however, solutions to this problem exist only for vertex-to-vertex distance distributions [60], with the case of vertex angles never addressed before.Since the focus of this article is to develop ISAC estimators, we leave this matter for a future contribution and meanwhile consider the proposed highly-accurate approximate model, as illustrated in Figure 4.

Probability Distribution of Scattering Angles in a Voxelated Grid
Scattering Angle, 3 It is visible in Figure 4 that for a sufficiently large N V , the distribution of scattering angles θ can be well modeled by a mixture of two scaled beta distributions, namely, Utilizing this empirical stochastic-geometric approach, the increasingly popular voxelated model utilized in various related works [47]- [49], [51]- [53] can be improved by the incorporation of random blockages of NLOS subpaths, 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.

D. Signal Model
Consider an uplink communication scenario between the group of N U UEs and a total of N A APs, under the models described above in Subsections II-A through II-C, with the N A APs connected to a central processing unit (CPU) via error-free fronthaul links of unlimited throughput, such that the received signals at all N A N R receive antennas are aggregated without loss of information or delay.
Then, the aggregated received signal matrix Y, over N T discrete transmission instances (symbol slots) is given by where G ∈ C N A N R ×N U is the effective channel matrix as described in Subsection II-B; X ∈ C N U ×N T is the transmit signal matrix collecting the symbols from all N U UEs, each drawn from the constellation X of cardinality N X ; and W ∈ C N A N R ×N T is the receive additive white Gaussian noise (AWGN) matrix with independent and identically distributed (i.i.d.) elements drawn from CN (0, N 0 ), where N 0 is the noise variance.
The transmit signal X comprises of a pilot block X P ∈ C N U ×N P and a data block where N P and N D denote the number of symbol slots allocated to the pilot and data sequences, respectively, with N T = N P + N D ; and where the pilot symbol matrix X P is assumed to be perfectly known at the CPU.
In view of equations ( 1), ( 4) and ( 5), the goal of the ISAC schemes to be hereafter presented can be concisely stated.The communication objective of the CPU is to estimate the unknown data symbol matrix X D , under the knowledge of only the pilot symbols in X P , 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 the said channel matrix G.

III. PROPOSED ISAC SOLUTION
By combining the channel decomposition model of equation ( 1), the received signal model of equation ( 4), and the transmit signal in model of equation ( 5), the overall system model becomes where the unknown variables of interest are the environment (voxel coefficients) vector v and the data symbol matrix X D .
Similar to related literature [47]- [49], [51]- [53], it is assumed hereafter that the LOS channel H, and the UE-to-voxel and voxel-to-AP subpath of components A and B in equation ( 6) are known, which still leaves an atypical relationship between the two variables diag(v) and X D .In particular, the latter unknowns are related, under equation ( 6), by an asymmetric bilinear system, requiring sophisticated algorithms to be either decoupled or jointly estimated [61]- [66].
In light of the above, we propose in the sequel two novel ISAC solutions for the joint estimation problem of the asymmetric bilinear system expressed by (6), by leveraging the Gaussian belief propagation (GaBP) MP framework.The first proposed method incorporates two separate linear estimation modules for each of the unknown variables v and X D , estimating them in an alternate fashion via feedback between the two modules.In turn, the second method utilizes only a single bilinear estimation module which enables the simultaneous extraction of both unknown variables.

A. Proposed Alternating Linear ISAC Algorithm (AL-ISAC)
The first proposed method, dubbed the "Alternating Linear ISAC (AL-ISAC)" algorithm, leverages two separate linear GaBP MP modules based on equation ( 6), which are respectively described as: 1) a linear GaBP module to estimate the environment vector v, given the transmit signal X, and conversely; and 2) a linear GaBP module to estimate the transmit signal matrix X, given the environment vector v.
The two linear GaBP modules and the constituting MP rules are derived in the next subsections, followed by the construction of the full ISAC algorithm encompassing the two derived modules.

1) Linear GaBP for Environment Vector v:
The linear GaBP algorithm 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 in ( 6) may be reformulated as where, since the channel matrix and BX ∈ C N V ×N T and are known, the described system in ( 7) is linear on v, to which a corresponding factor graph may be obtained as in Fig. 5.
corresponds to the factor nodes (square nodes) and each element v k of the unknown environment variable v, with k ∈ {1, • • • , N V }, correspond to the variable nodes (circular nodes).In turn, each (m, t)-th factor node on the factor graph has a corresponding soft-replica of each variable node element v k denoted by vk:m,t , with the corresponding mean-squared-error (MSE) given by < l a t e x i t s h a 1 _ b a s e 6 4 = " u j H / s y p i j Fig. 5: Factor graph of the linear system formulated for the estimation of v.
Utilizing the soft-replicas and their MSEs, the factor nodes perform soft-interference cancellation (IC) on the received signal y m,t for each variable v k , yielding the IC symbol where x n,t and w m,t are the (n, t)-th and (m, t)-th elements of X and W, with n ∈ {1, • • • , N U }; and c k,t ≜ N U u=1 b k,u x u,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 ȳv k:m,t can be modeled as with the corresponding conditional variance ν v k:m,t given by where is the noise variance.
The conditional variances for all v k are computed by all factor nodes, and the message is sent to the corresponding variable nodes.Consequently, the k-th variable node obtains the N A N R N T conditional variances from all factor nodes, from which the extrinsic belief ℓ v k is computed.In GaBP, self-interference is suppressed by excluding the conditional PDF of ȳk:m,t at the k-th variable node to yield ℓ v k:m,t with the PDF where the extrinsic mean µ v k:m,t and variance Ψ v k:m,t 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 v k , from which the updated soft-replica is obtained as where the normalizing factor in the denominator is the integrated updated posterior over C.
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 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 algorithm.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 μv k and variance Ψ v k expressed as which is consequently used to yield the final estimate by The above MP equations ( 8) to (20) fully describe the linear GaBP module to estimate the voxel environment v, which is summarized as a pseudocode in Algorithm 1.It is important to note that the signal matrix X is assumed given -i.e., X is not estimated by Algorithm 1and therefore is kept constant throughout the iterations, as seen by the pre-computation of the effective signals c k,t .The complementary block dedicated to the estimation of X given v will be described subsequently, leading to an alternate approach as previously mentioned.The algorithm also contains a damping update mechanism [67] using a damping factor η ∈ [0, 1] to prevent early convergence to a local optimum.
In order to derive the linear GaBP module to estimate the signal matrix X given v, the system model in ( 6) is first reduced to the one in (4), with the effective channel G ≜ H+Adiag(v)B , with the corresponding factor graph as illustrated in Fig. 6, where each element x n,t of the unknown signal matrix X with n ∈ {1, • • • , N U }, is the variable node (circular nodes).
Notice 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 corresponding to different 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 Algorithm 1.
. . . . . . . . .< l a t e x i t s h a 1 _ b a s e 6 4 = " 6 s y V i z 6 u L D V A m w S l I U i R p x X s J q 8 = " > A A A B 7 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B g 5 R E i n q S g h e P F e w H t K F s t p t 2 6 e 4 m 7 G 6 x < / l a t e x i t > x n,1 Fig. 6: The separated factor graphs of the linear system for X.
The soft-replica of the transmit signal matrix element x n,t to the (m, t)-th factor node is denoted by xm,t:n,t , with the corresponding MSE given by The soft-replica and the MSE are used in the soft-IC of the received signals for x n,t , following with the conditional PDF described by and the conditional variance ν x n,t:m,t is given by The conditional PDFs are combined with self-interference cancellation at the variable nodes to yield the extrinsic beliefs ℓ x n,t:m,t following with the extrinsic mean µ x n,t:m,t and variance Ψ x n,t:m,t respectively given by µ x n,t:m,t = Ψ x n,t:m,t In turn, since the symbols have a uniformly discrete prior from the symbol constellation X , with the symbol probability 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 ] denotes the average symbol power of the constellation X , and tanh(•) denotes the trigonometric hyperbolic tangent function.
Finally, the consensus PDF, which is taken after the iterations is given by with the consensus mean μx n,t and variance Ψ x n,t expressed as Equations ( 21) to (33), collected in the form of a pseudocode in Algorithm 2, fully describe the linear GaBP module for the estimation of the signal matrix X given the environment vector v, which together with the previously described module for the estimation of v given X, completes the proposed AL-ISAC scheme.All that remains is to describe the alternating procedure to estimate both unknown variables, which is addressed in the sequel.
Algorithm 2 : Linear GaBP Estimator for Signal Matrix X Inputs: Received signal matrix Y, channel matrices H, A, and B, environment vector v, noise variance N 0 , and prior distribution of transmit symbols P xn,t (x n,t ).Outputs: Estimated transmit signal matrix X. ▷ ∀n, t * The termination criteria can be set in accordance to Section IV-C and as described in Algorithm 1.

3) Combined Alternating Modular Structure:
With the two estimation modules given by Algorithm 1 and Algorithm 2, either of the two variables v or X may be estimated, assuming full information of the other variable.However, the very inherent problem of ISAC in equation ( 6), is that neither of the variables are fully known such that the linear modules may not be directly applied for estimation.
To address the problem, the proposed AL-ISAC algorithm successively applies the two linear GaBP modules to estimate the two sets of variables.To enable this, the received signal is separated into the blocks corresponding to the pilot phase and the data phase, as where First, by using only the pilot phase of the system (34a), Algorithm 1 is applied to estimate the initial environment vector ṽinit with the pilot block X P as known input signal matrix.
Next, by using the data phase of the system (34b), Algorithm 2 is applied to estimate the unknown data block XD using the initial environment estimate ṽinit as the known input environment vector.Finally, the environment vector is obtained by using Algorithm 2 again, but with the initial environment estimate ṽinit as the initialization value of the soft-replicas at all factor nodes, and [X P XD ] as the input transmit signal matrix.The described AL-ISAC algorithm is illustrated in a schematic form in Fig. 7, and summarized as pseudocode in Algorithm 3.
Algorithm 3 : Proposed Alternating Linear ISAC (AL-ISAC) Algorithm † Inputs: Received signal matrix Y, channel matrices H, A, and B, pilot matrix X P , noise variance N 0 , prior distribution of environment and transmit symbols P v k (v k ) and P xn,t (x n,t ).Outputs: Estimated environment vector ṽ and estimated data signal matrix XD .
Using only the block corresponding to the pilot sequence, i.e., t = {1, • • • , N P }: 1: Estimate initial environment ṽinit via Algorithm 1, using pilot X P as known signal input.
Using only the block corresponding to the data sequence, i.e., t = {N P + 1, • • • , N P + N D }: 2: Estimate data signal XD via Algorithm 2, using ṽinit as known environment input.
4: Output ṽ as the estimated environment vector, and XD as estimated data signal matrix.† Although our simulations indicate that a single iteration (as shown in Fig. 7) is sufficient, the linear GaBP modules in steps 2 and 3 can be iterated multiple times, with feedback at modular level and possibly adaptive MP denoising between feedback loops.These, and other potential improvements remain open points for a follow up work.

(Alg. 1)
< l a t e x i t s h a _ b a s e = " T H of .
(Alg. 1) < l a t e x i t s h a 1 _ b a s e 6 4 = " T H R H 1 i f P + y i l P w = < / l a t e x i t > ṽ < l a t e x i t s h a 1 _ b a s e 6 4 = " 0 N U A x S 6 q x H 6 7 P E f s 8 a W 2 o e W x p g k = " > < l a t e x i t s h a 1 _ b a s e 6 4 = " r s x H 6 0 Y 5 c 4 h + g P j 8 w c P y 5 W V < / l a t e x i t > XD < l a t e x i t s h a 1 _ b a s e 6 4 = " c z z 7 n k i y j c r B w 4 a t o X n T 7 S 6 d / x H 6 0 Y 5 c 4 h + g P j 8 w c P y 5 W V < / l a t e x i t >

Using only data phase
Using both pilot and data phase < l a t e x i t s h a 1 _ b a s e 6 4 = " c z z 7 n k i y j c r B w 4 a t o X n T 7 S 6 d / Despite having a potential complexity advantage, especially for scenarios with large numbers of UEs, as shown later in Subsection IV-A, the alternating approach of the AL-ISAC algorithm has the drawback of causing a heavy dependence on the length of the pilot sequence X P , which as shall be shown in Section IV, affects the performance of both environment and data signal estimation, in addition to the obvious trade-off with the total communication throughput.Aiming to circumvent this deficiency, in the next subsection we propose another ISAC method in which both v and X D are estimated simultaneously via a bilinear inference method.

B. Proposed Bilinear ISAC Algorithm (Bi-ISAC)
In this section, we develop a new ISAC algorithm in which the sensing and communication variables v and X D are estimated in parallel, by using a bilinear message passing technique which incorporates the uncertainty of both estimates at each iteration, thus requiring only a single estimation module to acquire both variables, as illustrated in Fig. 8.
We start by observing that the unique asymmetric bilinear relationship of v and X D , as per equation (7), prevents the application of recently discovered bilinear estimators, such as the bilinear generalized approximate message passing (BiGAMP) [61], which operates only on symmetric systems described by equations in the form Y = VX + W for the joint estimation of the unknowns V and X; or the parametric BiGAMP [64], [65], which works on systems with In contrast to these two examples, the problem dealt with here is, as described by equation ( 6), in neither of the aforementioned forms, nor can it be transformed to fit general bilinear forms, which implies that new, purpose-built bilinear Gaussian belief propagation (BiGaBP) [62], [63], [66] MP rules must be derived for its solution.Therefore, the BiGaBP message passing is performed on a tripartite factor graph as illustrated in Fig. 9, where the factor nodes (square Joint Est. of and .(Alg.4) q d 0 7 R H 1 i f P + y i l P w = < / l a t e x i t > ṽ < l a t e x i t s h a 1 _ b a s e 6 4 = " 0 N U A x S 6 q x H 6 7 P E f s 8 a W 2 o e W x p g k = " > A C k 3 m 6 X P 8 a l W h j i I p X 4 R 4 J n 6 e y M j o V L T 0 N e T R U i 1 6 B X i f 9 4 g h e D K z X i U p M A i O j 8 U p A J D j I s q 8 J B L R k F M N S F U c p 0 V 0 z G R h I I u r K Z L s B e / v E y 6 5 w 3 7 o t G 8 a 9 Z b z b K O K j p G J + g M 2 e g S t d A t a q M O o u g R P a N X 9 G Y 8 G S / G u / E x H 6 0 Y 5 c 4 h + g P j 8 w c P y 5 W V < / l a t e x i t > XD < l a t e x i t s h a 1 _ b a s e 6 4 = " c z z 7 n k i y j c r B w 4 a t o X n T 7 S 6 d / x < / l a t e x i t > x n,1 < l a t e x i t s h a 1 _ b a s e 6 4 = " u j H / s y p i j X 6 i E I Y A 6 T v h q n B / / i g = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 1 G P R i 8 e K 9 g P a U D b b T b t 0 s w m 7 k 0 I J / Q l e P C j i 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 s y V i z 6 u L D V A m w S l I U i R p x X s J q 8 = " > A A A B 7 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B g 5 R E i n q S g h e P F e w H t K F s t p t 2 6 e 4 m 7 G 6 Fig. 9: The tripartite factor graph of the bilinear system.nodes) are the received symbols, and the two sets of variable nodes (circular nodes) corresponds to the environment vector v and the signal matrix X, respectively.An important distinction is made between the two types of variable nodes, which is that a data variable node receives messages only from N A N R factor nodes corresponding to the same time instance t, while an environment variable node receives messages from all N A N R N T factor nodes.
We highlight the higher complexity of the factor graph in Fig. 9 compared to those of linear GaBP schemes shown in Figs. 5 and 6, due to the asymmetric and embedded system structure of equation ( 6) in relation to both variables together, as opposed to those of linear systems where one variable is considered at a time.Other than that, the messages transferred over the graph edges are the same information as in the linear GaBP in Section III-A, i.e., the soft-replicas, MSEs, and conditional PDFs of the variables of interest.However, since neither of the latter is known, except as soft-replicas, the corresponding calculation of the messages must incorporate the uncertainties in both variables, in the form of the respective MSEs.
In light of the above, similarly to Section III-A, the exchanged messages are constructed on the basis of soft-replicas of the variables.Since the entries in X corresponding to t ∈ {1, • • • , N P } are pilot symbols X P , the corresponding soft-replicas are set to their known values, i.e., xn,t:m,t = (x p ) n,t ∀t ∈ {1, • • • , N P }, with the corresponding MSE values set to 0. The remaining softreplicas and MSEs for t ∈ {N P + 1, • • • , N T } are as given in equations ( 8) and (21).
The complete description of the proposed bilinear ISAC (Bi-ISAC) algorithm is summarized in the form of a pseudocode in Algorithm 4, and the corresponding equations of the message passing rules are elaborated in the following.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, while our contributions are not limited to such a conditions.Instead, the proposed methods operate over fully dense received signals and with no support of IRSs.
Due to this distinction in system set-up, an adequate system parametrization must be considered for a fair comparison of our proposed algorithms against the SotA method of [47].Specifically, in the SCMA scheme of [47], each single-antenna UE transmits an M -bit code by utilizing d f out of R orthogonal frequency bands, which is received by a single AP with N R receive antennas.Since the considered system for our proposed algorithms is a single-frequency model (4), the diversity gain between each UE and the CPU is mimicked by setting the number of APs as N A = d f , such that N A N R = d f • N R , such that the number of nodes and edges of the final factor graph is the same in both systems.Fig. 11 illustrates the sensing and communication performances of the two proposed ISAC algorithms in comparison to the SotA ISAC algorithm of [47], in terms of the MSE and SERs, respectively.First, in Fig. 11a, the sensing performance, i.e., the MSE of the voxel occupancy coefficients, is evaluated.Under equivalent system parameters, in particular with a pilot ratio 7   of ρ = 0.5, the MSE of the proposed Bi-ISAC algorithm is found to significantly outperform that of the SotA 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.The value ρ = 0.5 is taken from [47] in order to enable direct comparison.
Then, in Fig. 11b, the communication performances of the ISAC systems are evaluated in terms of the SER of the estimated symbols.It can be seen that both proposed algorithms exhibit superior symbol estimation performances compared to the SotA, 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 claim that both proposed methods generally outperform the SotA method of [47] in both sensing and communication functionalities.

C. Convergence Behavior of the Proposed Algorithms
In view of the superior performance of the two proposed ISAC algorithms in comparison with the SotA [47], we proceed in this section to further analyze the two proposed ISAC algorithms in detail, aiming also to clarify advantages of each relative to the other.To that extent, we first study the convergence behavior of the proposed ISAC algorithms so as to obtain insight on the appropriate MP termination criteria and damping parameters to be utilized.It can be observed in Fig. 12a that while different values of η v all result in convergence of MSE, the convergent value is sensitive to the parameterization of η v , especially for the Bi-ISAC which suffers from high MSE with low damping, as opposed to the AL-ISAC which is less prone to such errors in both the initial (white circles) and the final (black circles) estimation.
Similar and consequent results in the SER are observed in Fig. 12b, where convergence is also achieved in all cases, but the convergent value is largely affected by the MSE performance and η v , while the effect of η x is not as prominent to the final result compared to η v .
In light of the above results, the parameterization η x = 0.5 and η v = 0.9 and the convergence criteria of λ = 100 iterations is selected in the following simulations to ensure a reliable convergence behavior 8 .

D. Robustness Analysis of the Proposed Algorithms
To that end, we shall utilize performance metrics for the sensing and communication functions of ISAC systems other than the MSE of the estimated voxel coefficients and the SER of the estimated communication symbols, which were used in [47]- [49] and the previous section.
For the sensing function in particular, we remark that metrics used for radar-based ISAC cannot be used directly due to the unique voxelated occupancy grid-based approach followed here.It is therefore sensible to instead introduce a new metric, referred to as the voxel-occupancy-errorrate (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 groundtruth, and the incorrect estimation of an empty voxel element in the presence of an occupied ground-truth, respectively.Mathematically, the VOER is therefore defined as E ||v − ṽ|| 0 /N V , where v is the ground truth, ṽ is the estimate vector, while ∥ • ∥ 0 denotes the ℓ 0 -norm of a vector.
Notice that for the trivial all-empty (or "blind") estimator, which returns ṽ = 0 N V ×1 , the VOER reduces to E v ≜ E[||v|| 0 ]/N V = VOER empty , which is the average sparsity of the environment.We can therefore utilize this figure as an absolute reference of performance, in the sense VOER ≪ VOER empty indicates a good sensing performance by a given ISAC method.
Finally, instead of the SER often used in related literature, we opt to evaluate the communication performance of the proposed ISAC schemes in terms of the more descriptive BER, defined as BER ≜ E[B e ]/B, where B e denotes the number of errorneously detected data bits of X D , and B is the total number of bits conveyed in X D .

Fig. 1 :
Fig. 1: Illustration of a voxelated grid map-based model of the environment of a given ROI.

P3 > 3 $
at h u n av ai la b le AP UE LOS UE-to-AP NLOS voxel-to-AP NLOS UE-to-voxel (a) Path unavailability due to obtuse scattering.

P
at h u n av ai la b le AP R e. ec te d w av e UE LOS UE-to-AP NLOS voxel-to-AP NLOS UE-to-voxel (b) Path unavailability due to skewed surface.

Fig. 3 :
Fig. 3: Illustration of physical phenomena leading to the unavailability of propagation paths.
and where γ is a weighing factor and the quantities a 1 , a 2 , b 1 , b 2 are shape parameters optimised to match the empirical data obtained by evaluating equation (3), with c U , c A and c V taken randomly within the voxelated grid.
j k T w u e n 6 H / S O L b d U 7 t 4 X S x U z u d x Z G E P 9 u E Q X D i D C l x C F e p A I I R 7 e I Q n S 1 o P 1 r P 1 M m v N W P O Z X f g G 6 / U D 9 / 2 Q h g = = < / l a t e x i t > x n,N T < l a t e x i t s h a 1 _ b a s e 6 4 = " V t 7 o h Y h 3 s N H S p A S 3 X f d t L w A G e F s = " > A A A B 7 3 i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g I e F O g l p Y B G w s I 5 g P S I 6 w t 9 l L l u z t n r t 7 Y j j y J 2 w s F L H 1 7 9 j 5 b 9 x L r t D E B w O P 9 2 a Y H T L W 2 j o d j d t l 4 R s + B d e P G i M V / + N N / + N X d i D g p M 0 m c y 8 l 8 4 b P x Z c G 9 f 9 d g p r 6 x u b W 8 X t 0 s 7 u 3 v 5 B + f C o p a N E M W y y S E S q 7 V O N g k t s G m 4 E t m O F N P Q F P v r j 2 8 x / n K D S P J I P Z h p j L 6 R D y Q

F 3 /
o 2 T N g t t P T B w O O d e 7 p k T J J w p 7 T j f 1 s b m 1 v b O b m W v u n 9 w e H R s n 9 S 6 K k 4 l h Q 6 N e S z 7 A V H A m Y C O Z p p D P 5 F A o o B D L 5 j e F X 5 v B l K x W D z q e Q J + R M a C h Y w S b a S h X f M 0 4 y P I v I j o S R B m s z w f 2 n W n 4 S y A 1 4 l b 7 q l 9 d J 5 6 r u X d c b D 4 1 a s 1 H U U Y Y z O I d L 8 O A G m n A P L W g D g w S e 4 R X e n N R 5 c d 6 d j 6 W 1 5 B Q 7 p / A H z u c P J z S R u g = = < / l a t e x i t > X P < l a t e x i t s h a 1 _ b a s e 6 4 = " r s 7 R f o 2 B n B 3 1 l T d o q O D A p e q U Q 9

y 6 5 w 3 7
o t G 8 a 9 Z b z b K O K j p G J + g M 2 e g S t d A t a q M O o u g R P a N X 9 G Y 8 G S / G u / E x H 6 0 Y 5 c 4 h + g P j 8 w c P y 5 W V < / l a t e x i t > XD < l a t e x i t s h a 1 _ b a s e 6 4 = " c z z 7 n k i y j c r B w 4 a t o X n T 7 S 6 d / / l a t e x i t > Y P < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 r J Y W + i 5 w F m x s l P D/ W Q 0 R b I p G B A = " > A A A B 8 3 i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s x I U Z c F X b i s Y B / S G U o m z b S h S W Z I M k I Z + h t u X C j i1 p 9 x 5 9 + Y a W e h r Q c C h 3 P u 5 Z 6 c M O F M G 9 f 9 d k p r 6 x u b W + X t y s 7 u 3 v 5 B 9 f C o o + N U E d o m M Y 9 V L 8 S a c i Z p 2 z D D a S 9 R F I u Q 0 2 4 4 u c n 9 7 h N V m s X y w U w T G g g 8 k i x i B B s r + b 7 A Z h x G 2 e N s c D u o 1 t y 6 O w d a J V 5 B a l C g N a h + + c O Y r K M I J n M I 5 e H A F d b i F B j S B g Y J n e I U 3 x z g v z r v z s R g t O P n O M f y B 8 / k D z N a Q + A = = < / l a t e x i t > Y < l a t e x i t s h a 1 _ b a s e 6 4 = " A 5 P j c 2 u v e v s H K P W w D 7 5 k c 4 x Z 0 q 1 p k I w 5 0 9 e N L 2 z h n n e a N 4 1 6 + 1 m G U c V H I F j c A p M c A H a 4 B Z 0 Q B d g 8 A i e w S t 4 0 5 6 0 F + 1 D + 5 y 1 V r R y 5 h D 8 k / b 9 A 7 k t p X g = < / l a t e x i t > [X P XD ] = X < l a t e x i t s h a 1 _ b a s e 6 4 = " p 4 2 7 z B g I O 8 w p U D C 4 / 1
7 q l 9 d J 5 6 r u X d c b D 4 1 a s 1 H U U Y Y z O I d L 8 O A G m n A P L W g D g w S e 4 R X e n N R 5 c d 6 d j 6 W 1 5 B Q 7 p / A H z u c P J z S R u g = = < / l a t e x i t > X P < l a t e x i t s h a 1 _ b a s e 6 4 = " r s 7 R f o 2 B n B 3 1 l T d o q O D A p e q U Q 9 1 b 3 6 h u 1 r a 2 d 3 b 3 z P 2 D r o p T S V m H x i K W f Z 8 o J n j E O s B B s H 4 i G Q l 9 w X r + 5 L r w e w 9 M K h 5 H 9 z B N m B u S U c Q D T g l o y T O P H O B i y D I n J D D 2 g 6 y f 5 9 6 N Z 9 a t h j U DX i Z 2 S e q o R N s z v 5 x h T N O Q R U A F U W p g W w m 4 G Z H A q W B 5 z U k V S w i d k B E b a B q R k C k 3 m 6 X P 8 a l W h j i I p X 4 R 4 J n 6 e y M j o V L T 0 N e T R U i 1 6 B X i f 9 4 g h e D K z X i U p M A i O j 8 U p A J D j I s q 8 J B L R k F M N S F U c p 0 V 0 z G R h I I u r K Z L s B e / v Ey 6 5 w 3 7 o t G 8 a 9 Z b z b K O K j p G J + g M 2 e g S t d A t a q M O o u g R P a N X 9 G Y 8 G S / G u / E x H 6 0 Y 5 c 4 h + g P j 8 w c P y 5 W V < / l a t e x i t > XD < l a t e x i t s h a 1 _ b a s e 6 4 = " K 4 + n V U C 3 e Z q K l v 2 C N B w k N N U H F + o = " > A A A B 8 X i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s x I U Z c F N y 4 r 2 I e 2 Q 8 m k m T Y 0 k x m S O 0 I Z + h d u X C j i 1 r 9 x 5 9 + Y t r P Q 1 g O B w z n 3 k n r K M I J n M I 5 e H A F d b i F B j S B g Y J n e I U 3 x z g v z r v z s R g t O P n O M f y B 8 / k D z N a Q + A = = < / l a t e x i t > Y < l a t e x i t s h a 1 _ b a s e 6 4 = " r s 7 R f o 2 B n B 3 1 l T d o q O D A p e q U Q 9
j k T w u e n 6 H / S O L b d U 7 t 4 X S x U z u d x Z G E P 9 u E Q X D i D C l x C F e p A I I R 7 e I Q n S 1 o P 1 r P 1 M m v N W P O Z X f g G 6 / U D 9 / 2 Q h g = = < / l a t e x i t > x n,N T < l a t e x i t s h a 1 _ b a s e 6 4 = " V t 7 o h Y h 3 s N H S p A S 3 X f d t L w A G e F s = " > A A A B 7 3 i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g I e F Og l p Y B G w s I 5 g P S I 6 w t 9 l L l u z t n r t 7 Y j j y J 2 w s F L H 1 7 9 j 5 b 9 x L r t D E B w O P 9 2 a Y r b / h z r 9 x + l h o 6 4 E L h 3 P u 5 d 5 7 o o x R p R 3 n 2 y o t L a + s r p X X K x u b W 9 s 7 9 u 5 e W 4 l c Y u J h w Y S 8 j 5 r b l 1 d w 6 0 S r y C 1 K B A a 1 j 9 G o w i k g g q D e F Y 6 7 7 n x s b P s D K M c D q r D B J N Y 0 y m e E z 7 l k o s q P a z+ b k z d G a V E Q o j Z U s a N F d / T 2 R Y a J 2 K w H Y K b C Z 6 2 c v F / 7 x + Y s I b P 2 M y T g y V Z L E o T D g y E c p / R y O m K D E 8 t Q Q T x e y t i E y w w s T Y h C o 2 B G / 5 5 V X S u a x 7 V / X G Q 6 P W v C 3 i K M M J n M I 5e H A N T b i H F r S B w B S e 4 R X e n N h 5 c d 6 d j 0 V r y S l m j u E P n M 8 f C a m P X w = = < / l a t e x i t > y m,1 . . . . . .< l a t e x i t s h a 1 _ b a s e 6 4 = " t W W o c n A l P F 8 o T h z W 8 + U J d R N y J n 0 = " > A A A B 8 H i c d V D L S s N A F L 3 x W e u r 6 t L N Y B F c S E i 0 1 a 6 k 4 M a V V O h L 2 h A m 0 0 k 7 d C Y J M x O h h H 6 F G x e K u P V z 3 P k 3 T h + C z w M X D u f c y 7 3 3 B A l n S j v O u 7 W w u L S 8 s p p b y 6 9 v b G 5 t F 3 Z 2 m y p O k o b 0 L 4 / B T 9 T 5 o n t n t m l 2 5 K x e r F P I 4 c 7 M M B H I E L 5 1 C F K 6 h B A w g I u I d H e L K k 9 W A 9 W y + z 1 g V r P r M H 3 2 C 9 f g D 4 m p C I < / l a t e x i t > y m,N T < l a t e x i t s h a 1 _ b a s e 6 4 = " B Y / 8 u t d P x P y N 4 Y J Q O 7 E b S l N B w r P 6 e S H G o 1 D A M T G e I d V / N e p n 4 n 9 d K d P e s n T I R J 5 o K M l n U T T j U E c z C g R 0 m K d F 8 a A g m k p l b I e l j i Y k 2 E d o m B G / 2 5 X n S O C 5 5 J 6 X y d b l

Figure 12
Figure12depicts the convergence behavior of the proposed ISAC algorithms in terms of the MSE of the voxel coefficient soft-replicas vk:m,t and the SER of the data symbol soft-replicas xn,t:m,t , for three selected combinations of their respective damping parameters η x ∈ [0, 1] and η v ∈ [0, 1], following the damped update rule given by[67] x(τ) n,t:m,t

Fig. 12 :
Fig. 12: Convergence behavior of the proposed AL-ISAC and Bi-ISAC algorithms with varying damping factors η x and η v with ρ = 0.1, E v = 1.5%, SNR = 15dB, as a function of MP iterations in a system with N U = 4, N A N R = 12, N V = 512, and N T = 100. ȳx Until termination criteria is satisfied * do 1: Compute the effective channel matrix G ≜ H + Adiag(v)B; 2: Initialize the soft-replica at all variable nodes as xn,t:m,t = E xn,t [x n,t ]; ▷ ∀m, n, t 3: Initialize the MSE at all variable nodes as ψ x n,t:m,t via (21); ▷ ∀m, n, t 11: Project xn,t to the symbol constellation X ; ▷ ∀n, t 12: Output projected xn,t as final hard estimate;