Functioning JTMSK decoder, within stand-alone program msk.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/jtms3@2519 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Joe Taylor 2012-07-17 19:31:01 +00:00
parent 71ddb3fd1e
commit 8b55eb3a28
17 changed files with 666 additions and 12 deletions

View File

@ -1,6 +1,6 @@
[Setup]
AppName=JTMS3
AppVerName=JTMS3 Version 0.2 r2516
AppVerName=JTMS3 Version 0.2 r2517
AppCopyright=Copyright (C) 2001-2012 by Joe Taylor, K1JT
DefaultDirName=c:\JTMS3
DefaultGroupName=JTMS3

View File

@ -6,10 +6,10 @@
QT += core gui network
CONFIG += qwt thread
#CONFIG += console
CONFIG += console
TARGET = jtms3
VERSION = 0.1
VERSION = 0.2
TEMPLATE = app
win32 {

View File

@ -17,7 +17,7 @@ CFLAGS = -I. -fbounds-check
%.o: %.F90
${FC} ${FFLAGS} -c $<
all: libm65.a ms3.exe t1.exe
all: libm65.a msk.exe
OBJS1 = trimlist.o display.o getdphi.o pctile.o ccf65.o \
decode1a.o sort.o filbig.o fil6521.o afc65b.o \
@ -32,15 +32,17 @@ OBJS1 = trimlist.o display.o getdphi.o pctile.o ccf65.o \
geocentric.o moon2.o toxyz.o dot.o dcoord.o f77_wisdom.o \
gen65.o chkmsg.o ptt.o astrosub.o astro0.o recvpkt.o symspec.o \
iqcal.o iqfix.o timf2.o s3avg.o genjtms3.o analytic.o \
db.o specjtms.o genmsk.o
db.o specjtms.o genmsk.o mskdf.o tweak1.o syncmsk.o \
lenmsk.o decodemsk.o ping.o makepings.o alignmsg.o match.o \
rtping.o jtmsk.o hipass.o setupmsk.o foldmsk.o
libm65.a: $(OBJS1)
ar cr libm65.a $(OBJS1)
ranlib libm65.a
OBJS3 = ms3.o
ms3.exe: $(OBJS3) libm65.a
$(FC) -o ms3.exe $(OBJS3) libm65.a ../libfftw3f_win.a \
OBJS3 = msk.o
msk.exe: $(OBJS3) libm65.a
$(FC) -o msk.exe $(OBJS3) libm65.a ../libfftw3f_win.a
OBJS2 = JT65code.o
JT65code.exe: $(OBJS2) libm65.a

View File

@ -2,9 +2,10 @@ subroutine analytic(d,npts,nfft,s,c)
! Convert real data to analytic signal
parameter (NFFTMAX=128*1024)
real d(npts)
real s(npts)
complex c(npts)
complex c(NFFTMAX)
nh=nfft/2
fac=2.0/nfft

39
libm65/decodemsk.f90 Normal file
View File

@ -0,0 +1,39 @@
subroutine decodemsk(cdat,npts,cw,i1,nchar,s2,msg)
! DF snd sync have been established, now decode the message
complex cdat(npts)
complex cw(168,0:63) !Complex waveforms for codewords
real s2(0:63,400)
character msg*400
complex z,zmax
character cc*64
! 1 2 3 4 5 6
! 0123456789012345678901234567890123456789012345678901234567890123
data cc/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ./?- _ @'/
msg=' '
do j=1,nchar !Find best match for each character
ia=i1 + (j-1)*168
smax=0.
do k=0,40
kk=k
if(k.eq.40) kk=57
z=0.
do i=1,168
z=z + cdat(ia+i)*conjg(cw(i,kk))
enddo
ss=abs(z)
s2(k,j)=ss
if(ss.gt.smax) then
smax=ss
zmax=z
kpk=kk
endif
enddo
msg(j:j)=cc(kpk+1:kpk+1)
if(kpk.eq.57) msg(j:j)=' '
enddo
return
end subroutine decodemsk

50
libm65/foldmsk.f90 Normal file
View File

@ -0,0 +1,50 @@
subroutine foldmsk(s2,msglen,nchar,mycall,msg,msg29)
! Fold the 2-d "goodness of fit" array s2 modulo message length,
! then decode the folded message.
real s2(0:63,400)
real fs2(0:63,29)
integer nfs2(29)
character mycall*12
character msg*400,msg29*29
character cc*64
! 1 2 3 4 5 6
! 0123456789012345678901234567890123456789012345678901234567890123
data cc/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ./?- _ @'/
fs2=0.
nfs2=0
do j=1,nchar !Fold s2 into fs2, modulo msglen
jj=mod(j-1,msglen)+1
nfs2(jj)=nfs2(jj)+1
do i=0,40
fs2(i,jj)=fs2(i,jj) + s2(i,j)
enddo
enddo
msg=' '
do j=1,msglen
smax=0.
do k=0,40
if(fs2(k,j).gt.smax) then
smax=fs2(k,j)
kpk=k
endif
enddo
if(kpk.eq.40) kpk=57
msg(j:j)=cc(kpk+1:kpk+1)
if(kpk.eq.57) msg(j:j)=' '
enddo
msg29=msg(1:msglen)
call alignmsg(' ',2,msg29,msglen,idone)
if(idone.eq.0) call alignmsg('CQ', 3,msg29,msglen,idone)
if(idone.eq.0) call alignmsg('QRZ', 3,msg29,msglen,idone)
if(idone.eq.0) call alignmsg(mycall,4,msg29,msglen,idone)
if(idone.eq.0) call alignmsg(' ', 1,msg29,msglen,idone)
msg29=adjustl(msg29)
return
end subroutine foldmsk

View File

@ -65,7 +65,6 @@ subroutine genmsk(msg28,iwave,nwave)
f0=48000.d0/nsps
dfgen=0.5d0*f0
foffset=1500.d0 - f0
print*,f0,dfgen,foffset
t=0.d0
k=0
phi=0.d0
@ -95,5 +94,9 @@ subroutine genmsk(msg28,iwave,nwave)
if(isrch.eq.0) iwave(k+1:)=0
nwave=k
! call makepings(iwave,nwave)
! write(71) iwave
! call flush(71)
return
end subroutine genmsk

86
libm65/jtmsk.f90 Normal file
View File

@ -0,0 +1,86 @@
subroutine jtmsk(dat,npts,cfile6,t2,mswidth,ndb,nrpt,Nfreeze, &
ntol,MouseDF,pick)
! Decode a JTMS ping
parameter (NZ=30*48000)
real dat(npts) !Raw data
complex cdat(NZ) !Analytic form of signal
character*6 cfile6 !FileID
logical pick
character*12 mycall,hiscall
real s(NZ) !Power spectrum
real s2(0:63,400)
real r(60000)
complex cw(168,0:63) !Complex waveforms for all codewords
complex cwb(168) !Complex waveform for <space>
logical first
character msg*400,msg29*29
character*90 line
common/ccom/nline,tping(100),line(100)
data first/.true./
save first,cw,cwb
save cdat !Fix its address, for four2
if(first) call setupmsk(cw,cwb) !Calculate waveforms for codewords
first=.false.
nsps=24 !Samples per symbol
f0=1000.0 !Nominal frequency for bit=0
n=log(float(npts))/log(2.0) + 1.0
nfft1=2**n !FFT length
call analytic(dat,npts,nfft1,s,cdat) !Convert to analytic signal
call mskdf(cdat,npts,t2,nfft1,f0,nfreeze,mousedf,ntol, &
dfx,snrsq2) !Get DF
print*,'b',dfx,snrsq2
sq2lim=7.0
if(pick) sq2lim=5.0
if(snrsq2.lt.sq2lim) go to 900 !Reject non-JTMS signals
call tweak1(cdat,npts,-dfx,cdat) !Mix to standard frequency
! DF is known, now establish character sync.
call syncmsk(cdat,npts,cwb,r,i1) !Get character sync
call lenmsk(r,npts,msglen) !Find message length
s2=0.
nchar=(npts-168+1-i1)/168
if(nchar.gt.400) nchar=400
call decodemsk(cdat,npts,cw,i1,nchar,s2,msg) !Decode the message
print*,'B',i1,msglen,nchar
print*,msg(1:nchar)
! ia=1
! if(nchar.ge.40) ia=min(nchar/3,nchar-28)
! ib=min(ia+28,nchar) !Can better limits ia, ib be found?
! print*,'A',ia,ib,nchar
! print*,msg(1:nchar)
! msg29=adjustl(msg(ia:ib))
msg=adjustl(msg)
ib=min(nchar,45)
ndf=nint(dfx)
nchk=max(20,nint(1.5*msglen))
if(msglen.eq.0 .or. nchar.lt.nchk .or. pick) then
if(nline.le.99) nline=nline+1
tping(nline)=t2
write(line(nline),1110) cfile6,t2,mswidth,ndb,nrpt,ndf,msg(1:45)
1110 format(a6,f5.1,i5,i3,1x,i2.2,i5,5x,a45)
endif
if(msglen.gt.0 .and. nchar.ge.nchk) then
call foldmsk(s2,msglen,nchar,mycall,msg,msg29) !Decode folded message
if(nline.le.99) nline=nline+1
tping(nline)=t2
write(line(nline),1120) cfile6,t2,mswidth,ndb,nrpt,ndf,msg29
1120 format(a6,f5.1,i5,i3,1x,i2.2,i5,5x,a29,11x,'*')
endif
900 continue
return
end subroutine jtmsk

56
libm65/lenmsk.f90 Normal file
View File

@ -0,0 +1,56 @@
subroutine lenmsk(r,npts,msglen)
! Determine length of the user message in a JTMS ping.
real r(60000)
real acf(4872)
integer np(9)
data np/5,7,9,11,13,17,19,23,29/ !Permissible message lengths
save acf !Why necessary? (But don't remove!)
msglen=0 !Use ACF to find msg length
if(npts.ge.8*168) then
r=r-sum(r(1:npts))/npts
acfmax=0.
acf0=dot_product(r(1:npts),r(1:npts))
kz=min(nint(0.75*npts),29*168)
do k=8,kz
fac=float(npts)/(npts-k)
acf(k)=fac*dot_product(r(1:npts),r(1+k:npts+k))/acf0
enddo
call hipass(acf(8),kz-7,50)
do k=8,kz !Find acfmax, kpk
if(acf(k).gt.acfmax) then
acfmax=acf(k)
kpk=k
endif
enddo
sumsq=0.
n=0
do k=8,kz !Find rms, skipping around kpk
if(abs(k-kpk).gt.10) then
sumsq=sumsq+acf(k)**2
n=n+1
endif
enddo
rms=sqrt(sumsq/n)
acf=acf/rms !Normalize the acf
amax2=0.
acflim=3.5
do i=1,9
k=168*np(i) !Check only the permitted lengths
if(k.gt.kz) go to 10
if(acf(k).gt.acflim .and. acf(k).gt.amax2) then
amax2=acf(k) !Save best value >3.5 sigma
msglen=np(i) !Save message length
kpk2=k
endif
enddo
10 continue
endif
return
end subroutine lenmsk

59
libm65/msk.f90 Normal file
View File

@ -0,0 +1,59 @@
program msk
! Starting code for a JTMSK decoder.
parameter (NSMAX=30*48000)
character*80 infile
character*6 cfile6
character*12 arg
real dat(NSMAX)
real x(NSMAX)
complex cx(0:NSMAX/2)
integer hdr(11)
integer*2 id
common/mscom/id(NSMAX),s1(215,703),s2(215,703)
nargs=iargc()
if(nargs.lt.1) then
print*,'Usage: msk <snr>'
go to 999
endif
call getarg(1,arg)
read(arg,*) snr
open(71,file='dat.71',form='unformatted',status='old')
read(71) id
cfile6='123400'
npts=30*48000
kstep=2048
minsigdb=6
mousedf=0
ntol=200
call random_number(x)
nfft=NSMAX
call four2a(x,nfft,1,-1,0)
df=48000.0/nfft
ia=nint(300.0/df)
ib=nint(2800.0/df)
cx(:ia)=0.
cx(ib:)=0.
call four2a(cx,nfft,1,1,-1)
x(1)=0.
sq=0.
do i=1,NSMAX
sq=sq + x(i)**2
enddo
rms=sqrt(sq/NSMAX)
x=x/rms
sig=(10.0**(0.05*snr))/32768.0
dat=sig*id + x
k=0
do iblk=1,npts/kstep
k=k+kstep
call rtping(dat,k,cfile6,MinSigdB,MouseDF,ntol)
enddo
999 end program msk

49
libm65/msk0.f90 Normal file
View File

@ -0,0 +1,49 @@
program msk
! Starting code for a JTMSK decoder.
parameter (NSMAX=30*48000)
character*80 infile
character*6 cfile6
real dat(NSMAX)
integer hdr(11)
integer*2 id
common/mscom/id(NSMAX),s1(215,703),s2(215,703)
nargs=iargc()
if(nargs.lt.1) then
print*,'Usage: msk file1 [file2 ...]'
print*,' Reads data from *.wav files.'
go to 999
endif
npts=30*48000
kstep=2048
minsigdb=6
mousedf=0
ntol=200
do ifile=1,nargs
call getarg(ifile,infile)
open(10,file=infile,access='stream',status='old',err=998)
read(10) hdr
read(10) id
close(10)
hdr(1)=hdr(2)
i1=index(infile,'.wav')
cfile6=infile(i1-6:i1-1)
dat=id
k=0
do iblk=1,npts/kstep
k=k+kstep
call rtping(dat,k,cfile6,MinSigdB,MouseDF,ntol)
enddo
enddo
go to 999
998 print*,'Cannot open file:'
print*,infile
999 end program msk

75
libm65/mskdf.f90 Normal file
View File

@ -0,0 +1,75 @@
subroutine mskdf(cdat,npts,t2,nfft1,f0,nfreeze,mousedf,ntol,dfx,snrsq2)
! Determine DF for a JTMS signal. Also find ferr, a measure of
! frequency differerence between 1st and 2nd harmonic.
! (Should be 0.000)
parameter (NZ=128*1024)
complex cdat(npts)
integer ntol
real sq(NZ)
real ccf(-2600:2600) !Correct limits?
real tmp(NZ)
complex c(NZ)
data nsps/8/
save c
df1=48000.0/nfft1
nh=nfft1/2
fac=1.0/(nfft1**2)
do i=1,npts
c(i)=fac*cdat(i)**2
enddo
c(npts+1:nfft1)=0.
call four2a(c,nfft1,1,-1,1)
! In the "doubled-frequencies" spectrum of squared cdat:
fa=2.0*(f0-400)
fb=2.0*(f0+400)
j0=nint(2.0*f0/df1)
ja=nint(fa/df1)
jb=nint(fb/df1)
jd=nfft1/nsps
do j=1,nh+1
sq(j)=real(c(j))**2 + aimag(c(j))**2
! if(j*df1.lt.6000.0) write(54,3009) j*df1,sq(j),db(sq(j))
!3009 format(3f12.3)
enddo
ccf=0.
do j=ja,jb
ccf(j-j0-1)=sq(j)+sq(j+jd)
enddo
call pctile(ccf(ja-j0-1),tmp,jb-ja+1,50,base)
ccf=ccf/base
if(NFreeze.gt.0) then
fa=2.0*(f0+MouseDF-ntol)
fb=2.0*(f0+MouseDF+ntol)
endif
ja=nint(fa/df1)
jb=nint(fb/df1)
! rewind 51
smax=0.
do j=ja,jb
k=j-j0-1
if(ccf(k).gt.smax) then
smax=ccf(k)
jpk=j
endif
f=0.5*k*df1
! write(51,3002) f,ccf(k)
!3002 format(2f12.3)
enddo
! call flush(51)
fpk=(jpk-1)*df1
dfx=0.5*fpk-f0
snrsq2=smax
return
end subroutine mskdf

42
libm65/ping.f90 Normal file
View File

@ -0,0 +1,42 @@
subroutine ping(s,nz,dtbuf,slim,wmin,pingdat,nping)
! Detect pings and make note of their start time, duration, and peak strength.
real*4 s(nz)
real*4 pingdat(3,100)
logical inside
nping=0
peak=0.
inside=.false.
snrlim=10.0**(0.1*slim) - 1.0
sdown=10.0*log10(0.25*snrlim+1.0)
i0=0 !Silence compiler warnings.
tstart=0.0
do i=2,nz
if(s(i).ge.slim .and. .not.inside) then
i0=i
tstart=i0*dtbuf
inside=.true.
peak=0.
endif
if(inside .and. s(i).gt.peak) then
peak=s(i)
endif
if(inside .and. (s(i).lt.sdown .or. i.eq.nz)) then
if(i.gt.i0) then
if(dtbuf*(i-i0).ge.wmin) then
if(nping.le.99) nping=nping+1
pingdat(1,nping)=tstart
pingdat(2,nping)=dtbuf*(i-i0)
pingdat(3,nping)=peak
endif
inside=.false.
peak=0.
endif
endif
enddo
return
end subroutine ping

103
libm65/rtping.f90 Normal file
View File

@ -0,0 +1,103 @@
subroutine rtping(dat,k,cfile6,MinSigdB,MouseDF,ntol)
!subroutine rtping(dat,jz,nz,MinSigdB,MinWidth,NFreeze,DFTolerance, &
! MouseDF,istart,pick,cfile6,mycall,hiscall,mode,ps0)
! Decode Multi-Tone FSK441 mesages.
parameter (NSMAX=30*48000)
parameter (NZMAX=NSMAX/2048)
real dat(NSMAX) !Raw audio data
logical pick
character*6 cfile6
real sig(NZMAX) !Sq-law detected signal, sampled at 43 ms
real sigdb(NZMAX) !Signal in dB, sampled at 43 ms
real work(NZMAX)
real pingdat(3,100)
! character msg*40,msg3*3
character*90 line
common/ccom/nline,tping(100),line(100)
data nping0/0/
save
slim=MinSigdB
! nf1=-ntol
! nf2=ntol
dt=1.0/48000.0
kstep=2048
! pick=.false.
istart=1
jz=k
! Find signal power
j=k/kstep
sig(j)=dot_product(dat(k-kstep+1:k),dat(k-kstep+1:k))/kstep
if(j.lt.10) return
! Remove baseline, compute signal level in dB
call pctile (sig,work,j,50,base1)
do i=1,j
sigdb(i)=db(sig(i)/base1)
if(j.eq.703) write(13,3001) i,sig(i),sigdb(i)
3001 format(i5,2e12.3)
enddo
dtbuf=kstep*dt
wmin=0.040
call ping(sigdb,j,dtbuf,slim,wmin,pingdat,nping)
! If this is a "mouse pick" and no ping was found, force a pseudo-ping
! at center of data.
! if(pick.and.nping.eq.0) then
! if(nping.le.99) nping=nping+1
! pingdat(1,nping)=0.5*jz*dt
! pingdat(2,nping)=0.16
! pingdat(3,nping)=1.0
! endif
do iping=1,nping
! Find starting place and length of data to be analyzed:
tstart=pingdat(1,iping)
width=pingdat(2,iping)
peak=pingdat(3,iping)
! mswidth=10*nint(100.0*width)
jj=(tstart-0.02)/dt
if(jj.lt.1) jj=1
jjz=nint((width+0.02)/dt)+1
jjz=min(jjz,jz+1-jj)
! Compute average spectrum of this ping.
! call spec441(dat(jj),jjz,ps,f0)
! Decode the message.
! msg=' '
! call longx(dat(jj),jjz,ps,DFTolerance,noffset,msg,msglen,bauderr)
! Assemble a signal report:
nwidth=0
if(width.ge.0.04) nwidth=1 !These might depend on NSPD
if(width.ge.0.12) nwidth=2
if(width.gt.1.00) nwidth=3
nstrength=6
if(peak.ge.11.0) nstrength=7
if(peak.ge.17.0) nstrength=8
if(peak.ge.23.0) nstrength=9
! nrpt=10*nwidth + nstrength
t2=tstart + dt*(istart-1)
jjzz=min(jjz,2*48000) !Max data size 2 s
!###
jjzz=14400
jj=jj-200
!###
if(nping.gt.nping0) then
print*,'a',jj,jjzz,jj*dt,jjzz*dt,t2,width
call jtmsk(dat(jj),jjzz,cfile6,t2,mswidth,int(peak),nrpt, &
nfreeze,DFTolerance,MouseDF,pick)
nping0=nping
endif
enddo
return
end subroutine rtping

51
libm65/setupmsk.f90 Normal file
View File

@ -0,0 +1,51 @@
subroutine setupmsk(cw,cwb)
! Calculate the JTMS character waveforms.
complex cw(168,0:63)
complex cwb(168)
integer nb(7)
! real*8 twopi,dt,f0,f1
character cc*64
! 1 2 3 4 5 6
! 0123456789012345678901234567890123456789012345678901234567890123
data cc/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ./?- _ @'/
nsps=24
twopi=8.d0*atan(1.d0)
dt=1.d0/48000.d0 !Sample interval
f0=1000.d0
f1=2000.d0
dphi0=twopi*dt*f0
dphi1=twopi*dt*f1
do i=0,63
k=0
m=0
do n=5,0,-1 !Each character gets 6+1 bits
k=k+1
nb(k)=iand(1,ishft(i,-n))
m=m+nb(k)
enddo
k=k+1
nb(k)=iand(m,1) !Insert parity bit
phi=0.
j=0
do k=1,7 !Generate the waveform
if(nb(k).eq.0) then
dphi=dphi0
else
dphi=dphi1
endif
do ii=1,nsps
j=j+1
phi=phi+dphi
cw(j,i)=cmplx(cos(phi),sin(phi))
enddo
enddo
enddo
cwb=cw(1:168,57)
return
end subroutine setupmsk

38
libm65/syncmsk.f90 Normal file
View File

@ -0,0 +1,38 @@
subroutine syncmsk(cdat,npts,cwb,r,i1)
! Establish character sync within a JTMS ping.
complex cdat(npts) !Analytic signal
complex cwb(168) !Complex waveform for 'space'
real r(60000)
real tmp(60000)
integer hist(168),hmax(1)
complex z
r=0.
jz=npts-168+1
do j=1,jz
z=0.
ss=0.
do i=1,168
ss=ss + abs(cdat(i+j-1)) !Total power
z=z + cdat(i+j-1)*conjg(cwb(i)) !Signal matching <space>
enddo
r(j)=abs(z)/ss !Goodness-of-fit to <space>
! write(52,3001) j/168.0,r(j),cdat(j)
!3001 format(4f12.3)
enddo
ncut=99.0*float(jz-10)/float(jz)
call pctile(r,tmp,jz,ncut,rlim)
hist=0
do j=1,jz
k=mod(j-1,168)+1
if(r(j).gt.rlim) hist(k)=hist(k)+1
enddo
hmax=maxloc(hist)
i1=hmax(1)
return
end subroutine syncmsk

View File

@ -1309,7 +1309,7 @@ void MainWindow::ba2msg(QByteArray ba, char message[]) //ba2msg()
{
bool eom;
eom=false;
for(int i=0;i<22; i++) {
for(int i=0;i<28; i++) {
if((int)ba[i] == 0) eom=true;
if(eom) {
message[i]=32;
@ -1317,7 +1317,7 @@ void MainWindow::ba2msg(QByteArray ba, char message[]) //ba2msg()
message[i]=ba[i];
}
}
message[22]=0;
message[28]=0;
}
void MainWindow::on_txFirstCheckBox_stateChanged(int nstate) //TxFirst