mirror of
https://github.com/saitohirga/WSJT-X.git
synced 2025-11-03 21:40:52 -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
6d2e429144
commit
c539863b6a
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, &
|
||||
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)
|
||||
ihsym=ihsym+1
|
||||
cx=w*cx00 !Apply window for 2nd forward FFT
|
||||
|
||||
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
|
||||
call four2a(cx,nfft3,1,1,1) !Second forward FFT (X)
|
||||
|
||||
n=ihsym
|
||||
do i=1,NFFT
|
||||
n=min(322,ihsym)
|
||||
do i=1,nfft3
|
||||
sx=real(cx(i))**2 + aimag(cx(i))**2
|
||||
ss(1,n,i)=sx ! Pol = 0
|
||||
ss(1,n,i)=sx
|
||||
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
|
||||
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
|
||||
|
||||
@ -1,4 +1,4 @@
|
||||
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
|
||||
@ -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
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user