【论文精读系列】之《Turbulence-Induced 2D Correlated Image Distortion》代码分析一

AAvariance.mAAcorrelation.mcalc_correlation_Cv.m

项目网址：Turbulence-Induced Image Distortion 论文代码：Matlab files for simulating the distortion vector field 注：论文作者代码使用MATLAB编写，在这里只进行展示注解。

AAvariance.m

function [ Variance, VarianceCoef ] = AAvariance( Cn_2, D, L, lamda, wave, CalcType ) %calculating the variance of the AOA. %If CalacType = 'full', then we ucalculate the varicance based on a formula %from the paper "Closed form approximations for the AOA variance of plane %and spherical waves prpagating through homogenous and isotoropic turbulence. %If not, the varianve is calculated according to the known equations, which %are good only if q>>1 beta = 0.5216; f = sqrt(lamda*L); q = D/f; if strcmp(CalcType, 'full') VarianceCoef = sqrt(3)/16*gamma(1/6)*gamma(8/3)*((beta/2)^(-1/3)); if strcmp(wave, 'plane') VarianceCoef =VarianceCoef*(1+6/5*((pi/2)^(1/6))*((beta*q)^(1/3))*((1+((pi/2)^2)*((beta*q)^4))^(5/12)... *sin(5/6*atan(2/(pi*(beta*q)^2))))); else %%sphericak wave F = hypergeom([1/6 17/6], 23/6, 1-0.5i*pi*(beta*q)^2); VarianceCoef = 3/8*VarianceCoef*(1+16/17*real((-2i/pi/((beta*q)^2))^(-1/6)*F)); end else %% requires q>>1 if strcmp(wave, 'sphere') VarianceCoef = 1.09; else %% plane wave VarianceCoef = 2.91; end end Variance = VarianceCoef*Cn_2*L*(D^(-1/3)); end

注： ① 代码实现计算光波前的 $AA$ 的方差，涉及论文公式 1,2,3,4,5,6 ② 函数：strcmp/gamma/atan/hypergeom/real ③ 代码分了 $q$ 是否远大于1，平面波还是球面波，共4种情况讨论

AAcorrelation.m

function b = AAcorrelation(gamma,wave,beta_max,nbeta,nkappa,nxi) %Based on: Eq. (4) in Belen'ki et.al. "Turbulence-induced edge image %waviness: theory and experiment." % %AA (AOA) (Angle Of Arrival) correlation scale from the AA correlation function %b(theta) = < phi(theta0)*phi(theta0+theta)> / sigmaA_2 %phi :the AA %sigmaA_2 : the AA variance %theta : angular separation between two sources %< > : statistical average % %Output: %Structure b with fields: %beta = theta/theta_D : vector of normalized angles %b_long : the correlation of the LONGITUDINAL AA, PARALLEL to angular displacement. %b_lat : the correlation of the LATERAL AA, TRANSVERSE to angular displacement. % %Input: %gamma = L0/D; %L0 outer scale %wave = 'plane' or 'sphere' [default = 'sphere'] %beta_max = largest value of beta, in units of gamma [default = 2000, so largest beta is 2000*gamma] %nbeta = number of normalized angle values [default = 1000] %nkappa = number of kappa values between 1e-2 and 1e2 [default = 2000] %nxi = number of xi values between 0 and 1 [default = 1000] % %Note: %Cn_2 : turbulence refractive index structure function % - horizontal path : we assume that Cn_2 is constant, not needed % Check inputs if(~exist('wave', 'var')), wave = 'sphere'; end if(~exist('beta_max', 'var')) switch wave case 'plane', beta_max = 16*gamma; case 'sphere', beta_max = 1000*gamma; end end if(~exist('nbeta', 'var')), nbeta = 1000; end if(~exist('nkappa', 'var')), nkappa = 2000; end if(~exist('nxi', 'var')), nxi = 1000; end %% %kappa %vector of the spectrum variable - integrates to inf %kappa = 2*pi/lambda %lambda = [8:0.01:13]*1e-6; %wavelength - spectral response of imager in meters % kappa = [0:0.05:50]; nkappa = 1000; kappa = logspace(-2,2,nkappa); nxi = 1000; xi = linspace(0,1,nxi); % Normalized by L beta = logspace(-2,log10(beta_max),nbeta); T = round(15e-9 * nbeta * nkappa * nxi); % Expected computing time switch wave case 'plane' Gk = ((2*gamma*kappa).^2 + 1).^(-11/6) .* (besselj(2, kappa).^2./kappa); I_minus = zeros(size(beta)); I_plus = zeros(size(beta)); for(t = 1:nbeta) J0 = besselj(0, 2*beta(t)*kappa); J2 = besselj(2, 2*beta(t)*kappa); I_minus(t) = trapz(kappa, Gk .* (J0 - J2)); I_plus(t) = trapz(kappa, Gk .* (J0 + J2)); end case 'sphere' Q = (1-xi).^(5/3); Gk = ((2*gamma*kappa).^2 + 1).^(-11/6) .* (besselj(2, kappa).^2./kappa); I_minus = zeros(size(beta)); I_plus = zeros(size(beta)); h = waitbar(0,['Computing AA correlation (~', num2str(T), ' min)']); for(t = 1:nbeta) J0 = besselj(0, 2*beta(t)*(xi .* (1-xi))'*kappa); J2 = besselj(2, 2*beta(t)*(xi .* (1-xi))'*kappa); inner_minus = trapz(kappa, repmat(Gk, [nxi 1]) .* (J0 - J2), 2); inner_plus = trapz(kappa, repmat(Gk, [nxi 1]) .* (J0 + J2), 2); I_minus(t) = trapz(xi, inner_minus); I_plus(t) = trapz(xi, inner_plus); waitbar(t / nbeta, h) end close(h) % Fix inaccuracies i = find(diff(beta > 1000)); I_minus(i:end) = 0; end b_long = I_minus / I_minus(1); b_lat = I_plus / I_plus(1); b = struct('beta', beta, 'long', b_long, 'lat', b_lat); return %% test gamma = 4; tic b = AAcorrelation(gamma, 'sphere', 1000); toc plot(b.beta, b.long, b.beta, b.lat), axis([0 b.beta(end) -0.2 1]) xlabel('beta = theta/theta_{D}'), ylabel('b'), legend('b_{long}','b_{lat}') semilogx(b.beta, b.long, b.beta, b.lat) xlabel('theta/theta\_D'), ylabel('b'), legend('long', 'lat') % Figure % b1 = AAcorrelation(1); b4 = AAcorrelation(4); b10 = AAcorrelation(10); % save('AA_gamma1', 'b1') % save('AA_gamma4', 'b4') load('AA_gamma1'); load('AA_gamma4');b4=b; load('AA_gamma10'); b10=b; %% figure(1), hold on; plot(b10.beta, b10.lat, 'r-', b10.beta, b10.long, 'r--','LineWidth',1.75); plot(b4.beta, b4.lat, 'b-', b4.beta, b4.long, 'b--','LineWidth',1.75); % plot(b1.beta, b1.lat, '-.', b1.beta, b1.long, '-.'); plot([0 100], [0 0], 'k-') legend('b_{\perp} ({\rho} = 10)','b_{||} ({\rho} = 10)', ... 'b_{\perp} ({\rho} = 4)','b_{||} ({\rho} = 4)'); % 'b_{lat} ({\it\gamma} = 1)','b_{long} ({\it\gamma} = 1)'); set(gca, 'fontsize', 16, 'LineWidth',1.2); box on; xlabel('{\it\eta} = {\it\theta} /{\it\theta_{D}}'), ylabel('b'); axis([0 400 -0.1 1.001]), hold off; % Use the same color for the same rho. Say blue for rho=10, red for rho=4. % % * Use solid line for the perp, and dashed line for the parallel, in each of the cases. % % * Use a line width much wider. Try 2pt. Now it is too deliacte to see convninetly without zoom. % % * enlarge the fonts on the axes so they are closer in size to the font-size of the caption. You can delute the x-axix, showing, say 0,50,100,200, 300, 400 % % * Use thicker lines in the axes. Try thickness of 2 points.

注： ① 代码实现计算 $AA(o)$ 和 $AA(o^{\prime })$ 的相关函数，涉及论文公式 12,14,15 ② 函数：exist/logspace/besselj/trapz/num2str/repmat/waitbar/diff/find/struct/semilogx ③ 代码分平面波还是球面波，共2种情况讨论 ④ 公式15涉及的积分，都采用梯形数值方法计算 ⑤ b_long代表 $b_{\parallel }$ ，b_lat代表 $b_{\perp }$ ，gamma其实是 $\rho$，beta其实是 $\eta$，Gk是 $H(\kappa )$

calc_correlation_Cv.m

function [C,X,Y] = calc_correlation_Cv(sz, b, f, theta_D) % Calculates the cross correlation function (matrix) % AS modified from Marina Alterman beta_vec = b.beta; b_lat_vec = b.lat; b_long_vec = b.long; % Grid in Physical Coordinates dx = 1; [X,Y] = meshgrid(-sz*dx:dx:(sz-1)*dx, -sz*dx:dx:(sz-1)*dx); r = sqrt(X.^2+Y.^2); beta_mat = r/(f*theta_D); %assuming small angles. Close to the center of the image. Relative to the middle pixel % Small correction to avoid exterpolation problems beta_mat = min(beta_mat, beta_vec(end)); % Interpolate correlation function on grid b_long_interp = interp1(beta_vec, b_long_vec, beta_mat(:)); b_long_interp_mat = reshape(b_long_interp, size(beta_mat)); b_lat_interp = interp1(beta_vec, b_lat_vec, beta_mat(:)); b_lat_interp_mat = reshape(b_lat_interp, size(beta_mat)); % Covariance Matrix c11 = b_lat_interp_mat - (b_lat_interp_mat - b_long_interp_mat)./(X.^2+Y.^2).*X.^2; c21 = -(b_lat_interp_mat - b_long_interp_mat)./(X.^2+Y.^2).*X.*Y; % c22 = b_lat_interp_mat - (b_lat_interp_mat - b_long_interp_mat)./(X.^2+Y.^2).*Y.^2; c22 = c11'; % Note that Y = X' c12 = c21; % Output C = cat(4, cat(3, c11, c21), cat(3, c12, c22)); o = sz+1; C(o, o, :, :) = eye(2); return %-------------------------------------------------------------------------- % Check C11 = C(:,:,1,1); C22 = C(:,:,2,2); detC = C(:,:,1,1).*C(:,:,2,2)-C(:,:,1,2).*C(:,:,2,1); [sum(C11(:) < 0), sum(C22(:) < 0), sum(detC(:) < 0)]

注： ① 代码实现计算自相关域 $\mathbf{C}(\mathbf{v})$ ，涉及论文公式 17,32,33 ② 函数：meshgrid/interp1/reshape/cat ③ 这几行代码是公式32展开写的结果：

c11 = b_lat_interp_mat - (b_lat_interp_mat - b_long_interp_mat)./(X.^2+Y.^2).*X.^2; c21 = -(b_lat_interp_mat - b_long_interp_mat)./(X.^2+Y.^2).*X.*Y; c22 = c11'; % Note that Y = X' c12 = c21;

且矩阵 $C_{xy}$ 与 $C_{yx}$ 一样，$C_{xx}$ 与 $C_{yy}$ 是转置关系，可以简化计算