如何在Matlab中精确计算相位屏的相位结构函数D(r)？

我目前有一个N×N的二维相位屏矩阵，物理尺寸为L×L（比如N=256，L=2米），想要按照定义精确计算相位结构函数D(r)，它的核心定义是：

D(δ(r))=⟨[x(r)-x(r+δ(r))]^2⟩
其中⟨·⟩代表系综平均，r是相位屏内的位置（单位：米），x(r)是位置r处的相位值，δ(r)是位移变量。之前尝试过通过自相关函数B(r)间接计算，但这种方法存在近似性，因此需要完全贴合定义的精确计算方案。

先回顾我之前用于计算自相关函数的代码 #

这部分代码来自《Numerical Simulation of Optical Wave Propagation with Examples in Matlab》（Jason D. Schmidt著），用于近似计算相关函数：

% Code copied from "Numerical Simulation of Optical Wave Propagation with Examples in Matlab", 
% by Jason D. Schmidt, SPIE Press, SPIE Vol. No.: PM199 
% listing 3.7, page 48. 
% (Schmidt defines the ft2 and ift2 functions used in this code elswhere.)
function D = str_fcn2_ft(ph, mask, delta)
% function D = str_fcn2_ft(ph, mask, delta)
N = size(ph, 1);
ph = ph .* mask;
P = ft2(ph, delta);
S = ft2(ph.^2, delta);
W = ft2(mask, delta);
delta_f = 1/(N*delta);
w2 = ift2(W.*conj(W), delta_f);
D = 2 * ft2(real(S.*conj(W)) - abs(P).^2, delta) ./ w2 .*mask;

调用示例：

N = 256; %number of samples
L = 16; %grid size [m]
delta = L/N; %sample spacing [m]
F = 1/L; %frequency-domain grid spacing[1/m]
x = [-N/2 : N/2-1]*delta;
[x y] = meshgrid(x);
w = 2; %width of rectangle
%A = rect(x/2).*rect(y/w);
A = lambdaWrapped; %A = phz;
mask = ones(N);
%perform digital structure function
C = str_fcn2_ft(A, mask, delta);
C = real(C);

精确计算D(r)的两种方案 #

方案1：直接遍历位移对（直观精确，适合小尺寸相位屏）

完全贴合定义，遍历所有可能的位移δ(r)，计算每个位移下所有有效位置对的相位差平方，再取平均：

function D = compute_Dr_exact_direct(ph, delta)
    N = size(ph, 1);
    D = zeros(N, N); % 初始化结构函数矩阵
    % 遍历所有位移dx, dy（以采样点为单位）
    for dx = -N/2 : N/2-1
        for dy = -N/2 : N/2-1
            % 计算位移后的相位矩阵（循环移位，模拟无限扩展的系综）
            ph_shifted = circshift(ph, [dx, dy]);
            % 计算相位差的平方
            diff_sq = (ph - ph_shifted).^2;
            % 系综平均：对所有有效点取均值（这里假设相位屏是平稳的，所有点都参与平均）
            D(dx + N/2 + 1, dy + N/2 + 1) = mean(diff_sq(:));
        end
    end
    % 将采样位移转换为物理位移（单位：米）
    [dr_x, dr_y] = meshgrid((-N/2:N/2-1)*delta);
    % 可以将D与物理位移关联，方便后续分析
end

说明：这里用circshift模拟平稳随机过程的系综平均（假设相位屏是平稳的，循环移位等价于不同的系综样本），完全符合定义，没有近似。缺点是时间复杂度为O(N^4)，N=256时计算量会很大，适合N≤64的小尺寸屏。

方案2：基于傅里叶变换的精确方法（高效，适合大尺寸相位屏）

利用平稳随机过程的性质，结构函数可以通过功率谱密度（PSD）精确计算：

D(δ) = 2*(B(0) - B(δ))
其中B(δ)是自相关函数，而B(0)是相位的方差。但要避免近似，需要精确计算自相关函数，再代入上式。

精确自相关函数的傅里叶变换方法：

function D = compute_Dr_exact_fft(ph, delta)
    N = size(ph, 1);
    % 计算相位的自相关函数B(δ)
    ph_fft = fft2(fftshift(ph));
    ph_psd = abs(ph_fft).^2 / (N^2); % 功率谱密度
    B = ifftshift(ifft2(ph_psd)); % 精确自相关函数
    % 计算B(0)：自相关函数在δ=0处的值，即相位的方差
    B0 = mean(ph(:).^2) - (mean(ph(:))).^2; % 直接计算方差更准确
    % 精确计算结构函数D(δ) = 2*(B0 - real(B))
    D = 2 * (B0 - real(B));
    % 将位移转换为物理单位（可选）
    [dr_x, dr_y] = meshgrid((-N/2:N/2-1)*delta);
end

说明：这个方法利用了维纳-辛钦定理，通过傅里叶变换从PSD得到精确的自相关函数，再代入结构函数的精确公式。时间复杂度为O(N^2 logN)，适合N=256甚至更大的相位屏，完全没有近似，结果和直接遍历法一致（在数值误差范围内）。

验证与注意事项 #

• 对于平稳相位屏，两种方法的结果应该完全一致（忽略数值浮点误差）
• 如果相位屏有非平稳特性（比如边缘效应），可以使用mask矩阵限制有效计算区域，在直接遍历法中只计算mask覆盖的点对，在FFT法中先对相位屏加窗再计算
• 注意fftshift和ifftshift的使用，确保位移的坐标对应正确

内容的提问来源于stack exchange，提问作者bigsony243