radar2d4.f90

基于matlab的探地雷达信号处理开源程序

  do j = 1,nz
    vmean(j) = 0.0
    do i = 1,ntr
      vmean(j) = vmean(j)+ v(j,i)/2.
    end do
    vmean(j) = vmean(j)/float(ntr)
  end do
end subroutine
! --------------------------------------------------
subroutine mean_q(qm, q, ntr, nz)
!*****  computes mean quality factor
implicit none
integer :: nz,ntr,i,j
real(4) :: qm(nz), q(nz, ntr)
! - - - begin - - -
do j = 1,nz
  qm(j) = 0.0
  do i = 1,ntr
    qm(j) = qm(j)+ q(j,i)
  end do
  qm(j) = qm(j)/float(ntr)
end do
end subroutine
! --------------------------------------------------
function source(tet1,tet2,tet)  result (output_4)
!*****  Computes Antenna radiation pattern
implicit none
! - - - arg types - - -
real(4) :: output_4                                       
! - - - local declarations - - -
real(8),parameter :: pi = 3.141592653589793d0
integer :: tet1,tet2
real(4) :: tet,src,alf1,alf2
! - - - begin - - -
alf1 = pi*tet1/180
alf2 = pi*tet2/180
if (tet <= alf1) then
   output_4 = (1.0)
   return
else
   if (tet <= alf2) then
      src = (0.42-0.5*cos((1+(tet-alf1)/(alf2-alf1))*pi)+  &
             0.08*cos((1+(tet-alf1)/(alf2-alf1))*2*pi))
      output_4 = src
      return
   else
      output_4 = 0.0
      return
   end if
end if
return
end function
! --------------------------------------------------
function npfft(nmin,nmax) result (output_4)
!*************************************************************************
!Return optimal n between nmin and nmax for prime factor fft.
!*************************************************************************
!Input:
!nmin		lower bound on returned value (see notes below)
!nmax		desired (but not guaranteed) upper bound on returned value
!Returned:	valid n for prime factor fft
!Notes:
!The returned n will be composed of mutually prime factors from
!the set {2,3,4,5,7,8,9,11,13,16}.  Because n cannot exceed
!720720 = 5*7*9*11*13*16, 720720 is returned if nmin exceeds 720720.
!If nmin does not exceed 720720, then the returned n will not be 
!less than nmin.  The optimal n is chosen to minimize the estimated
!cost of performing the fft, while satisfying the constraint, if
!possible, that n not exceed nmax.
!*************************************************************************
!Author:   Dave Hale, Colorado School of Mines
!F90 code: Andreas Tzanis, University of Athens
!*************************************************************************
implicit none
! - - - arg types - - -
integer :: output_4                                       
integer :: nmin,nmax
! - - - local declarations - - -
integer :: i,j
integer*4,parameter :: ntab=240
integer*4 :: nn(ntab)
real*4    :: c(ntab)
data nn/  &
     1,     2,     3,     4,     5,     6,     7,     8, &
     9,    10,    11,    12,    13,    14,    15,    16, &
    18,    20,    21,    22,    24,    26,    28,    30, &
    33,    35,    36,    39,    40,    42,    44,    45, &
    48,    52,    55,    56,    60,    63,    65,    66, &
    70,    72,    77,    78,    80,    84,    88,    90, &
    91,    99,   104,   105,   110,   112,   117,   120, &
   126,   130,   132,   140,   143,   144,   154,   156, &
   165,   168,   176,   180,   182,   195,   198,   208, &
   210,   220,   231,   234,   240,   252,   260,   264, &
   273,   280,   286,   308,   312,   315,   330,   336, &
   360,   364,   385,   390,   396,   420,   429,   440, &
   455,   462,   468,   495,   504,   520,   528,   546, &
   560,   572,   585,   616,   624,   630,   660,   693, &
   715,   720,   728,   770,   780,   792,   819,   840, &
   858,   880,   910,   924,   936,   990,  1001,  1008, &
  1040,  1092,  1144,  1155,  1170,  1232,  1260,  1287, &
  1320,  1365,  1386,  1430,  1456,  1540,  1560,  1584, &
  1638,  1680,  1716,  1820,  1848,  1872,  1980,  2002, &
  2145,  2184,  2288,  2310,  2340,  2520,  2574,  2640, &
  2730,  2772,  2860,  3003,  3080,  3120,  3276,  3432, &
  3465,  3640,  3696,  3960,  4004,  4095,  4290,  4368, &
  4620,  4680,  5005,  5040,  5148,  5460,  5544,  5720, &
  6006,  6160,  6435,  6552,  6864,  6930,  7280,  7920, &
  8008,  8190,  8580,  9009,  9240,  9360, 10010, 10296, &
 10920, 11088, 11440, 12012, 12870, 13104, 13860, 15015, &
 16016, 16380, 17160, 18018, 18480, 20020, 20592, 21840, &
 24024, 25740, 27720, 30030, 32760, 34320, 36036, 40040, &
 45045, 48048, 51480, 55440, 60060, 65520, 72072, 80080, &
 90090,102960,120120,144144,180180,240240,360360,720720/
data c/  &
 0.000052,0.000061,0.000030,0.000053,0.000066,0.000067,0.000071,0.000062, &
 0.000079,0.000080,0.000052,0.000069,0.000103,0.000123,0.000050,0.000086, &
 0.000108,0.000101,0.000098,0.000135,0.000090,0.000165,0.000084,0.000132, &
 0.000158,0.000138,0.000147,0.000207,0.000156,0.000158,0.000176,0.000171, &
 0.000185,0.000227,0.000242,0.000194,0.000215,0.000233,0.000288,0.000271, &
 0.000248,0.000247,0.000285,0.000395,0.000285,0.000209,0.000332,0.000321, &
 0.000372,0.000400,0.000391,0.000358,0.000440,0.000367,0.000494,0.000413, &
 0.000424,0.000549,0.000480,0.000450,0.000637,0.000497,0.000590,0.000626, &
 0.000654,0.000536,0.000656,0.000611,0.000730,0.000839,0.000786,0.000835, &
 0.000751,0.000826,0.000926,0.000991,0.000852,0.000820,0.001053,0.000987, &
 0.001152,0.000952,0.001299,0.001155,0.001270,0.001156,0.001397,0.001173, &
 0.001259,0.001471,0.001569,0.001767,0.001552,0.001516,0.002015,0.001748, &
 0.001988,0.001921,0.001956,0.002106,0.001769,0.002196,0.002127,0.002454, &
 0.002099,0.002632,0.002665,0.002397,0.002711,0.002496,0.002812,0.002949, &
 0.003571,0.002783,0.003060,0.003392,0.003553,0.003198,0.003726,0.003234, &
 0.004354,0.003800,0.004304,0.003975,0.004123,0.004517,0.005066,0.003902, &
 0.004785,0.005017,0.005599,0.005380,0.005730,0.005323,0.005112,0.006658, &
 0.005974,0.006781,0.006413,0.007622,0.006679,0.007032,0.007538,0.007126, &
 0.007979,0.007225,0.008961,0.008818,0.008427,0.009004,0.009398,0.010830, &
 0.012010,0.010586,0.012058,0.011673,0.011700,0.011062,0.014313,0.013021, &
 0.014606,0.013216,0.015789,0.016988,0.014911,0.016393,0.016741,0.018821, &
 0.018138,0.018892,0.018634,0.020216,0.022455,0.022523,0.026087,0.023474, &
 0.024590,0.025641,0.030303,0.025253,0.030364,0.031250,0.029412,0.034404, &
 0.037500,0.034091,0.040214,0.037221,0.042735,0.040214,0.042980,0.045872, &
 0.049505,0.049834,0.055762,0.057034,0.054945,0.056818,0.066667,0.065502, &
 0.068182,0.065217,0.075000,0.078534,0.087719,0.081081,0.084270,0.102740, &
 0.106383,0.105634,0.119048,0.123967,0.119048,0.137615,0.140187,0.154639, &
 0.168539,0.180723,0.180723,0.220588,0.241935,0.254237,0.254237,0.288462, &
 0.357143,0.357143,0.384615,0.384615,0.454545,0.517241,0.576923,0.625000, &
 0.833333,0.789474,1.153846,1.153846,1.875000,2.500000,3.750000,7.5000000/
! - - - begin - - -
i = 1
do while (i <= ntab .AND. nn(i) <= nmin)
   i = i+1
   j = i
   do while (j<=ntab .AND. nn(j) <= nmax)
      j = j+1
	  IF (c(j) < c(i) ) i = j
   enddo
enddo
output_4 = nn(i)
return
end function
! --------------------------------------------------
      subroutine fft(a,b,ntot,n,nspan,isn)
!  multivariate complex fourier transform, computed in place
!    using mixed-radix fast fourier transform algorithm.
!  by r. c. singleton, stanford research institute, sept. 1968
!  arrays a and b originally hold the real and imaginary
!    components of the data, and return the real and
!    imaginary components of the resulting fourier coefficients.
!  multivariate data is indexed according to the fortran
!    array element successor function, without limit
!    on the number of implied multiple subscripts.
!    the subroutine is called once for each variate.
!    the calls for a multivariate transform may be in any order.
!  ntot is the total number of complex data values.
!  n is the dimension of the current variable.
!  nspan/n is the spacing of consecutive data values
!    while indexing the current variable.
!  the sign of isn determines the sign of the complex
!    exponential, and the magnitude of isn is normally one.
!  a tri-variate transform with a(n1,n2,n3), b(n1,n2,n3)
!    is computed by
!      call fft(a,b,n1*n2*n3,n1,n1,1)
!      call fft(a,b,n1*n2*n3,n2,n1*n2,1)
!      call fft(a,b,n1*n2*n3,n3,n1*n2*n3,1)
!  for a single-variate transform,
!    ntot = n = nspan = (number of complex data values), e.g.
!      call fft(a,b,n,n,n,1)
!  the data can alternatively be stored in a single complex array c
!    in standard fortran fashion, i.e. alternating real and imaginary
!    parts. then with most fortran compilers, the complex array c can
!    be equivalenced to a real array a, the magnitude of isn changed
!    to two to give correct indexing increment, and a(1) and a(2) used
!    to pass the initial addresses for the sequences of real and
!    imaginary values, e.g.
!       complex c(ntot)
!       real    a(2*ntot)
!       equivalence (c(1),a(1))
!       call fft(a(1),a(2),ntot,n,nspan,2)
!  arrays at(maxf), ck(maxf), bt(maxf), sk(maxf), and np(maxp)
!    are used for temporary storage.  if the available storage
!    is insufficient, the program is terminated by a stop.
!    maxf must be .ge. the maximum prime factor of n.
!    maxp must be .gt. the number of prime factors of n.
!    in addition, if the square-free portion k of n has two or
!    more prime factors, then maxp must be .ge. k-1.
      dimension a(1),b(1)
!  array storage in nfac for a maximum of 15 prime factors of n.
!  if n has more than one square-free factor, the product of the
!    square-free factors must be .le. 210
      dimension nfac(11),np(209)
!  array storage for maximum prime factor of 23
      dimension at(23),ck(23),bt(23),sk(23)
      equivalence (i,ii)
!  the following two constants should agree with the array dimensions.
      maxp=209
!
! Date: Wed, 9 Aug 1995 09:38:49 -0400
! From: [e-mail]
      maxf=23
!
      if(n .lt. 2) return
      inc=isn
      c72=0.30901699437494742
      s72=0.95105651629515357
      s120=0.86602540378443865
      rad=6.2831853071796
      if(isn .ge. 0) go to 10
      s72=-s72
      s120=-s120
      rad=-rad
      inc=-inc
   10 nt=inc*ntot
      ks=inc*nspan
      kspan=ks
      nn=nt-inc
      jc=ks/n
      radf=rad*float(jc)*0.5
      i=0
      jf=0
!  determine the factors of n
      m=0
      k=n
      go to 20
   15 m=m+1
      nfac(m)=4
      k=k/16
   20 if(k-(k/16)*16 .eq. 0) go to 15
      j=3
      jj=9
      go to 30
   25 m=m+1
      nfac(m)=j
      k=k/jj
   30 if(mod(k,jj) .eq. 0) go to 25
      j=j+2
      jj=j**2
      if(jj .le. k) go to 30
      if(k .gt. 4) go to 40
      kt=m
      nfac(m+1)=k
      if(k .ne. 1) m=m+1
      go to 80
   40 if(k-(k/4)*4 .ne. 0) go to 50
      m=m+1
      nfac(m)=2
      k=k/4
   50 kt=m
      j=2
   60 if(mod(k,j) .ne. 0) go to 70
      m=m+1
      nfac(m)=j
      k=k/j
   70 j=((j+1)/2)*2+1
      if(j .le. k) go to 60
   80 if(kt .eq. 0) go to 100
      j=kt
   90 m=m+1
      nfac(m)=nfac(j)
      j=j-1
      if(j .ne. 0) go to 90
!  compute fourier transform
  100 sd=radf/float(kspan)
      cd=2.0*sin(sd)**2
      sd=sin(sd+sd)
      kk=1
      i=i+1
      if(nfac(i) .ne. 2) go to 400
!  transform for factor of 2 (including rotation factor)
      kspan=kspan/2
      k1=kspan+2
  210 k2=kk+kspan
      ak=a(k2)
      bk=b(k2)
      a(k2)=a(kk)-ak
      b(k2)=b(kk)-bk
      a(kk)=a(kk)+ak
      b(kk)=b(kk)+bk
      kk=k2+kspan
      if(kk .le. nn) go to 210
      kk=kk-nn
      if(kk .le. jc) go to 210
      if(kk .gt. kspan) go to 800
  220 c1=1.0-cd
      s1=sd
  230 k2=kk+kspan
      ak=a(kk)-a(k2)
      bk=b(kk)-b(k2)
      a(kk)=a(kk)+a(k2)
      b(kk)=b(kk)+b(k2)
      a(k2)=c1*ak-s1*bk
      b(k2)=s1*ak+c1*bk
      kk=k2+kspan
      if(kk .lt. nt) go to 230
      k2=kk-nt
      c1=-c1
      kk=k1-k2
      if(kk .gt. k2) go to 230
      ak=c1-(cd*c1+sd*s1)
      s1=(sd*c1-cd*s1)+s1
      c1=2.0-(ak**2+s1**2)
      s1=c1*s1
      c1=c1*ak
      kk=kk+jc
      if(kk .lt. k2) go to 230
      k1=k1+inc+inc
      kk=(k1-kspan)/2+jc
      if(kk .le. jc+jc) go to 220
      go to 100
!  transform for factor of 3 (optional code)
  320 k1=kk+kspan
      k2=k1+kspan
      ak=a(kk)
      bk=b(kk)
      aj=a(k1)+a(k2)
      bj=b(k1)+b(k2)
      a(kk)=ak+aj
      b(kk)=bk+bj
      ak=-0.5*aj+ak
      bk=-0.5*bj+bk
      aj=(a(k1)-a(k2))*s120
      bj=(b(k1)-b(k2))*s120
      a(k1)=ak-bj
      b(k1)=bk+aj
      a(k2)=ak+bj
      b(k2)=bk-aj
      kk=k2+kspan
      if(kk .lt. nn) go to 320
      kk=kk-nn
      if(kk .le. kspan) go to 320
      go to 700
!  transform for factor of 4
  400 if(nfac(i) .ne. 4) go to 600
      kspnn=kspan
      kspan=kspan/4
  410 c1=1.0
      s1=0
  420 k1=kk+kspan
      k2=k1+kspan
      k3=k2+kspan