![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/851050/image.png)
Why do I get a wrong plot?
5 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/855310/image.png)
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/855315/image.png)
Note that the OTFs are circularly symmetric in the frequency domain, and they are functions of
, where
and
are the spatial frequencies.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/855320/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/855325/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/855330/image.png)
- The atmospheric OTF:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/855335/image.png)
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/855340/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/855345/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/855350/image.png)
- The diffraction OTF:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/855355/image.png)
where
is the cut-off frequency given by
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/855360/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/855365/image.png)
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
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/855345/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/855350/image.png)
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/851020/image.jpeg)
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?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/855345/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/855350/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/855370/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/855375/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/855350/image.png)
3 comentarios
Star Strider
el 3 de En. de 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.
.
Más respuestas (0)
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!