Temporary, more-or-less working version of symspec.f90.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@2623 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Joe Taylor 2012-10-01 19:05:22 +00:00
parent e36f1be536
commit e72a766337
4 changed files with 184 additions and 174 deletions

46
jt9.txt
View File

@ -1,31 +1,31 @@
JT9 is a mode designed for amateur QSOs and beacon-like transmissions
at MF and LF. The mode uses the same 72-bit user messages as JT65.
Convolutional error-control coding (ECC) uses constraint length K=32,
rate r=1/2, and a zero tail, which leads to an encoded message length
of (72+31)*2 = 206 bits. Modulation is 9-FSK: 8 tones for data, one
for synchronization. Sixteen symbol intervals are used for
synchronization, so a transmission requires a total of 207/3 + 16 = 85
channel symbols. Symbol durations tsym are approximately
(TRperiod-10)/85, where TRperiod is the T/R sequence length in
seconds. Exact symbol lengths are chosen so that nsps, the number of
samples per symbol (at 12000 samples per second) is a number with no
prime factor greater than 7. This choice makes for efficient FFTs.
Tone spacing of the 9-FSK modulation is df=1/tsym=12000/nsps, the same
as the keying rate. The total occupied bandwidth is 9*df.
JT9 is a mode designed for amateur QSOs at MF and LF. The mode uses
the same 72-bit structured messages as JT65. Error control coding
(ECC) uses a convolutional code with constraint length K=32, rate
r=1/2, and a zero tail, leading to an encoded message length of
(72+31)*2 = 206 information-carrying bits. Modulation is 9-FSK: 8
tones for data, one for synchronization. Sixteen symbol intervals are
used for synchronization, so a transmission requires a total of 207/3
+ 16 = 85 channel symbols. Symbol durations tsym are approximately
(TRperiod-8)/85, where TRperiod is the T/R sequence length in seconds.
Exact symbol lengths are chosen so that nsps, the number of samples
per symbol (at 12000 samples per second) is a number with no prime
factor greater than 7. This choice makes for efficient FFTs. Tone
spacing of the 9-FSK modulation is df=1/tsym=12000/nsps, the same as
the keying rate. The total occupied bandwidth is 9*df.
Parameters of five JT9 sub-modes are summarized in the following
table, along with S/N thresholds measured by simulation on an AWGN
(additive white Gaussian noise) channel.
channel.
-------------------------------------------------------------------------
Mode nsps nsps2 df tsym BW S/N* Tdec Tfree Factors
12000 750 (Hz) (s) (Hz) (dB) (s) (s) of nsps
12000 1500 (Hz) (s) (Hz) (dB) (s) (s) of nsps
-------------------------------------------------------------------------
JT9-1 6912 432 1.736 0.58 15.6 -26.9 52.5 7.5 2^8 3^3
JT9-2 15360 960 0.781 1.28 7.0 -30.2 112.3 7.7 2^10 3 5
JT9-5 40960 2560 0.293 3.41 2.6 -34.4 293.6 6.4 2^13 5
JT9-10 82944 5184 0.145 6.91 1.3 -37.5 591.0 9.0 2^10 3^4
JT9-30 250880 15750 0.048 20.91 0.4 -42.3 1788.5 11.5 2^5 3^2 5^3 7
JT9-1 6912 864 1.736 0.58 15.6 -26.9 52.5 7.5 2^8 3^3
JT9-2 15360 1920 0.781 1.28 7.0 -30.2 112.3 7.7 2^10 3 5
JT9-5 40960 5120 0.293 3.41 2.6 -34.4 293.6 6.4 2^13 5
JT9-10 82944 10368 0.145 6.91 1.3 -37.5 591.0 9.0 2^10 3^4
JT9-30 252000 31500 0.048 21.00 0.4 -42.3 1788.5 11.5 2^5 3^2 5^3 7
-------------------------------------------------------------------------
* Noise power measured in a 2500 Hz bandwidth.
@ -43,8 +43,8 @@ Transmitting
Receiving
---------
1. Apply noise blanking with the timf2 method
2. Filter to 500 Hz bandwidth and downsample (1/16) to 750 Hz (FIR
filter with 61 taps). Complex data to array c0 (max 1.35M)
2. Filter to 1000 Hz bandwidth and downsample (1/8) to 1500 Hz, saving
complex data to array c0(1350000).
3. Compute symbol-length spectra at half-symbol steps. Use for
waterfall display s(15750) and save in ss(184,15750) and
savg(15750) for detecting sync vectors.

View File

@ -3,25 +3,21 @@ program jt9
! Decoder for JT9. Can run stand-alone, reading data from *.wav files;
! or as the back end of wsjt-x, with data placed in a shared memory region.
parameter (NSMAX=60*96000)
parameter (NFFT=32768)
integer*2 i2(4,87)
real*8 hsym
real*4 ssz5a(NFFT)
logical*1 lstrong(0:1023)
common/tracer/limtrace,lu
real*8 fc0,fcenter
character*80 arg,infile
character mycall*12,hiscall*12,mygrid*6,hisgrid*6,datetime*20
common/datcom/dd(4,5760000),ss(4,322,NFFT),savg(4,NFFT),fc0,nutc0,junk(34)
common/npar/fcenter,nutc,idphi,mousedf,mousefqso,nagain, &
ndepth,ndiskdat,neme,newdat,nfa,nfb,nfcal,nfshift, &
mcall3,nkeep,ntol,nxant,nrxlog,nfsample,nxpol,mode65, &
mycall,mygrid,hiscall,hisgrid,datetime
parameter (NMAX=1800*12000) !Total sample intervals per 30 minutes
parameter (NDMAX=1800*1500) !Sample intervals at 1500 Hz rate
parameter (NSMAX=22000) !Max length of saved spectra
integer*4 ihdr(11)
real*4 s(NSMAX)
logical*1 lstrong(0:1023)
integer*2 id2
complex c0
common/jt8com/id2(NMAX),ss(184,NSMAX),savg(NSMAX),c0(NDMAX),nutc,junk(20)
common/tracer/limtrace,lu
nargs=iargc()
if(nargs.lt.1) then
print*,'Usage: jt9 TRp file1 [file2 ...]'
print*,'Usage: jt9 TRperiod file1 [file2 ...]'
print*,' Reads data from *.wav files.'
print*,''
print*,' jt9 -s'
@ -30,101 +26,83 @@ program jt9
endif
call getarg(1,arg)
if(arg(1:2).eq.'-s') then
call m65a
call ftnquit
! call jt9a
! call ftnquit
go to 999
endif
nfsample=96000
nxpol=1
mode65=2
ifile1=1
if(arg.eq.'95238') then
nfsample=95238
call getarg(2,arg)
ifile1=2
endif
limtrace=0
read(arg,*) ntrperiod
ifile1=2
limtrace=10000
lu=12
nfa=100
nfb=162
nfshift=6
ndepth=2
nfcal=344
idphi=-50
ntol=500
nkeep=10
call timer('jt9 ',0) !###
call ftninit('.')
nfa=1000
nfb=2000
ntol=500
mousedf=0
mousefqso=1500
newdat=1
nb=0
nbslider=100
! call ftninit('.')
do ifile=ifile1,nargs
call getarg(ifile,infile)
open(10,file=infile,access='stream',status='old',err=998)
i1=index(infile,'.tf2')
read(10) ihdr
i1=index(infile,'.wav')
read(infile(i1-4:i1-1),*,err=1) nutc0
go to 2
1 nutc0=0
2 hsym=2048.d0*96000.d0/11025.d0 !Samples per half symbol
nhsym0=-999
k=0
fcenter=144.125d0
mousedf=0
mousefqso=125
newdat=1
mycall='K1JT'
2 nsps=0
if(ntrperiod.eq.1) nsps=6912
if(ntrperiod.eq.2) nsps=15360
if(ntrperiod.eq.5) nsps=40960
if(ntrperiod.eq.10) nsps=82944
if(ntrperiod.eq.30) nsps=252000
if(nsps.eq.0) stop 'Error: bad TRprtiod'
if(ifile.eq.ifile1) call timer('m65 ',0)
do irec=1,9999999
call timer('read_tf2',0)
read(10) i2
call timer('read_tf2',1)
call timer('float ',0)
do i=1,87
k=k+1
dd(1,k)=i2(1,i)
dd(2,k)=i2(2,i)
dd(3,k)=i2(3,i)
dd(4,k)=i2(4,i)
enddo
call timer('float ',1)
nhsym=(k-2048)/hsym
kstep=nsps/2
k=0
nhsym0=-999
npts=(60*ntrperiod-6)*12000
call timer('read_wav',0)
read(10) id2(1:npts)
call timer('read_wav',1)
! do i=1,npts
! id2(i)=100.0*sin(6.283185307*1046.875*i/12000.0)
! enddo
! if(ifile.eq.ifile1) call timer('jt9 ',0)
do iblk=1,npts/kstep
k=iblk*kstep
nhsym=(k-2048)/kstep
if(nhsym.ge.1 .and. nhsym.ne.nhsym0) then
ndiskdat=1
nb=0
! Emit signal readyForFFT
call timer('symspec ',0)
fgreen=-13.0
iqadjust=1
iqapply=1
nbslider=100
gainx=0.9962
gainy=1.0265
phasex=0.01426
phasey=-0.01195
call symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample,fgreen, &
iqadjust,iqapply,gainx,gainy,phasex,phasey,rejectx,rejecty, &
pxdb,pydb,ssz5a,nkhz,ihsym,nzap,slimit,lstrong)
call symspecx(k,ntrperiod,nsps,ndiskdat,nb,nbslider,pxdb,s, &
ihsym,nzap,slimit,lstrong)
call timer('symspec ',1)
nhsym0=nhsym
if(ihsym.ge.278) go to 10
if(ihsym.ge.184) go to 10
endif
enddo
10 continue
if(iqadjust.ne.0) write(*,3002) rejectx,rejecty
3002 format('Image rejection:',2f7.1,' dB')
nutc=nutc0
nstandalone=1
call decode0(dd,ss,savg,nstandalone,nfsample)
! call decode0(dd,ss,savg,nstandalone,nfsample)
enddo
call timer('m65 ',1)
call timer('m65 ',101)
call ftnquit
call timer('jt9 ',1)
call timer('jt9 ',101)
! call ftnquit
go to 999
998 print*,'Cannot open file:'
print*,infile
999 end program m65
999 end program jt9

View File

@ -1,51 +1,74 @@
subroutine symspecx(k,nsps,ndiskdat,nb,nbslider,pxdb,s,ihsym, &
subroutine symspecx(k,ntrperiod,nsps,ndiskdat,nb,nbslider,pxdb,s,ihsym, &
nzap,slimit,lstrong)
! k pointer to the most recent new data
! nsps samples per symbol (at 12000 Hz)
! ndiskdat 0/1 to indicate if data from disk
! nb 0/1 status of noise blanker (off/on)
! pxdb power (0-60 dB)
! s spectrum for waterfall display
! ihsym index number of this half-symbol (1-322)
! nzap number of samples zero'ed by noise blanker
! Input:
! k pointer to the most recent new data
! ntrperiod T/R sequence length, minutes
! nsps samples per symbol (12000 Hz)
! ndiskdat 0/1 to indicate if data from disk
! nb 0/1 status of noise blanker (off/on)
! nbslider NB setting, 0-100
! Output:
! pxdb power (0-60 dB)
! s spectrum for waterfall display
! ihsym index number of this half-symbol (1-322)
! nzap number of samples zero'ed by noise blanker
! slimit NB scale adjustment
! lstrong true if strong signal at this freq
parameter (NMAX=1800*12000) !Total sample intervals per 30 minutes
parameter (NSMAX=10000) !Max length of saved spectra
parameter (MAXFFT=262144) !Max length of FFTs
integer*2 id2
real*8 ts,hsym
real*8 fcenter
real*4 s(NFFT),w(NFFT)
parameter (NDMAX=1800*1500) !Sample intervals at 1500 Hz rate
parameter (NSMAX=22000) !Max length of saved spectra
parameter (NFFT1=1024)
parameter (NFFT2=1024,NFFT2A=NFFT2/8)
parameter (MAXFFT3=32768)
real*4 s(NSMAX),w(NFFT1),w3(MAXFFT3)
real*4 stmp(NFFT2/2)
real*4 x0(NFFT1),x1(NFFT1)
real*4 x2(NFFT2)
complex cx2(0:NFFT2/2)
complex cx2a(NFFT2A)
complex z,zfac
complex zsumx
complex cx(NFFT)
complex cx00(NFFT)
complex cx(MAXFFT3)
complex cx00(NFFT1)
complex cx0(0:1023),cx1(0:1023)
logical*1 lstrong(0:1023)
common/jt8com/id2(NMAX),ss(184,NSMAX),savg(NSMAX),fcenter,nutc,junk(20)
data rms/999.0/,k0/99999999/,ntrperiod0/0/
integer*2 id2
complex c0
common/jt8com/id2(NMAX),ss(184,NSMAX),savg(NSMAX),c0(NDMAX),nutc,junk(20)
equivalence (x2,cx2)
data rms/999.0/,k0/99999999/,ntrperiod0/0/,nfft3z/0/
save
nfft3=nsps
hsym=nsps/2
if(ntrperiod.eq.1) nfft3=1024
if(ntrperiod.eq.2) nfft3=2048
if(ntrperiod.eq.5) nfft3=6144
if(ntrperiod.eq.10) nfft3=12288
if(ntrperiod.eq.30) nfft3=32768
jstep=nsps/16
if(k.gt.NMAX) go to 999
if(k.lt.nfft3) then
ihsym=0
go to 999 !Wait for enough samples to start
go to 999 !Wait for enough samples to start
endif
if(k0.eq.99999999) then
if(nfft3.ne.nfft3z) then
pi=4.0*atan(1.0)
do i=1,nfft3
w(i)=(sin(i*pi/nfft3))**2 !Window for nfft3
w3(i)=(sin(i*pi/nfft3))**2 !Window for nfft3
enddo
stmp=0.
nfft3z=nfft3
endif
if(k.lt.k0) then
ts=1.d0 - hsym
ja=-2*jstep
savg=0.
ihsym=0
k1=0
k8=0
if(ndiskdat.eq.0) id2(k+1:)=0. !### Should not be needed ??? ###
endif
k0=k
@ -55,64 +78,73 @@ subroutine symspecx(k,nsps,ndiskdat,nb,nbslider,pxdb,s,ihsym, &
peaklimit=sigmas*max(10.0,rms)
faclim=3.0
px=0.
df2=12000.0/NFFT2
nwindow=2
! nwindow=0 !### No windowing ###
nfft1=1024
kstep=nfft1
if(nwindow.ne.0) kstep=nfft1/2
nblks=(k-k1)/kstep
! nwindow=2
nwindow=0 !### No windowing ###
kstep1=NFFT1
if(nwindow.ne.0) kstep1=NFFT1/2
fac=1.0/(NFFT1*NFFT2)
nblks=(k-k1)/kstep1
do nblk=1,nblks
j=k1+1
do i=0,nfft1-1
cx0(i)=cmplx(dd(1,j+i),dd(2,j+i))
do i=1,NFFT1
x0(i)=fac*id2(k1+i)
enddo
call timf2x(k,nfft1,nwindow,nb,peaklimit,faclim,cx0,cx1, &
slimit,lstrong,px,nzap)
! call timf2x(k,NFFT1,nwindow,nb,peaklimit,faclim,x0,x1, &
! slimit,lstrong,px,nzap)
x1=x0 !###
x2=x1
call four2a(x2,NFFT2,1,-1,0) !Second forward FFT, r2c
do i=0,kstep-1
dd(1,j+i)=real(cx1(i))
enddo
k1=k1+kstep
i0=nint(1000.0/df2) + 1
cx2a(1:NFFT2A/2)=cx2(i0:NFFT2A/2+i0-1)
cx2a(NFFT2A/2+1:NFFT2A)=cx2(i0-1-NFFT2A/2:i0-1)
call four2a(cx2a,NFFT2A,1,1,1)
c0(k8+1:k8+NFFT2A)=cx2a
!### Test for gliches at multiples of 128
! if(k8.lt.1000) then
! do i=k8+1,k8+NFFT2A
! write(82,4002) i,c0(i)
!4002 format(i8,2e12.3)
! enddo
! endif
!###
k1=k1+kstep1
k8=k8+kstep1/8
enddo
ts=ts+hsym
ja=ts !Index of first sample
jb=ja+nfft3-1 !Last sample
i=0
fac=0.0002
do j=ja,jb !Copy data into cx
x1=dd(1,j)
i=i+1
cx(i)=fac*cmplx(x1,x2)
ja=ja+jstep !Index of first sample
if(ja.lt.0) go to 999
do i=1,nfft3 !Copy data into cx
cx(i)=c0(ja+i)
enddo
if(nzap/178.lt.50 .and. (ndiskdat.eq.0 .or. ihsym.lt.280)) then
nsum=nblks*kstep - nzap
if(nsum.le.0) nsum=1
rmsx=sqrt(0.5*px/nsum)
rms=rmsx
endif
pxdb=0.
if(rmsx.gt.1.0) pxdb=20.0*log10(rmsx)
if(pxdb.gt.60.0) pxdb=60.0
cx00=cx
ihsym=ihsym+1
cx=w*cx00 !Apply window for 2nd forward FFT
call four2a(cx,nfft3,1,1,1) !Second forward FFT (X)
call four2a(cx,nfft3,1,-1,1) !Third forward FFT (X)
n=min(322,ihsym)
do i=1,nfft3
n=min(184,ihsym)
df3=1500.0/nfft3
iz=min(NSMAX,nint(1000.0/df3))
do i=1,iz
sx=real(cx(i))**2 + aimag(cx(i))**2
ss(1,n,i)=sx
savg(1,i)=savg(1,i) + sx
ssz5a(i)=sx
ss(n,i)=sx
savg(i)=savg(i) + sx
s(i)=sx
enddo
if(ihsym.eq.175) then
do i=1,iz
write(71,3001) i*df3,savg(i),10.0*log10(savg(i))
3001 format(f12.6,e12.3,f12.3)
enddo
endif
999 return
end subroutine symspec
end subroutine symspecx

View File

@ -1,4 +1,4 @@
//--------------------------------------------------------------- MainWindow
//-------------------------------------------------------------- MainWindow
#include "mainwindow.h"
#include "ui_mainwindow.h"
#include "devsetup.h"