下面把你的代码先“还原成数学”，再分别从解析解、结果图的物理意义、DFT 当作“全频带滤波器”的角度，以及 ICI（载波间干扰）强弱来对比 AF 与 AFfs。

---

1) 信号与两种处理在解析层面的等价式

为便于推导，把关键对象写成二维索引（快时域 $n=0\ldots N-1$、慢时域脉冲 $m=0\ldots M-1$）。

• 发送 OFDM（每个脉冲一符号，忽略 CP）

  $$x[n,m]=\frac{1}{\sqrt N}\sum_{k=0}^{N-1} X[k,m]\,e^{j2\pi kn/N}, \quad X[k,m]\in \text{16-PSK}.$$

  代码里 $x= \text{pskmod}(\cdot)$，时域符号为 $F'X$。

• 延时 $\tau$ 与“跨整帧”多普勒 $\nu$：你的指数 $\exp\{j2\pi \nu \, \ell/(2MN)\}$（$\ell$ 是把 $[n,m]$ 展成一维后的索引）在第二个 $N$ 点（你用的 y2(N+1:end,:)）可分解成

  $$e^{j2\pi \nu m/M}\;\underbrace{e^{j2\pi \frac{\nu}{2MN}n}}_{\text{同一 OFDM 符号内的相位旋转}}.$$

  于是接收

  $$r[n,m]=x[(n-\tau)\!\!\!\!\pmod N,m]\cdot e^{j2\pi \nu m/M}\cdot e^{j2\pi \frac{\nu}{2MN}n}.$$

---

A. 你的 AF（时间—频率连续格点的“全帧”模糊函数）

代码里对每个延迟 $i$ 与“多普勒栅格” $f_j=\tfrac{j}{2MN}$ 计算

$$A(i,f_j)=\sum_{\ell=0}^{2MN-1} r[\ell]\;s^*[\ell-i]\;e^{-j2\pi f_j \ell},$$

然后只取第二段 $N$ 点并把 $M$ 个脉冲求和（temp4(:) 全加）。这就是对整个 2MN 点序列的时间延迟 + 频移的最大似然匹配滤波（非可分离的 2D 核），Doppler 分辨率是 $1/(2MN)$（非常细），所以对分数多普勒也能精确“对上”。

---

B. 你的 AFfs（快时域相关 + 慢时域 FFT 的“可分离”处理）

把 AFfs = abs(fft((F' * ((F*y2(N+1:end,:)).*conj(x))).')).^2 按步骤展开：

1. 快时域 DFT 解调：
   $Y[k,m]=\sum_{n} r[n,m]e^{-j2\pi kn/N}$.
   代入 $r[n,m]$ 得

   $$Y[k,m]=e^{j2\pi \nu m/M}\sum_{p=0}^{N-1} X[p,m]e^{-j2\pi p\tau/N}\underbrace{\sum_{n=0}^{N-1}e^{j2\pi\frac{(p-k)+\varepsilon}{N}n}}_{\displaystyle D_N(p-k+\varepsilon)},$$

   其中

   $$\varepsilon \triangleq \frac{\nu}{2M},\qquad D_N(\delta)=\frac{\sin\big(\pi (\,\delta\,)\big)}{\sin\big(\pi(\,\delta/N\,)\big)}e^{j\phi}$$

   是离散 Dirichlet 核。只要 $\varepsilon\neq 0$，$Y[k,m]$ 就不只受 $X[k,m]$ 一项影响，出现ICI（各 $p\neq k$ 的串扰）。

2. 快时域匹配（去掉已知数据）： 乘以 $X^*[k,m]$ 只消除了“本子载波”的相位，无法消去 ICI 项。

3. 快时域 IDFT 得到延迟/距离相关：
   $g[n,m]=\frac{1}{\sqrt N}\sum_k X^*[k,m]Y[k,m]e^{j2\pi kn/N}$，这是对每个脉冲的圆相关（理想无 ICI 时等价于时域匹配滤波）。

4. 慢时域 Doppler FFT：
   $G[n,q]=\sum_{m} g[n,m] e^{-j2\pi qm/M}$。若 $\nu$ 是整数，慢时域不会泄漏，峰恰在 $q=\nu$；若非整数，出现慢时域泄漏（sinc 扩散）。

总之，AFfs 是“先快时域、再慢时域”的可分离 2D 处理；当符号内存在相位斜率 $e^{j2\pi \nu n/(2MN)}$ 时，它与真实最优核不一致，从而产生 ICI 与峰值衰减。

---

2) 结果图 ifftshift(20*log10(.)) 的物理意义差异

• 你画的是 AFfs：已是功率 $|G|^2$。从物理量纲讲，功率应使用 $10\log_{10}(\cdot)$；你用 $20\log_{10}$ 会把动态范围拉大 2 倍（dB 数字大 2×），但峰/旁瓣位置不变。实际绘图建议：20*log10(AFfs+eps) 或更规范 10*log10(AFfs+eps)，避免 $-\infty$。

• ifftshift/fftshift：
  这里你对一个 $N\times M$ 的 RD 图整体 ifftshift，等于把零 Doppler/零频率搬到图中心，同时也对“距离”维做了循环位移。严格说，应只对 Doppler 轴做 fftshift（或 ifftshift），距离轴是否平移视你的“零延迟”定义。对 AF（若你也画它）同理：对 Doppler 维做中心化，对延迟轴是否中心化取决于是线性相关还是圆相关。

• 直观差别：

  • AF（全帧、非可分离匹配）——在 $(i=\tau,\; f=\nu/(2MN))$ 处会出现更接近“拇指钉”的尖峰，分数 Doppler 时峰值也高、泄漏低。
  • AFfs（可分离 + 子载波匹配）——快时域内存在 ICI，峰值降低、旁瓣/基底抬高；若 $\nu$ 非整数，还会在 Doppler 轴出现 sinc 泄漏。

---

3) 从“DFT 矩阵 = 全频带滤波器”的视角看两者差别

• 在 AF 中，你直接用时域参考波形 $s[\ell-i]e^{-j2\pi f \ell}$ 对整个 2MN 点序列做全带宽匹配，完全等价于“把信号过一个与其共轭时间反转并带有线性频移的全频带滤波器”。

• 在 AFfs 中，你把全带宽分解成 $N$ 个矩形分析滤波器（DFT 子带），在子带上做点乘（与 $X^*$），再综合（IDFT）。矩形分析窗意味着子带旁瓣高，只要符号内有频偏/多普勒（哪怕很小），能量便在子带间扩散，形成 ICI；这本质上是用可分离、矩形原型滤波器组去逼近一个本应非可分离的最优核。

结论：把 DFT 当作“全频带滤波器”是对 AF 的一种近似分解。近似成立的前提是符号内多普勒近似常相（$\nu\approx 0$)，否则会偏离最优。

---

4) ICI 强弱的定量比较（给出可用的闭式近似）

在快时域解调阶段，主子载波（$\Delta=p-k=0$）的幅度系数相对理想值 $N$ 的衰减近似为

$$\underbrace{\frac{\sin(\pi \varepsilon)}{N\,\sin(\pi \varepsilon/N)}}_{\displaystyle \approx\;\frac{\sin(\pi \varepsilon)}{\pi \varepsilon}\;(N\text{ 足够大时})}, \qquad \varepsilon=\frac{\nu}{2M}.$$

因此主瓣功率占比近似

$$\left(\frac{\sin(\pi \varepsilon)}{\pi \varepsilon}\right)^{\!2},\quad\text{ICI 功率}\approx 1-\left(\frac{\sin(\pi \varepsilon)}{\pi \varepsilon}\right)^{\!2}.$$

• 你的数值：$M=16,\;\nu=15\Rightarrow \varepsilon=\tfrac{15}{32}\approx0.46875$。
  则

  $$\left(\frac{\sin(\pi \varepsilon)}{\pi \varepsilon}\right)^2 \approx 0.457,$$

  也就是仅约 46% 的能量留在“正确”的子载波上，超过 50% 被打到邻近子载波（ICI）。这正是你用 AFfs 时峰值降低、旁瓣/底噪抬高的根源。
  而 AF 直接在时域做“频移+延迟”的全帧匹配，不依赖子载波正交性，几乎不受快时域 ICI 的影响（只受 Doppler 栅格化误差影响，你的栅格是 $1/(2MN)$ 很细，误差很小）。

• 慢时域（脉冲维）：你的 $\exp\{j2\pi \nu m/M\}$ 里 $\nu=15$ 是整数，所以 AFfs 在 Doppler FFT 上无泄漏（峰在 $q=15$）。若 $\nu$ 非整数，AFfs 在 Doppler 轴也会出现 sinc 泄漏；AF 因为 Doppler 栅格细，也能更好对齐。

综述：在这段代码设定下（$\nu=15$）、ICI：AFfs $\gg$ AF。当 $\nu$ 小、$\varepsilon\to 0$ 时，两者收敛；当 $\nu$ 大或 $M$ 小（$\varepsilon=\nu/(2M)$ 增大）时，AFfs 的 ICI 急剧变重。

---

5) 复杂度与实现差异（工程观）

• AF： 你用双层循环扫 $i=1\ldots N$、$j=1\ldots 2MN$，本质是全搜索匹配滤波（近似 ML），复杂度高、但最稳健。
• AFfs： 两个 FFT（快时域 $N$、慢时域 $M$）+ 若干点乘，复杂度低得多，但在符号内存在多普勒时是近似方法，会带来 ICI 和峰值损失。

---

6) 代码层面的几点微调建议（不影响上述结论）

1. A-F(i,:) = row_AF; 应为 AF(i,:) = row_AF;。
2. F2 = dftmtx(M)./sqrt(dftmtx(M)); 没用而且不对。若要单位化慢时域 DFT：F_M = dftmtx(M)/sqrt(M);
3. 作图建议只对 Doppler 轴 fftshift：

   matlab
   imagesc( fftshift(10*log10(AFfs + eps), 2) ); axis xy;

4. 功率谱用 10*log10 更合规范；避免对 0 取 log：加 eps。

---

一句话总结

• AF = 对整帧做“延迟+频移”的非可分离匹配滤波 → 对分数多普勒鲁棒、ICI 最低、但算力最高。
• AFfs = “快时域相关 + 慢时域 FFT”的可分离近似 → 计算快，但只在符号内多普勒很小（$\varepsilon\approx 0$）时等价；否则快时域出现显著 ICI，峰值压低、旁瓣/底噪抬高。
• ifftshift(20*log10(.)) 只是把零频搬中心并转 dB 显示；物理峰位/旁瓣来自上述两种匹配核的本质差异，而非 shift 本身。

好的，我们就聚焦 AFfs，把“DFT＝矩形子带滤波器组”的处理一般化成“可设计窄带滤波器组（channelized / FBMC 式）”，并从解析解角度给出等价形式、ICI 的来源与可控性、以及检测器（相关函数）如何定义与融合。核心思想：用频域子空间投影替代“把所有子载波在快时域一次性累加（IDFT）”。

---

1) 把 AFfs 抽象成“滤波器组+相关”的通式

令单个脉冲的时域参考波形为 $s[n]$，接收为 $r[n]$（这里我们先固定慢时域脉冲 $m$，延迟 $i$ 与多普勒归一化为 $\varepsilon=\nu/(2M)$，与代码一致）。AFfs 的本质是：

• 先对 $r$ 做分析滤波器组（DFT 等价于 $N$ 个矩形子带）得到子带系数 $Y[k]$；
• 用已知调制 $X[k]$ 做“子带匹配”（乘以 $X^*[k]$）；
• 再对 $k$ 维做 IDFT（把所有子带相干叠加到快时域延迟格点），形成“快时域相关”；
• 最后对慢时域做 FFT 得到 Doppler 维。

把它一般化：取任意原型脉冲 $g[n]$（非矩形），构造分析滤波器

$$h_k[n]=g[n]\,e^{j2\pi kn/N},\quad k=0,\ldots,N-1,$$

其频响 $H_k(\omega)=G(\omega-\tfrac{2\pi k}{N})$。对每个脉冲：

• 子带输出： $Y[k]=\sum_n r[n]\,h_k^*[n]$；
• 参考子带： $S[k]=\sum_n s[n]\,h_k^*[n]$（或直接用已知 $X[k]$ 的等效参考）；
• 子带相关（快时域匹配的频域形式）： $Z[k]=Y[k]\;S^*[k]$。

若像 AFfs 一样把所有 $k$ 相干合成（IDFT）得到

$$c[i]=\frac{1}{\sqrt N}\sum_{k=0}^{N-1} Z[k]\,e^{j2\pi ki/N},$$

则这一步就是“把各子带的相关量在快时域延迟栅格上相位对齐并求和”。

结论：AFfs 是“分析滤波器组 $\to$ 子带匹配 $\to$ 子带相干合成”的可分离处理；当 $g[n]$ 取矩形窗时，就退化为 DFT（你的代码）。

---

2) ICI 的解析表达：由原型脉冲的交叉模糊函数决定

考虑单脉冲内的分数多普勒 $\varepsilon$，忽略慢时域相位；把 $r[n]$ 写成

$$r[n]=\sum_{p=0}^{N-1} X[p]\;s_p[n-\tau]\;e^{j2\pi (p+\varepsilon)n/N},$$

其中 $s_p[n]=g[n]e^{j2\pi pn/N}$（用原型 $g$ 成形的第 $p$ 个子带“子载波”）。

子带 $k$ 的分析输出

$$Y[k]=\sum_{p} X[p]\underbrace{\sum_n g[n-\tau]\;g[n]\;e^{j2\pi (p-k+\varepsilon)n/N}}_{\displaystyle A_g\big(\tau,\tfrac{p-k+\varepsilon}{N}\big)\;e^{-j2\pi p\tau/N}},$$

其中

$$A_g(\tau,\nu)\triangleq\sum_n g[n-\tau]\;g[n]\;e^{j2\pi \nu n}$$

是原型脉冲 $g$ 的离散模糊函数。于是子带相关量

$$Z[k]=Y[k]\;S^*[k]=\sum_{p} X[p]X^*[k]\;A_g\!\left(\tau,\tfrac{p-k+\varepsilon}{N}\right)\,e^{-j2\pi p\tau/N}.$$

关键结论

• ICI 系数正是 $A_g\big(\tau,\frac{p-k+\varepsilon}{N}\big)$ 的旁瓣：当 $p\neq k$ 或 $\varepsilon\neq 0$ 时，非零。
• 矩形窗 $g=\mathbf{1}_{[0,N-1]}$ 时，$A_g$ 化为 Dirichlet 核（旁瓣慢衰、ICI 重）；
• 选择频域定形更好的 $g$（如 Kaiser/Dolph-Chebyshev/PHYDYAS 等），可显著压低 $A_g$ 的离轴能量，从而减小 ICI。
• 但 $A_g$ 在 $\tau$ 与 $\nu$ 上满足海森堡不等式：频域更“尖” $\Rightarrow$ 时间域更“宽”（可能增大 ISI/旁瓣），这是根本权衡。

---

3) “不把各子带相干相加”的检测器：多子空间能量/相关的投影检波

你提出：用可设计窄带滤波器划分多个快时域子空间，分别计算相关函数，而不是像 AFfs 那样把它们相干求和（IDFT 叠加）。理论上可这样刻画：

• 记 $\mathbf{h}_{k,i}$ 为第 $k$ 子带、延迟 $i$ 的参考向量（已含 $g$ 与调制）；
• 全带匹配滤波对应的检测统计量（对单脉冲，为简洁忽略慢时域）： $$T_{\text{MF}}(i)=\big|\langle \mathbf{r},\sum_k \mathbf{h}_{k,i}\rangle\big|^2.$$
• 子空间分解 + 非相干融合（你的设想）： $$T_{\text{NC}}(i)=\sum_k \big|\langle \mathbf{r},\mathbf{h}_{k,i}\rangle\big|^2 =\mathbf{r}^\mathrm{H}\!\!\left(\sum_k \mathbf{h}_{k,i}\mathbf{h}_{k,i}^\mathrm{H}\right)\!\mathbf{r}.$$

二者差别在于核函数：

• $T_{\text{MF}}$ 用的是 $\big(\sum_k\mathbf{h}_{k,i}\big)\big(\sum_k\mathbf{h}_{k,i}\big)^\mathrm{H}$：跨子带存在相干交叉项；
• $T_{\text{NC}}$ 用的是 $\sum_k \mathbf{h}_{k,i}\mathbf{h}_{k,i}^\mathrm{H}$：没有跨子带交叉项，因此不会把不同子带的 ICI 相干“叠噪”。

因此，你的设计在理论上是合理的：它等价于把全带模板投影到一组近正交的频域子空间，并对每个投影能量做“相关—取模平方—求和”的非相干融合。这样可抑制由分数多普勒引起的跨子带相位错配所带来的“相干错误叠加”，从而降低 ICI 的有效积累。

---

4) 何时“非相干融合”优于“相干叠加”？（SNR/ICI 权衡）

设噪声为白噪声、子带近正交。

• 当模型无失配（$\varepsilon\approx 0$ 或用正确补偿）时，$T_{\text{MF}}$ 为最优（Neyman–Pearson 意义），而 $T_{\text{NC}}$ 会有最多约 3 dB 的非相干损失（因为把本可相干合并的能量做了功率叠加）。
• 当存在分数多普勒导致各子带的相位倾斜不同、相干合成错相时，$T_{\text{MF}}$ 的峰值被系统性削弱且 ICI 相干累积，而 $T_{\text{NC}}$不引入错相损失且不相干叠加 ICI，常常净收益为正（尤其当 $|\varepsilon|$ 大、或 $g$ 不是理想窄带时）。

提示：你也可以在每个子带内先做小尺度 Doppler 细补偿或估计子带相位，再做相干融合（“分波段相位校正 + 合并”），可在低 ICI 的同时逼近最优 SNR。

---

5) ICI 的“可设计”上界：用 $g$ 的模糊函数直接给界

根据上节推导，从子带 $p$ 泄到 $k$ 的 ICI 系数为

$$C_{p\to k}(\tau,\varepsilon)=A_g\!\left(\tau,\tfrac{p-k+\varepsilon}{N}\right).$$

因此可用以下设计准则抑制 ICI：

1. 频域旁瓣：选择 $g$ 使得 $|A_g(\tau,\nu)|$ 随 $|\nu|$ 快速衰减（典型：Kaiser、Dolph-Chebyshev、最小旁瓣窗；FBMC/PHYDYAS 原型）。
2. 子带宽度：减小每个子带等效带宽（增加子带数），使得当 $|p-k|\ge 1$ 时 $|A_g|$ 很小。
3. 时间展宽与 ISI：更窄的频带必然带来更长的 $g[n]$（重叠处理），$|A_g(\tau,\cdot)|$ 在 $\tau$ 维更宽 → 需要足够的循环前缀/保护间隔或重叠正交设计（如 OQAM/PR 滤波器组）以控制 ISI。
4. 帧间/子带正交性：优选**紧框架/完全重构（PR）**滤波器组，满足 $\sum_k |H_k(\omega)|^2=\text{常数}$，保证能量守恒与良好的数值稳定性。

---

6) 把慢时域（Doppler）合进来：子带 RD 图与融合

对每个脉冲 $m$ 重复上节步骤，得到 $Z_m[k]$。做延迟合成（可相干或非相干），再对慢时域做 Doppler FFT：

• 子带 RD 图：

  $$C_k[i,q]=\sum_{m=0}^{M-1} c_k[i,m]\;e^{-j2\pi qm/M},$$

  其中 $c_k[i,m]$ 是第 $k$ 子带、延迟 $i$ 的（快时域）相关输出。

• 融合方式：

  • 非相干融合：$\;C_{\text{NC}}[i,q]=\sum_k |C_k[i,q]|^2$（鲁棒抑制 ICI）；
  • 相干融合（带相位校正）：$\;C_{\text{CC}}[i,q]=\big|\sum_k \alpha_k(q)\,C_k[i,q]\big|^2$，其中 $\alpha_k(q)$ 是子带相位补偿（可依据 $g$ 的群时延与 $\varepsilon$ 近似给定，或数据驱动估计）。

物理图像：多普勒在快时域引入“频率内斜率”，跨宽频带会造成跨子带相位不一致；把频带分块并分别成像，可以避免把“彼此错相”的量强行相干叠加。

---

7) 与 AF（全帧匹配）的关系：投影—重构的等价性

设 $\mathcal{P}_k$ 是第 $k$ 子带的**（近）正交投影**算子（理想砖墙滤波器组就是正交投影），则

$$\sum_k \langle \mathcal{P}_k r,\;\mathcal{P}_k s_{\tau,\nu}\rangle = \left\langle r,\;\Big(\sum_k \mathcal{P}_k\Big)s_{\tau,\nu}\right\rangle.$$

若 $\sum_k \mathcal{P}_k=I$（PR 滤波器组），把各子带的相关“相干求和”就完全等价于全带匹配滤波（AF）；
而不相干求和则对应于核 $\sum_k \mathcal{P}_k s\,s^\mathrm{H}\mathcal{P}_k$，天然消除了跨子带的相干交叉项，这正是你想要降低 ICI 的数学依据。

---

8) 实用设计要点（从理论走向可实现）

1. 选择原型 $g[n]$：令 $|G(\omega)|$ 旁瓣低、过渡带陡（Kaiser/多项式窗/PHYDYAS），并使滤波器组满足紧框架/PR（如 CMFB/OFDM-OQAM 家族）。
2. 子带数与带宽：子带越窄，ICI 抑制越强，但时域展宽越大；用重叠处理（OLA）或增加保护间隔来消除 ISI。
3. 融合策略：先用非相干融合做稳健探测；若需要更高峰值/更低旁瓣，再在检测峰附近做子带相位估计 + 相干合并的二阶段处理。
4. 复杂度：分析/合成均可用多相结构实现，复杂度与 DFT 级别相当；非相干融合避免了全带 IDFT 的“错相叠加”，数值更稳定。

---

9) 小结（回答你的问题）

• 合理性：把快时域分解为可设计窄带子空间、分别计算相关、再做**（非）相干融合**，在理论上完全成立；它等价于对全带匹配滤波的子空间投影检波，可显著降低分数多普勒导致的 ICI 有效累积。
• 推导关键：ICI 系数完全由**原型窗的模糊函数 $A_g(\tau,\nu)$**决定；矩形窗给出 Dirichlet（重 ICI），优良原型则把 ICI 压到很低。
• 代价与条件：频域更干净 $\Rightarrow$ 时域更长、可能有 ISI；需 PR/紧框架的滤波器组、合适的保护策略与融合策略来平衡。

如果你愿意，下一步我可以把上面的通式写成与你代码变量一一对应的可替换实现（给出 g[n] 设计、分析/合成的矩阵形式、以及 AFfs 的“子带非相干/相干融合”两种版本），并附上可直接运行的 MATLAB 片段。

太好了！下面把“可设计窄带滤波器组 + 子空间相关”写成与你现有变量一一对应、并且可以直接替换 AFfs 的 MATLAB 代码。代码包含：

1. 原型窗 $g[n]$ 的设计（Kaiser/Dolph 可选）；
2. 分析/合成的矩阵（或等价实现）形式；
3. 两种融合版本：非相干融合 (NC) 与 相干融合 (CC)；
4. 保持你的变量命名（F, x, y2, N, M, ...），仅替换你原来的 AFfs = ... 那一行附近；
5. 代码注释把每一步与公式/变量做了对齐。

---

可直接替换的 MATLAB 片段

把这段放到你循环内、y2 已经得到且你正在计算 AFfs 的地方（即原来那行 AFfs = abs(fft((F'*((F*y2(N+1:end,:)).*conj(x))).')).^2; 的位置），用下面整段替换即可。

matlab
%% ======  子带滤波器组版 AFfs：可设计 g[n] + NC/CC 两种融合  ======
% 与你现有变量一致：N: 每脉冲采样点数, M: 脉冲数, x: N×M 子载波数据(频域), 
% y2(N+1:end, :) 是第二个OFDM符号块的时域数据 (N×M)

% ---------- 0) 基本矩阵 ----------
F  = dftmtx(N) / sqrt(N);      % N点单位化 DFT 矩阵
Fi = F';                       % IDFT

% ---------- 1) 原型窗 g[n] 设计（可选：Kaiser 或 Dolph-Chebyshev） ----------
use_dolph = false;     % 选 false 用 Kaiser，true 用 Dolph
if ~use_dolph
    % Kaiser 窗：旁瓣 ~ -60 dB（beta≈8.6），你可改小/大以调旁瓣/主瓣
    beta = 8.6; 
    g = kaiser(N, beta); 
else
    % Dolph-Chebyshev 窗：目标旁瓣 ~ -60 dB（可改），主瓣更窄
    ripple_db = 60; 
    g = chebwin(N, ripple_db);
end
g = g(:);                        % 列向量
Gdiag = diag(g);                 % 时间域对角矩阵（分析窗）

% ---------- 2) 取当前需要成像的时域块 r[n,m] 与参考 s[n,m] ----------
Rseg = y2(N+1:end,:);           % N×M，与你原代码保持一致
Sseg = Fi * x;                  % 参考的时域符号：s[n,m] = IDFT{X[k,m]}

% ---------- 3) 分析滤波器组输出（子带系数） ----------
% 子带输出 Y[k,m] = DFT{ g[n] * r[n,m] }_k
Y = F * (Gdiag * Rseg);         % N×M
% 参考子带 S[k,m] = DFT{ g[n] * s[n,m] }_k
S = F * (Gdiag * Sseg);         % N×M

% 子带“匹配”量：Z[k,m] = Y[k,m] * conj(S[k,m])   —— 与你AFfs里的 X^*Y 对齐
Z = Y .* conj(S);               % N×M

% ---------- 4) 频带划分：把 N 个频点分成 B 个子带（每带 L 个bin） ----------
% 说明：单个bin做非相干会丢延迟分辨率，所以用“子带”而非“单bin”
B = 8;                          % 子带数（可调：更大→更窄带→更强 ICI 抑制，但时域更长）
L = floor(N / B);               % 每带bin数
% 为保险起见，若 N 不是 B 的整数倍，把末带吃掉余数
N_eff = L * B;
Z = Z(1:N_eff, :);
% 预先分配
AFfs_nc = zeros(N, M);          % 最终 非相干 融合 RD 图
AFfs_cc = zeros(N, M);          % 最终 相干   融合 RD 图

% ---------- 5) 各子带成像（快时域相关 + 慢时域Doppler） ----------
% 思路：
%   对每个子带 b，只在该带频点上做 IDFT（等价：对Z做零填充后 ifft），
%   得到子带的快时域相关 c_b[n,m]；再对慢时域做 FFT 成像成 C_b[n,q]。
%   - 非相干融合 (NC)：在 Doppler 域对 |C_b|^2 做子带求和；
%   - 相干融合 (CC)：先在快时域对子带 c_b 复数相加，再做 Doppler FFT，取 |·|^2。
CB_list = cell(B,1);            % 存各子带的 c_b(n,m)

for b = 1:B
    k0 = (b-1)*L + 1; 
    k1 = b*L;
    % 取第 b 个子带的 Z_b[k,m]，并零填回 N 维以便 ifft
    Zb = zeros(N, M);
    Zb(k0:k1, :) = Z(k0:k1, :);

    % 快时域相关（仅由该子带贡献）：c_b[n,m] = IDFT_k{ Z_b[k,m] }
    %   注意：等价于你的 "F' * (X^* .* Y)" 但只保留该子带频点。
    cb = Fi * Zb;                        % N×M
    CB_list{b} = cb;                     % 保存以便后续融合

    % 子带 Doppler FFT：C_b[n,q]
    Cb = fft(cb, M, 2);                  % 沿慢时域 m 方向做 M点 FFT → N×M

    % 非相干融合累加：|C_b|^2
    AFfs_nc = AFfs_nc + abs(Cb).^2;
end

% 相干融合（可选：加入每子带相位校正 alpha_b(q) 后再相加）
%   这里先给最简单的“直接相加版”；如果需要更进一步，
%   你可以按q(或估计的 ε)设计 alpha_b(q) 做相位对齐后再相加。
c_sum = zeros(N, M);
for b = 1:B
    c_sum = c_sum + CB_list{b};          % 快时域先相干求和
end
AFfs_cc = abs( fft(c_sum, M, 2) ).^2;    % 再做慢时域 FFT，取功率

% ---------- 6) 可视化（与原风格一致；注意功率建议用 10*log10） ----------
figure; 
subplot(1,2,1);
imagesc( fftshift(10*log10(AFfs_nc + eps), 2) ); axis xy tight;
title(sprintf('NC: 子带非相干融合  (B=%d, L=%d)', B, L)); xlabel('Doppler (fftshift)'); ylabel('Delay');

subplot(1,2,2);
imagesc( fftshift(10*log10(AFfs_cc + eps), 2) ); axis xy tight;
title(sprintf('CC: 子带相干融合  (B=%d, L=%d)', B, L)); xlabel('Doppler (fftshift)'); ylabel('Delay');

% 若你只想要一个 AFfs 输出，可将 AFfs_nc 或 AFfs_cc 作为“新AFfs”
AFfs = AFfs_nc;   % 或 AFfs = AFfs_cc;

---

变量与公式的一一对齐（便于你审阅/二次推导）

• 原代码：
  AFfs = abs( fft( (F' * ((F*y2(N+1:end,:)).*conj(x))).' ) ).^2;
  等价于：对每个脉冲 $m$，先做 $Y=F\,r$，乘上 $X^*$，再做快时域 IDFT（圆相关），最后对慢时域做 FFT、取功率。矩形窗 → ICI 重。

• 现在：

  1. $g$（g/Gdiag）：分析窗（可设计），把 DFT 变成“可设计滤波器组”；
  2. $Y = F\,G\,R$，$S = F\,G\,S_{\text{time}}$（Y,S）；
  3. $Z = Y \odot S^*$（Z）；
  4. 频带分组 $\{k_0\!\!:\!k_1\}$（B,L）；
  5. 子带快时域相关 $c_b = F^\mathrm{H} Z_b$（cb = Fi*Zb）；
  6. 子带 Doppler 变换 $C_b = \text{FFT}_m\{c_b\}$；
  7. NC：$\sum_b |C_b|^2$（功率求和，抑制跨带相干 ICI）；
     CC：$|\text{FFT}_m\{\sum_b c_b\}|^2$（相干求和，SNR 高但对分数多普勒更敏感）。

---

使用与调参建议

• B/L（子带数/宽度）：
  • B↑（L↓）→ 子带更窄，分数多普勒下 ICI 抑制更强；但时域展宽、潜在 ISI 风险上升（尤其无 CP/保护时）。
• 窗型 $g$：
  • Kaiser：beta 控制旁瓣；beta≈8.6 ~ −60 dB 是个稳妥起点。
  • Dolph-Chebyshev：旁瓣等波动、主瓣更窄；对 ICI 抑制更狠，但时域更长。
• 融合模式：
  • 首选 NC（稳健、抗分数多普勒错相）；
  • 峰附近再做子带相位校正 + CC 合并可进一步逼近最优（可根据估计的 $\varepsilon$ / 目标 Doppler bin，对每带乘相位 $\alpha_b(q)$ 再相加）。

---

如果你想把 子带相位校正（$\alpha_b(q)$）也写成现成代码（例如“估计每个子带的线性相位斜率，再做 CC”），我可以把上面 c_sum 那段扩展成两阶段处理：先用 NC 找峰，再在峰附近估 $\varepsilon$ 并做校正合并。

没问题！下面是在你“子带滤波器组版 AFfs”代码基础上，**加入“两阶段：峰值定位→细 Doppler(δ) 估计→子带相位校正→相干合并”**的现成 MATLAB 片段。它能给出

• $\hat\delta$：慢时域的分数多普勒（相对 FFT bin 的小数部分，范围 $[-0.5,0.5)$）；
• $\hat\varepsilon$：符号内分数多普勒斜率（与你上文定义一致，$\hat\varepsilon=\frac{q_0+\hat\delta}{2M}$）。

用法：把这段紧接在你已经计算出 CB_list、AFfs_nc 的地方（也就是我之前给你的“子带版 AFfs”代码第 5 步之后）粘上即可运行。

matlab
%% ====== 两阶段：NC找峰 → 估计分数多普勒δ → 子带相位校正(CC) ======

% ---------- A) 在 NC 图上找全局峰 (i0, q0) ----------
[~, idx_max] = max(AFfs_nc(:));
[i0, q0] = ind2sub(size(AFfs_nc), idx_max);   % 延迟索引 i0 ∈ [1..N]，Doppler索引 q0 ∈ [1..M]

% ---------- B) 细 Doppler 估计：在 q0 附近用小栅格搜索 δ ∈ [-0.5, 0.5) ----------
Nd   = 41;                               % 栅格点数（奇数便于包含0）
delta_grid = linspace(-0.5, 0.5 - 1/Nd, Nd);  % 分数多普勒相对bin的偏移
Edelta = zeros(1, Nd);                   % 各δ的能量度量
m = 0:M-1;                               % 慢时域索引

for id = 1:Nd
    delta = delta_grid(id);
    phase_rot = exp(-1j*2*pi*delta*m/M);    % 对每个脉冲做细Doppler旋转
    Esum = 0;
    for b = 1:B
        cb = CB_list{b};                  % N×M：子带快时域相关
        cb_rot = cb .* phase_rot;         % 细补偿
        Cb_rot = fft(cb_rot, M, 2);       % Doppler FFT
        Esum = Esum + abs(Cb_rot(i0, q0)).^2;   % 只看峰值处的能量（也可用邻域能量更稳健）
    end
    Edelta(id) = Esum;
end

[~, id_best] = max(Edelta);
delta_hat = delta_grid(id_best);               % 估计的分数多普勒（慢时域，相对bin）
% 可选：也可对 delta 做抛物线插值细化，这里先用网格最大值

% ---------- C) 用 δ̂ 重新补偿每个子带，并在峰处估计“子带相位” ----------
phase_rot_hat = exp(-1j*2*pi*delta_hat*m/M);
Cb_hat_all = cell(B,1);
for b = 1:B
    cb = CB_list{b};
    cb_rot = cb .* phase_rot_hat;             % 应用细Doppler补偿
    Cb_hat_all{b} = fft(cb_rot, M, 2);        % N×M
end

% 在 (i0,q0) 处用一个小的延迟邻域来稳健估计每个子带的相位
winR = 1;                                     % 邻域半径：i ∈ [i0-winR .. i0+winR]
i_lo = max(1,   i0 - winR);
i_hi = min(N,   i0 + winR);

alpha = zeros(B,1);                           % 子带相位校正系数 α_b(q0)
for b = 1:B
    Cb = Cb_hat_all{b};
    % 邻域内复数求和作为相位参考（更抗噪）
    ref_phasor = sum(Cb(i_lo:i_hi, q0), 1);
    alpha(b) = conj(ref_phasor) ./ (abs(ref_phasor) + eps);  % 单位模对齐
end

% ---------- D) 子带相干合并（带相位校正）：得到 CC(phase-corrected) ----------
% 方案1（在 Doppler 域相干相加）：C_sum_corr[n,q] = Σ_b α_b * C_b_hat[n,q]
C_sum_corr = zeros(N, M);
for b = 1:B
    C_sum_corr = C_sum_corr + alpha(b) .* Cb_hat_all{b};
end
AFfs_cc_pc = abs(C_sum_corr).^2;               % 相干合并后的功率图

% 方案2（在快时域先合并再做Doppler）：等价近似，通常效果相当
% c_sum_corr = zeros(N, M);
% for b = 1:B
%     cb = CB_list{b} .* phase_rot_hat;       % 先做细Doppler补偿
%     % 用峰处相位对齐：等效于给该子带快时域乘常数相位
%     Cb = Cb_hat_all{b};
%     ref_phasor = sum(Cb(i_lo:i_hi, q0), 1);
%     alpha_b = conj(ref_phasor) ./ (abs(ref_phasor) + eps);
%     c_sum_corr = c_sum_corr + alpha_b .* cb;
% end
% AFfs_cc_pc = abs( fft(c_sum_corr, M, 2) ).^2;

% ---------- E) 输出/可视化 ----------
% 建议：同时对比 NC、原始 CC（未校正）、以及相位校正后的 CC
% 这里假设你已有 AFfs_cc (未校正) 与 AFfs_nc：
figure;
subplot(1,3,1);
imagesc( fftshift(10*log10(AFfs_nc + eps), 2) ); axis xy tight;
title('NC（非相干融合）'); xlabel('Doppler (fftshift)'); ylabel('Delay');

subplot(1,3,2);
imagesc( fftshift(10*log10(AFfs_cc + eps), 2) ); axis xy tight;
title('CC（未相位校正）'); xlabel('Doppler (fftshift)'); ylabel('Delay');

subplot(1,3,3);
imagesc( fftshift(10*log10(AFfs_cc_pc + eps), 2) ); axis xy tight;
title(sprintf('CC（相位校正）  \\deltâ=%.3f', delta_hat));
xlabel('Doppler (fftshift)'); ylabel('Delay');

% 若需要作为“最终 AFfs”参与后续PSR等计算：
AFfs = AFfs_cc_pc;

% 也给出 ε̂（符号内分数多普勒）：  ε̂ = (q0_bin_index + δ̂) / (2M)
% 注意 q0 是未shift的 FFT 索引（从1计数），对应离散频率 (q0-1)/M
eps_hat = ((q0 - 1) + delta_hat) / (2*M);
fprintf('Estimated fractional Doppler: delta_hat = %.4f (of a bin), epsilon_hat = %.6f\n', delta_hat, eps_hat);

说明与要点

• δ̂ 搜索：对每个候选 δ，对 所有子带 的 cb 统一做慢时域相位旋转，再做 Doppler FFT，并在全局峰 $(i_0,q_0)$ 处累加能量，选能量最大的 δ。这样既抑制了泄漏，又提高稳健性。
• 子带相位校正 $\alpha_b$：用峰值延迟附近的小邻域复数和估计每个子带的相位，再把各子带在 Doppler 域相干相加。这一步消除了“跨子带错相”，避免把 ICI 错相叠加成噪底。
• $\hat\varepsilon$：与你之前定义一致，可直接用于分析“符号内相位斜率”的影响：$\hat\varepsilon=\frac{q_0+\hat\delta}{2M}$。

如果你想把 δ̂ 的搜索换成抛物线插值（三点幅度对数抛物线）或改用二维邻域能量，也很好加；要做“自适应子带数量/窗参数”的网格搜索，也可以在外层 loop 上包起来，我可以再给你一版自动调参脚本。