Why do I get a wrong plot?

Question

Daigo 2022-1-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1621030-why-do-i-get-a-wrong-plot

评论： Minh 2023-3-21

I'd like to make a plot of a quantity called "PSF" but I cannot get a correct result for some reason. I will first explain the theory and then my code.

Background

Point spread function (PSF) is computed by the inverse-Fourier transform of the overall OTF:

The optical transfer function (OTF) for an imaging system in optical turbulence can be modeled as the product of the atmospheric OTF and the diffraction OTF:

Note that the OTFs are circularly symmetric in the frequency domain, and they are functions of

, where

and

are the spatial frequencies.

The atmospheric OTF:

where λ is the wavelength, ℓ is the camera focal length,

is the Fried parameter, and D is the aperture diameter. The parameter α works as a switch. When

, it models a long exposure, and when

, it models a short exposure.

The diffraction OTF:

where

is the cut-off frequency given by

.

Simulation

I used the following MATLAB code to compute the PSFs based on the theories above. I simply want to visualize the cross-sections of the PSFs for

and

.

clear all;

% Unit conversion parameters

nm2m = 1e-9;

m2um = 1e+6;

% Parameters

D = 0.2034; % [m], diameter of camera aperture

fl = 1.2; % [m], focal length

wvl = 525*nm2m; % [m], central wavelength

fn = fl/D; % f-number

rho_c = 1/(wvl*fn); % [Hz] optical cut-off frequency

delta_f = 1/(2*rho_c); % [m], Nyquist pixel spacing

r0 = 0.0478; % [m], Coherence length

% Sampling grid

N = 1024; % Number of sampling points

delta = delta_f/2; % Grid spacing (spatial domain)

n_vec = -N/2 : 1 : N/2-1;

df = 1 / (N*delta); % Grid spacing (frequency domain)

[fX,fY] = meshgrid(n_vec*df); % Sampling grid (frequency domain)

[~,rho] = cart2pol(fX, fY);

mask = circ(fX, fY, 2*rho_c);

% Diffraction-limited OTF

rho_n = rho/(2*rho_c);

Hdif = (2/pi)*(acos(rho_n)-rho_n.*sqrt(1-rho_n.^2));

figure;

cc = 0;

for alpha = [0 1]

cc = cc + 1;

% Atmospheric OTF

% (alpha = 0 -> long exposure, alpha = 1 -> short exposure)

numer = wvl*fl*rho;

term1 = -3.44*(numer/r0).^(5/3);

term2 = 1-alpha*(numer/D).^(1/3);

term_all = term1.*term2;

Hatm = exp(term_all);

Hatm = Hatm.*mask;

% Overall OTF

Hall = Hatm .* Hdif;

% PSF (theory)

PSF = ift2(Hall, df);

PSF = abs(PSF);

PSF = PSF/max(max(PSF));

PSF_slice = PSF(N/2+1, :);

% Plot

subplot(2,2,cc);

x_seq = delta*n_vec;

plot(x_seq/delta_f, PSF_slice), grid on; xlim([-30,30]);

xlabel('x (Nyquist pixels)'); ylabel('h(x,0)');

end

function g = ift2(G, delta_f)

% Function to compute 2D discrete inverse Fourier transform

N = size(G,1);

g = ifftshift(ifft2(ifftshift(G))) * (N * delta_f)^2;

end

function z = circ(x, y, D)

% Circular function

r = sqrt(x.^2 + y.^2);

z = double(r<D/2);

z(r==D/2) = 0.5;

end

Below is the comparison of the result from a paper and my results of this code:

Even though I see a perfect agreement of the results for

, I don't get a correct result for

. That means, the computation of

is probably correct but there is something wrong with

only when

. I am wondering if there is an undersampling problem or numerical error in my computation but I cannot identify the cause. Do you have any idea what I'm doing wrong?

3 个评论
显示 1更早的评论隐藏 1更早的评论

Daigo 2022-1-3

I also suspected that it looks like an envelope function. I double-checked the paper but it doesn't say so... The connected red dots are the ones obtained by something similar to a Monte Carlo simulation, and they are validated against the theoretical results (blue graphs). So what we want to reproduce is the blue graph.

Star Strider 2022-1-3

‘I double-checked the paper but it doesn't say so ...’

... which does not necessarily mean that it isn’t.

I’m not certain. I’m having difficulty following the code, so the only suggestion I have is to perhaps give ‘n_vec’ a finer resolution (smaller step increment), or zero-pad it to give finer resolution.

.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Daigo 2022-1-6

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1621030-why-do-i-get-a-wrong-plot#answer_869560

编辑：Daigo 2022-1-6

I contacted the author of the paper. It turned out that the equation for the diffraction OTF has a typo. The correct one is

After I modified my code to reflect this change, I got the correct result. After all, the problem was not related to MATLAB at all but thank you for your suggestions!