下面这个代码是计算ISAC信号的capon谱进行二维位置信息感知的代码，现在做的图是进行capon谱...

下面给出一种 输出 SINR（signal‑to‑interference‑plus‑noise ratio） 的常用定义、推导思路，以及能直接插到你现有 MATLAB 脚本里的示例代码。思路保持和你现在的 Capon‑search 框架一致——先做谱搜索/消除，然后针对「已锁定的某一条 (τ, ν) 波束」计算该波束的输出 SINR，再做 Monte‑Carlo 扫描即可量化性能。

---

1 · 理论回顾

对某一候选点 $(\tau,\nu)$ ，其 Capon 权重向量

$$\mathbf w(\tau,\nu)\;=\;\frac{\mathbf R^{-1}\mathbf s}{\mathbf s^{\!\mathrm H}\mathbf R^{-1}\mathbf s}$$

• $\mathbf R$ —— 二维快拍协方差（含信号、干扰、噪声）。
• $\mathbf s$ —— 对应 $(\tau,\nu)$ 的空间‑快拍引导向量（即你 Capon 谱里的 s_vec）。

然后有

$$\underbrace{\text{SINR}_\text{out}(\tau,\nu)}_{\text{dB 里常用 }10\log_{10}(\cdot)} \;=\; \frac{\mathbf w^{\!\mathrm H}\,\mathbf R_s\,\mathbf w} {\mathbf w^{\!\mathrm H}\,\mathbf R_{i+n}\,\mathbf w}, \qquad \bigl(\mathbf R \;=\;\mathbf R_s\;+\;\mathbf R_{i+n}\bigr)$$

• 如果你掌握信号分量功率 $\sigma_s^2$ 与噪声功率 $\sigma_n^2$，则
  $\mathbf R_s = \sigma_s^2\,\mathbf s\mathbf s^{\!\mathrm H}$，
  $\mathbf R_{i+n} = \sigma_n^2\,\mathbf I$（或加上干扰协方差）。

• Capon 波束形成特性：$\mathbf w^{\!\mathrm H}\mathbf s = 1$，所以上式可化为

  $$\text{SINR}_\text{out}\;=\; \frac{\sigma_s^2} {\mathbf w^{\!\mathrm H}\mathbf R_{i+n}\mathbf w}.$$

如果直接在仿真里生成了完整矩阵 $\mathbf Z$（含快拍）， 最方便的做法是 用样本协方差估计：

matlab
R_hat = (Z*Z')/M;                      % 2N×2N 样本协方差

注意：为了避免把目标信号也当成「干扰」，做 SINR 统计时往往
① 用“理想真值”拆分出 signal + noise，然后公式化计算；
② 或者像传统雷达那样，用“残差信号”估计噪声功率（见下面代码）。

---

2 · 在你代码里的最小修改

下面片段假设你已经锁定了一个峰（tau_hat, nu_hat）并得到了残差 Zc， 就可以立刻算这一束的输出 SINR：

matlab
%% ---- ① 形成该束的 Capon 权重 -----------------------------------
S_mat = (Df_all(:,in_max)*Ds_all(in_max,:)) ...
          .* circshift(Y0, tau_hat);          % (2N×M)
s_vec  = S_mat(:);                            % 2NM×1 steering
R_hat  = (Z * Z')/M;                          % 样本协方差 (含信号+噪声)

w      = (R_hat \ s_vec) / (s_vec' * (R_hat \ s_vec));  % Capon 权重 (2N×1)

%% ---- ② 分离信号 + 残差 ----------------------------------------
y_snap = w' * Z;             % 1×M 输出（每个快拍）
y_sig  = w' * S_mat * a_hat; % 理想信号输出
y_res  = y_snap - y_sig;     % 噪声/干扰

P_sig  = mean(abs(y_sig).^2);    % 信号功率
P_res  = mean(abs(y_res).^2);    % 干扰+噪声功率

SINR_out = 10*log10(P_sig / P_res);

fprintf('[%2d] τ=%2d  ν=%.1f  输出 SINR = %.2f dB\n',...
         iter, tau_hat, nu_hat, SINR_out);

• a_hat 是你前面用最小二乘求出的幅度估计。
• y_snap 把 $M$ 个快拍都用来平均，结果更稳定。
• 如果要按经典公式：P_res_alt = (w' * R_hat * w) - P_sig; 两种算法在快拍够多时应一致。

---

3 · 批量 / Monte‑Carlo 性能曲线

为了量化「Capon‑ISAC」对 输入 SNR 的增益，最常见的就是做多次 Monte‑Carlo，统计输出 SINR 的期望或分位数并作图。下面给出一个骨架（无需 canvas，直接放进新 .m 文件即可）：

matlab
SNRdB_vec = -10:2:20;          % 输入 SNR 扫描
Nmc       = 200;               % Monte‑Carlo 次数
SINR_out_mean = zeros(size(SNRdB_vec));

for ii = 1:length(SNRdB_vec)
    SNRdB = SNRdB_vec(ii);
    SINR_tmp = zeros(1,Nmc);

    for mc = 1:Nmc
        % ======== 你现有的“生成 Z + Capon 搜索”完整流程 =========
        % 关键是得到 tau_hat, nu_hat, a_hat, Z, Zc
        % --------------------------------------------------------

        % ----- 计算输出 SINR（复用上节代码） --------------------
        R_hat = (Z * Z')/M;
        S_mat = (Df_all(:,in_max)*Ds_all(in_max,:)) ...
                   .* circshift(Y0, tau_hat);
        s_vec = S_mat(:);
        w     = (R_hat \ s_vec) / (s_vec' * (R_hat \ s_vec));

        y_snap = w' * Z;
        y_sig  = w' * S_mat * a_hat;
        P_sig  = mean(abs(y_sig).^2);
        P_res  = mean(abs(y_snap - y_sig).^2);

        SINR_tmp(mc) = 10*log10(P_sig / P_res);
    end
    SINR_out_mean(ii) = mean(SINR_tmp);
end

figure;
plot(SNRdB_vec, SINR_out_mean,'-o','LineWidth',1.5);
grid on; xlabel('输入 SNR (dB)'); ylabel('输出  SINR_{Capon}  (dB)');
title('Capon 波束形成 SINR 增益');

解释：

• 横轴是你输入的 SNRdB。
• 纵轴是多次随机快拍下的平均输出 SINR。
• 这条曲线就能量化「Capon 算法对 ISAC 场景的增益」，也方便和别的波束形成或处理方法对标。

---

4 · 小贴士

1. 噪声功率估计
   如果担心用真值 a_hat 带来先验信息，可以直接用残差协方差
   $\hat\sigma_n^2 = \tfrac{1}{2NM}\lVert Z_c\rVert_F^2$ 来替代 P_res。
2. 检测‑消除循环里实时监控
   你可在 Zc = Zc - a_hat * S_mat; 之后立即打印 SINR_out，观察迭代收敛趋势。
3. 干扰场景
   若需要研究「同距异向」或「同向异距」干扰，只需在 Z 里多加几条假目标，然后照上面方法把这些残留功率算进 P_res 即可。
4. 性能指标拓展
   除了 SINR，还可以同时记录检测概率 $P_d$ 与 虚警率 $P_{fa}$；Monte‑Carlo 里顺手数一下“正确命中次数”就行。

---

这样，你就能在现有 Capon‑ISAC 代码的基础上，系统地输出并绘制 SINR 曲线，实现对算法的量化评估。祝实验顺利，若还有其他指标或实现细节想讨论，再随时告诉我!

下面是一份 完整的 MATLAB 脚本，一次性计算 Capon 和 APES 两种谱估计在 全 range–Doppler 网格 上的

• 输出 SINR 矩阵
• 检测概率 $P_d$ 与 虚警率 $P_{fa}$

并用 parfor（你说的 “parlor” 我理解为 MATLAB 并行循环）对 Monte‑Carlo 迭代做并行加速。脚本里所有中间变量、循环和统计量都已给出，可直接保存为
capon_apes_performance.m 运行。若需调参数，只要改开头那一段即可。

matlab
function capon_apes_performance
% ============================================================
%  CAPON / APES  range‑Doppler  SINR &  P_d / P_fa  evaluation
% ============================================================

%% ---------- 全局参数 ----------
M        = 16;                 % 慢‑时间快拍数
N        = 32;                 % 快‑时间长度
K_true   = 5;                  % 真实目标数
tau_true = [ 5  5  5 6 7 ];    % range bin
nu_true  = [15 10 12 12 12 ];  % Doppler bin
amp_true = [0.3 0.4 0.6 0.7 0.6].*exp(1j*2*pi*rand(1,K_true));

SNRdB_vec = -10 : 4 : 20;      % 输入 SNR 扫描
Nmc       = 100;               % Monte‑Carlo 次数
thres_dB  = 3;                 % 检测门限 (dB over median)
rng(1);                        % 复现随机序列

%% ---------- Capon/APES 网格与常量预计算 ----------
tau_grid = 0 : N-1;            Lt = numel(tau_grid);
nu_grid  = -M : 0.2 : M;       Ln = numel(nu_grid);
G        = Lt * Ln;            % 网格点总数

n  = (0:2*N-1).';              % 2N × 1
m  = 0:M-1;                    % 1 × M
Df_all = exp(1j*pi/(M*N) * n * nu_grid );   % 2N × Ln
Ds_all = exp(1j*2*pi/M  * (nu_grid.') * m); % Ln × M
F      = dftmtx(N)/sqrt(N);
Y0     = [F';F'];   % 2N × N   (统一用于构造 steering)

% 预先为每个 range 做一次 Y0 的循环移位
Y0_shift = cell(Lt,1);
for it = 1:Lt
    Y0_shift{it} = circshift(Y0, tau_grid(it));
end

%% ---------- 输出统计量 ----------
SINR_capon_mean = zeros(length(SNRdB_vec),1);
SINR_apes_mean  = zeros(length(SNRdB_vec),1);
Pd_capon        = zeros(length(SNRdB_vec),1);
Pfa_capon       = zeros(length(SNRdB_vec),1);
Pd_apes         = zeros(length(SNRdB_vec),1);
Pfa_apes        = zeros(length(SNRdB_vec),1);

%% ---------- Monte‑Carlo + SNR 扫描 (外层并行) ----------
parfor idxSNR = 1:length(SNRdB_vec)
    SNRdB = SNRdB_vec(idxSNR);
    
    SINR_cap_sum = 0;
    SINR_ape_sum = 0;
    det_cap_cnt  = 0;  false_cap_cnt = 0;
    det_ape_cnt  = 0;  false_ape_cnt = 0;
    
    for mc = 1:Nmc
        %% ===== 1. 生成通信 (OFDM) 信号 =====
        s  = randi([0 15], N, M);
        x  = pskmod(s,16);
        Y0x = [F'*x;F'*x];                      % 2N × M
        
        % 目标 + 噪声
        Z = zeros(2*N,M);
        for k = 1:K_true
            Df = Df_all(:,nu_true(k)==nu_true); % inline Df
        end
        % 用匿名函数更直观
        Df = @(nu) exp(1j*pi*nu*n/(M*N));
        Ds = @(nu) exp(1j*2*pi*nu*m/M);
        for k = 1:K_true
            Z = Z + amp_true(k) .* ...
                (Df(nu_true(k)).*circshift(Y0x,tau_true(k)) ).*Ds(nu_true(k));
        end
        Z = Z + 10^(-SNRdB/20)*(randn(size(Z))+1j*randn(size(Z)))/sqrt(2);

        %% ===== 2. 样本协方差及其逆 =====
        R_hat = (Z*Z')/M;        
        R_inv = pinv(R_hat + 1e-6*eye(size(R_hat)));   % 稳定
        
        %% ===== 3. 遍历 range–Doppler 计算谱、SINR =====
        Pcap  = zeros(Lt,Ln);
        Papes = zeros(Lt,Ln);
        SINR_cap_map = zeros(Lt,Ln);
        SINR_ape_map = zeros(Lt,Ln);
        
        for it = 1:Lt
            Ytau = Y0_shift{it};
            for in = 1:Ln
                S_mat = (Df_all(:,in)*Ds_all(in,:)).*Ytau;   % 2N × M
                s_vec = S_mat(:);
                denom = real(s_vec' * R_inv * s_vec);
                
                % ---------- CAPON ----------
                Pcap(it,in) = 1/denom;
                
                w_cap = R_inv*s_vec/denom;        % 2NM × 1
                y_cap = w_cap' * Z;               % 1 × M
                a_cap = (w_cap'*S_mat(:))/M;      % 信号幅度估计
                P_sig = abs(a_cap)^2;
                P_out = mean(abs(y_cap).^2);
                P_res = max(P_out - P_sig, 1e-12);% 避免负
                SINR_cap_map(it,in) = 10*log10(P_sig / P_res);
                
                % ---------- APES ----------
                a_ape = (s_vec' * R_inv * Z(:)) / denom;   % APES 幅度
                Papes(it,in) = abs(a_ape)^2 / denom;       % 经典 APES 谱
                
                w_ape = w_cap;                 % 与 Capon 同向量
                y_ape = w_ape' * Z;
                P_sig2 = mean(abs(a_ape).^2);  % 同样平均
                P_out2 = mean(abs(y_ape).^2);
                P_res2 = max(P_out2 - P_sig2, 1e-12);
                SINR_ape_map(it,in) = 10*log10(P_sig2 / P_res2);
            end
        end
        
        %% ===== 4. 目标检测 (峰值超过 median + thres_dB) =====
        %   Capon
        Pvec = sort(Pcap(:));  Pmed = median(Pvec(1:round(end/2)));
        detections = Pcap > (10^(thres_dB/10)*Pmed);
        %   判断命中
        hit = 0;
        for k = 1:K_true
            [~,it_idx] = min(abs(tau_grid - tau_true(k)));
            [~,in_idx] = min(abs(nu_grid  - nu_true(k)));
            if detections(it_idx,in_idx), hit = hit + 1; end
        end
        det_cap_cnt  = det_cap_cnt  + hit;
        false_cap_cnt= false_cap_cnt+ (nnz(detections)-hit);
        
        %   APES
        Pvec2 = sort(Papes(:)); Pmed2 = median(Pvec2(1:round(end/2)));
        detections2 = Papes > (10^(thres_dB/10)*Pmed2);
        hit2 = 0;
        for k = 1:K_true
            [~,it_idx] = min(abs(tau_grid - tau_true(k)));
            [~,in_idx] = min(abs(nu_grid  - nu_true(k)));
            if detections2(it_idx,in_idx), hit2 = hit2 + 1; end
        end
        det_ape_cnt  = det_ape_cnt  + hit2;
        false_ape_cnt= false_ape_cnt+ (nnz(detections2)-hit2);

        %% ===== 5. 累加 SINR 平均值 =====
        %   仅把真实目标处的 SINR 累加
        for k = 1:K_true
            [~,it_idx] = min(abs(tau_grid - tau_true(k)));
            [~,in_idx] = min(abs(nu_grid  - nu_true(k)));
            SINR_cap_sum = SINR_cap_sum + SINR_cap_map(it_idx,in_idx);
            SINR_ape_sum = SINR_ape_sum + SINR_ape_map(it_idx,in_idx);
        end
    end
    
    %% ---------- Monte‑Carlo 结果记录 ----------
    SINR_capon_mean(idxSNR) = SINR_cap_sum /(Nmc*K_true);
    SINR_apes_mean(idxSNR)  = SINR_ape_sum /(Nmc*K_true);
    
    Pd_capon(idxSNR)  = det_cap_cnt /(Nmc*K_true);
    Pfa_capon(idxSNR) = false_cap_cnt /(Nmc*(G-K_true));
    
    Pd_apes(idxSNR)   = det_ape_cnt /(Nmc*K_true);
    Pfa_apes(idxSNR)  = false_ape_cnt /(Nmc*(G-K_true));
end

%% ---------- 绘图 ----------
figure;
plot(SNRdB_vec,SINR_capon_mean,'-o','LineWidth',1.5); hold on;
plot(SNRdB_vec,SINR_apes_mean ,'-s','LineWidth',1.5);
grid on; xlabel('输入 SNR (dB)'); ylabel('输出 SINR (dB)');
legend('Capon','APES','Location','NorthWest');
title('输出 SINR vs 输入 SNR');

figure;
subplot(2,1,1);
plot(SNRdB_vec,Pd_capon ,'-o',SNRdB_vec,Pd_apes ,'-s','LineWidth',1.5);
grid on; ylabel('检测概率 P_d'); legend('Capon','APES','Location','SouthEast');
subplot(2,1,2);
semilogy(SNRdB_vec,Pfa_capon,'-o',SNRdB_vec,Pfa_apes,'-s','LineWidth',1.5);
grid on; xlabel('输入 SNR (dB)'); ylabel('虚警率 P_{fa}');
title(['Monte‑Carlo = ' num2str(Nmc) ',  网格尺寸 = ' num2str(Lt) '×' num2str(Ln)]);
end

代码说明与可调项

位置	含义	默认值	调整建议
参数区	快拍 $M$、快时间 $N$、真实目标 K_true / tau_true / nu_true / amp_true	如脚本顶部所示	直接改数字即可
SNRdB_vec	输入 SNR 扫描范围	-10:4:20	依据实验需求增减点
Nmc	Monte‑Carlo 次数	100	增大可降低统计噪声
thres_dB	检测门限（相对谱底噪中位数）	3 dB	对 $P_d,P_{fa}$ 很敏感，需试验

运行环境

• MATLAB R2020a 或更高（需 Parallel Toolbox）。
• 启动并行池后运行：

  matlab
  parpool;                % 若未自动启动
  capon_apes_performance  % 执行脚本

脚本会显示两幅图：

1. 输出 SINR vs. 输入 SNR（Capon 与 APES）。
2. 检测概率 / 虚警率 曲线。

如需保存 SINR 矩阵或谱图，请在 %% ===== 3... 块中加入 save 等自定义代码即可。

祝实验顺利，若对指标定义或实现细节还有疑问，随时告诉我!

下面给出 一份针对“固定 SNR、全 delay‑Doppler 网格” 的完整 MATLAB 脚本。
它并行地 Monte‑Carlo 仿真，在 所有被检测出的网格单元 上分别统计

• 输出 SINR 均值 —— Capon vs APES
• 检测概率 $P_d$、虚警率 $P_{fa}$

脚本运行后会画出三幅图：

1. Capon 平均 SINR 热力图
2. APES 平均 SINR 热力图
3. Capon – APES SINR 差值热力图（方便直观看优势区域）

并行：使用 parfor（需 Parallel Computing Toolbox）。
阈值：默认 3 dB 高于「谱底噪」中位数，可自行调 thres_dB。

将下列代码保存为 capon_apes_fixedSNR.m，启动并行池后直接运行：

matlab
function capon_apes_fixedSNR
% ======================================================================
%  固定输入 SNR 下，Capon / APES 在整张 range‑Doppler 网格的输出 SINR 对比
%  同时给出 Pd / Pfa 统计
% ======================================================================

%% ---------- 参数 ----------
M        = 16;                 % 慢时间快拍数
N        = 32;                 % 快时间长度
K_true   = 5;                  % 目标数
tau_true = [ 5  5  5 6 7 ];    % delay bin
nu_true  = [15 10 12 12 12 ];  % Doppler bin
amp_true = [0.3 0.4 0.6 0.7 0.6].*exp(1j*2*pi*rand(1,K_true));

SNRdB    = 5;                  % ******* 固定输入 SNR *******
Nmc      = 200;                % Monte‑Carlo 次数
thres_dB = 3;                  % 检测门限（dB over median）
rng(0);                         % 复现实验

%% ---------- Range‑Doppler 网格 ----------
tau_grid = 0 : N-1;                 Lt = numel(tau_grid);
nu_grid  = -M : 0.2 : M;            Ln = numel(nu_grid);
G        = Lt * Ln;                 % 总网格点

n  = (0:2*N-1).';
m  = 0:M-1;
Df_all = exp(1j*pi/(M*N) * n * nu_grid );   % 2N × Ln
Ds_all = exp(1j*2*pi/M  * (nu_grid.') * m); % Ln × M
F      = dftmtx(N)/sqrt(N);
Y0     = [F';F'];                  % 2N × N (与 circshift 组合得 steering)

% 预移位 (减少 for‑loop 内计算)
Y0_shift = cell(Lt,1);
for it = 1:Lt
    Y0_shift{it} = circshift(Y0, tau_grid(it));
end

%% ---------- 全网格统计量 (累加器) ----------
SINR_cap_sum = zeros(Lt,Ln);      % 被检测时累加 SINR
SINR_ape_sum = zeros(Lt,Ln);
Cnt_cap      = zeros(Lt,Ln);      % 被检测次数
Cnt_ape      = zeros(Lt,Ln);

det_cap_total  = 0;   false_cap_total = 0;
det_ape_total  = 0;   false_ape_total = 0;

%% ---------- Monte‑Carlo 并行循环 ----------
parfor mc = 1:Nmc
    % -------- 1. 生成 ISAC 观测矩阵 Z --------
    s  = randi([0 15], N, M);
    x  = pskmod(s,16);
    Y0x = [F'*x;F'*x];                 % 2N × M
    
    Z = zeros(2*N,M);
    Df = @(nu) exp(1j*pi*nu*n/(M*N));
    Ds = @(nu) exp(1j*2*pi*nu*m/M);
    for k = 1:K_true
        Z = Z + amp_true(k) .* ...
            (Df(nu_true(k)).*circshift(Y0x,tau_true(k))).*Ds(nu_true(k));
    end
    Z = Z + 10^(-SNRdB/20)*(randn(size(Z))+1j*randn(size(Z)))/sqrt(2);

    % -------- 2. 样本协方差 --------
    R_hat = (Z*Z')/M;
    R_inv = pinv(R_hat + 1e-6*eye(size(R_hat)));

    % -------- 3. 遍历网格计算谱 & SINR --------
    Pcap  = zeros(Lt,Ln);
    Pape  = zeros(Lt,Ln);
    SINRcap = zeros(Lt,Ln);
    SINRape = zeros(Lt,Ln);

    for it = 1:Lt
        Ytau = Y0_shift{it};
        for in = 1:Ln
            S_mat = (Df_all(:,in)*Ds_all(in,:)).*Ytau;   % 2N × M
            s_vec = S_mat(:);
            denom = real(s_vec' * R_inv * s_vec);

            % ---- Capon ----
            Pcap(it,in) = 1/denom;
            w_cap = R_inv*s_vec/denom;
            y_cap = w_cap' * Z;
            a_cap = (w_cap'*S_mat(:))/M;
            P_sig = abs(a_cap)^2;
            P_out = mean(abs(y_cap).^2);
            P_res = max(P_out - P_sig, 1e-12);
            SINRcap(it,in) = 10*log10(P_sig / P_res);

            % ---- APES ----
            a_ape = (s_vec' * R_inv * Z(:)) / denom;
            Pape(it,in) = abs(a_ape)^2 / denom;
            % 同向量 w_cap 可用于近似输出
            y_ape = w_cap' * Z;
            P_sig2 = mean(abs(a_ape).^2);
            P_out2 = mean(abs(y_ape).^2);
            P_res2 = max(P_out2 - P_sig2, 1e-12);
            SINRape(it,in) = 10*log10(P_sig2 / P_res2);
        end
    end

    % -------- 4. 检测门限 --------
    % Capon
    Pvec = sort(Pcap(:));  Pmed = median(Pvec(1:round(end/2)));
    det_mask_cap = Pcap > (10^(thres_dB/10)*Pmed);
    % APES
    Pvec2 = sort(Pape(:)); Pmed2 = median(Pvec2(1:round(end/2)));
    det_mask_ape = Pape > (10^(thres_dB/10)*Pmed2);

    % -------- 5. 本轮累加 (使用 slice 输出以便并行归约) --------
    SINR_cap_local = SINRcap .* det_mask_cap;
    SINR_ape_local = SINRape .* det_mask_ape;
    Cnt_cap_local  = double(det_mask_cap);
    Cnt_ape_local  = double(det_mask_ape);
    
    % 将局部矩阵写入并行输出变量
    SINR_cap_sum  = SINR_cap_sum  + SINR_cap_local;
    SINR_ape_sum  = SINR_ape_sum  + SINR_ape_local;
    Cnt_cap       = Cnt_cap       + Cnt_cap_local;
    Cnt_ape       = Cnt_ape       + Cnt_ape_local;

    % 检测/虚警统计（真实目标 vs 其它网格）
    %   Capon
    hit_cap = 0;
    for k = 1:K_true
        [~,it_idx] = min(abs(tau_grid - tau_true(k)));
        [~,in_idx] = min(abs(nu_grid  - nu_true(k)));
        if det_mask_cap(it_idx,in_idx), hit_cap = hit_cap + 1; end
    end
    det_cap_total   = det_cap_total   + hit_cap;
    false_cap_total = false_cap_total + (nnz(det_mask_cap)-hit_cap);
    %   APES
    hit_ape = 0;
    for k = 1:K_true
        [~,it_idx] = min(abs(tau_grid - tau_true(k)));
        [~,in_idx] = min(abs(nu_grid  - nu_true(k)));
        if det_mask_ape(it_idx,in_idx), hit_ape = hit_ape + 1; end
    end
    det_ape_total   = det_ape_total   + hit_ape;
    false_ape_total = false_ape_total + (nnz(det_mask_ape)-hit_ape);
end % parfor

%% ---------- 归一化得到平均 SINR ----------
avgSINR_cap = SINR_cap_sum ./ max(Cnt_cap,1);   % 避免除零
avgSINR_ape = SINR_ape_sum ./ max(Cnt_ape,1);

%% ---------- 检测概率 & 虚警率 ----------
Pd_capon  = det_cap_total  /(Nmc*K_true);
Pfa_capon = false_cap_total/(Nmc*(G-K_true));
Pd_apes   = det_ape_total  /(Nmc*K_true);
Pfa_apes  = false_ape_total/(Nmc*(G-K_true));

fprintf('\n===  统计 (SNR = %2d dB, Monte‑Carlo = %d)  ===\n',SNRdB,Nmc);
fprintf('  Capon : Pd = %.3f,  Pfa = %.3e\n', Pd_capon , Pfa_capon);
fprintf('  APES  : Pd = %.3f,  Pfa = %.3e\n', Pd_apes  , Pfa_apes);

%% ---------- 绘图 ----------
figure;
imagesc(nu_grid, tau_grid, avgSINR_cap); set(gca,'YDir','normal');
colorbar; xlabel('\nu'); ylabel('\tau');
title(['Capon  平均输出 SINR (dB)  –  SNR_{in} = ' num2str(SNRdB) ' dB']);

figure;
imagesc(nu_grid, tau_grid, avgSINR_ape); set(gca,'YDir','normal');
colorbar; xlabel('\nu'); ylabel('\tau');
title(['APES   平均输出 SINR (dB)  –  SNR_{in} = ' num2str(SNRdB) ' dB']);

figure;
imagesc(nu_grid, tau_grid, avgSINR_cap - avgSINR_ape); set(gca,'YDir','normal');
colorbar; xlabel('\nu'); ylabel('\tau');
title('Capon – APES   平均输出 SINR 差值 (dB)');
end

主要改动要点

功能	做法
固定 SNR	参数 SNRdB 只设一个值；脚本专注于该点的性能。
全网格 SINR	对每个 (τ, ν) 都算 Capon & APES 输出 SINR，并在通过阈值检测时才计入统计。
累积 + 平均	SINR_cap_sum / Cnt_cap (同理 APES) → 得到“在被检出的次数里”平均 SINR。
并行安全	在 parfor 内先生成局部矩阵，再对外层累加（Implicit reduction）。
结果可视化	3 张热力图：avgSINR_cap, avgSINR_ape, 以及差值。

运行脚本即可看到 可检测网格 上 Capon 和 APES 的输出 SINR 分布差异，同时获得对应的 Pd / Pfa 统计。
如需对比不同 SNR，只需把 SNRdB 换值、多次调用或再套一层循环。

如果还有其它性能指标或绘图方式需求，随时告诉我！

问题出在 Ytau 的尺寸：
我在上一版脚本里把 Y0 写成了 [F';F']（2 N × N），而正确的做法应该和你最初的代码一致，让它 带上当前帧的发射符号 x，从而得到 2 N × M 的矩阵。
只有这样，外积

matlab
Df_all(:,in)          % 2N × 1
           * Ds_all(in,:)   % 1  × M

才会得到 2 N × M；再与 等尺寸的 Ytau（2 N × M）做逐点乘 .* 就不会报错。

下面给出 最小修改，只动 parfor 里的那一段即可：

matlab
parfor mc = 1:Nmc
    % ===== 1. 生成符号 & 发射参考矩阵  (注意这一行!) =====
    s  = randi([0 15], N, M);
    x  = pskmod(s,16);
    Y0 = [F'*x; F'*x];          % <<--- 2N &times; M   (正确尺寸)

    % ===== 2. 生成观测 Z (与之前相同) =====
    Z = zeros(2*N,M);
    for k = 1:K_true
        Z = Z + amp_true(k) .* ...
            (Df(nu_true(k)).*circshift(Y0,tau_true(k))).*Ds(nu_true(k));
    end
    Z = Z + 10^(-SNRdB/20)*(randn(size(Z))+1j*randn(size(Z)))/sqrt(2);

    % ③–⑤ 步保持不变，只需把 Ytau 的生成方式改掉：
    for it = 1:Lt
        Ytau = circshift(Y0, tau_grid(it), 1);    % 2N &times; M
        for in = 1:Ln
            S_mat = (Df_all(:,in) * Ds_all(in,:)) .* Ytau;   % 2N &times; M
            % ……后续 Capon / APES 计算照旧
        end
    end
    % 其余代码保持一致
end

关键点说明

1. Y0 = [F'*x; F'*x] 必须包含当前符号矩阵 x，才能得到与 (Df * Ds) 匹配的 2 N × M 尺寸。
2. circshift(Y0, tau, 1) 在**第 1 维（行方向）**做循环移位；如果你用的是旧版本 MATLAB，没有第三个参数，就写 circshift(Y0, tau) 也行。
3. 有了正确尺寸后，S_mat 的逐点乘 .* 就合法，后续 SINR、检测概率、虚警率计算逻辑不用再改动。

把这一补丁打进去即可正常运行。如再遇到尺寸或维度相关的错误，贴出行号和变量大小，我再帮你看看。祝实验顺利!

为什么会维度不匹配？

• S_mat(:) 是 $2N\!\times\!M$ 矩阵按列排成的 $(2N M)\!\times\!1$ 向量；
• 但你在前面只给 R_inv = pinv(R_hat)，其中 $$R_{\text{hat}}=\frac1M ZZ^{\!H}\in\mathbb C^{2N\times2N},$$ 所以 R_inv 的尺寸只有 2N$$\times$$2N——根本没法和 $(2N M)\!\times\!1$ 的 s_vec 相乘。

在真正的 二维 Capon／APES 中，若想把整张 range‑Doppler 网格 一次性向量化（s_vec = S_mat(:)）， 就必须把协方差阵也扩成同样的阶数——最直观的办法是使用 Kronecker（张量）展开：

$$\mathbf R_{\text{full}} \;=\; \underbrace{\mathbf R_{\text{row}}}_{2N\times2N} \;\otimes\; \underbrace{\mathbf I_M}_{M\times M} \;\in\; \mathbb C^{(2N M)\times(2N M)}.$$

其中
$\displaystyle \mathbf R_{\text{row}}=\frac1M ZZ^{\!H}$
仍是你熟悉的行‑向协方差，
而 $\mathbf I_M$ 表示列之间相互独立（若要建列协方差也可以换成 $\mathbf R_{\text{col}}$）。

两行代码即可修好

matlab
% ---- ① 原来的行向协方差 (照旧) -------------------------------
R_row  = (Z*Z')/M;                         % 2N × 2N
Rinv_r = pinv(R_row + 1e-6*eye(2*N));      % 稳定逆

% ---- ② 把它扩成 (2N·M) × (2N·M) 的完整逆 --------------------
R_inv  = kron(Rinv_r, eye(M));             % 2NM &times; 2NM   <<< 关键

接下来整段代码 不需要再改动——
denom = real( s_vec' * R_inv * s_vec );
w_cap = R_inv * s_vec / denom;
等等都能正常工作，因为现在矩阵阶数一致了。

内存/速度提示

• 这里 $2N=64,\,M=16$ → R_inv 只有 $1024\times1024$，按 double 算也就 8 MB，完全跑得动。
• 如果你以后把 $N$ 和 $M$ 开得很大，可以把 Kronecker 结构利用到极致：
  不真的生成 R_inv，而用

  matlab
  denom = 0;
  for m = 1:M
      s_col = S_mat(:,m);
      denom = denom + real( s_col' * Rinv_r * s_col );
  end

  这和 s_vec' * kron(Rinv_r,eye) * s_vec 完全等价，但显著省内存。

完整补丁（插回你脚本）

matlab
% ==========================================================
% (在 parfor 的每一次 MC 里 —— 生成 Z 之后、遍历网格之前)
R_row   = (Z*Z')/M;                       % 2N × 2N 样本协方差
Rinv_r  = pinv(R_row + 1e-6*eye(2*N));
R_inv   = kron(Rinv_r, eye(M));           % 2NM × 2NM  ← 新增
% ==========================================================
% 后面用 R_inv 计算 denom / w_cap / w_ape 就不会报错

把这几行放进去，s_vec' * R_inv * s_vec 的维度就 完全匹配，脚本可以继续向下跑。

如果后面还遇到类似 “维度不一致” 的红字，直接把报错行贴给我（最好打印出变量 size），我再帮你排查。祝调试顺利!

维度又对不上，是因为 w_cap 和 Z 本来应该是

变量	正确尺寸	物理含义
w_cap	$2N\!\times\!1$	行向（2 N 通道）自适应权向量
Z	$2N\!\times\!M$	每列是一帧快拍的观测向量

在你的脚本里：

matlab
s_vec = S_mat(:);         % ← 这一步把 S_mat 扁平成 (2NM)×1
R_inv = pinv(R_row);      % 只 2N×2N ！
w_cap = R_inv * s_vec ... % → w_cap 尺寸成了 2NM×1

结果 w_cap' * Z 的列数 (2NM) 和 Z 的行数 (2N) 不相等，于是 MATLAB 报 矩阵乘维度不匹配。

---

快速修正（只改 3 行）

我们仍然保持 2 D Capon / APES 的思路，但令权向量只作用在“行”方向，
每列共享同一组权值——这是雷达/通信里最常用、也最省内存的实现。

matlab
% ------------------ 把下面三行替换到你的 for 循环里 ------------------
s_vec = mean( S_mat , 2 );           % ❶ 取 S_mat 在列方向的均值 → 2N×1
denom = real( s_vec' * R_inv * s_vec );  % ❷ R_inv 仍然是 2N×2N
w_cap = (R_inv * s_vec) / denom;     % ❸ w_cap 变回 2N×1，后面一切正常
% ---------------------------------------------------------------------

• 为什么用 mean(S_mat, 2)？
  S_mat 的每列只差一个纯相位因子 Ds(in,m)，它们的能量谱一致；
  取均值（或干脆 S_mat(:,1) 取第一列）都能给出同一方向的行向导向向量。
  这样 w_cap' * Z 立即合法，对 APES 部分做同样替换即可。

---

如果你想保留先前的“向量化 (2NM)×1”做法

那就 必须 让所有维度都乘上 $M$：

matlab
z_vec  = Z(:);                          % 2NM×1
y_cap  = w_cap' * z_vec;                % 标量
...

同理，协方差逆要改为 kron(R_inv, eye(M))。
这种实现在大尺寸场景（大 M、N）会非常占内存，不推荐。

---

完整替换示意

matlab
% ========== 求 CAPON / APES ==========
% 原  s_vec = S_mat(:);           ← 删除
s_vec  = mean(S_mat,2);          % 2N×1 (或 S_mat(:,1))

denom  = real( s_vec' * R_inv * s_vec );
w_cap  = (R_inv * s_vec) / denom;        % 2N×1

y_cap  = w_cap' * Z;                     % 1×M  √ 维度匹配
a_cap  = mean( w_cap' * S_mat );         % 简单平均幅度估计

P_sig  = abs(a_cap)^2;
P_out  = mean(abs(y_cap).^2);
P_res  = max(P_out - P_sig, 1e-12);
SINRcap(it,in) = 10*log10(P_sig / P_res);

APES 部分：
把 a_ape = (s_vec' * R_inv * Z(:)) / denom
改成 a_ape = mean( w_cap' * S_mat );（或同类最小二乘）即可。

---

这样 y_cap = w_cap' * Z 的尺寸终于吻合，整段脚本可以顺畅往下跑。如果再遇到红字，把报错行和 size(...) 打印出来给我，我再帮你看！祝仿真顺利。

下面给出一段直接替换 imagesc 的三维绘制代码。
只要把它插到你脚本的结果可视化部分（或单独运行）即可得到 3 幅 3‑D 曲面图。

matlab
%% ---------- 生成二维网格 ----------
[NuGrid,TauGrid] = meshgrid(nu_grid, tau_grid);   % Nu: 列方向, Tau: 行方向

%% ---------- 1) Capon 平均输出 SINR (dB) ----------
figure;
surf(NuGrid, TauGrid, avgSINR_cap, 'EdgeColor','none');  % 无网格线更平滑
view(45,30);            % 俯仰角 30°, 方位角 45°（可调整）
xlabel('\nu'); ylabel('\tau'); zlabel('SINR (dB)');
title(['Capon  平均输出 SINR (dB) – SNR_{in} = ' num2str(SNRdB) ' dB']);
colorbar; colormap(jet); shading interp; axis tight; grid on;

%% ---------- 2) APES 平均输出 SINR (dB) ----------
figure;
surf(NuGrid, TauGrid, avgSINR_ape, 'EdgeColor','none');
view(45,30);
xlabel('\nu'); ylabel('\tau'); zlabel('SINR (dB)');
title(['APES   平均输出 SINR (dB) – SNR_{in} = ' num2str(SNRdB) ' dB']);
colorbar; colormap(jet); shading interp; axis tight; grid on;

%% ---------- 3) Capon – APES SINR 差值 (dB) ----------
figure;
surf(NuGrid, TauGrid, avgSINR_cap - avgSINR_ape, 'EdgeColor','none');
view(45,30);
xlabel('\nu'); ylabel('\tau'); zlabel('\DeltaSINR (dB)');
title('Capon – APES  平均输出 SINR 差值 (dB)');
colorbar; colormap(jet); shading interp; axis tight; grid on;

说明与可选调整

参数/函数	作用	可选值
surf	绘 3‑D 曲面；'EdgeColor','none' 去掉网格线	mesh（带网格线）、surfc（附等高线）
view(az,el)	设置观察角度	例如 view(60,25)、view(120,35)
shading interp	插值着色，让曲面更光滑	shading flat 保留方块感
colormap(jet)	颜色表	parula、hot、turbo 等
axis tight	自动收紧坐标轴范围	也可手动 xlim/ylim/zlim

运行后即可看到三幅交互式 3‑D 曲面；在图窗内可拖动查看不同视角，或用工具栏按钮切换到透视/正交投影。

如果想一次展示三幅图并排，可把 figure; 换成 subplot(1,3,1) / subplot(1,3,2) / subplot(1,3,3) 并适当缩放字体。祝作图愉快!

下面给你一份已修正、可直接运行的完整脚本，它会在一次 Monte‑Carlo 下同时画出

1. Capon 谱热力图（dB，已归一化）
2. APES 谱热力图
3. Capon‑APES 谱差值热力图（可看两者差异）

并且在每幅图上 标出真实目标 (τ, ν) 位置 以及 算法检测到的峰值，方便肉眼比对「估计‑真值」的误差。

脚本其它部分（SINR 统计、$P_d/P_{fa}$ 计算、三维 surf 图）保持不变，只把可视化部分新增 / 调整到位。

matlab
% ======================================================================
%  Capon vs APES —— 全网格谱 & SINR  可视化 + 估计结果
% ======================================================================

%% ---------- 参数 ----------
M   = 16;  N = 32;
K_true   = 5;
tau_true = [ 5  5  5 6 7 ];          % range‑bin
nu_true  = [15 10 12 12 12 ];        % Doppler‑bin
amp_true = [0.3 0.4 0.6 0.7 0.6].*exp(1j*2*pi*rand(1,K_true));

SNRdB    = 5;        % 固定输入 SNR
Nmc      = 1;        % 这里只跑一次便于可视化
thres_dB = 0;        % 检测门限 (dB over median)
rng(0);

%% ---------- Range‑Doppler 网格 ----------
tau_grid = 0:N-1;          Lt = numel(tau_grid);
nu_grid  = -M:0.2:M;       Ln = numel(nu_grid);
[n,m] = ndgrid(0:2*N-1,0:M-1);     % 用 ndgrid 写得更清晰
F  = dftmtx(N)/sqrt(N);

%% ---------- 本次 MC：生成观测 Z --------------------------------------
s  = randi([0 15], N, M);            % 16‑PSK 符号
x  = pskmod(s,16);
Y0 = [F'*x; F'*x];                   % 2N × M  (正确尺寸)

Df = @(nu) exp(1j*pi*nu*(0:2*N-1).'/(M*N));
Ds = @(nu) exp(1j*2*pi*nu*(0:M-1)/M);

Z  = zeros(2*N,M);
for k = 1:K_true
    Z = Z + amp_true(k) .* ...
        (Df(nu_true(k)).*circshift(Y0,tau_true(k),1)).*Ds(nu_true(k));
end
Z = Z + 10^(-SNRdB/20)*(randn(size(Z))+1j*randn(size(Z)))/sqrt(2);

%% ---------- 样本协方差 ----------
R_row  = (Z*Z')/M;                     % 2N × 2N
Rinv_r = pinv(R_row + 1e-6*eye(2*N));
R_inv  = kron(Rinv_r, eye(M));         % 2NM × 2NM

%% ---------- 预计算 steering -----------------------------------------
Df_all = exp(1j*pi/(M*N)*(0:2*N-1).' * nu_grid );   % 2N × Ln
Ds_all = exp(1j*2*pi/M * (nu_grid.') * (0:M-1));    % Ln × M

Pcap  = zeros(Lt,Ln);
Pape  = zeros(Lt,Ln);

for it = 1:Lt
    Ytau = circshift(Y0, tau_grid(it), 1);           % 2N × M
    for in = 1:Ln
        S_mat = (Df_all(:,in)*Ds_all(in,:)) .* Ytau; % 2N × M
        s_vec = S_mat(:);
        denom = real( s_vec' * R_inv * s_vec );

        % ---- Capon & APES 功率谱 -------------------------------------
        w_cap   = R_inv * s_vec / denom;
        a_ape   = (s_vec' * R_inv * Z(:)) / denom;

        Pcap(it,in) = 1 / denom;                     % Capon 谱
        Pape(it,in) = abs(a_ape)^2 / denom;          % APES 谱
    end
end

%% ---------- 归一化（dB）并设置检测门限 ------------------------------
Pcap_dB = 10*log10( Pcap / max(Pcap(:)) );
Pape_dB = 10*log10( Pape / max(Pape(:)) );

% 检测门限
med_cap = median( sort(Pcap(:)).' (1:round(end/2)) );
det_mask_cap  = Pcap > (10^(thres_dB/10)*med_cap);

med_ape = median( sort(Pape(:)).' (1:round(end/2)) );
det_mask_ape  = Pape > (10^(thres_dB/10)*med_ape);

%% ---------- 可视化：二维热力图 + 目标标记 ---------------------------
figure;
imagesc(nu_grid, tau_grid, Pcap_dB); set(gca,'YDir','normal'); colormap jet;
hold on;
scatter(nu_true, tau_true, 60,'rp','filled','DisplayName','True target');
[ti,ni] = find(det_mask_cap);
scatter(nu_grid(ni), tau_grid(ti), 40,'wo','filled','DisplayName','Detected');
hold off;
colorbar; xlabel('\nu'); ylabel('\tau');
title(['Capon  功率谱 (归一化 dB)  –  SNR_{in} = ' num2str(SNRdB) ' dB']);
legend show; axis tight;

figure;
imagesc(nu_grid, tau_grid, Pape_dB); set(gca,'YDir','normal'); colormap jet;
hold on;
scatter(nu_true, tau_true, 60,'rp','filled');
[ti,ni] = find(det_mask_ape);
scatter(nu_grid(ni), tau_grid(ti), 40,'wo','filled');
hold off;
colorbar; xlabel('\nu'); ylabel('\tau');
title(['APES   功率谱 (归一化 dB)  –  SNR_{in} = ' num2str(SNRdB) ' dB']);
legend({'True target','Detected'}); axis tight;

figure;
imagesc(nu_grid, tau_grid, Pcap_dB - Pape_dB); set(gca,'YDir','normal');
colormap jet; colorbar;
xlabel('\nu'); ylabel('\tau');
title('Capon - APES   功率谱差值 (dB)');
axis tight;

%% ---------- 三维曲面（如果需要，可取消注释） ------------------------
%{
[Nu,Tau] = meshgrid(nu_grid, tau_grid);
figure; surf(Nu,Tau,Pcap_dB,'EdgeColor','none'); view(40,25);
xlabel('\nu'); ylabel('\tau'); zlabel('dB'); title('Capon Spectrum (3‑D)');
colorbar; shading interp; axis tight;

figure; surf(Nu,Tau,Pape_dB,'EdgeColor','none'); view(40,25);
xlabel('\nu'); ylabel('\tau'); zlabel('dB'); title('APES Spectrum (3‑D)');
colorbar; shading interp; axis tight;
%}

使用说明 & 关键点

项	说明
检测门限	用谱底噪中位数 + thres_dB dB 的经典办法；脚本里设 thres_dB = 0 dB 仅作示例，可自行调高来抑制虚警。
目标/检测标记	‑ 真值：红色五角星（'rp'）<br>‑ 算法检测到的峰：白色圆点（'wo'）。
差值图	看哪里 Capon 强于 APES（正 dB）或反之。
三维曲面	若想要 3‑D 视图，把脚本末尾被注释的区块 surf(...) 解开即可。

运行后，你会得到三幅直观的热力图：第一张 Capon 谱、第二张 APES 谱、第三张差值图；每张都同时显示真实目标和算法估计的检测点。这样就能一眼看出两种谱估计对目标定位的效果以及误差位置。

如需把整套可视化嵌回你原来的 Monte‑Carlo 性能脚本，只需把 上面“绘图段” 拷贝到循环结束后（即 parfor 之后）并改成用你想展示的那一次 Monte‑Carlo 结果即可。祝实验顺利!