# GJB6289命中精度纠偏型估计算法的matlab代码实现

• JOHNCLOCK
了解作者
• matlab
开发工具
• 4KB
文件大小
• zip
文件格式
• 0
收藏次数
• 5 积分
下载积分
• 0
下载次数
• 2022-10-07 21:08
上传日期
GJB6289命中精度纠偏型估计算法的matlab代码实现，可以用于正态/非正态命中数据的估计与评定
GJB6289-2008.zip
• func_L1jk.m
383B
• getRhat_GJB6289.m
7.1KB
• func_L2ijk.m
412B

function [Rhat, Rchat] =getRhat_GJB6289(Xi, Zi) if size(Xi, 1) ~= 1 error('input Xi must be row vec'); elseif size(Zi, 1) ~= 1 error('input Xi must be row vec'); end N = numel(Xi); mu_x_hat = 1/N*sum(Xi); mu_z_hat = 1/N*sum(Zi); sig_x_hat = sqrt(1/(N-1)*sum((Xi-mu_x_hat).^2)); sig_z_hat = sqrt(1/(N-1)*sum((Zi-mu_z_hat).^2)); p_hat = sum((Xi-mu_x_hat).*(Zi-mu_z_hat))/sqrt(sum((Xi-mu_x_hat).^2).*sum((Zi-mu_z_hat).^2)); bxi = 1/2*atan(2*p_hat*sig_x_hat*sig_z_hat... /(sig_x_hat^2-sig_z_hat^2)); mu_xp_hat = mu_x_hat*cos(bxi) + mu_z_hat*sin(bxi); mu_zp_hat = mu_z_hat*cos(bxi) - mu_x_hat*sin(bxi); C1 = (sig_z_hat*cos(bxi))^2 + (sig_x_hat*sin(bxi))^2 -... p_hat*sig_x_hat*sig_z_hat*sin(2*bxi); C2 = (sig_x_hat*cos(bxi))^2 + (sig_z_hat*sin(bxi))^2 +... p_hat*sig_x_hat*sig_z_hat*sin(2*bxi); sig_xp_hat = sig_x_hat*sig_z_hat*sqrt((1-p_hat^2)/C1); sig_zp_hat = sig_x_hat*sig_z_hat*sqrt((1-p_hat^2)/C2); u1 = mu_x_hat; u2 = mu_z_hat; s1 = sig_x_hat; s2 = sig_z_hat; [~, pval] = corrcoef(Xi', Zi'); if pval(1,2) < 0.1 u1 = mu_xp_hat; u2 = mu_zp_hat; s1 = sig_xp_hat; s2 = sig_zp_hat; end a = 1/4*(1./s2^2 - 1./s1^2); b = 1/4*(1./s2^2 + 1./s1^2); c = 1/(s1*s2)*exp(-1/2*((u1/s1).^2 + (u2/s2).^2)); sR = sqrt(u1^2+u2^2); func = @(x,y) exp(-1/2 * ((x-u1).^2./s1.^2 + (y-u2).^2./s2^2)); CEP_a = 0; Rhat = sR; while (abs(CEP_a - 0.5) > 1e-6) ry = @(x) sqrt(Rhat.^2-x.^2); rym = @(x) -sqrt(Rhat.^2-x.^2); CEP_a = 1/(2*pi*s1*s2)*integral2(func, -Rhat, Rhat, rym, ry); if CEP_a > 0.5 Rhat = Rhat - Rhat*abs(CEP_a - 0.5); elseif CEP_a < 0.5 Rhat = Rhat + Rhat*abs(CEP_a - 0.5); end end theta1 = u1; theta2 = u2; theta3 = s1^2; theta4 = s2^2; v1 = theta3; v2 = theta4; v3 = 2*theta3^2; v4 = 2*theta4^2; w1 = 0; w2 = 0; w3 = 8*theta3^3; w4 = 8*theta4^3; R = Rhat; % psi(x) = log(x), frist order; % psi(x) = x, C-F L1_00 = func_L1jk(0, 0, R, u1, s1, u2, s2); T = L1_00*c*R; L1_20 = func_L1jk(2, 0, R, u1, s1, u2, s2); L1_02 = func_L1jk(0, 2, R, u1, s1, u2, s2); L1_01 = func_L1jk(0, 1, R, u1, s1, u2, s2); L1_10 = func_L1jk(1, 0, R, u1, s1, u2, s2); L2_210 = func_L2ijk(2, 1, 0, R, u1, s1, u2, s2); Tsup1 = L2_210/s1^2; L2_201 = func_L2ijk(2, 0, 1, R, u1, s1, u2, s2); Tsup2 = L2_201/s2^2; L2_320 = func_L2ijk(3, 2, 0, R, u1, s1, u2, s2); Tsup3 = 1/(2*s1^4)*(L2_320 - 2*u1*L2_210); L2_302 = func_L2ijk(3, 0, 2, R, u1, s1, u2, s2); Tsup4 = 1/(2*s2^4)*(L2_302 - 2*u2*L2_201); L2_311 = func_L2ijk(3, 1, 1, R, u1, s1, u2, s2); L2_430 = func_L2ijk(4, 3, 0, R, u1, s1, u2, s2); L2_403 = func_L2ijk(4, 0, 3, R, u1, s1, u2, s2); L2_412 = func_L2ijk(4, 1, 2, R, u1, s1, u2, s2); L2_421 = func_L2ijk(4, 2, 1, R, u1, s1, u2, s2); L2_540 = func_L2ijk(5, 4, 0, R, u1, s1, u2, s2); L2_504 = func_L2ijk(5, 0, 4, R, u1, s1, u2, s2); L2_522 = func_L2ijk(5, 2, 2, R, u1, s1, u2, s2); c1 = -u1*c/(s1^2); c2 = -u2*c/(s2^2); c3 = c/(2*s1^2)*(u1^2/s1^2 - 1); c4 = c/(2*s2^2)*(u2^2/s2^2 - 1); c11 = c/s1^2*(u1^2/s1^2 - 1); c12 = u1*u2*c/((s1*s2)^2); c13 = u1*c/(2*(s1^4))*(3-u1^2/s1^2); c14 = u1*c/(2*(s1*s2)^2)*(1-u2^2/s2^2); c21 = c12; c22 = c/s2^2*(u2^2/s2^2 - 1); c23 = u2*c/(2*s1^2*s2^2)*(1-u1^2/s1^2); c24 = u2*c/(2*s2^4)*(3-u2^2/s2^2); c31 = c13; c32 = c23; c33 = c/(4*s1^4)*(3 - 6*u1^2/s1^2 + u1^4/s1^4); c34 = c/(4*s1^2*s2^2)*(1 - u1^2/s1^2 - u2^2/s2^2 + u1^2*u2^2/(s1^2*s2^2)); c41 = c14; c42 = c24; c43 = c34; c44 = c/(4*s2^4)*(3- 6*u2^2/s2^2 + u2^4/s2^4); cij = [c11, c12, c13, c14; c21, c22, c23, c24; c31, c32, c33, c34; c41, c42, c43, c44]; R1 = -1/T * (c1/(2*c)+c*Tsup1); R2 = -1/T * (c2/(2*c)+c*Tsup2); R3 = -1/T * (c3/(2*c)+c*Tsup3); R4 = -1/T * (c4/(2*c)+c*Tsup4); T1 = (c1/c + R1/R - R1*R/s2^2)*T + 4*a*c*R1*R^2*L1_20 + c*R/s1^2*(R+u1*R1)*L1_10 + c*u2/s2^2*R1*R*L1_01; T2 = (c2/c + R2/R - R2*R/s2^2)*T + 4*a*c*R2*R^2*L1_20 + c*R/s2^2*(R+u2*R2)*L1_01 + c*u1/s1^2*R2*R*L1_10; T3 = (c3/c + R3/R - R3*R/s2^2)*T + 4*a*c*R3*R^2*L1_20 + c*R^3/(2*s1^4)*L1_20 + c*u1*R/s1^2*(R3-R/s1^2)*L1_10 ... +c*u2/s2^2*R3*R*L1_01; T4 = (c4/c + R4/R - R4*R/s2^2)*T + 4*a*c*R4*R^2*L1_20 + c*R^3/(2*s2^4)*L1_02 + c*u2*R/s2^2*(R4-R/s2^2)*L1_01 ... +c*u1/s1^2*R4*R*L1_10; Tk = [T1; T2; T3; T4]; T1s1 = R1*R^2/s1^2 * L1_10 + 1/(s1^4)*L2_320; T2s1 = R2*R^2/s1^2 * L1_10 + 1/(s1^2*s2^2)*L2_311; T3s1 = R3*R^2/s1^2 * L1_10 - 1/(s1^4)*L2_210 + 1/(2*s1^6)*L2_430 - u1/s1^6*L2_320; T4s1 = R4*R^2/s1^2 * L1_10 + 1/(2*s1^2*s2^4)*L2_412 - u2/(s1^2*s2^4)*L2_311; T1s2 = R1*R^2/s2^2*L1_01 + 1/(s1^2*s2^2)*L2_311; T2s2 = R2*R^2/s2^2*L1_01 + 1/s2^4*L2_302; T3s2 = R3*R^2/s2^2*L1_01 + 1/(2*s1^4*s2^2)*L2_421 - u1/(s1^4*s2^2)*L2_311; T4s2 = R4*R^2/s2^2*L1_01 - 1/(s2^4)*L2_201 + 1/(2*s2^6)*L2_403 - u2/(s2^6)*L2_302; T1s3 = R1*R^3/(2*s1^4)*L1_20 - u1*R1*R^2/(s1^4)*L1_10 - 1/(s1^4)*L2_210 ... +1/(2*s1^6)*(L2_430 - 2*u1*L2_320); T2s3 = R2*R^3/(2*s1^4)*L1_20 - u1*R2*R^2/(s1^4)*L1_10 + 1/(2*s1^4*s2^2) ... *(L2_421 - 2*u1*L2_311); T3s3 = R3*R^3/(2*s1^4)*L1_20 - u1*R3*R^2/(s1^4)*L1_10 - 1/(s1^6)*(L2_320 - 2*u1*L2_210) ... +1/(4*s1^8)*(L2_540 - 4*u1*L2_430 + 4*u1^2*L2_320); T4s3 = R4*R^3/(2*s1^4)*L1_20 - u1*R4*R^2/(s1^4)*L1_10 + 1/(4*s1^4*s2^4)... *(L2_522 - 2*u2*L2_421 - 2*u1*L2_412 + 4*u1*u2*L2_311); T1s4 = R1*R^3/(2*s2^4)*L1_02 - u2*R1*R^2/(s2^4)*L1_01 + 1/(2*s1^2*s2^4)*L2_412 ... -u2/(s1^2*s2^4)*L2_311; T2s4 = R2*R^3/(2*s2^4)*L1_02 - u2*R2*R^2/(s2^4)*L1_01 - 1/(s2^4)*L2_201 + 1/(2*s2^6)*(L2_403 - 2*u2*L2_302); T3s4 = R3*R^3/(2*s2^4)*L1_02 - u2*R3*R^2/(s2^4)*L1_01 + 1/(4*s1^4*s2^4)... *(L2_522 - 2*u1*L2_412 - 2*u2*L2_421 + 4*u1*u2*L2_311); T4s4 = R4*R^3/(2*s2^4)*L1_02 - u2*R4*R^2/(s2^4)*L1_01 - 1/(s2^6)*(L2_302 - 2*u2*L2_201) ... +1/(4*s2^8)*(L2_504 - 4*u2*L2_403 + 4*u2^2*L2_302); ck = [c1;c2;c3;c4]; Tjsi = [T1s1, T1s2, T1s3, T1s4; T2s1, T2s2, T2s3, T2s4; T3s1, T3s2, T3s3, T3s4; T4s1, T4s2, T4s3, T4s4]; Tsupk = [Tsup1; Tsup2; Tsup3; Tsup4]; Rk = [R1; R2; R3; R4]; Rij = zeros(4, 4); psi_ij = zeros(4, 4); for i = 1:4 for j = 1:4 Rij(i,j) = -1/T * (1/2*(cij(i,j)/c - ck(i)*ck(j)/c^2) + Rk(i)*Tk(j) + ... Tsupk(i)*ck(j) + c*Tjsi(j, i)); % psi_ij(i, j) = Rij(i,j)/R - 1/R^2 * Rk(i) * Rk(j); psi_ij(i, j) = Rij(i,j); end end psi1 = R1; psi2 = R2; psi3 = R3; psi4 = R4; vk = [v1;v2;v3;v4]; wk = [w1;w2;w3;w4]; Hpsi = 1./sqrt(psi1^2*v1 + psi2^2*v2 + psi3^2*v3 + psi4^2*v4); psik = [psi1; psi2; psi3; psi4]; % pvptheta_ij = diag([1,1,4*theta3,4*theta4]); pvptheta_ij = [0,0,1,0; 0,0,0,1; 0,0,4*theta3,0; 0,0,0,4*theta4]; tmp = zeros(4,4); for i = 1:4 for j = 1:4 tmp(i,j) = pvptheta_ij(i,j)*psik(i)^2 + 2*psik(i)*vk(i)*psi_ij(i,j); end end pHpsi_pthetak = -Hpsi^3/2 * sum(tmp, 1)'; Ghatpsi1_p1 = 1/6*Hpsi^3*sum(wk.*(psik.^3)); tmp1 = 0; for i = 1:4 for j = 1:4 tmp1 = tmp1 + vk(i)*vk(j)*psik(i)*psi_ij(i,j)*psik(j); end end Ghatpsi1_p2 = 1/2*Hpsi^3*tmp1; Ghatpsi1_p3 = -1/2*Hpsi*sum(vk.*(diag(psi_ij))); Ghatpsi1 = Ghatpsi1_p1 + Ghatpsi1_p2 + Ghatpsi1_p3; Ghatpsi3_p1 = -1/6*Hpsi^5*sum(wk.*(psik.^3)); Ghatpsi3_p2 = -1/2*Hpsi^5*tmp1; Ghatpsi3_p3 = -Hpsi^2*sum(vk.*pHpsi_pthetak.*psik); Ghatpsi3 = Ghatpsi3_p1 + Ghatpsi3_p2 + Ghatpsi3_p3; Rchat = Rhat + 1/N*(Ghatpsi3/Hpsi^3 + Ghatpsi1/Hpsi); end

相关推荐
• libiconv-1.1.tar.gz
字符集转换程序
• VisualC++.rar
vc++数字图像处理 ,是一本很不错的介绍数字图像方面的书籍,这里有本书的全部源码
• sharewareluncher.zip
执行和去除共享软件日期限制的程序
• opctkit.rar
opc client 的开发工具
• 一个XML留言本源代码(C#).rar
用C#编写的XML源程序