请深入学习一下下面的推导过程，用中文文本的方式讲给我，并用符合IEEE TSP风格的推导过程用latex代码重新整理一下：
\documentclass[11pt]{article}

\usepackage[a4paper,margin=1in]{geometry}
\usepackage{amsmath,amssymb,amsthm,bm}
\usepackage{mathtools}
\usepackage{hyperref}

% ---------- Macros ----------
\newcommand{\C}{\mathbb{C}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\E}{\mathbb{E}}
\newcommand{\ii}{\mathrm{i}}

\newcommand{\ww}{w}
\newcommand{\Cw}{\mathbf C_{\ww}}
\newcommand{\Hk}{\mathbf H_k}
\newcommand{\Mv}{\mathbf M_v}
\DeclareMathOperator{\Ree}{Re}
\DeclareMathOperator{\Imm}{Im}
\newcommand{\indic}[1]{\mathbb{1}{#1}}

\newtheorem{assumption}{Assumption}

\title{Fourth-Moment Analysis of CE-OFDM Signals in a Single-User ISAC System}
\author{}
\date{}

\begin{document}
\maketitle

\section{System Model and Assumptions}

Consider the fourth-moment analysis of Constant Envelope (CE)-OFDM signals in a single-user Integrated Sensing and Communication (ISAC) system.
The core objective is to evaluate the fourth moment (equivalently expressed later as $\E\{\|F x\|_4^4\}$).

\subsection{Transmitted Symbols and Constellation Characteristics}

Let the transmitted symbol vector $\bm{s}\in\C^{N_d\times 1}$ be independent and identically distributed (i.i.d.), and its constellation $\mathcal{S}$ satisfies:

\begin{assumption}[Unit power normalization]
Constellation points are normalized to unit power:
\begin{equation}
\E(|s|^2)=1,\quad \forall s\in\mathcal{S}.
\end{equation}
\end{assumption}

\begin{assumption}[Zero mean and zero pseudo-variance]
The mean and pseudo-variance of constellation points are both zero:
\begin{equation}
\E(s)=0,\qquad \E(s^2)=0,\quad \forall s\in\mathcal{S}.
\end{equation}
\end{assumption}

Define the (normalized) kurtosis of the constellation as
\begin{equation}
\mu_4 \triangleq
\frac{\E{|s-\E(s)|^4}}{\big(\E{|s-\E(s)|^2}\big)^2}
= \E(|s|^4),
\end{equation}
where the last equality uses Assumptions 1--2 and $\E(|s|^2)=1$.

\subsection{Frequency-/Time-Domain Signal Modeling}

\paragraph{1) Hermitian-symmetric frequency-domain vector.}
To obtain a real-valued time-domain OFDM symbol, map $\bm{s}$ to a Hermitian-symmetric frequency-domain vector $\bm{v}\in\C^{N\times 1}$:
\begin{equation}
\bm{v}

\big[,0,\ s_1,\ s_2,\ \ldots,\ s_{N_d},\ \bm{z}p^{\mathsf{T}},\ 0,\ s{N_d}^,\ \ldots,\ s_2^,\ s_1^,\big]^{\mathsf{T}},
\end{equation}
where $\bm{z}_p\in\R^{N_z\times 1}$ is the zero-padding vector, and the total length is
\begin{equation}
N = 2(N_d+1)+N_z.
\end{equation}
The Hermitian symmetry is
\begin{equation}
v[0]=v[N/2]=0,\qquad v[N/2+k]=v^[N/2-k],\ \ 1\le k\le N/2-1.
\end{equation}

\paragraph{2) Time-domain signal generation.}
Let $F_N\in\C^{N\times N}$ be the normalized Discrete Fourier Transform (DFT) matrix. The time-domain signal is generated via inverse DFT:
\begin{equation}
\tilde{\bm{v}} = F_N^{\mathsf{H}}\bm{v}\in\R^{N\times 1}.
\end{equation}
Using Hermitian symmetry, the $n$-th element can be written as
\begin{equation}
\tilde{v}[n]
= \frac{2}{N}\sum_{k=1}^{N/2-1}
\Big(
\Ree{s[k]}\cos!\big(\tfrac{2\pi k n}{N}\big)
-\Imm{s[k]}\sin!\big(\tfrac{2\pi k n}{N}\big)
\Big).
\end{equation}

\paragraph{3) CE-OFDM signal definition.}
Define the CE-OFDM signal (elementwise exponential) as

Let $\mathbf{x} \in \mathbb C^{N\times 1}$ be defined as
\begin{equation}
\mathbf{x} = \tilde{A}\exp!\left(j \tilde{h}\mathbf{f}n^{H}\bm v \right)\in\C^{N\times 1},
\end{equation}
where the exponential is applied element-wise.
The $k$-th frequency-domain autocorrelation sample is defined as
\begin{align}
r_k&=\mathbf{x}^H\mathbf{J}k\mathbf{x}\notag\
& = \sqrt{N}\mathbf{x}^H\mathbf{F}N^H\text{Diag}(\mathbf{f}{N-k+1})\mathbf{F}N\mathbf{x}\notag\
&=\sum{n=1}^N\big|\mathbf{f}n\mathbf{x}\big|^2\exp!\left(j\frac{2\pi k (n-1)}{N}\right)\notag\
&=\sum{n=1}^N\left(\sum{p=1}^N\sum{q=1}^N\mathbf{f}_n[p]\mathbf{f}n^*[q]\exp\left(j\left(\mathbf{f}p^H\mathbf{v}-\mathbf{f}q^H\mathbf{v}\right)\right)\right)\exp!\left(j\frac{2\pi k (n-1)}{N}\right)\notag\
& = \frac{1}{N} \sum{n=1}^N\exp!\left(j\frac{2\pi k (n-1)}{N}\right)\sum{p=1}^N\sum{q=1}^N\exp\bigg[j(\mathbf{f}_p^H-\mathbf{f}_q^H)\mathbf{v})\bigg]\exp\left(-j\frac{2\pi(n-1)(p-q)}{N}\right)
\end{align}

where $\mathbf{f}_n^{H}$ denotes the $n$-th row of the unitary DFT matrix $\mathbf{f}_n$.

\begin{align}
|r_k|^2 &= \sum_{n=1}^N\big|\mathbf{f}n\mathbf{x}\big|^2\exp!\left(j\frac{2\pi k (n-1)}{N}\right)\sum{m=1}^N\big|\mathbf{f}m\mathbf{x}\big|^2\exp!\left(-j\frac{2\pi k (m-1)}{N}\right)\notag\
& = \sum{n=1}^N\sum_{m=1}^N \big|\mathbf{f}_n\mathbf{x}\big|^2\big|\mathbf{f}_m\mathbf{x}\big|^2\exp!\left(j\frac{2\pi k (n-m)}{N}\right)
\end{align}

where
\begin{align}
&\big|\mathbf{f}n\mathbf{x}\big|^2\big|\mathbf{f}m\mathbf{x}\big|^2 \notag\
&=\frac{1}{N^2}\sum{p=1}^N\sum{q=1}^N\exp\bigg[j(\mathbf{f}p^H-\mathbf{f}q^H)\mathbf{v})\bigg]\exp\left(-j\frac{2\pi(n-1)(p-q)}{N}\right)\sum{p'=1}^N\sum{q'=1}^N\exp\bigg[-j(\mathbf{f}{p'}^H-\mathbf{f}{q'}^H)\mathbf{v})\bigg]\exp\left(j\frac{2\pi(m-1)(p'-q')}{N}\right)\notag\
&= \frac{1}{N^2}\sum_{p,q,p'q'}\exp(j\mathbf{a}^H\mathbf{v})\exp\left(-j\frac{2\pi(n-1)(p-q)}{N}\right)\exp\left(j\frac{2\pi(m-1)(p'-q')}{N}\right)
\end{align}
with
\begin{equation}
\mathbf{a}^H = \mathbf{f}p^H-\mathbf{f}q^H-\mathbf{f}{p'}^H+\mathbf{f}{q'}^H
\end{equation}

\begin{align}
|r_k|^2 & = \sum_{n=1}^N\sum_{m=1}^N \big|\mathbf{f}n\mathbf{x}\big|^2\big|\mathbf{f}m\mathbf{x}\big|^2\exp!\left(j\frac{2\pi k (n-m)}{N}\right)\notag\
& =\frac{1}{N^2}\sum{n=1}^N\sum{m=1}^N\exp!\left(j\frac{2\pi k (n-m)}{N}\right)\sum_{p,q,p'q'}\exp\left(-j\frac{2\pi(n-1)(p-q)}{N}\right)\exp\left(j\frac{2\pi(m-1)(p'-q')}{N}\right)\exp(j\mathbf{a}^H\mathbf{v})
\end{align}

\begin{align}
\mathbb{E}{|r_k|^2}& =\frac{1}{N^2}\mathbb{E}\left{\sum_{n=1}^N\sum_{m=1}^N\exp!\left(j\frac{2\pi k (n-m)}{N}\right)\sum_{p,q,p'q'}\exp\left(-j\frac{2\pi(n-1)(p-q)}{N}\right)\exp\left(j\frac{2\pi(m-1)(p'-q')}{N}\right)\exp(j\mathbf{a}^H\mathbf{v})\right}\notag\
& =\frac{1}{N^2}\sum_{n=1}^N\sum_{m=1}^N\exp!\left(j\frac{2\pi k (n-m)}{N}\right)\sum_{p,q,p'q'}\exp\left(-j\frac{2\pi(n-1)(p-q)}{N}\right)\exp\left(j\frac{2\pi(m-1)(p'-q')}{N}\right)\notag\
& \times \mathbb{E}\left{\exp(j\mathbf{a}^H\mathbf{v})\right}\notag\
\end{align}

\subsection{Characteristic function approximation (Lindeberg-CLT)}

We now replace the previous definition of $\bm a$ by
\begin{equation}
\bm a = \tilde h,(\bm f_p-\bm f_q-\bm f_{p'}+\bm f_{q'}),\qquad
\bm a^{\mathsf H}=\tilde h,(\bm f_p^{\mathsf H}-\bm f_q^{\mathsf H}-\bm f_{p'}^{\mathsf H}+\bm f_{q'}^{\mathsf H}).
\end{equation}
(If the new model does not include $\tilde h$, simply set $\tilde h=1$ in what follows.)

Using the Hermitian symmetry of $\bm v$ and the conjugate property of DFT columns, one still has
$a[N-r]=a^*[r]$. Hence,
\begin{equation}
\bm a^{\mathsf H}\bm v=\sum_{r=1}^{N_d} z_r,
\qquad
z_r \triangleq a[r]^* s_r + a[r] s_r^*
=2,\Re{a[r]^* s_r}.
\end{equation}
The random variables $\{z_r\}$ are independent but generally not identically distributed.

Under the Lindeberg condition
\begin{equation}
\max_{r}\frac{|a[r]|^2}{\sum_{r=1}^{N_d}|a[r]|^2}\to 0,
\qquad \mu_4<\infty,
\end{equation}
a cumulant expansion truncated to fourth order gives the approximation
\begin{equation}\label{eq:charfun_approx_new}
\E!\left[\exp!\big(\ii,\bm a^{\mathsf H}\bm v\big)\right]
\approx
\exp!\left(
-\sum_{r=1}^{N_d}|a[r]|^2
+\frac{\mu_4-2}{4}\sum_{r=1}^{N_d}|a[r]|^4
\right).
\end{equation}

After summing over $k$ (DFT orthogonality), the new sign pattern yields the modular constraint
\begin{equation}
p-q-p'+q' \equiv 0\ (\mathrm{mod}\ N)
\qquad\Leftrightarrow\qquad
p+q'\equiv q+p'\ (\mathrm{mod}\ N).
\end{equation}
A convenient re-parameterization is
\begin{equation}
(p,q,p',q')=(u+v,\ v,\ u,\ 0)\quad(\mathrm{mod}\ N).
\end{equation}

Let $\omega\triangleq 2\pi/N$. Define
\begin{align}
U_r(u,v)
&\triangleq e^{j\omega r(u+v)}-e^{j\omega rv}-e^{j\omega ru}+1 \nonumber\
&= (1-e^{j\omega ru})(1-e^{j\omega rv}).
\end{align}
With the normalized DFT convention, this implies
\begin{equation}
a[r]=\frac{\tilde h}{N}U_r(u,v).
\end{equation}
Introduce the deterministic sums
\begin{equation}
S_2(u,v)\triangleq \sum_{r=1}^{N_d}|U_r(u,v)|^2,\qquad
S_4(u,v)\triangleq \sum_{r=1}^{N_d}|U_r(u,v)|^4.
\end{equation}
Then
\begin{equation}
\sum_{r=1}^{N_d}|a[r]|^2=\frac{\tilde h^2}{N^2}S_2(u,v),
\qquad
\sum_{r=1}^{N_d}|a[r]|^4=\frac{\tilde h^4}{N^4}S_4(u,v),
\end{equation}
and \eqref{eq:charfun_approx_new} becomes
\begin{equation}
\E!\left[\exp!\big(j\bm a^{\mathsf H}\bm v\big)\right]
\approx
\exp!\left(
-\frac{\tilde h^2}{N^2}S_2(u,v)
+\frac{\mu_4-2}{4}\frac{\tilde h^4}{N^4}S_4(u,v)
\right).
\end{equation}

\section{Closed-Form Expressions for $S_2$ and $S_4$}

Define the auxiliary function
\begin{equation}
\Psi(d)\triangleq \sum_{r=1}^{N_d}\cos!\Big(\frac{2\pi}{N} r d\Big),
\qquad \Psi(d)=\Psi(-d).
\end{equation}
Then the standard closed form is
\begin{equation}
\boxed{
\Psi(d)=
\begin{cases}
N_d, & d\equiv 0\ (\mathrm{mod}\ N),\
\dfrac{\sin!\left(\frac{N_d\omega d}{2}\right)\cos!\left(\frac{(N_d+1)\omega d}{2}\right)}
{\sin!\left(\frac{\omega d}{2}\right)}, & \text{otherwise}.
\end{cases}}
\end{equation}

\subsection{Second-order sum (new index/sign pattern)}
For $U_r(p,q,p',q') = e^{j\omega rp}-e^{j\omega rq}-e^{j\omega rp'}+e^{j\omega rq'}$
(one may equivalently write $U_r=e^{j\omega rp}+e^{j\omega rq'}-e^{j\omega rq}-e^{j\omega rp'}$),
the closed form becomes
\begin{equation}
\boxed{
S_2(p,q,p',q')

(4N_d+2)\Big[
\Psi(p-q')+\Psi(q-p')
-\Psi(p-q)-\Psi(p-p')-\Psi(q'-q)-\Psi(q'-p')
\Big].}
\end{equation}

\subsection{Fourth-order sum (updated difference vector)}
Let $(c_0,c_1,c_2)=\big(\tfrac{3}{2},-2,\tfrac{1}{2}\big)$ and define the difference vector
\begin{equation}
\bm g=\big[g_1,\ldots,g_6\big]^{\mathsf{T}}

\big[
p-q',\ q-p',\ p-q,\ p-p',\ q'-q,\ q'-p'
\big]^{\mathsf{T}}.
\end{equation}
A compact expression consistent with the provided formula is
\begin{align}
S_4(p,q,p',q')
&=
(16N_d+16)\Big[
\Psi(p-q')+\Psi(q-p')
-\Psi(p-q)-\Psi(p-p')-\Psi(q'-q)-\Psi(q'-p')
\Big]\nonumber\
&\quad
+2\sum_{m=0}^{2}\sum_{n=0}^{2} c_m c_n
\sum_{i=1}^{6}\sum_{j'=1}^{6}
\Big[
\Psi(m g_i - n g_{j'}) + \Psi(m g_i + n g_{j'})
\Big].
\end{align}

\subsection{Kernel-matrix form (updated $\mathbf b$ and $\mathbf n(u,v)$)}
Define

$$\mathbf b\triangleq[1,-1,-1,1]^T,\qquad \mathbf n(u,v)\triangleq[u+v,\ v,\ u,\ 0]^T.$$

Define the $4\times 4$ kernel matrix
\begin{equation}
[\mathbf K_2(u,v)]{ij}\triangleq \Psi(n_i-n_j),
\qquad i,j\in{1,2,3,4},
\end{equation}
then
\begin{equation}
S_2(u,v)=\mathbf b^T\mathbf K_2(u,v)\mathbf b.
\end{equation}
Let $\mathbf w\triangleq \mathbf b\otimes \mathbf b\in\mathbb R^{16}$ and index pairs $(i,j)$ as one index in $\{1,\dots,16\}$.
Define the $16\times16$ kernel
\begin{equation}
[\mathbf K_4(u,v)]{(i,j),(k,\ell)}\triangleq
\Psi\big((n_i-n_j)+(n_k-n_\ell)\big),
\end{equation}
then
\begin{equation}
S_4(u,v)=\mathbf w^T\mathbf K_4(u,v)\mathbf w.
\end{equation}
Finally,
\begin{equation}
\E!\left[\exp!\big(j\bm a^{\mathsf{H}}\bm v\big)\right]
\approx
\exp!\left(
-\frac{\tilde{h}^2}{N^2}\mathbf b^T\mathbf K_2(u,v)\mathbf b
+\frac{\mu_4-2}{4}\frac{\tilde{h}^4}{N^4}(\mathbf b\otimes \mathbf b)^T\mathbf K_4(u,v)(\mathbf b\otimes \mathbf b)
\right).
\end{equation}

% 2nd-order truncation (discard S4):
\begin{equation}
\mathbb E!\left[\exp(j\bm a^H\bm v)\right]\approx
\exp!\left(-\frac{\tilde h^2}{N^2}S_2(p,q,p',q')\right).
\end{equation}

\begin{equation}
\boxed{
\mathbb E{|r_k|^2}\approx
\frac{1}{N^3}\sum_{n=1}^{N}\sum_{m=1}^{N} g_n g_m,e^{j\frac{2\pi}{N}k(n-m)}
\sum_{p,q,p',q'}
e^{-j\frac{2\pi}{N}\left((n-1)(q-p)+(m-1)(p'-q')\right)}
\exp!\left(-\frac{\tilde h^2}{N^2}S_2(p,q,p',q')\right).
}
\end{equation}

% ===== Iceberg extraction: terms with a = 0 (pairings) =====
% For a^H = f_p^H - f_q^H - f_{p'}^H + f_{q'}^H:
% Pairing-1: p=q and q'=p'
% Pairing-2: p=p' and q=q'
% These two pairings contribute (2N^2 - N) quadruples in (p,q,p',q').

\begin{equation}
\mathcal I_k(\bm g)

\frac{1}{N^3}(2N^2-N)\sum_{n=1}^{N}\sum_{m=1}^{N} g_n g_m e^{j\frac{2\pi}{N}k(n-m)}

\left(2-\frac{1}{N}\right)\big|\bm f_k^H\bm g\big|^2.
\end{equation}

\begin{equation}
\boxed{
\mathbb E{|r_k|^2}\approx \mathcal I_k(\bm g)+\mathcal S_k(\bm g),
\qquad
\mathcal I_k(\bm g)=\left(2-\frac{1}{N}\right)|\bm f_k^H\bm g|^2,
}
\end{equation}

\begin{equation}
\mathcal S_k(\bm g)\triangleq
\frac{1}{N^3}\sum_{n=1}^{N}\sum_{m=1}^{N} g_n g_m,e^{j\frac{2\pi}{N}k(n-m)}
\sum_{p,q,p',q'\notin \text{pairings}}
e^{-j\frac{2\pi}{N}\left((n-1)(q-p)+(m-1)(p'-q')\right)}
\exp!\left(-\frac{\tilde h^2}{N^2}S_2(p,q,p',q')\right).
\end{equation}

% ============================================================
% TSP-style derivation starting from the 4-index summation
% (no (u,v,t) re-parameterization), and reusing K_2(·) in 4 args
% ============================================================

\paragraph{Starting point (4-index form).}
Let $\omega \triangleq 2\pi/N$ and $\alpha \triangleq \tilde h^2/N^2$. From
\begin{equation}\label{eq:Ek_rk2_start}
\boxed{
\mathbb E{|r_k|^2}\approx
\frac{1}{N^3}\sum_{n=1}^{N}\sum_{m=1}^{N} g_n g_m,e^{j\omega k(n-m)}
\sum_{p,q,p',q'}
e^{-j\omega\left((n-1)(q-p)+(m-1)(p'-q')\right)}
\exp!\left(-\alpha,S_2(p,q,p',q')\right),
}
\end{equation}
we re-introduce the \emph{four-parameter} kernel representation of $S_2$.
Define the index vector and kernel matrix

$$\mathbf n(p,q,p',q')\triangleq [p,\ q,\ p',\ q']^T\ (\mathrm{mod}\ N),\qquad [\mathbf K_2(p,q,p',q')]_{ij}\triangleq \Psi(n_i-n_j),\ i,j\in\{1,2,3,4\},$$

with $\Psi(d)\triangleq \sum_{r=1}^{N_d}\cos(\omega r d)$, and set

$$\mathbf b\triangleq[1,-1,-1,1]^T.$$

Then
\begin{equation}\label{eq:S2_K2_4args}
\boxed{
S_2(p,q,p',q')=\mathbf b^T\mathbf K_2(p,q,p',q')\mathbf b,
}
\end{equation}
so that the weight in \eqref{eq:Ek_rk2_start} becomes $\exp\!\big(-\alpha\,\mathbf b^T\mathbf K_2(p,q,p',q')\mathbf b\big)$ without introducing any new 1D kernel.

\paragraph{Pairings and iceberg term $\mathcal I_k(\mathbf g)$.}
Under the CE-OFDM model $x_j=\exp(-j\tilde h\,\mathbf f_j^H\mathbf v)$, the 4th-order moment entering the derivation has the characteristic-function form
$\mathbb E[x_qx_{p'}x_p^*x_{q'}^*]=\mathbb E[\exp(j\mathbf a^H\mathbf v)]$ with

$$\mathbf a^H=\tilde h\big(\mathbf f_p^H-\mathbf f_q^H-\mathbf f_{p'}^H+\mathbf f_{q'}^H\big).$$

Hence $\mathbf a=\mathbf 0$ (equivalently $S_2=0$ and $\exp(-\alpha S_2)=1$) occurs if and only if one of the following \emph{pairings} holds:

$$\text{(A)}\ p=q\ \text{and}\ p'=q',\qquad \text{(B)}\ p=p'\ \text{and}\ q=q'.$$

These contribute $|\text{A}|+|\text{B}|-|\text{A}\cap\text{B}|=N^2+N^2-N=2N^2-N$ quadruples. Extracting them from \eqref{eq:Ek_rk2_start} yields
\begin{equation}\label{eq:Ik_def}
\mathcal I_k(\mathbf g)

\frac{1}{N^3}(2N^2-N)\sum_{n=1}^{N}\sum_{m=1}^{N} g_n g_m e^{j\omega k(n-m)}.
\end{equation}
Using

$$\sum_{n=1}^{N}\sum_{m=1}^{N} g_n g_m e^{j\omega k(n-m)} =\left|\sum_{n=1}^N g_n e^{j\omega k(n-1)}\right|^2 = N\,|\mathbf f_k^H\mathbf g|^2,$$

we obtain the closed form
\begin{equation}\label{eq:Ik_closed}
\boxed{
\mathcal I_k(\mathbf g)=\left(2-\frac{1}{N}\right)|\mathbf f_k^H\mathbf g|^2.
}
\end{equation}

\paragraph{Sea-level term $\mathcal S_k(\mathbf g)$ and elimination of the $(n,m)$ sums.}
Let the pairing set be

$$\mathcal P \triangleq \{(p,q,p',q'):\ (p=q\land p'=q')\ \text{or}\ (p=p'\land q=q')\}.$$

Define

\mathbb E|r_k|^2 \approx \mathcal I_k(\mathbf g)+\mathcal S_k(\mathbf g), \qquad \mathcal S_k(\mathbf g)\ \text{is \eqref{eq:Ek_rk2_start} restricted to}\ (p,q,p',q')\notin\mathcal P.

Next, combine the exponential factors in \eqref{eq:Ek_rk2_start} to separate the $n$- and $m$-dependent parts:

$$e^{j\omega k(n-m)}e^{-j\omega(n-1)(q-p)}=e^{j\omega(n-1)(k+p-q)},\qquad e^{-j\omega k(n-m)}e^{-j\omega(m-1)(p'-q')}=e^{j\omega(m-1)(-k-p'+q')}.$$

Therefore,

$$\sum_{n=1}^N g_n e^{j\omega(n-1)(k+p-q)}=\sqrt N\,\mathbf f_{\modN{k+p-q}}^H\mathbf g,\qquad \sum_{m=1}^N g_m e^{j\omega(m-1)(-k-p'+q')}=\sqrt N\,\mathbf f_{\modN{-k-p'+q'}}^H\mathbf g,$$

where $\modN{\cdot}$ denotes reduction modulo $N$ into $\{0,\ldots,N-1\}$. Substituting back gives the \emph{pure 4-index} sea-level representation
\begin{equation}\label{eq:Sk_4idx_DFT}
\boxed{
\mathcal S_k(\mathbf g)=
\frac{1}{N^2}!!\sum_{(p,q,p',q')\notin\mathcal P}!!
\big(\mathbf f_{\modN{k+p-q}}^H\mathbf g\big),
\big(\mathbf f_{\modN{-k-p'+q'}}^H\mathbf g\big),
\exp!\left(-\alpha,\mathbf b^T\mathbf K_2(p,q,p',q')\mathbf b\right).
}
\end{equation}

\begin{equation}
\boxed{
\mathcal Q_k=
\frac{1}{N^2}!!\sum_{(p,q,p',q')\notin\mathcal P}!!
\big(\mathbf f_{\modN{k+p-q}}^H\big),
\big(\mathbf f_{\modN{-k-p'+q'}}^H\big),
\exp!\left(-\alpha,\mathbf b^T\mathbf K_2(p,q,p',q')\mathbf b\right).
}
\end{equation}

\begin{equation}
\boxed{
\mathcal Q_k=\mathcal Q_k^{SP}+\mathcal Q_k^{bulk}
}
\end{equation}

\paragraph{Structured hyperplanes (single-pair layers) and matrix forms.}
We partition the non-pairing set into structured hyperplanes and a residual bulk:

$$\text{(H1)}\ p'=q',\qquad \text{(H2)}\ p=q,\qquad \text{(H3)}\ q=q',\qquad \text{(H4)}\ p=p',\qquad \text{bulk: complement}.$$

On each hyperplane, $\mathbf b^T\mathbf K_2(p,q,p',q')\mathbf b$ reduces (by shift-invariance of $\Psi(\cdot)$) to a function of a single modular difference $d\in\ZN$. We encode this \emph{using $K_2$ itself} via the 1D weight
\begin{equation}\label{eq:kappa_def}
\kappa[d]\triangleq \exp!\left(-\alpha,\mathbf b^T\mathbf K_2(d,0,0,0)\mathbf b\right),\qquad d\in\ZN ,
\end{equation}
which is merely a shorthand (no new model) since it is generated by the original four-argument kernel.

% ===== kappa[d] in closed form and vectorization =====

% Parameters

$$\omega \triangleq \frac{2\pi}{N},\qquad \alpha\triangleq \frac{\tilde h^2}{N^2},\qquad \mathbf b \triangleq [\,1,-1,-1,1\,]^T.$$

% Psi(d)

$$\Psi(d)\triangleq \sum_{r=1}^{N_d}\cos(\omega r d),\qquad \Psi(d)=\Psi(-d).$$

% Key simplification for the special representative (d,0,0,0)

$$\boxed{ \mathbf b^T \mathbf K_2(d,0,0,0)\mathbf b =2\big(\Psi(0)-\Psi(d)\big) =2\big(N_d-\Psi(d)\big). }$$

% kappa[d] explicit (sum form)

$$\boxed{ \kappa[d]\triangleq \exp\!\left(-\alpha\,\mathbf b^T\mathbf K_2(d,0,0,0)\mathbf b\right) =\exp\!\left(-2\alpha\big(N_d-\Psi(d)\big)\right) =\exp\!\left( -\frac{2\tilde h^2}{N^2}\Big(N_d-\sum_{r=1}^{N_d}\cos(\omega r d)\Big) \right). }$$

% ===== kappa as a length-N vector =====

$$\boldsymbol\kappa \triangleq \big[\kappa[0],\kappa[1],\ldots,\kappa[N-1]\big]^T \in \mathbb R^{N}, \qquad (\boldsymbol\kappa)_d = \kappa[d],\ \ d\in\mathbb Z_N.$$

% (Optional) elementwise definition

$$\boxed{ (\boldsymbol\kappa)_d =\exp\!\left(-2\alpha\big(N_d-\Psi(d)\big)\right), \qquad d=0,1,\ldots,N-1. }$$

\emph{Low-rank layers:} the hyperplanes (H1) and (H2) lead to two (Hermitianized) rank-$\le2$ quadratic forms

$$\mathbf Q_{k,ut}^{(0)}\triangleq \frac{1}{N}\sum_{d=1}^{N-1}\kappa[d]\ \mathbf f_{-k}^{*}\mathbf f_{k+d}^{H}, \qquad \mathbf Q_{k,ut}\triangleq \tfrac12\big(\mathbf Q_{k,ut}^{(0)}+\mathbf Q_{k,ut}^{(0)H}\big),$$ $$\mathbf Q_{k,v0}^{(0)}\triangleq \frac{1}{N}\sum_{d=1}^{N-1}\kappa[d]\ \mathbf f_{k}^{*}\mathbf f_{d-k}^{H}, \qquad \mathbf Q_{k,v0}\triangleq \tfrac12\big(\mathbf Q_{k,v0}^{(0)}+\mathbf Q_{k,v0}^{(0)H}\big).$$

% ===== Normalized DFT convention =====

$$\omega \triangleq \frac{2\pi}{N},\qquad [\mathbf F]_{n,d}=\frac{1}{\sqrt N}e^{-j\omega (n-1)d}, \quad n=1,\ldots,N,\ d=0,\ldots,N-1,$$

and let $\mathbf f_d$ denote the $d$-th column of $\mathbf F$ (indices are taken mod $N$).

% ===== Pack kappa into a length-N vector and take one FFT =====

$$\tilde\kappa[0]\triangleq 0,\qquad \tilde\kappa[d]\triangleq \kappa[d],\ d=1,\ldots,N-1, \qquad \tilde{\boldsymbol\kappa}\triangleq[\tilde\kappa[0],\ldots,\tilde\kappa[N-1]]^T\in\mathbb R^N,$$ $$\mathbf s \triangleq \mathbf F\,\tilde{\boldsymbol\kappa} =\sum_{d=1}^{N-1}\kappa[d]\mathbf f_d\in\mathbb C^N.$$

% ===== Key identity: elementwise product of DFT columns =====
For normalized DFT columns, the elementwise products satisfy

$$\mathbf f_k^*\odot \mathbf f_d=\frac{1}{\sqrt N}\mathbf f_{d-k}, \qquad \mathbf f_k\odot \mathbf f_d=\frac{1}{\sqrt N}\mathbf f_{d+k}, \qquad (\,\mathrm{mod}\ N\,).$$

Hence,

$$\sum_{d=1}^{N-1}\kappa[d]\mathbf f_{d+k} =\sqrt N\,(\mathbf f_{-k}\odot \mathbf s), \qquad \sum_{d=1}^{N-1}\kappa[d]\mathbf f_{d-k} =\sqrt N\,(\mathbf f_{k}\odot \mathbf s).$$

% ===== Compact forms for Q_{k,ut}^{(0)} and Q_{k,v0}^{(0)} without D_k =====
Using $\mathbf f_{-k}^*=\mathbf f_k$ and $\mathbf f_k^*=\mathbf f_{-k}$,

$$\boxed{ \mathbf Q_{k,ut}^{(0)} =\frac{1}{N}\,\mathbf f_k\Big(\sum_{d=1}^{N-1}\kappa[d]\mathbf f_{d+k}\Big)^H =\frac{\sqrt N}{N}\,\mathbf f_k\,(\mathbf f_{-k}\odot \mathbf s)^H, }$$ $$\boxed{ \mathbf Q_{k,v0}^{(0)} =\frac{1}{N}\,\mathbf f_{-k}\Big(\sum_{d=1}^{N-1}\kappa[d]\mathbf f_{d-k}\Big)^H =\frac{\sqrt N}{N}\,\mathbf f_{-k}\,(\mathbf f_{k}\odot \mathbf s)^H, } \qquad \mathbf s=\mathbf F\tilde{\boldsymbol\kappa}.$$

% ===== Hermitian symmetrization =====

$$\mathbf Q_{k,ut}\triangleq \tfrac12\big(\mathbf Q_{k,ut}^{(0)}+\mathbf Q_{k,ut}^{(0)H}\big), \qquad \mathbf Q_{k,v0}\triangleq \tfrac12\big(\mathbf Q_{k,v0}^{(0)}+\mathbf Q_{k,v0}^{(0)H}\big).$$

% ===== Final low-rank sum (no D_k) =====

$$\boxed{ \mathbf Q_{k,ut}+\mathbf Q_{k,v0} = \frac{\sqrt N}{2N}\Big( \mathbf f_k(\mathbf f_{-k}\odot \mathbf s)^H + (\mathbf f_{-k}\odot \mathbf s)\mathbf f_k^H + \mathbf f_{-k}(\mathbf f_{k}\odot \mathbf s)^H + (\mathbf f_{k}\odot \mathbf s)\mathbf f_{-k}^H \Big), \quad \mathbf s=\mathbf F\tilde{\boldsymbol\kappa}. }$$

\emph{Diagonal layer:} the hyperplanes (H3) and (H4) collapse (via the DFT shift/modulation property) to a diagonal matrix

$$\mathbf D_d \triangleq \mathrm{diag}\!\big(e^{-j\omega d(n-1)}\big)_{n=1}^N,\qquad \mathbf Q_{\mathrm{diag}}\triangleq \frac{1}{N^2}\sum_{d=1}^{N-1}\kappa[d]\big(\mathbf D_d+\mathbf D_{-d}\big).$$

% ===== Key idea: represent the diagonal operator by its diagonal vector =====
Since $\mathbf Q_{\mathrm{diag}}$ is diagonal, define its diagonal vector

$$\mathbf q_{\mathrm{diag}}\triangleq \mathbf Q_{\mathrm{diag}}\mathbf 1 \in \mathbb R^N,$$

% ===== Compute q_diag explicitly =====
For each $n=1,\ldots,N$,
\begin{align*}
(\mathbf q_{\mathrm{diag}})n
&=\frac{1}{N^2}\sum{d=1}^{N-1}\kappa[d]\Big(e^{-j\omega d(n-1)}+e^{+j\omega d(n-1)}\Big)\
&=\frac{2}{N^2}\sum_{d=1}^{N-1}\kappa[d]\cos!\big(\omega d(n-1)\big).
\end{align*}

% ===== Pack kappa into a length-N vector and express via DFT columns =====
Define the length-$N$ vector $\tilde{\boldsymbol\kappa}\in\mathbb R^N$ by

$$\tilde\kappa[0]\triangleq 0,\qquad \tilde\kappa[d]\triangleq \kappa[d],\ d=1,\ldots,N-1, \qquad \tilde{\boldsymbol\kappa}=[\tilde\kappa[0],\ldots,\tilde\kappa[N-1]]^T.$$

Then

$$\mathbf s \triangleq \mathbf F\tilde{\boldsymbol\kappa} = \sum_{d=1}^{N-1}\kappa[d]\mathbf f_d,$$

and

$$(\mathbf q_{\mathrm{diag}})_n =\frac{1}{N^2}\Big(\sqrt N\,s_n+\sqrt N\,s_n^*\Big) =\frac{2}{N^{3/2}}\Re\{s_n\}.$$

Therefore,

$$\boxed{ \mathbf q_{\mathrm{diag}} =\frac{2}{N^{3/2}}\Re\{\mathbf F\tilde{\boldsymbol\kappa}\} =\frac{2}{N^{3/2}}\Re\left\{\sum_{d=1}^{N-1}\kappa[d]\mathbf f_d\right\}, }$$

and the diagonal operator can be applied without any $\mathrm{diag}(\cdot)$ as

$$\mathbf Q_{\mathrm{diag}}= \left(\frac{2}{N^{3/2}}\mathbf F\tilde{\boldsymbol\kappa}\right).$$

Thus the structured part is
\begin{equation}\label{eq:Sk_sp_quadratic}
\boxed{
\mathcal S_k^{\mathrm{sp}}(\mathbf g)=\mathbf g^H\mathbf Q_k^{\mathrm{sp}}\mathbf g,\qquad
\mathbf Q_k^{\mathrm{sp}}=\mathbf Q_{\mathrm{diag}}+\mathbf Q_{k,ut}+\mathbf Q_{k,v0}.
}
\end{equation}

% ===== Normalized DFT =====

$$\omega \triangleq \frac{2\pi}{N},\qquad [\mathbf F]_{n,d}=\frac{1}{\sqrt N}e^{-j\omega (n-1)d}, \quad n=1,\ldots,N,\ d=0,\ldots,N-1,$$

and let $\mathbf f_d$ be the $d$-th column of $\mathbf F$ (indices mod $N$).

% ===== kappa vector and one FFT =====

$$\tilde\kappa[0]\triangleq 0,\qquad \tilde\kappa[d]\triangleq \kappa[d],\ d=1,\ldots,N-1, \qquad \tilde{\boldsymbol\kappa}\triangleq[\tilde\kappa[0],\ldots,\tilde\kappa[N-1]]^T,$$ $$\mathbf s\triangleq \mathbf F\tilde{\boldsymbol\kappa}.$$

% ===== s-factor matrix shared by all terms =====

$$\mathbf S(\mathbf s)\triangleq \mathbf 1\,\mathbf s^H+\mathbf s\,\mathbf 1^T, \qquad \text{so } [\mathbf S(\mathbf s)]_{n,n}=s_n^*+s_n=2\Re\{s_n\}.$$

% ===== Q_diag without diag(·): diagonal picked by Hadamard with I =====

$$\boxed{ \mathbf Q_{\mathrm{diag}} =\frac{2}{N^{3/2}}\operatorname{diag}(\Re\{\mathbf s\}) =\frac{1}{N^{3/2}}\big(\mathbf I\odot \mathbf S(\mathbf s)\big). }$$

% ===== Low-rank sum expressed with the SAME S(s) =====

$$\mathbf A_k\triangleq \mathbf f_k\mathbf f_{-k}^H+\mathbf f_{-k}\mathbf f_k^H,$$ $$\boxed{ \mathbf Q_{k,ut}+\mathbf Q_{k,v0} =\frac{\sqrt N}{2N}\big(\mathbf A_k\odot \mathbf S(\mathbf s)\big). }$$

% Assume normalized DFT F and its columns f_k.
% Let z = F \tilde{kappa}. To make the Hadamard product well-defined, lift z to an NxN matrix:

$$\mathbf z \triangleq \mathbf F\tilde{\boldsymbol\kappa},\qquad \tilde{\mathbf S}\triangleq \mathbf 1\,\mathbf z^H+\mathbf z\,\mathbf 1^T .$$

Note that $\mathrm{diag}(\tilde{\mathbf S})=2\Re\{\mathbf z\}=2\Re\{\mathbf F\tilde{\boldsymbol\kappa}\}$.

Define

$$\mathbf A_k\triangleq \mathbf f_k\mathbf f_{-k}^H+\mathbf f_{-k}\mathbf f_k^H, \qquad \mathbf Q_k^{\mathrm{sp}} =\Big(\tfrac{1}{N^{3/2}}\mathbf I+\tfrac{\sqrt N}{2N}\mathbf A_k\Big)\odot \tilde{\mathbf S}.$$

If $\mathbf Q_k^{\mathrm{sp}}$ is Hermitian PSD, then

$$\boxed{ \mathcal S_k^{\mathrm{sp}}(\mathbf g) =\mathbf g^H\mathbf Q_k^{\mathrm{sp}}\mathbf g =\big\|(\mathbf Q_k^{\mathrm{sp}})^{1/2}\mathbf g\big\|_2^2 = \left\| \left( \Big(\tfrac{1}{N^{3/2}}\mathbf I+\tfrac{\sqrt N}{2N}\mathbf A_k\Big)\odot \tilde{\mathbf S} \right)^{1/2}\mathbf g \right\|_2^2 . }$$

\paragraph{Bulk term as a vec--matrix--vec product (still 4 indices).}
Define the DFT coefficient vector $\boldsymbol\gamma\triangleq \mathbf F^H\mathbf g$ so that $\gamma_\ell=\mathbf f_\ell^H\mathbf g$.
Introduce the two $N\times N$ circulant-by-difference matrices

$$[\mathbf H_k]_{p,q}\triangleq \mathbf f_{\modN{k+p-q}}^H\mathbf g=\gamma_{\modN{k+p-q}}, \qquad [\mathbf G_k]_{p',q'}\triangleq \mathbf f_{\modN{-k-p'+q'}}^H\mathbf g=\gamma_{\modN{-k+(q'-p')}}.$$

Let $\mathbf c_k[d]\triangleq \gamma_{\modN{k+d}}$ and $\mathbf c_{-k}[d]\triangleq \gamma_{\modN{-k+d}}$. Then

$$\mathbf H_k=\mathrm{circ}(\mathbf c_k),\qquad \mathbf G_k=\mathrm{circ}(\mathbf c_{-k})^{T}.$$

Next define the $N^2\times N^2$ kernel matrix and an explicit bulk mask:

$$[\mathbf K]_{(p,q),(p',q')}\triangleq \exp\!\left(-\alpha\,\mathbf b^T\mathbf K_2(p,q,p',q')\mathbf b\right),$$ $$[\mathbf P^{\mathrm{bulk}}]_{(p,q),(p',q')}\triangleq \mathbf 1\{p\neq q,\ p'\neq q',\ p\neq p',\ q\neq q'\}\qquad(\mathrm{mod}\ N),$$

which removes all pairings and the structured hyperplanes. The bulk contribution becomes
\begin{equation}\label{eq:Sk_bulk_vec}
\boxed{
\mathcal S_k^{(\mathrm{bulk})}(\mathbf g)=
\frac{1}{N^2},
\operatorname{vec}(\mathbf H_k)^{T},
\big(\mathbf K\odot \mathbf P^{\mathrm{bulk}}\big),
\operatorname{vec}(\mathbf G_k).
}
\end{equation}

\paragraph{Final TSP-style decomposition.}
Combining \eqref{eq:Ik_closed}, \eqref{eq:Sk_sp_quadratic}, and \eqref{eq:Sk_bulk_vec}, we obtain

$$\boxed{ \mathbb E|r_k|^2 \approx \underbrace{\left(2-\frac1N\right)|\mathbf f_k^H\mathbf g|^2}_{\mathcal I_k(\mathbf g)} + \underbrace{\mathbf g^H\mathbf Q_k^{\mathrm{sp}}\mathbf g}_{\mathcal S_k^{\mathrm{sp}}(\mathbf g)} + \underbrace{\frac{1}{N^2}\operatorname{vec}(\mathbf H_k)^{T}\big(\mathbf K\odot \mathbf P^{\mathrm{bulk}}\big)\operatorname{vec}(\mathbf G_k)}_{\mathcal S_k^{(\mathrm{bulk})}(\mathbf g)}. }$$

\end{document}