logoChat01.ai

你现在是无线通信的专家。我们现在一起来分析一个问题，就是OFDM在分析AF时候的性能。我们都知道IC...

創建於：2026年1月21日

使用 GPT-5.2 Thinking 回答 Chat01

問題

你现在是无线通信的专家。我们现在一起来分析一个问题，就是OFDM在分析AF时候的性能。我们都知道ICI对OFDM造成的影响很大，原因是快时间doppler的影响。我们下面的代码是分别分析相干积分前后的2D-FFT的AF和逐点匹配的AF。我们现在想对于相干积分前的AF，直接对每一个doppler的点进行快时间部分的多普勒补偿，来消除这个干扰。请在下面这个代码里面实现，并与补偿前的进行对比，我们要对比的就是相干积分前的AFfs进行补偿前后的对比：clc
clear

format long

times = 1;
M = 16; % slow time index 脉冲数量 16
N = 32*times; % fast time index 每个脉冲的采样点数 32
Nmc = 10;
rng(1545);
tau = 4;
range1 = -M:M;
rec = zeros(length(range1),1);
count = 0;

for nu = 8
count = count + 1;
cAF = zeros(N,MN2); % 初始化常规模糊函数 cAF
for mc = 1:Nmc

text
    %OFDM gen
    s = randi(16,[N,M])-1; 
    x = pskmod(s,16); 
    F = dftmtx(N)./sqrt(N);
    y = repmat(F'*x,[2,1]);
    
    % delay doppler OFDM
    yy1 = y(:); 
    yy2 = circshift(yy1,tau); 
    yy3 = yy2.*exp(1i*2*pi*(0:M*N*2-1)/M/N/2*nu).'; 
    
    % transmit OFDM
    xx = repmat(F'*x,[2,1]);
    xx = xx(:); 

    y2 = reshape(yy3,N*2,M); 
    z = y2(:); 

    AF = zeros(N,M*N*2); % 初始化常规模糊函数 AF
    parfor i = 1:N
            row_AF = zeros(1, M*N*2);  
            for j = 1:M*N*2
                temp1 = z .* exp(-1i*2*pi*(0:M*N*2-1)'/(M*N*2)*j);
                temp2 = temp1 .* circshift(conj(xx), i);
                temp3 = reshape(temp2, N*2, M);
                temp4 = temp3(N+1:end, :);
                row_AF(j) = sum(temp4(:));
            end
            AF(i, :) = row_AF;  
    end
    cAF = cAF + abs(AF).^2; 
end

AFfs = zeros(N,M); % 初始化修改后的快慢时间模糊函数
F2 = dftmtx(M)./sqrt(dftmtx(M)); 
AFfs = abs(fft((F'*((F*y2(N+1:end,:)).*conj(x))).')).^2; % 快慢时间模糊函数的功率谱计算

 %figure;
 %surf(ifftshift(20*log10(AFfs)));

 % PSR calc
if nu~=0
    rec(count) = 10*log10(max(max(AFfs))./ (sum(sum(AFfs))-max(max(AFfs))).*(M*N+1)); 
else
    [~, idx] = max(AFfs(:));
    AFfs_filtered = AFfs;
    AFfs_filtered(idx) = 0;
    rec(count) = 10*log10(max(AFfs, [], 'all') / sum(AFfs_filtered, 'all') * (M*N+1));
end

end

figure;
plot(range1,rec-max(rec),'o-');
xlabel('Doppler Shift Index')
ylim([-350,0])
title('N = 16*32')
grid on

figure;
plot(range1,rec,'o-');
xlabel('Doppler Shift Index')
title('N = 16*32')
grid on

%% 非相干累积
AF3 = zeros(N,N*2);
for i = 1:M
AF3 = AF3 + cAF(:,i:M:end);
end
AF3 = AF3 / M;

figure
surf(ifftshift(10*log10(AF3./max(max(AF3)))));
shading interp;
zlabel('Ambiguity Function/dB')
xlabel('Doppler Index');
ylabel('Delay Index');
title('After Non-Coherent Intergration')

% 绘图
figure;
imagesc(ifftshift(10log10(AF3./max(max(AF3)))));
colormap jet;
colorbar;
clim([-30 0]);
xlabel('Doppler Index (Cyclic Shifted)');
ylabel('Delay Index');
title('After Non-Coherent Intergration')
axis xy;
set(gca, 'XTick', [1, N, 2N]);
set(gca, 'XTickLabel', {'-N', '0', 'N'});
set(gca, 'YTick', [1, N/2, N]);
set(gca, 'YTickLabel', {'0', 'N/2', 'N'});

%% 未进行非相干累积
AF2 = zeros(N, N2); % 修正为正确的行维度
for i = 1:M
AF2 = AF2 + abs(AF(:, N2*(i-1)+1:N2i)).^2; % 逐块累加
end
AF2 = AF2 ./ M; % 平均操作

figure
surf(ifftshift(10*log10(AF2./max(max(AF2)))));
shading interp;
zlabel('Ambiguity Function/dB')
xlabel('Doppler Index');
ylabel('Delay Index');
title('Before Non-Coherent Intergration')

figure;
imagesc(ifftshift(10log10(AF2./max(max(AF2)))));
colormap jet;
colorbar;
clim([-30 0]);
xlabel('Doppler Index (Cyclic Shifted)');
ylabel('Delay Index');
title('Before Non-Coherent Intergration')
axis xy;
set(gca, 'XTick', [1, N, 2N]);
set(gca, 'XTickLabel', {'-N', '0', 'N'});
set(gca, 'YTick', [1, N/2, N]);
set(gca, 'YTickLabel', {'0', 'N/2', 'N'});

%% zero delay and doppler cut

% zero Doppler
zero_doppler_AF2 = AF2(:, 1);
zero_doppler_AF3 = AF3(:, 1);

figure;
plot(10log10(abs(zero_doppler_AF2./max(zero_doppler_AF2))), 'r', 'LineWidth', 1.5);
hold on;
plot(10log10(abs(zero_doppler_AF3./max(zero_doppler_AF3))), 'b', 'LineWidth', 1.5);
xlabel('Delay Index');
legend('AF2', 'AF3');
title('Zero Doppler Cut')
grid on

% zero delay
doppler_index = -(N-1):N;
zero_delay_AF2 = AF2(tau, :);
zero_delay_AF3 = AF3(tau, :);
figure;
plot(doppler_index,ifftshift(10log10(abs(zero_delay_AF2./max(zero_delay_AF2)))), 'r', 'LineWidth', 1.5);
hold on;
plot(doppler_index,ifftshift(10log10(abs(zero_delay_AF3./max(zero_delay_AF3)))), 'b', 'LineWidth', 1.5);
xlabel('Doppler Index');
legend('AF2', 'AF3');
title('Zero Delay Cut')
grid on

思考

Modifying MATLAB for Doppler compensation

I need to adjust the MATLAB code to implement fast-time Doppler compensation for each Doppler bin before coherent integration. Given the formula used for AFfs, I should apply frequency correction for each Doppler shift. Specifically, this involves multiplying by exp(-j2πn/N * fd_hat) for each Doppler point and verifying improvements to AFfs.

Parsing the AFfs formula for Doppler compensation

The term for Doppler phase increment is clarified as π*nu/(MN). It affects both fast and slow time. In the AFfs coherent method, I transform to the frequency domain for the fast time, match with x, and then perform an FFT across the slow time. The AFfs formula involves FFTs and matched filtering, then a 2D FFT of the output.

Compensating fast-time Doppler per bin

To compensate for fast-time Doppler per Doppler bin, I'll apply a phase ramp for each pulse before the FFT, with the focus on the fast-time component. We'll hypothesize Doppler, de-rotate the time-domain samples, and multiply by exp(-j2π fD n) for each pulse to correct the fast-time shift, leaving slow-time Doppler handled by the slow-time FFT.

Compensating Doppler for entire AFfs

To compute AFfs for each Doppler bin after slow-time FFT, I need to apply compensation for each bin based on fD for that bin. This means applying a fast-time derotation with fD mapped across both fast and slow time. The mapping uses fD = (k - floor(M/2)) / (MN2), adjusting for the full signal length L = MN2.

Deriving Doppler compensation for AFfs

To apply Doppler compensation for each Doppler bin, I need to correct the phase shift for each pulse and match the slow-time FFT bins. The correct Doppler frequency fD for a given bin k is derived as k/(MN2). This matches the physical Doppler frequency by considering the phase change across both fast and slow time.

Implementing Doppler compensation in AFfs

To implement fast-time compensation for each Doppler bin k, I’ll multiply the fast-time samples by exp(-j2π fD n) before performing the FFT. For each bin, calculate fD based on the Doppler index k, where fD = (k-1)/(MN2). Then, I can compute the compensated and raw AFfs and compare them.

Conducting Doppler compensation and range processing

For each Doppler bin k, I’ll apply the Doppler compensation to all pulses in the y2 matrix by multiplying it with a phase shift. Then, I’ll compute the range output after applying the matched filter via FFT. The computed range matrix is then used for Doppler extraction using the FFT, after flipping dimensions as needed. Ultimately, the result is assigned into the AFfs matrix.

Fixing dimensions in FFT processing

It seems like there’s confusion in the dimensions for AFfs. The code uses FFT without dimension specification, which leads to dimension mismatches. I’ll adjust the FFT to apply across the correct dimension (slow-time) along columns. This should correct the transposition error and ensure minimal changes to the existing structure.

Implementing code to fix dimension mismatch and doppler compensation

I need to ensure that the user can keep their original formula while adding the corrected version. It looks like there might be a dimension mismatch with their current code, so we should fix that first. I'll provide an updated line for the range matrix and FFT, and then implement Doppler compensation. The code will adjust each Doppler bin (k), applying the compensation to each pulse.

Doppler compensation and FFT processing

For each Doppler bin, I need to apply compensation based on the corresponding Doppler frequency (fd). After we derotate with that fd, the FFT should focus on the residual phase at the corresponding bin. The method would involve FFT over range while keeping the Doppler match from fast-time derotation. Then, we’ll correlate across pulses.

Doppler compensation and AFfs comparison

To compute Doppler output for each bin, I apply Doppler compensation based on the corresponding frequency shift for that bin. After compensating for this, I can compute the Doppler output using standard FFT. I'll compare the results before and after compensation by normalizing the outputs and visualizing them using imagesc to highlight any changes. Additionally, I'll calculate the PSR for both raw and compensated results.

Adjusting AFfs calculation for MC loop and comparison

It seems they compute AFfs outside the MC loop, which could be the issue. To match the desired comparison, I plan to compute AFfs for each MC run and average it over all loops. I'll accumulate AFfs_raw and AFfs_comp, then show modified MATLAB code that includes this averaging approach. Then, I’ll plot both the raw and compensated AFfs for comparison.

Thought for 2m 30s

回答

下面给你一个**“逐 Doppler 点做快时间 Doppler 补偿”的实现方式，并且把补偿前/补偿后**的 **AFfs（相干积分前 2D-FFT AF）**直接对比画出来。

核心想法是：

  • 你现在的 AFfs 做法等价于：先做快时间 FFT（子载波域），再乘 conj(x) 匹配滤波，再做快时间 IFFT 得到延迟像，最后沿慢时间做 FFT 得到 Doppler。
  • 但 Doppler 相位是沿整个快时间采样点变化的，所以会在快时间 FFT 时产生 ICI。
  • 因此对每一个 Doppler bin（慢时间 FFT 的第 k 个点），先对每个脉冲内部的快时间序列乘上一个反向旋转： e−j2πfd(k) ne^{-j2\pi f_d(k)\,n} 其中 nn 是快时间采样索引，fd(k)f_d(k) 是该 Doppler bin 对应的“每采样的归一化频率”（cycles/sample）。
  • 注意：我们只补偿“快时间内的相位滚动”，不补偿跨脉冲的相位项（那一项本来就由慢时间 FFT 来处理）。

1) 直接可用的代码修改（在你代码里替换/插入）

你原来 AFfs = ... 那一段有维度不一致的隐患（.' 之后是 M×N，但你定义了 AFfs=zeros(N,M)）。我这里给出一个更明确的写法：沿慢时间维（第2维）做 FFT，得到 N×M 的 AFfs。

把下面这段放到你每次 mc 里、得到 y2 之后（即 y2 = reshape(yy3,N*2,M); 之后），并把你原来的 AFfs = ... 替换掉即可。

matlab
% ====== 在 mc 循环内：y2 已经得到之后加入/替换 ======
y_fast = y2(N+1:end, :);              % 取后半段: N×M

% --- (A) 补偿前：原始 AFfs（2D-FFT 形式，先快时间处理，再慢时间FFT） ---
range_mat = F' * ((F*y_fast) .* conj(x));     % N×M
AFfs_raw_mc = abs( fft(range_mat, M, 2) ).^2; % 沿慢时间(第2维)做FFT => N×M


% --- (B) 逐 Doppler bin 的快时间 Doppler 补偿 AFfs ---
AFfs_comp_mc = zeros(N, M);

n_fast = (0:N-1).';                   % 快时间采样索引（每个脉冲内重置）
p_idx  = (0:M-1).';                   % 慢时间脉冲索引

for k = 1:M
    % 把 k 映射到“居中”的 Doppler bin（便于同时覆盖正负频率）
    k0 = k-1;
    if k0 > M/2
        k0 = k0 - M;
    end

    % 该 Doppler bin 对应的“每采样”归一化频率 cycles/sample
    % 推导：慢时间FFT bin 满足 exp(j2π f_d * (N*2) * p) <-> exp(j2π k p/M)
    % => f_d = k/(M*N*2)
    fd = k0 / (M*N*2);

    % 只做“快时间内”的去旋转（消ICI）
    rot_fast = exp(-1i*2*pi*fd*n_fast);        % N×1
    y_fast_comp = y_fast .* rot_fast;          % N×M (隐式扩展)

    % 快时间处理（子载波域匹配 + IFFT 到 delay）
    range_mat_comp = F' * ((F*y_fast_comp) .* conj(x));  % N×M

    % 该 Doppler 点的慢时间相干和（等价于取慢时间FFT的第k个bin）
    steering = exp(-1i*2*pi*p_idx*(k-1)/M);    % M×1
    AFfs_comp_mc(:, k) = abs(range_mat_comp * steering).^2; % N×1
end

2) 为了对比：在 mc 外做累加平均 + 画图

你原代码里 AFfs 没有在 Nmc 上平均（cAF 平均了，但 AFfs 没有）。为了公平对比补偿前后，建议你在 nu 循环里加两行累加器：

在 for nu = 8 内，for mc = 1:Nmc 前面加：

matlab
AFfs_raw_acc  = zeros(N, M);
AFfs_comp_acc = zeros(N, M);

然后在每次 mc 结束时累加：

matlab
AFfs_raw_acc  = AFfs_raw_acc  + AFfs_raw_mc;
AFfs_comp_acc = AFfs_comp_acc + AFfs_comp_mc;

mc 循环结束后得到平均：

matlab
AFfs_raw  = AFfs_raw_acc  / Nmc;
AFfs_comp = AFfs_comp_acc / Nmc;

3) 补偿前后 AFfs 的对比绘图（相干积分前）

把下面画图代码接在 AFfs_raw / AFfs_comp 得到之后即可：

matlab
% 归一化到各自峰值，转 dB
AFfs_raw_dB  = 10*log10(AFfs_raw  ./ max(AFfs_raw(:)));
AFfs_comp_dB = 10*log10(AFfs_comp ./ max(AFfs_comp(:)));

figure;
imagesc(ifftshift(AFfs_raw_dB, 2));  % 只对Doppler维做shift（第2维）
axis xy; colormap jet; colorbar;
clim([-30 0]);
xlabel('Doppler Bin (shifted)');
ylabel('Delay Bin');
title('AFfs BEFORE fast-time Doppler compensation (coherent, pre-NCI)');

figure;
imagesc(ifftshift(AFfs_comp_dB, 2));
axis xy; colormap jet; colorbar;
clim([-30 0]);
xlabel('Doppler Bin (shifted)');
ylabel('Delay Bin');
title('AFfs AFTER fast-time Doppler compensation (per-Doppler, coherent, pre-NCI)');

figure;
imagesc(ifftshift(AFfs_comp_dB - AFfs_raw_dB, 2));
axis xy; colormap jet; colorbar;
xlabel('Doppler Bin (shifted)');
ylabel('Delay Bin');
title('AFfs Improvement (AFfs\_comp - AFfs\_raw) in dB');

4) 可选：用一个数值指标比较（PSR）

如果你也想像原来那样用 PSR 比较一眼看提升，可以加：

matlab
peak_raw = max(AFfs_raw(:));
psr_raw = 10*log10( peak_raw / (sum(AFfs_raw(:)) - peak_raw) * (M*N+1) );

peak_comp = max(AFfs_comp(:));
psr_comp = 10*log10( peak_comp / (sum(AFfs_comp(:)) - peak_comp) * (M*N+1) );

fprintf('PSR raw  = %.3f dB\n', psr_raw);
fprintf('PSR comp = %.3f dB\n', psr_comp);
fprintf('Gain     = %.3f dB\n', psr_comp - psr_raw);

你会看到什么现象（通常）

  • AFfs_raw：主峰会被 ICI“糊开”，旁瓣/杂散抬高。
  • AFfs_comp：主峰更集中，杂散显著下降（尤其是 Doppler 不小的时候），PSR 上升。

如果你接下来想更“严格”地对齐你逐点匹配的 AF（你 parfor 那个），我们还可以把 fd 的映射和 Doppler 轴定义完全统一（比如和你 exp(-j2π n/(M*N*2)*j) 的 j 对齐），这样补偿会更精确。

問題

给我一个完整的，对比这两个方案的代码

思考

Providing MATLAB code for compensation and comparison

The user asked for MATLAB code to calculate and compare two AFfs schemes: one with compensation before coherent integration and one after. Their code computes AFfs directly, but not averaged. I'll modify it to compute averaged AFfs per Doppler bin, integrate it within the loop, and ensure performance is acceptable while preserving the original code structure, fixing potential errors along the way.

Clarifying AFfs computation and compensation logic

In MATLAB, dftmtx provides an unnormalized DFT matrix, and multiplying by F'*x gives the IDFT. For AFfs, I compute Yk as F*y_fast, then multiply by conj(x), and use F'* to get the delay profile. For compensation, each Doppler bin has a frequency fd=k0M∗N∗2f_d = \frac{k_0}{M*N*2}, affecting y_fast via a phase rotation across the fast-time index.

Understanding slow-time and Doppler bin mapping

In slow-time FFT, the index p corresponds to blocks of N2 samples, and the frequency per bin is calculated as k/(MN2). This works out for Doppler bins, where the phase advances for each pulse by exp(j2π fd N2). For AFfs, we map Doppler indices between -M/2 and M/2-1 using ifftshift. The computational cost is manageable with N=32 and averaging across Nmc=10, making looping feasible.

Generating the full script with AFfs and comparison

The user wants a full script comparing the two schemes for compensation, both before and after coherent integration, and asked for plots and a PSR printout. I'll include both schemes and ensure the original AF2 and AF3 code remains intact, while avoiding overwriting from the parfor loop. The results will be averaged over pulses but not across multiple MCs, which might be fine for their purposes.

Refining AF2, AF3, and PSR calculation logic

I’ll fix the inconsistency between AF2 and AF3 by using cAF consistently. For AF2, I'll modify the code to use cAF for averaging across both MCs and M. This addresses possible confusion from using the last MC's AF. For PSR, I'll keep the user’s setup but may offer an option to scan across the full Doppler range for comparison. I'll focus on the requested compensation comparison but also provide insights with potential adjustments for PSR.

Including PSR and nu sweep for comparison

For a more robust comparison, I’ll compute the PSR curves for both raw and compensated AFfs, and plot them against the range. It might increase computation time, but it's manageable with small N and M sizes. I'll keep the nu sweep optional, defaulting to 8, and provide the full MATLAB code with both methods for the user’s reference.

Reorganizing code structure for clarity

I’ll structure the code with parameters like times, M, N, Nmc, rng, and tau. I’ll define F using a DFT matrix and set nu_list to range1. During each loop over nu, I’ll compute AFfs_raw and AFfs_comp, then calculate the PSR for each scheme. Afterward, I'll plot the 2D maps for specific nu values and include AF2 and AF3 calculations.

Addressing MATLAB code adjustments and clarifications

I’ll update the code by adjusting averages for consistency (dividing by M*Nmc rather than just M) and annotating these changes. The code also needs certain adjustments for MATLAB, like ensuring compatibility with parfor, matrix multiplication, and identifying ifftshift and temp1 definitions. I'll add a new plot for AFfs improvement and preserve consistent scaling and output formats.

Thought for 1m 21s

回答

下面是一份完整可运行的 MATLAB 脚本：同时计算并对比

  1. 补偿前的相干（2D-FFT）AFfs_raw
  2. 逐 Doppler 点做快时间 Doppler 补偿后的 AFfs_comp

并给出：

  • AFfs_raw / AFfs_comp 的 2D 图对比（以 nu_show 为例）
  • PSR 随 nu 变化的曲线对比（两个方案）

你只要把下面整段脚本保存为一个 .m 文件运行即可。

matlab
clc; clear;
format long

%% ===================== 参数区 =====================
times = 1;
M = 16;                 % slow time pulses
N = 32*times;           % fast time samples per pulse (useful part)
Nmc = 10;               % Monte Carlo runs
rng(1545);
tau = 4;                % delay (samples)
range1 = -M:M;          % Doppler index scan (for PSR curve)

nu_list = range1;       % 你也可以改成 nu_list = 8; 只看一个 Doppler
nu_show = 8;            % 用于画 AFfs 2D 对比图的 Doppler（必须在 nu_list 里）

F = dftmtx(N)./sqrt(N); % unitary DFT matrix

% 记录 PSR
rec_raw  = zeros(length(nu_list),1);
rec_comp = zeros(length(nu_list),1);

% 保存 nu_show 下的结果用于画图
AFfs_raw_show  = [];
AFfs_comp_show = [];
cAF_show = [];          % 也保存逐点匹配的 cAF（用于你原来的 AF2/AF3 部分对照）

%% ===================== 主循环：扫 Doppler =====================
for idx_nu = 1:length(nu_list)
    nu = nu_list(idx_nu);

    % 逐点匹配 AF 的非相干累积容器（和你原来一致）
    cAF = zeros(N, M*N*2);

    % 2D-FFT AFfs 的 MC 平均容器（新增：raw & comp）
    AFfs_raw_acc  = zeros(N, M);
    AFfs_comp_acc = zeros(N, M);

    for mc = 1:Nmc
        %% -------- OFDM gen --------
        s = randi(16,[N,M]) - 1;
        x = pskmod(s,16);                 % N×M (subcarrier symbols)
        y = repmat(F'*x,[2,1]);           % time-domain (N×M)->(2N×M), repeat twice

        %% -------- delay + Doppler (fast-time Doppler causes ICI) --------
        yy1 = y(:);
        yy2 = circshift(yy1, tau);

        L = M*N*2; % total samples
        % Doppler across global sample index
        yy3 = yy2 .* exp(1i*2*pi*(0:L-1)'/L * nu);

        %% -------- transmit reference --------
        xx = repmat(F'*x,[2,1]);
        xx = xx(:);

        %% -------- reshape received --------
        y2 = reshape(yy3, N*2, M);        % (2N)×M
        z  = y2(:);

        %% ============================================================
        % (1) 逐点匹配的 AF（和你原来一致）-> 累积到 cAF
        % ============================================================
        AF = zeros(N, M*N*2);
        parfor i = 1:N
            row_AF = zeros(1, M*N*2);
            for j = 1:M*N*2
                temp1 = z .* exp(-1i*2*pi*(0:L-1)'/L * j);
                temp2 = temp1 .* circshift(conj(xx), i);
                temp3 = reshape(temp2, N*2, M);
                temp4 = temp3(N+1:end, :);
                row_AF(j) = sum(temp4(:));
            end
            AF(i,:) = row_AF;
        end
        cAF = cAF + abs(AF).^2;

        %% ============================================================
        % (2) 相干(2D-FFT) AFfs：补偿前 & 逐 Doppler 点快时间补偿后
        % ============================================================

        % 取后半段（对应你原来 temp4 的那段）
        y_fast = y2(N+1:end, :);          % N×M

        % ---------- (A) 补偿前：原始 AFfs_raw ----------
        range_mat = F' * ((F*y_fast) .* conj(x));       % N×M (delay×pulse)
        AFfs_raw_mc = abs( fft(range_mat, M, 2) ).^2;   % 沿慢时间(第2维)FFT -> N×M
        AFfs_raw_acc = AFfs_raw_acc + AFfs_raw_mc;

        % ---------- (B) 逐 Doppler bin 的快时间 Doppler 补偿：AFfs_comp ----------
        AFfs_comp_mc = zeros(N, M);

        n_fast = (0:N-1).';               % fast-time sample index within each pulse
        p_idx  = (0:M-1).';               % pulse index

        for k = 1:M
            % 将 k 映射到居中的 Doppler bin（覆盖正负频率）
            k0 = k-1;
            if k0 > M/2
                k0 = k0 - M;
            end

            % slow-time FFT bin k0 对应的“每采样”归一化频率 cycles/sample
            % fd = k0 / (M * N * 2)
            fd = k0 / (M*N*2);

            % 只补偿“快时间内”的相位滚动，用于抑制 ICI
            rot_fast = exp(-1i*2*pi*fd*n_fast);         % N×1
            y_fast_comp = y_fast .* rot_fast;           % N×M (隐式扩展)

            % 快时间 OFDM 匹配（子载波域乘 conj(x)）+ IFFT 到 delay
            range_mat_comp = F' * ((F*y_fast_comp) .* conj(x));  % N×M

            % 取慢时间 FFT 的第 k 个 bin（等价相干合成该 Doppler 点）
            steering = exp(-1i*2*pi*p_idx*(k-1)/M);     % M×1
            AFfs_comp_mc(:,k) = abs(range_mat_comp * steering).^2;
        end

        AFfs_comp_acc = AFfs_comp_acc + AFfs_comp_mc;
    end

    %% -------- MC 平均 --------
    AFfs_raw  = AFfs_raw_acc  / Nmc;
    AFfs_comp = AFfs_comp_acc / Nmc;

    %% -------- PSR（raw vs comp）--------
    peak_raw = max(AFfs_raw(:));
    rec_raw(idx_nu) = 10*log10( peak_raw / (sum(AFfs_raw(:)) - peak_raw) * (M*N+1) );

    peak_comp = max(AFfs_comp(:));
    rec_comp(idx_nu) = 10*log10( peak_comp / (sum(AFfs_comp(:)) - peak_comp) * (M*N+1) );

    %% -------- 保存 nu_show 的 2D 图数据（用于画图）--------
    if nu == nu_show
        AFfs_raw_show  = AFfs_raw;
        AFfs_comp_show = AFfs_comp;
        cAF_show = cAF;
    end
end

%% ===================== 画 PSR 曲线对比 =====================
figure;
plot(nu_list, rec_raw - max(rec_raw), 'o-', 'LineWidth', 1.2); hold on;
plot(nu_list, rec_comp - max(rec_comp), 's-', 'LineWidth', 1.2);
xlabel('Doppler Shift Index');
ylabel('PSR (normalized to each max) / dB');
title('PSR vs Doppler: BEFORE vs AFTER fast-time Doppler compensation');
grid on;
legend('Raw AFfs', 'Compensated AFfs', 'Location', 'best');

figure;
plot(nu_list, rec_raw, 'o-', 'LineWidth', 1.2); hold on;
plot(nu_list, rec_comp, 's-', 'LineWidth', 1.2);
xlabel('Doppler Shift Index');
ylabel('PSR / dB');
title('PSR vs Doppler (absolute): BEFORE vs AFTER compensation');
grid on;
legend('Raw AFfs', 'Compensated AFfs', 'Location', 'best');

%% ===================== 画 nu_show 下的 AFfs 2D 对比 =====================
if ~isempty(AFfs_raw_show)
    AFfs_raw_dB  = 10*log10(AFfs_raw_show  ./ max(AFfs_raw_show(:)));
    AFfs_comp_dB = 10*log10(AFfs_comp_show ./ max(AFfs_comp_show(:)));

    figure;
    imagesc(ifftshift(AFfs_raw_dB,2));
    axis xy; colormap jet; colorbar; caxis([-30 0]);
    xlabel('Doppler Bin (shifted)');
    ylabel('Delay Bin');
    title(['AFfs BEFORE compensation (nu = ', num2str(nu_show), ')']);

    figure;
    imagesc(ifftshift(AFfs_comp_dB,2));
    axis xy; colormap jet; colorbar; caxis([-30 0]);
    xlabel('Doppler Bin (shifted)');
    ylabel('Delay Bin');
    title(['AFfs AFTER compensation (nu = ', num2str(nu_show), ')']);

    figure;
    imagesc(ifftshift(AFfs_comp_dB - AFfs_raw_dB,2));
    axis xy; colormap jet; colorbar;
    xlabel('Doppler Bin (shifted)');
    ylabel('Delay Bin');
    title(['AFfs improvement (comp - raw) in dB (nu = ', num2str(nu_show), ')']);
end

%% ===================== （可选）保留你原来的 AF2/AF3 对照：用 nu_show 的 cAF_show =====================
if ~isempty(cAF_show)
    % 注意：cAF_show 是对 mc 累加后的 abs(AF)^2，因此平均时建议除以 Nmc
    cAF_avg = cAF_show / Nmc;

    %% 非相干累积（AF3）
    AF3 = zeros(N, N*2);
    for i = 1:M
        AF3 = AF3 + cAF_avg(:, i:M:end);
    end
    AF3 = AF3 / M;

    figure;
    imagesc(ifftshift(10*log10(AF3./max(AF3(:)))));
    axis xy; colormap jet; colorbar; caxis([-30 0]);
    xlabel('Doppler Index (Cyclic Shifted)');
    ylabel('Delay Index');
    title(['After Non-Coherent Integration (nu = ', num2str(nu_show), ')']);

    %% 未进行非相干累积（AF2）：按脉冲块累加
    AF2 = zeros(N, N*2);
    for i = 1:M
        AF2 = AF2 + cAF_avg(:, N*2*(i-1)+1:N*2*i);
    end
    AF2 = AF2 / M;

    figure;
    imagesc(ifftshift(10*log10(AF2./max(AF2(:)))));
    axis xy; colormap jet; colorbar; caxis([-30 0]);
    xlabel('Doppler Index (Cyclic Shifted)');
    ylabel('Delay Index');
    title(['Before Non-Coherent Integration (nu = ', num2str(nu_show), ')']);

    %% zero Doppler cut
    zero_doppler_AF2 = AF2(:, 1);
    zero_doppler_AF3 = AF3(:, 1);
    figure;
    plot(10*log10(abs(zero_doppler_AF2./max(zero_doppler_AF2))), 'r', 'LineWidth', 1.5); hold on;
    plot(10*log10(abs(zero_doppler_AF3./max(zero_doppler_AF3))), 'b', 'LineWidth', 1.5);
    xlabel('Delay Index'); grid on;
    legend('AF2', 'AF3');
    title(['Zero Doppler Cut (nu = ', num2str(nu_show), ')']);

    %% zero delay cut
    doppler_index = -(N-1):N;
    zero_delay_AF2 = AF2(tau, :);
    zero_delay_AF3 = AF3(tau, :);
    figure;
    plot(doppler_index, ifftshift(10*log10(abs(zero_delay_AF2./max(zero_delay_AF2)))), 'r', 'LineWidth', 1.5); hold on;
    plot(doppler_index, ifftshift(10*log10(abs(zero_delay_AF3./max(zero_delay_AF3)))), 'b', 'LineWidth', 1.5);
    xlabel('Doppler Index'); grid on;
    legend('AF2', 'AF3');
    title(['Zero Delay Cut (nu = ', num2str(nu_show), ')']);
end

如果你想进一步“更精准地对齐”你逐点匹配 AF 里用的 Doppler 栅格（比如把 k 映射成和 j 一样的循环频率定义），我也可以帮你把 fd 的定义改成严格同一套频率采样，这样补偿会更贴近你 point-wise AF 的 Doppler 轴。

問題

逐点的af不用留着 就比较我让你比较的就好

思考

Delivering MATLAB script for AFfs comparison

The user wants a complete MATLAB script comparing two schemes, excluding pointwise AF components. The script should compute and plot AFfs for a given nu_show value and calculate PSR vs nu. I'll ensure the dimensions are correct, use built-in fft/ifft, avoid F2 line, and allow flexible parameter scanning. I'll keep it straightforward and clear.

Rewriting AFfs computation and compensation

The initial AFfs formula seems off. I'll reframe it by using range_mat = F' * ((F*y_fast) .* conj(x)) and compute AFfs_raw with fft(range_mat, M, 2). For compensation, we rotate y_fast for each Doppler bin and compute range_mat_comp. We will vectorize for efficiency. Fast-time rotation needs to use the correct sample indices.

Fast-time index correction and phase handling

When selecting y2(N+1:end,:), I'm choosing fast-time indices n from N..2N-1, but the compensation uses n_fast from 0..N-1. The shift results in a constant phase factor exp(-j2π fd N), which is the same across pulses. This doesn't affect Doppler combination, just the global phase. Each pulse has a consistent phase, so the de-rotation doesn't need to account for it.

Mapping Doppler bins and slow-time FFT

I’m mapping Doppler bins to fd_target = nu/(M*2N), and the compensation involves using nu/L. For slow-time FFT, the frequency has a phase increment exp(j2π fd *2N). The corresponding slow-time normalized frequency is fd*2N, meaning the slow-time FFT bin aligns with nu. If nu is an integer, this matches the FFT bin index perfectly.

Aligning FFT bins and compensation

The M-point FFT bins correspond to k0/M cycles per pulse, and if nu lies within [-M/2, M/2-1], it directly maps. For compensation, fd uses k0/(M*2N) with centered bins, and shifting nu correctly matches this. I’ll include PSR plots, 2D images, and optional loops for efficiency. No need for parfor now.

Thought for 51s

回答

下面是完整精简版脚本：只计算并对比你要的 相干积分前 AFfs（2D-FFT AF） 的两种方案：

  • 方案A（raw）：不做快时间 Doppler 补偿，直接 2D-FFT 得到 AFfs_raw
  • 方案B（comp）：对每个 Doppler bin，先做快时间 Doppler 去旋转（抑制 ICI），再得到 AFfs_comp

并输出：

  1. PSR 随 nu 扫描的曲线对比
  2. 指定 nu_show 下的 AFfs_raw / AFfs_comp 2D 图对比 + 提升量图

直接复制保存为 .m 运行即可（逐点 AF 已删除）。

matlab
clc; clear;
format long

%% ===================== 参数区 =====================
times = 1;
M = 16;                 % slow-time pulses
N = 32*times;           % fast-time samples per pulse (useful part)
Nmc = 10;               % Monte Carlo runs
rng(1545);

tau = 4;                % delay (samples)
nu_list = -M:M;         % Doppler index scan（和你原来 range1 一致）
nu_show = 8;            % 用于画 2D AFfs 对比的 Doppler（必须在 nu_list 里）

F = dftmtx(N)./sqrt(N); % unitary DFT (N×N)

% 记录 PSR
rec_raw  = zeros(length(nu_list),1);
rec_comp = zeros(length(nu_list),1);

% 保存 nu_show 下的 2D AFfs 用于画图
AFfs_raw_show  = [];
AFfs_comp_show = [];

%% ===================== 主循环：扫 Doppler =====================
for idx_nu = 1:length(nu_list)
    nu = nu_list(idx_nu);

    % 2D-FFT AFfs 的 MC 累加容器（raw & comp）
    AFfs_raw_acc  = zeros(N, M);
    AFfs_comp_acc = zeros(N, M);

    for mc = 1:Nmc
        %% -------- OFDM gen --------
        s = randi(16,[N,M]) - 1;     % N×M
        x = pskmod(s,16);            % N×M (subcarrier symbols)

        % time-domain OFDM (repeat twice -> 2N×M)
        y = repmat(F' * x, [2,1]);   % (2N)×M

        %% -------- delay + Doppler（快时间 Doppler 导致 ICI）--------
        yy1 = y(:);
        yy2 = circshift(yy1, tau);

        L = M * N * 2; % total samples
        yy3 = yy2 .* exp(1i*2*pi*(0:L-1)'/L * nu);  % global fast-time Doppler

        %% -------- reshape received --------
        y2 = reshape(yy3, N*2, M);   % (2N)×M

        % 取后半段（与你原来 temp4 / y2(N+1:end,:) 一致）
        y_fast = y2(N+1:end, :);     % N×M

        %% ============================================================
        % 方案A：补偿前 AFfs_raw（2D-FFT）
        % 处理链：快时间FFT->乘conj(x)->快时间IFFT(=F') 得delay像 -> 慢时间FFT
        %% ============================================================
        range_mat = F' * ((F*y_fast) .* conj(x));      % N×M (delay×pulse)
        AFfs_raw_mc = abs( fft(range_mat, M, 2) ).^2;  % N×M (delay×DopplerBin)
        AFfs_raw_acc = AFfs_raw_acc + AFfs_raw_mc;

        %% ============================================================
        % 方案B：逐 Doppler bin 快时间 Doppler 补偿 AFfs_comp
        % 对每个 Doppler bin k：
        %   先 y_fast(n,p) * exp(-j2π fd(k) n) 去旋转（只补偿快时间相位滚动）
        %   再做与方案A同样的快时间处理
        %   最后取慢时间FFT的第k个bin（steering 相干合成）
        %% ============================================================
        AFfs_comp_mc = zeros(N, M);

        n_fast = (0:N-1).';          % fast-time index within each pulse
        p_idx  = (0:M-1).';          % pulse index

        for k = 1:M
            % centered bin index k0 in [-M/2, M/2-1]
            k0 = k-1;
            if k0 > M/2
                k0 = k0 - M;
            end

            % slow-time bin k0 对应的“每采样”归一化频率 cycles/sample
            % 因为 pulse-to-pulse 间隔是 2N 个采样点：
            %   fd = k0 / (M * 2N)
            fd = k0 / (M * N * 2);

            % 快时间 Doppler 去旋转（抑制 ICI）
            rot_fast = exp(-1i*2*pi*fd*n_fast);    % N×1
            y_fast_comp = y_fast .* rot_fast;      % N×M (implicit expand)

            % 与方案A同样的快时间处理得到 delay×pulse
            range_mat_comp = F' * ((F*y_fast_comp) .* conj(x));  % N×M

            % 取慢时间 FFT 第 k 个 Doppler bin（相干合成）
            steering = exp(-1i*2*pi*p_idx*(k-1)/M); % M×1
            AFfs_comp_mc(:,k) = abs(range_mat_comp * steering).^2;
        end

        AFfs_comp_acc = AFfs_comp_acc + AFfs_comp_mc;
    end

    %% -------- MC 平均 --------
    AFfs_raw  = AFfs_raw_acc  / Nmc;
    AFfs_comp = AFfs_comp_acc / Nmc;

    %% -------- PSR（raw vs comp）--------
    peak_raw = max(AFfs_raw(:));
    rec_raw(idx_nu) = 10*log10( peak_raw / (sum(AFfs_raw(:)) - peak_raw) * (M*N+1) );

    peak_comp = max(AFfs_comp(:));
    rec_comp(idx_nu) = 10*log10( peak_comp / (sum(AFfs_comp(:)) - peak_comp) * (M*N+1) );

    %% -------- 保存 nu_show 用于 2D 对比图 --------
    if nu == nu_show
        AFfs_raw_show  = AFfs_raw;
        AFfs_comp_show = AFfs_comp;
    end
end

%% ===================== PSR 曲线对比 =====================
figure;
plot(nu_list, rec_raw,  'o-', 'LineWidth', 1.2); hold on;
plot(nu_list, rec_comp, 's-', 'LineWidth', 1.2);
xlabel('Doppler Shift Index \nu');
ylabel('PSR / dB');
title('PSR vs Doppler: AFfs BEFORE vs AFTER fast-time Doppler compensation');
grid on;
legend('Raw AFfs', 'Compensated AFfs', 'Location', 'best');

figure;
plot(nu_list, rec_raw - max(rec_raw),  'o-', 'LineWidth', 1.2); hold on;
plot(nu_list, rec_comp - max(rec_comp),'s-', 'LineWidth', 1.2);
xlabel('Doppler Shift Index \nu');
ylabel('PSR (normalized) / dB');
title('PSR vs Doppler (normalized): Raw vs Compensated');
grid on;
legend('Raw AFfs', 'Compensated AFfs', 'Location', 'best');

%% ===================== nu_show 下的 AFfs 2D 对比 =====================
if ~isempty(AFfs_raw_show)
    AFfs_raw_dB  = 10*log10(AFfs_raw_show  ./ max(AFfs_raw_show(:)));
    AFfs_comp_dB = 10*log10(AFfs_comp_show ./ max(AFfs_comp_show(:)));

    figure;
    imagesc(ifftshift(AFfs_raw_dB, 2));  % 只对 Doppler 维做 shift
    axis xy; colormap jet; colorbar; caxis([-30 0]);
    xlabel('Doppler Bin (shifted)');
    ylabel('Delay Bin');
    title(['AFfs BEFORE compensation (nu = ', num2str(nu_show), ')']);

    figure;
    imagesc(ifftshift(AFfs_comp_dB, 2));
    axis xy; colormap jet; colorbar; caxis([-30 0]);
    xlabel('Doppler Bin (shifted)');
    ylabel('Delay Bin');
    title(['AFfs AFTER compensation (nu = ', num2str(nu_show), ')']);

    figure;
    imagesc(ifftshift(AFfs_comp_dB - AFfs_raw_dB, 2));
    axis xy; colormap jet; colorbar;
    xlabel('Doppler Bin (shifted)');
    ylabel('Delay Bin');
    title(['AFfs improvement (comp - raw) in dB (nu = ', num2str(nu_show), ')']);
end

如果你希望把 “逐 Doppler 点补偿”的 fd 定义改成严格按你原来 Doppler index（nu）那套去对齐（包括循环频率/shift 方式/频率标定），我也可以把 Doppler 轴从 bin 映射到 -(M-1):M 或者按你 nu_list 的定义重建，这样对比会更“同一坐标系”。

問題

我们希望的psr应该是越低越好么

回答

不是。按你代码里这个 PSR 定义，PSR 越大越好（越高表示主峰越“尖”、旁瓣/杂波越小）。

你这里算的是（本质上）：

PSR(dB)=10log⁡10(Pmax⁡∑P−Pmax⁡⋅(MN+1))\text{PSR(dB)}=10\log_{10}\Bigg(\frac{P_{\max}}{\sum P - P_{\max}}\cdot (MN+1)\Bigg)
  • Pmax⁡P_{\max}：AF 图里的最大峰值（主峰）
  • ∑P−Pmax⁡\sum P - P_{\max}：除主峰外的总能量（旁瓣/杂散/噪声底）
  • 乘的 (MN+1)(MN+1) 是个尺度因子（把比值做个归一化，便于不同尺寸比较）

因此：

  • 主峰更高、旁瓣/杂散更低 ⇒ 分子/分母更大 ⇒ PSR 更大（更好）
  • ICI 更严重时通常会把能量“摊开”，让 ∑P−Pmax⁡\sum P - P_{\max} 变大、Pmax⁡P_{\max} 变小 ⇒ PSR 变小（更差）

你图里如果出现 rec-max(rec) 全是负值，那只是做了“相对最大值归一化”，0 dB 那个点是最好，越往下越差；但在“绝对 PSR（rec）”那张图里，仍然是数值越大越好。

如果你想让“越低越好”，那通常用的是 ISLR/PSLR 这类“旁瓣比”指标（例如旁瓣能量/主瓣能量），那种确实是越低越好；但你现在这个 PSR 不是那个定义。

問題

image.pngimage.png

可是从结果来看 补偿后的效果更差 是为什么

思考

Considering compensation and de-rotation in fast and slow time

Okay, so slow-time bin index k maps to nu mod M. The mapping works well, but I’m wondering—there seems to be a potential issue with double-compensation. We rotate in both fast time and slow time, but we’re not aligning fast and slow rotations at the same time. I'll need to reassess how these two interact for better accuracy.

Adjusting de-rotation application for better PSR

The lower PSR might be due to how we computed the range matrix with fast-time compensation, but without adjusting for the global indices in each pulse. I think the de-rotation should happen before the FFT, and might need to apply to the time-domain before CP removal—especially for the second half of the pulse, where offsets should be considered.

Addressing pulse offset and phase factors

Okay, the constant offset produces a phase factor that affects the entire pulse, and it stays consistent across subcarriers — no real issue unless I misapply Doppler shifts later. The key seems to be accounting for pulse index in de-rotation. The Doppler term should take into account the full pulse index, not just the fast time.

Investigating PSR drop with de-rotation

When de-rotating by exp(-j2π fd n), we remove intra-pulse Doppler, leaving a phase that’s slow-time linear with p, which the FFT handles well. But the drop in PSR could be due to "focusing" non-matching bins, increasing sidelobes. The peak-to-sum ratio (PSR) depends on this, so either the peak could decrease or the sidelobe energy could increase.

Considering trade-offs in compensation and PSR

The bin-by-bin compensation could increase the noise floor by breaking the orthogonality between bins, leading to higher sidelobe energy. This alters the overall energy distribution, which affects PSR. While evaluating the true Doppler bin might show a sharper peak, comparing across bins could result in a worse PSR due to increased sidelobes. It seems evaluating PSR per bin might be a better approach.

Clarifying compensation behavior and mapping issues

The mapping and compensation seem tricky, especially at Nyquist (nu=8), where handling aliasing can cause ambiguity. It looks like this results in a reduction of peak values, even at nu=0, as compensation suppresses energy. The current approach might deviate from raw results when fd=0, but ideally, they should match if handled correctly.

Considering compensation effects on PSR

For bin 1, the compensation matches the raw result; however, PSR differs due to how other bins are processed. Fast-time de-rotation for bins with non-zero fd may increase sidelobe energy. To keep comparisons fair, we should apply compensation across all doppler bins or try compensating in the frequency domain.

Explaining PSR decrease with compensation

The PSR measurement is global, so compensation can change processing across bins, affecting energy across them. The true global Doppler for compensation might need adjustment. Our compensation introduces additional sidelobe energy, as seen in the decreased PSR. This may be due to mismatched processing across bins, increasing total energy.

Addressing PSR and frequency-domain compensation

The columns in AFfs_comp don't form a unitary transform, so energy may get skewed across bins, affecting PSR. To assess improvement, it's better to compare at the target Doppler bin or use a focused window for PSR. Compensation for ICI should be done in the frequency domain, correcting for CFO once instead of per Doppler bin.

Improving PSR and compensation approach

To avoid increasing sidelobes and unfair results, it’s better to perform CFO compensation for the specific Doppler bin (nu) first, only if it's known. If unknown, do a search for max PSR. Instead of calculating across the whole map, compare PSR at the center or mainlobe.

Refining compensation approach for PSR improvement

The key issue lies in how compensation adjusts Doppler bins. While targeting a specific bin with compensation is beneficial, it seems we may introduce extra energy in other bins, lowering PSR. To improve, we should focus on coherent summation at zero Doppler after compensation without overemphasizing steering, which may unfocus high bins.

Addressing Doppler compensation and PSR issues

It seems the main issue lies in how Doppler compensation is applied—fast-time rotation is being mismatched with slow-time compensation, resulting in aliasing and increased sidelobe energy. To fix this, applying compensation in the subcarrier domain before fast-time rotation may help, and we should compute PSR focused on the correct Doppler or adjust for aliasing effects.

Exploring Doppler Compensation and PSR Behavior

The high PSR peaks at nu=0 and ±1 suggest ICI due to CFO, affecting energy calculations. A correction method using fast-time compensation is proposed but concerns arise because columns aren’t orthogonal post-compensation. To evaluate better, we might compute PSR per Doppler hypothesis and compare per-bin results, normalizing total energy.

Proposing improvements for ICI and PSR calculations

It seems the compensation per Doppler bin is causing increased energy due to mismatched CFO corrections. To address that, we could propose measuring PSR differently—using PSLR (peak-to-second-peak ratio) or ISLR (sidelobe-to-mainlobe ratio) instead of summing across all bins. A single CFO correction based on true Doppler might also improve outcomes.

Thought for 4m 50s

回答

PSR 按你现在这套定义（主峰 / 其余能量）是 越大越好。你现在看到“补偿后 PSR 更低”，通常不是“补偿把主峰做坏了”，而是补偿把整张 AF 图的“其余能量”抬高了，导致分母变大 → PSR 变小。

下面是最常见、也最符合你这张曲线形状的原因（按重要性排序）：

1) 你做的“逐 Doppler 点补偿”会破坏 Doppler 维的正交性（能量不再守恒）

原始 AFfs_raw 是对 同一个 range_mat 做慢时间 FFT 得到所有 Doppler bin：

  • 这是一个线性、（近似）正交的变换（FFT）
  • “不属于该 bin 的能量”会比较干净地被压到别的地方

但你的 AFfs_comp 是：

每个 Doppler bin 都先用不同的快时间去旋转，再各算一列。

这等价于做了一个“多假设匹配滤波器组”，每一列用的都是不同预处理后的数据，于是：

  • 很多“本来在 raw 里会被正交性压下去的东西”，在 comp 里可能被某些假设 意外对齐，产生更大的输出
  • 结果就是：整张 AFfs_comp 的“旁瓣/杂散总能量”变大
    即使主峰没怎么变，PSR 也会下降

这也是你图里“补偿后整条曲线都更低”的典型现象。

2) 你现在的 PSR 是对“整张 2D 图”做的，全图分母会被“错误 Doppler 假设”污染

你 PSR 的分母是 sum(AFfs(:)) - peak（全图除主峰外的所有能量）。

但你逐 Doppler 点补偿后，错误 Doppler 假设对应的列可能被抬高很多；这些能量全部进入分母，PSR 自然变差。

换句话说：你现在的 PSR 更像是在惩罚“整张图是不是很稀疏”。
而逐 Doppler 点补偿本质是“每列都做一次匹配”，会让图变得没那么稀疏。

3) 你扫的 nu_list = -M:M 本身在慢时间（M 点）上是混叠的

慢时间只有 M=16 个 bin，但你扫了 33 个 nu。
nu 在慢时间上会 mod M 混叠，但它在快时间造成的频偏（ICI 强弱）并不等价，所以你用“bin 对应的 fd(k)”去补偿时，某些 nu 会出现“补偿不对味”。

这会让“补偿后曲线更平滑、更低”更明显。

怎么验证“补偿其实在正确 Doppler 上是有效的”

最直接的验证方法：不要用全图 PSR，而是只在“目标 Doppler 列”上算一个 delay 维的 PSR（或 PSLR）。

对每个真值 nu，它对应的慢时间 Doppler bin 是：

matlab
k_true = mod(nu, M);
if k_true < 0, k_true = k_true + M; end
k_true = k_true + 1; % MATLAB index

然后只用这一列算“列内 PSR”：

matlab
col_raw  = AFfs_raw(:,  k_true);
col_comp = AFfs_comp(:, k_true);

psr_col_raw  = 10*log10( max(col_raw)  / (sum(col_raw)  - max(col_raw))  * (N+1) );
psr_col_comp = 10*log10( max(col_comp) / (sum(col_comp) - max(col_comp)) * (N+1) );

你大概率会看到：列内指标 comp ≥ raw（尤其是 Doppler 不小的时候），说明快时间补偿确实在抑制 ICI；只是“全图 PSR”被错误假设的列能量抬高给拖累了。

如果你坚持要“全图 PSR”也变好，该怎么改？

想让全图 PSR 也提升，你需要保证 Doppler 维仍然是对同一个数据矩阵做一个一致的变换。两条常用路线：

  1. 先估计/假设一个 Doppler（或 CFO），对整段快时间统一补偿一次，再正常做 2D-FFT（这会真正减少 ICI 且不引入“每列不同预处理”的副作用）。
  2. 做真正的 ICI-aware 等化/解耦（例如 ICI 矩阵补偿），而不是“每个 bin 独立旋转”。

如果你愿意，我可以在你那份“完整代码”里直接加上 (a) 全图 PSR + (b) 目标列 PSR 两套曲线一起画出来，你就能一眼看到：全图 PSR 为什么掉、但目标列为什么其实提升了。

問題

我们这样 不在原来的代码里面补偿，在我下面这个为基础来补偿：function [AFfs, tau_idx, nu_idx, peak_vals] = calc_AFfs_multi( ...
N, M, Nmc, tau, nu, flag, K, alpha, Kest, guard_r, guard_c)
% calc_AFfs_multi (N-structure, no 2N)
% 多目标 Monte-Carlo 期望 |RDM|^2，标准 2D-FFT 处理

text
tau   = tau(:).';
nu    = nu(:).';
if nargin < 8 || isempty(alpha),   alpha   = 1; end
if nargin < 9 || isempty(Kest),    Kest    = max(1, numel(tau)); end
if nargin < 10 || isempty(guard_r), guard_r = 2; end
if nargin < 11 || isempty(guard_c), guard_c = 2; end
alpha = alpha(:).';

AF_sum = zeros(N,M);

rng(1545,'twister');
seeds = randi(2^31-1, Nmc, 1, 'uint32');

parfor mc = 1:Nmc
    rng(double(seeds(mc)),'twister');

    [Y, Xtd, Xf] = simulate_ofdm_echo_multi(N,M,tau,nu,flag,K,alpha);

    % 快时间 FFT -> 频域匹配滤波 -> 快时间 IFFT
    Yf = fft(Y, [], 1);                    % N×M
    R  = ifft( Yf .* conj(Xf), [], 1 );    % N×M，range 压缩

    % 慢时间 FFT -> Doppler 轴
    RDM = fft(R, [], 2);                   % N×M

    AF_trial = abs(RDM).^2;
    AF_sum   = AF_sum + AF_trial;
end

AFfs = AF_sum / Nmc;

[tau_idx, nu_idx, peak_vals] = topK_peaks(AFfs, Kest, guard_r, guard_c);

end

%% ---------- Top-K 峰值 + 保护单元 ----------

function [r_idx, c_idx, vals] = topK_peaks(A, K, guard_r, guard_c)
Awork = A;
r_idx = zeros(K,1);
c_idx = zeros(K,1);
vals = zeros(K,1);

text
for k = 1:K
    [val, lin] = max(Awork(:));
    if ~isfinite(val) || val <= 0
        r_idx = r_idx(1:k-1);
        c_idx = c_idx(1:k-1);
        vals  = vals(1:k-1);
        break;
    end
    [r, c] = ind2sub(size(Awork), lin);
    r_idx(k) = r;
    c_idx(k) = c;
    vals(k)  = val;

    r1 = max(1, r-guard_r); r2 = min(size(Awork,1), r+guard_r);
    c1 = max(1, c-guard_c); c2 = min(size(Awork,2), c+guard_c);
    Awork(r1:r2, c1:c2) = -Inf;
end

end

function [Y, Xtd, Xf] = simulate_ofdm_echo_multi(N,M,tau_vec,nu_vec,flag,K,amp)
tau_vec = tau_vec(:).';
nu_vec = nu_vec(:).';
amp = amp(:).';

text
P = max([numel(tau_vec), numel(nu_vec), numel(amp)]);
if numel(tau_vec)==1, tau_vec = repmat(tau_vec,1,P); end
if numel(nu_vec)==1,  nu_vec  = repmat(nu_vec, 1,P); end
if numel(amp)==1,     amp     = repmat(amp,    1,P); end

tau_vec = mod(round(tau_vec), N);
nu_vec  = mod(round(nu_vec),  M);

s = randi(K,[N,M]) - 1;
if flag == 1
    Xf = pskmod(s, K);
else
    Xf = qammod(s, K);
end

Xtd = sqrt(N) * ifft(Xf, [], 1);

k = (0:N-1).';
m = 0:M-1;

Y = complex(zeros(N,M));
for p = 1:P
    tau_p = tau_vec(p);
    nu_p  = nu_vec(p);

    X_shift = circshift(Xtd, tau_p, 1);

    dop_fast = exp(1j*2*pi*nu_p * k/(M*N));
    dop_slow = exp(1j*2*pi*nu_p * m/M);
    Dop = dop_fast * dop_slow;

    Y = Y + amp(p) * (X_shift .* Dop);
end

end
看看应该补偿在哪一行是比较合理的，我认为是要和DAFT补偿的逻辑一样的

思考

Considering Doppler compensation and DAFT logic

For per-Doppler bin compensation, you multiply Y by a fast-time correction term, then perform the FFT and match filter. But applying different Y for each bin can break orthogonality, affecting PSR. Regarding DAFT, it could refer to a method in radar dealing with Doppler-delay coupling, like "Doppler-Aware Fourier Transform" or "Discrete Affine Fourier Transform."

Explaining DAFT-like compensation

To compensate similarly to DAFT, we multiply by a phase term before the FFT due to the linear phase introduced by fast-time Doppler. After FFT, compensating with a circular shift in frequency is another approach. This shift is tied to how the phase causes Inter-Carrier Interference (ICI), especially when ν is not a multiple of M.

Explaining compensation for Doppler processing

To compensate for frequency offset, you multiply by exp(-j2π ν_h k/(M N)) after simulating the echo. This requires broadcasting across the columns in the matrix. For preserving slow-time Doppler, we don’t remove the slow-phase term. To compute the RDM per Doppler bin, apply per-bin compensation before FFT.

Considering fast-time compensation and DAFT-style logic

Fast-time compensation should be applied to Y before the fast FFT, requiring looping over l for separate RDM slices. This results in a matrix with separate compensations. The user is likely asking for comparisons before coherent integration (before the slow FFT). The best place to apply DAFT-style compensation is after simulating the OFDM echo, before the FFT. The compensation can be done by pre-multiplying in the time domain to restore orthogonality.

Proposing two compensation versions and code modifications

The code will define a phase matrix based on the Doppler shift and multiply it by Y. For compensation, map nu_h to a centered range and apply the phase shift. We can implement two versions:

  1. Single compensation with an estimated nu_hat for global PSR improvement.
  2. Per-bin compensation to compute each Doppler bin output.

We’ll modify the function to accept a 'comp_mode' flag for different compensation modes, ensuring flexibility for the user’s needs.

Providing compensation function and comparison

To compare both compensation schemes, I’ll create a new function calc_AFfs_multi_comp that returns both raw and compensated AFfs outputs. The function handles fast-time compensation for each Doppler bin and computes the RDM accordingly, adjusting for Doppler shifts with a steering matrix. We'll also emphasize using k as the fast-time index.

Addressing issue with nu mapping and propose changes

The simulation uses mod(round(nu_vec), M) but this doesn't allow for negative nu values. To fix this, we should modify the code to handle negative nu properly using centered mapping like nu_vec = mod(round(nu_vec) + floor(M/2), M) - floor(M/2). This ensures the Doppler axis is consistent for simulation and compensation. We'll also propose changes to handle signed nu values and add compensation logic.

Explaining frequency domain compensation and updates

Compensation in the frequency domain could involve circularly shifting Yf. However, Doppler compensation in the time domain is simpler. I’ll show an updated function, calc_AFfs_multi_compare, that returns both the raw and compensated results. The compensation targets the dop_fast term, and I’ll suggest evaluating the PSR after applying single compensation or per-bin compensation. We’ll include both approaches in the final code.

Explaining compensation adjustments and mapping

The user asked where compensation should be applied. After simulating, compensation can be inserted after the line where we define Y.

For consistency with simulation, nu_h = c-1 should be used, without centered mapping. If negative nu_h is needed, it should be mapped properly. Adjust the simulation logic to allow signed values or use modulo, ensuring consistency across internal operations and the compensation process.

Providing options for Doppler compensation approach

They want compensation for each Doppler point, but need to decide how to handle negative values. Two approaches:

  1. Keep simulation wrapped (modulo) as is, using nu_h = k - 1.
  2. Switch to signed mapping, adjusting both the simulation and compensation.

I'll provide both options in the calc_AFfs_multi_compare function, keeping it minimally modified. Also, I'll compute PSR separately and include performance optimizations like precomputing the steering matrix.

Thought for 2m 30s

回答

你这个版本里，最合理、也最“像 DAFT/预调制再FFT”的补偿位置就是：

在 simulate_ofdm_echo_multi 得到 Y 之后、做快时间 FFT (Yf = fft(Y,[],1)) 之前，对 Y 做快时间去旋转（de-rotation）。

原因很简单：ICI 是在你对含有快时间线性相位的信号做 FFT 时产生的（子载波正交性被破坏），所以要在 FFT 之前把这段线性相位“拿掉”，才能抑制 ICI。

并且你在仿真里快时间 Doppler 项就是这句生成的：

matlab
dop_fast = exp(1j*2*pi*nu_p * k/(M*N));

所以补偿项就应该是它的共轭（反相位）：

matlab
exp(-1j*2*pi*nu_h * k/(M*N))

下面给你一份基于你函数的完整“对比两种方案”的实现：

  • AFfs_raw：原始标准 2D-FFT
  • AFfs_comp：对每个 Doppler bin 单独做快时间补偿（你想要的方式）

⚠️ 关键点：你 simulate_ofdm_echo_multi 里把 nu 做了 mod(...,M)，也就是 Doppler 是 0…M-1 的“环形索引”。
所以下面补偿时 也用同样的 0…M-1 去补偿（否则会对不上，补偿会变差）。

完整代码：对比 raw vs per-Doppler 快时间补偿

matlab
function [AFfs_raw, AFfs_comp, peaks_raw, peaks_comp] = calc_AFfs_multi_compare( ...
        N, M, Nmc, tau, nu, flag, K, alpha, Kest, guard_r, guard_c)
% calc_AFfs_multi_compare
% 输出两套：标准2D-FFT的AFfs_raw，以及“逐Doppler点快时间补偿”的AFfs_comp

    tau   = tau(:).';
    nu    = nu(:).';
    if nargin < 8 || isempty(alpha),    alpha   = 1; end
    if nargin < 9 || isempty(Kest),     Kest    = max(1, numel(tau)); end
    if nargin < 10 || isempty(guard_r), guard_r = 2; end
    if nargin < 11 || isempty(guard_c), guard_c = 2; end
    alpha = alpha(:).';

    AF_sum_raw  = zeros(N, M);
    AF_sum_comp = zeros(N, M);

    % 用于 Doppler bin 提取的慢时间 steering（与 fft(R,[],2) 同相位定义）
    m = (0:M-1).';
    W = exp(-1j * 2*pi/M * (m * (0:M-1)));   % M×M，列对应 Doppler bin 0..M-1

    % 快时间索引
    k = (0:N-1).';

    rng(1545,'twister');
    seeds = randi(2^31-1, Nmc, 1, 'uint32');

    parfor mc = 1:Nmc
        rng(double(seeds(mc)),'twister');

        [Y, ~, Xf] = simulate_ofdm_echo_multi(N, M, tau, nu, flag, K, alpha); % Y: N×M

        %% ======================
        % (A) RAW：标准 2D-FFT
        %% ======================
        Yf  = fft(Y, [], 1);                 % N×M
        R   = ifft(Yf .* conj(Xf), [], 1);   % N×M (range)
        RDM = fft(R, [], 2);                 % N×M (doppler)
        AF_trial_raw = abs(RDM).^2;

        %% ==========================================
        % (B) COMP：逐 Doppler bin 的快时间补偿
        % 对每个 Doppler bin q=0..M-1:
        %   先 Y .* exp(-j2π q k/(M N))
        %   再走同样的“快时间FFT-匹配-IFFT”
        %   最后只取该 bin（用 steering 等价 fft 的第 q 列）
        %% ==========================================
        AF_trial_comp = zeros(N, M);

        for q = 0:M-1
            rot_fast = exp(-1j * 2*pi * q * k/(M*N));   % N×1
            Yc = Y .* rot_fast;                         % N×M（隐式扩展；旧版可用bsxfun）

            Yf_c = fft(Yc, [], 1);
            R_c  = ifft(Yf_c .* conj(Xf), [], 1);

            % 取 Doppler bin q（等价于 RDM(:,q+1)）
            AF_trial_comp(:, q+1) = abs(R_c * W(:, q+1)).^2;
        end

        AF_sum_raw  = AF_sum_raw  + AF_trial_raw;
        AF_sum_comp = AF_sum_comp + AF_trial_comp;
    end

    AFfs_raw  = AF_sum_raw  / Nmc;
    AFfs_comp = AF_sum_comp / Nmc;

    % 可选：输出 topK 峰值信息，方便你后续算PSR/定位
    [r1,c1,v1] = topK_peaks(AFfs_raw,  Kest, guard_r, guard_c);
    [r2,c2,v2] = topK_peaks(AFfs_comp, Kest, guard_r, guard_c);

    peaks_raw  = struct('tau_idx', r1, 'nu_idx', c1, 'vals', v1);
    peaks_comp = struct('tau_idx', r2, 'nu_idx', c2, 'vals', v2);
end


%% ---------- Top-K 峰值 + 保护单元 ----------
function [r_idx, c_idx, vals] = topK_peaks(A, K, guard_r, guard_c)
    Awork = A;
    r_idx = zeros(K,1);
    c_idx = zeros(K,1);
    vals  = zeros(K,1);

    for k = 1:K
        [val, lin] = max(Awork(:));
        if ~isfinite(val) || val <= 0
            r_idx = r_idx(1:k-1);
            c_idx = c_idx(1:k-1);
            vals  = vals(1:k-1);
            break;
        end
        [r, c] = ind2sub(size(Awork), lin);
        r_idx(k) = r;
        c_idx(k) = c;
        vals(k)  = val;

        r1 = max(1, r-guard_r); r2 = min(size(Awork,1), r+guard_r);
        c1 = max(1, c-guard_c); c2 = min(size(Awork,2), c+guard_c);
        Awork(r1:r2, c1:c2) = -Inf;
    end
end


function [Y, Xtd, Xf] = simulate_ofdm_echo_multi(N,M,tau_vec,nu_vec,flag,K,amp)
    tau_vec = tau_vec(:).';
    nu_vec  = nu_vec(:).';
    amp     = amp(:).';

    P = max([numel(tau_vec), numel(nu_vec), numel(amp)]);
    if numel(tau_vec)==1, tau_vec = repmat(tau_vec,1,P); end
    if numel(nu_vec)==1,  nu_vec  = repmat(nu_vec, 1,P); end
    if numel(amp)==1,     amp     = repmat(amp,    1,P); end

    tau_vec = mod(round(tau_vec), N);
    nu_vec  = mod(round(nu_vec),  M);   % 注意：这里把 Doppler 变成 0..M-1 的环形索引

    s = randi(K,[N,M]) - 1;
    if flag == 1
        Xf = pskmod(s, K);
    else
        Xf = qammod(s, K);
    end

    Xtd = sqrt(N) * ifft(Xf, [], 1);

    k = (0:N-1).';
    m = 0:M-1;

    Y = complex(zeros(N,M));
    for p = 1:P
        tau_p = tau_vec(p);
        nu_p  = nu_vec(p);

        X_shift = circshift(Xtd, tau_p, 1);

        dop_fast = exp(1j*2*pi*nu_p * k/(M*N));
        dop_slow = exp(1j*2*pi*nu_p * m/M);
        Dop = dop_fast * dop_slow;

        Y = Y + amp(p) * (X_shift .* Dop);
    end
end

你要补偿“在哪一行”？

就在你原函数里这两行之间：

matlab
[Y, Xtd, Xf] = simulate_ofdm_echo_multi(...);

Yf = fft(Y, [], 1);

补偿就是插在中间（对 Y 乘 exp(-j2π nu_h k/(MN))），这就是最标准的“DAFT式：先预相位/去相位，再做FFT”。

一个很重要的提醒（避免你再遇到“补偿反而更差”）

如果你用的是全图 PSR（主峰 / 其余所有能量），那么“逐 Doppler 假设”这种做法会让很多错误假设的列也被某些东西“部分对齐”，分母会变大，所以 PSR 可能看起来更差。

更公平的评价方式是：

  • 只在“真 Doppler 对应的那一列”算 1D-PSR/PSLR（列内），或
  • 用 topK_peaks 找到主峰后，在其保护单元外算 ISLR/PSLR。

如果你把你现在算 PSR 的方式贴出来（你希望用全图还是列内），我可以把对比指标那部分也给你改成“更符合你想看的 ICI 改善”的版本。

問題

把你现在这个思路 写成一个m文件，而不是函数

思考

Developing MATLAB script with Doppler compensation

I need to create a MATLAB script instead of a function and write it in Chinese. I'll include local functions at the end of the script. The goal is to compare raw and compensated AFfs values, compute per-bin Doppler compensation, and implement plotting. I'll use their variables, such as N, M, and Nmc, and apply the Doppler compensation inside the loop per bin.

Creating efficient MATLAB script for Doppler and PSR

To improve efficiency, instead of recomputing the full FFT each time, I’ll use steering. I'll write a script that sets parameters, defines a true nu value, and compares performance for various Doppler values by scanning through the range -M:M. I'll include two variants for PSR calculations: full map and column PSR.

Preparing MATLAB script with PSR metrics and parfor

The user asked for a script instead of a function, so I’ll include both PSR metrics to show a comparison between methods. The script will have local functions like simulate_ofdm_echo_multi, topK_peaks, and compute_psr. I’ll use a mixture of for loops and parfor to keep things efficient, ensuring proper variable slicing.

Implementing PSR with parfor and plotting

In the script, I'll use parfor with reduction on numeric arrays for AF_sum_raw and AF_sum_comp. I'll ensure deterministic results by replicating their entry for random seeds, as in their function. I’ll compute both full and column PSR for each nu_true, and use imagesc for plotting.

Thought for 48s

回答

下面给你一份完整的 MATLAB 脚本（.m 文件），不是函数；直接运行即可对比两套方案：

  • Raw（不补偿）：Y -> FFT快时间 -> 频域匹配 -> IFFT快时间 -> FFT慢时间
  • Comp（逐 Doppler bin 快时间补偿）：对每个 Doppler bin q，在快时间 FFT 之前做
      Y(k,m)←Y(k,m)exp⁡(−j2πqk/(MN))  \;Y(k,m)\leftarrow Y(k,m)\exp(-j2\pi q k/(MN))\;（DAFT/去旋转逻辑），再走同样链路，并只取该 bin 输出。

同时我给了两种 PSR：

  • PSR_full：你之前那种“全图 PSR”（容易出现你看到的“补偿后更差”现象）
  • PSR_col：只在“真值 Doppler 对应列”上算的 PSR（更能反映 ICI 抑制是否有效）

把下面保存为例如：demo_AFfs_fastDopplerComp.m 直接运行。

matlab
clc; clear; close all;
format long;

%% ===================== 参数区 =====================
N   = 32;        % fast-time samples per pulse
M   = 16;        % slow-time pulses
Nmc = 10;        % Monte-Carlo runs

tau0 = 4;        % delay
flag = 1;        % 1: PSK, else: QAM
Kmod = 16;       % modulation order
alpha = 1;       % target amplitude (single target)

nu_list = -M:M;  % 扫描的 Doppler index（注意：仿真里 nu 会 mod M）
nu_show = 8;     % 画2D图用的 Doppler index（必须在 nu_list 里）

% 保护单元（如果你想找峰值也可用；本脚本主要画图/算PSR）
guard_r = 2;
guard_c = 2;

% 快时间/慢时间索引
k_fast = (0:N-1).';
m_slow = (0:M-1).';

% slow-time FFT 的 steering（与 fft(R,[],2) 的相位定义一致）
W = exp(-1j * 2*pi/M * (m_slow * (0:M-1)));  % M×M

% 记录 PSR：两种口径（全图/真列）
PSR_full_raw  = zeros(numel(nu_list),1);
PSR_full_comp = zeros(numel(nu_list),1);
PSR_col_raw   = zeros(numel(nu_list),1);
PSR_col_comp  = zeros(numel(nu_list),1);

% 保存用于展示的 2D AFfs（nu_show）
AFfs_raw_show  = [];
AFfs_comp_show = [];

%% ===================== 主循环：扫 Doppler =====================
rng(1545,'twister');                       % 固定随机性（可复现）
seeds = randi(2^31-1, Nmc, 1, 'uint32');   % parfor 可复现

for idx = 1:numel(nu_list)
    nu0 = nu_list(idx);

    AF_sum_raw  = zeros(N,M);
    AF_sum_comp = zeros(N,M);

    parfor mc = 1:Nmc
        rng(double(seeds(mc)),'twister');

        % ---------- 生成回波（N-structure：Y 是 N×M） ----------
        [Y, ~, Xf] = simulate_ofdm_echo_multi(N, M, tau0, nu0, flag, Kmod, alpha);

        %% ============ 方案 A：RAW（标准2D-FFT） ============
        Yf  = fft(Y, [], 1);                  % 快时间 FFT
        R   = ifft(Yf .* conj(Xf), [], 1);    % 匹配滤波 + 快时间 IFFT（range 压缩）
        RDM = fft(R, [], 2);                  % 慢时间 FFT（Doppler）
        AF_trial_raw = abs(RDM).^2;

        %% ============ 方案 B：逐 Doppler bin 快时间补偿 ============
        % 关键：补偿位置 = 在快时间 FFT 之前，对 Y 乘 e^{-j2π q k/(MN)}
        AF_trial_comp = zeros(N,M);

        for q = 0:M-1
            rot_fast = exp(-1j * 2*pi * q * k_fast/(M*N));  % N×1
            % 兼容老版本 MATLAB：用 bsxfun 做维度扩展
            Yc = bsxfun(@times, Y, rot_fast);               % N×M

            Yf_c = fft(Yc, [], 1);
            R_c  = ifft(Yf_c .* conj(Xf), [], 1);

            % 只取 Doppler bin = q（等价于 fft(R_c,[],2) 的第 q+1 列）
            AF_trial_comp(:, q+1) = abs(R_c * W(:, q+1)).^2;
        end

        % parfor reduction
        AF_sum_raw  = AF_sum_raw  + AF_trial_raw;
        AF_sum_comp = AF_sum_comp + AF_trial_comp;
    end

    AFfs_raw  = AF_sum_raw  / Nmc;
    AFfs_comp = AF_sum_comp / Nmc;

    %% ===================== PSR 计算（两种口径） =====================
    % (1) 全图 PSR（你之前用的那种：peak / (sum-peak)）
    peak_raw  = max(AFfs_raw(:));
    peak_comp = max(AFfs_comp(:));

    PSR_full_raw(idx)  = 10*log10( peak_raw  / (sum(AFfs_raw(:))  - peak_raw )  * (M*N+1) );
    PSR_full_comp(idx) = 10*log10( peak_comp / (sum(AFfs_comp(:)) - peak_comp) * (M*N+1) );

    % (2) 真值 Doppler 对应列的 PSR（更贴近“ICI抑制是否有效”）
    % 注意：simulate_ofdm_echo_multi 内部 nu 会 mod M，所以真值列为：
    q_true = mod(round(nu0), M);         % 0..M-1
    col_raw  = AFfs_raw(:,  q_true+1);
    col_comp = AFfs_comp(:, q_true+1);

    pk1 = max(col_raw);
    pk2 = max(col_comp);
    PSR_col_raw(idx)  = 10*log10( pk1 / (sum(col_raw)  - pk1) * (N+1) );
    PSR_col_comp(idx) = 10*log10( pk2 / (sum(col_comp) - pk2) * (N+1) );

    %% 保存展示用 2D 图
    if nu0 == nu_show
        AFfs_raw_show  = AFfs_raw;
        AFfs_comp_show = AFfs_comp;
    end
end

%% ===================== 画 PSR 曲线对比 =====================
figure;
plot(nu_list, PSR_full_raw,  'o-', 'LineWidth', 1.2); hold on;
plot(nu_list, PSR_full_comp, 's-', 'LineWidth', 1.2);
grid on; xlabel('Doppler index \nu'); ylabel('PSR\_full / dB');
title('PSR\_full (全图口径)：Raw vs Per-bin Fast-time Compensation');
legend('Raw', 'Comp', 'Location', 'best');

figure;
plot(nu_list, PSR_col_raw,  'o-', 'LineWidth', 1.2); hold on;
plot(nu_list, PSR_col_comp, 's-', 'LineWidth', 1.2);
grid on; xlabel('Doppler index \nu'); ylabel('PSR\_col / dB');
title('PSR\_col（真值列口径）：Raw vs Per-bin Fast-time Compensation');
legend('Raw', 'Comp', 'Location', 'best');

%% ===================== 画 2D AFfs 对比（nu_show） =====================
if ~isempty(AFfs_raw_show)
    AF_raw_dB  = 10*log10(AFfs_raw_show  ./ max(AFfs_raw_show(:)));
    AF_comp_dB = 10*log10(AFfs_comp_show ./ max(AFfs_comp_show(:)));

    figure;
    imagesc(fftshift(AF_raw_dB, 2)); axis xy; colormap jet; colorbar; caxis([-30 0]);
    xlabel('Doppler bin (fftshift)'); ylabel('Delay bin');
    title(['AFfs RAW (nu = ', num2str(nu_show), ')']);

    figure;
    imagesc(fftshift(AF_comp_dB, 2)); axis xy; colormap jet; colorbar; caxis([-30 0]);
    xlabel('Doppler bin (fftshift)'); ylabel('Delay bin');
    title(['AFfs COMP (nu = ', num2str(nu_show), ')']);

    figure;
    imagesc(fftshift(AF_comp_dB - AF_raw_dB, 2)); axis xy; colormap jet; colorbar;
    xlabel('Doppler bin (fftshift)'); ylabel('Delay bin');
    title(['AFfs improvement (COMP - RAW) dB (nu = ', num2str(nu_show), ')']);
end

%% ===================== 本脚本用到的本地函数 =====================

function [Y, Xtd, Xf] = simulate_ofdm_echo_multi(N,M,tau_vec,nu_vec,flag,K,amp)
    tau_vec = tau_vec(:).';
    nu_vec  = nu_vec(:).';
    amp     = amp(:).';

    P = max([numel(tau_vec), numel(nu_vec), numel(amp)]);
    if numel(tau_vec)==1, tau_vec = repmat(tau_vec,1,P); end
    if numel(nu_vec)==1,  nu_vec  = repmat(nu_vec, 1,P); end
    if numel(amp)==1,     amp     = repmat(amp,    1,P); end

    tau_vec = mod(round(tau_vec), N);
    nu_vec  = mod(round(nu_vec),  M);   % 这里把 Doppler index 做了 mod M（环形）

    s = randi(K,[N,M]) - 1;
    if flag == 1
        Xf = pskmod(s, K);
    else
        Xf = qammod(s, K);
    end

    Xtd = sqrt(N) * ifft(Xf, [], 1);

    k = (0:N-1).';
    m = 0:M-1;

    Y = complex(zeros(N,M));
    for p = 1:P
        tau_p = tau_vec(p);
        nu_p  = nu_vec(p);

        X_shift = circshift(Xtd, tau_p, 1);

        % 你原模型：快时间 dop_fast + 慢时间 dop_slow
        dop_fast = exp(1j*2*pi*nu_p * k/(M*N));
        dop_slow = exp(1j*2*pi*nu_p * m/M);
        Dop = dop_fast * dop_slow;

        Y = Y + amp(p) * (X_shift .* Dop);
    end
end

你关心的“补偿应该加在哪一行？”

就在脚本里这段的位置（与 DAFT 的“先去相位再 FFT”完全一致）：

matlab
[Y, ~, Xf] = simulate_ofdm_echo_multi(...);

% ——就在这里对 Y 做快时间去旋转（de-rotation）——
% 然后再做 Yf = fft(Y,[],1)
Yf = fft(Y, [], 1);

如果你运行后仍然发现 PSR_full 变差、但 PSR_col 变好，这是正常的：逐 bin 补偿本质是“多假设匹配”，会把很多“错误 bin”也抬起来，从而拉大全图分母；但对“真值 bin”会更聚焦、更能体现 ICI 被压下去。你要我把指标改成 PSLR/ISLR（更常用于 AF 旁瓣评价）也可以。

問題

不用imagesc来画图 而是用surf来对比

回答

下面给你一份完整脚本 .m（不是函数），和我上一版逻辑一致，只是把 2D 图对比的绘图全部改成 surf（Raw / Comp / 差值），并且用了 fftshift 把 Doppler 维居中显示。直接复制保存为 demo_AFfs_fastDopplerComp_surf.m 运行即可。

matlab
clc; clear; close all;
format long;

%% ===================== 参数区 =====================
N   = 32;        % fast-time samples per pulse
M   = 16;        % slow-time pulses
Nmc = 10;        % Monte-Carlo runs

tau0 = 4;        % delay
flag = 1;        % 1: PSK, else: QAM
Kmod = 16;       % modulation order
alpha = 1;       % target amplitude (single target)

nu_list = -M:M;  % 扫描的 Doppler index（注意：仿真里 nu 会 mod M）
nu_show = 8;     % 用于画 2D surf 对比的 Doppler index（必须在 nu_list 里）

% 快时间/慢时间索引
k_fast = (0:N-1).';
m_slow = (0:M-1).';

% slow-time FFT 的 steering（与 fft(R,[],2) 的相位定义一致）
W = exp(-1j * 2*pi/M * (m_slow * (0:M-1)));  % M×M

% 记录 PSR：两种口径（全图/真列）
PSR_full_raw  = zeros(numel(nu_list),1);
PSR_full_comp = zeros(numel(nu_list),1);
PSR_col_raw   = zeros(numel(nu_list),1);
PSR_col_comp  = zeros(numel(nu_list),1);

% 保存用于展示的 2D AFfs（nu_show）
AFfs_raw_show  = [];
AFfs_comp_show = [];

%% ===================== 主循环：扫 Doppler =====================
rng(1545,'twister');                       % 固定随机性（可复现）
seeds = randi(2^31-1, Nmc, 1, 'uint32');   % parfor 可复现

for idx = 1:numel(nu_list)
    nu0 = nu_list(idx);

    AF_sum_raw  = zeros(N,M);
    AF_sum_comp = zeros(N,M);

    parfor mc = 1:Nmc
        rng(double(seeds(mc)),'twister');

        % ---------- 生成回波（N-structure：Y 是 N×M） ----------
        [Y, ~, Xf] = simulate_ofdm_echo_multi(N, M, tau0, nu0, flag, Kmod, alpha);

        %% ============ 方案 A：RAW（标准2D-FFT） ============
        Yf  = fft(Y, [], 1);                  % 快时间 FFT
        R   = ifft(Yf .* conj(Xf), [], 1);    % 匹配滤波 + 快时间 IFFT（range 压缩）
        RDM = fft(R, [], 2);                  % 慢时间 FFT（Doppler）
        AF_trial_raw = abs(RDM).^2;

        %% ============ 方案 B：逐 Doppler bin 快时间补偿 ============
        % 补偿位置：在快时间 FFT 之前，对 Y 乘 e^{-j2π q k/(MN)}
        AF_trial_comp = zeros(N,M);

        for q = 0:M-1
            rot_fast = exp(-1j * 2*pi * q * k_fast/(M*N));  % N×1
            Yc = bsxfun(@times, Y, rot_fast);               % N×M

            Yf_c = fft(Yc, [], 1);
            R_c  = ifft(Yf_c .* conj(Xf), [], 1);

            % 只取 Doppler bin = q（等价于 fft(R_c,[],2) 的第 q+1 列）
            AF_trial_comp(:, q+1) = abs(R_c * W(:, q+1)).^2;
        end

        AF_sum_raw  = AF_sum_raw  + AF_trial_raw;
        AF_sum_comp = AF_sum_comp + AF_trial_comp;
    end

    AFfs_raw  = AF_sum_raw  / Nmc;
    AFfs_comp = AF_sum_comp / Nmc;

    %% ===================== PSR 计算（两种口径） =====================
    % (1) 全图 PSR
    peak_raw  = max(AFfs_raw(:));
    peak_comp = max(AFfs_comp(:));
    PSR_full_raw(idx)  = 10*log10( peak_raw  / (sum(AFfs_raw(:))  - peak_raw )  * (M*N+1) );
    PSR_full_comp(idx) = 10*log10( peak_comp / (sum(AFfs_comp(:)) - peak_comp) * (M*N+1) );

    % (2) 真值 Doppler 对应列 PSR（更能反映 ICI 抑制）
    q_true = mod(round(nu0), M);          % 0..M-1（因为仿真里 nu 会 mod M）
    col_raw  = AFfs_raw(:,  q_true+1);
    col_comp = AFfs_comp(:, q_true+1);
    pk1 = max(col_raw);
    pk2 = max(col_comp);
    PSR_col_raw(idx)  = 10*log10( pk1 / (sum(col_raw)  - pk1) * (N+1) );
    PSR_col_comp(idx) = 10*log10( pk2 / (sum(col_comp) - pk2) * (N+1) );

    %% 保存展示用 2D 图
    if nu0 == nu_show
        AFfs_raw_show  = AFfs_raw;
        AFfs_comp_show = AFfs_comp;
    end
end

%% ===================== PSR 曲线对比 =====================
figure;
plot(nu_list, PSR_full_raw,  'o-', 'LineWidth', 1.2); hold on;
plot(nu_list, PSR_full_comp, 's-', 'LineWidth', 1.2);
grid on; xlabel('Doppler index \nu'); ylabel('PSR\_full / dB');
title('PSR\_full (全图口径)：Raw vs Per-bin Fast-time Compensation');
legend('Raw', 'Comp', 'Location', 'best');

figure;
plot(nu_list, PSR_col_raw,  'o-', 'LineWidth', 1.2); hold on;
plot(nu_list, PSR_col_comp, 's-', 'LineWidth', 1.2);
grid on; xlabel('Doppler index \nu'); ylabel('PSR\_col / dB');
title('PSR\_col（真值列口径）：Raw vs Per-bin Fast-time Compensation');
legend('Raw', 'Comp', 'Location', 'best');

%% ===================== 用 surf 画 2D AFfs 对比（nu_show） =====================
if ~isempty(AFfs_raw_show)
    % 归一化到各自峰值并转 dB
    AF_raw_dB  = 10*log10(AFfs_raw_show  ./ max(AFfs_raw_show(:)));
    AF_comp_dB = 10*log10(AFfs_comp_show ./ max(AFfs_comp_show(:)));
    AF_impr_dB = AF_comp_dB - AF_raw_dB;

    % Doppler 维居中显示（只 shift 第2维）
    AF_raw_dB_s  = fftshift(AF_raw_dB,  2);
    AF_comp_dB_s = fftshift(AF_comp_dB, 2);
    AF_impr_dB_s = fftshift(AF_impr_dB, 2);

    % 轴
    delay_axis   = 0:N-1;
    doppler_axis = (-floor(M/2)):(ceil(M/2)-1);  % 对应 fftshift 后的列顺序
    [DOP, DEL]   = meshgrid(doppler_axis, delay_axis);

    % --- RAW surf ---
    figure;
    surf(DOP, DEL, AF_raw_dB_s);
    shading interp; colorbar;
    xlabel('Doppler bin (shifted)');
    ylabel('Delay bin');
    zlabel('Normalized Power / dB');
    title(['AFfs RAW (surf)  nu = ', num2str(nu_show)]);
    zlim([-40 0]);  % 你想看更深可改

    % --- COMP surf ---
    figure;
    surf(DOP, DEL, AF_comp_dB_s);
    shading interp; colorbar;
    xlabel('Doppler bin (shifted)');
    ylabel('Delay bin');
    zlabel('Normalized Power / dB');
    title(['AFfs COMP (surf)  nu = ', num2str(nu_show)]);
    zlim([-40 0]);

    % --- Improvement surf ---
    figure;
    surf(DOP, DEL, AF_impr_dB_s);
    shading interp; colorbar;
    xlabel('Doppler bin (shifted)');
    ylabel('Delay bin');
    zlabel('Delta / dB');
    title(['AFfs Improvement (COMP - RAW) (surf)  nu = ', num2str(nu_show)]);
end

%% ===================== 本脚本用到的本地函数 =====================
function [Y, Xtd, Xf] = simulate_ofdm_echo_multi(N,M,tau_vec,nu_vec,flag,K,amp)
    tau_vec = tau_vec(:).';
    nu_vec  = nu_vec(:).';
    amp     = amp(:).';

    P = max([numel(tau_vec), numel(nu_vec), numel(amp)]);
    if numel(tau_vec)==1, tau_vec = repmat(tau_vec,1,P); end
    if numel(nu_vec)==1,  nu_vec  = repmat(nu_vec, 1,P); end
    if numel(amp)==1,     amp     = repmat(amp,    1,P); end

    tau_vec = mod(round(tau_vec), N);
    nu_vec  = mod(round(nu_vec),  M);   % Doppler index 做 mod M（环形）

    s = randi(K,[N,M]) - 1;
    if flag == 1
        Xf = pskmod(s, K);
    else
        Xf = qammod(s, K);
    end

    Xtd = sqrt(N) * ifft(Xf, [], 1);

    k = (0:N-1).';
    m = 0:M-1;

    Y = complex(zeros(N,M));
    for p = 1:P
        tau_p = tau_vec(p);
        nu_p  = nu_vec(p);

        X_shift = circshift(Xtd, tau_p, 1);

        dop_fast = exp(1j*2*pi*nu_p * k/(M*N));
        dop_slow = exp(1j*2*pi*nu_p * m/M);
        Dop = dop_fast * dop_slow;

        Y = Y + amp(p) * (X_shift .* Dop);
    end
end

如果你还想把 surf 的视角（view）、或者把 z 轴改成线性幅度而不是 dB、或者把 zlim 做成自动/一致（raw 和 comp 用同一个峰值归一化）我也可以顺手帮你改。

問題

image.pngimage.png
image.pngimage.png

我现在设置N = 64; % fast-time samples per pulse
M = 64; % slow-time pulses
Nmc = 100; % Monte-Carlo runs
可是当我目标doppler放在M/2的位置上，显示出来的结果在边上，当目标doppler最大的时候，目标又跑到相应图的中央了。分析一下是什么原因并改正

思考

Clarifying Doppler axis mapping after fftshift

I’m working through the mapping of Doppler bins after fftshift. For M=64, the bins range from -32 to 31, but the shift causes the ordering to adjust. After fftshift, the center bin corresponds to the maximum Doppler, while the edges show the target Doppler at M/2. I’ll need to validate the modulo mapping that’s also likely involved here.

Clarifying Doppler center alignment

When nu0 is set to M/2 (32), it wraps to 32, which is the Nyquist frequency and appears at the edge after fftshift. This is expected for even M. However, the user expects M/2 to appear at the center. To fix this, I'll standardize the representation to a centered range [-M/2, M/2-1], removing any mod issues.

Addressing user confusion with axis labels after fftshift

The user might expect nu=32 to be at the center, but in the current setup, fftshift places the Nyquist frequency (bin 32) at the edge. For a centered axis, nu=32 could correspond to -32 (the Nyquist). To place nu=32 at the center, we may need to adjust the shifting or labeling convention.

Resolving alignment issues with Doppler frequency plot

It seems the user wants the target Doppler at M/2 to appear at the center of an unshifted axis (0..M-1) instead of as a frequency. When they set nu=M-1 (maximum Doppler), it ends up at the center due to fftshift. To address this, I’ll either avoid fftshift or adjust it based on the user’s expectations.

Addressing Doppler bin labeling and plotting issues

I’ll explain to the user that when using fftshift, bin M/2 (Nyquist) maps to the edge, and bin M-1 maps to -1, which causes it to appear near the center. I’ll propose two options: keeping the axis as 0..M-1 without shift or adjusting the nu mapping for centered representation. I’ll also correct the figure labels with proper nu values.

Thought for 1m 2s

回答

你这个现象不是算法坏了，而是Doppler bin 的“循环频率 + fftshift 显示方式”导致的，尤其 M 为偶数时 M/2 那个 bin 是 Nyquist 点，在 fftshift 后必然跑到边界。

1) 为什么 “nu = M/2 在边上、nu = M-1 在中间”？

你仿真里 Doppler 被写成了 循环索引：

matlab
nu_vec  = mod(round(nu_vec),  M);   % 0..M-1

而 fft(R,[],2) 的 Doppler bin 本来就是 0..M-1。
你画图时又做了：

matlab
AF = fftshift(AF, 2);
doppler_axis = -M/2 : M/2-1;

对于 M=64：

  • q = 63（你说的“最大”）在“居中坐标”里等价于 -1
    → fftshift 后会出现在 0 附近（中间），所以你看到它“跑到中央”。
  • q = 32 = M/2 是 Nyquist bin（等价于 -32）
    → 对偶数长度 fftshift 会把 Nyquist 分量放到第一列（最边上），所以你看到它“跑到边上”。

这完全是 FFT/fftshift 的标准行为。

2) 你要“nu=0..M-1 按索引显示在 0..M-1 上”，怎么改？

最简单：画图不要 fftshift，横轴直接用 0:M-1。
这样：

  • nu = 32 就会在横轴中间附近
  • nu = 63 就会在横轴末端（不会跑回中间）

把你脚本里 surf 画图那一段替换成下面这个（RAW/COMP/提升都改好了）：

matlab
% ===================== 用 surf 画 2D AFfs 对比（不 fftshift） =====================
if ~isempty(AFfs_raw_show)
    % 归一化到各自峰值并转 dB
    AF_raw_dB  = 10*log10(AFfs_raw_show  ./ max(AFfs_raw_show(:)));
    AF_comp_dB = 10*log10(AFfs_comp_show ./ max(AFfs_comp_show(:)));
    AF_impr_dB = AF_comp_dB - AF_raw_dB;

    % 轴（注意：这里是 0..M-1，不做 shift）
    delay_axis   = 0:N-1;
    doppler_axis = 0:M-1;
    [DOP, DEL]   = meshgrid(doppler_axis, delay_axis);

    figure;
    surf(DOP, DEL, AF_raw_dB);
    shading interp; colorbar;
    xlabel('Doppler bin (0..M-1)');
    ylabel('Delay bin');
    zlabel('Normalized Power / dB');
    title(['AFfs RAW (surf)  nu\_in = ', num2str(nu_show), ...
           ', nu\_bin = ', num2str(mod(round(nu_show),M))]);
    zlim([-40 0]);

    figure;
    surf(DOP, DEL, AF_comp_dB);
    shading interp; colorbar;
    xlabel('Doppler bin (0..M-1)');
    ylabel('Delay bin');
    zlabel('Normalized Power / dB');
    title(['AFfs COMP (surf)  nu\_in = ', num2str(nu_show), ...
           ', nu\_bin = ', num2str(mod(round(nu_show),M))]);
    zlim([-40 0]);

    figure;
    surf(DOP, DEL, AF_impr_dB);
    shading interp; colorbar;
    xlabel('Doppler bin (0..M-1)');
    ylabel('Delay bin');
    zlabel('\Delta / dB');
    title(['AFfs Improvement (COMP-RAW) (surf)  nu\_in = ', num2str(nu_show)]);
end

3) 如果你坚持要用 fftshift 的“物理频率轴 [-M/2, M/2-1]”怎么办？

那就要接受一个事实：偶数点 FFT 的 Nyquist（±M/2）在 shift 后就是在边上。
也就是说：nu = M/2 显示在边界是“正确的”，不是错。

你要我可以给你同时画两张：

  • 一张 0..M-1（工程 bin 索引）
  • 一张 fftshift（正负频率）
    这样就不会混淆“索引”和“频率”。

你现在更偏向哪种坐标系：0..M-1（按 bin）还是 [-M/2..M/2-1]（按正负频率）？我可以把整份脚本统一成你要的那一种。

問題

image.pngimage.png

刚才这个没啥问题了，再看看真值列，结果好像有点一般,因为真值的补偿的感觉很好了。真值列要不我们扩大一个口径，比如周边几列这种？然后同时在下面的代码里面再加上zero delay 和zero doppler的slice，以及实际值的delay和doppler的slice：clc; clear; close all;
format long;

%% ===================== 参数区 =====================
N = 64; % fast-time samples per pulse
M = 64; % slow-time pulses
Nmc = 100; % Monte-Carlo runs

tau0 = 4; % delay
flag = 1; % 1: PSK, else: QAM
Kmod = 16; % modulation order
alpha = 1; % target amplitude (single target)

nu_list = -M:M; % 扫描的 Doppler index（注意：仿真里 nu 会 mod M）
nu_show = M-1; % 用于画 2D surf 对比的 Doppler index（必须在 nu_list 里）

% 快时间/慢时间索引
k_fast = (0:N-1).';
m_slow = (0:M-1).';

% slow-time FFT 的 steering（与 fft(R,[],2) 的相位定义一致）
W = exp(-1j * 2*pi/M * (m_slow * (0:M-1))); % M×M

% 记录 PSR：两种口径（全图/真列）
PSR_full_raw = zeros(numel(nu_list),1);
PSR_full_comp = zeros(numel(nu_list),1);
PSR_col_raw = zeros(numel(nu_list),1);
PSR_col_comp = zeros(numel(nu_list),1);

% 保存用于展示的 2D AFfs（nu_show）
AFfs_raw_show = [];
AFfs_comp_show = [];

%% ===================== 主循环：扫 Doppler =====================
rng(1545,'twister'); % 固定随机性（可复现）
seeds = randi(2^31-1, Nmc, 1, 'uint32'); % parfor 可复现

for idx = 1:numel(nu_list)
nu0 = nu_list(idx);

text
AF_sum_raw  = zeros(N,M);
AF_sum_comp = zeros(N,M);

parfor mc = 1:Nmc
    rng(double(seeds(mc)),'twister');

    % ---------- 生成回波（N-structure：Y 是 N×M） ----------
    [Y, ~, Xf] = simulate_ofdm_echo_multi(N, M, tau0, nu0, flag, Kmod, alpha);

    %% ============ 方案 A：RAW（标准2D-FFT） ============
    Yf  = fft(Y, [], 1);                  % 快时间 FFT
    R   = ifft(Yf .* conj(Xf), [], 1);    % 匹配滤波 + 快时间 IFFT（range 压缩）
    RDM = fft(R, [], 2);                  % 慢时间 FFT（Doppler）
    AF_trial_raw = abs(RDM).^2;

    %% ============ 方案 B：逐 Doppler bin 快时间补偿 ============
    % 补偿位置：在快时间 FFT 之前，对 Y 乘 e^{-j2π q k/(MN)}
    AF_trial_comp = zeros(N,M);

    for q = 0:M-1
        rot_fast = exp(-1j * 2*pi * q * k_fast/(M*N));  % N×1
        Yc = bsxfun(@times, Y, rot_fast);               % N×M

        Yf_c = fft(Yc, [], 1);
        R_c  = ifft(Yf_c .* conj(Xf), [], 1);

        % 只取 Doppler bin = q（等价于 fft(R_c,[],2) 的第 q+1 列）
        AF_trial_comp(:, q+1) = abs(R_c * W(:, q+1)).^2;
    end

    AF_sum_raw  = AF_sum_raw  + AF_trial_raw;
    AF_sum_comp = AF_sum_comp + AF_trial_comp;
end

AFfs_raw  = AF_sum_raw  / Nmc;
AFfs_comp = AF_sum_comp / Nmc;

%% ===================== PSR 计算（两种口径） =====================
% (1) 全图 PSR
peak_raw  = max(AFfs_raw(:));
peak_comp = max(AFfs_comp(:));
PSR_full_raw(idx)  = 10*log10( peak_raw  / (sum(AFfs_raw(:))  - peak_raw )  * (M*N+1) );
PSR_full_comp(idx) = 10*log10( peak_comp / (sum(AFfs_comp(:)) - peak_comp) * (M*N+1) );

% (2) 真值 Doppler 对应列 PSR（更能反映 ICI 抑制）
q_true = mod(round(nu0), M);          % 0..M-1（因为仿真里 nu 会 mod M）
col_raw  = AFfs_raw(:,  q_true+1);
col_comp = AFfs_comp(:, q_true+1);
pk1 = max(col_raw);
pk2 = max(col_comp);
PSR_col_raw(idx)  = 10*log10( pk1 / (sum(col_raw)  - pk1) * (N+1) );
PSR_col_comp(idx) = 10*log10( pk2 / (sum(col_comp) - pk2) * (N+1) );

%% 保存展示用 2D 图
if nu0 == nu_show
    AFfs_raw_show  = AFfs_raw;
    AFfs_comp_show = AFfs_comp;
end

end

%% ===================== PSR 曲线对比 =====================
figure;
plot(nu_list, PSR_full_raw, 'o-', 'LineWidth', 1.2); hold on;
plot(nu_list, PSR_full_comp, 's-', 'LineWidth', 1.2);
grid on; xlabel('Doppler index \nu'); ylabel('PSR_full / dB');
title('PSR_full (全图口径)：Raw vs Per-bin Fast-time Compensation');
legend('Raw', 'Comp', 'Location', 'best');

figure;
plot(nu_list, PSR_col_raw, 'o-', 'LineWidth', 1.2); hold on;
plot(nu_list, PSR_col_comp, 's-', 'LineWidth', 1.2);
grid on; xlabel('Doppler index \nu'); ylabel('PSR_col / dB');
title('PSR_col（真值列口径）：Raw vs Per-bin Fast-time Compensation');
legend('Raw', 'Comp', 'Location', 'best');

%% ===================== 用 surf 画 2D AFfs 对比（nu_show） =====================
% ===================== 用 surf 画 2D AFfs 对比（不 fftshift） =====================
if ~isempty(AFfs_raw_show)
% 归一化到各自峰值并转 dB
AF_raw_dB = 10log10(AFfs_raw_show ./ max(AFfs_raw_show(:)));
AF_comp_dB = 10log10(AFfs_comp_show ./ max(AFfs_comp_show(:)));
AF_impr_dB = AF_comp_dB - AF_raw_dB;

text
% 轴（注意：这里是 0..M-1，不做 shift）
delay_axis   = 0:N-1;
doppler_axis = 0:M-1;
[DOP, DEL]   = meshgrid(doppler_axis, delay_axis);

figure;
surf(DOP, DEL, AF_raw_dB);
shading interp; colorbar;
xlabel('Doppler bin (0..M-1)');
ylabel('Delay bin');
zlabel('Normalized Power / dB');
title(['AFfs RAW (surf)  nu\_in = ', num2str(nu_show), ...
       ', nu\_bin = ', num2str(mod(round(nu_show),M))]);
%zlim([-60 0]);

figure;
surf(DOP, DEL, AF_comp_dB);
shading interp; colorbar;
xlabel('Doppler bin (0..M-1)');
ylabel('Delay bin');
zlabel('Normalized Power / dB');
title(['AFfs COMP (surf)  nu\_in = ', num2str(nu_show), ...
       ', nu\_bin = ', num2str(mod(round(nu_show),M))]);
%zlim([-60 0]);

figure;
surf(DOP, DEL, AF_impr_dB);
shading interp; colorbar;
xlabel('Doppler bin (0..M-1)');
ylabel('Delay bin');
zlabel('\Delta / dB');
title(['AFfs Improvement (COMP-RAW) (surf)  nu\_in = ', num2str(nu_show)]);

end

%% ===================== 本脚本用到的本地函数 =====================
function [Y, Xtd, Xf] = simulate_ofdm_echo_multi(N,M,tau_vec,nu_vec,flag,K,amp)
tau_vec = tau_vec(:).';
nu_vec = nu_vec(:).';
amp = amp(:).';

text
P = max([numel(tau_vec), numel(nu_vec), numel(amp)]);
if numel(tau_vec)==1, tau_vec = repmat(tau_vec,1,P); end
if numel(nu_vec)==1,  nu_vec  = repmat(nu_vec, 1,P); end
if numel(amp)==1,     amp     = repmat(amp,    1,P); end

tau_vec = mod(round(tau_vec), N);
nu_vec  = mod(round(nu_vec),  M);   % Doppler index 做 mod M（环形）

s = randi(K,[N,M]) - 1;
if flag == 1
    Xf = pskmod(s, K);
else
    Xf = qammod(s, K);
end

Xtd = sqrt(N) * ifft(Xf, [], 1);

k = (0:N-1).';
m = 0:M-1;

Y = complex(zeros(N,M));
for p = 1:P
    tau_p = tau_vec(p);
    nu_p  = nu_vec(p);

    X_shift = circshift(Xtd, tau_p, 1);

    dop_fast = exp(1j*2*pi*nu_p * k/(M*N));
    dop_slow = exp(1j*2*pi*nu_p * m/M);
    Dop = dop_fast * dop_slow;

    Y = Y + amp(p) * (X_shift .* Dop);
end

end

思考

Updating code for expanded aperture and slices

Looks like the goal is to modify the code to expand the aperture for the true value column by including neighboring Doppler bins and delay values. They want to compute PSR not just for a single Doppler but across a band of columns. I should define half-width for both Doppler and delay and ensure zero-delay and zero-Doppler slices are correctly implemented. I'll need to clarify the slicing logic too.

Plotting slices and expanding Doppler aperture

For zero Doppler, bin 0 corresponds to Doppler = 0. I’ll generate plots for both raw and computed values based on the delay slice (along Doppler) at tau0 and the Doppler slice (along delay) at q_true, also considering surf plots. Additionally, I’ll calculate the aperture PSR by summing columns in the Doppler neighborhood.

Implementing banded PSR and slices

I'll implement "banded column PSR" (PSR_band) by summing selected columns and adding PSR_band arrays. I also need to save AFfs matrices for nu_show to plot slices. Since nu_show is a Doppler input, I’ll calculate q_true_show and save for slices. I'll focus on slices for nu_show, using AFfs_raw_show and comp for plotting. Additionally, I’ll introduce Bdop as a parameter with a value of 1 or 2.

Calculating PSR_band in loop

After computing AFfs, I'll calculate PSR_band_raw and PSR_band_comp within the loop. First, I'll determine the Doppler bin q_true and its surrounding band based on Bdop. Then, I'll sum the raw and computed AFfs values across the selected Doppler columns. The resulting PSR is calculated using the maximum of the summed values, with the possibility of using averaging instead of summing, which wouldn’t affect PSR's scaling.

Computing PSR_band2D with slices and plots

I’ll compute PSR_band2D using a subset of bins and then create plots of nu_list vs. PSR_band. For slices, I’ll include zero Doppler and Delay slices, along with actual slices based on tau0 and q_true_show. For Doppler, I’ll consider whether to use a centered or unshifted axis, keeping it consistent to avoid confusion. I'll use fftshift for signed axis but will keep it minimal, as required.

Thought for 57s

回答

你现在的 **PSR_col（只取真值 Doppler 那一列）**之所以“看起来一般”，最常见原因是：

  • 即使补偿后 ICI 变小，目标能量仍可能分布在真值列附近的几列（例如 ±1、±2 个 Doppler bin），尤其当 Doppler 接近 Nyquist、或者存在离栅（非整数 bin）/有限长度泄漏时。
  • 你只看一列会“漏掉”这些能量 → 指标提升被低估。

所以可以把“真值列口径”扩大成 真值 Doppler 附近 ±B 列的能量合成（band integration），再算一个 PSR_band，更符合你想观察的“补偿后能量回聚”。

另外你要的切片我也加上了（Raw vs Comp 同图对比）：

  • Zero Doppler slice：Doppler=0 列（沿 delay）
  • Zero Delay slice：Delay=0 行（沿 doppler）
  • True Doppler slice：真值 Doppler 列（沿 delay）
  • True Delay slice：真值 delay 行（沿 doppler）

下面是你给的脚本基础上直接改好的完整 .m 文件（仍然用 surf 画 2D，对切片用 plot）：

matlab
clc; clear; close all;
format long;

%% ===================== 参数区 =====================
N   = 64;        % fast-time samples per pulse
M   = 64;        % slow-time pulses
Nmc = 100;       % Monte-Carlo runs

tau0 = 4;        % true delay (samples)
flag = 1;        % 1: PSK, else: QAM
Kmod = 16;       % modulation order
alpha = 1;       % target amplitude (single target)

nu_list = -M:M;     % scan input Doppler index (will mod M inside sim)
nu_show = M-1;      % choose one nu for 2D/slices visualization

Bdop = 2;           % >>> 真值 Doppler 周边口径：±Bdop 列合成（建议 1~3）

% 索引
k_fast = (0:N-1).';
m_slow = (0:M-1).';

% slow-time FFT steering (consistent with fft(.,[],2))
W = exp(-1j * 2*pi/M * (m_slow * (0:M-1)));  % M×M, column q+1 is bin q

% 记录 PSR：全图、真值列、真值邻域列合成
PSR_full_raw   = zeros(numel(nu_list),1);
PSR_full_comp  = zeros(numel(nu_list),1);
PSR_col_raw    = zeros(numel(nu_list),1);
PSR_col_comp   = zeros(numel(nu_list),1);
PSR_band_raw   = zeros(numel(nu_list),1);
PSR_band_comp  = zeros(numel(nu_list),1);

% 保存用于展示的 2D AFfs（nu_show）
AFfs_raw_show  = [];
AFfs_comp_show = [];
q_true_show    = [];

%% ===================== 主循环：扫 Doppler =====================
rng(1545,'twister');                       % reproducible
seeds = randi(2^31-1, Nmc, 1, 'uint32');   % parfor reproducible

for idx = 1:numel(nu_list)
    nu0 = nu_list(idx);

    AF_sum_raw  = zeros(N,M);
    AF_sum_comp = zeros(N,M);

    parfor mc = 1:Nmc
        rng(double(seeds(mc)),'twister');

        % ---------- 生成回波（N-structure：Y 是 N×M） ----------
        [Y, ~, Xf] = simulate_ofdm_echo_multi(N, M, tau0, nu0, flag, Kmod, alpha);

        %% ============ 方案 A：RAW（标准2D-FFT） ============
        Yf  = fft(Y, [], 1);                  % fast-time FFT
        R   = ifft(Yf .* conj(Xf), [], 1);    % match + IFFT (range)
        RDM = fft(R, [], 2);                  % slow-time FFT (Doppler)
        AF_trial_raw = abs(RDM).^2;

        %% ============ 方案 B：逐 Doppler bin 快时间补偿 ============
        % 补偿位置：在快时间 FFT 之前，对 Y 乘 e^{-j2π q k/(MN)}
        AF_trial_comp = zeros(N,M);

        for q = 0:M-1
            rot_fast = exp(-1j * 2*pi * q * k_fast/(M*N));  % N×1
            Yc = bsxfun(@times, Y, rot_fast);               % N×M

            Yf_c = fft(Yc, [], 1);
            R_c  = ifft(Yf_c .* conj(Xf), [], 1);

            % 只取 Doppler bin = q
            AF_trial_comp(:, q+1) = abs(R_c * W(:, q+1)).^2;
        end

        AF_sum_raw  = AF_sum_raw  + AF_trial_raw;
        AF_sum_comp = AF_sum_comp + AF_trial_comp;
    end

    AFfs_raw  = AF_sum_raw  / Nmc;
    AFfs_comp = AF_sum_comp / Nmc;

    %% ===================== PSR 计算 =====================
    % (1) 全图 PSR
    peak_raw  = max(AFfs_raw(:));
    peak_comp = max(AFfs_comp(:));
    PSR_full_raw(idx)  = 10*log10( peak_raw  / (sum(AFfs_raw(:))  - peak_raw )  * (M*N+1) );
    PSR_full_comp(idx) = 10*log10( peak_comp / (sum(AFfs_comp(:)) - peak_comp) * (M*N+1) );

    % 真值 Doppler bin（因为仿真内部 nu 做了 mod M）
    q_true = mod(round(nu0), M);   % 0..M-1

    % (2) 真值列口径 PSR_col（只看 q_true 那一列）
    col_raw  = AFfs_raw(:,  q_true+1);
    col_comp = AFfs_comp(:, q_true+1);
    pk1 = max(col_raw);  pk2 = max(col_comp);
    PSR_col_raw(idx)  = 10*log10( pk1 / (sum(col_raw)  - pk1) * (N+1) );
    PSR_col_comp(idx) = 10*log10( pk2 / (sum(col_comp) - pk2) * (N+1) );

    % (3) 扩大口径：真值 Doppler 周边 ±Bdop 列合成（band integration）
    q_band = mod((q_true-Bdop):(q_true+Bdop), M);     % wrap-around
    cols   = q_band + 1;

    band_raw  = sum(AFfs_raw(:,  cols), 2);           % N×1 (delay profile)
    band_comp = sum(AFfs_comp(:, cols), 2);

    pkb1 = max(band_raw);  pkb2 = max(band_comp);
    PSR_band_raw(idx)  = 10*log10( pkb1 / (sum(band_raw)  - pkb1) * (N+1) );
    PSR_band_comp(idx) = 10*log10( pkb2 / (sum(band_comp) - pkb2) * (N+1) );

    %% 保存展示用 2D 图
    if nu0 == nu_show
        AFfs_raw_show  = AFfs_raw;
        AFfs_comp_show = AFfs_comp;
        q_true_show    = q_true;
    end
end

%% ===================== PSR 曲线对比 =====================
figure;
plot(nu_list, PSR_full_raw,  'o-', 'LineWidth', 1.2); hold on;
plot(nu_list, PSR_full_comp, 's-', 'LineWidth', 1.2);
grid on; xlabel('Doppler index \nu'); ylabel('PSR\_full / dB');
title('PSR\_full（全图口径）：Raw vs Comp');
legend('Raw', 'Comp', 'Location', 'best');

figure;
plot(nu_list, PSR_col_raw,  'o-', 'LineWidth', 1.2); hold on;
plot(nu_list, PSR_col_comp, 's-', 'LineWidth', 1.2);
grid on; xlabel('Doppler index \nu'); ylabel('PSR\_col / dB');
title('PSR\_col（真值列口径）：Raw vs Comp');
legend('Raw', 'Comp', 'Location', 'best');

figure;
plot(nu_list, PSR_band_raw,  'o-', 'LineWidth', 1.2); hold on;
plot(nu_list, PSR_band_comp, 's-', 'LineWidth', 1.2);
grid on; xlabel('Doppler index \nu'); ylabel('PSR\_band / dB');
title(['PSR\_band（真值±', num2str(Bdop), '列合成口径）：Raw vs Comp']);
legend('Raw', 'Comp', 'Location', 'best');

%% ===================== 用 surf 画 2D AFfs 对比（nu_show） =====================
if ~isempty(AFfs_raw_show)
    % 归一化到各自峰值并转 dB
    AF_raw_dB  = 10*log10(AFfs_raw_show  ./ max(AFfs_raw_show(:)));
    AF_comp_dB = 10*log10(AFfs_comp_show ./ max(AFfs_comp_show(:)));
    AF_impr_dB = AF_comp_dB - AF_raw_dB;

    delay_axis   = 0:N-1;
    doppler_axis = 0:M-1;
    [DOP, DEL]   = meshgrid(doppler_axis, delay_axis);

    figure;
    surf(DOP, DEL, AF_raw_dB);
    shading interp; colorbar;
    xlabel('Doppler bin (0..M-1)'); ylabel('Delay bin'); zlabel('Norm Power / dB');
    title(['AFfs RAW (surf)  nu\_in=', num2str(nu_show), ', nu\_bin=', num2str(mod(round(nu_show),M))]);

    figure;
    surf(DOP, DEL, AF_comp_dB);
    shading interp; colorbar;
    xlabel('Doppler bin (0..M-1)'); ylabel('Delay bin'); zlabel('Norm Power / dB');
    title(['AFfs COMP (surf) nu\_in=', num2str(nu_show), ', nu\_bin=', num2str(mod(round(nu_show),M))]);

    figure;
    surf(DOP, DEL, AF_impr_dB);
    shading interp; colorbar;
    xlabel('Doppler bin (0..M-1)'); ylabel('Delay bin'); zlabel('\Delta / dB');
    title(['AFfs Improvement (COMP-RAW) (surf)  nu\_in=', num2str(nu_show)]);

    %% ===================== slices：zero / true =====================
    % 为了 slice 显示更直观：Doppler 维用 signed 轴（fftshift）
    dop_signed = (-floor(M/2)):(ceil(M/2)-1);

    % 工具：安全归一化转 dB
    to_dB = @(v) 10*log10( max(v, eps) ./ max(max(v), eps) );

    % ---- Zero Doppler slice (Doppler bin = 0, column 1), along delay ----
    zeroD_raw  = AFfs_raw_show(:, 1);
    zeroD_comp = AFfs_comp_show(:,1);

    figure;
    plot(delay_axis, to_dB(zeroD_raw),  'o-', 'LineWidth', 1.2); hold on;
    plot(delay_axis, to_dB(zeroD_comp), 's-', 'LineWidth', 1.2);
    grid on; xlabel('Delay bin'); ylabel('Norm Power / dB');
    title('Zero Doppler slice (bin=0): Raw vs Comp');
    legend('Raw','Comp','Location','best');

    % ---- True Doppler slice (column q_true_show+1), along delay ----
    trueD_raw  = AFfs_raw_show(:,  q_true_show+1);
    trueD_comp = AFfs_comp_show(:, q_true_show+1);

    figure;
    plot(delay_axis, to_dB(trueD_raw),  'o-', 'LineWidth', 1.2); hold on;
    plot(delay_axis, to_dB(trueD_comp), 's-', 'LineWidth', 1.2);
    grid on; xlabel('Delay bin'); ylabel('Norm Power / dB');
    title(['True Doppler slice (bin=', num2str(q_true_show), '): Raw vs Comp']);
    legend('Raw','Comp','Location','best');

    % ---- Zero Delay slice (Delay bin = 0, row 1), along doppler ----
    zeroTau_raw  = AFfs_raw_show(1, :);
    zeroTau_comp = AFfs_comp_show(1,:);

    figure;
    plot(dop_signed, to_dB(fftshift(zeroTau_raw)),  'o-', 'LineWidth', 1.2); hold on;
    plot(dop_signed, to_dB(fftshift(zeroTau_comp)), 's-', 'LineWidth', 1.2);
    grid on; xlabel('Doppler bin (signed, fftshift)'); ylabel('Norm Power / dB');
    title('Zero Delay slice (tau=0): Raw vs Comp');
    legend('Raw','Comp','Location','best');

    % ---- True Delay slice (Delay bin = tau0, row tau0+1), along doppler ----
    if tau0+1 <= N
        trueTau_raw  = AFfs_raw_show(tau0+1, :);
        trueTau_comp = AFfs_comp_show(tau0+1,:);

        figure;
        plot(dop_signed, to_dB(fftshift(trueTau_raw)),  'o-', 'LineWidth', 1.2); hold on;
        plot(dop_signed, to_dB(fftshift(trueTau_comp)), 's-', 'LineWidth', 1.2);
        grid on; xlabel('Doppler bin (signed, fftshift)'); ylabel('Norm Power / dB');
        title(['True Delay slice (tau=', num2str(tau0), '): Raw vs Comp']);
        legend('Raw','Comp','Location','best');
    end

    % ---- 扩大口径的 true Doppler band（±Bdop）也给一个 delay slice ----
    q_band_show = mod((q_true_show-Bdop):(q_true_show+Bdop), M);
    cols_show   = q_band_show + 1;
    band_raw_show  = sum(AFfs_raw_show(:,  cols_show), 2);
    band_comp_show = sum(AFfs_comp_show(:, cols_show), 2);

    figure;
    plot(delay_axis, to_dB(band_raw_show),  'o-', 'LineWidth', 1.2); hold on;
    plot(delay_axis, to_dB(band_comp_show), 's-', 'LineWidth', 1.2);
    grid on; xlabel('Delay bin'); ylabel('Norm Power / dB');
    title(['True Doppler band slice (bin ', num2str(q_true_show), ' ±', num2str(Bdop), '): Raw vs Comp']);
    legend('Raw','Comp','Location','best');
end

%% ===================== 本脚本用到的本地函数 =====================
function [Y, Xtd, Xf] = simulate_ofdm_echo_multi(N,M,tau_vec,nu_vec,flag,K,amp)
    tau_vec = tau_vec(:).';
    nu_vec  = nu_vec(:).';
    amp     = amp(:).';

    P = max([numel(tau_vec), numel(nu_vec), numel(amp)]);
    if numel(tau_vec)==1, tau_vec = repmat(tau_vec,1,P); end
    if numel(nu_vec)==1,  nu_vec  = repmat(nu_vec, 1,P); end
    if numel(amp)==1,     amp     = repmat(amp,    1,P); end

    tau_vec = mod(round(tau_vec), N);
    nu_vec  = mod(round(nu_vec),  M);   % Doppler index -> 0..M-1 cyclic

    s = randi(K,[N,M]) - 1;
    if flag == 1
        Xf = pskmod(s, K);
    else
        Xf = qammod(s, K);
    end

    Xtd = sqrt(N) * ifft(Xf, [], 1);

    k = (0:N-1).';
    m = 0:M-1;

    Y = complex(zeros(N,M));
    for p = 1:P
        tau_p = tau_vec(p);
        nu_p  = nu_vec(p);

        X_shift = circshift(Xtd, tau_p, 1);

        dop_fast = exp(1j*2*pi*nu_p * k/(M*N));
        dop_slow = exp(1j*2*pi*nu_p * m/M);
        Dop = dop_fast * dop_slow;

        Y = Y + amp(p) * (X_shift .* Dop);
    end
end

你可以先这么看效果

  1. 把 Bdop 从 0/1/2/3 试一下：
  • Bdop=0 就退化成你原来的 PSR_col
  • Bdop=1~2 往往能明显把“补偿后能量回聚”的提升体现出来
  1. 看 True Doppler band slice 那张 delay 曲线：
    补偿好时应该 主峰更尖、底更低。

如果你还想更“雷达味儿”的指标（比如 PSLR/ISLR，用保护单元排除主瓣后再算旁瓣），我也能按你现在的 peak/guard 逻辑把它加进去，会比 PSR 更稳定。

問題

我们的比较方法就不用PSR了 这个方法很奇怪。然后我们现在从解析解的角度，把这种每一个doppler都补偿的滤波器组的方法，和afdm的方法进行对比，从理论上分析一下哪一种对于ICI的抑制效果更好。下面是一个可以参考的思路，但是不一定对：
\documentclass{article}
\usepackage{amsmath,amssymb}
\usepackage{geometry}
\geometry{a4paper,margin=1in}

\title{ICI Suppression Performance Comparison: AFDM vs. Doppler-Compensated/Uncompensated OFDM}
\author{Wireless Communication Analysis}
\date{}

\begin{document}
\maketitle

\begin{abstract}
This paper presents an analytical derivation of Inter-Carrier Interference (ICI) for uncompensated OFDM, Doppler-compensated OFDM, and Affine Frequency Division Multiplexing (AFDM) over double-selective channels. The ICI power expressions are derived to quantify the suppression performance, revealing the inherent advantage of AFDM's slant compensation over OFDM's linear compensation.
\end{abstract}

\section{1. Assumptions and Channel Model}
To ensure a fair comparison, a unified double-selective channel model (frequency-selective and time-varying, i.e., Doppler spread + delay spread) is established with the following key assumptions:
\begin{enumerate}
\item The channel is a multipath fading channel with PP independent paths. For the pp-th path, the delay is τp\tau_p, Doppler shift is fd,pf_{d,p}, and complex gain is hph_p (hp∼CN(0,σh2)h_p \sim \mathcal{CN}(0,\sigma_h^2)). It satisfies τ0<τ1<⋯<τP−1≤Tcp\tau_0 < \tau_1 < \dots < \tau_{P-1} \leq T_{\text{cp}}, where TcpT_{\text{cp}} is the length of Cyclic Prefix (CP) or Chirp Periodic Prefix (CPP) to eliminate Inter-Symbol Interference (ISI).
\item The number of subcarriers is NN, symbol period is TsT_s, subcarrier spacing is Δf=1/Ts\Delta f = 1/T_s, and normalized Doppler shift is νp=fd,pTs\nu_p = f_{d,p} T_s (characterizing the Doppler impact on symbol period).
\item Doppler-compensated OFDM can only perfectly compensate for integer multiple Doppler shifts, leaving a residual fractional Doppler shift ϵp=νp−⌊νp⌋\epsilon_p = \nu_p - \lfloor \nu_p \rfloor (0≤ϵp<10 \leq \epsilon_p < 1), which is the core source of ICI.
\item AFDM adopts Discrete Affine Fourier Transform (DAFT): Inverse DAFT (IDAFT) for modulation and DAFT for demodulation. The chirp parameters are c1c_1 (time-domain chirp) and c2c_2 (frequency-domain chirp), satisfying the optimization condition c1=2(⌊νmax⌋+ξ)2Nc_1 = \frac{2(\lfloor \nu_{\text{max}} \rfloor + \xi)}{2N} (νmax\nu_{\text{max}} is the maximum Doppler shift, ξ\xi is the guard parameter) and c2≤1/(2N)c_2 \leq 1/(2N).
\end{enumerate}

The channel impulse response (Linear Time-Varying, LTV model) is:

g(t,τ)=∑p=0P−1hpej2πfd,ptδ(t−τp)g(t,\tau) = \sum_{p=0}^{P-1} h_p e^{j2\pi f_{d,p} t} \delta(t - \tau_p)

\section{2. ICI Analysis of Uncompensated OFDM}
\subsection{2.1 Transmit and Receive Signal Model}
The OFDM transmit signal after IDFT modulation is:

xOFDM(n)=1N∑k=0N−1Xkej2πkn/N,n=0,1,…,N−1x_{\text{OFDM}}(n) = \frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} X_k e^{j2\pi k n / N}, \quad n=0,1,\dots,N-1

After inserting CP, the time-domain received signal over the double-selective channel is:

rOFDM,uncomp(n)=∑p=0P−1hpej2πfd,pnTsxOFDM(n−τp′)+w(n)r_{\text{OFDM,uncomp}}(n) = \sum_{p=0}^{P-1} h_p e^{j2\pi f_{d,p} n T_s} x_{\text{OFDM}}(n - \tau_p') + w(n)

where τp′=⌊τp/(Ts/N)⌋\tau_p' = \lfloor \tau_p / (T_s / N) \rfloor (delay quantized to sampling points), and w(n)w(n) is Additive White Gaussian Noise (AWGN).

\subsection{2.2 FFT Demodulation and ICI Power}
The received signal is demodulated by FFT. For the kk-th subcarrier, the received signal is:

Yk,uncompOFDM=1N∑n=0N−1rOFDM,uncomp(n)e−j2πkn/NY_{k,\text{uncomp}}^{\text{OFDM}} = \frac{1}{\sqrt{N}} \sum_{n=0}^{N-1} r_{\text{OFDM,uncomp}}(n) e^{-j2\pi k n / N}

Substituting xOFDM(n−τp′)x_{\text{OFDM}}(n - \tau_p') and separating the desired signal, ICI, and noise:

Yk,uncompOFDM=HkXk+ICIOFDM,uncomp,k+WkY_{k,\text{uncomp}}^{\text{OFDM}} = H_k X_k + \text{ICI}_{\text{OFDM,uncomp},k} + W_k

where Hk=∑p=0P−1hpe−j2πkτp′/NH_k = \sum_{p=0}^{P-1} h_p e^{-j2\pi k \tau_p' / N} (channel frequency response), WkW_k is the noise component, and the ICI term is:

ICIOFDM,uncomp,k=1N∑m≠kXm∑p=0P−1hpe−j2πmτp′/N⋅sin⁡(π(νp+k−m))π(νp+k−m)\text{ICI}_{\text{OFDM,uncomp},k} = \frac{1}{N} \sum_{m \neq k} X_m \sum_{p=0}^{P-1} h_p e^{-j2\pi m \tau_p' / N} \cdot \frac{\sin(\pi (\nu_p + k - m))}{\pi (\nu_p + k - m)}

The ICI power is the second moment of the ICI term:

PICI,OFDM,uncomp=E[∣ICIOFDM,uncomp,k∣2]=∑p=0P−1∣hp∣2⋅1N2∑m≠k∣Xm∣2⋅∣sin⁡(π(νp+k−m))π(νp+k−m)∣2P_{\text{ICI,OFDM,uncomp}} = \mathbb{E}[|\text{ICI}_{\text{OFDM,uncomp},k}|^2] = \sum_{p=0}^{P-1} |h_p|^2 \cdot \frac{1}{N^2} \sum_{m \neq k} |X_m|^2 \cdot \left| \frac{\sin(\pi (\nu_p + k - m))}{\pi (\nu_p + k - m)} \right|^2

\section{3. ICI Analysis of Doppler-Compensated OFDM}
\subsection{3.1 Doppler Compensation and Received Signal}
The Doppler compensation operation multiplies the received signal by e−j2π⌊νp⌋n/Ne^{-j2\pi \lfloor \nu_p \rfloor n / N} to eliminate integer multiple Doppler shifts. The compensated signal is:

rOFDM,comp(n)=∑p=0P−1hpej2πϵpn/NxOFDM(n−τp′)+w(n)e−j2π⌊νp⌋n/Nr_{\text{OFDM,comp}}(n) = \sum_{p=0}^{P-1} h_p e^{j2\pi \epsilon_p n / N} x_{\text{OFDM}}(n - \tau_p') + w(n) e^{-j2\pi \lfloor \nu_p \rfloor n / N}

\subsection{3.2 ICI Power Expression}
After FFT demodulation, the ICI term for the kk-th subcarrier is:

ICIOFDM,comp,k=1N∑m≠kXm∑p=0P−1hpe−j2πmτp′/N⋅sin⁡(π(ϵp+k−m))π(ϵp+k−m)\text{ICI}_{\text{OFDM,comp},k} = \frac{1}{N} \sum_{m \neq k} X_m \sum_{p=0}^{P-1} h_p e^{-j2\pi m \tau_p' / N} \cdot \frac{\sin(\pi (\epsilon_p + k - m))}{\pi (\epsilon_p + k - m)}

The corresponding ICI power is:

PICI,OFDM,comp=E[∣ICIOFDM,comp,k∣2]=∑p=0P−1∣hp∣2⋅1N2∑m≠k∣Xm∣2⋅∣sin⁡(π(ϵp+k−m))π(ϵp+k−m)∣2P_{\text{ICI,OFDM,comp}} = \mathbb{E}[|\text{ICI}_{\text{OFDM,comp},k}|^2] = \sum_{p=0}^{P-1} |h_p|^2 \cdot \frac{1}{N^2} \sum_{m \neq k} |X_m|^2 \cdot \left| \frac{\sin(\pi (\epsilon_p + k - m))}{\pi (\epsilon_p + k - m)} \right|^2

\section{4. ICI Analysis of AFDM}
\subsection{4.1 Transmit and Receive Signal Model}
The AFDM transmit signal after IDAFT modulation (core of slant compensation with chirp parameters) is:

xAFDM(n)=1N∑k=0N−1Xkej2π(c1n2+c2k2+kn/N),n=0,1,…,N−1x_{\text{AFDM}}(n) = \frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} X_k e^{j2\pi (c_1 n^2 + c_2 k^2 + k n / N)}, \quad n=0,1,\dots,N-1

After inserting CPP, the time-domain received signal is:

rAFDM(n)=∑p=0P−1hpej2πfd,pnTsxAFDM(n−τp′)+w(n)r_{\text{AFDM}}(n) = \sum_{p=0}^{P-1} h_p e^{j2\pi f_{d,p} n T_s} x_{\text{AFDM}}(n - \tau_p') + w(n)

The DAFT demodulation kernel at the receiver is:

DAFT{r}(k)=1N∑n=0N−1rAFDM(n)e−j2π(c1n2+c2k2+kn/N)\text{DAFT}\{r\}(k) = \frac{1}{\sqrt{N}} \sum_{n=0}^{N-1} r_{\text{AFDM}}(n) e^{-j2\pi (c_1 n^2 + c_2 k^2 + k n / N)}

\subsection{4.2 ICI Suppression Mechanism and Power}
Substituting the received signal into DAFT and simplifying using the optimized c1c_1, let θp=2π(c2k2−fd,pτp′Ts)\theta_p = 2\pi (c_2 k^2 - f_{d,p} \tau_p' T_s), the demodulated signal is:

YkAFDM=∑p=0P−1hpejθpXk−⌊νp⌋⋅1N∑n=0N−1ej2πϵpn/N+ICIAFDM,k+WkY_k^{\text{AFDM}} = \sum_{p=0}^{P-1} h_p e^{j\theta_p} X_{k - \lfloor \nu_p \rfloor} \cdot \frac{1}{N} \sum_{n=0}^{N-1} e^{j2\pi \epsilon_p n / N} + \text{ICI}_{\text{AFDM},k} + W_k

Due to the sparse response of AFDM in the DAFT domain, ICI is confined to a narrow bandwidth. The ICI power is:

PICI,AFDM=E[∣ICIAFDM,k∣2]=∑p=0P−1∣hp∣2⋅1N2∑∣m−k∣>Δ∣Xm∣2⋅∣∑n=0N−1ej2π(c1n2+ϵpn/N)∣2P_{\text{ICI,AFDM}} = \mathbb{E}[|\text{ICI}_{\text{AFDM},k}|^2] = \sum_{p=0}^{P-1} |h_p|^2 \cdot \frac{1}{N^2} \sum_{|m - k| > \Delta} |X_m|^2 \cdot \left| \sum_{n=0}^{N-1} e^{j2\pi (c_1 n^2 + \epsilon_p n / N)} \right|^2

where Δ≤2\Delta \leq 2 is the interference bandwidth, and the exponential sum approaches zero when ∣m−k∣>Δ|m - k| > \Delta.

\section{5. Performance Comparison Conclusion}

  1. Uncompensated OFDM suffers from severe ICI due to full Doppler shifts, with ICI spreading to all subcarriers via sinc side lobes.
  2. Doppler-compensated OFDM reduces ICI by eliminating integer Doppler shifts, but residual fractional Doppler still causes non-negligible ICI.
  3. AFDM achieves time-frequency 2D joint compensation via DAFT and chirp parameters, confining ICI to a narrow bandwidth, leading to significantly lower ICI power than both OFDM schemes.

\end{document}

思考

Evaluating fractional compensation in OFDM

I’m working through how to minimize a residual in OFDM. The key idea involves handling fractional parts of the frequency term by compensating using integer adjustments. By introducing smaller grid hypotheses (non-integer values), I can reduce the residual. Would it require non-FFT demod? After FFT, the coefficient simplifies.

Comparing OFDM and AFDM for Doppler interference reduction

With grid adjustments, residuals can be reduced but not eliminated for all paths. In AFDM, using DAFT chirps transforms Doppler shifts into delay shifts. This results in limited interference for off-grid Doppler due to chirp pre- and post-processing, achieving near-zero interference for integer Doppler with fractional still reducing it effectively.

Comparing AFDM and filter bank for ICI suppression

For multi-path scenarios with varying Doppler shifts, AFDM outperforms, as it handles both delay and Doppler at once, leading to sparse responses and minimizing global ICI. In contrast, the filter bank method requires O(M) operations. AFDM uses a DAFT per symbol and incorporates equalization in the delay-Doppler domain. The filter bank acts as a matched filter, preserving energy using Parseval’s theorem.

Exploring ICI suppression and AFDM metrics

When considering ICI suppression, rotating the signal changes its magnitude but doesn't conserve energy due to multiple de-rotations across hypotheses. I think a metric like the off-diagonal energy of the effective channel matrix (E||H - diag(H)||^2) is needed. I’ll craft a more structured, analytical solution involving both OFDM and AFDM techniques.

Comparing ICI suppression in filter bank and AFDM

For suppressing ICI, the filter bank chooses the best q for each path, but multiple Dopplers require complex multi-branch combining. ICI suppression is best with a dominant Doppler (q chosen near ν_dominant). In AFDM, the channel is sparse after DAFT demod, and off-grid Doppler causes limited spreading, resulting in lower ICI compared to OFDM’s broader interference.

Thought for 1m 27s

回答

可以，把“PSR”先放一边，用解析解来对比更干净。下面我用一个统一的离散基带模型，把：

  • (A) OFDM + 每个 Doppler 都补偿的滤波器组（Doppler-hypothesis bank）
  • (B) AFDM（DAFT/IDAFT + chirp 参数）

放在同一框架下比较 ICI（子载波间干扰）的“结构”和“能量规模”，并说明在什么条件下谁更强。

1) 统一的离散 LTV 信道模型（去 CP/CPP 后）

取一帧长度 NN（快时间采样），慢时间索引先不引入（单个 OFDM 符号内看 ICI 就够了）：

r[n]=∑p=0P−1hp ej2πνpnN x[(n−τp)N]+w[n],n=0,…,N−1r[n]=\sum_{p=0}^{P-1} h_p\,e^{j2\pi \nu_p \frac{n}{N}}\,x[(n-\tau_p)_N]+w[n],\quad n=0,\dots,N-1
  • τp\tau_p 为整数采样延迟（CP/CPP 足够长保证循环卷积）
  • νp\nu_p 为归一化 Doppler（以子载波间隔为单位），可为非整数
  • 关键：ICI 的根源就是这个 ej2πνpn/Ne^{j2\pi \nu_p n/N} 破坏了 FFT 的正交性

2) OFDM 的 ICI：本质是 Dirichlet 核（而不是简单 sinc）

OFDM 发射：

x[n]=1N∑m=0N−1X[m]ej2πmnNx[n]=\frac{1}{\sqrt{N}}\sum_{m=0}^{N-1}X[m]e^{j2\pi \frac{mn}{N}}

接收 FFT：

Y[k]=1N∑n=0N−1r[n]e−j2πknNY[k]=\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}r[n]e^{-j2\pi \frac{kn}{N}}

代入后得到一个标准形式（把“期望项”和“ICI项”统一写成矩阵）：

Y[k]=∑m=0N−1Hk,mOFDM X[m]+W[k]Y[k]=\sum_{m=0}^{N-1} H_{k,m}^{\text{OFDM}}\,X[m] + W[k]

其中耦合系数为

Hk,mOFDM=∑php e−j2πmτpN  CN(νp+k−m)H_{k,m}^{\text{OFDM}}=\sum_{p} h_p\,e^{-j2\pi \frac{m\tau_p}{N}}\;C_N(\nu_p+k-m)

而

CN(a)≜1N∑n=0N−1ej2πanN=ejπaN−1N  sin⁡(πa)Nsin⁡(πa/N)C_N(a)\triangleq \frac{1}{N}\sum_{n=0}^{N-1}e^{j2\pi a \frac{n}{N}} = e^{j\pi a\frac{N-1}{N}}\;\frac{\sin(\pi a)}{N\sin(\pi a/N)}

结论 1（非常关键）：

  • 当 aa 非整数时，∣CN(a)∣|C_N(a)| 的旁瓣是Dirichlet 核，会“长尾”地铺到很多子载波（不是只影响邻近几个）。
  • 你的参考文档里用的 sin⁡(π(⋅))π(⋅)\frac{\sin(\pi(\cdot))}{\pi(\cdot)} 是大 NN 的近似（sinc），但严格离散情形应是上面的 Dirichlet 形式；这会影响对“干扰扩散范围/能量规模”的判断。

于是 ICI 功率（以每个子载波平均、符号独立同分布）可用“非对角能量”刻画：

PICIOFDM∝∑m≠kE∣Hk,mOFDM∣2≈∑p∣hp∣2∑l≠0∣CN(νp+l)∣2P_{\text{ICI}}^{\text{OFDM}} \propto \sum_{m\neq k}\mathbb{E}|H_{k,m}^{\text{OFDM}}|^2 \approx \sum_p |h_p|^2 \sum_{l\neq 0}|C_N(\nu_p+l)|^2

这里 l=k−ml=k-m。

3) “每个 Doppler 都补偿”的滤波器组：等价于把 νp\nu_p 平移成 νp−q\nu_p-q

你现在做的方案（每个 Doppler bin 都做一次快时间去旋转）本质是：

对某个假设 qq（可以是整数 bin，也可以做细网格）：

rq[n]=r[n]  e−j2πqnNr_q[n]=r[n]\;e^{-j2\pi q \frac{n}{N}}

再做 FFT。

那么耦合系数直接变成：

Hk,m(q)=∑php e−j2πmτpN  CN(νp−q+k−m)H_{k,m}^{(q)}=\sum_p h_p\,e^{-j2\pi \frac{m\tau_p}{N}}\;C_N(\nu_p-q+k-m)

也就是把每条径的 Doppler 从 νp\nu_p 变成了残差

δp(q)≜νp−q\delta_p(q)\triangleq \nu_p-q

因此，对“某个固定的 qq”而言，ICI 能量就是：

PICI(q)≈∑p∣hp∣2∑l≠0∣CN(δp(q)+l)∣2P_{\text{ICI}}^{(q)} \approx \sum_p |h_p|^2 \sum_{l\neq 0}|C_N(\delta_p(q)+l)|^2

关键性质

  • 如果所有路径 Doppler 都相同且为整数：选 q=νq=\nu，则 δp=0\delta_p=0，ICI 理论上可为 0（只剩多径相位和延迟造成的频域响应，仍是对角的）。
  • 如果是单径或“一个主径 Doppler 占绝对主导”：选 q≈νdomq \approx \nu_{\text{dom}}，主径残差很小，主径 ICI 会显著下降。
  • 但如果多径具有不同 νp\nu_p：你只能选一个 qq，必然存在很多路径残差 δp(q)\delta_p(q) 不小 → 这些路径仍产生“长尾 Dirichlet”型 ICI，无法同时压住。

所以滤波器组的“上限”取决于：信道是不是“近似单 Doppler”。

你代码里“把每个 qq 都算一遍并组成整张 RDM”的含义

这其实是在算一个“多假设匹配滤波输出集合”。
它不是一个能量守恒的正交变换（每一列都用不同预处理），所以“全图能量/旁瓣指标”会显得怪——这也是你之前用 PSR 会很纠结的根源之一。

真正公平的使用方式是：用估计/检测选出最可能的 qq，然后只拿这个 qq 的输出做后续解调/成像。

4) AFDM（DAFT）的核心：把“Doppler 引起的子载波耦合”变成“近似稀疏的移位”

AFDM 的调制/解调可抽象成：

x=AHs,s^=Arx = \mathbf{A}^H s,\quad \hat{s}=\mathbf{A} r

其中 A\mathbf{A} 是 DAFT 矩阵（含 chirp 参数 c1,c2c_1,c_2 的二次相位核）。

在满足典型 AFDM 设计条件（CPP 足够长、c1c_1 选取覆盖最大 Doppler、并满足论文里的“可分离/可稀疏化”条件）时，AFDM 在 DAFT 域把每条径近似映射成二维移位（delay-Doppler 的稀疏 tap）：

  • 对“整数 Doppler/整数 delay”的理想情况：DAFT 域响应接近单点（或极窄带），ICI/IDE（干扰扩散）天然很小。
  • 对“非整数 Doppler”的情况：扩散不再是 Dirichlet 核那种长尾，而是由二次相位和（高斯和/二次指数和）决定，其典型性质是能量更集中、旁瓣衰减更快/更可控，并且扩散宽度可以通过 guard/参数设计约束在较窄邻域（这是 AFDM 的卖点之一）。

你参考里写的：

∑n=0N−1ej2π(c1n2+ϵn/N)\sum_{n=0}^{N-1} e^{j2\pi (c_1 n^2 + \epsilon n/N)}

就是这种二次指数和的影子。它与 OFDM 的线性相位和（Dirichlet 核）相比，一个关键差别是：

  • OFDM 的干扰核是 “线性相位和” → Dirichlet（长尾）
  • AFDM 的干扰核是 “二次相位和” → 更强的相消，干扰更局部/可控

（严格的界可以用高斯和性质给出：在不少参数取值下，离开主支持区的幅度上界可做到比 Dirichlet 更“瘦”。）

5) 谁对 ICI 抑制更好？给一个清晰结论（按信道类型分）

情况 A：单径 / 所有路径共享同一 Doppler（“单 Doppler 信道”）

  • 滤波器组选对 q=νq=\nu（且 ν\nu 为整数）：OFDM 的 ICI 可以做到几乎为 0
  • AFDM 也能很好，但这时它的优势不明显（甚至你还要付 chirp/DAFT 的实现复杂度）

✅ 结论：单 Doppler 场景，滤波器组（选对 qq）不输 AFDM，甚至可“零 ICI”。

情况 B：多径、且不同路径 Doppler 不同（“多 Doppler 双选信道”）

  • 滤波器组：只能用一个 qq 去旋转，最多照顾一个主径；其它路径残差 Doppler 仍产生 Dirichlet 长尾 ICI
  • AFDM：目标就是把“多径+多 Doppler”在 DAFT 域做成稀疏/窄带扩散，更适合“多 Doppler 叠加”

✅ 结论：多 Doppler 多径时，AFDM 通常比“单 qq”滤波器组更能整体压住 ICI。
如果你允许滤波器组做“多 qq 分支联合估计 + 多径分离/联合检测”，那能逼近 ML，但复杂度会飙升，且本质上你在做“二维稀疏估计/多假设检测”，就和 AFDM 的优势区间重叠了。

情况 C：Doppler 非整数（离栅）

  • 滤波器组如果只允许整数 qq：残差 δ∈[−0.5,0.5)\delta\in[-0.5,0.5) 仍带来 ICI
  • 如果允许更细 Doppler 网格（分数 qq）：可以进一步减小残差，但你需要更多分支（复杂度线性增长）
  • AFDM：通过参数设计让离栅 Doppler 的扩散更局部，通常更“经济”

✅ 结论：离栅 Doppler 下，AFDM 的“局部扩散 + 稀疏结构”往往更占优；滤波器组要靠更细网格和更复杂的选择/联合检测来追。

6) 建议你们用什么“解析对比指标”（替代 PSR）

如果要纯理论比较 ICI 抑制，我建议直接比有效信道矩阵的非对角能量比（ODR / ICI ratio）：

ρ≜∑k∑m≠kE∣Hk,m∣2∑kE∣Hk,k∣2\rho \triangleq \frac{\sum_{k}\sum_{m\neq k}\mathbb{E}|H_{k,m}|^2}{\sum_{k}\mathbb{E}|H_{k,k}|^2}
  • OFDM：Hk,mH_{k,m} 由 Dirichlet 核决定，长尾 → ρ\rho 往往不小（尤其多 Doppler）
  • 滤波器组：ρ(q)\rho(q) 取决于残差 δp(q)\delta_p(q)。你可以比较
    • 最佳分支：min⁡qρ(q)\min_q \rho(q)（代表“选对 qq”的能力上限）
    • 或主径分支：q=round(νdom)q=\text{round}(\nu_{\text{dom}})
  • AFDM：HH 在 DAFT 域更稀疏/更窄带扩散，ρ\rho 通常更小，且对多径多 Doppler 更稳定

如果你愿意，我也可以把上面这个 ρ\rho（或“主邻域外能量占比”）直接接到你们的仿真里当评价指标，比 PSR 更贴近“ICI 抑制”。

一句话总结

  • 滤波器组（逐 Doppler 补偿）：本质是把 OFDM 的 ICI 核从 ν\nu 平移到残差 ν−q\nu-q。对“单 Doppler/主径强”的信道非常有效，选对 qq 可接近零 ICI；但对“多径多 Doppler”无法同时对齐所有路径，ICI 仍呈 Dirichlet 长尾扩散。
  • AFDM：通过 chirp/DAFT 把双选信道在变换域做成稀疏/局部扩散，整体上更适合多 Doppler 多径，通常能得到更低、更稳定的 ICI（尤其离栅 Doppler 时优势更明显）。

如果你把你们要对比的“AFDM 参数设计”（c1,c2,ξc_1,c_2,\xi 的取法、是否 CPP、最大 Doppler 假设）给我，我可以把上面这套解析推导进一步“落到你们的符号体系”，写成一段更像论文的推导（包括从 CN(⋅)C_N(\cdot) 到 ICI 功率闭式/上界，以及 AFDM 二次指数和的界）。

問題

2512.08352v1.pdf

c_1,c_2先不用给具体的值 我们都先带符号进去，而不是带着具体的数。让我们从2D-FFT的角度，一步一步来分析没有补偿的OFDM，补偿的OFDM和AFDM的有效信道矩阵的模糊函数，具体来说，这个是2D-FFT之后打模平方的操作。我们先来参考下面文章的建模思路来推导这个解析解

問題

2512.08352v1.pdf

c_1,c_2先不用给具体的值 我们都先带符号进去，而不是带着具体的数。让我们从2D-FFT的角度，一步一步来分析没有补偿的OFDM，补偿的OFDM和AFDM的有效信道矩阵的模糊函数，具体来说，这个是2D-FFT之后打模平方的操作。我们先来参考下面文章的建模思路来推导这个解析解，尽快给我答案，不要停止思考。这篇文章其实是忽略了ICI的影响的，我们先基于这个框架，建立一个考虑ICI影响的框架解析解，即模糊函数的解析解

分享這個問答