pd_swerling4.m
上传用户:szahd2008
上传日期:2020-09-25
资源大小:1275k
文件大小:2k
源码类别:

传真(Fax)编程

开发平台:

Matlab

  1. function pd = pd_swerling4 (nfa, np, snrbar)
  2. % This function is used to calculate the probability of detection
  3. % for Swerling 4 targets.
  4. format long
  5. snrbar = 10.0^(snrbar/10.);
  6. eps = 0.00000001;
  7. delmax = .00001;
  8. delta =10000.;
  9. % Calculate the threshold Vt
  10. pfa =  np * log(2) / nfa;
  11. sqrtpfa = sqrt(-log10(pfa));
  12. sqrtnp = sqrt(np); 
  13. vt0 = np - sqrtnp + 2.3 * sqrtpfa * (sqrtpfa + sqrtnp - 1.0);
  14. vt = vt0;
  15. while (abs(delta) >= vt0)
  16.    igf = incomplete_gamma(vt0,np);
  17.    num = 0.5^(np/nfa) - igf;
  18.    temp = (np-1) * log(vt0+eps) - vt0 - factor(np-1);
  19.    deno = exp(temp);
  20.    vt = vt0 + (num / (deno+eps));
  21.    delta = abs(vt - vt0) * 10000.0; 
  22.    vt0 = vt;
  23. end
  24. h8 = snrbar /2.0;
  25. beta = 1.0 + h8;
  26. beta2 = 2.0 * beta^2 - 1.0;
  27. beta3 = 2.0 * beta^3;
  28. if (np >= 50)
  29.    temp1 = 2.0 * beta -1;
  30.    omegabar = sqrt(np * temp1);
  31.    c3 = (beta3 - 1.) / 3.0 / beta2 / omegabar;
  32.    c4 = (beta3 * beta3 - 1.0) / 4. / np /beta2 /beta2;;
  33.    c6 = c3 * c3 /2.0;
  34.    V = (vt - np * (1.0 + snrbar)) / omegabar;
  35.    Vsqr = V *V;
  36.    val1 = exp(-Vsqr / 2.0) / sqrt( 2.0 * pi);
  37.    val2 = c3 * (V^2 -1.0) + c4 * V * (3.0 - V^2) - ... 
  38.       c6 * V * (V^4 - 10. * V^2 + 15.0);
  39.    q = 0.5 * erfc (V/sqrt(2.0));
  40.    pd =  q - val1 * val2;
  41.    return
  42. else
  43.    snr = 1.0;
  44.    gamma0 = incomplete_gamma(vt/beta,np);
  45.    a1 = (vt / beta)^np / (exp(factor(np)) * exp(vt/beta));
  46.    sum = gamma0;
  47.    for i = 1:1:np
  48.       temp1 = 1;
  49.       if (i == 1)
  50.          ai = a1;
  51.       else
  52.          ai = (vt / beta) * a1 / (np + i -1);
  53.       end
  54.       a1 = ai;
  55.       gammai = gamma0 - ai;
  56.       gamma0 = gammai;
  57.       a1 = ai;
  58.       for ii = 1:1:i
  59.          temp1 = temp1 * (np + 1 - ii);
  60.       end
  61.       term = (snrbar /2.0)^i * gammai * temp1 / exp(factor(i));
  62.       sum = sum + term;
  63.    end
  64.    pd = 1.0 - sum / beta^np;
  65. end
  66. pd = max(pd,0.);
  67. return
  68.