More work on msk40.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@6968 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Steven Franke 2016-07-30 00:44:37 +00:00
parent 5aa3657fc6
commit 6e45134550
3 changed files with 401 additions and 4 deletions

View File

@ -192,9 +192,9 @@ subroutine detectmsk144(cbig,n,pchk_file,lines,nmessages,nutc,ntol,t00)
enddo enddo
cc=cc1+cc2 cc=cc1+cc2
dd=abs(cc1)*abs(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)) 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 enddo
ccr=ccr1+ccr2 ccr=ccr1+ccr2
ddr=abs(ccr1)*abs(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)) crmax=maxval(abs(ccr))
ishort=0 ishort=0
if( crmax .gt. cmax ) ishort=1 if( crmax .gt. cmax ) ishort=1
if( ishort .eq. 1 ) write(*,*) 'Short message detected'
! Find 6 largest peaks ! Find 6 largest peaks
do ipk=1, 6 do ipk=1, 6
iloc=maxloc(abs(cc)) iloc=maxloc(abs(cc))

396
lib/detectmsk40.f90 Normal file
View File

@ -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

View File

@ -83,7 +83,7 @@ subroutine genmsk144(msg0,ichk,msgsent,i4tone,itype,pchk_file)
if(message(1:1).eq.'<') then if(message(1:1).eq.'<') then
call genmsk40(message,msgsent,ichk,i4tone,itype) call genmsk40(message,msgsent,ichk,i4tone,itype)
if(itype.lt.0) go to 999 if(itype.lt.0) go to 999
i4tone(39)=-40 i4tone(41)=-40
go to 999 go to 999
endif endif