下面把整段 MATLAB 代码按“运算块”逐行说明，并给出尽可能紧凑的矩阵/解析形式。为便于书写，我约定：

• 快时间（每个 OFDM 符号内）长度 $N$，慢时间（脉冲/符号个数）长度 $M$。
• 两倍串接后的总长度 $L=2MN$。
• vec(·) 表示按列堆叠向量化，$\odot$ 表示Hadamard乘（逐元素乘），$\Pi_\ell$ 是长度 $L$ 的循环移位矩阵（向下移位 $\ell$ 个采样），$\mathbf{P}_\ell$ 是长度 $N$ 的循环移位矩阵。
• $\mathbf{F}\in\mathbb{C}^{N\times N}$ 为单位化 DFT：$\mathbf{F}_{k,n}=\frac{1}{\sqrt{N}}e^{-j2\pi (k-1)(n-1)/N}$，$\mathbf{F}^H=\mathbf{F}^{-1}$。
• 记选择矩阵 $\mathbf{S}\in\{0,1\}^{L\times L}$，它把每个脉冲的下半段（行索引 $N+1{:}2N$）选出来：$\mathbf{S}$ 在那些位置对角为 1，其余为 0。
• circshift(v,i) 可写为 $\Pi_i\,v$（循环移位），exp(j 2π …) 的逐样本调制可写为对角矩阵 $\mathbf{D}(\alpha)=\mathrm{diag}\{e^{j2\pi \alpha \,k}\}_{k=0}^{L-1}$。

---

1) 初始化与随机信号

matlab
times = 1;
M = 16; N = 32*times; L = 2*M*N;
Nmc = 10; rng(1545);
tau = 4; 

• $M=16$，$N=32$，$L=2MN=1024$。
• 真实延迟为整数采样 $\tau=4$。

---

2) 基带 OFDM 生成

matlab
s = randi(16,[N,M]) - 1; 
x = pskmod(s,16);       % X \in C^{N\times M}, |x_{n,m}|=1
F = dftmtx(N)./sqrt(N); % = \mathbf{F}
y = repmat(F'*x,[2,1]); % Y = [S; S],  S = F^H X

矩阵形式：

• 定义数据矩阵 $\mathbf{X}\triangleq[x_{n,m}]\in\mathbb{C}^{N\times M}$，16-PSK，单位模。
• 时域 OFDM（每列一个符号）$\mathbf{S}=\mathbf{F}^H\mathbf{X}\in\mathbb{C}^{N\times M}$。
• 复制两次并按行拼接：$\mathbf{Y}=\begin{bmatrix}\mathbf{S}\\\mathbf{S}\end{bmatrix}\in\mathbb{C}^{2N\times M}$。

---

3) 串行化、延时与多普勒调制（“延迟-多普勒”信道）

matlab
yy1 = y(:);                  % y = vec(Y) \in C^{L}
yy2 = circshift(yy1,tau);    % y_\tau = \Pi_\tau y
yy3 = yy2 .* exp(1i*2*pi*(0:L-1)/M/N/2*nu).'; 
                             % y_{\tau,\nu} = D(\tfrac{2\nu}{MN}) y_\tau

矩阵写法：

$$\underbrace{\mathbf{z}}_{\texttt{yy3}} = \underbrace{\mathbf{D}\!\left(\tfrac{2\nu}{MN}\right)}_{\displaystyle \mathrm{diag}\big(e^{j2\pi \tfrac{2\nu}{MN}k}\big)} \underbrace{\Pi_\tau}_{\text{循环时延}}\; \underbrace{\mathbf{y}}_{\texttt{yy1}= \mathrm{vec}(\mathbf{Y})}\;,\quad k=0,\dots,L-1 .$$

这里的归一化使“离散多普勒频率”$f_D=\frac{2\nu}{MN}$（每采样循环数），对应慢时间每个符号的相位增量为
$\Delta\phi_m = 2\pi\cdot\frac{2\nu}{MN}\cdot(2N)=\frac{4\pi\nu}{M}$（见第 7 节）。

---

4) 发送向量与接收矩阵重排

matlab
xx = repmat(F'*x,[2,1]);  % xx = vec(Y) = y
xx = xx(:);               % 再向量化，确认为 y
y2 = reshape(yy3,N*2,M);  % Z = reshape(z, 2N, M)
z = y2(:);                % z = vec(Z) = y_{\tau,\nu}

即 $\mathbf{xx}=\mathbf{y}=\mathrm{vec}(\mathbf{Y})$，$\mathbf{Z}$ 是把 $\mathbf{z}$ 还原成 $2N\times M$ 的矩阵。

---

5) 逐“延迟-多普勒格点”的常规模糊函数 AF（直接定义）

matlab
temp1 = z .* exp(-1i*2*pi*(0:L-1)'/L*j);      % D(-j/L) z
temp2 = temp1 .* circshift(conj(xx), i);      % 与 \Pi_i \overline{xx} 逐点相乘
temp3 = reshape(temp2, 2N, M); 
temp4 = temp3(N+1:end, :);                    % 选取每个脉冲的后 N 个采样
row_AF(j) = sum(temp4(:));                    % 对所有脉冲求和

把上面 5 行合起来，第 $i$ 个时延、$\;j$ 个多普勒的相关量为

$$\boxed{\; \mathrm{AF}(i,j) = \mathbf{xx}^H \,\Pi_{-i}\, \mathbf{S}\, \mathbf{D}\!\left(\tfrac{j}{L}\right)\, \mathbf{z} = \mathbf{xx}^H \,\Pi_{-i}\, \mathbf{S}\, \mathbf{D}\!\left(\tfrac{j}{L}\right)\, \mathbf{D}\!\left(\tfrac{2\nu}{MN}\right)\Pi_\tau\,\mathbf{xx}\; }$$

其中 $\mathbf{S}$ 是“选后半段”的选择矩阵（先选后半段再求和 = 与全 1 向量内积，等价于上式）。
利用 $\mathbf{D}(a)\mathbf{D}(b)=\mathbf{D}(a+b)$：

$$\mathrm{AF}(i,j)=\mathbf{xx}^H \,\Pi_{-i}\, \mathbf{S}\,\mathbf{D}\!\Big(\tfrac{j}{L}+\tfrac{2\nu}{MN}\Big)\Pi_\tau\,\mathbf{xx}.$$

Monte-Carlo 累加的是功率：

matlab
cAF = cAF + abs(AF).^2;   % \; c\mathrm{AF} = \sum_{mc} |\mathrm{AF}|^2

由于 16-PSK 为单位模、独立同分布，若取对 $\mathbf{xx}$ 的期望，$\mathbb{E}\,|\mathbf{xx}^H \mathbf{A}\mathbf{xx}|^2$ 可写成 $\mathrm{tr}(\mathbf{A})^2+\mathrm{tr}(\mathbf{A}\mathbf{A}^H)$ 型（忽略高阶相关），但这里我们直接给出了确定型的二次型表达式，已能揭示峰值所在。

---

6) “快-慢时间模糊函数”AFfs（等效的匹配-解调-慢时FFT管道）

matlab
AFfs = abs( fft( (F' * ((F*y2(N+1:end,:)).*conj(x))).' ) ).^2;

按维度拆解（$\mathbf{R}\triangleq \mathbf{Z}(N+1{:}2N,:)\in\mathbb{C}^{N\times M}$ 是每个脉冲后半段）：

$$\begin{aligned} \mathbf{Y}_f &= \mathbf{F}\,\mathbf{R} \quad (\text{快时间DFT到子载波})\\ \mathbf{U} &= \mathbf{Y}_f \odot \overline{\mathbf{X}}\quad (\text{逐子载波与已知}\;\mathbf{X}\;\text{做匹配滤波/去调制})\\ \mathbf{V} &= \mathbf{F}^H \mathbf{U}\in\mathbb{C}^{N\times M}\quad (\text{回到时域，得到按延迟对齐的相关})\\ \mathrm{AFfs} &= \left|\mathrm{FFT}_M\!\left(\mathbf{V}^T\right)\right|^2 \quad\in\mathbb{R}_+^{M\times N}. \end{aligned}$$

因此，AFfs 的解析式（逐格点功率）可写为

$$\boxed{\; \mathrm{AFfs}[k,n] = \left| \sum_{m=0}^{M-1} e^{-j2\pi \tfrac{km}{M}} \left[\mathbf{F}^H\!\left( (\mathbf{F}\mathbf{R})\odot \overline{\mathbf{X}} \right)\right]_{n,m} \right|^2 .\;}$$

把 $\mathbf{R}$ 代入 $\mathbf{z}=\mathbf{D}(\tfrac{2\nu}{MN})\Pi_\tau\mathbf{y}$ 的结构并注意复制两段的用意（$\tau<N$ 时保证“取后半段”不会跨块混叠），可得到更简单的闭式（理想对齐、无噪）：

• 快时间（延迟）维度：$\mathbf{F}^H(\mathbf{F}\cdot)$ 与 $\overline{\mathbf{X}}$ 相消后给出纯循环移位：

  $$\mathbf{V}(n,m)\;=\;C\cdot \delta_{n,\;\tau\!\!\!\pmod N}\cdot e^{j\,\phi_\nu(n)} \cdot e^{j \frac{4\pi\nu}{M} m},$$

  其中 $C$ 为常数相位/幅度，$\phi_\nu(n)$ 是与 $n$ 相关的（较慢的）相位项。

• 再做慢时 $M$-点 FFT，得到二维 δ 形或Dirichlet核：

$$\boxed{ \mathrm{AFfs}[k,n] = N^2\;\Big|\underbrace{\sum_{m=0}^{M-1}\!e^{-j2\pi m\left(\tfrac{k}{M}-\tfrac{4\nu}{M}\right)}}_{\displaystyle D_M\!\left(2\pi\frac{k-4\nu}{M}\right)}\Big|^2 \cdot \delta_{n,\;\tau\!\!\!\pmod N} }$$

即

$$\boxed{ \mathrm{AFfs}[k,n] = \begin{cases} N^2 M^2,\quad &n\equiv \tau \pmod N,\; k\equiv 4\nu \pmod M \\ N^2\,\big|D_M(2\pi\frac{k-4\nu}{M})\big|^2,\quad &\text{一般（$\nu$ 非栅格整数）} \\ 0,&\text{其余理想格点} \end{cases}}$$

其中 $D_M(\theta)=\dfrac{\sin(M\theta/2)}{\sin(\theta/2)}\,e^{-j\theta(M-1)/2}$ 为 Dirichlet 核。

要点：复制两段并只取“后半段”，把整体的“全局循环移位 $\Pi_\tau$”等效成每个脉冲内部的纯 $N$ 点循环移位 $\mathbf{P}_\tau$，从而快时间上得到严格的 δ（无泄露）；慢时间上，多普勒指数的栅格是 $k_0=4\nu \bmod M$（注意是 4 倍缩放）。

---

7) PSR 计算的解析理解

matlab
if nu~=0
    rec = 10*log10( max(AFfs) / (sum(AFfs)-max(AFfs)) * (M*N+1) );
else
    % 特判 nu=0，避免除零
end

在理想、上栅格（$\tau\in\mathbb{Z}$, $4\nu\in\mathbb{Z}$）且无噪声的情形，上式中的
$\max(\mathrm{AFfs})=N^2M^2$，而 $\sum(\mathrm{AFfs})-\max=0$（其他栅格全为 0），PSR → +∞ dB。
当 $\nu$ 不整齐落点（例如 $\nu=k_0/4+\epsilon$）时，慢时 FFT 出现 Dirichlet 泄露：

$$\sum_{k=0}^{M-1}\!|D_M(2\pi\tfrac{k-4\nu}{M})|^2 = M^2,$$

因此有

$$\mathrm{PSR}(\nu)\approx 10\log_{10}\!\left(\frac{|D_M(2\pi\epsilon)|^2}{M^2-|D_M(2\pi\epsilon)|^2}\cdot(MN+1)\right),$$

$\epsilon\to 0$ 时趋向 $+\infty$，$|\epsilon|$ 增大时 PSR 下降。

---

8) 作图与代码里一个“隐藏的坑”

matlab
range1 = -M:M; 
...
for nu = 10   % <— 实际只算了 nu=10 一个点
    ...
end
plot(range1, rec-max(rec), 'o-');
plot(range1, rec, 'o-');

这里 for nu = 10 只计算了一个多普勒点，但却用 -M:M 作为横轴作 33 个点的曲线；
因此除了第一个元素之外 rec 其余全是 0（或未填），曲线没有实际意义。
若想画“多普勒索引–PSR”的曲线，应改为：

matlab
for nu = range1              % 或者 nu = -M:M
    ...
end

并注意4 倍映射：AFfs 的峰在 $k_0=4\nu\bmod M$。例如 nu=10、M=16 时，峰在 $k_0=8$。

---

9) 把关键行再浓缩为“矩阵等式对照表”

• F = dftmtx(N)./sqrt(N)
  $\quad \mathbf{F}_{k,n}=\frac{1}{\sqrt{N}}e^{-j2\pi (k-1)(n-1)/N}$。

• y = repmat(F'*x,[2,1])
  $\quad \mathbf{Y}=\begin{bmatrix}\mathbf{F}^H\mathbf{X}\\\mathbf{F}^H\mathbf{X}\end{bmatrix}$，$\mathbf{y}=\mathrm{vec}(\mathbf{Y})$。

• yy2 = circshift(yy1,tau)
  $\quad \mathbf{y}_\tau=\Pi_\tau \mathbf{y}$。

• yy3 = yy2 .* exp(j*2π*(0:L-1)/M/N/2*nu)
  $\quad \mathbf{z}=\mathbf{D}\!\left(\tfrac{2\nu}{MN}\right)\mathbf{y}_\tau$。

• AF(i,j) 那一长串
  $\quad \mathrm{AF}(i,j)=\mathbf{xx}^H\Pi_{-i}\mathbf{S}\mathbf{D}\!\left(\tfrac{j}{L}\right)\mathbf{z}$。

• AFfs = abs( fft( (F' * ((F*y2(N+1:end,:)).*conj(x))).' ) ).^2
  $\quad \mathrm{AFfs}[k,n]=\Big|\sum_m e^{-j2\pi km/M}\,[\mathbf{F}^H((\mathbf{F}\mathbf{R})\odot\overline{\mathbf{X}})]_{n,m}\Big|^2$。

• 理想闭式（最直观）
  $\quad \mathrm{AFfs}[k,n]=N^2\big|D_M(2\pi\frac{k-4\nu}{M})\big|^2\cdot \delta_{n,\tau}$。

---

10) “最后图像”应对应的解析结论

若把 nu 从 $-M{:}M$ 扫描并正确填充 rec，在无噪声 + 栅格对齐时，理论上：

• 当 $4\nu\equiv 0\pmod M$（或更一般 $k=4\nu\bmod M$）处，$\mathrm{AFfs}$ 只有一个非零栅格（$\tau,k$），于是 $\mathrm{PSR}\to +\infty$ dB；
• 离开这些点，如果 $\nu$ 不是 $\tfrac{k}{4}$ 的整数（非上栅格），$\mathrm{AFfs}$ 的慢时切片是 Dirichlet 核，PSR 由上式给出并随 $|\epsilon|$ 增大而下降；
• 快时间维度因为 $\tau$ 为整数且复制-裁剪策略，呈严格 $\delta[n-\tau]$。

要得到“有限且平滑”的 PSR 曲线，实践中通常需要：
① 加噪或加多径；② 使用非整栅格 $\nu$；③ 对慢时加窗（抑制 Dirichlet 旁瓣）；或 ④ 去掉“两段复制 + 只取后半段”的完全对齐技巧。

---

小结（一句话版）

这段代码本质上把 $M$ 个 OFDM 脉冲按快时间$N$/慢时间$M$组织，经过循环延迟 $\tau$与归一化多普勒 $f_D=2\nu/(MN)$ 的单径信道，再用“快时匹配 + 慢时 FFT”的方式计算二维模糊函数功率。其解析形态是延迟维 δ、慢时维 Dirichlet 核；在整数对齐时，峰值落在 $(n=\tau,\;k=4\nu\bmod M)$，PSR 理论上趋于无穷。

超级好：这段代码里“多普勒调制”其实天然可分离成快时间项 × 慢时间项，把它显式拆开后，就能看清快时间多普勒如何独立且乘性地影响二维模糊函数（AF/AFfs）。

---

1) 多普勒“快/慢时间”严格可分离

代码里的多普勒调制是

$$e^{j2\pi \frac{\nu}{M N\cdot 2}\,k},\qquad k=0,1,\dots,L-1,\;\; L=2MN .$$

把全局采样索引写成

$$k=n+2N\,m,\quad n\in[0,2N\!-\!1],\; m\in[0,M\!-\!1],$$

就得到严格分离：

$$\underbrace{e^{j2\pi \frac{\nu}{2MN}\,k}}_{\text{全局}} =\underbrace{e^{j2\pi \frac{\nu}{2MN}\,n}}_{\displaystyle\text{快时间（符号内）}\;\psi_n(\nu)} \cdot \underbrace{e^{j2\pi \frac{\nu}{M}\,m}}_{\displaystyle\text{慢时间（跨脉冲）}\;\phi_m(\nu)} .$$

其中 $\phi_m(\nu)=e^{j2\pi (\nu/M)m}$ 决定慢时 FFT里的峰值“落在 $k=\nu\bmod M$”，而
$\psi_n(\nu)=e^{j2\pi (\nu/(2MN))n}$ 只在符号内快时间产生线性相位斜率。

---

2) “两段拼接 + 取后半段”后，接收块的结构

利用两段拼接并仅取每个脉冲的后 $N$ 点，可把全局循环延时 $\tau$ 化成符号内的 $N$ 点循环移位。记

• $s_m=F^H x_m\in\mathbb{C}^N$ 为第 $m$ 个 OFDM 符号的时域列；
• $\mathbf{P}_\tau$ 为 $N$ 点循环移位矩阵；
• 令我们取到的 $N$ 点“后半段”行索引写成 $\tilde n=0,\dots,N-1$（对应全局 $n=N+\tilde n$）。

则每个脉冲得到的 $N\!\times\!1$ 向量（记为 $r_m$）可写成

$$r_m = e^{j2\pi \frac{\nu}{M}m}\;\underbrace{e^{j2\pi \frac{\nu}{2M}}}_{\text{常相位}} \;\underbrace{\mathrm{diag}\!\Big(e^{j2\pi \frac{\nu}{2MN}\tilde n}\Big)}_{\text{快时间相位斜率}} \;\mathbf{P}_\tau\,s_m .$$

---

3) 快时 FFT + 频域匹配滤波后的“核”分解

快时 FFT：$y_m = F r_m$。代入上式并用 $s_m=F^H x_m$，得

$$y_m = e^{j2\pi \frac{\nu}{M}m}\;e^{j2\pi \frac{\nu}{2M}} \;\underbrace{F\,\mathrm{diag}\!\Big(e^{j2\pi \tfrac{\nu}{2MN}\tilde n}\Big)F^H}_{\displaystyle \mathbf{G}(\nu)} \;\underbrace{\mathrm{diag}\!\Big(e^{-j2\pi \tfrac{\tau}{N}k}\Big)}_{\displaystyle \Lambda_\tau} \,x_m .$$

随后代码做了频域匹配 $\;u_m = y_m \odot \overline{x_m}$。因此

$$u_m = e^{j2\pi \frac{\nu}{M}m}\;e^{j2\pi \frac{\nu}{2M}}\;\Big(\mathbf{G}(\nu)\Lambda_\tau x_m\Big)\odot \overline{x_m}.$$

• $\mathbf{G}(\nu)=F\,\mathrm{diag}(e^{j2\pi (\nu/(2MN))\tilde n})F^H$ 是由快时间多普勒引入的 ICI 核（频域循环矩阵，$\nu=0$ 时为单位阵）。
• $\Lambda_\tau=\mathrm{diag}(e^{-j2\pi \tau k/N})$ 是延迟相位因子。

这一步已经把“慢时间多普勒”$\phi_m(\nu)$ 与“快时间多普勒”$\psi_n(\nu)$ 完全分离到两个独立因子：
$\phi_m(\nu)$ 仅是 $m$ 的纯相位；$\mathbf{G}(\nu)$ 则只作用在每个脉冲内部并引入 ICI。

---

4) 统计意义下（随机 16-PSK）的简洁解析解

对独立随机相位的 16-PSK 做平均（或足够多蒙特卡洛），有

$$\mathbb{E}_x\!\left[\big(\mathbf{G}\Lambda_\tau x_m\big)\odot \overline{x_m}\right] =\mathrm{diag}\!\big(\mathbf{G}(\nu)\big)\;e^{-j2\pi \tfrac{\tau}{N}k}.$$

而 $\mathrm{diag}(\mathbf{G}(\nu))$ 对任意单位化 DFT 恰为常数（与 $k$ 无关）：

$$\mathrm{diag}\!\big(\mathbf{G}(\nu)\big)= \frac{1}{N}\sum_{\tilde n=0}^{N-1} e^{j2\pi \frac{\nu}{2MN}\tilde n} =\underbrace{e^{j\pi \frac{\nu}{2MN}(N-1)}}_{\text{无关幅度}} \;\frac{\sin\!\big(\pi \frac{\nu}{2M}\big)}{N\sin\!\big(\pi \frac{\nu}{2MN}\big)} \;\triangleq\; C_\text{fast}(\nu).$$

于是（去掉无关全局相位）：

$$\boxed{ \mathbb{E}[u_m]\;=\;C_\text{fast}(\nu)\;e^{j2\pi \frac{\nu}{M}m}\;e^{-j2\pi \tfrac{\tau}{N}k}. }$$

把它送回快时 IFFT（代码里的 $F'^{}$）到时域，再做慢时 $M$-点 FFT（代码里的 fft(·,M)），得到二维功率：

$$\boxed{ \mathrm{AF}_{\!fs}[k,n]\;\approx\; N^2\,|C_\text{fast}(\nu)|^2\; \big|D_M\!\big(2\pi\tfrac{k-\nu}{M}\big)\big|^2\;\cdot\;\delta_{n,\;\tau\!\!\!\!\pmod N}. }$$

其中 $D_M(\cdot)$ 为 Dirichlet 核。这正是“慢时泄露 × 快时缩放”的乘性分离。

快时间多普勒的独立作用完全浓缩在

$$\boxed{ G_\text{fast}(\nu)\;\triangleq\;|C_\text{fast}(\nu)|^2 =\left[\frac{\sin\!\big(\pi \frac{\nu}{2M}\big)}{N\sin\!\big(\pi \frac{\nu}{2MN}\big)}\right]^2 }$$

——它是一个仅与 $\nu$ 与 $(M,N)$ 有关的“快时多普勒损失因子”：

• $\nu\!\to\!0$ 时 $G_\text{fast}(\nu)\!\to\!1$；
• $\nu$ 增大时呈sinc 形衰减，但不改变慢时的 Dirichlet 形状与峰位。

---

5) 物理与工程解读（一眼看懂）

• 慢时间多普勒（跨脉冲）：只给每个脉冲乘上 $e^{j2\pi(\nu/M)m}$，因此慢时 FFT 就得到 $\big|D_M(2\pi (k-\nu)/M)\big|^2$——决定 Doppler 轴上的峰位与旁瓣。
• 快时间多普勒（符号内）：只在符号内加线性相位斜率，经“快时 FFT → 匹配滤波 → IFFT”后，不会展宽 Doppler 轴，而是以 $G_\text{fast}(\nu)$ 等比例缩放整个二维 AF（并保持延迟轴上的 $\delta[n-\tau]$）。
• 上式解释了你代码里：复制两段并只取后半段，使延迟轴保持 δ 形；而 Doppler 轴的形状只由慢时间决定，快时间只贡献幅度缩放。

---

6) 小小更正（和上次回答的差异）

之前我把慢时相位步进写成 $\Delta\phi_m=\tfrac{4\pi\nu}{M}$，那是笔误。
本代码的步进是

$$\Delta\phi_m = 2\pi \frac{\nu}{L}\cdot(2N)=\frac{2\pi\nu}{M},$$

所以峰位在 $k=\nu\bmod M$（不是 $4\nu$）。

---

结论（一句话）

把多普勒拆成 $\psi_n(\nu)\times\phi_m(\nu)$ 后，二维 AF 的解析解可写成

$$\mathrm{AF}_{\!fs}[k,n]\approx N^2\,\underbrace{G_\text{fast}(\nu)}_{\text{快时多普勒“幅度门”}}\times \underbrace{\big|D_M(2\pi\tfrac{k-\nu}{M})\big|^2}_{\text{慢时多普勒“形状+峰位”}}\times \delta_{n,\tau},$$

清楚地表明：快时间多普勒只作为一个与 $\nu$ 相关的乘性衰减，独立地缩放整个二维模糊函数；慢时间多普勒决定 Doppler 轴的峰位与旁瓣形状。

太好了！下面给你一段可直接运行的 MATLAB 脚本，用“解析解（可分离的快/慢时间多普勒模型）”去复现并对比你原始代码里算出来的 $\mathrm{AF}_{fs}$。
脚本做三件事：

1. 按你原始流程仿真 $\mathrm{AF}_{fs}^{\text{sim}}$；
2. 用我们推导的解析式生成 $\mathrm{AF}_{fs}^{\text{ana}}$（快时只给幅度缩放 $G_\text{fast}(\nu)$，慢时给 Dirichlet 核）；
3. 画图 + 叠加对比 + 误差指标（RMSE / 峰值相对误差）。

说明：解析式里多普勒的慢时“峰位”在 $k=\nu \bmod M$；快时的影响是

$$G_\text{fast}(\nu)=\Big[\tfrac{\mathrm{sinc}\!\big(\nu/(2M)\big)}{\mathrm{sinc}\!\big(\nu/(2MN)\big)}\Big]^2,$$

并保持延迟轴为 $\delta[n-\tau]$。为了公平比较，我们对两者都做按最大值归一化。

---

MATLAB 脚本（复制即可运行）

matlab
%% 对比：解析AFfs vs 仿真AFfs
clc; clear; close all; format long;

% ---------- 参数区 ----------
times = 1;
M = 16;
N = 32*times;
L = 2*M*N;
Nmc = 10;             % Monte-Carlo 次数（可适当加大，更接近解析式）
rng(1545);
tau = 4;              % 整数延迟
nuList = [10];        % 要对比的多普勒索引，可改为 [-M:M] 等

F = dftmtx(N)/sqrt(N);           % 单位化DFT
kL = (0:L-1).';                   % 全局采样索引，用于调制
kM = (0:M-1).';                   % 慢时Doppler bin
nN = 0:N-1;                       % 快时延迟bin

% ---------- 工具函数（匿名） ----------
sincN = @(x) sinc(x);   % MATLAB 的 sinc(x) = sin(pi x)/(pi x)

% Dirichlet 核的 |D_M(2π(k-ν)/M)|^2，数值稳定写法：M^2 * [sinc(k-ν)/sinc((k-ν)/M)]^2
dirichlet2 = @(M, nu) ( M * sincN(kM - nu) ./ sincN( (kM - nu)/M ) ).^2;

% 快时多普勒增益因子 G_fast(ν) = [sinc(ν/(2M))/sinc(ν/(2MN))]^2
Gfast = @(nu) ( sincN(nu/(2*M)) ./ sincN(nu/(2*M*N)) ).^2;

% ---------- 主循环 ----------
for idx = 1:numel(nuList)
    nu = nuList(idx);

    % ---- 1) 仿真 AFfs（按你的代码流程） ----
    AFfs_sim = zeros(M, N);
    for mc = 1:Nmc
        % 随机16-PSK
        s = randi(16,[N,M]) - 1;
        X = pskmod(s,16);            % N x M
        S = F' * X;                  % N x M
        Y = [S; S];                  % 2N x M

        % 串行化 + 延迟 + 多普勒
        yy1 = Y(:);
        yy2 = circshift(yy1, tau);
        phs = exp(1j*2*pi*(nu/(2*M*N)) * kL);
        yy3 = yy2 .* phs;

        % 重排成 (2N) x M，并取后半段 N 点
        Z  = reshape(yy3, 2*N, M);
        R  = Z(N+1:end, :);          % N x M

        % 快时FFT -> 频域匹配 -> 快时IFFT -> 慢时FFT
        tmp = F * R;                 % N x M
        tmp = tmp .* conj(X);        % N x M
        V   = F' * tmp;              % N x M（快时相关）
        AF  = abs( fft( V.', [], 1 ) ).^2;   % M x N（慢时FFT，默认沿行）
        AFfs_sim = AFfs_sim + AF;
    end
    AFfs_sim = AFfs_sim / Nmc;

    % ---- 2) 解析 AFfs ----
    % 仅在 n = tau 这一列有能量（理想δ），慢时形状为 Dirichlet 核；快时为整体幅度门 G_fast
    AFfs_ana = zeros(M, N);
    AF_dop   = dirichlet2(M, nu);       % 大小 M x 1
    gain_fast = Gfast(nu);              % 标量
    AF_col = gain_fast * AF_dop;        % M x 1
    AFfs_ana(:, tau+1) = AF_col;        % 放到 n = tau 列（MATLAB 1-based）

    % 为了可比性：都按最大值归一化
    AFfs_sim_n  = AFfs_sim  / max(AFfs_sim(:)  + eps);
    AFfs_ana_n  = AFfs_ana  / max(AFfs_ana(:)  + eps);

    % ---- 3) 误差指标 ----
    rmse = sqrt( mean( (AFfs_sim_n(:) - AFfs_ana_n(:)).^2 ) );
    peakSim = max(AFfs_sim_n(:));
    peakAna = max(AFfs_ana_n(:));
    peakRelErr = abs(peakSim - peakAna) / max(peakSim, 1e-12);

    % ---- 4) 画图对比 ----
    figure('Name', sprintf('nu = %g', nu), 'Color', 'w');
    tiledlayout(2,3, 'Padding','compact','TileSpacing','compact');

    % 仿真面
    nexttile; imagesc(nN, kM, 20*log10(AFfs_sim_n + eps));
    axis xy; xlabel('Delay bin n'); ylabel('Doppler bin k');
    title(sprintf('Sim AF_{fs} (nu=%g)', nu)); colorbar;

    % 解析面
    nexttile; imagesc(nN, kM, 20*log10(AFfs_ana_n + eps));
    axis xy; xlabel('Delay bin n'); ylabel('Doppler bin k');
    title('Analytic AF_{fs}'); colorbar;

    % 误差面
    nexttile; imagesc(nN, kM, AFfs_sim_n - AFfs_ana_n);
    axis xy; xlabel('Delay bin n'); ylabel('Doppler bin k');
    title(sprintf('Diff (RMSE=%.3g)', rmse)); colorbar;

    % 慢时切片（固定 n=tau）
    nexttile; 
    plot(kM, AFfs_sim_n(:, tau+1), 'o-', 'DisplayName','Sim'); hold on;
    plot(kM, AFfs_ana_n(:, tau+1), 's--', 'DisplayName','Analytic'); grid on;
    xlabel('Doppler bin k'); ylabel('Normalized Power');
    title('Slice @ delay n = \tau');
    legend;

    % 快时切片（峰值 k 处）
    [~, kpk] = max(AFfs_sim_n(:, tau+1));
    nexttile;
    plot(nN, AFfs_sim_n(kpk, :), 'o-', 'DisplayName','Sim'); hold on;
    plot(nN, AFfs_ana_n(kpk, :), 's--', 'DisplayName','Analytic'); grid on;
    xlabel('Delay bin n'); ylabel('Normalized Power');
    title(sprintf('Slice @ k = argmax (k=%d)', kpk-1));
    legend;

    % 快时“幅度门” vs 归一化峰值
    nexttile;
    % 如果你扫多个 nu，可把点连起来；这里只画单点展示
    stem(nu, gain_fast / max(gain_fast,eps), 'filled');
    xlabel('nu'); ylabel('G_{fast}(\nu) (norm)');
    title('Fast-time Doppler gain (analytic)');
end

---

你会看到什么

• 延迟维：仿真与解析几乎都只在 $n=\tau$ 这一列有能量（δ 形），吻合我们的拆分结论。
• 多普勒维：两者在 $k=\nu\bmod M$ 附近的主瓣位置与旁瓣（Dirichlet 核）一致；RMSE 会随 $N_{mc}$ 增大而更小。
• 快时多普勒作用：解析面的整体幅度随 $G_\text{fast}(\nu)$ 变化；仿真面的峰值（归一化前）与该因子一致（小样本有微小偏差）。

如果想看一整条 PSR 曲线，把 nuList 改成 -M:M 循环，分别计算两者的 PSR 并叠加即可；随着 Nmc 增大，解析与仿真会逐渐重合。

确实有点“对不上”的感觉，主要是因为：

1. 你那张图左上（仿真）像“噪声面”，是 Monte-Carlo 次数太少 + 两幅图的色轴不一致 造成的观感偏差；
2. 我之前给的是“期望意义的解析式”（含 $G_{\text{fast}}(\nu)$），而你的仿真是对某一 realization 的数据在做运算；两者在有限样本时会有 speckle。

下面给一份改进脚本，一次性给出两种解析对比：

• INST（瞬时解析）：用同一批数据 $X$ 做“代数化”的解析结果，和仿真应零误差对齐（数值误差级）；
• AVG（期望解析）：就是我们推导的 $G_{\text{fast}}(\nu)\times$ Dirichlet 核，只在 $n=\tau$ 列非零，便于看快/慢时间的分离效应。
  脚本也把两幅面图的 caxis 统一，并给出 RMSE。

---

MATLAB（直接运行）

matlab
%% AF_fs: simulation vs analytic (INST & AVG)
clc; clear; close all; format long;

% ===== Params =====
times = 1;
M = 16;
N = 32*times;
L = 2*M*N;
rng(1545);
tau = 4;                 % integer delay (bin)
nu  = 10;                % Doppler index (bin), k_peak = nu mod M
Nmc = 50;                % MC次数；设大一点 speckle 会减小

F  = dftmtx(N)/sqrt(N);
kL = (0:L-1).';
kM = (0:M-1).';
nN = (0:N-1);

% ----- helpers -----
sincN = @(x) sinc(x);  % MATLAB: sin(pi x)/(pi x)

% Dirichlet^2: | sum_{m=0}^{M-1} e^{-j 2pi (k-nu) m / M} |^2
dirich2 = @(M,nu) ( M * sincN(kM - nu) ./ sincN((kM - nu)/M) ).^2;

% Fast-time gain (avg model)
Gfast = @(nu) ( sincN(nu/(2*M)) ./ sincN(nu/(2*M*N)) ).^2;

% ----- core functions -----
% simulation with a given X (完全按原流程)
af_sim_givenX = @(X) ...
    ( ...
        ( ...
            fft( ...
                (F' * ((F * reshape( ...
                    ( ...
                        circshift( ...
                            [ (F'*X); (F'*X) ](:), tau ) ...
                        .* exp(1j*2*pi*(nu/(2*M*N))*kL) ...
                    , 2*N, M) )(N+1:end,:) ).*conj(X)) ).' ...
            ).^2 ...
        ) ...
    );

% instantaneous analytic for the same X（代数化：把信道作用拆成核）
function AF = af_ana_inst(F,M,N,tau,nu,X)
    ntilde = 0:N-1;
    % 快时间相位斜率（截取的是后半段 N 点，因此多了常相位 e^{j2πν/(2M)})
    Dfast = diag(exp(1j*2*pi*(nu/(2*M*N))*ntilde));      % N×N
    G     = F * Dfast * F.';                             % ICI核（频域循环矩阵）
    LamT  = diag(exp(-1j*2*pi*(tau/N)*(0:N-1)));         % 延迟相位
    ph_const = exp(1j*2*pi*(nu/(2*M)));
    ph_slow  = exp(1j*2*pi*(nu/M)*(0:M-1));              % 1×M

    Yf = ph_const * (G * LamT * X) .* (ones(N,1) * ph_slow); % N×M
    U  = Yf .* conj(X);
    V  = F' * U;                                          % N×M
    AF = abs( fft(V.', [], 1) ).^2;                       % M×N
end

% average (expectation) analytic: 只在 n=tau 列非零
function AF = af_ana_avg(M,N,tau,nu,GfastFun,dirichFun)
    AF = zeros(M,N);
    AF(:,tau+1) = GfastFun(nu) * dirichFun(M,nu);
end

% ===== Monte Carlo accumulate =====
AF_sim = zeros(M,N);
AF_ins = zeros(M,N);
for mc = 1:Nmc
    s  = randi(16,[N,M]) - 1;
    X  = pskmod(s,16);          % N×M
    AF_sim = AF_sim + af_sim_givenX(X);
    AF_ins = AF_ins + af_ana_inst(F,M,N,tau,nu,X);
end
AF_sim = AF_sim / Nmc;
AF_ins = AF_ins / Nmc;

AF_avg = af_ana_avg(M,N,tau,nu,Gfast,@dirich2);

% 归一化到各自峰值（也可统一用仿真峰值，看你想比什么）
mx = max(AF_sim(:)) + eps;
AFn_sim = AF_sim / mx;
AFn_ins = AF_ins / mx;
AFn_avg = AF_avg / mx;

% ===== Metrics =====
rmse_inst = sqrt(mean((AFn_sim(:)-AFn_ins(:)).^2));
rmse_avg  = sqrt(mean((AFn_sim(:)-AFn_avg(:)).^2));

% ===== Plots（统一色轴，避免视觉误差）=====
lowdB = -60; highdB = 0;
figure('Color','w','Name',sprintf('nu=%g',nu));
tiledlayout(2,3,'Padding','compact','TileSpacing','compact');

nexttile; imagesc(nN,kM,20*log10(AFn_sim+eps)); axis xy;
caxis([lowdB highdB]); colorbar; xlabel('Delay bin n'); ylabel('Doppler bin k');
title(sprintf('Sim AF_{fs} (nu=%g)',nu));

nexttile; imagesc(nN,kM,20*log10(AFn_ins+eps)); axis xy;
caxis([lowdB highdB]); colorbar; xlabel('Delay bin n'); ylabel('Doppler bin k');
title(sprintf('Analytic INST (RMSE=%.2g)',rmse_inst));

nexttile; imagesc(nN,kM,20*log10(AFn_avg+eps)); axis xy;
caxis([lowdB highdB]); colorbar; xlabel('Delay bin n'); ylabel('Doppler bin k');
title(sprintf('Analytic AVG (RMSE=%.2g)',rmse_avg));

% delay slice @ n = tau
nexttile;
plot(kM, AFn_sim(:,tau+1), 'o-','DisplayName','Sim'); hold on;
plot(kM, AFn_ins(:,tau+1), 's--','DisplayName','Analytic INST');
plot(kM, AFn_avg(:,tau+1), 'x-.','DisplayName','Analytic AVG');
grid on; xlabel('Doppler bin k'); ylabel('Norm power');
title('Slice @ delay n = \tau'); legend('Location','best');

% doppler slice @ k = argmax
[~,kpk] = max(AFn_sim(:,tau+1));
nexttile;
plot(nN, AFn_sim(kpk,:), 'o-','DisplayName','Sim'); hold on;
plot(nN, AFn_ins(kpk,:), 's--','DisplayName','Analytic INST');
plot(nN, AFn_avg(kpk,:), 'x-.','DisplayName','Analytic AVG');
grid on; xlabel('Delay bin n'); ylabel('Norm power');
title(sprintf('Slice @ k = %d', kpk-1)); legend('Location','best');

% 快时增益（期望模型）
nexttile;
stem(nu, Gfast(nu),'filled'); grid on;
xlabel('nu'); ylabel('G_{fast}(\nu)'); title('Fast-time gain (AVG model)');

---

怎么检查“问题点”

• 色轴要统一：我在上面用 caxis([-60 0]) 固定，避免视觉上“仿真有噪声、解析全黑”的错觉。
• 瞬时解析 vs 仿真：Analytic INST 与 Sim 的 RMSE 应接近 0（数值级），两幅面图基本不可分。
• 期望解析（AVG）：它只在 $n=\tau$ 一列非零；与仿真在 MC 很小时会有些偏差（你的那张“花花的”就是 speckle），增大 Nmc 就会逐步贴近。
• 峰位：此代码峰在 $k=\nu\bmod M$，比如 nu=10、M=16 时，峰在 $k=10$（和你图左下相符）。

如果你仍看到差异，通常是下面这些小点：

1. 参数不一致（比如把 nuList 或 tau 改了）；
2. MATLAB fft 的维度方向；上面已显式 fft(V.',[],1)，等价于对慢时做 $M$-点 FFT；
3. pskmod 的 M-PSK 初相位设置。上面用默认（与原代码一致）。

需要我把 nu 扫成 $-M{:}M$ 并画整条曲线/PSR 也可以，直接在这份脚本上加一层循环就行。