From e72a7663372d305801a3ed0bcc58bcbb5a088d65 Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Mon, 1 Oct 2012 19:05:22 +0000 Subject: [PATCH] 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 --- jt9.txt | 46 ++++++------- libm65/jt9.f90 | 144 +++++++++++++++++---------------------- libm65/symspec.f90 | 166 +++++++++++++++++++++++++++------------------ mainwindow.cpp | 2 +- 4 files changed, 184 insertions(+), 174 deletions(-) diff --git a/jt9.txt b/jt9.txt index 5aeb0dee1..05699bd46 100644 --- a/jt9.txt +++ b/jt9.txt @@ -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. diff --git a/libm65/jt9.f90 b/libm65/jt9.f90 index ced01e62a..2de0cfe6c 100644 --- a/libm65/jt9.f90 +++ b/libm65/jt9.f90 @@ -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 diff --git a/libm65/symspec.f90 b/libm65/symspec.f90 index c8c8a3198..7bb512b33 100644 --- a/libm65/symspec.f90 +++ b/libm65/symspec.f90 @@ -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 diff --git a/mainwindow.cpp b/mainwindow.cpp index b4b246613..48d16902e 100644 --- a/mainwindow.cpp +++ b/mainwindow.cpp @@ -1,4 +1,4 @@ -//--------------------------------------------------------------- MainWindow +//-------------------------------------------------------------- MainWindow #include "mainwindow.h" #include "ui_mainwindow.h" #include "devsetup.h"