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源程序
    • VB6MultiThread.rar
      VB多线程源代码
    • ScrollBitmap_demo.zip
      Displaying a large bitmap file on a dialog box, in its original size, is quite difficult in the VC++ environment. However, it is possible to display a large bitmap to a predefined area of the dialog by using the StretchBlt( ) function.The major disadvantage of this is that the clarity of the image will be lost. Check out this article for displaying large bitmaps into the desired area of your dialog box in its original size with a scrolling technique used to show the entire bitmap. 滚动显示位图 在VC++环境下,在一个对话框中显示一个原始尺寸的大小的位图文件相当是困难的。然而,通过使用 StretchBlt()函数一个给定的区域显示一个大的位图是可能的。主要的缺点是图像将会失真。看了这篇通过卷动技术显示整个位图技术的文章,你将能够以它的原始尺寸在给定对话框的区域内显示一个大位图。 来源: http://www.codeguru.com/bitmap/ScrollBitmap.html