Clear Filters
Clear Filters

Why do I get a wrong plot?

7 views (last 30 days)
Daigo
Daigo on 2 Jan 2022
Commented: Minh on 21 Mar 2023
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 Comments
Daigo
Daigo on 3 Jan 2022
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
Star Strider on 3 Jan 2022
‘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.
.

Sign in to comment.

Accepted Answer

Daigo
Daigo on 6 Jan 2022
Edited: Daigo on 6 Jan 2022
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!

More Answers (0)

Products


Release

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!