diff --git a/libm65/jt9.f90 b/libm65/jt9.f90 new file mode 100644 index 000000000..ced01e62a --- /dev/null +++ b/libm65/jt9.f90 @@ -0,0 +1,130 @@ +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 + + nargs=iargc() + if(nargs.lt.1) then + print*,'Usage: jt9 TRp file1 [file2 ...]' + print*,' Reads data from *.wav files.' + print*,'' + print*,' jt9 -s' + print*,' Gets data from shared memory region.' + go to 999 + endif + call getarg(1,arg) + if(arg(1:2).eq.'-s') then + call m65a + 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 + lu=12 + nfa=100 + nfb=162 + nfshift=6 + ndepth=2 + nfcal=344 + idphi=-50 + ntol=500 + nkeep=10 + + 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(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' + + 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 + 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 timer('symspec ',1) + nhsym0=nhsym + if(ihsym.ge.278) 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) + enddo + + call timer('m65 ',1) + call timer('m65 ',101) + call ftnquit + go to 999 + +998 print*,'Cannot open file:' + print*,infile + +999 end program m65 diff --git a/libm65/jt9a.f90 b/libm65/jt9a.f90 new file mode 100644 index 000000000..2f5db5562 --- /dev/null +++ b/libm65/jt9a.f90 @@ -0,0 +1,97 @@ +subroutine m65a + +! NB: this interface block is required by g95, but must be omitted +! for gfortran. (????) + +#ifndef UNIX + interface + function address_m65() + end function address_m65 + end interface +#endif + + integer*1 attach_m65,lock_m65,unlock_m65 + integer size_m65 + integer*1, pointer :: address_m65,p_m65 + character*80 cwd + logical fileExists + common/tracer/limtrace,lu + + call getcwd(cwd) + call ftninit(trim(cwd)) + limtrace=0 + lu=12 + i1=attach_m65() + +10 inquire(file=trim(cwd)//'/.lock',exist=fileExists) + if(fileExists) then + call sleep_msec(100) + go to 10 + endif + + inquire(file=trim(cwd)//'/.quit',exist=fileExists) + if(fileExists) then + call ftnquit + i=detach_m65() + go to 999 + endif + + nbytes=size_m65() + if(nbytes.le.0) then + print*,'m65a: Shared memory mem_m65 does not exist.' + print*,'Program m65a should be started automatically from within map65.' + go to 999 + endif + p_m65=>address_m65() + call m65b(p_m65,nbytes) + + write(*,1010) +1010 format('') + flush(6) + +100 inquire(file=trim(cwd)//'/.lock',exist=fileExists) + if(fileExists) go to 10 + call sleep_msec(100) + go to 100 + +999 return +end subroutine m65a + +subroutine m65b(m65com,nbytes) + integer*1 m65com(0:nbytes-1) + kss=4*4*60*96000 + ksavg=kss+4*4*322*32768 + kfcenter=ksavg+4*4*32768 + call m65c(m65com(0),m65com(kss),m65com(ksavg),m65com(kfcenter)) + return +end subroutine m65b + +subroutine m65c(dd,ss,savg,nparams0) + integer*1 detach_m65 + real*4 dd(4,5760000),ss(4,322,32768),savg(4,32768) + real*8 fcenter + integer nparams0(37),nparams(37) + character*12 mycall,hiscall + character*6 mygrid,hisgrid + character*20 datetime + 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 + equivalence (nparams,fcenter) + + nparams=nparams0 !Copy parameters into common/npar/ + npatience=1 + if(iand(nrxlog,1).ne.0) then + write(21,1000) datetime(:17) +1000 format(/'UTC Date: 'a17/78('-')) + flush(21) + endif + if(iand(nrxlog,2).ne.0) rewind 21 + if(iand(nrxlog,4).ne.0) rewind 26 + + nstandalone=0 + if(sum(nparams).ne.0) call decode0(dd,ss,savg,nstandalone) + + return +end subroutine m65c diff --git a/libm65/symspec.f90 b/libm65/symspec.f90 index a317d99fe..c8c8a3198 100644 --- a/libm65/symspec.f90 +++ b/libm65/symspec.f90 @@ -1,55 +1,52 @@ -subroutine 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) +subroutine symspecx(k,nsps,ndiskdat,nb,nbslider,pxdb,s,ihsym, & + nzap,slimit,lstrong) ! k pointer to the most recent new data -! nxpol 0/1 to indicate single- or dual-polarization +! nsps samples per symbol (at 12000 Hz) ! ndiskdat 0/1 to indicate if data from disk -! nb 0/1 status of noise blanker -! idphi Phase correction for Y channel, degrees -! nfsample sample rate (Hz) -! fgreen Frequency of green marker in I/Q calibrate mode (-48.0 to +48.0 kHz) -! iqadjust 0/1 to indicate whether IQ adjustment is active -! iqapply 0/1 to indicate whether to apply I/Q calibration -! pxdb power in x channel (0-60 dB) -! pydb power in y channel (0-60 dB) -! ssz5a polarized spectrum, for waterfall display -! nkhz integer kHz portion of center frequency, e.g., 125 for 144.125 +! 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 - parameter (NSMAX=60*96000) !Total sample intervals per minute - parameter (NFFT=32768) !Length of FFTs + 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 - common/datcom/dd(4,5760000),ss(4,322,NFFT),savg(4,NFFT),fcenter,nutc,junk(34) - real*4 ssz5a(NFFT),w(NFFT) + real*4 s(NFFT),w(NFFT) complex z,zfac - complex zsumx,zsumy - complex cx(NFFT),cy(NFFT) + complex zsumx + complex cx(NFFT) + complex cx00(NFFT) complex cx0(0:1023),cx1(0:1023) - complex cy0(0:1023),cy1(0:1023) logical*1 lstrong(0:1023) - data rms/999.0/,k0/99999999/,nadjx/0/,nadjy/0/ + common/jt8com/id2(NMAX),ss(184,NSMAX),savg(NSMAX),fcenter,nutc,junk(20) + data rms/999.0/,k0/99999999/,ntrperiod0/0/ save - if(k.gt.5751000) go to 999 - if(k.lt.NFFT) then + nfft3=nsps + hsym=nsps/2 + if(k.gt.NMAX) go to 999 + if(k.lt.nfft3) then ihsym=0 go to 999 !Wait for enough samples to start endif if(k0.eq.99999999) then pi=4.0*atan(1.0) - do i=1,NFFT - w(i)=(sin(i*pi/NFFT))**2 + do i=1,nfft3 + w(i)=(sin(i*pi/nfft3))**2 !Window for nfft3 enddo endif + if(k.lt.k0) then ts=1.d0 - hsym savg=0. ihsym=0 k1=0 - if(ndiskdat.eq.0) dd(1:4,k+1:5760000)=0. !### Should not be needed ??? ### + if(ndiskdat.eq.0) id2(k+1:)=0. !### Should not be needed ??? ### endif k0=k @@ -58,138 +55,64 @@ subroutine symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample,fgreen, & peaklimit=sigmas*max(10.0,rms) faclim=3.0 px=0. - py=0. - iqapply0=0 - iqadjust0=0 - if(iqadjust.ne.0) iqapply0=0 nwindow=2 -! nwindow=0 !### No wondowing ### - nfft2=1024 - kstep=nfft2 - if(nwindow.ne.0) kstep=nfft2/2 +! nwindow=0 !### No windowing ### + nfft1=1024 + kstep=nfft1 + if(nwindow.ne.0) kstep=nfft1/2 nblks=(k-k1)/kstep do nblk=1,nblks j=k1+1 - do i=0,nfft2-1 + do i=0,nfft1-1 cx0(i)=cmplx(dd(1,j+i),dd(2,j+i)) - if(nxpol.ne.0) cy0(i)=cmplx(dd(3,j+i),dd(4,j+i)) enddo - call timf2(k,nxpol,nfft2,nwindow,nb,peaklimit,iqadjust0,iqapply0,faclim, & - cx0,cy0,gainx,gainy,phasex,phasey,cx1,cy1,slimit,lstrong, & - px,py,nzap) + call timf2x(k,nfft1,nwindow,nb,peaklimit,faclim,cx0,cx1, & + slimit,lstrong,px,nzap) do i=0,kstep-1 dd(1,j+i)=real(cx1(i)) - dd(2,j+i)=aimag(cx1(i)) - if(nxpol.ne.0) then - dd(3,j+i)=real(cy1(i)) - dd(4,j+i)=aimag(cy1(i)) - endif enddo k1=k1+kstep enddo - hsym=2048.d0*96000.d0/11025.d0 !Samples per JT65 half-symbol - if(nfsample.eq.95238) hsym=2048.d0*95238.1d0/11025.d0 - npts=NFFT !Samples used in each half-symbol FFT - - ihsym=ihsym+1 - ja=ts+hsym !Index of first sample - jb=ja+npts-1 !Last sample - ts=ts+hsym + ja=ts !Index of first sample + jb=ja+nfft3-1 !Last sample + i=0 fac=0.0002 - dphi=idphi/57.2957795 - zfac=fac*cmplx(cos(dphi),sin(dphi)) - do j=ja,jb !Copy data into cx, cy + do j=ja,jb !Copy data into cx x1=dd(1,j) - x2=dd(2,j) - if(nxpol.ne.0) then - x3=dd(3,j) - x4=dd(4,j) - else - x3=0. - x4=0. - endif i=i+1 cx(i)=fac*cmplx(x1,x2) - cy(i)=zfac*cmplx(x3,x4) !NB: cy includes dphi correction 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) - rmsy=sqrt(0.5*py/nsum) rms=rmsx - if(nxpol.ne.0) rms=sqrt((px+py)/(4.0*nsum)) endif pxdb=0. - pydb=0. if(rmsx.gt.1.0) pxdb=20.0*log10(rmsx) - if(rmsy.gt.1.0) pydb=20.0*log10(rmsy) if(pxdb.gt.60.0) pxdb=60.0 - if(pydb.gt.60.0) pydb=60.0 - cx=w*cx !Apply window for 2nd forward FFT - if(nxpol.ne.0) cy=w*cy + cx00=cx - call four2a(cx,NFFT,1,1,1) !Second forward FFT - if(iqadjust.eq.0) nadjx=0 - if(iqadjust.ne.0 .and. nadjx.lt.50) call iqcal(nadjx,cx,NFFT,gainx,phasex, & - zsumx,ipkx,rejectx0) - if(iqapply.ne.0) call iqfix(cx,NFFT,gainx,phasex) - - if(nxpol.ne.0) then - call four2a(cy,NFFT,1,1,1) - if(iqadjust.eq.0) nadjy=0 - if(iqadjust.ne.0 .and. nadjy.lt.50) call iqcal(nadjy,cy,NFFT,gainy,phasey,& - zsumy,ipky,rejecty) - if(iqapply.ne.0) call iqfix(cy,NFFT,gainy,phasey) - endif - - n=ihsym - do i=1,NFFT - sx=real(cx(i))**2 + aimag(cx(i))**2 - ss(1,n,i)=sx ! Pol = 0 - savg(1,i)=savg(1,i) + sx + ihsym=ihsym+1 + cx=w*cx00 !Apply window for 2nd forward FFT - if(nxpol.ne.0) then - z=cx(i) + cy(i) - s45=0.5*(real(z)**2 + aimag(z)**2) - ss(2,n,i)=s45 ! Pol = 45 - savg(2,i)=savg(2,i) + s45 + call four2a(cx,nfft3,1,1,1) !Second forward FFT (X) - sy=real(cy(i))**2 + aimag(cy(i))**2 - ss(3,n,i)=sy ! Pol = 90 - savg(3,i)=savg(3,i) + sy - - z=cx(i) - cy(i) - s135=0.5*(real(z)**2 + aimag(z)**2) - ss(4,n,i)=s135 ! Pol = 135 - savg(4,i)=savg(4,i) + s135 - - z=cx(i)*conjg(cy(i)) - q=sx - sy - u=2.0*real(z) - ssz5a(i)=0.707*sqrt(q*q + u*u) !Spectrum of linear polarization -! Leif's formula: -! ssz5a(i)=0.5*(sx+sy) + (real(z)**2 + aimag(z)**2 - sx*sy)/(sx+sy) - else - ssz5a(i)=sx - endif + n=min(322,ihsym) + do i=1,nfft3 + sx=real(cx(i))**2 + aimag(cx(i))**2 + ss(1,n,i)=sx + savg(1,i)=savg(1,i) + sx + ssz5a(i)=sx enddo - if(ihsym.eq.278) then - if(iqadjust.ne.0 .and. ipkx.ne.0 .and. ipky.ne.0) then - rejectx=10.0*log10(savg(1,1+nfft-ipkx)/savg(1,1+ipkx)) - rejecty=10.0*log10(savg(3,1+nfft-ipky)/savg(3,1+ipky)) - endif - endif - nkhz=nint(1000.d0*(fcenter-int(fcenter))) - if(fcenter.eq.0.d0) nkhz=125 999 return end subroutine symspec diff --git a/libm65/timf2x.f90 b/libm65/timf2x.f90 index cf0ccdf67..5a8328a70 100644 --- a/libm65/timf2x.f90 +++ b/libm65/timf2x.f90 @@ -1,13 +1,13 @@ -subroutine timf2x(k,nfft,ntrperiod,nwindow,nb,peaklimit,faclim,cx0,cx1, & +subroutine timf2(k,nfft,nwindow,nb,peaklimit,faclim,cx0,cx1, & slimit,lstrong,px,nzap) ! Sequential processing of time-domain I/Q data, using Linrad-like ! "first FFT" and "first backward FFT". -! cx0 - complex input data -! nfft - length of FFTs -! nwindow - 0 for no window, 2 for sin^2 window -! cx1 - output data +! cx0 - complex input data +! nfft - length of FFTs +! nwindow - 0 for no window, 2 for sin^2 window +! cx1 - output data ! Non-windowed processing means no overlap, so kstep=nfft. ! Sin^2 window has 50% overlap, kstep=nfft/2. @@ -15,24 +15,23 @@ subroutine timf2x(k,nfft,ntrperiod,nwindow,nb,peaklimit,faclim,cx0,cx1, & ! Frequencies with strong signals are identified and separated. Back ! transforms are done separately for weak and strong signals, so that ! noise blanking can be applied to the weak-signal portion. Strong and -! weak signals are finally re-combined in the time domain. +! weak are finally re-combined, in the time domain. - parameter (MAXFFT=32768,MAXNH=MAXFFT/2) + parameter (MAXFFT=1024,MAXNH=MAXFFT/2) parameter (MAXSIGS=100) complex cx0(0:nfft-1),cx1(0:nfft-1) complex cx(0:MAXFFT-1),cxt(0:MAXFFT-1) complex cxs(0:MAXFFT-1),covxs(0:MAXNH-1) !Strong X signals complex cxw(0:MAXFFT-1),covxw(0:MAXNH-1) !Weak X signals - complex cxw2(0:8191) - complex cxs2(0:8191) real*4 w(0:MAXFFT-1) - real*4 s(0:MAXFFT-1) + real*4 s(0:MAXFFT-1),stmp(0:MAXFFT-1) logical*1 lstrong(0:MAXFFT-1),lprev integer ia(MAXSIGS),ib(MAXSIGS) + complex h,u,v logical first data first/.true./ data k0/99999999/ - save + save w,covxs,covxw,s,ntc,ntot,nh,kstep,fac,first,k0 if(first) then pi=4.0*atan(1.0) @@ -40,10 +39,9 @@ subroutine timf2x(k,nfft,ntrperiod,nwindow,nb,peaklimit,faclim,cx0,cx1, & w(i)=(sin(i*pi/nfft))**2 enddo s=0. + ntc=0 + ntot=0 nh=nfft/2 - nfft2=nfft/4 - if(ntrperiod.ge.300) nfft2=nfft/32 - nh2=nfft2/2 kstep=nfft if(nwindow.eq.2) kstep=nh fac=1.0/nfft @@ -59,17 +57,28 @@ subroutine timf2x(k,nfft,ntrperiod,nwindow,nb,peaklimit,faclim,cx0,cx1, & cx(0:nfft-1)=cx0 if(nwindow.eq.2) cx(0:nfft-1)=w(0:nfft-1)*cx(0:nfft-1) - call four2a(cx,nfft,1,-1,0) !First forward FFT, r2c + call four2a(cx,nfft,1,1,1) !First forward FFT + cxt(0:nfft-1)=cx(0:nfft-1) ! Identify frequencies with strong signals, copy frequency-domain ! data into array cs (strong) or cw (weak). + ntot=ntot+1 + if(mod(ntot,128).eq.5) then + call pctile(s,stmp,1024,50,xmedian) + slimit=faclim*xmedian + endif + + if(ntc.lt.96000/nfft) ntc=ntc+1 + uu=1.0/ntc + smax=0. do i=0,nfft-1 - s(i)=real(cxt(i))**2 + aimag(cxt(i))**2 + p=real(cxt(i))**2 + aimag(cxt(i))**2 + s(i)=(1.0-uu)*s(i) + uu*p + lstrong(i)=(s(i).gt.slimit) + if(s(i).gt.smax) smax=s(i) enddo - ave=sum(s(0:nfft-1))/nfft - lstrong(0:nfft-1)=s(0:nfft-1).gt.10.0*ave nsigs=0 lprev=.false. @@ -105,26 +114,19 @@ subroutine timf2x(k,nfft,ntrperiod,nwindow,nb,peaklimit,faclim,cx0,cx1, & cxs(i)=fac*cxt(i) cxw(i)=0. else - cxs(i)=0. cxw(i)=fac*cxt(i) + cxs(i)=0. endif enddo - df=12000.0/nfft - i0=nint(1500.0/df) - cxw2(0:nh2)=cxw(i0:i0+nh2) - cxw2(nfft2-nh2:nfft2-1)=cxw(i0-nh2:i0-1) - cxs2(0:nh2)=cxs(i0:i0+nh2) - cxs2(nfft2-nh2:nfft2-1)=cxs(i0-nh2:i0-1) - - call four2a(cxw2,nfft2,1,1,1) !Transform weak and strong X - call four2a(cxs2,nfft2,1,1,1) !back to time domain, separately + call four2a(cxw,nfft,1,-1,1) !Transform weak and strong X + call four2a(cxs,nfft,1,-1,1) !back to time domain, separately if(nwindow.eq.2) then - cxw2(0:nh2-1)=cxw2(0:nh2-1)+covxw(0:nh2-1) !Add prev segment's 2nd half - covxw(0:nh2-1)=cxw2(nh2:nfft2-1) !Save 2nd half - cxs2(0:nh2-1)=cxs2(0:nh2-1)+covxs(0:nh2-1) !Ditto for strong signals - covxs(0:nh2-1)=cxs2(nh2:nfft2-1) + cxw(0:nh-1)=cxw(0:nh-1)+covxw(0:nh-1) !Add previous segment's 2nd half + covxw(0:nh-1)=cxw(nh:nfft-1) !Save 2nd half + cxs(0:nh-1)=cxs(0:nh-1)+covxs(0:nh-1) !Ditto for strong signals + covxs(0:nh-1)=cxs(nh:nfft-1) endif ! Apply noise blanking to weak data @@ -132,19 +134,18 @@ subroutine timf2x(k,nfft,ntrperiod,nwindow,nb,peaklimit,faclim,cx0,cx1, & do i=0,kstep-1 peak=abs(cxw(i)) if(peak.gt.peaklimit) then - cxw2(i)=0. + cxw(i)=0. nzap=nzap+1 endif enddo endif ! Compute power levels from weak data only - px=0. do i=0,kstep-1 - px=px + real(cxw2(i))**2 + aimag(cxw2(i))**2 + px=px + real(cxw(i))**2 + aimag(cxw(i))**2 enddo - cx1(0:kstep-1)=cxw2(0:kstep-1) + cxs2(0:kstep-1) !Weak + strong + cx1(0:kstep-1)=cxw(0:kstep-1) + cxs(0:kstep-1) !Recombine weak + strong return -end subroutine timf2x +end subroutine timf2