基于相干衍射成像模拟的matlab仿真（完整源码）.rar

function output = Propagate(input,propagator,dx,lamda,z) % 波前传播 %%%%%%%%%%%%%%%%%%%%%%%%%%% % input: the wavefront to propagate % propagator: one of 'fourier','fresnel' or 'angular spectrum' % dx: the pixel spacing of input wavefront % wavelength: the wavelength of the illumination % z: the distance to propagate % output: the distance to propagate % 频域坐标 k = 2*pi/lamda; [ysize, xsize] = size(input); x = -xsize/2:xsize/2-1; y = -ysize/2:ysize/2-1; fx = x./(dx*xsize); fy = y./(dx*ysize); [fx, fy] = meshgrid(fx,fy); switch propagator case 'fourier' if z>0 output = fftshift(fft2(fftshift(input))); % output = fft2(input); else output = ifftshift(ifft2(ifftshift(input))); % output = ifft2(input); end case 'angular spectrum' % Calculate phase distribution for each plane wave component w = sqrt(1/lamda^2 - fx.^2 - fy.^2); % exclude evanescent waves notEvanescent = imag(w) == 0; % Compute FFT of input F = fftshift(fft2(fftshift(input))); % multiply FFT by phase-shift and inverse transform output = ifftshift(ifft2(ifftshift(F.*exp(2i*pi*z*w).*notEvanescent))); case 'fresnel' % Calculate approx phase distribution for each plane wave component w = fx.^2 + fy.^2; % Comput FFT F = fftshift(fft2(fftshift(input))); % multiply by phase-shift and inverse transform output = ifftshift(ifft2(ifftshift(F.*exp(-1i*pi*z*lamda*w)))); end end