From fba5846b49ba0032b8a73ba953895101390b7a3e Mon Sep 17 00:00:00 2001
From: Joe Taylor <k1jt@arrl.org>
Date: Fri, 28 Sep 2012 18:07:07 +0000
Subject: [PATCH] Program jt9sim is now basically functional.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@2603 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
---
 jt9.txt                | 128 ++++++++++++++++++++---------------------
 libm65/genjt9.f90      |  94 ++++++++++++++++--------------
 libm65/graycode.f      |  10 ----
 libm65/graycode.f90    |   9 +++
 libm65/igray.c         |   4 --
 libm65/interleave9.f90 |  39 +++++++++++++
 libm65/jt9sim.f90      |  79 ++++++++++++++-----------
 mainwindow.cpp         |   2 +-
 8 files changed, 209 insertions(+), 156 deletions(-)
 delete mode 100644 libm65/graycode.f
 create mode 100644 libm65/graycode.f90
 create mode 100644 libm65/interleave9.f90

diff --git a/jt9.txt b/jt9.txt
index 7e554e864..a4e11239e 100644
--- a/jt9.txt
+++ b/jt9.txt
@@ -1,65 +1,63 @@
-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,
-augmented by a 12-bit cyclic redundancy check (CRC).  Error-control
-coding uses a convolutional code with constraint length K=16, rate
-r=1/2, and a zero tail, leading to an encoded message length of
-(72+12+15)*2 = 198 bits.  Modulation is 9-FSK -- 8 tones for data, and
-one for a synchronization vector.  With 15 symbol intervals used for
-synchronization, a transmission requires a total of 198/3 + 15 = 81
-channel symbols.  Tone spacing df of the 9-FSK modulation is equal to
-the keying rate; symbol duration tsym = 1/df, and the total occupied
-bandwidth is 9*df.  Transmission length is slightly less than the T/R
-sequence time, to allow for message decoding at time Tdec, a time 
-Tfree seconds before the next transmission starts.
-
-Parameters of the 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.  Parameter nsps1 is the number
-of samples per symbol at sample rate 12000 Hz; nsps is the same 
-quantity at 750 Hz.
-
-------------------------------------------------------------------------
-Mode    nsps1  nsps   df    tsym  BW    S/N*   Tdec  Tfree   Factors 
-        12000   750  (Hz)   (s)  (Hz)   (dB)   (s)    (s)    of nsps
-------------------------------------------------------------------------
-JT9-1    7168   448  1.674  0.60 15.1  -26.9   51.9   8.1   2^10 7
-JT9-2   16000  1000  0.750  1.33  6.8  -30.2  111.5   8.5   2^7 5^3
-JT9-5   42336  2646  0.283  3.53  2.5  -34.4  289.4  10.6   2^5 3^3 7^2
-JT9-10  86400  5400  0.139  7.20  1.3  -37.5  586.7  13.3   2^7 3^3 5^2
-JT9-30 262144 16384  0.046 21.85  0.4  -42.3 1773.4  26.6   2^18
-------------------------------------------------------------------------
-* Noise power measured in a 2500 Hz bandwidth.
-
-
-Transmitting
-------------
-1. Source encode the structured message to 72 bits
-2. Add 12-bit CRC
-3. Convolutionally encode (K=16, r=1/2) to (72+12+15)*2 = 198 bits
-4. Interleave to scramble the bit order
-5. Assemble 3-bit groups to make 198/3 = 66 symbols
-6. Apply Gray code to the symbol values
-7. Insert 15 sync symbols ==> 66+15=81 channel symbols with values 0-8
-
-
-Receiving
----------
-1.  Apply noise blanking with the timf2 method
-2.  Filter to 500 Hz bandwidth, downsample (1/16) to 750 Hz complex
-    using FIR filter with 61 taps.  Data in array c0(1800*750) (1.35M)
-3.  Compute symbol-length spectra at half-symbol steps.  Use for
-    waterfalls s(16384) and save in ss(180,16384), savg(16384)
-4.  At time Tdec, look for sync vectors in ss() to get  estimates of DF, DT
-5.  Do full-length FFTs of length NFFT1=96*nsps
-6.  For each candidate JT9 signal, do inverse FFT of length 1536.  This 
-    gives 16 complex samples per symbol, sync tone should be close to f=0.
-7.  Use afc65b method to get improved values of DF, DT.
-8.  Tweak freq and time offset to 0.
-9.  Compute 8-bin spectra of 66 data symbols: s2(8,66).  Re-order the
-    bins by removing Gray code.
-10. Compute soft symbols for 198 bits
-11. Re-order the bits by removing interleaving.
-12. Pack bits into bytes, send to Viterbi decoder ==> 72-bit message, 
-    12-bit CRC, and metric.
-13. Declare erasure if CRC check fails
-14. If CRC is good, remove source encoding and display user message.
+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 set to approximately
+(TRperiod-10)/85, where TRperiod is the T/R sequence length in
+seconds.  The exact symbol lengths are chosen so as to make nsps, the
+number of samples at 12000 samples per second, a number rich in
+factors less than 7.  This choice makes for efficient FFTs.  Tone
+spacing of the 9-FSK modulation is df=1/tsym, the same as the keying
+rate.  The total occupied bandwidth is 9*df.
+
+Parameters of the 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.  
+
+-------------------------------------------------------------------------
+Mode     nsps nsps2   df    tsym  BW    S/N*   Tdec  Tfree   Factors 
+        12000   750  (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
+-------------------------------------------------------------------------
+* Noise power measured in a 2500 Hz bandwidth.
+
+
+Transmitting
+------------
+1. Source encode the structured message to 72 bits
+2. Convolutionally encode (K=32, r=1/2) to (72+31)*2 = 206 bits
+3. Interleave to scramble the bit order
+4. Assemble 3-bit groups to make (206+1)/3 = 69 symbols
+5. Gray-code the symbol values
+6. Insert 16 sync symbols ==> 69+15=81 channel symbols, values 0-8
+
+
+Receiving
+---------
+1.  Apply noise blanking with the timf2 method
+2.  Filter to 500 Hz bandwidth, downsample (1/16) to 750 Hz complex
+    (FIR filter with 61 taps).  Complex data to array c0 (max 1.35M)
+3.  Compute symbol-length spectra at half-symbol steps.  Use for
+    waterfalls s(15750) and detecting sync vectors; save in ss(184,15750)
+    and savg(15750)
+4.  At time Tdec, look for sync vectors in ss(); get estimates of DF, DT
+5.  Do full-length FFTs, NFFT1=96*nsps2
+6.  For each candidate signal, do inverse FFT of length 1536.  This 
+    will provide 16 complex samples per symbol, and sync tone should be 
+    close to f=0.
+7.  Use afc65b method to get improved values of DF, DT.
+8.  Tweak freq and time offset to 0.
+9.  Compute 8-bin spectra of 69 data symbols: s2(8,69).  Re-order the
+    bins by removing Gray code.
+10. Compute soft symbols for 206 bits
+11. Re-order the bits by removing interleaving.
+12. Pack bits into bytes, send to Fano decoder
+13. If Fano succeeds, remove source encoding and display user message.
diff --git a/libm65/genjt9.f90 b/libm65/genjt9.f90
index a1b60f64a..d87a5a2a5 100644
--- a/libm65/genjt9.f90
+++ b/libm65/genjt9.f90
@@ -1,77 +1,86 @@
-subroutine genjt9(message,minutes,lwave,msgsent,d7,iwave,nwave)
+subroutine genjt9(message,minutes,lwave,msgsent,d8,iwave,nwave)
 
-! Encodes a "JT9-minutes" message and returns array d7(81) of tone
+! Encodes a "JT9-minutes" message and returns array d7(85) of tone
 ! values in the range 0-8.  If lwave is true, also returns a generated
 ! waveform in iwave(nwave), at 12000 samples per second.
 
   parameter (NMAX=1800*12000)   !Max length of wave file
   character*22 message          !Message to be generated
   character*22 msgsent          !Message as it will be received
-  character*3 cok
   real*8 dt,phi,f,f0,dfgen,dphi,twopi
   integer*2 iwave(NMAX)         !Generated wave file
-  integer   d0(14)              !72-bit message + 12-bit CRC, as 6-bit bytes
-  integer*1 d1(84)              !72+12 = 84 single bits 
-  integer   d2(11)              !72+12 bits, as 8-bit words
-  integer*1 d3(11)              !72+12 bits, as 8-bit bytes
-  integer*1 d4(198)             !Encoded information-carrying bits
-  integer*1 d5(198)             !Bits from d4, after interleaving
-  integer*1 d6(66)              !Symbols from d5, values 0-7
-  integer*1 d7(66)              !Gray-coded symbols, values 0-7
-  integer*1 d8(81)              !Channel symbols including sync, values 0-8
+
+  integer*4 d0(12)              !72-bit message as 6-bit words
+  integer*1 d1(72)              !72 single bits 
+  integer*4 d2(9)               !72 bits as 8-bit words
+  integer*1 d3(9)               !72 bits as 8-bit bytes
+  integer*1 d4(206)             !Encoded information-carrying bits
+  integer*1 d5(206)             !Bits from d4, after interleaving
+  integer*4 d6(69)              !Symbols from d5, values 0-7
+  integer*4 d7(69)              !Gray-coded symbols, values 0-7
+  integer*4 d8(85)              !Channel symbols including sync, values 0-8
 
   logical lwave
-  integer isync(15)
-  data isync/1,6,11,16,21,26,31,39,45,51,57,63,69,75,81/
+  integer isync(85)
+  data isync/                                    &
+       1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,  &
+       1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,  &
+       0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,  &
+       0,0,1,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,  &
+       1,0,0,0,1/
   data twopi/6.283185307179586476d0/
   save
 
-  call chkmsg(message,cok,nspecial,flip)
   call packmsg(message,d0)            !Pack message into 12 6-bit bytes
-  d0(11)=0                            !### Temporary CRC=0 ###
-  d0(12)=0
-  call unpackbits(d0,14,6,d1)         !Unpack into 84 bits
-  nbits=84
-  call packbits(d1,nbits,8,d2)        !Pack into 11 8-bit words
-  nbytes=(nbits+7)/8
-  do i=1,nbytes
+  call unpackbits(d0,12,6,d1)         !Unpack into 72 bits
+
+  nbits=72
+  nbytes=9
+  call packbits(d1,nbits,8,d2)        !Pack into 9 8-bit words
+  do i=1,9
      if(d2(i).lt.128) d3(i)=d2(i)
      if(d2(i).ge.128) d3(i)=d2(i)-256
   enddo
 
-  call enc216(d3,nbits,d4,nsym2,16,2)   !Convolutional code, K=16, r=1/2
-! NB: Should give nsym2 = (72+12+15)*2 = 198
-
-!  call interleavejt9(d4,1,d5)
-  d5=d4                             !### TEMP ###
-
+  call encode232(d3,nbytes,d4)        !Convolutional code, K=32, r=1/2
+  nsym2=206
+  call interleave9(d4,1,d5)
   call packbits(d5,nsym2,3,d6)
-!  call graycode(d6,nsym2,1,d7)
+  call graycode(d6,69,1,d7)
 
-! Now insert sync symbols and add 1 to the tone numbers.
+! Insert sync symbols (ntone=0) and add 1 to the data-tone numbers.
+  j=0
+  do i=1,85
+     if(isync(i).eq.1) then
+        d8(i)=0
+     else
+        j=j+1
+        d8(i)=d7(j)+1
+     endif
+  enddo
 
   nwave=0
 
   if(lwave) then
-     nsps1=0
-     if(minutes.eq.1)  nsps1=7168
-     if(minutes.eq.2)  nsps1=16000
-     if(minutes.eq.5)  nsps1=42336
-     if(minutes.eq.10) nsps1=86400
-     if(minutes.eq.30) nsps1=262144
-     if(nsps1.eq.0) stop 'Bad value for minutes in genjt9.'
+     nsps=0
+     if(minutes.eq.1)  nsps=6912
+     if(minutes.eq.2)  nsps=15360
+     if(minutes.eq.5)  nsps=40960
+     if(minutes.eq.10) nsps=82944
+     if(minutes.eq.30) nsps=250880
+     if(nsps.eq.0) stop 'Bad value for minutes in genjt9.'
 
 ! Set up necessary constants
      dt=1.d0/12000.d0
      f0=1500.d0
-     dfgen=12000.d0/nsps1
+     dfgen=12000.d0/nsps
      phi=0.d0
      i=0
      k=0
-     do j=1,81
+     do j=1,85
         f=f0 +d7(j)*dfgen
         dphi=twopi*dt*f
-        do i=1,nsps1
+        do i=1,nsps
            phi=phi+dphi
            if(phi.gt.twopi) phi=phi-twopi
            xphi=phi
@@ -79,9 +88,10 @@ subroutine genjt9(message,minutes,lwave,msgsent,d7,iwave,nwave)
            iwave(k)=32767.0*sin(xphi)
         enddo
      enddo
-     nwave=81*nsps1
+     nwave=85*nsps
   endif
-  call unpackmsg(dgen,msgsent)
+
+  call unpackmsg(d0,msgsent)
 
   return
 end subroutine genjt9
diff --git a/libm65/graycode.f b/libm65/graycode.f
deleted file mode 100644
index 24c03e943..000000000
--- a/libm65/graycode.f
+++ /dev/null
@@ -1,10 +0,0 @@
-      subroutine graycode(dat,n,idir)
-
-      integer dat(n)
-      do i=1,n
-         dat(i)=igray(dat(i),idir)
-      enddo
-
-      return
-      end
-
diff --git a/libm65/graycode.f90 b/libm65/graycode.f90
new file mode 100644
index 000000000..21c1f9086
--- /dev/null
+++ b/libm65/graycode.f90
@@ -0,0 +1,9 @@
+subroutine graycode(ia,n,idir,ib)
+
+  integer ia(n),ib(n)
+  do i=1,n
+     ib(i)=igray(ia(i),idir)
+  enddo
+
+  return
+end subroutine graycode
diff --git a/libm65/igray.c b/libm65/igray.c
index 395f79712..4646898b5 100644
--- a/libm65/igray.c
+++ b/libm65/igray.c
@@ -1,8 +1,4 @@
-#ifdef CVF
-extern int __stdcall IGRAY(int *n0, int *idir)
-#else
 int igray_(int *n0, int *idir)
-#endif
 {
   int n;
   unsigned long sh;
diff --git a/libm65/interleave9.f90 b/libm65/interleave9.f90
new file mode 100644
index 000000000..0626cf878
--- /dev/null
+++ b/libm65/interleave9.f90
@@ -0,0 +1,39 @@
+subroutine interleave9(ia,ndir,ib)
+  integer*1 ia(0:205),ib(0:205)
+  integer j0(0:205)
+  logical first
+  data first/.true./
+  save first,j0
+
+  if(first) then
+     k=-1
+     do i=0,255
+        m=i
+        n=iand(m,1)
+        n=2*n + iand(m/2,1)
+        n=2*n + iand(m/4,1)
+        n=2*n + iand(m/8,1)
+        n=2*n + iand(m/16,1)
+        n=2*n + iand(m/32,1)
+        n=2*n + iand(m/64,1)
+        n=2*n + iand(m/128,1)
+        if(n.le.205) then
+           k=k+1
+           j0(k)=n
+        endif
+     enddo
+!     first=.false.
+  endif
+
+  if(ndir.gt.0) then
+     do i=0,205
+        ib(j0(i))=ia(i)
+     enddo
+  else
+     do i=0,205
+        ib(i)=ia(j0(i))
+     enddo
+  endif
+
+  return
+end subroutine interleave9
diff --git a/libm65/jt9sim.f90 b/libm65/jt9sim.f90
index f409ace8f..a4dc719df 100644
--- a/libm65/jt9sim.f90
+++ b/libm65/jt9sim.f90
@@ -1,19 +1,18 @@
 program jt9sim
 
-! Generate simulated data for testing of MAP65
+! Generate simulated data for testing of WSJT-X
 
   parameter (NMAX=1800*12000)
-  real*4 d4(NMAX)                     !Floating-point data
   integer ihdr(11)
-  integer*2 id2(NMAX)                 !i*2 data
-  integer*2 iwave(NMAX)               !Generated waveform (no noise)
-  real*8 f0,f,dt,twopi,phi,dphi,baud
+  integer*2 iwave                     !Generated waveform (no noise)
+  real*8 f0,f,dt,twopi,phi,dphi,baud,fspan
   character msg0*22,message*22,msgsent*22,arg*8,fname*11
   logical lwave
-  integer*1 d7(81)
+  integer*4 d8(85)
+  common/acom/dat(NMAX),iwave(NMAX)
 
   nargs=iargc()
-  if(nargs.ne.9) then
+  if(nargs.ne.6) then
      print*,'Usage: jt9sim "message" fspan nsigs minutes SNR nfiles'
      print*,'Example:  "CQ K1ABC FN42" 200  20      2    -28    1'
      print*,' '
@@ -41,21 +40,21 @@ program jt9sim
   fsample=12000.d0                   !Sample rate (Hz)
   dt=1.d0/fsample                    !Sample interval (s)
   twopi=8.d0*atan(1.d0)
-  rad=360.0/twopi
   npts=12000*(60*minutes-6)
   lwave=.false.
-  nsps1=0
-  if(minutes.eq.1)  nsps1=7168
-  if(minutes.eq.2)  nsps1=16000
-  if(minutes.eq.5)  nsps1=42336
-  if(minutes.eq.10) nsps1=86400
-  if(minutes.eq.30) nsps1=262144
-  if(nsps1.eq.0) stop 'Bad value for minutes.'
+  nsps=0
+  if(minutes.eq.1)  nsps=6912
+  if(minutes.eq.2)  nsps=15360
+  if(minutes.eq.5)  nsps=40960
+  if(minutes.eq.10) nsps=82944
+  if(minutes.eq.30) nsps=250880
+  if(nsps.eq.0) stop 'Bad value for minutes.'
+  ihdr=0                             !Temporary ###
   
   open(12,file='msgs.txt',status='old')
 
   write(*,1000)
-1000 format('File  N   freq     S/N  Message'/    &
+1000 format('File  N    freq      S/N  Message'/    &
             '---------------------------------------------------')
 
   do ifile=1,nfiles
@@ -66,10 +65,14 @@ program jt9sim
 1002 format('000000_',2i2.2)
      open(10,file=fname//'.wav',access='stream',status='unknown')
 
-     call noisegen(d4,npts)                      !Generate Gaussian noise
+     if(snrdb.lt.90) then
+        call noisegen(dat,npts)                 !Generate Gaussian noise
+     else
+        dat(1:npts)=0.
+     endif
 
      if(msg0.ne.'                      ') then
-        call genjt9(message,minutes,lwave,msgsent,d7,iwave,nwave)
+        call genjt9(message,minutes,lwave,msgsent,d8,iwave,nwave)
      endif
 
      rewind 12
@@ -78,40 +81,48 @@ program jt9sim
         if(msg0.eq.'                      ') then
            read(12,1004) message
 1004       format(a22)
-           call genjt9(message,minutes,lwave,msgsent,d7,iwave,nwave)
+           call genjt9(message,minutes,lwave,msgsent,d8,iwave,nwave)
         endif
-           
-        f=1500.d0
-        if(nsigs.gt.1) f=1500.0 - 0.5* fspan + fspan*(isig-1.0)/(nsigs-1.0)
+
+        f=f0
+        if(nsigs.gt.1) f=f0 - 0.5d0*fspan + fspan*(isig-1.d0)/(nsigs-1.d0)
         snrdbx=snrdb
-        if(snrdb.ge.-1.0) snrdbx=-15.0 - 15.0*(isig-1.0)/nsigs
-        sig=sqrt(2.2*2500.0/96000.0) * 10.0**(0.05*snrdbx)
-        write(*,1020) ifile,isig,0.001*f,snrdbx,msgsent
-1020    format(i3,i4,f8.3,f7.1,2x,a22)
+!        if(snrdb.ge.-1.0) snrdbx=-15.0 - 15.0*(isig-1.0)/nsigs
+        sig=sqrt(2500.0/12000.0) * 10.0**(0.05*snrdbx)
+        write(*,1020) ifile,isig,f,snrdbx,msgsent
+1020    format(i3,i4,f10.3,f7.1,2x,a22)
 
         phi=0.
-        baud=12000.0/nsps1
+        baud=12000.0/nsps
         k=12000                             !Start at t = 1 s
-        do isym=1,81
-           freq=f + d7(isym)*baud
-           dphi=twopi*freq*dt + 0.5*twopi
-           do i=1,nsps1
+        do isym=1,85
+           freq=f + d8(isym)*baud
+           dphi=twopi*freq*dt
+           do i=1,nsps
               phi=phi + dphi
               if(phi.lt.-twopi) phi=phi+twopi
               if(phi.gt.twopi) phi=phi-twopi
               xphi=phi
               k=k+1
-              d4(k)=d4(k) + sin(phi)
+              dat(k)=dat(k) + sig*sin(xphi)
            enddo
         enddo
      enddo
 
      do i=1,npts
-        id2(i)=nint(rms*d4(i))
+        iwave(i)=nint(rms*dat(i))
+!        iwave(i)=nint(100.0*dat(i))
      enddo
 
-     write(10) ihdr,id2(1:npts)
+     write(10) ihdr,iwave(1:npts)
      close(10)
+
+     
+     do i=1,30000
+        write(71,3001) i,iwave(12000+i)
+3001    format(2i8)
+     enddo
+
   enddo
 
 999 end program jt9sim
diff --git a/mainwindow.cpp b/mainwindow.cpp
index 2df44ab0f..d1d67f2ab 100644
--- a/mainwindow.cpp
+++ b/mainwindow.cpp
@@ -1,4 +1,4 @@
-//---------------------------------------------------------------- MainWindow
+//--------------------------------------------------------------- MainWindow
 #include "mainwindow.h"
 #include "ui_mainwindow.h"
 #include "devsetup.h"