Why do I get a wrong plot?

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?

[Star Strider]

It could be that the ‘correct results’ for are the envelope of that function, in which case there is acceptable agreement. There is not enough detail in the ‘correct results’ plot to determine what the connected dots represent.

[Daigo]

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]

‘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]

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!

[Minh]

What was the paper called?