Permutation-Driven AFDM DFE: From Sparse Detection Ordering to Kendall-Tau Structural Diversity

1. Introduction

Affine frequency division multiplexing (AFDM) provides an effective chirp-based multicarrier solution for next-generation wireless communications resilient to doubly dispersive fading [1]. By matching chirp sweep rates to the delay-Doppler profile of the channel through the discrete affine Fourier transform (DAFT), AFDM maps dense time-varying multipath into a quasi-sparse representation, thereby mitigating the inter-carrier interference (ICI) that destroys the subcarrier orthogonality of conventional orthogonal frequency division multiplexing (OFDM) [2]. This resilience to doubly selective dispersion is critical for high-speed railway and vehicular links, unmanned aerial vehicle (UAV)/air-to-ground and non-terrestrial connections, and millimeter-wave/sub-terahertz (sub-THz) integrated sensing and communications (ISAC) deployments [3,4,5]. Within these scenarios, recent studies position AFDM as a compact alternative to orthogonal time frequency space (OTFS) [6,7]. Furthermore, AFDM achieves full diversity gain with fast Fourier transform (FFT)-compatible complexity [8], resists hardware impairments in multiple-input multiple-output (MIMO) systems [9], and supports higher spectral efficiency through generalized spatial modulation (GSM) [10,11] or relaxed orthogonality constraints [12]. These transmission benefits position AFDM as a promising candidate waveform for future networks, but they also make reliable low-complexity detection over doubly selective channels a central receiver-design challenge.

Although recent channel-estimation advances [13] can provide accurate channel state information (CSI), reliable high-spectral-efficiency AFDM detection remains challenging because rich multipath and fractional delay and Doppler shifts induce severe inter-symbol interference despite DAFT-domain sparsity. Linear equalizers, such as zero-forcing (ZF) and linear minimum mean-square error (LMMSE) [14,15], therefore often remain sub-optimal, trading noise enhancement against residual interference under deep fading. Iterative detectors have therefore been investigated, but with distinct limitations. Iterative parallel interference cancellation (I-PIC) schemes [16] feed back block-wise tentative decisions to iteratively subtract interference, but suffer severe block-wise error propagation across iterations, where localized decision errors contaminate the entire block in subsequent feedback steps. Probabilistic message-passing paradigms offer statistical error control, yet low-complexity approximations such as the message passing algorithm (MPA) [17] lack robustness under high-order constellations and require stringent damping control. Conversely, exact matrix-solving counterparts such as orthogonal approximate message passing (OAMP) [18] achieve robust stability but incur high latency from numerous internal iterations, precluding low-latency deployment.

To circumvent the computational burden of iterative schemes, conventional DFE exploits structured channel interference for non-iterative reception [19,20], and the DAFT-domain symbol separation of AFDM naturally aligns with its successive interference cancellation (SIC) mechanism. However, DFE faces two practical bottlenecks. First, under doubly selective channels and high-order modulations, causal cancellation propagates early decision errors to subsequent stages. Second, single-branch DFE exhibits uneven symbol-wise reliability, weakening its synergy with channel coding schemes such as low-density parity-check (LDPC) codes [21]. Bidirectional DFE (BDFE) using time reversal [22,23,24] and bidirectional arbitrated decision-feedback equalization (BAD) [25] introduce branch diversity to mitigate error propagation, but they do not structurally decouple errors or optimize trajectories.

To overcome these bottlenecks, we introduce a permutation-driven DFE (P-DFE) methodology. Traditional single-carrier DFE (SC-DFE) [20] is restricted by chronological causality. Block-based waveforms such as AFDM and OTFS [26] remove this temporal constraint, yet existing block DFE and multi-branch BDFE receivers [22,24,27] underuse this freedom by confining detection to predefined natural or time-reversed sequences. Inspired by MIMO detection [28,29], we formulate sequence permutation as an explicit optimization variable. The resulting framework optimizes single-branch reliability and supports multi-branch structural diversity.

Permutation thus serves a dual role. First, it shapes interference by reshaping the signal-to-interference-plus-noise ratio (SINR) profile to suppress initial errors. Second, it uses the resulting single-branch reliability to generate parallel structural diversity. Using the Kendall–Tau (KT) distance [30] to quantify structural dissimilarity, we formulate a criterion for diverse yet quality-preserving decision trajectories. This reduces shared error-propagation paths, improves demodulated soft information, and integrates with modern channel coding.

The proposed framework differs from prior art in its use of detection order. Linear ZF/LMMSE equalizers lack decision feedback and are limited by residual interference or noise enhancement; conventional DFE, BDFE, and BAD-type receivers exploit feedback but remain tied to natural, reversed, or predefined trajectories; and iterative MPA/OAMP-type receivers improve reliability through repeated inference or matrix solving, with higher latency and cost. In contrast, P-DFE optimizes the detection permutation to shape single-branch reliability, while the proposed permutation-driven quad-branch DFE (P-QDFE) uses structurally diverse permutations to suppress correlated error propagation without iterative detection. The main contributions are summarized as follows.

1)

We establish the P-DFE framework by formulating an exact mapping between the physical symbol domain and the sequential decision domain. Using this mapping, we reconstruct the interference model and derive recursions for the feedforward disturbance and propagated feedback error. This decomposition identifies the feedforward disturbance variance as the primary bottleneck, motivating permutation-based optimization.

2)

Formulating permutation as a mechanism to shape the feedforward disturbance variance profile, we apply majorization theory under a genie-aided assumption to establish Sorted QR Decomposition (SQRD) [29] and Column-Pivoted QR (CPQR) [31] as statistical performance references for the proposed P-DFE. Guided by these references, we develop the Energy Graph (EG) ordering algorithm. By exploiting the sparse AFDM interference topology, EG approximates the SQRD performance while bypassing the complexity of repeated matrix factorization.

3)

We extend permutation to parallel branch generation and propose permutation-driven multi-branch DFE (P-MBDFE). P-QDFE uses a quality-preserving shift/reversal family and a conservative KT harmonic-mean criterion. Starting from the physical symbol-wise error correlation, matched-ancestor and block-average relaxations yield the KT-based correlation model, calibrated by a linear scaling factor, for soft-output fusion; with LDPC coding, P-QDFE matches the evaluated iterated OAMP receiver at the reported pre-decoding bit error rate (BER) thresholds, while support-truncated incomplete Cholesky factorization (ICHOL) gives LMMSE-level critical-path complexity.

The paper is organized as follows. Section II introduces the AFDM transceiver and structured channel model. Section III develops P-DFE via permutation-based interference reconstruction, SINR, and error decomposition. Section IV establishes the genie-aided references and EG ordering algorithm. Section V introduces the non-iterative parallel P-MBDFE scheme, focusing on quality-preserving branch generation and KT optimization. Section VI presents simulation results, and Section VII concludes the paper.

Notations: Scalars, vectors, and matrices are denoted by lowercase, bold lowercase, and bold uppercase letters, respectively. For matrix entries, subscripts denote row and column. The operators ${(·)}^T$, ${(·)}^H$, and ${(·)}^{*}$ denote transpose, Hermitian transpose, and complex conjugate. ${\mathbf{I}}_N$ is the $N\times N$ identity matrix, and $\mathbf{F}$ is the normalized discrete Fourier transform (DFT) matrix. $\mathrm{diag}(·)$ forms a diagonal matrix. $\mathcal{CN}(\mu ,{\sigma }^2)$ represents a complex Gaussian distribution with mean $\mu$ and variance ${\sigma }^2$. We use $Q(·)$ for the Gaussian Q-function, while $\mathcal{Q}(·):\mathbb{C}\to \mathcal{A}$ denotes hard slicing to the nearest constellation point. Here ${[·]}_N$ denotes modulo-N reduction, $\lceil ·\rceil$ the ceiling operation, $u_n^{\downarrow }$ the n-th largest component of a vector $\mathbf{u}$, and $\mathbb{I}\{·\}$ the binary indicator function. $\mathcal{A}$ denotes the signal constellation.

2. AFDM System and Channel Model

AFDM employs the DAFT to map data symbols onto orthogonal chirps, transforming a dense doubly selective channel into a quasi-sparse effective channel. For zero-mean equiprobable symbols $\mathbf{x}\in {\mathcal{A}}^N$ with $\mathbb{E}\left[\mathbf{x}{\mathbf{x}}^H\right]=E_s{\mathbf{I}}_N$, the inverse DAFT gives $\mathbf{s}={\mathbf{A}}^H\mathbf{x}$. The forward DAFT matrix is $\mathbf{A}={\mathrm{\Lambda }}_{c_2}\mathbf{F}{\mathrm{\Lambda }}_{c_1}$, where ${\left[\mathbf{F}\right]}_{n,m}=N^{-1/2}{e}^{-j2\pi nm/N}$ and ${\mathrm{\Lambda }}_c=\mathrm{diag}{\left(e^{-j2\pi c{n}^2}\right)}_{n=0}^{N-1}$. Both chirp matrices are unitary for real c; $c_1$ controls delay-Doppler separation and $c_2$ the secondary chirp phase. Frequently used quantities are summarized in Table 1; local indices and auxiliary variables are defined at their first occurrence.

A chirp-periodic prefix (CPP) of length $L_{\mathrm{cp}}\ge \lceil {\mathcal{l}}_{\max }\rceil$ enables a circular block model, with ${\mathcal{l}}_{\max }$ the maximum path delay in samples. To accommodate both integer- and fractional-delay paths, define the circulant delay operator

$${\mathrm{\Pi }}_{\mathcal{l}}={\mathbf{F}}^H\mathrm{diag}{\left(e^{-j2\pi k\mathcal{l}/N}\right)}_{k=0}^{N-1}\mathbf{F},\mathcal{l}\in \mathbb{R},$$

(1)

which reduces to the usual cyclic-shift power ${\mathrm{\Pi }}^{\mathcal{l}}$, with $\mathrm{\Pi }\triangleq {\mathrm{\Pi }}_1$, when ℓ is integer. The signal propagates through a linear time-varying (LTV) channel modeled as

$${\mathbf{H}}_{\mathrm{tv}}={\displaystyle \sum _{i=1}^L}{h}_i{\mathrm{\Delta }}_{{\nu }_i}{\mathrm{\Pi }}_{{\mathcal{l}}_i},$$

(2)

where L is the path count, and $h_i$, ${\mathcal{l}}_i$, and ${\nu }_i$ are the complex gain, normalized delay in sampling intervals, and normalized Doppler shift of the i-th path, respectively. Here, ${\mathcal{l}}_i$ and ${\nu }_i$ may both be fractional, and ${\mathrm{\Delta }}_{\nu }=\mathrm{diag}{\left(e^{j2\pi \nu n/N}\right)}_{n=0}^{N-1}$. After CPP removal, $\mathbf{r}={\mathbf{H}}_{\mathrm{tv}}\mathbf{s}+\mathbf{w}$, with $\mathbf{w}\sim \mathcal{CN}(\mathbf{0},{\sigma }^2{\mathbf{I}}_N)$ independent of $\mathbf{x}$. Unless otherwise specified, expectations below are conditioned on the realized ${\mathbf{H}}_{\mathrm{eff}}$. Applying the forward DAFT yields

$$\mathbf{y}=\mathbf{A}\mathbf{r}={\mathbf{H}}_{\mathrm{eff}}\mathbf{x}+\tilde{\mathbf{w}},$$

(3)

where $\tilde{\mathbf{w}}=\mathbf{A}\mathbf{w}$ remains white Gaussian noise and ${\mathbf{H}}_{\mathrm{eff}}\triangleq \mathbf{A}{\mathbf{H}}_{\mathrm{tv}}{\mathbf{A}}^H\in {\mathbb{C}}^{N\times N}$ is the effective DAFT-domain channel. For the integer-delay specialization, substituting the DAFT definition and (2) into (3) gives the compact $(n,m)$-th entry

$${\left[{\mathbf{H}}_{\mathrm{eff}}\right]}_{n,m}={\displaystyle \sum _{i=1}^L}{\tilde{h}}_i(n,m){{\rm Y}}_N\left(q_i(n,m)\right),$$

(4)

where $n,m\in \{0,\dots ,N-1\}$, $q_i(n,m)\triangleq {\nu }_i-{[n-m]}_N-{\xi }_i$, ${\xi }_i\triangleq 2c_{1}N{\mathcal{l}}_i$, and

$$\begin{array}{cc}\hfill {\tilde{h}}_i(n,m)\triangleq h_i\exp \left\{j2\pi \right[& c_2(m^2-n^2)-{\displaystyle \frac{m{\mathcal{l}}_i}{N}}+c_1{\mathcal{l}}_i^2\hfill \\ \hfill & +q_i(n,m)\left(1-{\displaystyle \frac{1}{N}}\right)\left]\right\}.\hfill \end{array}$$

(5)

The Dirichlet kernel ${{\rm Y}}_N\left(x\right)={\displaystyle \frac{\sin \left(\pi x\right)}{N\sin (\pi x/N)}}{e}^{-j\pi x(1-1/N)}$ is interpreted by continuous extension at removable singularities. For fractional ${\mathcal{l}}_i$, (3) remains exact through (1), whereas (4)–(5) is the integer-delay specialization and fractional paths additionally exhibit delay-domain leakage. To separate the dominant integer-grid bins, $c_1$ is tuned as [1]

$$c_1={\displaystyle \frac{2{\alpha }_{\max }+1}{2N}},\mathrm{s}.\mathrm{t}.(2{\alpha }_{\max }+1)\lceil {\mathcal{l}}_{\max }\rceil +2{\alpha }_{\max }<N,$$

(6)

where ${\nu }_{\max }\triangleq {\max }_i\left|{\nu }_i\right|$ and ${\alpha }_{\max }=\lceil {\nu }_{\max }\rceil$. Thus, ${\mathbf{H}}_{\mathrm{eff}}$ is structurally sparse on the integer grid. When both delay and Doppler lie on that grid, each path occupies one dominant diagonal. Fractional Doppler and delay cause leakage; in that case (6) separates the dominant bins rather than eliminating leakage, whose energy nevertheless remains concentrated around those diagonals. This quasi-sparsity underpins the proposed equalizer.

3. Permutation-Driven Interference Reconstruction and P-DFE Scheme

Exploiting the block-wise reception of AFDM, P-DFE permutes the symbols within each received block before successive detection, thereby replacing the fixed symbol-detection order of conventional DFE with a configurable one. This ordering affects both its feedforward disturbance and its exposure to propagated decision errors. This section establishes the physical-to-stage mapping and derives the order-dependent reliability model used for permutation design. Figure 1 organizes the receiver into the four processing steps followed below.

3.1. P-DFE Model and Interference Reconstruction

Once permutation is introduced, the physical index of a symbol no longer coincides with its detection stage. We therefore distinguish the physical symbol domain from the detection-stage domain as

$$\mathcal{S}\triangleq \{1,\dots ,N\},\mathcal{K}\triangleq \{1,\dots ,N\},$$

(7)

where $\mathcal{S}$ indexes the transmitted information symbols and $\mathcal{K}$ indexes their processing stages after permutation. The designed permutation matrix then defines a complete mapping between these two domains.

Let $\pi ={[{\pi }_1,\dots ,{\pi }_N]}^T$ denote the stage-order vector, where ${\pi }_k$ is the physical symbol assigned to triangular stage k. Since backward substitution proceeds from N to 1, the chronological decision trajectory is $[{\pi }_N,\dots ,{\pi }_1]$. The inverse rank map and associated permutation matrix are defined by

$${\varphi }_{\pi }\left({\pi }_k\right)=k,{\left[{\mathbf{P}}_{\pi }\right]}_{k,m}=\mathbb{I}\{m={\pi }_k\},$$

(8)

where ${\varphi }_{\pi }:\mathcal{S}\to \mathcal{K}$. The receiver flow is then expressed in the four steps of Figure 1.

1) Permutation ($\mathbf{y}\to {\mathbf{y}}_{\pi }$, $\mathbf{G}\to {\mathbf{G}}_{\pi }$): The permutation block maps the received model to the stage domain as

$${\mathbf{y}}_{\pi }={\mathbf{H}}_{\pi }{\mathbf{x}}_{\pi }+{\tilde{\mathbf{w}}}_{\pi },$$

(9)

where ${\mathbf{x}}_{\pi }={\mathbf{P}}_{\pi }\mathbf{x}$, ${\mathbf{y}}_{\pi }={\mathbf{P}}_{\pi }\mathbf{y}$, ${\mathbf{H}}_{\pi }={\mathbf{P}}_{\pi }{\mathbf{H}}_{\mathrm{eff}}{\mathbf{P}}_{\pi }^T$, and ${\tilde{\mathbf{w}}}_{\pi }={\mathbf{P}}_{\pi }\tilde{\mathbf{w}}\sim \mathcal{CN}(\mathbf{0},{\sigma }^2{\mathbf{I}}_N)$. Let $\lambda \triangleq {\sigma }^2/E_s$ denote the noise-to-symbol-energy ratio. The corresponding noise-regularized Gram matrices are

$$\mathbf{G}\triangleq {\mathbf{H}}_{\mathrm{eff}}^H{\mathbf{H}}_{\mathrm{eff}}+\lambda {\mathbf{I}}_N,{\mathbf{G}}_{\pi }={\mathbf{P}}_{\pi }\mathbf{G}{\mathbf{P}}_{\pi }^T={\mathbf{H}}_{\pi }^H{\mathbf{H}}_{\pi }+\lambda {\mathbf{I}}_N.$$

(10)

This operation does not change the physical channel; it changes the triangular couplings and pivot sequence seen by the detector.

2) Forward Filtering (${\mathbf{y}}_{\pi }\to {\mathbf{z}}_{\pi }$): The feedforward (FF) equalizer in Figure 1 computes the Cholesky factorization of the Hermitian positive-definite matrix ${\mathbf{G}}_{\pi }$ with its diagonal scaling separated:

$${\mathbf{G}}_{\pi }={\mathbf{L}}_{\pi }{\mathbf{D}}_{\pi }{\mathbf{L}}_{\pi }^H,$$

(11)

where ${\mathbf{L}}_{\pi }$ is unit lower triangular and ${\mathbf{D}}_{\pi }$ is positive diagonal; equivalently, the conventional Cholesky factor is ${\mathbf{L}}_{\pi }{\mathbf{D}}_{\pi }^{1/2}$. These factors are computed recursively as

$$\begin{array}{cc}\hfill d_k^{\left(\pi \right)}& =g_{k,k}^{\left(\pi \right)}-{\displaystyle \sum _{j=1}^{k-1}}{\left|l_{k,j}^{\left(\pi \right)}\right|}^2{d}_j^{\left(\pi \right)},\hfill \end{array}$$

(12a)

$$\begin{array}{cc}\hfill l_{i,k}^{\left(\pi \right)}& ={\displaystyle \frac{g_{i,k}^{\left(\pi \right)}-{\sum }_{j=1}^{k-1}{l}_{i,j}^{\left(\pi \right)}{d}_j^{\left(\pi \right)}{l}_{k,j}^{\left(\pi \right)*}}{d_k^{\left(\pi \right)}}},i>k,\hfill \end{array}$$

(12b)

where $g_{i,k}^{\left(\pi \right)}={\left[{\mathbf{G}}_{\pi }\right]}_{i,k}$. The FF equalizer output is

$$\begin{array}{cc}\hfill {\mathbf{z}}_{\pi }& ={\mathbf{L}}_{\pi }^{-1}{\mathbf{H}}_{\pi }^H{\mathbf{y}}_{\pi }={\mathbf{D}}_{\pi }{\mathbf{L}}_{\pi }^H{\mathbf{x}}_{\pi }+{\tilde{\mathbf{n}}}_{\pi },\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\tilde{\mathbf{n}}}_{\pi }& ={\mathbf{L}}_{\pi }^{-1}\left({\mathbf{H}}_{\pi }^H{\tilde{\mathbf{w}}}_{\pi }-\lambda {\mathbf{x}}_{\pi }\right).\hfill \end{array}$$

(14)

3) Backward Substitution (${\mathbf{z}}_{\pi }\to {\tilde{\mathbf{x}}}_{\pi }$): The feedback (FB) equalizer in Figure 1 performs backward substitution for $k=N,N-1,\dots ,1$. Equations (11) and (14) give $\mathbb{E}\left[\right|{\tilde{n}}_{\pi ,k}{|}^2]={\sigma }^2{d}_k^{\left(\pi \right)}$ and $\mathbb{E}\left[{\tilde{n}}_{\pi ,k}{x}_{\pi ,k}^{*}\right]=-{\sigma }^2$. Its scalar input from the FF equalizer is

$$z_{\pi ,k}=d_k^{\left(\pi \right)}{x}_{\pi ,k}+{\displaystyle \sum _{i=k+1}^N}{d}_k^{\left(\pi \right)}{l}_{i,k}^{\left(\pi \right)*}{x}_{\pi ,i}+{\tilde{n}}_{\pi ,k}.$$

(15)

For previously detected stages $i>k$, let ${\widehat{x}}_{\pi ,i}=\mathcal{Q}\left({\tilde{x}}_{\pi ,i}\right)$ and $\mathrm{\Delta }{x}_{\pi ,i}\triangleq x_{\pi ,i}-{\widehat{x}}_{\pi ,i}$. Substitution into (15) gives the pre-normalized soft output

$$\begin{array}{cc}\hfill {\breve{x}}_{\pi ,k}& \triangleq {\displaystyle \frac{z_{\pi ,k}-{\sum }_{i=k+1}^N{d}_k^{\left(\pi \right)}{l}_{i,k}^{\left(\pi \right)*}{\widehat{x}}_{\pi ,i}}{d_k^{\left(\pi \right)}}}\hfill \\ \hfill & =x_{\pi ,k}+{\displaystyle \frac{{\tilde{n}}_{\pi ,k}}{d_k^{\left(\pi \right)}}}+{\displaystyle \sum _{i=k+1}^N}{l}_{i,k}^{\left(\pi \right)*}\mathrm{\Delta }{x}_{\pi ,i}.\hfill \end{array}$$

(16)

Under correct prior decisions, these moments give $\mathbb{E}\left[{\breve{x}}_{\pi ,k}{x}_{\pi ,k}^{*}\right]/E_s={\mu }_{\pi ,k}$, where ${\mu }_{\pi ,k}\triangleq 1-\lambda /d_k^{\left(\pi \right)}=(d_k^{\left(\pi \right)}-\lambda )/d_k^{\left(\pi \right)}>0$. The unbiased soft output is therefore

$$\begin{array}{cc}\hfill {\tilde{x}}_{\pi ,k}\triangleq {\displaystyle \frac{{\breve{x}}_{\pi ,k}}{{\mu }_{\pi ,k}}}& =x_{\pi ,k}+\underset{\triangleq n_{\mathrm{f},\pi ,k}}{\underbrace{{\displaystyle \frac{{\tilde{n}}_{\pi ,k}+\lambda x_{\pi ,k}}{d_k^{\left(\pi \right)}-\lambda }}}}\hfill \\ \hfill & +\underset{\triangleq n_{\mathrm{b},\pi ,k}}{\underbrace{{\displaystyle \frac{d_k^{\left(\pi \right)}}{d_k^{\left(\pi \right)}-\lambda }}{\sum }_{i=k+1}^N{l}_{i,k}^{\left(\pi \right)*}\mathrm{\Delta }{x}_{\pi ,i}}}.\hfill \end{array}$$

(17)

Thus the slicer operates on an unbiased output whose error is decomposed into an orthogonal FF disturbance and a propagated FB error.

4) Remapping (${\tilde{\mathbf{x}}}_{\pi }\to \tilde{\mathbf{x}}$): Collecting the scalar outputs gives ${\tilde{\mathbf{x}}}_{\pi }$, and the inverse permutation $\tilde{\mathbf{x}}={\mathbf{P}}_{\pi }^H{\tilde{\mathbf{x}}}_{\pi }={\mathbf{P}}_{\pi }^T{\tilde{\mathbf{x}}}_{\pi }$ returns them to the physical symbol domain. Hence $e_{\pi ,k}\triangleq {\tilde{x}}_{\pi ,k}-x_{\pi ,k}=n_{\mathrm{f},\pi ,k}+n_{\mathrm{b},\pi ,k}$, and performance is governed by the stage-wise SINR profile.

3.2. Order-Dependent Reliability Model

3.2.1. FF Disturbance Variance

The FF disturbance $n_{\mathrm{f},\pi ,k}$ is the stochastic trigger for decision errors. Using the two moments above in (17), its variance is

$${\sigma }_{\mathrm{f},\pi ,k}^2\triangleq \mathbb{E}\left[\right|n_{\mathrm{f},\pi ,k}{|}^2]={\displaystyle \frac{{\sigma }^2}{d_k^{\left(\pi \right)}-\lambda }}.$$

(18)

Substituting (12a) into (18) gives

$${\sigma }_{\mathrm{f},\pi ,k}^2={\sigma }^2{\left(g_{k,k}^{\left(\pi \right)}-{\displaystyle \sum _{j=1}^{k-1}}{\left|l_{k,j}^{\left(\pi \right)}\right|}^2{d}_j^{\left(\pi \right)}-\lambda \right)}^{-1}.$$

(19)

Equation (19) shows that this variance is conditioned on the preceding factorization stages rather than being an intrinsic property of physical symbol ${\pi }_k$. Different permutations therefore produce different residual projections. Since $\lambda$ is permutation independent, maximizing $d_k^{\left(\pi \right)}$ remains equivalent to minimizing ${\sigma }_{\mathrm{f},\pi ,k}^2$, making the diagonal pivot a direct surrogate for FF disturbance ordering.

3.2.2. FB Error Propagation

The FB error variance is governed jointly by the triangular coupling and the squared hard-decision error. To distinguish the latter from symbol error probability, define

$${\mathcal{M}}_{\mathrm{\Delta }}(\mathcal{A},v)\triangleq {\displaystyle \frac{1}{\left|\mathcal{A}\right|}}\sum _{a\in \mathcal{A}}{\mathbb{E}}_{u\sim \mathcal{CN}(0,v)}\left[{|a-\mathcal{Q}(a+u)|}^2\right].$$

(20)

In the high-reliability nearest-neighbor regime, ${\mathcal{M}}_{\mathrm{\Delta }}(\mathcal{A},v)\approx P_s(\mathcal{A},v)d_{\min }^2$. Under the scalar equivalent-Gaussian propagation model and neglecting cross-correlations among hard-decision errors, (17) gives

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{b},\pi ,k}^2& \triangleq \mathbb{E}\left[\right|n_{\mathrm{b},\pi ,k}{|}^2]\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill & \approx {\displaystyle \frac{1}{{\mu }_{\pi ,k}^2}}{\displaystyle \sum _{i=k+1}^N}{\left|l_{i,k}^{\left(\pi \right)}\right|}^2{\mathcal{M}}_{\mathrm{\Delta }}\left(\mathcal{A},{\sigma }_{\mathrm{f},\pi ,i}^2+{\sigma }_{\mathrm{b},\pi ,i}^2\right).\hfill \end{array}$$

(22)

Substituting (12) and ${\mu }_{\pi ,k}=(d_k^{\left(\pi \right)}-\lambda )/d_k^{\left(\pi \right)}$ into (22) yields

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{b},\pi ,k}^2\approx {\displaystyle \sum _{i=k+1}^N}& {\left|{\displaystyle \frac{g_{i,k}^{\left(\pi \right)}-{\sum }_{j=1}^{k-1}{l}_{i,j}^{\left(\pi \right)}{d}_j^{\left(\pi \right)}{l}_{k,j}^{\left(\pi \right)*}}{d_k^{\left(\pi \right)}-\lambda }}\right|}^2\hfill \\ \hfill & \times {\mathcal{M}}_{\mathrm{\Delta }}\left(\mathcal{A},{\sigma }_{\mathrm{f},\pi ,i}^2+{\sigma }_{\mathrm{b},\pi ,i}^2\right).\hfill \end{array}$$

(23)

Equation (23) shows that permutation jointly reshapes the residual triangular coupling, the effective denominator $d_k^{\left(\pi \right)}-\lambda$, and the reliability of later decisions. The FB recursion starts from ${\sigma }_{\mathrm{b},\pi ,N}^2=0$, and ${\mathcal{M}}_{\mathrm{\Delta }}$ may be stored in a look-up table.

3.2.3. Stage-Wise SINR

Under the same scalar auxiliary model, $n_{\mathrm{f},\pi ,k}$ and $n_{\mathrm{b},\pi ,k}$ in (17) are treated as zero mean, mutually uncorrelated, and uncorrelated with $x_{\pi ,k}$. Using the FF disturbance and FB error variances in (18) and (22), the effective stage-wise SINR is modeled as

$${\mathrm{SINR}}_{\pi ,k}={\displaystyle \frac{E_s}{{\sigma }_{\mathrm{f},\pi ,k}^2+{\sigma }_{\mathrm{b},\pi ,k}^2}}.$$

(24)

Direct optimization is difficult because ${\sigma }_{\mathrm{b},\pi ,k}^2$ recursively depends on later decisions. We therefore use the tractable FF disturbance profile as a surrogate: reducing its vulnerable stages lowers the corresponding decision-error moments, while the net FB effect remains ordering dependent through (23).

4. Shaping the FF Disturbance Profile via Permutation Design

Building on the reconstructed interference model, this section utilizes permutation to explicitly shape the FF disturbance variance profile. By analyzing the sequential Cholesky factors under a genie-aided assumption, we apply majorization theory to establish SQRD and CPQR as statistical performance references. Consequently, we exploit the sparse AFDM interference topology to develop the low-complexity EG ordering, which approximates the SQRD reference without repeated full factorization.

4.1. Genie-Aided Analysis: SQRD and CPQR Bounds

Under genie-aided (GA) FB decisions [32], ${\sigma }_{\mathrm{b},\pi ,k}^2=0$ and (24) becomes ${\mathrm{SINR}}_{\pi ,k}^{\mathrm{GA}}=d_k^{\left(\pi \right)}/\lambda -1$. The determinant invariant therefore applies to $1+{\mathrm{SINR}}_{\pi ,k}^{\mathrm{GA}}$ rather than to the physical SINR itself. Its geometric mean is

$$\overline{\gamma }\triangleq {\left({\displaystyle \prod _{k=1}^N}\left[1+{\mathrm{SINR}}_{\pi ,k}^{\mathrm{GA}}\right]\right)}^{1/N}={\displaystyle \frac{\det {\left(\mathbf{G}\right)}^{1/N}}{\lambda }}.$$

(25)

Hence $1+{\mathrm{SINR}}_{\pi ,k}^{\mathrm{GA}}=\overline{\gamma }{e}^{{\delta }_k}$, where the ordering-dependent fluctuation satisfies

$${\delta }_k\triangleq \ln \left({\displaystyle \frac{d_k^{\left(\pi \right)}}{\lambda \overline{\gamma }}}\right),{\displaystyle \sum _{k=1}^N}{\delta }_k=0.$$

(26)

where ${\delta }_{\pi }={[{\delta }_1,\dots ,{\delta }_N]}^T$ denotes the profile generated by permutation $\pi$.

Let $P_b\left(\gamma \right)$ denote the uncoded additive white Gaussian noise (AWGN) bit-error function of the adopted constellation and define $f\left({\delta }_k\right)\triangleq P_b(\overline{\gamma }{e}^{{\delta }_k}-1)$. A second-order Taylor expansion around ${\delta }_k=0$ motivates the proxy

$${\widetilde{\mathrm{BER}}}_{\pi }^{\mathrm{GA}}\triangleq f\left(0\right)+{\displaystyle \frac{f^{\prime \prime }\left(0\right)}{2}}{\displaystyle \frac{1}{N}}{\displaystyle \sum _{k=1}^N}{\delta }_k^2,$$

(27)

which approximates $N^{-1}{\sum }_k{P}_b\left({\mathrm{SINR}}_{\pi ,k}^{\mathrm{GA}}\right)$; the first-order term vanishes because ${\sum }_k{\delta }_k=0$. Let $C_1\triangleq f\left(0\right)$ and $C_2\triangleq f^{\prime \prime }\left(0\right)/2$. Because $\overline{\gamma }$ varies across channel realizations, the moments in (28) and (31) are conditioned on a fixed $\overline{\gamma }$, with this conditioning suppressed for notational consistency with the original formulation. The expectation becomes

$$\mathbb{E}\left[{\widetilde{\mathrm{BER}}}_{\pi }^{\mathrm{GA}}\right]=C_1+C_2\left(\underset{\triangleq \mathrm{\Psi }\left(\mathbb{E}\left[{\delta }_{\pi }\right]\right)}{\underbrace{{\displaystyle \frac{1}{N}}\sum _{k=1}^N{\left(\mathbb{E}\left[{\delta }_k\right]\right)}^2}}+\underset{\triangleq \mathrm{\Phi }\left({\delta }_{\pi }\right)}{\underbrace{{\displaystyle \frac{1}{N}}\sum _{k=1}^N\mathrm{Var}\left({\delta }_k\right)}}\right).$$

(28)

Here, $\mathrm{\Psi }$ captures the mean-profile distortion and $\mathrm{\Phi }$ captures its channel-to-channel fluctuation.

Definition 1

(Majorization and Schur-Convexity [33]). For $\mathbf{x},\mathbf{y}\in {\mathbb{R}}^N$ satisfying ${\sum }_{n=1}^N{x}_n={\sum }_{n=1}^N{y}_n$, $\mathbf{x}$ majorizes $\mathbf{y}$ (denoted by $\mathbf{x}\succ \mathbf{y}$) if $S_k\left(\mathbf{x}\right)\ge S_k\left(\mathbf{y}\right)$ for all $k=1,\dots ,N-1$, where $S_k\left(\mathbf{u}\right)\triangleq {\sum }_{n=1}^k{u}_n^{\downarrow }$. A function $\mathrm{\Psi }:{\mathbb{R}}^N\to \mathbb{R}$ is Schur-convex if $\mathbf{x}\succ \mathbf{y}$ implies $\mathrm{\Psi }\left(\mathbf{x}\right)\ge \mathrm{\Psi }\left(\mathbf{y}\right)$.

To derive the statistical performance references of P-DFE, SQRD [29] and CPQR [31] construct $\pi$ iteratively. At iteration k, let ${\mathcal{P}}_{k-1}\subset \mathcal{S}$ contain the physical symbols selected in the first $k-1$ factorization stages. In the notation of (12), the candidate pivot of a remaining symbol p is

$$d_{p\mid {\mathcal{P}}_{k-1}}=g_{p,p}-\sum _{q\in {\mathcal{P}}_{k-1}}{\left|l_{p,q}\right|}^2{d}_q,$$

(29)

where $l_{p,q}$ and $d_q$ denote the partial Cholesky quantities associated with the current selected set; equivalently, (29) is the Schur complement associated with ${\mathcal{P}}_{k-1}$. The next index is chosen over $p\notin {\mathcal{P}}_{k-1}$ as

$${\pi }_k^{\mathrm{SQRD}}=\arg \underset{p}{\min }{d}_{p\mid {\mathcal{P}}_{k-1}},{\pi }_k^{\mathrm{CPQR}}=\arg \underset{p}{\max }{d}_{p\mid {\mathcal{P}}_{k-1}}.$$

(30)

where ties are resolved by the smallest physical index. Since the FB recursion proceeds from stage N to stage 1, SQRD reserves stronger pivots for the earliest decisions, whereas CPQR applies the opposite allocation.

To investigate the statistical properties of SQRD and CPQR, Figure 2 compares the cumulative partial sums and profile variations. The nonintersecting expected-gap curves motivate the following statistical majorization hypothesis:

$$\mathbb{E}\left[{\delta }_{\mathrm{SQRD}}\right]\prec \mathbb{E}\left[{\delta }_{\pi }\right]\prec \mathbb{E}\left[{\delta }_{\mathrm{CPQR}}\right].$$

(31)

The annotations $(\mathrm{\Psi },\mathrm{\Phi },\mathrm{\Sigma })$ in Figure 2 are the corresponding un profile statistics, with $\mathrm{\Sigma }=\mathrm{\Psi }+\mathrm{\Phi }$. They provide descriptive motivation for (31), rather than a realization-wise proof, because sorting and expectation do not commute. The same trend motivates $\mathrm{\Phi }\left({\delta }_{\mathrm{SQRD}}\right)\le \mathrm{\Phi }\left({\delta }_{\pi }\right)\le \mathrm{\Phi }\left({\delta }_{\mathrm{CPQR}}\right)$.

Proposition 1

(Statistical Proxy Bounds). Suppose (31) and the stated variance ordering hold, with $C_2>0$ for almost every $\overline{\gamma }$. Averaging the resulting inequalities over $\overline{\gamma }$ gives

$$\mathbb{E}\left[{\widetilde{\mathrm{BER}}}_{\mathrm{SQRD}}^{\mathrm{GA}}\right]\le \mathbb{E}\left[{\widetilde{\mathrm{BER}}}_{\pi }^{\mathrm{GA}}\right]\le \mathbb{E}\left[{\widetilde{\mathrm{BER}}}_{\mathrm{CPQR}}^{\mathrm{GA}}\right].$$

(32)

Proof. 

See Appendix A.    □

Thus SQRD and CPQR serve as the statistical performance references of P-DFE under the stated proxy model, rather than universal channel-by-channel BER bounds.

4.2. Fast Graph-Based Ordering Approximation

The complexity of standard SQRD is on the order of $\mathcal{O}\left(N^3\right)$, which scales poorly with large block sizes. To mitigate this, we propose the EG ordering algorithm shown in Figure 3. EG replaces dense QR pivoting with reverse elimination on a pruned energy graph. As shown by the three steps in Figure 3, the procedure consists of energy-graph and heap initialization, node selection, and neighbor-key update.

4.2.1. Energy-Graph and Heap Initialization

Let ${\mathbf{G}}_0={\mathbf{H}}_{\mathrm{eff}}^H{\mathbf{H}}_{\mathrm{eff}}$ denote the unregularized Gram matrix over $\mathcal{S}$, with ${\left[{\mathbf{G}}_0\right]}_{p,p}>0$ for each active symbol. For $p\ne q$, define the normalized coupling weight ${\alpha }_{p,q}\triangleq {\left|{\left[{\mathbf{G}}_0\right]}_{p,q}\right|}^2/\left({\left[{\mathbf{G}}_0\right]}_{p,p}{\left[{\mathbf{G}}_0\right]}_{q,q}\right)$. The pruned energy graph is

$${\mathcal{G}}_{\mathrm{EG}}\left({\delta }_{\mathrm{EG}}\right)\triangleq \left(\mathcal{S},\mathcal{E}\left({\delta }_{\mathrm{EG}}\right)\right),$$

(33)

where $(q,p)\in \mathcal{E}\left({\delta }_{\mathrm{EG}}\right)$ iff $q\ne p$, $|{\left[{\mathbf{G}}_0\right]}_{p,q}|>0$, and ${\alpha }_{p,q}>{\delta }_{\mathrm{EG}}$; the edge means that q is retained as a dominant interference source for p. Thus ${\delta }_{\mathrm{EG}}=0$ retains nonzero couplings rather than all $N(N-1)$ pairs, while fractional leakage requires a positive threshold or explicit support truncation for sparsity. The retained energy and average directed degree are

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{E}}\left({\delta }_{\mathrm{EG}}\right)& ={\displaystyle \frac{{\sum }_{(q,p)\in \mathcal{E}\left({\delta }_{\mathrm{EG}}\right)}{\left|{\left[{\mathbf{G}}_0\right]}_{p,q}\right|}^2}{{\sum }_{\begin{array}{c}q,p\in \mathcal{S}\\ q\ne p\end{array}}{\left|{\left[{\mathbf{G}}_0\right]}_{p,q}\right|}^2}},\hfill \\ \hfill D_{\mathrm{avg}}\left({\delta }_{\mathrm{EG}}\right)& ={\displaystyle \frac{\left|\mathcal{E}\right({\delta }_{\mathrm{EG}}\left)\right|}{N}}.\hfill \end{array}$$

(34)

where ${\rho }_{\mathrm{E}}$ is used only when its denominator is nonzero. Once ${\delta }_{\mathrm{EG}}$ is fixed, we write ${\mathcal{G}}_{\mathrm{EG}}=(\mathcal{S},\mathcal{E})$ and define ${\mathcal{N}}_{\mathrm{in}}\left(p\right)=\{q:(q,p)\in \mathcal{E}\}$ and ${\mathcal{N}}_{\mathrm{out}}\left(p\right)=\{q:(p,q)\in \mathcal{E}\}$. For a remaining symbol set ${\mathcal{S}}_{\mathrm{rem},t}\subseteq \mathcal{S}$, the current residual interference power and heap key are

$$J_p^{\left(t\right)}\triangleq \sum _{\begin{array}{c}q\in {\mathcal{N}}_{\mathrm{in}}\left(p\right)\\ q\in {\mathcal{S}}_{\mathrm{rem},t}\end{array}}{\displaystyle \frac{|{\left[{\mathbf{G}}_0\right]}_{p,q}{|}^2}{{\left[{\mathbf{G}}_0\right]}_{p,p}}},{\vartheta }_p^{\left(t\right)}\triangleq {\displaystyle \frac{E_s{\left[{\mathbf{G}}_0\right]}_{p,p}}{{\sigma }^2+E_s{J}_p^{\left(t\right)}}}.$$

(35)

The heap is initialized with ${\mathcal{S}}_{\mathrm{rem},1}=\mathcal{S}$ and keys ${\left\{{\vartheta }_p^{\left(1\right)}\right\}}_{p\in \mathcal{S}}$.

4.2.2. Node Selection

At iteration t, the residual graph is the induced graph ${\mathcal{G}}_t=({\mathcal{S}}_{\mathrm{rem},t},{\mathcal{E}}_t)$, where ${\mathcal{E}}_t=\{(q,p)\in \mathcal{E}:q,p\in {\mathcal{S}}_{\mathrm{rem},t}\}$. EG extracts the maximum-key symbol from the heap and writes it into the reverse detection order:

$$m_t=\arg \underset{p\in {\mathcal{S}}_{\mathrm{rem},t}}{\max }{\vartheta }_p^{\left(t\right)},{\pi }_{\star ,N-t+1}^{\mathrm{EG}}=m_t,$$

(36)

then removes it by setting ${\mathcal{S}}_{\mathrm{rem},t+1}={\mathcal{S}}_{\mathrm{rem},t}\setminus \left\{m_t\right\}$. Ties are resolved by the smallest physical index.

4.2.3. Neighbor Update

After $m_t$ is removed, only its remaining out-neighbors lose one residual interferer. For $q\in {\mathcal{N}}_{\mathrm{out}}\left(m_t\right)\cap {\mathcal{S}}_{\mathrm{rem},t+1}$,

$$J_q^{(t+1)}=\max \left\{0,J_q^{\left(t\right)}-{\displaystyle \frac{|{\left[{\mathbf{G}}_0\right]}_{q,m_t}{|}^2}{{\left[{\mathbf{G}}_0\right]}_{q,q}}}\right\}.$$

(37)

and the updated key is ${\vartheta }_q^{(t+1)}=E_s{\left[{\mathbf{G}}_0\right]}_{q,q}/({\sigma }^2+E_s{J}_q^{(t+1)})$. The procedure then returns to Step 2) Node Selection and repeats until all symbols are removed. The complete EG ordering flow is summarized in Algorithm 1.

Algorithm 1 Energy-Graph (EG) Ordering Algorithm

Input: Effective channel ${\mathbf{H}}_{\mathrm{eff}}$, noise variance ${\sigma }^2$, symbol energy $E_s$, pruning threshold ${\delta }_{\mathrm{EG}}$. Output: EG order ${\pi }_{\star }^{\mathrm{EG}}$ and permutation matrix ${\mathbf{P}}_{\star }$.   1: ${\mathbf{G}}_0\leftarrow {\mathbf{H}}_{\mathrm{eff}}^H{\mathbf{H}}_{\mathrm{eff}}$   2: Build ${\mathcal{G}}_{\mathrm{EG}}=(\mathcal{S},\mathcal{E})$ and neighbor lists using (33)   3: for all $p\in \mathcal{S}$ do   4:     $J_p\leftarrow {\sum }_{q\in {\mathcal{N}}_{\mathrm{in}}\left(p\right)}{\left|{\left[{\mathbf{G}}_0\right]}_{p,q}\right|}^2/{\left[{\mathbf{G}}_0\right]}_{p,p}$   5:     ${\vartheta }_p\leftarrow E_s{\left[{\mathbf{G}}_0\right]}_{p,p}/({\sigma }^2+E_s{J}_p)$   6: end for   7: Insert all $(p,{\vartheta }_p)$ into a max-priority queue $\mathcal{H}$   8: ${\mathcal{S}}_{\mathrm{rem}}\leftarrow \mathcal{S}$   9: for $t=1$ to N do 10:     $m_t\leftarrow \mathrm{ExtractMax}\left(\mathcal{H}\right)$ (smallest-index tie break) 11:     ${\pi }_{\star ,N-t+1}^{\mathrm{EG}}\leftarrow m_t$, ${\mathcal{S}}_{\mathrm{rem}}\leftarrow {\mathcal{S}}_{\mathrm{rem}}\setminus \left\{m_t\right\}$ 12:     for all $q\in {\mathcal{N}}_{\mathrm{out}}\left(m_t\right)\cap {\mathcal{S}}_{\mathrm{rem}}$ do 13:        $J_q\leftarrow \max \{0,J_q-|{\left[{\mathbf{G}}_0\right]}_{q,m_t}{|}^2/{\left[{\mathbf{G}}_0\right]}_{q,q}\}$ 14:        ${\vartheta }_q\leftarrow E_s{\left[{\mathbf{G}}_0\right]}_{q,q}/({\sigma }^2+E_s{J}_q)$ 15:        ${key}_{\mathcal{H}}\left(q\right)\leftarrow {\vartheta }_q$ 16:     end for 17: end for 18: ${\mathbf{P}}_{\star }\leftarrow {\left[\mathbb{I}\{u={\pi }_{\star ,k}^{\mathrm{EG}}\}\right]}_{k,u}$ 19: return ${\pi }_{\star }^{\mathrm{EG}}$ and ${\mathbf{P}}_{\star }$

Given ${\mathbf{G}}_0$ and the retained edge lists, initialization requires $\mathcal{O}(N+|\mathcal{E}\left|\right)$ operations, heap extraction requires $\mathcal{O}\left(N\log N\right)$, and neighbor updates require $\mathcal{O}\left(\right|\mathcal{E}\left|\log N\right)$. Each node is extracted exactly once and each directed edge triggers at most one key update. Therefore,

$${\mathcal{C}}_{\mathrm{order}}=\mathcal{O}\left((N+|\mathcal{E}\left|\right)\log N\right).$$

(38)

Figure 4 shows that, in the evaluated setting, EG tracks the SQRD reference across FF disturbance, FB error, and reconstructed-SINR profiles while avoiding repeated global pivot searches.

5. Decoupling Error Propagation via Parallel P-MBDFE Scheme

Although ordering enhances the performance of a single-branch equalizer, the P-DFE remains highly vulnerable to error propagation under severe channel conditions or high-order modulations. To address this limitation, we extend reordering beyond the optimization of a single-branch perspective; instead, we leverage permutation to construct a family of multiple branches that simultaneously exhibit high performance and reduced shared feedback topology. Guided by the structural symmetry of AFDM channels and the KT distance metric, we propose the P-MBDFE architecture to achieve these design objectives.

5.1. P-MBDFE Roadmap

Figure 5 summarizes the proposed P-MBDFE receiver as a continuation of the single-branch ordering framework developed in Section 3 and Section 4.2. The EG module first provides a high-quality base order ${\pi }_{\star }^{\mathrm{EG}}$ for a single P-DFE branch. Starting from this base order, the multi-branch construction module generates a set of candidate permutations by exploiting AFDM structural symmetries, while the offline KT-based selection path identifies the optimal branch set ${\mathcal{B}}^{\star }$. Each selected element is a permutation applied to an otherwise identical P-DFE detector, rather than a new detector state.

This roadmap also clarifies the organization of the remainder of this section. Section 5.2 first derives the fusion objective and shows that P-MBDFE must simultaneously preserve branch-wise reliability and suppress residual inter-branch correlation. Section 5.3 addresses the first requirement by constructing a quality-preserving shift/reversal branch family, and Section 5.4 addresses the second requirement by selecting branches with large KT distance. The selected branches are then executed in parallel, fused by maximum-ratio combining (MRC), and converted into bit log-likelihood ratios (LLRs) for decoding.

5.2. Fusion Gain of P-MBDFE and Optimization Objective

To suppress feedback sensitivity, P-MBDFE introduces B branches with distinct permutation matrices ${\left\{{\mathbf{P}}_{{\pi }_b}\right\}}_{b=1}^B$. For branch b, let $e_m^{\left(b\right)}\triangleq {\tilde{x}}_m^{\left(b\right)}-x_m$ and ${ϵ}_{b,m}^2\triangleq \mathbb{E}\left[\right|e_m^{\left(b\right)}{|}^2]$; in implementation, this branch mean-square error (MSE) is estimated by ${\sigma }_{\mathrm{f},{\pi }_b,k}^2+{\sigma }_{\mathrm{b},{\pi }_b,k}^2$ at $k={\varphi }_{{\pi }_b}\left(m\right)$ using (18) and (22). Since the full cross-branch covariance is unavailable online, the branch MSEs form the inverse-variance weights

$${\omega }_{b,m}\triangleq {\displaystyle \frac{{ϵ}_{b,m}^{-2}}{{\sum }_{c=1}^B{ϵ}_{c,m}^{-2}}},{\tilde{x}}_m={\displaystyle \sum _{b=1}^B}{\omega }_{b,m}{\tilde{x}}_m^{\left(b\right)}.$$

(39)

For two branches, define the physical symbol-wise error correlation

$${\rho }_m({\pi }_b,{\pi }_c)\triangleq {\displaystyle \frac{Re\left\{\mathbb{E}\left[e_m^{\left(b\right)}{\left(e_m^{\left(c\right)}\right)}^{*}\right]\right\}}{{ϵ}_{b,m}{ϵ}_{c,m}}}\in [-1,1].$$

(40)

With $S_m\triangleq {\sum }_b{ϵ}_{b,m}^{-2}$, expanding $\mathbb{E}\left[\right|{\sum }_b{\omega }_{b,m}{e}_m^{\left(b\right)}{|}^2]$ and substituting (39) and (40) gives

$${ϵ}_{\mathrm{fuse},m}^2={\displaystyle \frac{1}{S_m}}\left[1+\underset{\triangleq {\mathrm{\Delta }}_{\mathrm{corr},m}}{\underbrace{{\displaystyle \frac{2}{S_m}}{\sum }_{b<c}{\displaystyle \frac{{\rho }_m({\pi }_b,{\pi }_c)}{{ϵ}_{b,m}{ϵ}_{c,m}}}}}\right].$$

(41)

For equiprobable constellation labels, the matched Max-Log LLR under the standard Gaussian residual model is

$${\mathrm{\Lambda }}_{m,j}={\displaystyle \frac{1}{{ϵ}_{\mathrm{fuse},m}^2}}\left(\underset{a\in {\mathcal{A}}_j^{\left(1\right)}}{\min }|{\tilde{x}}_m{-a|}^2-\underset{a\in {\mathcal{A}}_j^{\left(0\right)}}{\min }{|{\tilde{x}}_m-a|}^2\right),$$

(42)

where ${\mathcal{A}}_j^{\left(0\right)}$ and ${\mathcal{A}}_j^{\left(1\right)}$ denote the subsets of constellation points with the j-th bit equal to 0 and 1, respectively. A smaller actual variance in (41) makes the fused observation more informative and improves the resulting LLR quality. The actual ${\rho }_m$, however, is unavailable online, so the matched variance in (42) cannot be evaluated directly; merely replacing it by an underestimated value would inflate $|{\mathrm{\Lambda }}_{m,j}|$ and produce overconfident soft information.

Improving the fused observation and calibrating its LLRs are distinct but successive objectives addressed in three stages. First, Section 5.3 preserves a large $S_m$ through quality-preserving branch generation. Second, Section 5.4 links shared feedback ancestry to KT overlap and maximizes harmonic KT separation, enabling branch selection without unavailable symbol-wise correlations and yielding a provable envelope for the correlation penalty. Third, after branch selection, a single scale factor sets the correlation magnitude used in the fusion-variance estimate and LLR computation.

5.3. Quality-Preserving Branch Generation

We first preserve branch-wise reliability before optimizing cross-branch diversity. Equation (41) shows that this requires a large $S_m={\sum }_b{ϵ}_{b,m}^{-2}$, with each reliable branch contributing through its inverse MSE. Permutation diversity alone does not ensure this: an arbitrary detection order can create weak stages and severe error propagation, yielding a large ${ϵ}_{b,m}^2$ and therefore only a small contribution to $S_m$. At $k={\varphi }_{{\pi }_b}\left(m\right)$, the fusion model approximates the branch MSE as ${ϵ}_{b,m}^2\approx {\sigma }_{\mathrm{f},{\pi }_b,k}^2+{\sigma }_{\mathrm{b},{\pi }_b,k}^2$. Combining this approximation with (24) gives ${ϵ}_{b,m}^{-2}\approx {\mathrm{SINR}}_{{\pi }_b,k}/E_s$, so branches with higher stage-wise SINR contribute more strongly to $S_m$.

The EG base order therefore serves as the reliability reference for branch generation. To obtain distinct decision trajectories without sacrificing this reliability, we exploit the near-circulant structure of its ordered Gram matrix. Cyclic shifts and reversal approximately preserve the stage-wise SINR profile and hence the branch contributions to $S_m$, yielding the candidate family

$$\mathcal{F}\left({\mathbf{P}}_{\star }\right)\triangleq \left\{{\mathrm{\Pi }}^s{\mathbf{P}}_{\star },\mathbf{J}{\mathrm{\Pi }}^s{\mathbf{P}}_{\star }\mid s\in \{0,\dots ,N-1\}\right\},$$

(43)

where $\mathrm{\Pi }$ is the cyclic-shift matrix and ${\left[\mathbf{J}\right]}_{k,\mathcal{l}}=\mathbb{I}\{k+\mathcal{l}=N+1\}$ for one-based indices. Each element of $\mathcal{F}\left({\mathbf{P}}_{\star }\right)$ is a complete branch permutation matrix. The following proposition formalizes its SINR quasi-invariance relative to the EG base order.

Proposition 2

(SINR Quasi-Invariance via Structure Preservation). Let the base-ordered regularized Gram matrix be ${\mathbf{G}}_{\star }={\mathbf{P}}_{\star }\mathbf{G}{\mathbf{P}}_{\star }^T={\mathbf{C}}_{\star }+{\mathbf{R}}_{\star }$, where ${\mathbf{C}}_{\star }$ is the dominant Hermitian circulant component and ${\mathbf{R}}_{\star }$ is the residual. For any $s\in \{0,\dots ,N-1\}$, define ${\mathbf{P}}_{\mathrm{shift}}={\mathrm{\Pi }}^s{\mathbf{P}}_{\star }$ and ${\mathbf{P}}_{\mathrm{rev}}=\mathbf{J}{\mathrm{\Pi }}^s{\mathbf{P}}_{\star }$. Their Gram matrices satisfy ${\mathbf{G}}_{\mathrm{shift}}={\mathbf{G}}_{\star }+{\mathbf{V}}_{\mathrm{shift}}$ and ${\mathbf{G}}_{\mathrm{rev}}={\mathbf{G}}_{\star }^{*}+{\mathbf{V}}_{\mathrm{rev}}$, with $\max \{\parallel {\mathbf{V}}_{\mathrm{shift}}{\parallel }_F,\parallel {\mathbf{V}}_{\mathrm{rev}}{\parallel }_F\}\le 2\parallel {\mathbf{R}}_{\star }{\parallel }_F$. If ${\lambda }_{\min }\left({\mathbf{G}}_{\star }\right)\ge {\lambda }_0>0$ and $\parallel {\mathbf{R}}_{\star }{\parallel }_F$ is sufficiently small, the stage-wise GA SINR profiles satisfy

$$\underset{\begin{array}{c}b\in \{\mathrm{shift},\mathrm{rev}\}\\ 1\le k\le N\end{array}}{\max }\left|{\mathrm{SINR}}_{b,k}^{\mathrm{GA}}-{\mathrm{SINR}}_{\star ,k}^{\mathrm{GA}}\right|=\mathcal{O}\left(\parallel {\mathbf{R}}_{\star }{\parallel }_F\right).$$

(44)

Proof. 

See Appendix B.    □

Thus, (44) establishes SINR-profile quasi-invariance, and Figure 6 confirms that the shifted and reversed branches closely track the EG base-order profile.

5.4. Diversity Selection via Kendall–Tau Distance

After preserving branch quality, the remaining design objective in (41) is to reduce ${\mathrm{\Delta }}_{\mathrm{corr},m}$. Because the exact symbol-wise correlations are unavailable during offline branch selection, we seek a permutation-only characterization of their order dependence.

To construct this characterization, we first express each branch error in the physical-symbol domain. Figure 7 illustrates the detected/undetected partition induced by each order through the rank function ${\varphi }_{\pi }$.

For a physical symbol $m\in \mathcal{S}$ under permutation $\pi$, define the undetected set${\mathcal{U}}_{\pi }\left(m\right)$ (symbols not yet detected when m is processed) and the detected set${\mathcal{D}}_{\pi }\left(m\right)$ (symbols already detected and used for cancellation when m is processed) as

$$\begin{array}{cc}\hfill {\mathcal{U}}_{\pi }\left(m\right)& \triangleq \{p\in \mathcal{S}\setminus \left\{m\right\}\mid {\varphi }_{\pi }\left(p\right)<{\varphi }_{\pi }\left(m\right)\},\hfill \end{array}$$

(45a)

$$\begin{array}{cc}\hfill {\mathcal{D}}_{\pi }\left(m\right)& \triangleq \{p\in \mathcal{S}\setminus \left\{m\right\}\mid {\varphi }_{\pi }\left(p\right)>{\varphi }_{\pi }\left(m\right)\}.\hfill \end{array}$$

(45b)

Based on (17), the soft-decision error $e_{\pi ,k}\triangleq {\tilde{x}}_{\pi ,k}-x_{\pi ,k}$ at stage k equals $n_{\mathrm{f},\pi ,k}$ plus weighted contributions from previously detected decision errors:

$$e_{\pi ,k}=n_{\mathrm{f},\pi ,k}+{\displaystyle \frac{1}{{\mu }_{\pi ,k}}}{\displaystyle \sum _{i=k+1}^N}{l}_{i,k}^{\left(\pi \right)*}\mathrm{\Delta }{x}_{\pi ,i}.$$

(46)

For a physical symbol $m\in \mathcal{S}$, let $k={\varphi }_{\pi }\left(m\right)$. Remapping (46) with the detected set in (45b) gives

$$e_m^{\left(\pi \right)}=n_{\mathrm{f},m}^{\left(\pi \right)}+\sum _{p\in {\mathcal{D}}_{\pi }\left(m\right)}{\zeta }_{m,p}^{\left(\pi \right)}\mathrm{\Delta }{x}_p^{\left(\pi \right)},$$

(47)

where $n_{\mathrm{f},m}^{\left(\pi \right)}$ denotes $n_{\mathrm{f},\pi ,k}$ after remapping to physical symbol m, ${\zeta }_{m,p}^{\left(\pi \right)}\triangleq l_{{\varphi }_{\pi }\left(p\right),{\varphi }_{\pi }\left(m\right)}^{\left(\pi \right)*}/{\mu }_{\pi ,{\varphi }_{\pi }\left(m\right)}$ is the order-dependent transfer coefficient from $\mathrm{\Delta }{x}_p^{\left(\pi \right)}$ to $e_m^{\left(\pi \right)}$, and $\mathrm{\Delta }{x}_p^{\left(\pi \right)}$ is the hard-decision error at physical symbol p. Expanding (47) for two branches $\pi$ and ${\pi }^{\prime }$, while neglecting cross terms between the $n_{\mathrm{f},m}$ terms and the decision errors, gives

$$\begin{array}{cc}\hfill C_m(\pi ,{\pi }^{\prime })& \triangleq \mathbb{E}\left[e_m^{\left(\pi \right)}{\left(e_m^{\left({\pi }^{\prime }\right)}\right)}^{*}\right]\hfill \\ \hfill & =C_m^{\left(0\right)}+\sum _{p\in {\mathcal{D}}_{\pi }\left(m\right)}\sum _{q\in {\mathcal{D}}_{{\pi }^{\prime }}\left(m\right)}{\zeta }_{m,p}^{\left(\pi \right)}{\left({\zeta }_{m,q}^{\left({\pi }^{\prime }\right)}\right)}^{*}{C}_{p,q}^{\mathrm{err}}\hfill \\ \hfill & \approx C_m^{\left(0\right)}+\sum _{p\in {\mathcal{D}}_m^{\cap }}{\zeta }_{m,p}^{\left(\pi \right)}{\left({\zeta }_{m,p}^{\left({\pi }^{\prime }\right)}\right)}^{*}{C}_p^{\mathrm{err}},\hfill \end{array}$$

(48)

where $C_m^{\left(0\right)}\triangleq \mathbb{E}\left[n_{\mathrm{f},m}^{\left(\pi \right)}{\left(n_{\mathrm{f},m}^{\left({\pi }^{\prime }\right)}\right)}^{*}\right]$ is the cross-covariance between the two remapped $n_{\mathrm{f},m}$ terms, $C_{p,q}^{\mathrm{err}}\triangleq \mathbb{E}\left[\mathrm{\Delta }{x}_p^{\left(\pi \right)}{\left(\mathrm{\Delta }{x}_q^{\left({\pi }^{\prime }\right)}\right)}^{*}\right]$, $C_p^{\mathrm{err}}\equiv C_{p,p}^{\mathrm{err}}$, and ${\mathcal{D}}_m^{\cap }\triangleq {\mathcal{D}}_{\pi }\left(m\right)\cap {\mathcal{D}}_{{\pi }^{\prime }}\left(m\right)$. The last line retains contributions from decision errors associated with symbols detected under both orders and neglects cross-symbol decision-error moments. Normalizing its real part according to (40) gives

$${\rho }_m({\pi }_b,{\pi }_c)\approx {\rho }_{0,m}^{(b,c)}+\sum _{p\in {\mathcal{D}}_m^{\cap }}{a}_{m,p}^{(b,c)},$$

(49)

where ${\rho }_{0,m}^{(b,c)}\triangleq Re\left\{C_m^{\left(0\right)}\right\}/\left({ϵ}_{b,m}{ϵ}_{c,m}\right)$ and $a_{m,p}^{(b,c)}\triangleq Re\left\{{\zeta }_{m,p}^{\left({\pi }_b\right)}{\left({\zeta }_{m,p}^{\left({\pi }_c\right)}\right)}^{*}{C}_p^{\mathrm{err}}\right\}/\left({ϵ}_{b,m}{ϵ}_{c,m}\right)$. For any unordered symbol pair $\{u,v\}$, the two detection orders either preserve or reverse its relative order. In the former case, exactly one orientation of the pair appears in the detected-set intersections; in the latter, neither orientation appears and the pair contributes one inversion. This one-to-one correspondence motivates the Kendall–Tau distance.

Definition 2

(Kendall–Tau Distance). For two permutations π and ${\pi }^{\prime }$, define

$$\begin{array}{cc}\hfill d_{\mathrm{KT}}(\pi ,{\pi }^{\prime })& \triangleq \sum _{u<v}\mathbb{I}\{[{\varphi }_{\pi }\left(u\right)-{\varphi }_{\pi }\left(v\right)]\hfill \\ & \times [{\varphi }_{{\pi }^{\prime }}\left(u\right)-{\varphi }_{{\pi }^{\prime }}\left(v\right)]<0\},\hfill \end{array}$$

which counts their pairwise order inversions [30].

Because permutations have no rank ties, this correspondence gives ${\sum }_{m\in \mathcal{S}}|{\mathcal{D}}_{\pi }\left(m\right)\cap {\mathcal{D}}_{{\pi }^{\prime }}\left(m\right)|=\left({\displaystyle \genfrac{}{}{0pt}{}{N}{2}}\right)-d_{\mathrm{KT}}(\pi ,{\pi }^{\prime })$. Dividing by the total number of symbol pairs gives

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{KT}}(\pi ,{\pi }^{\prime })& \triangleq {\displaystyle \frac{{\sum }_{m\in \mathcal{S}}|{\mathcal{D}}_{\pi }\left(m\right)\cap {\mathcal{D}}_{{\pi }^{\prime }}\left(m\right)|}{\left(\genfrac{}{}{0pt}{}{N}{2}\right)}}\hfill \\ \hfill & =1-{\displaystyle \frac{d_{\mathrm{KT}}(\pi ,{\pi }^{\prime })}{\left(\genfrac{}{}{0pt}{}{N}{2}\right)}}.\hfill \end{array}$$

(50)

Using ${\rho }_{\mathrm{KT}}$ in the correlation penalty of (41), while retaining the branch-specific inverse-MSE weights, gives

$${\mathrm{\Delta }}_{\mathrm{KT},m}\left(\mathcal{B}\right)\triangleq {\displaystyle \frac{2}{S_m}}\sum _{b<c}{\displaystyle \frac{{\rho }_{\mathrm{KT}}({\pi }_b,{\pi }_c)}{{ϵ}_{b,m}{ϵ}_{c,m}}}.$$

(51)

Equation (51) uses ${\rho }_{\mathrm{KT}}$ in place of the unavailable pairwise correlations while preserving the inverse-MSE weights in (41).

To penalize highly similar branch pairs, we use the harmonic mean of their KT distances. Its relation to ${\mathrm{\Delta }}_{\mathrm{KT},m}$ is established below.

Proposition 3

(Harmonic KT Branch Selection). For a B-branch set $\mathcal{B}$, let $M\triangleq \left({\displaystyle \genfrac{}{}{0pt}{}{B}{2}}\right)$, $C_N\triangleq \left({\displaystyle \genfrac{}{}{0pt}{}{N}{2}}\right)$, and $d_{bc}\triangleq d_{\mathrm{KT}}({\pi }_b,{\pi }_c)$. Define

$$d_{\mathrm{H}}\left(\mathcal{B}\right)\triangleq {\displaystyle \frac{M}{{\displaystyle \sum _{b<c}{d}_{bc}^{-1}}}},$$

(52)

with $d_{\mathrm{H}}\left(\mathcal{B}\right)=0$ if any $d_{bc}=0$, and let $d_{\mathrm{A}}\left(\mathcal{B}\right)\triangleq M^{-1}{\sum }_{b<c}{d}_{bc}$. Then $d_{\mathrm{H}}\left(\mathcal{B}\right)\le d_{\mathrm{A}}\left(\mathcal{B}\right)$, with equality if and only if all pairwise distances are equal, and, for every symbol m,

$${\mathrm{\Delta }}_{\mathrm{KT},m}\left(\mathcal{B}\right)\le M\left(1-{\displaystyle \frac{d_{\mathrm{H}}\left(\mathcal{B}\right)}{C_N}}\right).$$

Hence maximizing $d_{\mathrm{H}}$ exactly minimizes the aggregate inverse-distance risk ${\sum }_{b<c}{d}_{bc}^{-1}$, balances the weakest branch pairs, and minimizes the stated reliability-independent envelope of ${\mathrm{\Delta }}_{\mathrm{KT},m}$.

Proof. 

See Appendix C.    □

Proposition 3 establishes $d_{\mathrm{H}}$ as the robust criterion for reducing ${\mathrm{\Delta }}_{\mathrm{KT},m}$. The offline selected subset is therefore

$${\mathcal{B}}_B^{\star }=\arg \underset{\begin{array}{c}\mathcal{B}\subset \mathcal{F}\left({\mathbf{P}}_{\star }\right)\\ \left|\mathcal{B}\right|=B\end{array}}{\max }{d}_{\mathrm{H}}\left(\mathcal{B}\right).$$

(53)

Figure 8 supports the inverse relation between $d_{\mathrm{H}}$ and the actual cross-branch correlation over the evaluated branch subsets. For the block-level approximation, we neglect the average of ${\rho }_{0,m}^{(b,c)}$ in (49) and replace $a_{m,p}^{(b,c)}$ by a common value $\overline{a}>0$ for $p\in {\mathcal{D}}_m^{\cap }$. Let ${\overline{\rho }}_{bc}\triangleq {\displaystyle \frac{1}{N}}{\sum }_m{\rho }_m({\pi }_b,{\pi }_c)$. Equations (49) and (50) then give

$${\overline{\rho }}_{bc}\approx \underset{\theta }{\underbrace{{\displaystyle \frac{N-1}{2}}\overline{a}}}{\rho }_{\mathrm{KT}}({\pi }_b,{\pi }_c).$$

(54)

Applying this block-level approximation to the correlation penalty gives

$${\widehat{\mathrm{\Delta }}}_{\mathrm{corr},m}\left(\theta \right)\triangleq {\displaystyle \frac{2}{S_m}}\sum _{b<c}{\displaystyle \frac{\theta {\rho }_{\mathrm{KT}}({\pi }_b,{\pi }_c)}{{ϵ}_{b,m}{ϵ}_{c,m}}}=\theta {\mathrm{\Delta }}_{\mathrm{KT},m}.$$

The variance supplied to the soft demodulator is therefore

$${\widehat{ϵ}}_{\mathrm{fuse},m}^2\triangleq {\displaystyle \frac{1}{S_m}}\left(1+\theta {\mathrm{\Delta }}_{\mathrm{KT},m}\right).$$

(55)

Thus, $d_{\mathrm{H}}$ selects branches with less similar detection orders, whereas $\theta$ sets the correlation magnitude used in the fusion variance and hence the LLR scale.

With the offline branch set in (53) and the calibrated fusion variance in (55), the complete P-MBDFE receiver operates in three phases: 1) Branch Generation: Constructing ${\mathcal{B}}_B^{\star }$ from the EG-ordered base permutation; 2) Parallel Detection: Executing P-DFE for each branch to obtain soft symbols and reliability metrics; and 3) Fusion: Combining the outputs via the MRC rule in (39), evaluating the fusion variance via (51) and (55), and computing final LLRs for decoding via (42).

As a concrete implementation example of this framework with $B=4$, as derived in Appendix D, a KT-optimal P-QDFE configuration is

$${\mathcal{B}}_4^{\mathrm{P}-\mathrm{QDFE}}=\left\{{\mathbf{P}}_{\star },{\mathrm{\Pi }}^{{\displaystyle \frac{N}{2}}}{\mathbf{P}}_{\star },\mathbf{J}{\mathbf{P}}_{\star },\mathbf{J}{\mathrm{\Pi }}^{{\displaystyle \frac{N}{2}}}{\mathbf{P}}_{\star }\right\}.$$

(56)

This specific instance, referred to as the P-QDFE scheme, comprises the base order, its semi-circle shift, and their time-reversed counterparts, as verified in Figure 8. The complete P-MBDFE procedure is summarized in Algorithm 2.

Algorithm 2 Proposed P-MBDFE Receiver

Input: ${\mathbf{H}}_{\mathrm{eff}}$, $\mathbf{y}$, ${\sigma }^2$, $E_s$, $\mathcal{A}$, B, ${\delta }_{\mathrm{EG}}$, ${\delta }_{\mathrm{ICHOL}}$, and $\theta$. Output: $\tilde{\mathbf{x}}$, $\widehat{\mathbf{x}}$, and $\left\{{\mathrm{\Lambda }}_{m,j}\right\}$.   1: $({\pi }_{\star }^{\mathrm{EG}},{\mathbf{P}}_{\star })\leftarrow EG({\mathbf{H}}_{\mathrm{eff}},{\sigma }^2,E_s,{\delta }_{\mathrm{EG}})$ (Alg. 1)   2: Construct $\mathcal{F}\left({\mathbf{P}}_{\star }\right)$ using (43)   3: Obtain ${\mathcal{B}}_B^{\star }$ from the offline solution of (53)   4: Compute ${\left\{{\rho }_{\mathrm{KT}}({\pi }_b,{\pi }_c)\right\}}_{b<c}$ using (50)   5: $\lambda \leftarrow {\sigma }^2/E_s$, $\mathbf{G}\leftarrow {\mathbf{H}}_{\mathrm{eff}}^H{\mathbf{H}}_{\mathrm{eff}}+\lambda {\mathbf{I}}_N$, $\mathbf{r}\leftarrow {\mathbf{H}}_{\mathrm{eff}}^H\mathbf{y}$   6: for all ${\mathbf{P}}_{{\pi }_b}\in {\mathcal{B}}_B^{\star }$ in parallel do   7:    ${\mathbf{G}}_{{\pi }_b}\leftarrow {\mathbf{P}}_{{\pi }_b}\mathbf{G}{\mathbf{P}}_{{\pi }_b}^T$   8:    ${\mathbf{r}}_{{\pi }_b}\leftarrow {\mathbf{P}}_{{\pi }_b}\mathbf{r}={\mathbf{H}}_{{\pi }_b}^H{\mathbf{y}}_{{\pi }_b}$   9:    $({\mathbf{L}}_{{\pi }_b},{\mathbf{D}}_{{\pi }_b})\leftarrow \mathrm{ICHOL}({\mathbf{G}}_{{\pi }_b};{\delta }_{\mathrm{ICHOL}})$ 10:    ${\mathbf{z}}_{{\pi }_b}\leftarrow {\mathbf{L}}_{{\pi }_b}^{-1}{\mathbf{r}}_{{\pi }_b}$ 11:    for $k=N$ down to 1 do 12:       ${\breve{x}}_{{\pi }_b,k}\leftarrow z_{{\pi }_b,k}/d_k^{\left({\pi }_b\right)}-{\sum }_{i=k+1}^N{l}_{i,k}^{\left({\pi }_b\right)*}{\widehat{x}}_{{\pi }_b,i}$ 13:       ${\mu }_{{\pi }_b,k}\leftarrow 1-\lambda /d_k^{\left({\pi }_b\right)}$, ${\tilde{x}}_{{\pi }_b,k}\leftarrow {\breve{x}}_{{\pi }_b,k}/{\mu }_{{\pi }_b,k}$ 14:       ${\widehat{x}}_{{\pi }_b,k}\leftarrow \mathcal{Q}\left({\tilde{x}}_{{\pi }_b,k}\right)$, $m\leftarrow {\pi }_{b,k}$, ${\tilde{x}}_m^{\left(b\right)}\leftarrow {\tilde{x}}_{{\pi }_b,k}$ 15:       ${ϵ}_{b,m}^2\leftarrow {\sigma }_{\mathrm{f},{\pi }_b,k}^2+{\sigma }_{\mathrm{b},{\pi }_b,k}^2$ from (18) and (22) 16:    end for 17: end for 18: for $m=1$ to N do 19:    Compute ${\left\{{\omega }_{b,m}\right\}}_{b=1}^B$ and ${\tilde{x}}_m$ from (39) 20:    Compute ${\mathrm{\Delta }}_{\mathrm{KT},m}$ from (51) and ${\widehat{ϵ}}_{\mathrm{fuse},m}^2$ from (55) 21:    for $j=1$ to ${\log }_2\left|\mathcal{A}\right|$ do 22:      Compute ${\mathrm{\Lambda }}_{m,j}$ from (42) 23:    end for 24:    ${\widehat{x}}_m\leftarrow \mathcal{Q}\left({\tilde{x}}_m\right)$ 25: end for

5.5. Computational Complexity

The P-MBDFE scheme comprises four stages: shared initialization, base ordering, parallel detection, and fusion. The initialization exploits the sparse effective channel to compute the Gram matrix and matched filter, bounded by $\mathcal{O}\left(N{L}^2\right)$. The EG base ordering relies primarily on priority queue updates and has the complexity in (38). The proposed receiver uses a configurable ICHOL factorization [34] for each branch. Setting ${\delta }_{\mathrm{ICHOL}}=0$ retains the complete elimination graph and recovers the full $LD{L}^H$ factorization in (11), whereas ${\delta }_{\mathrm{ICHOL}}>0$ drops weak fill entries. Let $D_{\mathrm{chol}}\equiv D_{\mathrm{chol}}\left({\delta }_{\mathrm{ICHOL}}\right)$ denote the root-mean-square (RMS) number of retained nonzeros per factor column. The B branches then require $\mathcal{O}\left(BN{D}_{\mathrm{chol}}^2\right)$ operations for factorization and FF/FB processing. Finally, MRC fusion and Max-Log LLR extraction require $\mathcal{O}\left(N\right|\mathcal{A}\left|\right)$ operations. Since ${\rho }_{\mathrm{KT}}$ and $\theta$ are predefined, the LLR variance calibration in (55) introduces no additional computational complexity.

Therefore, the total sequential complexity is $\mathcal{O}\left(N\left(L^2+(1+{\displaystyle \frac{\left|\mathcal{E}\right|}{N}})\log N+B{D}_{\mathrm{chol}}^2+\left|\mathcal{A}\right|\right)\right)$. Moreover, the B-branch parallel execution confines the critical path to shared preprocessing, single-branch evaluation, and fusion. This bounds the parallel computational complexity order to $\mathcal{O}\left(N\left(L^2+(1+{\displaystyle \frac{\left|\mathcal{E}\right|}{N}})\log N+D_{\mathrm{chol}}^2+\left|\mathcal{A}\right|\right)\right)$, decoupling the execution latency from the number of branches.

6. Simulation Results

In this section, we evaluate the proposed equalizers under two channel models: multipath Rayleigh fading channels for stress testing, and the 3rd Generation Partnership Project (3GPP) tapped-delay-line C (TDL-C) profile [35] for dense high-mobility fractional delay-Doppler validation. For the Rayleigh fading channels, each path has equal average power $\mathbb{E}\left[\right|h_i{|}^2]=1/L$, where L is the number of paths, and the maximum normalized delay shift is ${\mathcal{l}}_{\max }=10$. Under integer Doppler, the path Doppler shifts ${\nu }_i$ are chosen uniformly from the integer set $[-{\nu }_{\max },{\nu }_{\max }]$ with ${\nu }_{\max }=4$. Under fractional Doppler, they follow the Jakes model ${\nu }_i={\nu }_{\max }\cos \left({\theta }_i\right)$ with ${\theta }_i\sim \mathcal{U}[-\pi ,\pi ]$ and ${\nu }_{\max }=4.0$. Unless otherwise stated, we set $N=512$ and $L_{\mathrm{cp}}=16$. The two-branch permutation-driven BDFE (P-BDFE) utilizes the EG base order and its reversal, whereas the four-branch P-QDFE utilizes the optimal configuration in (56). Coded curves use the 3GPP New Radio (NR) LDPC base graph 2 (BG2) protograph [36], lifting size $Z_c=256$, and rate matching.

Perfect CSI is assumed except for the TDL-C robustness test, where the receiver reconstructs ${\widehat{\mathbf{H}}}_{\mathrm{eff}}$ from perturbed path gains and delay-Doppler parameters [37,38]. For the imperfect-CSI sweep, define $\mathbf{h}={[h_1,\dots ,h_L]}^T$, ${\mathbf{e}}_h={[e_{h,1},\dots ,e_{h,L}]}^T$, and the path-gain normalized mean-square error (NMSE) $\eta \triangleq \mathbb{E}[\parallel {\mathbf{e}}_h{\parallel }_2^2\mid {\mathbf{h}]/\parallel \mathbf{h}\parallel }_2^2$, reported as $10{\log }_{10}\eta$ dB. Setting $P_{\mathrm{H}}={\parallel \mathbf{h}\parallel }_2^2$, generate ${\widehat{h}}_p=h_p+e_{h,p}$ with $e_{h,p}\sim \mathcal{CN}(0,\eta P_{\mathrm{H}}/L)$, ${\widehat{\mathcal{l}}}_p=\max (0,{\mathcal{l}}_p+\mathrm{\Delta }{\mathcal{l}}_p)$, and ${\widehat{\nu }}_p={\nu }_p+\mathrm{\Delta }{\nu }_p$. Here $\mathrm{\Delta }{\mathcal{l}}_p={\sigma }_{\mathcal{l}}\sqrt{\eta }{w}_{\mathcal{l},p}/a_p$, $\mathrm{\Delta }{\nu }_p={\sigma }_{\nu }\sqrt{\eta }{w}_{\nu ,p}/a_p$, $a_p=\left|h_p\right|/\sqrt{P_{\mathrm{H}}}$, and $w_{\mathcal{l},p},w_{\nu ,p}\sim \mathcal{N}(0,1)$.

Figure 9 evaluates the single-branch DFE under various orderings against their genie-aided bounds under 4-path fractional delay-Doppler channels. SQRD and CPQR bound the ensemble performance. We observe that: first, the low-complexity EG tightly tracks SQRD, yielding SNR gains over natural ordering of $0.57$ dB for QPSK at a BER of $1.0\times {10}^{-4}$ and $0.51$ dB for 16QAM at a BER of $1.7\times {10}^{-3}$ by deferring high-interference symbols to later stages to shape the stage-wise SINR profile; second, the gap between the actual and genie-aided EG curves widens under higher-order modulations, reflecting their stronger sensitivity to error propagation and motivating the multi-branch design evaluated next.

Figure 9b extends the same setting to P-BDFE and P-QDFE across constellations and orderings. Branch-diversity gain grows with constellation order, with all gains evaluated relative to the Bi-Natural baseline at the maximum evaluated SNR of the respective modulation. Under QPSK, the large Euclidean distance limits error propagation, yielding nearly identical performance across different equalizers. Under 8-ary phase-shift keying (8PSK), at the Bi-Natural BER of $1.20\times {10}^{-4}$ at 18 dB, Bi-EG and Quad-EG improve by $0.31$ dB and $0.53$ dB, respectively. For 16QAM, at the Bi-Natural BER of $3.18\times {10}^{-4}$ at 20 dB, Bi-EG and Quad-EG provide $1.20$ dB and $1.96$ dB gains. For 32-ary amplitude and phase-shift keying (32APSK), at the Bi-Natural BER of $6.44\times {10}^{-3}$ at 22 dB, the gains rise to $1.34$ dB and $3.32$ dB. This confirms KT-guided diversity and near-SQRD EG ordering.

Figure 9c evaluates the EG pruning rule in (33) using the BER permutation gain $G_{\mathrm{perm}}=({\mathrm{BER}}_{\mathrm{Nat}}-{\mathrm{BER}}_{\mathrm{EG}})/{\mathrm{BER}}_{\mathrm{Nat}}\times 100\%$, together with the retained coupling-energy ratio ${\rho }_{\mathrm{E}}\left({\delta }_{\mathrm{EG}}\right)$ and average graph degree $D_{\mathrm{avg}}\left({\delta }_{\mathrm{EG}}\right)$ defined in (34). For QPSK, the unpruned single-branch $G_{\mathrm{perm}}$ stays in a narrow range from $15.2\%$ to $17.6\%$ as L increases from 4 to 16, showing that the ordering gain is stable under low-order modulation. For 16QAM, $G_{\mathrm{perm}}$ rises from $39.2\%$ at $L=4$ and stays around $60\%$-$62\%$ for $L=8$-16. This indicates that EG remains robust and becomes more beneficial in richer multipath channels. At ${\delta }_{\mathrm{EG}}=0.01$, $D_{\mathrm{avg}}\left({\delta }_{\mathrm{EG}}\right)$ drops from 511 to $12.0$-$24.1$ while ${\rho }_{\mathrm{E}}\left({\delta }_{\mathrm{EG}}\right)$ still retains $52.5\%$-$72.7\%$ of the coupling energy. The resulting $G_{\mathrm{perm}}$ remains positive, ranging from $7.9\%$ to $14.2\%$ for QPSK and from $16.3\%$ to $38.9\%$ for 16QAM. Therefore, ${\delta }_{\mathrm{EG}}=0.01$ is used in the remaining simulations as a practical trade-off between graph sparsity and BER gain.

Figure 10 compares the uncoded BER performance of the proposed equalizers against DFE baselines and iterative benchmarks under ideal sparse integer delay-Doppler channels. Under QPSK, P-DFE yields a $0.48$ dB gain over DFE at $E_b/N_0=16$ dB, where the baseline BER is $6.15\times {10}^{-6}$. Meanwhile, due to the low modulation order, the multi-branch equalizers P-BDFE, P-QDFE, and BDFE exhibit similar performance, closely approaching the OAMP benchmark. Under 16QAM, the low-complexity MPA benchmark with $I=20$ iterations fails even in this ideal sparse channel, with damping factor $\beta =0.7$ stagnating near $3.5\times {10}^{-1}$ and $\beta =0.3$ still near $2.7\times {10}^{-2}$ at $E_b/N_0=20$ dB. In contrast, the proposed equalizers yield significant gains. P-DFE provides a $1.32$ dB gain over DFE at $E_b/N_0=20$ dB, where the baseline BER is $9.49\times {10}^{-4}$. Similarly, P-BDFE and P-QDFE achieve $1.31$ dB and $2.62$ dB gains over BDFE at $E_b/N_0=20$ dB, where the baseline BDFE BER is $8.46\times {10}^{-5}$, with P-QDFE reaching $7.8\times {10}^{-6}$ to closely approach the best evaluated OAMP configuration using $I=8$ and $\beta =0.3$, which achieves $4.4\times {10}^{-6}$. This confirms the efficacy of the proposed equalizers in ideal sparse channels.

Furthermore, Figure 11 investigates the equalizers with 16QAM under the dense fractional-delay and fractional-Doppler TDL-C channel and imperfect CSI settings summarized in Table 2. In Figure 11a, under perfect CSI, P-DFE achieves a BER of $2.12\times {10}^{-4}$ at $E_b/N_0=20$ dB. For multi-branch equalizers, P-BDFE yields a $0.93$ dB gain over OAMP using $I=4$ and $\beta =0.3$ at $E_b/N_0=20$ dB, where the baseline OAMP BER is $2.98\times {10}^{-5}$. Similarly, P-QDFE achieves a $0.47$ dB gain over OAMP using $I=8$ and $\beta =0.3$ at $E_b/N_0=20$ dB, where the baseline OAMP BER is $7.28\times {10}^{-6}$. In Figure 11b, the EG-ordered receivers retain their gain over the natural DFE/BDFE baselines up to a nominal path-gain NMSE of about $10{\log }_{10}\eta =-22$ dB, indicating that the EG ordering remains useful under moderate delay-Doppler mismatch for 16QAM. Across the full CSI-perturbation sweep, P-QDFE achieves the best BER, showing that the quad-branch architecture further improves robustness under non-ideal CSI.

Prior to coded BER evaluation, Table 3 and Table 4 validate the KT-based LLR variance calibration using the standard bit-wise normalized generalized mutual information (NGMI) computed from the LLRs in (42). As one practical selection criterion, Table 3 chooses $\theta =2.0$ by average NGMI, while Table 4 shows a $0.1579$ NGMI gain at $E_b/N_0=8$ dB over the uncorrected variance; the resulting P-QDFE soft output also exceeds unbiased LMMSE at every evaluated point.

Figure 12 illustrates the uncoded and coded BER performance under 16QAM over a 5-path fractional Doppler channel, utilizing a code rate of $R=0.625$, yielding a spectral efficiency of 2.50 bits/s/Hz, and the belief propagation (BP) decoder [39]. Under coded transmission, P-QDFE yields a $1.56$ dB gain over unbiased LMMSE at $E_b/N_0=15$ dB, where the baseline unbiased LMMSE BER is $1.99\times {10}^{-5}$. Notably, while the OAMP configuration with $I=15$, $\beta =0.7$, and an OAMP-specific LLR scaling factor of $0.2$, empirically optimized via parameter sweeping, achieves a lower uncoded BER, P-QDFE exhibits a steeper coded waterfall curve, yielding a $0.64$ dB gain over OAMP at $E_b/N_0=15$ dB, where the baseline OAMP BER is $2.81\times {10}^{-6}$. Although this scaling may not represent OAMP’s performance limit under a joint design, the comparison demonstrates the coded robustness of P-QDFE.

Figure 13 illustrates the coded BER-complexity trade-off under ICHOL, where ${\overline{D}}_{\mathrm{chol}}$ is the arithmetic-mean number of retained nonzeros per factor column; the asymptotic operation term uses the RMS degree defined above. At code rate $R=5/12$ and $E_b/N_0=11$ dB, increasing ${\delta }_{\mathrm{ICHOL}}$ reduces the average filter complexity. Under integer Doppler, ${\overline{D}}_{\mathrm{chol}}$ decreases from about $25.0$ to $14.4$ at ${\delta }_{\mathrm{ICHOL}}=1.58\times {10}^{-3}$, where P-QDFE and unbiased LMMSE yield BERs $1.05\times {10}^{-5}$ and $8.65\times {10}^{-5}$; at ${\delta }_{\mathrm{ICHOL}}={10}^{-4}$ their BERs are $8.32\times {10}^{-6}$ and $5.10\times {10}^{-5}$. Under fractional Doppler, spectral leakage raises the initial ${\overline{D}}_{\mathrm{chol}}$ to $236.3$, while the same threshold reduces it to $38.5$, corresponding to a $97.3\%$ complexity reduction. P-QDFE and unbiased LMMSE then yield $2.80\times {10}^{-4}$ and $9.11\times {10}^{-3}$, respectively. Furthermore, Table 5 and Figure 14 compare the critical path complexity for 16QAM ($L=5$, $d_{\mathrm{avg}}=25$). Iterative algorithms scale rapidly: MPA scales linearly with constellation size ($4.1\times {10}^5$ operations at $I=10$), while OAMP requires $\mathcal{O}\left(N^3\right)$ dense matrix operations ($1.4\times {10}^8$). In contrast, parallel P-QDFE decouples computational depth from the branch count, requiring only $3.7\times {10}^5$ operations, matching LMMSE ($3.4\times {10}^5$) and achieving multi-branch diversity without iterative latency scaling.

7. Conclusion

In this paper, a permutation-driven decision feedback equalization framework was proposed for AFDM, which unifies single-branch reliability and multi-branch diversity. In the single-branch P-DFE, permutation served as an interference-shaping mechanism. Guided by derived statistical bounds, a low-complexity EG algorithm was developed to approximate the optimal reliability sorting. Moreover, to decouple error propagation in high order constellations, permutation was utilized to generate structural diversity for the parallel P-MBDFE, where the optimal four-branch P-QDFE was configured via KT distance maximization. Simulations confirmed that the proposed P-QDFE suppressed error floors in high-order modulations. By utilizing LDPC coding, the non-iterative P-QDFE can attain nearly identical BER performance compared to the iterative OAMP equalizer, preventing convergence stagnation at pre-decoding thresholds. Finally, by integrating ICHOL, this parallel scheme provides an efficient receiver with promising complexity vs. BER trade-offs.