Remove some dead code.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7432 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Steven Franke 2016-12-30 20:26:59 +00:00
parent f30311f324
commit 2c2d31ea8d
7 changed files with 3 additions and 1219 deletions

View File

@ -329,8 +329,6 @@ set (wsjt_FSRCS
lib/deg2grid.f90
lib/degrade_snr.f90
lib/demod64a.f90
lib/detectmsk144.f90
lib/detectmsk40.f90
lib/determ.f90
lib/downsam9.f90
lib/encode232.f90
@ -409,7 +407,6 @@ set (wsjt_FSRCS
lib/moondopjpl.f90
lib/morse.f90
lib/move.f90
lib/msk144d.f90
lib/msk144d2.f90
lib/msk40decodeframe.f90
lib/msk144decodeframe.f90
@ -425,7 +422,6 @@ set (wsjt_FSRCS
lib/msk144sim.f90
lib/mskrtd.f90
lib/options.f90
lib/opdetmsk144.f90
lib/packjt.f90
lib/pctile.f90
lib/peakdt9.f90
@ -1074,9 +1070,6 @@ target_link_libraries (ldpcsim144 wsjt_fort wsjt_cxx)
add_executable (msk144sim lib/msk144sim.f90 wsjtx.rc)
target_link_libraries (msk144sim wsjt_fort wsjt_cxx)
add_executable (msk144d lib/msk144d.f90 wsjtx.rc)
target_link_libraries (msk144d wsjt_fort wsjt_cxx)
add_executable (msk144d2 lib/msk144d2.f90 wsjtx.rc)
target_link_libraries (msk144d2 wsjt_fort wsjt_cxx)

View File

@ -1,425 +0,0 @@
subroutine detectmsk144(cbig,n,lines,nmessages,nutc,ntol,t00)
use timer_module, only: timer
parameter (NSPM=864, NPTS=3*NSPM, MAXSTEPS=1700, NFFT=NSPM, MAXCAND=16)
character*22 msgreceived,allmessages(20)
character*80 lines(100)
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 cc(NPTS)
complex ccr(NPTS)
complex cc1(NPTS)
complex cc2(NPTS)
complex ccr1(NPTS)
complex ccr2(NPTS)
complex bb(6)
integer s8(8),s8r(8),hardbits(144)
integer, dimension(1) :: iloc
integer*1 decoded(80)
integer indices(MAXSTEPS)
integer ipeaks(10)
logical ismask(NFFT)
real cbi(42),cbq(42)
real detmet(-2:MAXSTEPS+3)
real detmet2(-2:MAXSTEPS+3)
real detfer(MAXSTEPS)
real rcw(12)
real dd(NPTS)
! real ddr(NPTS)
real ferrs(MAXCAND)
real pp(12) !Half-sine pulse shape
real snrs(MAXCAND)
real times(MAXCAND)
real tonespec(NFFT)
real*8 dt, df, fs, pi, twopi
real softbits(144)
real*8 lratio(128)
real llr(128)
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,fs,pi,twopi,dt,s8,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)
first=.false.
endif
! fill the detmet, detferr arrays
nstep=(n-NPTS)/216 ! 72ms/4=18ms steps
detmet=0
detmet2=0
detfer=-999.99
do istp=1,nstep
ns=1+216*(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
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)
ahavp=(sum(tonespec,ismask)-ah)/count(ismask)
trath=ah/(ahavp+0.01)
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)
alavp=(sum(tonespec,ismask)-al)/count(ismask)
tratl=al/(alavp+0.01)
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)
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. 4.0) ) exit
if( abs(detfer(il)) .le. ntol ) then
ndet=ndet+1
times(ndet)=((il-1)*216+NSPM/2)*dt
ferrs(ndet)=detfer(il)
snrs(ndet)=12.0*log10(detmet(il))/2-9.0
endif
! detmet(max(1,il-1):min(nstep,il+1))=0.0
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
times(ndet)=((il-1)*216+NSPM/2)*dt
ferrs(ndet)=detfer(il)
snrs(ndet)=12.0*log10(detmet2(il))/2-9.0
endif
! detmet2(max(1,il-1):min(nstep,il+1))=0.0
detmet2(il)=0.0
enddo
endif
nmessages=0
allmessages=char(0)
lines=char(0)
nshort=0
do ip=1,ndet ! Try to sync/demod/decode each candidate.
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=2*nint(snrs(ip)/2.0)
if( nsnr .lt. -4 ) nsnr=-4
if( nsnr .gt. 24 ) nsnr=24
! remove coarse freq error - should now be within a few Hz
call tweak1(cdat,NPTS,-(1500+ferr),cdat)
! attempt frame synchronization
! correlate with sync word waveforms
cc=0
ccr=0
cc1=0
cc2=0
ccr1=0
ccr2=0
do i=1,NPTS-(56*6+41)
cc1(i)=sum(cdat(i:i+41)*conjg(cb))
cc2(i)=sum(cdat(i+56*6:i+56*6+41)*conjg(cb))
enddo
cc=cc1+cc2
dd=abs(cc1)*abs(cc2)
! 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)
cmax=maxval(abs(cc))
! crmax=maxval(abs(ccr))
! if( crmax .gt. cmax ) then
! nshort=nshort+1
! endif
! Find 6 largest peaks
do ipk=1, 6
iloc=maxloc(abs(cc))
ic1=iloc(1)
iloc=maxloc(dd)
ic2=iloc(1)
ipeaks(ipk)=ic2
dd(max(1,ic2-7):min(NPTS-56*6-41,ic2+7))=0.0
! ipeaks(ipk)=ic1
! cc(max(1,ic1-7):min(NPTS-56*6-41,ic1+7))=0.0
!write(*,*) ipk,ic2
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,3 ! 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+864
! Estimate fine frequency error.
! Should a larger separation be used when frames are averaged?
cca=sum(cdat(ic:ic+41)*conjg(cb))
if( ic+56*6+41 .le. NPTS ) then
ccb=sum(cdat(ic+56*6:ic+56*6+41)*conjg(cb))
cfac=ccb*conjg(cca)
ferr2=atan2(imag(cfac),real(cfac))/(twopi*56*6*dt)
else
ccb=sum(cdat(ic-88*6:ic-88*6+41)*conjg(cb))
cfac=cca*conjg(ccb)
ferr2=atan2(imag(cfac),real(cfac))/(twopi*88*6*dt)
endif
! Final estimate of the carrier frequency - returned to the calling program
fest=1500+ferr+ferr2
do idf=0,4 ! frequency jitter
if( idf .eq. 0 ) then
deltaf=0.0
elseif( mod(idf,2) .eq. 0 ) then
deltaf=idf
else
deltaf=-(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,8 ! Hopefully we can eliminate some of these after looking at more examples
if( iav .eq. 1 ) then
c=cdat2(NSPM+1:2*NSPM)
elseif( iav .eq. 2 ) then
c=cdat2(NSPM-431:NSPM+432)
c=cshift(c,-432)
elseif( iav .eq. 3 ) then
c=cdat2(2*NSPM-431:2*NSPM+432)
c=cshift(c,-432)
elseif( iav .eq. 4 ) then
c=cdat2(1:NSPM)
elseif( iav .eq. 5 ) then
c=cdat2(2*NSPM+1:NPTS)
elseif( iav .eq. 6 ) then
c=cdat2(1:NSPM)+cdat2(NSPM+1:2*NSPM)
elseif( iav .eq. 7 ) then
c=cdat2(NSPM+1:2*NSPM)+cdat2(2*NSPM+1:NPTS)
elseif( iav .eq. 8 ) then
c=cdat2(1:NSPM)+cdat2(NSPM+1:2*NSPM)+cdat2(2*NSPM+1:NPTS)
endif
! Estimate final frequency error and carrier phase.
cca=sum(c(1:1+41)*conjg(cb))
ccb=sum(c(1+56*6:1+56*6+41)*conjg(cb))
cfac=ccb*conjg(cca)
ffin=atan2(imag(cfac),real(cfac))/(twopi*56*6*dt)
phase0=atan2(imag(cca+ccb),real(cca+ccb))
do ipha=1,1
if( ipha.eq.2 ) phase0=phase0+30*pi/180.0
if( ipha.eq.3 ) phase0=phase0-30*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
! sample to get softsamples
do i=1,72
softbits(2*i-1)=imag(c(1+(i-1)*12))
softbits(2*i)=real(c(7+(i-1)*12))
enddo
else
! matched filter -
! how much mismatch does the RX/TX/analytic filter cause?, how rig (pair) dependent is this loss?
softbits(1)=sum(imag(c(1:6))*pp(7:12))+sum(imag(c(864-5:864))*pp(1:6))
softbits(2)=sum(real(c(1:12))*pp)
do i=2,72
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
! sync word hard error weight is a good discriminator for
! frames that have reasonable probability of decoding
hardbits=0
do i=1,144
if( softbits(i) .ge. 0.0 ) then
hardbits(i)=1
endif
enddo
nbadsync1=(8-sum( (2*hardbits(1:8)-1)*s8 ) )/2
nbadsync2=(8-sum( (2*hardbits(1+56:8+56)-1)*s8 ) )/2
nbadsync=nbadsync1+nbadsync2
if( nbadsync .gt. 4 ) cycle
! normalize the softsymbols before submitting to decoder
sav=sum(softbits)/144
s2av=sum(softbits*softbits)/144
ssig=sqrt(s2av-sav*sav)
softbits=softbits/ssig
sigma=0.75
lratio(1:48)=softbits(9:9+47)
lratio(49:128)=softbits(65:65+80-1)
llr=2.0*lratio/(sigma*sigma)
lratio=exp(2.0*lratio/(sigma*sigma))
max_iterations=10
max_dither=1
call timer('bpdec144 ',0)
call bpdecode144(llr,max_iterations,decoded,niterations)
call timer('bpdec144 ',1)
if( niterations .ge. 0.0 ) then
call extractmessage144(decoded,msgreceived,nhashflag)
if( nhashflag .gt. 0 ) then ! CRCs match, so print it
ndupe=0
do im=1,nmessages
if( allmessages(im) .eq. msgreceived ) ndupe=1
enddo
if( ndupe .eq. 0 ) then
nmessages=nmessages+1
allmessages(nmessages)=msgreceived
write(lines(nmessages),1020) nutc,nsnr,t0,nint(fest),msgreceived
1020 format(i6.6,i4,f5.1,i5,' & ',a22)
endif
goto 999
else
msgreceived=' '
ndither=-99 ! -99 is bad hash flag
! write(78,1001) nutc,t0,nsnr,ip,ipk,is,idf,iav,ipha,deltaf,fest,ferr,ferr2, &
! ffin,bba,bbp,nbadsync1,nbadsync2, &
! phase0,niterations,ndither,msgreceived
endif
endif
enddo ! phase dither
enddo ! frame averaging loop
enddo ! frequency dithering loop
enddo ! sample-time dither loop
enddo ! peak loop - could be made more efficient by working harder to find good peaks
msgreceived=' '
ndither=-98
999 continue
if( nmessages .ge. 1 ) then
! write(78,1001) nutc,t0,nsnr,ip,ipk,is,idf,iav,ipha,deltaf,fest,ferr,ferr2, &
! ffin,bba,bbp,nbadsync1,nbadsync2, &
! phase0,niterations,ndither,msgreceived
! call flush(78)
!1001 format(i6.6,f8.2,i5,i5,i5,i5,i5,i5,i5,f6.2,f8.2,f8.2,f8.2,f8.2,f11.1,f8.2,i5,i5,f8.2,i5,i5,2x,a22)
exit
endif
enddo
return
end subroutine detectmsk144

View File

@ -1,425 +0,0 @@
subroutine detectmsk40(cbig,n,mycall,hiscall,lines,nmessages, &
nutc,ntol,t00)
use timer_module, only: timer
parameter (NSPM=240, NPTS=3*NSPM, MAXSTEPS=7500, NFFT=3*NSPM, MAXCAND=15)
character*4 rpt(0:15)
character*6 mycall,hiscall,mycall0,hiscall0
character*22 hashmsg,msgreceived
character*80 lines(100)
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:15)
integer indices(MAXSTEPS)
integer ipeaks(10)
integer*1 cw(32)
integer*1 decoded(16)
integer*1 testcw(32)
logical ismask(NFFT)
real cbi(42),cbq(42)
real detmet(-2:MAXSTEPS+3)
real detmet2(-2:MAXSTEPS+3)
real detfer(MAXSTEPS)
real rcw(12)
real ddr(NPTS)
real ferrs(MAXCAND)
real llr(32)
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 mycall0/'dummy'/,hiscall0/'dummy'/
data rpt/"-03 ","+00 ","+03 ","+06 ","+10 ","+13 ","+16 ", &
"R-03","R+00","R+03","R+06","R+10","R+13","R+16", &
"RRR ","73 "/
data s8/0,1,1,1,0,0,1,0/
data s8r/1,0,1,1,0,0,0,1/
! codeword for the message <K9AN K1JT> RRR
data testcw/0,1,0,0,0,1,0,1,0,0,1,1,1,1,0,1,0,1,1,1,1,1,0,0,0,0,0,0,1,1,1,0/
save df,first,cb,cbr,fs,pi,twopi,dt,s8,s8r,rcw,pp,nmatchedfilter,rpt,mycall0,hiscall0,ihash
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)
first=.false.
endif
if(mycall.ne.mycall0 .or. hiscall.ne.hiscall0) then
! do i=0,15
! hashmsg=trim(mycall)//' '//trim(hiscall)//' '//rpt(i)
hashmsg=trim(mycall)//' '//trim(hiscall)
call fmtmsg(hashmsg,iz)
call hash(hashmsg,22,ihash)
! nhashes(i)=iand(ihash,4095)
ihash=iand(ihash,4095)
! enddo
mycall0=mycall
hiscall0=hiscall
endif
! Fill the detmet, detferr arrays
nstepsize=60 ! 5ms steps
nstep=(n-NPTS)/nstepsize
detmet=0
detmet2=0
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)
ahavp=(sum(tonespec,ismask)-ah)/count(ismask)
trath=ah/(ahavp+0.01)
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)
alavp=(sum(tonespec,ismask)-al)/count(ismask)
tratl=al/(alavp+0.01)
fdiff=(ihpk+deltah-ilpk-deltal)*df
i2000=nint(2000/df)+1
i4000=nint(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)
detmet2(istp)=max(trath,tratl)
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),detmet2(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
if( ndet .lt. 3 ) then
do ip=1,MAXCAND-ndet ! Find candidates
iloc=maxloc(detmet2(1:nstep))
il=iloc(1)
if( (detmet2(il) .lt. 20.0) ) exit
if( abs(detfer(il)) .le. ntol ) then
ndet=ndet+1
times(ndet)=((il-1)*nstepsize+NSPM/2)*dt
ferrs(ndet)=detfer(il)
snrs(ndet)=12.0*log10(detmet2(il))/2-9.0
endif
detmet2(max(1,il-1):min(nstep,il+1))=0.0
! detmet2(il)=0.0
enddo
endif
! do ip=1,ndet
! write(*,'(i5,f7.2,f7.2,f7.2)') ip,times(ip),snrs(ip),ferrs(ip)
! enddo
nmessages=0
lines=char(0)
ncalls=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)
xsnr=snrs(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(15,*) i,abs(ccr(i)),ddr(i),abs(cdat(i))
!enddo
! Find 6 largest peaks
do ipk=1,6
iloc=maxloc(ddr)
ic1=iloc(1)
ipeaks(ipk)=ic1
ddr(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,4
! 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)
!write(*,*) abs(bb),bbp
if( ibb .le. 3 ) ibb=ibb-1
if( ibb .gt. 3 ) ibb=ibb-7
do id=1,3 ! 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(cbr))
if( ic+40*6+41 .le. NPTS ) then
ccb=sum(cdat(ic+40*6:ic+40*6+41)*conjg(cbr))
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(cbr))
cfac=cca*conjg(ccb)
ferr2=atan2(imag(cfac),real(cfac))/(twopi*40*6*dt)
endif
do idf=0,2 ! frequency jitter
if( idf .eq. 0 ) then
deltaf=0.0
elseif( mod(idf,2) .eq. 0 ) then
deltaf=2.5*idf
else
deltaf=-2.5*(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(cbr))
phase0=atan2(imag(cca),real(cca))
do ipha=1,3
if( ipha.eq.2 ) phase0=phase0-30*pi/180.0
if( ipha.eq.3 ) phase0=phase0+30*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
! nerr=0
! do i=1,32
! if( testcw(i) .ne. hardbits(i+8) ) nerr=nerr+1
! enddo
! 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
if(xsnr.lt.0.0) sigma=0.75-0.0875*xsnr
lratio(1:32)=exp(2.0*softbits(9:40)/(sigma*sigma)) ! Use this for Radford Neal's routines
llr(1:32)=2.0*softbits(9:40)/(sigma*sigma) ! Use log likelihood for bpdecode40
max_iterations=5
max_dither=1
call bpdecode40(llr,max_iterations, decoded, niterations)
ncalls=ncalls+1
nhashflag=0
if( niterations .ge. 0 ) then
call encode_msk40(decoded,cw)
! call ldpc_encode(decoded,cw)
nhammd=0
cord=0.0
do i=1,32
if( cw(i) .ne. hardbits(i+8) ) then
nhammd=nhammd+1
cord=cord+abs(softbits(i+8))
endif
enddo
imsg=0
do i=1,16
imsg=ishft(imsg,1)+iand(1,decoded(17-i))
enddo
nrxrpt=iand(imsg,15)
nrxhash=(imsg-nrxrpt)/16
! if( nhammd .le. 5 .and. cord .lt. 1.7 .and. nrxhash .eq. nhashes(nrxrpt) ) then
if( nhammd .le. 5 .and. cord .lt. 1.7 .and. nrxhash .eq. ihash ) then
fest=1500+ferr+ferr2+deltaf
!write(14,'(i6.6,11i6,f7.1,f7.1)') nutc,ip,ipk,id,idf,iav,ipha,niterations,nbadsync,nrxrpt,ncalls,nhammd,cord,xsnr
nhashflag=1
msgreceived=' '
nmessages=1
write(msgreceived,'(a1,a,1x,a,a1,1x,a4)') "<",trim(mycall), &
trim(hiscall),">",rpt(nrxrpt)
write(lines(nmessages),1020) nutc,nsnr,t0,nint(fest),msgreceived
1020 format(i6.6,i4,f5.1,i5,' & ',a22)
return
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
return
end subroutine detectmsk40

View File

@ -1,4 +1,4 @@
subroutine fast_decode(id2,narg,ntrperiod,bShMsgs,line,mycall_12, &
subroutine fast_decode(id2,narg,ntrperiod,line,mycall_12, &
hiscall_12)
parameter (NMAX=30*12000)
@ -6,7 +6,6 @@ subroutine fast_decode(id2,narg,ntrperiod,bShMsgs,line,mycall_12, &
integer*2 id2a(NMAX)
integer*2 id2b(NMAX)
integer narg(0:14)
logical*1 bShMsgs
real dat(30*12000)
complex cdat(262145),cdat2(262145)
real psavg(450)
@ -48,42 +47,6 @@ subroutine fast_decode(id2,narg,ntrperiod,bShMsgs,line,mycall_12, &
if(nmode.eq.102) then
call fast9(id2,narg,line)
go to 900
else if(nmode.eq.104) then
! MSK144 mode
if(newdat.eq.1) then
id2b=id2a !Data for lower panel
id2a=id2 !Data for upper panel
nutcb=nutca
nutca=nutc
endif
ia=max(1,nint(t0*12000.0))
ib=nint(t1*12000.0)
if(ib.gt.ntrperiod*12000) ib=ntrperiod*12000
nz=ib-ia+1
if(newdat.eq.1) then
! Full sequence of new data
call msk144_decode(id2a(ia),nz,nutca,0,mycall,hiscall, &
bShMsgs,ntol,t0,line)
go to 100
endif
if(npick.eq.1) then
! Pick-decode from upper panel
call msk144_decode(id2(ia),nz,nutc,0,mycall,hiscall, &
bShMsgs,ntol,t0,line)
go to 100
endif
if(npick.eq.2) then
! Pick-decode from lower panel
! write(*,3001) newdat,npick,nutca
call msk144_decode(id2b(ia),nz,nutca,0,mycall,hiscall, &
bShMsgs,ntol,t0,line)
endif
100 continue
go to 900
endif
if(newdat.eq.1) then

View File

@ -1,97 +0,0 @@
program msk144d
! Test the msk144 decoder for WSJT-X
use options
use timer_module, only: timer
use timer_impl, only: init_timer
use readwav
character c
character*80 line(100)
character*512 pchk_file
logical :: display_help=.false.
logical*1 bShMsgs
type(wav_header) :: wav
integer*2 id2(30*12000)
character*500 infile
character*12 mycall,hiscall
character(len=500) optarg
type (option) :: long_options(6) = [ &
option ('dxcall',.true.,'d','hiscall',''), &
option ('evemode',.true.,'e','',''), &
option ('help',.false.,'h','Display this help message',''), &
option ('mycall',.true.,'m','mycall',''), &
option ('nftol',.true.,'n','nftol',''), &
option ('short',.false.,'s','enable Sh','') &
]
t0=0.0
ntol=100
mycall=''
hiscall=''
bShMsgs=.false.
do
call getopt('d:ehm:n:s',long_options,c,optarg,narglen,nstat,noffset,nremain,.true.)
if( nstat .ne. 0 ) then
exit
end if
select case (c)
case ('d')
read (optarg(:narglen), *) hiscall
case ('e')
t0=1e-4
case ('h')
display_help = .true.
case ('m')
read (optarg(:narglen), *) mycall
case ('n')
read (optarg(:narglen), *) ntol
case ('s')
bShMsgs=.true.
end select
end do
if(display_help .or. nstat.lt.0 .or. nremain.lt.1) then
print *, ''
print *, 'Usage: msk144d [OPTIONS] file1 [file2 ...]'
print *, ''
print *, ' msk144 decode pre-recorded .WAV file(s)'
print *, ''
print *, 'OPTIONS:'
do i = 1, size (long_options)
call long_options(i) % print (6)
end do
go to 999
endif
call init_timer ('timer.out')
call timer('msk144 ',0)
pchk_file='./peg-128-80-reg3.pchk'
ndecoded=0
do ifile=noffset+1,noffset+nremain
call get_command_argument(ifile,optarg,narglen)
infile=optarg(:narglen)
call timer('read ',0)
call wav%read (infile)
i1=index(infile,'.wav')
if( i1 .eq. 0 ) i1=index(infile,'.WAV')
read(infile(i1-6:i1-1),*,err=998) nutc
inquire(FILE=infile,SIZE=isize)
npts=min((isize-216)/2,360000)
read(unit=wav%lun) id2(1:npts)
close(unit=wav%lun)
call timer('read ',1)
call msk144_decode(id2,npts,nutc,1,mycall,hiscall,bShMsgs,ntol,t0,line)
enddo
call timer('msk144 ',1)
call timer('msk144 ',101)
go to 999
998 print*,'Cannot read from file:'
print*,infile
999 continue
end program msk144d

View File

@ -1,225 +0,0 @@
subroutine opdetmsk144(cbig,n,lines,nmessages,nutc,ntol,t00)
use timer_module, only: timer
! NSPM: number of samples per message frame
! NAVG: number of frames to coherently average
parameter (NSPM=864, NAVG=7, NPTS=NAVG*NSPM, NSTEP=6000)
character*22 msgreceived,allmessages(20)
character*80 lines(100)
complex cbig(n)
complex cdat(NPTS) !Analytic signal
complex cdat2(NPTS)
complex c(NSPM)
complex cr(NSPM,NAVG)
complex ct(NSPM)
complex cs(NSPM)
complex cb(42) !Complex waveform for sync word
complex cfac,cca,ccb
complex cc(0:NSPM-1)
complex csum
integer s8(8),hardbits(144)
integer, dimension(1) :: iloc
integer*1 decoded(80)
integer ipeaks(10)
real cbi(42),cbq(42)
real rcw(12)
real ccm(0:NSPM-1)
real ccms(0:NSPM-1)
real pp(12) !Half-sine pulse shape
real*8 dt, fs, pi, twopi
real softbits(144)
real*8 lratio(128)
real llr(128)
logical first
data first/.true./
data s8/0,1,1,1,0,0,1,0/
save first,cb,fs,pi,twopi,dt,s8,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
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)
first=.false.
endif
nmessages=0
allmessages=char(0)
lines=char(0)
nshort=0
trec=NPTS/12000.0 ! Duration of the data record
! write(*,*) "number of points in opdetmsk144",n
if( n .lt. NPTS .or. n .gt. 181000) return
nsteps=2*n/NSTEP-1 ! Hardwired for 0.5s steps
! write(*,*) 'nsteps ',nsteps
nsnr=-4
do istep=1,nsteps
ib=(istep-1)*NPTS/2+1
ie=ib+NPTS-1
if( ie .gt. n ) then
ie=n
ib=n-NPTS+1
endif
t0=t00+(ib-1)/12000.0+trec/2
cdat=cbig(ib:ie)
xmax=0.0
bestf=0.0
do ifr=-ntol,ntol ! search for frequency that maximizes sync correlation
ferr=ifr
! call timer('opdetavg ',0)
call tweak1(cdat,NPTS,-(1500+ferr),cdat2)
cr=reshape(cdat2,(/NSPM,NAVG/))
c=sum(cr,2)
! call timer('opdetavg ',1)
! call timer('opdetsync',0)
cc=0
do ish=0,NSPM-1
ct=cshift(c,ish)
! cc(ish)=sum(ct(1:42)*conjg(cb))+sum(ct(56*6:56*6+41)*conjg(cb))
! cc(ish)=sum(ct(1:42)*conjg(cb)+ ct(56*6:56*6+41)*conjg(cb))
csum=0
do j=1,42
csum=csum+(ct(j)+ct(56*6+j-1))*conjg(cb(j))
enddo
cc(ish)=csum
enddo
ccm=abs(cc)
xb=maxval(ccm)
! call timer('opdetsync',1)
if( xb .gt. xmax ) then
xmax=xb
bestf=ferr
cs=c
ccms=ccm
endif
enddo
fest=1500+bestf
! write(*,*) istep,fest,xmax
c=cs
ccm=ccms
! Find 2 largest peaks
do ipk=1, 2
iloc=maxloc(ccm)
ic2=iloc(1)
ipeaks(ipk)=ic2
ccm(max(0,ic2-7):min(NSPM-1,ic2+7))=0.0
enddo
do ipk=1,2
do is=1,3
ic0=ipeaks(ipk)
if( is.eq.2 ) ic0=max(1,ic0-1)
if( is.eq.3 ) ic0=min(NSPM,ic0+1)
ct=cshift(c,ic0-1)
! Estimate final frequency error and carrier phase.
cca=sum(ct(1:1+41)*conjg(cb))
ccb=sum(ct(1+56*6:1+56*6+41)*conjg(cb))
cfac=ccb*conjg(cca)
ffin=atan2(imag(cfac),real(cfac))/(twopi*56*6*dt)
phase0=atan2(imag(cca+ccb),real(cca+ccb))
! Remove phase error - want constellation rotated so that sample points lie on I/Q axes
cfac=cmplx(cos(phase0),sin(phase0))
ct=ct*conjg(cfac)
if( nmatchedfilter .eq. 0 ) then
! sample to get softsamples
do i=1,72
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(ct(1:6))*pp(7:12))+sum(imag(ct(864-5:864))*pp(1:6))
softbits(2)=sum(real(ct(1:12))*pp)
do i=2,72
softbits(2*i-1)=sum(imag(ct(1+(i-1)*12-6:1+(i-1)*12+5))*pp)
softbits(2*i)=sum(real(ct(7+(i-1)*12-6:7+(i-1)*12+5))*pp)
enddo
endif
! sync word hard error weight is a good discriminator for
! frames that have reasonable probability of decoding
hardbits=0
do i=1,144
if( softbits(i) .ge. 0.0 ) then
hardbits(i)=1
endif
enddo
nbadsync1=(8-sum( (2*hardbits(1:8)-1)*s8 ) )/2
nbadsync2=(8-sum( (2*hardbits(1+56:8+56)-1)*s8 ) )/2
nbadsync=nbadsync1+nbadsync2
!write(*,*) 'peak index ',ipk,' pk loc: ',ic0,' nbadsync ',nbadsync
if( nbadsync .gt. 4 ) cycle
! normalize the softsymbols before submitting to decoder
sav=sum(softbits)/144
s2av=sum(softbits*softbits)/144
ssig=sqrt(s2av-sav*sav)
softbits=softbits/ssig
sigma=0.75
lratio(1:48)=softbits(9:9+47)
lratio(49:128)=softbits(65:65+80-1)
llr=2.0*lratio/(sigma*sigma)
lratio=exp(2.0*lratio/(sigma*sigma))
max_iterations=10
max_dither=1
call timer('bpdec144 ',0)
call bpdecode144(llr,max_iterations,decoded,niterations)
call timer('bpdec144 ',1)
if( niterations .ge. 0.0 ) then
call extractmessage144(decoded,msgreceived,nhashflag)
if( nhashflag .gt. 0 ) then ! CRCs match, so print it
ndupe=0
do im=1,nmessages
if( allmessages(im) .eq. msgreceived ) ndupe=1
enddo
if( ndupe .eq. 0 ) then
nmessages=nmessages+1
allmessages(nmessages)=msgreceived
write(lines(nmessages),1020) nutc,nsnr,t0,nint(fest),msgreceived
1020 format(i6.6,i4,f5.1,i5,' ^ ',a22)
endif
goto 999
else
msgreceived=' '
ndither=-99 ! -99 is bad hash flag
endif
endif
enddo ! slicer dither
enddo ! peak loop
enddo
msgreceived=' '
ndither=-98
999 continue
return
end subroutine opdetmsk144

View File

@ -108,7 +108,7 @@ extern "C" {
float* level, float* sigdb, float* snr, float* dfreq,
float* width);
void fast_decode_(short id2[], int narg[], int* ntrperiod, bool* bShMsgs,
void fast_decode_(short id2[], int narg[], int* ntrperiod,
char msg[], char mycall[], char hiscall[],
int len1, int len2, int len3);
void degrade_snr_(short d2[], int* n, float* db, float* bandwidth);
@ -2383,7 +2383,7 @@ void MainWindow::decode() //decode()
narg[14]=m_config.aggressive();
memcpy(d2b,dec_data.d2,2*360000);
watcher3.setFuture (QtConcurrent::run (std::bind (fast_decode_,&d2b[0],
&narg[0],&m_TRperiod,&m_bShMsgs,&m_msg[0][0],
&narg[0],&m_TRperiod,&m_msg[0][0],
dec_data.params.mycall,dec_data.params.hiscall,8000,12,12)));
} else {
memcpy(to, from, qMin(mem_jt9->size(), size));