1subroutine msk144spd(cbig,n,ntol,nsuccess,msgreceived,fc,fret,tret,navg,ct,   &
2  softbits)
3
4! MSK144 short-ping-decoder
5
6  use packjt77
7  use timer_module, only: timer
8
9  parameter (NSPM=864, MAXSTEPS=100, NFFT=NSPM, MAXCAND=5, NPATTERNS=6)
10  character*37 msgreceived
11  complex cbig(n)
12  complex cdat(3*NSPM)                    !Analytic signal
13  complex c(NSPM)
14  complex ct(NSPM)
15  complex ctmp(NFFT)
16  integer, dimension(1) :: iloc
17  integer indices(MAXSTEPS)
18  integer npkloc(10)
19  integer navpatterns(3,NPATTERNS)
20  integer navmask(3)
21  integer nstart(MAXCAND)
22  logical ismask(NFFT)
23  real detmet(-2:MAXSTEPS+3)
24  real detmet2(-2:MAXSTEPS+3)
25  real detfer(MAXSTEPS)
26  real rcw(12)
27  real ferrs(MAXCAND)
28  real snrs(MAXCAND)
29  real softbits(144)
30  real tonespec(NFFT)
31  real tpat(NPATTERNS)
32  real*8 dt, df, fs, pi, twopi
33  logical first
34  data first/.true./
35  data navpatterns/ &
36       0,1,0, &
37       1,0,0, &
38       0,0,1, &
39       1,1,0, &
40       0,1,1, &
41       1,1,1/
42  data tpat/1.5,0.5,2.5,1.0,2.0,1.5/
43
44  save df,first,fs,pi,twopi,dt,tframe,rcw
45
46  if(first) then
47     nmatchedfilter=1
48! define half-sine pulse and raised-cosine edge window
49     pi=4d0*datan(1d0)
50     twopi=8d0*datan(1d0)
51     fs=12000.0
52     dt=1.0/fs
53     df=fs/NFFT
54     tframe=NSPM/fs
55
56     do i=1,12
57       angle=(i-1)*pi/12.0
58       rcw(i)=(1-cos(angle))/2
59     enddo
60
61     first=.false.
62  endif
63
64  ! fill the detmet, detferr arrays
65  nstep=(n-NSPM)/216  ! 72ms/4=18ms steps
66  detmet=0
67  detmet2=0
68  detfer=-999.99
69  nfhi=2*(fc+500)
70  nflo=2*(fc-500)
71  ihlo=nint((nfhi-2*ntol)/df) + 1
72  ihhi=nint((nfhi+2*ntol)/df) + 1
73  illo=nint((nflo-2*ntol)/df) + 1
74  ilhi=nint((nflo+2*ntol)/df) + 1
75  i2000=nint(nflo/df) + 1
76  i4000=nint(nfhi/df) + 1
77  do istp=1,nstep
78    ns=1+216*(istp-1)
79    ne=ns+NSPM-1
80    if( ne .gt. n ) exit
81    ctmp=cmplx(0.0,0.0)
82    ctmp(1:NSPM)=cbig(ns:ne)
83
84! Coarse carrier frequency sync - seek tones at 2000 Hz and 4000 Hz in
85! squared signal spectrum.
86
87    ctmp=ctmp**2
88    ctmp(1:12)=ctmp(1:12)*rcw
89    ctmp(NSPM-11:NSPM)=ctmp(NSPM-11:NSPM)*rcw(12:1:-1)
90    call four2a(ctmp,NFFT,1,-1,1)
91    tonespec=abs(ctmp)**2
92
93    ismask=.false.
94    ismask(ihlo:ihhi)=.true.  ! high tone search window
95    iloc=maxloc(tonespec,ismask)
96    ihpk=iloc(1)
97    deltah=-real( (ctmp(ihpk-1)-ctmp(ihpk+1)) / (2*ctmp(ihpk)-ctmp(ihpk-1)-ctmp(ihpk+1)) )
98    ah=tonespec(ihpk)
99    ahavp=(sum(tonespec,ismask)-ah)/count(ismask)
100    trath=ah/(ahavp+0.01)
101    ismask=.false.
102    ismask(illo:ilhi)=.true.   ! window for low tone
103    iloc=maxloc(tonespec,ismask)
104    ilpk=iloc(1)
105    deltal=-real( (ctmp(ilpk-1)-ctmp(ilpk+1)) / (2*ctmp(ilpk)-ctmp(ilpk-1)-ctmp(ilpk+1)) )
106    al=tonespec(ilpk)
107    alavp=(sum(tonespec,ismask)-al)/count(ismask)
108    tratl=al/(alavp+0.01)
109    fdiff=(ihpk+deltah-ilpk-deltal)*df
110    ferrh=(ihpk+deltah-i4000)*df/2.0
111    ferrl=(ilpk+deltal-i2000)*df/2.0
112    if( ah .ge. al ) then
113      ferr=ferrh
114    else
115      ferr=ferrl
116    endif
117    detmet(istp)=max(ah,al)
118    detmet2(istp)=max(trath,tratl)
119    detfer(istp)=ferr
120  enddo  ! end of detection-metric and frequency error estimation loop
121
122  call indexx(detmet(1:nstep),nstep,indices) !find median of detection metric vector
123  xmed=detmet(indices(nstep/4))
124  detmet=detmet/xmed ! noise floor of detection metric is 1.0
125  ndet=0
126
127  do ip=1,MAXCAND ! Find candidates
128    iloc=maxloc(detmet(1:nstep))
129    il=iloc(1)
130    if( (detmet(il) .lt. 3.0) ) exit
131    if( abs(detfer(il)) .le. ntol ) then
132      ndet=ndet+1
133      nstart(ndet)=1+(il-1)*216+1
134      ferrs(ndet)=detfer(il)
135      snrs(ndet)=12.0*log10(detmet(il))/2-9.0
136    endif
137    detmet(il)=0.0
138  enddo
139
140  if( ndet .lt. 3 ) then
141    do ip=1,MAXCAND-ndet ! Find candidates
142      iloc=maxloc(detmet2(1:nstep))
143      il=iloc(1)
144      if( (detmet2(il) .lt. 12.0) ) exit
145      if( abs(detfer(il)) .le. ntol ) then
146        ndet=ndet+1
147        nstart(ndet)=1+(il-1)*216+1
148        ferrs(ndet)=detfer(il)
149        snrs(ndet)=12.0*log10(detmet2(il))/2-9.0
150      endif
151      detmet2(il)=0.0
152    enddo
153  endif
154
155  nsuccess=0
156  msgreceived=' '
157  npeaks=2
158  ntol0=8
159  deltaf=2.0
160  do icand=1,ndet  ! Try to sync/demod/decode each candidate.
161    ib=max(1,nstart(icand)-NSPM)
162    ie=ib-1+3*NSPM
163    if( ie .gt. n ) then
164      ie=n
165      ib=ie-3*NSPM+1
166    endif
167    cdat=cbig(ib:ie)
168    fo=fc+ferrs(icand)
169    do iav=1,NPATTERNS
170      navmask=navpatterns(1:3,iav)
171      call msk144sync(cdat,3,ntol0,deltaf,navmask,npeaks,fo,fest,npkloc,   &
172           nsyncsuccess,xmax,c)
173
174      if( nsyncsuccess .eq. 0 ) cycle
175
176      do ipk=1,npeaks
177        do is=1,3
178          ic0=npkloc(ipk)
179          if( is.eq.2) ic0=max(1,ic0-1)
180          if( is.eq.3) ic0=min(NSPM,ic0+1)
181          ct=cshift(c,ic0-1)
182          call msk144decodeframe(ct,softbits,msgreceived,ndecodesuccess)
183          if( ndecodesuccess .gt. 0 ) then
184            tret=(nstart(icand)+NSPM/2)/fs
185            fret=fest
186            navg=sum(navmask)
187            nsuccess=1
188            return
189          endif
190        enddo
191      enddo
192    enddo
193  enddo       ! candidate loop
194
195  return
196end subroutine msk144spd
197