Convert savg(4,NFFT) to one-dimensional savg(NFFT).

This commit is contained in:
Joe Taylor 2022-12-11 17:10:40 -05:00
parent fed64b6d9f
commit 6ccd2a290b
10 changed files with 24 additions and 358 deletions

View File

@ -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

View File

@ -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

View File

@ -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

View File

@ -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

View File

@ -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, &

View File

@ -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)

View File

@ -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)

View File

@ -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

View File

@ -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)

View File

@ -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