From 6e4513455081c4c1e1f7a102dcea9dd9289ca896 Mon Sep 17 00:00:00 2001 From: Steven Franke Date: Sat, 30 Jul 2016 00:44:37 +0000 Subject: [PATCH] More work on msk40. git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@6968 ab8295b8-cf94-4d9e-aec4-7959e3be5d79 --- lib/detectmsk144.f90 | 7 +- lib/detectmsk40.f90 | 396 +++++++++++++++++++++++++++++++++++++++++++ lib/genmsk144.f90 | 2 +- 3 files changed, 401 insertions(+), 4 deletions(-) create mode 100644 lib/detectmsk40.f90 diff --git a/lib/detectmsk144.f90 b/lib/detectmsk144.f90 index 5c004d41d..e78ba2d65 100644 --- a/lib/detectmsk144.f90 +++ b/lib/detectmsk144.f90 @@ -192,9 +192,9 @@ subroutine detectmsk144(cbig,n,pchk_file,lines,nmessages,nutc,ntol,t00) enddo cc=cc1+cc2 dd=abs(cc1)*abs(cc2) - do i=1,NPTS-(32*6+41) + do i=1,NPTS-(40*6+41) ccr1(i)=sum(cdat(i:i+41)*conjg(cbr)) - ccr2(i)=sum(cdat(i+32*6:i+32*6+41)*conjg(cbr)) + ccr2(i)=sum(cdat(i+40*6:i+40*6+41)*conjg(cbr)) enddo ccr=ccr1+ccr2 ddr=abs(ccr1)*abs(ccr2) @@ -202,7 +202,8 @@ subroutine detectmsk144(cbig,n,pchk_file,lines,nmessages,nutc,ntol,t00) crmax=maxval(abs(ccr)) ishort=0 if( crmax .gt. cmax ) ishort=1 - +if( ishort .eq. 1 ) write(*,*) 'Short message detected' + ! Find 6 largest peaks do ipk=1, 6 iloc=maxloc(abs(cc)) diff --git a/lib/detectmsk40.f90 b/lib/detectmsk40.f90 new file mode 100644 index 000000000..362d8fb9a --- /dev/null +++ b/lib/detectmsk40.f90 @@ -0,0 +1,396 @@ +subroutine detectmsk40(cbig,n,mycall,partnercall,lines,nmessages,nutc,ntol,t00) + use timer_module, only: timer + + parameter (NSPM=240, NPTS=3*NSPM, MAXSTEPS=7500, NFFT=3*NSPM, MAXCAND=10) + character*4 rpt(0:63) + character*6 mycall,partnercall + character*22 hashmsg,msgreceived + character*80 lines(100) + character*40 pchk_file,gen_file + complex cbig(n) + complex cdat(NPTS) !Analytic signal + complex cdat2(NPTS) + complex c(NSPM) + complex ctmp(NFFT) + complex cb(42) !Complex waveform for sync word + complex cbr(42) !Complex waveform for reversed sync word + complex cfac,cca,ccb + complex ccr(NPTS) + complex ccr1(NPTS) + complex ccr2(NPTS) + complex bb(6) + integer s8(8),s8r(8),hardbits(40) + integer, dimension(1) :: iloc + integer nhashes(0:63) + integer indices(MAXSTEPS) + integer ipeaks(10) + integer*1 decoded(16) + logical ismask(NFFT) + real cbi(42),cbq(42) + real detmet(-2:MAXSTEPS+3) + real detfer(MAXSTEPS) + real rcw(12) + real ddr(NPTS) + real ferrs(MAXCAND) + real pp(12) !Half-sine pulse shape + real snrs(MAXCAND) + real softbits(40) + real times(MAXCAND) + real tonespec(NFFT) + real*8 dt, df, fs, pi, twopi + real*8 lratio(32) + logical first + data first/.true./ + data s8/0,1,1,1,0,0,1,0/ + data s8r/1,0,1,1,0,0,0,1/ + save df,first,cb,cbr,fs,nhashes,pi,twopi,dt,s8,s8r,rcw,pp,nmatchedfilter + + 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 + + do i=1,12 + angle=(i-1)*pi/12.0 + pp(i)=sin(angle) + rcw(i)=(1-cos(angle))/2 + enddo + +! define the sync word waveforms + s8=2*s8-1 + cbq(1:6)=pp(7:12)*s8(1) + cbq(7:18)=pp*s8(3) + cbq(19:30)=pp*s8(5) + cbq(31:42)=pp*s8(7) + cbi(1:12)=pp*s8(2) + cbi(13:24)=pp*s8(4) + cbi(25:36)=pp*s8(6) + cbi(37:42)=pp(1:6)*s8(8) + cb=cmplx(cbi,cbq) + 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) + cbr=cmplx(cbi,cbq) + + do i=0,30 + if( i.lt.5 ) then + write(rpt(i),'(a1,i2.2,a1)') '-',abs(i-5) + write(rpt(i+31),'(a2,i2.2,a1)') 'R-',abs(i-5) + else + write(rpt(i),'(a1,i2.2,a1)') '+',i-5 + write(rpt(i+31),'(a2,i2.2,a1)') 'R+',i-5 + endif + enddo + rpt(62)='RRR ' + rpt(63)='73 ' + + do i=0,63 + hashmsg=trim(mycall)//' '//trim(partnercall)//' '//rpt(i) + call fmtmsg(hashmsg,iz) + call hash(hashmsg,22,ihash) + nhashes(i)=iand(ihash,1023) + enddo + + first=.false. + endif + +! Temporarily hardwire filenames and init on every call + pchk_file="peg-32-16-reg3.pchk" + i=index(pchk_file,".pchk") + gen_file=pchk_file(1:i-1)//".gen" + call init_ldpc(trim(pchk_file)//char(0),trim(gen_file)//char(0)) + +! Fill the detmet, detferr arrays + nstepsize=48 ! 4ms steps + nstep=(n-NPTS)/nstepsize + detmet=0 + detmax=-999.99 + detfer=-999.99 + do istp=1,nstep + ns=1+nstepsize*(istp-1) + ne=ns+NPTS-1 + if( ne .gt. n ) exit + ctmp=cmplx(0.0,0.0) + ctmp(1:NPTS)=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(NPTS-11:NPTS)=ctmp(NPTS-11:NPTS)*rcw(12:1:-1) + call four2a(ctmp,NFFT,1,-1,1) + tonespec=abs(ctmp)**2 + + ihlo=(4000-2*ntol)/df+1 + ihhi=(4000+2*ntol)/df+1 + 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) + illo=(2000-2*ntol)/df+1 + ilhi=(2000+2*ntol)/df+1 + 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) + fdiff=(ihpk+deltah-ilpk-deltal)*df + i2000=2000/df+1 + i4000=4000/df+1 + 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) + detfer(istp)=ferr +! write(*,*) istp,ilpk,ihpk,ah,al + 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 i=1,nstep +!write(77,*) i,detmet(i),detfer(i) +!enddo + + do ip=1,MAXCAND ! find candidates + iloc=maxloc(detmet(1:nstep)) + il=iloc(1) + if( (detmet(il) .lt. 4.2) ) exit + if( abs(detfer(il)) .le. ntol ) then + ndet=ndet+1 + times(ndet)=((il-1)*nstepsize+NPTS/2)*dt + ferrs(ndet)=detfer(il) + snrs(ndet)=12.0*log10(detmet(il)-1)/2-8.0 + endif + detmet(max(1,il-3):min(nstep,il+3))=0.0 +! detmet(il)=0.0 + enddo + +! ndet=15 + do ip=1,ndet +! times(ip)=ip+0.012 +! snrs(ip)=-3.0 +! ferrs(ip)=0.0 + write(*,'(i5,f7.2,f7.2,f7.2)') ip,times(ip),snrs(ip),ferrs(ip) + enddo + + nmessages=0 + lines=char(0) + + do ip=1,ndet !run through the candidates and try to sync/demod/decode + imid=times(ip)*fs + if( imid .lt. NPTS/2 ) imid=NPTS/2 + if( imid .gt. n-NPTS/2 ) imid=n-NPTS/2 + t0=times(ip) + t00 + cdat=cbig(imid-NPTS/2+1:imid+NPTS/2) + ferr=ferrs(ip) + nsnr=nint(snrs(ip)) + if( nsnr .lt. -5 ) nsnr=-5 + if( nsnr .gt. 25 ) nsnr=25 + +! remove coarse freq error + call tweak1(cdat,NPTS,-(1500+ferr),cdat) + +! attempt frame synchronization +! correlate with sync word waveforms + ccr=0 + ccr1=0 + ccr2=0 + do i=1,NPTS-(40*6+41) + ccr1(i)=sum(cdat(i:i+41)*conjg(cbr)) + ccr2(i)=sum(cdat(i+40*6:i+40*6+41)*conjg(cbr)) + enddo + ccr=ccr1+ccr2 + ddr=abs(ccr1)*abs(ccr2) + crmax=maxval(abs(ccr)) + +do i=1,NPTS +write(14,*) i,abs(ccr(i)),ddr(i) +enddo + + +! Find 6 largest peaks + do ipk=1,6 + iloc=maxloc(abs(ccr)) + ic1=iloc(1) + iloc=maxloc(ddr) + ic2=iloc(1) + ipeaks(ipk)=ic1 + ccr(max(1,ic1-7):min(NPTS-40*6-41,ic1+7))=0.0 + enddo +do i=1,6 +write(*,*) i,ipeaks(i) +enddo + do ipk=1,6 + +! we want ic to be the index of the first sample of the frame + ic0=ipeaks(ipk) + +! fine adjustment of sync index + do i=1,6 + if( ic0+11+NSPM .le. NPTS ) then + bb(i) = sum( ( cdat(ic0+i-1+6:ic0+i-1+6+NSPM:6) * conjg( cdat(ic0+i-1:ic0+i-1+NSPM:6) ) )**2 ) + else + bb(i) = sum( ( cdat(ic0+i-1+6:NPTS:6) * conjg( cdat(ic0+i-1:NPTS-6:6) ) )**2 ) + endif + enddo + iloc=maxloc(abs(bb)) + ibb=iloc(1) + bba=abs(bb(ibb)) + bbp=atan2(-imag(bb(ibb)),-real(bb(ibb)))/(2*twopi*6*dt) + if( ibb .le. 3 ) ibb=ibb-1 + if( ibb .gt. 3 ) ibb=ibb-7 + + do id=1,1 ! slicer dither. + if( id .eq. 1 ) is=0 + if( id .eq. 2 ) is=-1 + if( id .eq. 3 ) is=1 + +! Adjust frame index to place peak of bb at desired lag + ic=ic0+ibb+is + if( ic .lt. 1 ) ic=ic+NSPM + +! Estimate fine frequency error. + cca=sum(cdat(ic:ic+41)*conjg(cb)) + if( ic+40*6+41 .le. NPTS ) then + ccb=sum(cdat(ic+40*6:ic+40*6+41)*conjg(cb)) + cfac=ccb*conjg(cca) + ferr2=atan2(imag(cfac),real(cfac))/(twopi*40*6*dt) + else + ccb=sum(cdat(ic-40*6:ic-40*6+41)*conjg(cb)) + cfac=cca*conjg(ccb) + ferr2=atan2(imag(cfac),real(cfac))/(twopi*40*6*dt) + endif + +! Final estimate of the carrier frequency - returned to the calling program + fest=1500+ferr+ferr2 + + do idf=0,0 ! frequency jitter + if( idf .eq. 0 ) then + deltaf=0.0 + elseif( mod(idf,2) .eq. 0 ) then + deltaf=2*idf + else + deltaf=-2*(idf+1) + endif + +! Remove fine frequency error + call tweak1(cdat,NPTS,-(ferr2+deltaf),cdat2) + +! place the beginning of frame at index NSPM+1 + cdat2=cshift(cdat2,ic-(NSPM+1)) + + do iav=1,4 ! Frame averaging patterns + if( iav .eq. 1 ) then + c=cdat2(NSPM+1:2*NSPM) + elseif( iav .eq. 2 ) then + c=cdat2(1:NSPM)+cdat2(NSPM+1:2*NSPM) + elseif( iav .eq. 3 ) then + c=cdat2(NSPM+1:2*NSPM)+cdat2(2*NSPM+1:3*NSPM) + elseif( iav .eq. 4 ) then + c=cdat2(1:NSPM)+cdat2(NSPM+1:2*NSPM)+cdat2(2*NSPM+1:3*NSPM) + endif + +! Estimate final frequency error and carrier phase. + cca=sum(c(1:1+41)*conjg(cb)) + phase0=atan2(imag(cca),real(cca)) + + do ipha=1,1 + if( ipha.eq.2 ) phase0=phase0-20*pi/180.0 + if( ipha.eq.3 ) phase0=phase0+20*pi/180.0 + +! Remove phase error - want constellation rotated so that sample points lie on I/Q axes + cfac=cmplx(cos(phase0),sin(phase0)) + c=c*conjg(cfac) + + if( nmatchedfilter .eq. 0 ) then + do i=1, 20 + softbits(2*i-1)=imag(c(1+(i-1)*12)) + softbits(2*i)=real(c(7+(i-1)*12)) + enddo + else ! matched filter + softbits(1)=sum(imag(c(1:6))*pp(7:12))+sum(imag(c(NSPM-5:NSPM))*pp(1:6)) + softbits(2)=sum(real(c(1:12))*pp) + do i=2,20 + softbits(2*i-1)=sum(imag(c(1+(i-1)*12-6:1+(i-1)*12+5))*pp) + softbits(2*i)=sum(real(c(7+(i-1)*12-6:7+(i-1)*12+5))*pp) + enddo + endif + + hardbits=0 ! use sync word hard error weight to decide whether to send to decoder + do i=1, 40 + if( softbits(i) .ge. 0.0 ) then + hardbits(i)=1 + endif + enddo + nbadsync1=(8-sum( (2*hardbits(1:8)-1)*s8r ) )/2 + nbadsync=nbadsync1 + if( nbadsync .gt. 3 ) cycle +write(*,*) ip,id,nbadsync + ! normalize the softsymbols before submitting to decoder + sav=sum(softbits)/40 + s2av=sum(softbits*softbits)/40 + ssig=sqrt(s2av-sav*sav) + softbits=softbits/ssig + + sigma=0.75 + lratio(1:32)=exp(2.0*softbits(9:40)/(sigma*sigma)) + + max_iterations=5 + max_dither=1 + call ldpc_decode(lratio,decoded,max_iterations,niterations,max_dither,ndither) + + nhashflag=0 + if( niterations .ge. 0 ) then + imsg=0 + do i=1,16 + imsg=ishft(imsg,1)+iand(1,decoded(i)) + enddo + nrxrpt=iand(imsg,63) + nrxhash=(imsg-nrxrpt)/64 +write(*,*) niterations,decoded,nrxrpt,nrxhash,nhashes(nrxrpt) + if( nrxhash .eq. nhashes(nrxrpt) ) then + nhashflag=1 + goto 999 + endif + endif + enddo ! phase loop + enddo ! frame averaging loop + enddo ! frequency dithering loop + enddo ! slicer dither loop + enddo ! time-sync correlation-peak loop + enddo ! candidate loop +999 continue + msgreceived=' ' + if( nhashflag.eq.1 ) then + nmessages=1 + write(msgreceived,'(a1,a,1x,a,a1,1x,a4)') "<",trim(mycall), & + trim(partnercall),">",rpt(nrxrpt) + write(lines(nmessages),1020) nutc,nsnr,t0,nint(fest),msgreceived +1020 format(i6.6,i4,f5.1,i5,' & ',a22) + + endif + + return +end subroutine detectmsk40 diff --git a/lib/genmsk144.f90 b/lib/genmsk144.f90 index 7ea737d42..bac51cda9 100644 --- a/lib/genmsk144.f90 +++ b/lib/genmsk144.f90 @@ -83,7 +83,7 @@ subroutine genmsk144(msg0,ichk,msgsent,i4tone,itype,pchk_file) if(message(1:1).eq.'<') then call genmsk40(message,msgsent,ichk,i4tone,itype) if(itype.lt.0) go to 999 - i4tone(39)=-40 + i4tone(41)=-40 go to 999 endif