mirror of
https://github.com/saitohirga/WSJT-X.git
synced 2024-11-25 21:58:38 -05:00
Starting work toward decoding program jt9.
git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@2622 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Joe Taylor
parent
dbbb7a1907
commit
e36f1be536
130
libm65/jt9.f90
Normal file
130
libm65/jt9.f90
Normal file
@ -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
|
||||||
97
libm65/jt9a.f90
Normal file
97
libm65/jt9a.f90
Normal file
@ -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('<m65aFinished>')
|
||||||
|
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
|
||||||
@ -1,55 +1,52 @@
|
|||||||
subroutine symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample,fgreen, &
|
subroutine symspecx(k,nsps,ndiskdat,nb,nbslider,pxdb,s,ihsym, &
|
||||||
iqadjust,iqapply,gainx,gainy,phasex,phasey,rejectx,rejecty, &
|
nzap,slimit,lstrong)
|
||||||
pxdb,pydb,ssz5a,nkhz,ihsym,nzap,slimit,lstrong)
|
|
||||||
|
|
||||||
! k pointer to the most recent new data
|
! 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
|
! ndiskdat 0/1 to indicate if data from disk
|
||||||
! nb 0/1 status of noise blanker
|
! nb 0/1 status of noise blanker (off/on)
|
||||||
! idphi Phase correction for Y channel, degrees
|
! pxdb power (0-60 dB)
|
||||||
! nfsample sample rate (Hz)
|
! s spectrum for waterfall display
|
||||||
! 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
|
|
||||||
! ihsym index number of this half-symbol (1-322)
|
! ihsym index number of this half-symbol (1-322)
|
||||||
! nzap number of samples zero'ed by noise blanker
|
! nzap number of samples zero'ed by noise blanker
|
||||||
|
|
||||||
parameter (NSMAX=60*96000) !Total sample intervals per minute
|
parameter (NMAX=1800*12000) !Total sample intervals per 30 minutes
|
||||||
parameter (NFFT=32768) !Length of FFTs
|
parameter (NSMAX=10000) !Max length of saved spectra
|
||||||
|
parameter (MAXFFT=262144) !Max length of FFTs
|
||||||
|
integer*2 id2
|
||||||
real*8 ts,hsym
|
real*8 ts,hsym
|
||||||
real*8 fcenter
|
real*8 fcenter
|
||||||
common/datcom/dd(4,5760000),ss(4,322,NFFT),savg(4,NFFT),fcenter,nutc,junk(34)
|
real*4 s(NFFT),w(NFFT)
|
||||||
real*4 ssz5a(NFFT),w(NFFT)
|
|
||||||
complex z,zfac
|
complex z,zfac
|
||||||
complex zsumx,zsumy
|
complex zsumx
|
||||||
complex cx(NFFT),cy(NFFT)
|
complex cx(NFFT)
|
||||||
|
complex cx00(NFFT)
|
||||||
complex cx0(0:1023),cx1(0:1023)
|
complex cx0(0:1023),cx1(0:1023)
|
||||||
complex cy0(0:1023),cy1(0:1023)
|
|
||||||
logical*1 lstrong(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
|
save
|
||||||
|
|
||||||
if(k.gt.5751000) go to 999
|
nfft3=nsps
|
||||||
if(k.lt.NFFT) then
|
hsym=nsps/2
|
||||||
|
if(k.gt.NMAX) go to 999
|
||||||
|
if(k.lt.nfft3) then
|
||||||
ihsym=0
|
ihsym=0
|
||||||
go to 999 !Wait for enough samples to start
|
go to 999 !Wait for enough samples to start
|
||||||
endif
|
endif
|
||||||
if(k0.eq.99999999) then
|
if(k0.eq.99999999) then
|
||||||
pi=4.0*atan(1.0)
|
pi=4.0*atan(1.0)
|
||||||
do i=1,NFFT
|
do i=1,nfft3
|
||||||
w(i)=(sin(i*pi/NFFT))**2
|
w(i)=(sin(i*pi/nfft3))**2 !Window for nfft3
|
||||||
enddo
|
enddo
|
||||||
endif
|
endif
|
||||||
|
|
||||||
if(k.lt.k0) then
|
if(k.lt.k0) then
|
||||||
ts=1.d0 - hsym
|
ts=1.d0 - hsym
|
||||||
savg=0.
|
savg=0.
|
||||||
ihsym=0
|
ihsym=0
|
||||||
k1=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
|
endif
|
||||||
k0=k
|
k0=k
|
||||||
|
|
||||||
@ -58,138 +55,64 @@ subroutine symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample,fgreen, &
|
|||||||
peaklimit=sigmas*max(10.0,rms)
|
peaklimit=sigmas*max(10.0,rms)
|
||||||
faclim=3.0
|
faclim=3.0
|
||||||
px=0.
|
px=0.
|
||||||
py=0.
|
|
||||||
|
|
||||||
iqapply0=0
|
|
||||||
iqadjust0=0
|
|
||||||
if(iqadjust.ne.0) iqapply0=0
|
|
||||||
nwindow=2
|
nwindow=2
|
||||||
! nwindow=0 !### No wondowing ###
|
! nwindow=0 !### No windowing ###
|
||||||
nfft2=1024
|
nfft1=1024
|
||||||
kstep=nfft2
|
kstep=nfft1
|
||||||
if(nwindow.ne.0) kstep=nfft2/2
|
if(nwindow.ne.0) kstep=nfft1/2
|
||||||
nblks=(k-k1)/kstep
|
nblks=(k-k1)/kstep
|
||||||
do nblk=1,nblks
|
do nblk=1,nblks
|
||||||
j=k1+1
|
j=k1+1
|
||||||
do i=0,nfft2-1
|
do i=0,nfft1-1
|
||||||
cx0(i)=cmplx(dd(1,j+i),dd(2,j+i))
|
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
|
enddo
|
||||||
call timf2(k,nxpol,nfft2,nwindow,nb,peaklimit,iqadjust0,iqapply0,faclim, &
|
call timf2x(k,nfft1,nwindow,nb,peaklimit,faclim,cx0,cx1, &
|
||||||
cx0,cy0,gainx,gainy,phasex,phasey,cx1,cy1,slimit,lstrong, &
|
slimit,lstrong,px,nzap)
|
||||||
px,py,nzap)
|
|
||||||
|
|
||||||
do i=0,kstep-1
|
do i=0,kstep-1
|
||||||
dd(1,j+i)=real(cx1(i))
|
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
|
enddo
|
||||||
k1=k1+kstep
|
k1=k1+kstep
|
||||||
enddo
|
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
|
ts=ts+hsym
|
||||||
|
ja=ts !Index of first sample
|
||||||
|
jb=ja+nfft3-1 !Last sample
|
||||||
|
|
||||||
i=0
|
i=0
|
||||||
fac=0.0002
|
fac=0.0002
|
||||||
dphi=idphi/57.2957795
|
do j=ja,jb !Copy data into cx
|
||||||
zfac=fac*cmplx(cos(dphi),sin(dphi))
|
|
||||||
do j=ja,jb !Copy data into cx, cy
|
|
||||||
x1=dd(1,j)
|
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
|
i=i+1
|
||||||
cx(i)=fac*cmplx(x1,x2)
|
cx(i)=fac*cmplx(x1,x2)
|
||||||
cy(i)=zfac*cmplx(x3,x4) !NB: cy includes dphi correction
|
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
if(nzap/178.lt.50 .and. (ndiskdat.eq.0 .or. ihsym.lt.280)) then
|
if(nzap/178.lt.50 .and. (ndiskdat.eq.0 .or. ihsym.lt.280)) then
|
||||||
nsum=nblks*kstep - nzap
|
nsum=nblks*kstep - nzap
|
||||||
if(nsum.le.0) nsum=1
|
if(nsum.le.0) nsum=1
|
||||||
rmsx=sqrt(0.5*px/nsum)
|
rmsx=sqrt(0.5*px/nsum)
|
||||||
rmsy=sqrt(0.5*py/nsum)
|
|
||||||
rms=rmsx
|
rms=rmsx
|
||||||
if(nxpol.ne.0) rms=sqrt((px+py)/(4.0*nsum))
|
|
||||||
endif
|
endif
|
||||||
pxdb=0.
|
pxdb=0.
|
||||||
pydb=0.
|
|
||||||
if(rmsx.gt.1.0) pxdb=20.0*log10(rmsx)
|
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(pxdb.gt.60.0) pxdb=60.0
|
||||||
if(pydb.gt.60.0) pydb=60.0
|
|
||||||
|
|
||||||
cx=w*cx !Apply window for 2nd forward FFT
|
cx00=cx
|
||||||
if(nxpol.ne.0) cy=w*cy
|
|
||||||
|
|
||||||
call four2a(cx,NFFT,1,1,1) !Second forward FFT
|
ihsym=ihsym+1
|
||||||
if(iqadjust.eq.0) nadjx=0
|
cx=w*cx00 !Apply window for 2nd forward FFT
|
||||||
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
|
|
||||||
|
|
||||||
if(nxpol.ne.0) then
|
call four2a(cx,nfft3,1,1,1) !Second forward FFT (X)
|
||||||
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
|
n=min(322,ihsym)
|
||||||
ss(3,n,i)=sy ! Pol = 90
|
do i=1,nfft3
|
||||||
savg(3,i)=savg(3,i) + sy
|
sx=real(cx(i))**2 + aimag(cx(i))**2
|
||||||
|
ss(1,n,i)=sx
|
||||||
z=cx(i) - cy(i)
|
savg(1,i)=savg(1,i) + sx
|
||||||
s135=0.5*(real(z)**2 + aimag(z)**2)
|
ssz5a(i)=sx
|
||||||
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
|
|
||||||
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
|
|
||||||
|
|
||||||
999 return
|
999 return
|
||||||
end subroutine symspec
|
end subroutine symspec
|
||||||
|
|||||||
@ -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)
|
slimit,lstrong,px,nzap)
|
||||||
|
|
||||||
! Sequential processing of time-domain I/Q data, using Linrad-like
|
! Sequential processing of time-domain I/Q data, using Linrad-like
|
||||||
! "first FFT" and "first backward FFT".
|
! "first FFT" and "first backward FFT".
|
||||||
|
|
||||||
! cx0 - complex input data
|
! cx0 - complex input data
|
||||||
! nfft - length of FFTs
|
! nfft - length of FFTs
|
||||||
! nwindow - 0 for no window, 2 for sin^2 window
|
! nwindow - 0 for no window, 2 for sin^2 window
|
||||||
! cx1 - output data
|
! cx1 - output data
|
||||||
|
|
||||||
! Non-windowed processing means no overlap, so kstep=nfft.
|
! Non-windowed processing means no overlap, so kstep=nfft.
|
||||||
! Sin^2 window has 50% overlap, kstep=nfft/2.
|
! 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
|
! Frequencies with strong signals are identified and separated. Back
|
||||||
! transforms are done separately for weak and strong signals, so that
|
! transforms are done separately for weak and strong signals, so that
|
||||||
! noise blanking can be applied to the weak-signal portion. Strong and
|
! 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)
|
parameter (MAXSIGS=100)
|
||||||
complex cx0(0:nfft-1),cx1(0:nfft-1)
|
complex cx0(0:nfft-1),cx1(0:nfft-1)
|
||||||
complex cx(0:MAXFFT-1),cxt(0:MAXFFT-1)
|
complex cx(0:MAXFFT-1),cxt(0:MAXFFT-1)
|
||||||
complex cxs(0:MAXFFT-1),covxs(0:MAXNH-1) !Strong X signals
|
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 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 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
|
logical*1 lstrong(0:MAXFFT-1),lprev
|
||||||
integer ia(MAXSIGS),ib(MAXSIGS)
|
integer ia(MAXSIGS),ib(MAXSIGS)
|
||||||
|
complex h,u,v
|
||||||
logical first
|
logical first
|
||||||
data first/.true./
|
data first/.true./
|
||||||
data k0/99999999/
|
data k0/99999999/
|
||||||
save
|
save w,covxs,covxw,s,ntc,ntot,nh,kstep,fac,first,k0
|
||||||
|
|
||||||
if(first) then
|
if(first) then
|
||||||
pi=4.0*atan(1.0)
|
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
|
w(i)=(sin(i*pi/nfft))**2
|
||||||
enddo
|
enddo
|
||||||
s=0.
|
s=0.
|
||||||
|
ntc=0
|
||||||
|
ntot=0
|
||||||
nh=nfft/2
|
nh=nfft/2
|
||||||
nfft2=nfft/4
|
|
||||||
if(ntrperiod.ge.300) nfft2=nfft/32
|
|
||||||
nh2=nfft2/2
|
|
||||||
kstep=nfft
|
kstep=nfft
|
||||||
if(nwindow.eq.2) kstep=nh
|
if(nwindow.eq.2) kstep=nh
|
||||||
fac=1.0/nfft
|
fac=1.0/nfft
|
||||||
@ -59,17 +57,28 @@ subroutine timf2x(k,nfft,ntrperiod,nwindow,nb,peaklimit,faclim,cx0,cx1, &
|
|||||||
|
|
||||||
cx(0:nfft-1)=cx0
|
cx(0:nfft-1)=cx0
|
||||||
if(nwindow.eq.2) cx(0:nfft-1)=w(0:nfft-1)*cx(0:nfft-1)
|
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)
|
cxt(0:nfft-1)=cx(0:nfft-1)
|
||||||
|
|
||||||
! Identify frequencies with strong signals, copy frequency-domain
|
! Identify frequencies with strong signals, copy frequency-domain
|
||||||
! data into array cs (strong) or cw (weak).
|
! 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
|
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
|
enddo
|
||||||
ave=sum(s(0:nfft-1))/nfft
|
|
||||||
lstrong(0:nfft-1)=s(0:nfft-1).gt.10.0*ave
|
|
||||||
|
|
||||||
nsigs=0
|
nsigs=0
|
||||||
lprev=.false.
|
lprev=.false.
|
||||||
@ -105,26 +114,19 @@ subroutine timf2x(k,nfft,ntrperiod,nwindow,nb,peaklimit,faclim,cx0,cx1, &
|
|||||||
cxs(i)=fac*cxt(i)
|
cxs(i)=fac*cxt(i)
|
||||||
cxw(i)=0.
|
cxw(i)=0.
|
||||||
else
|
else
|
||||||
cxs(i)=0.
|
|
||||||
cxw(i)=fac*cxt(i)
|
cxw(i)=fac*cxt(i)
|
||||||
|
cxs(i)=0.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
df=12000.0/nfft
|
call four2a(cxw,nfft,1,-1,1) !Transform weak and strong X
|
||||||
i0=nint(1500.0/df)
|
call four2a(cxs,nfft,1,-1,1) !back to time domain, separately
|
||||||
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
|
|
||||||
|
|
||||||
if(nwindow.eq.2) then
|
if(nwindow.eq.2) then
|
||||||
cxw2(0:nh2-1)=cxw2(0:nh2-1)+covxw(0:nh2-1) !Add prev segment's 2nd half
|
cxw(0:nh-1)=cxw(0:nh-1)+covxw(0:nh-1) !Add previous segment's 2nd half
|
||||||
covxw(0:nh2-1)=cxw2(nh2:nfft2-1) !Save 2nd half
|
covxw(0:nh-1)=cxw(nh:nfft-1) !Save 2nd half
|
||||||
cxs2(0:nh2-1)=cxs2(0:nh2-1)+covxs(0:nh2-1) !Ditto for strong signals
|
cxs(0:nh-1)=cxs(0:nh-1)+covxs(0:nh-1) !Ditto for strong signals
|
||||||
covxs(0:nh2-1)=cxs2(nh2:nfft2-1)
|
covxs(0:nh-1)=cxs(nh:nfft-1)
|
||||||
endif
|
endif
|
||||||
|
|
||||||
! Apply noise blanking to weak data
|
! 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
|
do i=0,kstep-1
|
||||||
peak=abs(cxw(i))
|
peak=abs(cxw(i))
|
||||||
if(peak.gt.peaklimit) then
|
if(peak.gt.peaklimit) then
|
||||||
cxw2(i)=0.
|
cxw(i)=0.
|
||||||
nzap=nzap+1
|
nzap=nzap+1
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
endif
|
endif
|
||||||
|
|
||||||
! Compute power levels from weak data only
|
! Compute power levels from weak data only
|
||||||
px=0.
|
|
||||||
do i=0,kstep-1
|
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
|
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
|
return
|
||||||
end subroutine timf2x
|
end subroutine timf2
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user