From 6ccd2a290b5027df3747f9f9691b511a556933b5 Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Sun, 11 Dec 2022 17:10:40 -0500 Subject: [PATCH] Convert savg(4,NFFT) to one-dimensional savg(NFFT). --- q65w/commons.h | 4 +- q65w/libm65/decode0.f90 | 2 +- q65w/libm65/m65.f90 | 195 ---------------------------------- q65w/libm65/m65a.f90 | 107 ------------------- q65w/libm65/m65c.f90 | 3 +- q65w/libm65/map65a.f90 | 14 +-- q65w/libm65/recvpkt.f90 | 2 +- q65w/libm65/symspec.f90 | 37 +------ q65w/libm65/wideband_sync.f90 | 16 +-- q65w/soundin.cpp | 2 +- 10 files changed, 24 insertions(+), 358 deletions(-) delete mode 100644 q65w/libm65/m65.f90 delete mode 100644 q65w/libm65/m65a.f90 diff --git a/q65w/commons.h b/q65w/commons.h index 4f2b0dbda..792a8adaf 100644 --- a/q65w/commons.h +++ b/q65w/commons.h @@ -8,7 +8,7 @@ extern "C" { extern struct { //This is "common/datcom/..." in Fortran float d4[4*5760000]; //Raw I/Q data from Linrad float ss[4*322*NFFT]; //Half-symbol spectra at 0,45,90,135 deg pol - float savg[4*NFFT]; //Avg spectra at 0,45,90,135 deg pol + float savg[NFFT]; //Avg spectra at 0,45,90,135 deg pol double fcenter; //Center freq from Linrad (MHz) int nutc; //UTC as integer, HHMM int idphi; //Phase correction for Y pol'n, degrees @@ -47,7 +47,7 @@ extern struct { //This is "common/datcom/..." in Fortran extern struct { //This is "common/datcom/..." in Fortran float d4[4*5760000]; //Raw I/Q data from Linrad float ss[4*322*NFFT]; //Half-symbol spectra at 0,45,90,135 deg pol - float savg[4*NFFT]; //Avg spectra at 0,45,90,135 deg pol + float savg[NFFT]; //Avg spectra at 0,45,90,135 deg pol double fcenter; //Center freq from Linrad (MHz) int nutc; //UTC as integer, HHMM int idphi; //Phase correction for Y pol'n, degrees diff --git a/q65w/libm65/decode0.f90 b/q65w/libm65/decode0.f90 index 7f1ba2855..b3c017a33 100644 --- a/q65w/libm65/decode0.f90 +++ b/q65w/libm65/decode0.f90 @@ -3,7 +3,7 @@ subroutine decode0(dd,ss,savg,nstandalone) use timer_module, only: timer parameter (NSMAX=60*96000) - real*4 dd(4,NSMAX),ss(4,322,NFFT),savg(4,NFFT) + real*4 dd(4,NSMAX),ss(4,322,NFFT),savg(NFFT) real*8 fcenter integer hist(0:32768) logical ldecoded diff --git a/q65w/libm65/m65.f90 b/q65w/libm65/m65.f90 deleted file mode 100644 index c93357c7e..000000000 --- a/q65w/libm65/m65.f90 +++ /dev/null @@ -1,195 +0,0 @@ -program m65 - -! Decoder for map65. Can run stand-alone, reading data from *.tf2 files; -! or as the back end of map65, with data placed in a shared memory region. - -! Fortran logical units -! -! 10 binary input data, *.tf2 files -! 11 prefixes.txt -! 12 wb_w65.txt -! 13 map65.log -! 14 -! 15 -! 16 tquick log -! 17 saved *.tf2 files -! 18 test file to be transmitted (wsjtgen.f90) -! 19 livecq.txt -! 20 -! 21 map65_rx.log -! 22 -! 23 CALL3.TXT -! 24 -! 25 -! 26 tmp26.txt - - use timer_module, only: timer - use timer_impl, only: init_timer, fini_timer - - include 'njunk.f90' - parameter (NFFT=32768) - parameter (NSMAX=60*96000) - parameter (NREAD=2048) - integer*2 i2(NREAD) - real*8 hsym - real*4 ssz5a(NFFT) - logical*1 lstrong(0:1023),ldecoded,eof - 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(NJUNK) - common/npar/fcenter,nutc,idphi,mousedf,mousefqso,nagain, & - ndepth,ndiskdat,neme,newdat,nfa,nfb,nfcal,nfshift, & - mcall3,nkeep,ntol,nxant,nrxlog,nfsample,nxpol,nmode, & - ndop00,nsave,max_drift,nhsym,mycall,mygrid,hiscall,hisgrid,datetime - common/early/nhsym1,nhsym2,ldecoded(32768) - - nargs=iargc() - if(nargs.ne.1 .and. nargs.lt.5) then - print*,'Usage: m65 Jsub Qsub Xpol <95238|96000> file1 [file2 ...]' - print*,'Examples: m65 B A X 96000 *.tf2' - print*,' m65 C C N 96000 *.iq' - print*,'' - print*,' m65 -s' - print*,' (Gets data from MAP65, via shared memory region.)' - go to 999 - endif - nstandalone=1 - nhsym1=280 - nhsym2=302 - call getarg(1,arg) - if(arg(1:2).eq.'-s') then - call m65a - go to 999 - endif - n=1 - if(arg(1:1).eq.'0') n=0 - if(arg(1:1).eq.'A') n=1 - if(arg(1:1).eq.'B') n=2 - if(arg(1:1).eq.'C') n=3 - - call getarg(2,arg) - m=1 - if(arg(1:1).eq.'0') m=0 - if(arg(1:1).eq.'A') m=1 - if(arg(1:1).eq.'B') m=2 - if(arg(1:1).eq.'C') m=3 - if(arg(1:1).eq.'D') m=4 - if(arg(1:1).eq.'E') m=5 - nmode=10*m + n - - call getarg(3,arg) - nxpol=0 - if(arg(1:1).eq.'X') nxpol=1 - - call getarg(4,arg) - nfsample=96000 - if(arg.eq.'95238') nfsample=95238 - - ifile1=5 - -! Some default parameters for command-line execution, in early testing. - mycall='K1JT' - mygrid='FN20QI' - hiscall='K9AN' - hisgrid='EN50' - nfa=100 !144.100 - nfb=162 !144.162 - ntol=100 - nkeep=10 !??? - mousefqso=140 !For IK4WLV in 210220_1814.tf2 - mousedf=0 - nfcal=0 - nkhz_center=125 - - if(nxpol.eq.0) then - nfa=55 !For KA1GT files - nfb=143 - mousefqso=69 !W2HRO signal - nkhz_center=100 - endif - - call ftninit('.') - call init_timer('timer.out') - call timer('m65 ',0) - - do ifile=ifile1,nargs - call getarg(ifile,infile) - open(10,file=infile,access='stream',status='old',err=998) - i1=index(infile,'.tf2') - if(i1.lt.1) i1=index(infile,'.iq') - 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 - read(10) fcenter - newdat=1 - nhsym0=-999 - k=0 - - nch=2 - if(nxpol.eq.1) nch=4 - eof=.false. - do irec=1,9999999 - if(.not.eof) read(10,end=4) i2 - go to 6 -4 eof=.true. -6 if(eof) i2=0 - do i=1,NREAD,nch - k=k+1 - if(k.gt.60*96000) exit - dd(1,k)=i2(i) - dd(2,k)=i2(i+1) - if(nxpol.eq.1) then - dd(3,k)=i2(i+2) - dd(4,k)=i2(i+3) - endif - enddo - nhsym=(k-2048)/hsym - if(nhsym.ge.1 .and. nhsym.ne.nhsym0) then - ndiskdat=1 - nb=0 -! Emit signal readyForFFT - fgreen=-13.0 - iqadjust=0 - iqapply=0 - nbslider=100 - gainx=0.9962 - gainy=1.0265 - phasex=0.01426 - phasey=-0.01195 - call timer('symspec ',0) - 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 - - nutc=nutc0 - if(nhsym.eq.nhsym1) call decode0(dd,ss,savg,nstandalone) - if(nhsym.eq.nhsym2) then - call decode0(dd,ss,savg,nstandalone) - exit - endif - endif - enddo ! irec - - if(iqadjust.ne.0) write(*,3002) rejectx,rejecty -3002 format('Image rejection:',2f7.1,' dB') - enddo ! ifile - - call timer('m65 ',1) - call timer('m65 ',101) - go to 999 - -998 print*,'Cannot open file:' - print*,infile - -999 call fini_timer() - if(arg(1:2).eq.'-s') then - write(21,1999) datetime(:17) -1999 format('Subprocess m65 terminated normally at UTC ',a17) - close(21) - endif - -end program m65 diff --git a/q65w/libm65/m65a.f90 b/q65w/libm65/m65a.f90 deleted file mode 100644 index ffeb176a3..000000000 --- a/q65w/libm65/m65a.f90 +++ /dev/null @@ -1,107 +0,0 @@ -subroutine m65a - - use timer_module, only: timer - use timer_impl, only: init_timer !, limtrace - use, intrinsic :: iso_c_binding, only: C_NULL_CHAR - use FFTW3 - - interface - function address_m65() - integer*1, pointer :: address_m65 - end function address_m65 - end interface - - integer*1 attach_m65 - integer size_m65 - integer*1, pointer :: p_m65 - character*80 cwd - character wisfile*256 - logical fileExists - - call getcwd(cwd) - call ftninit(trim(cwd)) - call init_timer (trim(cwd)//'/timer.out') - - 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 timer('decode0 ',101) - i=detach_m65() - ! Save FFTW wisdom and free memory - wisfile=trim(cwd)//'/m65_wisdom.dat'// C_NULL_CHAR - if(len(trim(wisfile)).gt.0) iret=fftwf_export_wisdom_to_filename(wisfile) - call four2a(a,-1,1,1,1) - call filbig(a,-1,1,0.0,0,0,0,0,0) !used for FFT plans - call fftwf_cleanup_threads() - call fftwf_cleanup() - 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) - call sleep_msec(500) ! wait for .lock to be recreated - go to 10 - -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) - - include 'njunk.f90' - real*4 dd(4,5760000),ss(4,322,32768),savg(4,32768) - real*8 fcenter - integer nparams0(NJUNK+2),nparams(NJUNK+2) - logical ldecoded - 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,nmode, & - ndop00,nsave,max_drift,nhsym,mycall,mygrid,hiscall,hisgrid, & - datetime,junk1,junk2 - common/early/nhsym1,nhsym2,ldecoded(32768) - equivalence (nparams,fcenter) - - nparams=nparams0 !Copy parameters into common/npar/ - npatience=1 - if(nhsym.eq.nhsym1 .and. 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) then - if(nhsym.eq.nhsym1) rewind(26) - if(nhsym.eq.nhsym2) backspace(26) - endif - - nstandalone=0 - if(sum(nparams).ne.0) call decode0(dd,ss,savg,nstandalone) - - return -end subroutine m65c diff --git a/q65w/libm65/m65c.f90 b/q65w/libm65/m65c.f90 index 7853ec507..be3914bc1 100644 --- a/q65w/libm65/m65c.f90 +++ b/q65w/libm65/m65c.f90 @@ -9,7 +9,6 @@ subroutine m65c(itimer) parameter (NFFT=32768) include 'njunk.f90' -! real*4 dd(4,5760000),ss(4,322,32768),savg(4,32768) real*8 fcenter integer nparams0(NJUNK+3),nparams(NJUNK+3) logical ldecoded,first @@ -17,7 +16,7 @@ subroutine m65c(itimer) character*6 mygrid,hisgrid character*20 datetime - common/datcom2/dd(4,5760000),ss(4,322,NFFT),savg(4,NFFT),nparams0 + common/datcom2/dd(4,5760000),ss(4,322,NFFT),savg(NFFT),nparams0 !### REMEMBER that /npar/ is not updated until nparams=nparams0 is executed. ### common/npar/fcenter,nutc,idphi,mousedf,mousefqso,nagain, & diff --git a/q65w/libm65/map65a.f90 b/q65w/libm65/map65a.f90 index 756624e65..e33cff612 100644 --- a/q65w/libm65/map65a.f90 +++ b/q65w/libm65/map65a.f90 @@ -12,7 +12,7 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, & parameter (MAXMSG=1000) !Size of decoded message list parameter (NSMAX=60*96000) real dd(4,NSMAX) - real*4 ss(4,322,NFFT),savg(4,NFFT) + real*4 ss(4,322,NFFT),savg(NFFT) real tavg(-50:50) !Temp for finding local base level real base(4) !Local basel level at 4 pol'ns real sig(MAXMSG,30) !Parameters of detected signals @@ -119,7 +119,7 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, & ntry=0 short=0. !Zero the whole short array jpz=1 - if(xpol) jpz=4 +! if(xpol) jpz=4 ! First steps for JT65 decoding do i=ia,ib !Search over freq range @@ -130,7 +130,7 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, & do ii=-50,50 iii=i+ii if(iii.ge.1 .and. iii.le.32768) then - tavg(ii)=savg(jp,iii) + tavg(ii)=savg(iii) else write(13,*) 'Error in iii:',iii,ia,ib,fa,fb flush(13) @@ -144,8 +144,8 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, & ! Find max signal at this frequency smax=0. do jp=1,jpz - if(savg(jp,i)/base(jp).gt.smax) then - smax=savg(jp,i)/base(jp) + if(savg(i)/base(jp).gt.smax) then + smax=savg(i)/base(jp) jpmax=jp endif enddo @@ -197,7 +197,7 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, & sig(km,9)=0 sig(km,10)=0 ! sig(km,11)=rms0 - sig(km,12)=savg(ipol2,i) + sig(km,12)=savg(i) sig(km,13)=0 sig(km,14)=0 sig(km,15)=0 @@ -257,7 +257,7 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, & sig(km,9)=nkv sig(km,10)=qual ! sig(km,11)=idphi - sig(km,12)=savg(ipol,i) + sig(km,12)=savg(i) sig(km,13)=a(1) sig(km,14)=a(2) sig(km,15)=a(3) diff --git a/q65w/libm65/recvpkt.f90 b/q65w/libm65/recvpkt.f90 index 88d51cc59..031ac8652 100644 --- a/q65w/libm65/recvpkt.f90 +++ b/q65w/libm65/recvpkt.f90 @@ -12,7 +12,7 @@ subroutine recvpkt(nsam,nblock2,userx_no,k,buf4,buf8,buf16) integer*2 jd(4),kd(2),nblock2 real*4 xd(4),yd(2) real*8 fcenter - common/datcom/dd(4,5760000),ss(4,322,NFFT),savg(4,NFFT),fcenter,nutc, & + common/datcom/dd(4,5760000),ss(4,322,NFFT),savg(NFFT),fcenter,nutc, & junk(NJUNK) equivalence (kd,d4) equivalence (jd,d8,yd) diff --git a/q65w/libm65/symspec.f90 b/q65w/libm65/symspec.f90 index e4160cfd1..55aa3110c 100644 --- a/q65w/libm65/symspec.f90 +++ b/q65w/libm65/symspec.f90 @@ -23,7 +23,7 @@ subroutine symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample, & parameter (NFFT=32768) !Length of FFTs real*8 ts,hsym real*8 fcenter - common/datcom/dd(4,5760000),ss(4,322,NFFT),savg(4,NFFT),fcenter,nutc, & + common/datcom/dd(4,5760000),ss(4,322,NFFT),savg(NFFT),fcenter,nutc, & junk(NJUNK) real*4 ssz5a(NFFT),w(NFFT),w2a(NFFT),w2b(NFFT) complex z,zfac @@ -181,42 +181,11 @@ subroutine symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample, & 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 - - 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 - - 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 + savg(i)=savg(i) + sx + ssz5a(i)=sx enddo 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 diff --git a/q65w/libm65/wideband_sync.f90 b/q65w/libm65/wideband_sync.f90 index d99391cd3..bdbdef326 100644 --- a/q65w/libm65/wideband_sync.f90 +++ b/q65w/libm65/wideband_sync.f90 @@ -34,7 +34,7 @@ subroutine get_candidates(ss,savg,xpol,jz,nfa,nfb,nts_jt65,nts_q65,cand,ncand) ! excised. Candidates are returned in the structure array cand(). parameter (MAX_PEAKS=100) - real ss(4,322,NFFT),savg(4,NFFT) + real ss(4,322,NFFT),savg(NFFT) real pavg(-20:20) integer indx(NFFT) logical xpol,skip,ldecoded @@ -122,8 +122,8 @@ subroutine wb_sync(ss,savg,xpol,jz,nfa,nfb) parameter (NFFT=32768) parameter (LAGMAX=30) real ss(4,322,NFFT) - real savg(4,NFFT) - real savg_med(4) + real savg(NFFT) + real savg_med real a(3) logical first,xpol integer isync(22) @@ -164,7 +164,7 @@ subroutine wb_sync(ss,savg,xpol,jz,nfa,nfb) if(ia.lt.1) ia=1 if(ib.gt.NFFT-1) ib=NFFT-1 - call pctile(savg(1,ia:ib),ib-ia+1,50,savg_med(1)) + call pctile(savg(ia:ib),ib-ia+1,50,savg_med) lagbest=0 ipolbest=1 @@ -181,7 +181,7 @@ subroutine wb_sync(ss,savg,xpol,jz,nfa,nfb) ccf4=ccf4 + ss(1,k,i+1) + ss(1,k+1,i+1) & + ss(1,k+2,i+1) enddo - ccf4=ccf4 - savg(1,i+1)*3*22/float(jz) + ccf4=ccf4 - savg(i+1)*3*22/float(jz) ccf=ccf4 ipol=1 if(ccf.gt.ccfmax) then @@ -198,7 +198,7 @@ subroutine wb_sync(ss,savg,xpol,jz,nfa,nfb) k=jsync0(j) + lag ccf4=ccf4 + ss(1,k,i+1) + ss(1,k+1,i+1) enddo - ccf4=ccf4 - savg(1,i+1)*2*63/float(jz) + ccf4=ccf4 - savg(i+1)*2*63/float(jz) ccf=ccf4 ipol=1 if(ccf.gt.ccfmax) then @@ -215,7 +215,7 @@ subroutine wb_sync(ss,savg,xpol,jz,nfa,nfb) k=jsync1(j) + lag ccf4=ccf4 + ss(1,k,i+1) + ss(1,k+1,i+1) enddo - ccf4=ccf4 - savg(1,i+1)*2*63/float(jz) + ccf4=ccf4 - savg(i+1)*2*63/float(jz) ccf=ccf4 ipol=1 if(ccf.gt.ccfmax) then @@ -235,7 +235,7 @@ subroutine wb_sync(ss,savg,xpol,jz,nfa,nfb) sync(i)%ipol=ipolbest sync(i)%iflip=flip sync(i)%birdie=.false. - if(ccfmax/(savg(ipolbest,i)/savg_med(ipolbest)).lt.3.0) sync(i)%birdie=.true. + if(ccfmax/(savg(i)/savg_med).lt.3.0) sync(i)%birdie=.true. enddo ! i (frequency bin) call pctile(sync(ia:ib)%ccfmax,ib-ia+1,50,base) diff --git a/q65w/soundin.cpp b/q65w/soundin.cpp index 33a9740b4..42d473fd1 100644 --- a/q65w/soundin.cpp +++ b/q65w/soundin.cpp @@ -16,7 +16,7 @@ extern "C" { double d8[2*60*96000]; //This is "common/datcom/..." in fortran float ss[4*322*NFFT]; - float savg[4*NFFT]; + float savg[NFFT]; double fcenter; int nutc; int idphi; //Phase correction for Y pol'n, degrees