First steps toward realtime short ping decoder for MSK40.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7139 ab8295b8-cf94-4d9e-aec4-7959e3be5d79

lib/msk40spd.f90

@@ -0,0 +1,194 @@
+subroutine msk40spd(cbig,n,ntol,nsuccess,msgreceived,fc,fret,tret)
+! msk40 short-ping-decoder
+
+  use timer_module, only: timer
+
+  parameter (NSPM=240, MAXSTEPS=100, NFFT=NSPM, MAXCAND=5, NPATTERNS=6)
+  character*22 msgreceived
+  complex cbig(n)
+  complex cdat(3*NSPM)                    !Analytic signal
+  complex c(NSPM)
+  complex ct(NSPM)
+  complex ctmp(NFFT)                  
+  integer, dimension(1) :: iloc
+  integer indices(MAXSTEPS)
+  integer npkloc(10)
+  integer navpatterns(3,NPATTERNS)
+  integer navmask(3)
+  integer nstart(MAXCAND)
+  logical ismask(NFFT)
+  real detmet(-2:MAXSTEPS+3)
+  real detmet2(-2:MAXSTEPS+3)
+  real detfer(MAXSTEPS)
+  real rcw(12)
+  real ferrs(MAXCAND)
+  real snrs(MAXCAND)
+  real tonespec(NFFT)
+  real tpat(NPATTERNS)
+  real*8 dt, df, fs, pi, twopi
+  logical first
+  data first/.true./
+  data navpatterns/ &
+       0,1,0, &
+       1,0,0, &
+       0,0,1, &
+       1,1,0, &
+       0,1,1, &
+       1,1,1/
+  data tpat/1.5,0.5,2.5,1.0,2.0,1.5/
+
+  save df,first,fs,pi,twopi,dt,tframe,rcw
+
+  if(first) then
+     nmatchedfilter=1
+! define half-sine pulse and raised-cosine edge window
+     pi=4d0*datan(1d0)
+     twopi=8d0*datan(1d0)
+     fs=12000.0
+     dt=1.0/fs
+     df=fs/NFFT
+     tframe=NSPM/fs
+
+     do i=1,12
+       angle=(i-1)*pi/12.0
+       rcw(i)=(1-cos(angle))/2
+     enddo
+
+     first=.false.
+  endif
+
+
+  ! fill the detmet, detferr arrays
+  nstep=(n-NSPM)/60  ! 20ms/4=5ms steps
+  detmet=0
+  detmet2=0
+  detfer=-999.99
+  nfhi=2*(fc+500)
+  nflo=2*(fc-500)
+  ihlo=(nfhi-2*ntol)/df+1
+  ihhi=(nfhi+2*ntol)/df+1
+  illo=(nflo-2*ntol)/df+1
+  ilhi=(nflo+2*ntol)/df+1
+  i2000=nflo/df+1
+  i4000=nfhi/df+1
+  do istp=1,nstep
+    ns=1+60*(istp-1)
+    ne=ns+NSPM-1
+    if( ne .gt. n ) exit
+    ctmp=cmplx(0.0,0.0)
+    ctmp(1:NSPM)=cbig(ns:ne)
+
+! Coarse carrier frequency sync - seek tones at 2000 Hz and 4000 Hz in 
+! squared signal spectrum.
+! search range for coarse frequency error is +/- 100 Hz
+
+    ctmp=ctmp**2
+    ctmp(1:12)=ctmp(1:12)*rcw
+    ctmp(NSPM-11:NSPM)=ctmp(NSPM-11:NSPM)*rcw(12:1:-1)
+    call four2a(ctmp,NFFT,1,-1,1)
+    tonespec=abs(ctmp)**2
+
+    ismask=.false.
+    ismask(ihlo:ihhi)=.true.  ! high tone search window
+    iloc=maxloc(tonespec,ismask)
+    ihpk=iloc(1)
+    deltah=-real( (ctmp(ihpk-1)-ctmp(ihpk+1)) / (2*ctmp(ihpk)-ctmp(ihpk-1)-ctmp(ihpk+1)) )
+    ah=tonespec(ihpk)
+    ahavp=(sum(tonespec,ismask)-ah)/count(ismask)
+    trath=ah/(ahavp+0.01)
+    ismask=.false.
+    ismask(illo:ilhi)=.true.   ! window for low tone
+    iloc=maxloc(tonespec,ismask)
+    ilpk=iloc(1)
+    deltal=-real( (ctmp(ilpk-1)-ctmp(ilpk+1)) / (2*ctmp(ilpk)-ctmp(ilpk-1)-ctmp(ilpk+1)) )
+    al=tonespec(ilpk)
+    alavp=(sum(tonespec,ismask)-al)/count(ismask)
+    tratl=al/(alavp+0.01)
+    fdiff=(ihpk+deltah-ilpk-deltal)*df
+    ferrh=(ihpk+deltah-i4000)*df/2.0
+    ferrl=(ilpk+deltal-i2000)*df/2.0
+    if( ah .ge. al ) then
+      ferr=ferrh
+    else
+      ferr=ferrl
+    endif
+    detmet(istp)=max(ah,al)
+    detmet2(istp)=max(trath,tratl)
+    detfer(istp)=ferr
+  enddo  ! end of detection-metric and frequency error estimation loop
+
+  call indexx(detmet(1:nstep),nstep,indices) !find median of detection metric vector
+  xmed=detmet(indices(nstep/4))
+  detmet=detmet/xmed ! noise floor of detection metric is 1.0
+  ndet=0
+
+  do ip=1,MAXCAND ! Find candidates
+    iloc=maxloc(detmet(1:nstep))
+    il=iloc(1)
+    if( (detmet(il) .lt. 3.0) ) exit 
+    if( abs(detfer(il)) .le. ntol ) then 
+      ndet=ndet+1
+      nstart(ndet)=1+(il-1)*60+1
+      ferrs(ndet)=detfer(il)
+      snrs(ndet)=12.0*log10(detmet(il))/2-9.0
+    endif
+    detmet(il)=0.0
+  enddo
+
+  if( ndet .lt. 3 ) then  
+    do ip=1,MAXCAND-ndet ! Find candidates
+      iloc=maxloc(detmet2(1:nstep))
+      il=iloc(1)
+      if( (detmet2(il) .lt. 12.0) ) exit 
+      if( abs(detfer(il)) .le. ntol ) then 
+        ndet=ndet+1
+        nstart(ndet)=1+(il-1)*60+1
+        ferrs(ndet)=detfer(il)
+        snrs(ndet)=12.0*log10(detmet2(il))/2-9.0
+      endif
+      detmet2(il)=0.0
+    enddo
+  endif
+
+  nsuccess=0
+  msgreceived=' '
+  npeaks=2
+  ntol0=29
+  deltaf=7.2
+  do icand=1,ndet  ! Try to sync/demod/decode each candidate.
+    ib=max(1,nstart(icand)-NSPM)
+    ie=ib-1+3*NSPM
+    if( ie .gt. n ) then
+      ie=n
+      ib=ie-3*NSPM+1
+    endif
+    cdat=cbig(ib:ie) 
+    fo=fc+ferrs(icand)
+    do iav=1,NPATTERNS
+      navmask=navpatterns(1:3,iav) 
+      call msk40sync(cdat,3,ntol0,deltaf,navmask,npeaks,fo,fest,npkloc,nsyncsuccess,c)
+
+      if( nsyncsuccess .eq. 0 ) cycle
+
+      do ipk=1,npeaks
+        do is=1,3
+          ic0=npkloc(ipk)
+          if( is.eq.2) ic0=max(1,ic0-1)
+          if( is.eq.3) ic0=min(NSPM,ic0+1)
+          ct=cshift(c,ic0-1)
+          call msk40decodeframe(ct,msgreceived,ndecodesuccess)
+
+          if( ndecodesuccess .gt. 0 ) then
+            tret=(nstart(icand)+NSPM/2)/fs
+            fret=fest
+!write(*,*) icand, iav, ipk, is, tret, fret, msgreceived
+            nsuccess=1
+            return
+          endif 
+        enddo
+      enddo
+    enddo
+  enddo       ! candidate loop
+
+  return
+end subroutine msk40spd

lib/msk40sync.f90

@@ -0,0 +1,113 @@
+subroutine msk40sync(cdat,nframes,ntol,delf,navmask,npeaks,fc,fest,   &
+     npklocs,nsuccess,c)
+
+  !$ use omp_lib
+
+  parameter (NSPM=240)
+  complex cdat(NSPM*nframes)
+  complex cdat2(NSPM*nframes,8)
+  complex c(NSPM)                    !Coherently averaged complex data
+  complex cs(NSPM,8)
+  complex cb(42)                     !Complex waveform for sync word 
+
+!  integer*8 count0,count1,clkfreq
+  integer s8r(8)
+  integer iloc(1)
+  integer npklocs(npeaks)
+  integer navmask(nframes)                 ! defines which frames to average
+
+  real cbi(42),cbq(42)
+  real pkamps(npeaks)
+  real xcc(0:NSPM-1)
+  real xccs(0:NSPM-1,8)
+  real xm(8)
+  real bf(8)
+  real pp(12)                        !Half-sine pulse shape
+  logical first
+  data first/.true./
+  data s8r/0,1,0,0,1,1,1,0/
+  save first,cb,fs,pi,twopi,dt,s8,pp
+
+!  call system_clock(count0,clkfreq)
+  if(first) then
+     pi=4.0*atan(1.0)
+     twopi=8.0*atan(1.0)
+     fs=12000.0
+     dt=1.0/fs
+
+     do i=1,12                       !Define half-sine pulse
+       angle=(i-1)*pi/12.0
+       pp(i)=sin(angle)
+     enddo
+
+! Define the sync word waveforms
+     s8r=2*s8r-1  
+     cbq(1:6)=pp(7:12)*s8r(1)
+     cbq(7:18)=pp*s8r(3)
+     cbq(19:30)=pp*s8r(5)
+     cbq(31:42)=pp*s8r(7)
+     cbi(1:12)=pp*s8r(2)
+     cbi(13:24)=pp*s8r(4)
+     cbi(25:36)=pp*s8r(6)
+     cbi(37:42)=pp(1:6)*s8r(8)
+     cb=cmplx(cbi,cbq)
+
+     first=.false.
+  endif
+
+  nfreqs=2*nint(ntol/delf) + 1
+  xm=0.0
+  bf=0.0
+  nthreads=1
+  !$ nthreads=min(8,int(OMP_GET_MAX_THREADS(),4))
+  nstep=nfreqs/nthreads
+
+  !$OMP PARALLEL NUM_THREADS(nthreads) PRIVATE(id,if1,if2)
+  id=1
+  !$ id=OMP_GET_THREAD_NUM() + 1            !Thread id = 1,2,...
+  if1=-nint(ntol/delf) + (id-1)*nstep
+  if2=if1+nstep-1
+  if(id.eq.nthreads) if2=nint(ntol/delf)
+  call msk40_freq_search(cdat,fc,if1,if2,delf,nframes,navmask,cb,    &
+       cdat2(1,id),xm(id),bf(id),cs(1,id),xccs(1,id))
+!  write(73,3002) id,if1,if2,nfreqs,nthreads,bf(id),xm(id)
+!3002 format(5i5,2f10.3)
+  !$OMP END PARALLEL
+
+  xmax=xm(1)
+  fest=fc+bf(1)
+  c=cs(1:NSPM,1)
+  xcc=xccs(0:NSPM-1,1)
+  if(nthreads.gt.1) then
+     do i=2,nthreads
+        if(xm(i).gt.xmax) then
+           xmax=xm(i)
+           fest=fc+bf(i)
+           c=cs(1:NSPM,i)
+           xcc=xccs(0:NSPM-1,i)
+        endif
+     enddo
+  endif
+
+! Find npeaks largest peaks
+  do ipk=1,npeaks
+     iloc=maxloc(xcc)
+     ic2=iloc(1)
+     npklocs(ipk)=ic2
+     pkamps(ipk)=xcc(ic2-1)
+     xcc(max(0,ic2-7):min(NSPM-1,ic2+7))=0.0
+  enddo
+
+  if( xmax .lt. 0.7 ) then
+    nsuccess=0
+  else
+    nsuccess=1
+  endif
+
+!  call system_clock(count1,clkfreq)
+!  t=float(count1-count0)/clkfreq
+!  write(72,3001) nfreqs,OMP_GET_MAX_THREADS(),nthreads,t
+!3001 format(3i6,f8.3)
+
+  return
+end subroutine msk40sync