matlab - 根据定义绘制相屏的相结构函数
问题描述
我已经有一个相位屏幕(一个 2-D NxN 矩阵和 LxL 的大小比例,例如:N = 256,L = 2 米)。
我想找到相位结构函数 - 由 D(delta(r)) = <[x(r)-x(r+delta(r))]^2> (<.> 定义的 D(r) 是整体平均, r 是相位屏中的位置(以米为单位),x 是相位屏中某点的相位值,delta(r) 在 Matlab 程序中是可变的且不固定)。你对我的目的有什么建议吗?
P/S:我试图通过自相关(定义为 B(r))计算 D(r),但这个计算仍然保留了一些近似值。因此,我想精确计算 D(r) 的结果。请您查看这张图片以更好地理解 D(r) 和 B(r) 的定义。下面是我计算 B(r) 的函数代码。
% 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;`
%fire run
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) 的一种方法是通过随机采样:您在屏幕上选择两个随机点,确定它们的距离和相位差的平方,然后更新累加器:
phi = rand(256,256)*(2*pi); % the data, phase
N = size(phi,1); % number of samples
L = 16; % grid size [m]
delta = L/N; % sample spacing [m]
D = zeros(1,sqrt(2)*N); % output function
count = D; % for computing mean
for n = 1:1e6 % find a good amount of points here, the more points the better the estimate
coords = randi(N,2,2);
r = round(norm(coords(1,:) - coords(2,:)));
if r<1
continue % skip if the two coordinates are the same
end
d = phi(coords(1,1),coords(1,2)) - phi(coords(2,1),coords(2,2));
d = mod(abs(d),pi); % you might not need this, depending on how A is constructed
D(r) = D(r) + d.^2;
count(r) = count(r) + 1;
end
I = count > 0;
D(I) = D(I) ./ count(I); % do not divide by 0, some bins might not have any samples
I = count < 100;
D(I) = 0; % ignore poor estimates
r = (1:length(D)) * delta;
plot(r,D)
如果您需要更高的精度,请考虑插值。将随机坐标计算为浮点值,并对相位进行插值以获得样本之间的值。D
然后需要更长的时间,索引为round(r*10)
或类似的东西。您将需要更多的随机样本来填充更大的累加器。
推荐阅读
- javascript - Pouchdb 使用 ID 最小化文档修订如何仅查询主文档?
- excel - 使用公式从不同的工作表中获取过滤后的唯一值列表
- jquery - 通过 AJAX (Laravel) 将数据从控制器传递到模态
- mysql - SQL 包含或喜欢
- excel - Excel 表查找值
- .net-core - 将 NuGet 包还原到缓存
- inheritance - Kotlin 列表是……“可实例化”的接口吗?
- python - 按 Pandas 中提到的 TimeFrame 对 Dataframe 进行分组
- react-native - React native expo - 以编程方式获取android中的默认字体系列
- reactjs - 打字稿不允许将变量作为值放入样式数组