mirror of https://github.com/saitohirga/WSJT-X.git
QRA64 decoding basically alive within MAP65.
git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/map65@7487 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Joe Taylor
parent
7e54259ce5
commit
b4145e9b38
|
|
@ -60,6 +60,8 @@ set (FSRCS
|
|||
astro.f90
|
||||
astro0.f90
|
||||
astrosub.f90
|
||||
averms.f90
|
||||
badmsg.f90
|
||||
ccf2.f90
|
||||
ccf65.f90
|
||||
cgen65.f90
|
||||
|
|
@ -79,6 +81,7 @@ set (FSRCS
|
|||
encode65.f90
|
||||
extract.f90
|
||||
fchisq.f90
|
||||
fchisq0.f90
|
||||
fil6521.f90
|
||||
filbig.f90
|
||||
fmtmsg.f90
|
||||
|
|
@ -100,6 +103,7 @@ set (FSRCS
|
|||
iqfix.f90
|
||||
jt65code.f90
|
||||
k2grid.f90
|
||||
lorentzian.f90
|
||||
map65a.f90
|
||||
moon2.f90
|
||||
moondop.f90
|
||||
|
|
@ -107,24 +111,31 @@ set (FSRCS
|
|||
noisegen.f90
|
||||
packjt.f90
|
||||
pctile.f90
|
||||
pctile2.f90
|
||||
pfxdump.f90
|
||||
qra64b.f90
|
||||
qra64c.f90
|
||||
recvpkt.f90
|
||||
rfile3a.f90
|
||||
s3avg.f90
|
||||
sec_midn.f90
|
||||
set.f90
|
||||
setup65.f90
|
||||
shell.f90
|
||||
sleep_msec.f90
|
||||
smo.f90
|
||||
sort.f90
|
||||
spec64.f90
|
||||
sun.f90
|
||||
symspec.f90
|
||||
sync64.f90
|
||||
timer.f90
|
||||
timf2.f90
|
||||
tm2.f90
|
||||
toxyz.f90
|
||||
trimlist.f90
|
||||
twkfreq.f90
|
||||
twkfreq2.f90
|
||||
zplot.f90
|
||||
|
||||
f77_wisdom.f
|
||||
|
|
@ -141,6 +152,15 @@ set (CSRCS
|
|||
tmoonsub.c
|
||||
usleep.c
|
||||
wrapkarn.c
|
||||
|
||||
../../wsjtx/lib/qra/qra64/qra64.c
|
||||
../../wsjtx/lib/qra/qra64/qra64_subs.c
|
||||
../../wsjtx/lib/qra/qracodes/npfwht.c
|
||||
../../wsjtx/lib/qra/qracodes/pdmath.c
|
||||
../../wsjtx/lib/qra/qracodes/qra12_63_64_irr_b.c
|
||||
../../wsjtx/lib/qra/qracodes/qra13_64_64_irr_e.c
|
||||
../../wsjtx/lib/qra/qracodes/qracodes.c
|
||||
../../wsjtx/lib/qra/qracodes/normrnd.c
|
||||
)
|
||||
|
||||
if (WIN32)
|
||||
|
|
|
|||
|
|
@ -0,0 +1,20 @@
|
|||
subroutine averms(x,n,nskip,ave,rms)
|
||||
real x(n)
|
||||
integer ipk(1)
|
||||
|
||||
ns=0
|
||||
s=0.
|
||||
sq=0.
|
||||
ipk=maxloc(x)
|
||||
do i=1,n
|
||||
if(abs(i-ipk(1)).gt.nskip) then
|
||||
s=s + x(i)
|
||||
sq=sq + x(i)**2
|
||||
ns=ns+1
|
||||
endif
|
||||
enddo
|
||||
ave=s/ns
|
||||
rms=sqrt(sq/ns - ave*ave)
|
||||
|
||||
return
|
||||
end subroutine averms
|
||||
|
|
@ -0,0 +1,46 @@
|
|||
subroutine badmsg(irc,dat,nc1,nc2,ng2)
|
||||
|
||||
! Get rid of a few QRA64 false decodes that cannot be correct messages.
|
||||
|
||||
integer dat(12) !Decoded message (as 12 integers)
|
||||
|
||||
ic1=ishft(dat(1),22) + ishft(dat(2),16) + ishft(dat(3),10)+ &
|
||||
ishft(dat(4),4) + iand(ishft(dat(5),-2),15)
|
||||
|
||||
! Test for "......" or "CQ 000"
|
||||
if(ic1.eq.262177560 .or. ic1.eq.262177563) then
|
||||
irc=-1
|
||||
return
|
||||
endif
|
||||
|
||||
ic2=ishft(iand(dat(5),3),26) + ishft(dat(6),20) + &
|
||||
ishft(dat(7),14) + ishft(dat(8),8) + ishft(dat(9),2) + &
|
||||
iand(ishft(dat(10),-4),3)
|
||||
|
||||
ig=ishft(iand(dat(10),15),12) + ishft(dat(11),6) + dat(12)
|
||||
|
||||
! Test for blank, -01 to -30, R-01 to R-30, RO, RRR, 73
|
||||
if(ig.ge.32401 .and. ig.le.32464) return
|
||||
|
||||
if(ig.ge.14220 .and. ig.le.14229) return !-41 to -50
|
||||
if(ig.ge.14040 .and. ig.le.14049) return !-31 to -40
|
||||
|
||||
if(ig.ge.13320 .and. ig.le.13329) return !+00 to +09
|
||||
if(ig.ge.13140 .and. ig.le.13149) return !+10 to +19
|
||||
if(ig.ge.12960 .and. ig.le.12969) return !+20 to +29
|
||||
if(ig.ge.12780 .and. ig.le.12789) return !+30 to +39
|
||||
if(ig.ge.12600 .and. ig.le.12609) return !+40 to +49
|
||||
|
||||
if(ig.ge.12420 .and. ig.le.12429) return !R-41 to R-50
|
||||
if(ig.ge.12240 .and. ig.le.12249) return !R-31 to R-40
|
||||
|
||||
if(ig.ge.11520 .and. ig.le.11529) return !R+00 to R+09
|
||||
if(ig.ge.11340 .and. ig.le.11349) return !R+10 to R+19
|
||||
if(ig.ge.11160 .and. ig.le.11169) return !R+20 to R+29
|
||||
if(ig.ge.10980 .and. ig.le.10989) return !R+30 to R+39
|
||||
if(ig.ge.10800 .and. ig.le.10809) return !R+40 to R+49
|
||||
|
||||
if(ic1.eq.nc1 .and. ic2.eq.nc2 .and. ng2.ne.32401 .and. ig.ne.ng2) irc=-1
|
||||
|
||||
return
|
||||
end subroutine badmsg
|
||||
|
|
@ -0,0 +1,23 @@
|
|||
real function fchisq0(y,npts,a)
|
||||
|
||||
real y(npts),a(4)
|
||||
|
||||
! rewind 51
|
||||
chisq = 0.
|
||||
do i=1,npts
|
||||
x=i
|
||||
z=(x-a(3))/(0.5*a(4))
|
||||
yfit=a(1)
|
||||
if(abs(z).lt.3.0) then
|
||||
d=1.0 + z*z
|
||||
yfit=a(1) + a(2) * (1.0/d - 0.1)
|
||||
endif
|
||||
chisq=chisq + (y(i) - yfit)**2
|
||||
! write(51,3001) i,y(i),yfit,y(i)-yfit
|
||||
!3001 format(i5,3f10.4)
|
||||
enddo
|
||||
fchisq0=chisq
|
||||
|
||||
return
|
||||
end function fchisq0
|
||||
|
||||
|
|
@ -0,0 +1,102 @@
|
|||
subroutine lorentzian(y,npts,a)
|
||||
|
||||
! Input: y(npts); assume x(i)=i, i=1,npts
|
||||
! Output: a(5)
|
||||
! a(1) = baseline
|
||||
! a(2) = amplitude
|
||||
! a(3) = x0
|
||||
! a(4) = width
|
||||
! a(5) = chisqr
|
||||
|
||||
real y(npts)
|
||||
real a(5)
|
||||
real deltaa(4)
|
||||
|
||||
a=0.
|
||||
df=12000.0/8192.0 !df = 1.465 Hz
|
||||
width=0.
|
||||
ipk=0
|
||||
ymax=-1.e30
|
||||
do i=1,npts
|
||||
if(y(i).gt.ymax) then
|
||||
ymax=y(i)
|
||||
ipk=i
|
||||
endif
|
||||
! write(50,3001) i,i*df,y(i)
|
||||
!3001 format(i6,2f12.3)
|
||||
enddo
|
||||
! base=(sum(y(ipk-149:ipk-50)) + sum(y(ipk+51:ipk+150)))/200.0
|
||||
base=(sum(y(1:20)) + sum(y(npts-19:npts)))/40.0
|
||||
stest=ymax - 0.5*(ymax-base)
|
||||
ssum=y(ipk)
|
||||
do i=1,50
|
||||
if(ipk+i.gt.npts) exit
|
||||
if(y(ipk+i).lt.stest) exit
|
||||
ssum=ssum + y(ipk+i)
|
||||
enddo
|
||||
do i=1,50
|
||||
if(ipk-i.lt.1) exit
|
||||
if(y(ipk-i).lt.stest) exit
|
||||
ssum=ssum + y(ipk-i)
|
||||
enddo
|
||||
ww=ssum/y(ipk)
|
||||
width=2
|
||||
t=ww*ww - 5.67
|
||||
if(t.gt.0.0) width=sqrt(t)
|
||||
a(1)=base
|
||||
a(2)=ymax-base
|
||||
a(3)=ipk
|
||||
a(4)=width
|
||||
|
||||
! Now find Lorentzian parameters
|
||||
|
||||
deltaa(1)=0.1
|
||||
deltaa(2)=0.1
|
||||
deltaa(3)=1.0
|
||||
deltaa(4)=1.0
|
||||
nterms=4
|
||||
|
||||
! Start the iteration
|
||||
chisqr=0.
|
||||
chisqr0=1.e6
|
||||
do iter=1,5
|
||||
do j=1,nterms
|
||||
chisq1=fchisq0(y,npts,a)
|
||||
fn=0.
|
||||
delta=deltaa(j)
|
||||
10 a(j)=a(j)+delta
|
||||
chisq2=fchisq0(y,npts,a)
|
||||
if(chisq2.eq.chisq1) go to 10
|
||||
if(chisq2.gt.chisq1) then
|
||||
delta=-delta !Reverse direction
|
||||
a(j)=a(j)+delta
|
||||
tmp=chisq1
|
||||
chisq1=chisq2
|
||||
chisq2=tmp
|
||||
endif
|
||||
20 fn=fn+1.0
|
||||
a(j)=a(j)+delta
|
||||
chisq3=fchisq0(y,npts,a)
|
||||
if(chisq3.lt.chisq2) then
|
||||
chisq1=chisq2
|
||||
chisq2=chisq3
|
||||
go to 20
|
||||
endif
|
||||
|
||||
! Find minimum of parabola defined by last three points
|
||||
delta=delta*(1./(1.+(chisq1-chisq2)/(chisq3-chisq2))+0.5)
|
||||
a(j)=a(j)-delta
|
||||
deltaa(j)=deltaa(j)*fn/3.
|
||||
! write(*,4000) iter,j,a,chisq2
|
||||
!4000 format(i1,i2,4f10.4,f11.3)
|
||||
enddo
|
||||
chisqr=fchisq0(y,npts,a)
|
||||
! write(*,4000) 0,0,a,chisqr
|
||||
if(chisqr/chisqr0.gt.0.99) exit
|
||||
chisqr0=chisqr
|
||||
enddo
|
||||
a(5)=chisqr
|
||||
|
||||
return
|
||||
end subroutine lorentzian
|
||||
|
||||
|
|
@ -33,7 +33,7 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, &
|
|||
data nfile/0/,nutc0/-999/,nid/0/,ip000/1/,ip001/1/,mousefqso0/-999/
|
||||
save
|
||||
|
||||
bqra64=nfast.ge.100
|
||||
bqra64=nfast.ge.100
|
||||
nfast=mod(nfast,100)
|
||||
mcall3a=mcall3b
|
||||
mousefqso0=mousefqso
|
||||
|
|
@ -218,7 +218,8 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, &
|
|||
a,dt,pol,nkv,nhist,nsum,nsave,qual,decoded)
|
||||
call timer('decode1a',1)
|
||||
if(nqd.eq.2) then
|
||||
call qra64b(nutc,ikhz)
|
||||
call qra64b(nutc,nqd,ikhz,mousedf,ntol,xpol,mycall, &
|
||||
hiscall,hisgrid)
|
||||
cycle
|
||||
endif
|
||||
|
||||
|
|
@ -348,6 +349,7 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, &
|
|||
1013 format('<QuickDecodeDone>',2i4)
|
||||
flush(6)
|
||||
endif
|
||||
if(nqd.eq.2) exit !### TESTING: do only QRA64
|
||||
if(nagain.eq.1 .and. nqd.eq.1) go to 999
|
||||
enddo
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,22 @@
|
|||
subroutine pctile2(x,npts,npct,xpct)
|
||||
|
||||
parameter (NMAX=100000)
|
||||
real*4 x(npts)
|
||||
real*4 tmp(NMAX)
|
||||
|
||||
if(npts.le.0) then
|
||||
xpct=1.0
|
||||
go to 900
|
||||
endif
|
||||
if(npts.gt.NMAX) stop
|
||||
|
||||
tmp(1:npts)=x
|
||||
call shell(npts,tmp)
|
||||
j=nint(npts*0.01*npct)
|
||||
if(j.lt.1) j=1
|
||||
if(j.gt.npts) j=npts
|
||||
xpct=tmp(j)
|
||||
|
||||
900 continue
|
||||
return
|
||||
end subroutine pctile2
|
||||
|
|
@ -1,14 +1,14 @@
|
|||
subroutine qra64b(nutc,ikhz)
|
||||
subroutine qra64b(nutc,nqd,ikhz,idf,ntol,xpol,mycall_12,hiscall_12, &
|
||||
hisgrid_6)
|
||||
|
||||
parameter (NFFT1=5376000) !56*96000
|
||||
parameter (NFFT2=336000) !56*6000 (downsampled by 1/16)
|
||||
complex cx(0:NFFT2-1),cy(0:NFFT2-1)
|
||||
complex ca(NFFT1),cb(NFFT1) !FFTs of raw x,y data
|
||||
logical xpol
|
||||
character*12 mycall_12,hiscall_12
|
||||
character*6 hisgrid_6
|
||||
common/cacb/ca,cb
|
||||
|
||||
!###
|
||||
if(nutc.ne.-999) return
|
||||
!###
|
||||
|
||||
df=96000.0/NFFT1
|
||||
k0=(ikhz-75.7)*1000.0/df
|
||||
|
|
@ -25,7 +25,9 @@ subroutine qra64b(nutc,ikhz)
|
|||
call four2a(cx,NFFT2,1,-1,1)
|
||||
call four2a(cy,NFFT2,1,-1,1)
|
||||
|
||||
write(67) nutc,cx,cy
|
||||
! write(67) nutc,cx,cy
|
||||
call qra64c(cx,cy,nutc,nqd,ikhz,idf,ntol,xplo,mycall_12, &
|
||||
hiscall_12,hisgrid_6)
|
||||
|
||||
return
|
||||
end subroutine qra64b
|
||||
|
|
|
|||
|
|
@ -0,0 +1,181 @@
|
|||
subroutine qra64c(cx,cy,nutc,nqd,ikhz,nfqso,ntol,xpol,mycall_12, &
|
||||
hiscall_12,hisgrid_6)
|
||||
|
||||
use packjt
|
||||
parameter (NFFT2=336000) !56*6000 (downsampled by 1/16)
|
||||
parameter (NMAX=60*12000,LN=1152*63)
|
||||
|
||||
! Required input data:
|
||||
! nutc,cx,cy,nf1,nf2,nfqso,ntol,mode64,minsync,ndepth,emedelay,
|
||||
! mycall_12,hiscall_12,hisgrid_6
|
||||
|
||||
character decoded*22
|
||||
character*12 mycall_12,hiscall_12
|
||||
character*6 mycall,hiscall,hisgrid_6
|
||||
character*4 hisgrid
|
||||
character*1 cp
|
||||
logical xpol,ltext
|
||||
complex cx(0:NFFT2-1),cy(0:NFFT2-1)
|
||||
complex c00(0:720000) !Complex spectrum of dd()
|
||||
complex c0(0:720000) !Complex data for dd()
|
||||
real a(3)
|
||||
real s3(LN) !Symbol spectra
|
||||
real s3a(LN) !Symbol spectra
|
||||
integer dat4(12) !Decoded message (as 12 integers)
|
||||
integer dat4x(12)
|
||||
integer nap(0:11)
|
||||
data nap/0,2,3,2,3,4,2,3,6,4,6,6/
|
||||
data nc1z/-1/,nc2z/-1/,ng2z/-1/,maxaptypez/-1/
|
||||
save
|
||||
|
||||
! For now:
|
||||
nf1=200
|
||||
nf2=2200
|
||||
! nfqso=808
|
||||
! ntol=1000
|
||||
mode64=1
|
||||
minsync=-1
|
||||
ndepth=3
|
||||
emedelay=2.5
|
||||
|
||||
irc=-1
|
||||
decoded=' '
|
||||
nft=99
|
||||
mycall=mycall_12(1:6)
|
||||
hiscall=hiscall_12(1:6)
|
||||
hisgrid=hisgrid_6(1:4)
|
||||
call packcall(mycall,nc1,ltext)
|
||||
call packcall(hiscall,nc2,ltext)
|
||||
call packgrid(hisgrid,ng2,ltext)
|
||||
nSubmode=0
|
||||
if(mode64.eq.2) nSubmode=1
|
||||
if(mode64.eq.4) nSubmode=2
|
||||
if(mode64.eq.8) nSubmode=3
|
||||
if(mode64.eq.16) nSubmode=4
|
||||
b90=1.0
|
||||
nFadingModel=1
|
||||
maxaptype=4
|
||||
if(iand(ndepth,64).ne.0) maxaptype=5
|
||||
if(nc1.ne.nc1z .or. nc2.ne.nc2z .or. ng2.ne.ng2z .or. &
|
||||
maxaptype.ne.maxaptypez) then
|
||||
do naptype=0,maxaptype
|
||||
if(naptype.eq.2 .and. maxaptype.eq.4) cycle
|
||||
call qra64_dec(s3,nc1,nc2,ng2,naptype,1,nSubmode,b90, &
|
||||
nFadingModel,dat4,snr2,irc)
|
||||
enddo
|
||||
nc1z=nc1
|
||||
nc2z=nc2
|
||||
ng2z=ng2
|
||||
maxaptypez=maxaptype
|
||||
endif
|
||||
naptype=maxaptype
|
||||
npts2=NFFT2
|
||||
|
||||
!1 read(67,end=999) nutc,cx,cy
|
||||
c00(0:NFFT2-1)=conjg(cy)
|
||||
call sync64(c00,nf1,nf2,nfqso,ntol,mode64,emedelay,dtx,f0,jpk0,sync, &
|
||||
sync2,width)
|
||||
|
||||
nfreq=nint(f0)
|
||||
if(mode64.eq.1 .and. minsync.ge.0 .and. (sync-7.0).lt.minsync) go to 900
|
||||
a=0.
|
||||
a(1)=-f0
|
||||
call twkfreq2(c00,c0,npts2,6000.0,a)
|
||||
|
||||
irc=-99
|
||||
s3lim=20.
|
||||
itz=11
|
||||
if(mode64.eq.4) itz=9
|
||||
if(mode64.eq.2) itz=7
|
||||
if(mode64.eq.1) itz=5
|
||||
|
||||
LL=64*(mode64+2)
|
||||
NN=63
|
||||
napmin=99
|
||||
do itry0=1,5
|
||||
idt=itry0/2
|
||||
if(mod(itry0,2).eq.0) idt=-idt
|
||||
jpk=jpk0 + 750*idt
|
||||
call spec64(c0,npts2,mode64,jpk,s3a,LL,NN)
|
||||
call pctile2(s3a,LL*NN,40,base)
|
||||
s3a=s3a/base
|
||||
where(s3a(1:LL*NN)>s3lim) s3a(1:LL*NN)=s3lim
|
||||
do iter=itz,0,-2
|
||||
b90=1.728**iter
|
||||
if(b90.gt.230.0) cycle
|
||||
if(b90.lt.0.15*width) exit
|
||||
s3(1:LL*NN)=s3a(1:LL*NN)
|
||||
call qra64_dec(s3,nc1,nc2,ng2,naptype,0,nSubmode,b90, &
|
||||
nFadingModel,dat4,snr2,irc)
|
||||
if(irc.eq.0) go to 10
|
||||
if(irc.gt.0) call badmsg(irc,dat4,nc1,nc2,ng2)
|
||||
iirc=max(0,min(irc,11))
|
||||
if(irc.gt.0 .and. nap(iirc).lt.napmin) then
|
||||
dat4x=dat4
|
||||
b90x=b90
|
||||
snr2x=snr2
|
||||
napmin=nap(iirc)
|
||||
irckeep=irc
|
||||
dtxkeep=jpk/6000.0 - 1.0
|
||||
itry0keep=itry0
|
||||
iterkeep=iter
|
||||
endif
|
||||
enddo
|
||||
if(irc.eq.0) exit
|
||||
enddo
|
||||
|
||||
if(napmin.ne.99) then
|
||||
dat4=dat4x
|
||||
b90=b90x
|
||||
snr2=snr2x
|
||||
irc=irckeep
|
||||
dtx=dtxkeep
|
||||
itry0=itry0keep
|
||||
iter=iterkeep
|
||||
endif
|
||||
10 decoded=' '
|
||||
|
||||
if(irc.ge.0) then
|
||||
call unpackmsg(dat4,decoded) !Unpack the user message
|
||||
call fmtmsg(decoded,iz)
|
||||
if(index(decoded,"000AAA ").ge.1) then
|
||||
! Suppress a certain type of garbage decode.
|
||||
decoded=' '
|
||||
irc=-1
|
||||
endif
|
||||
nft=100 + irc
|
||||
nsnr=nint(snr2)
|
||||
else
|
||||
snr2=0.
|
||||
endif
|
||||
|
||||
900 if(irc.lt.0) then
|
||||
sy=max(1.0,sync)
|
||||
if(nSubmode.eq.0) nsnr=nint(10.0*log10(sy)-35.0) !A
|
||||
if(nSubmode.eq.1) nsnr=nint(10.0*log10(sy)-34.0) !B
|
||||
if(nSubmode.eq.2) nsnr=nint(10.0*log10(sy)-29.0) !C
|
||||
if(nSubmode.eq.3) nsnr=nint(10.0*log10(sy)-29.0) !D
|
||||
if(nSubmode.eq.4) nsnr=nint(10.0*log10(sy)-24.0) !E
|
||||
endif
|
||||
|
||||
! write(*,1011) nutc/100,nsnr,dtx,nfreq,decoded
|
||||
!1011 format(i4.4,i4,f5.1,i5,1x,2x,1x,a22)
|
||||
|
||||
nkhz=108
|
||||
npol=0
|
||||
cp='H'
|
||||
nsync=sync
|
||||
ntxpol=0
|
||||
if(irc.ge.0) then
|
||||
write(*,1010) nkHz,nfreq,npol,nutc/100,dtx,nsnr,decoded,irc,ntxpol,cp
|
||||
!1010 format('!',i3,i5,i4,i7.6,f5.1,i4,2x,a22,i2,i5,i5,1x,a1)
|
||||
!1010 format(i3,i5,i4,i5.4,f5.1,i5,2x,a22,i2,i5,1x,a1)
|
||||
1010 format('!',i3,i5,i4,i6.4,f5.1,i5,2x,a22,i2,i5,1x,a1)
|
||||
else
|
||||
write(*,1010) nkHz,nfreq,npol,nutc/100,dtx,nsync
|
||||
endif
|
||||
|
||||
! goto 1
|
||||
|
||||
999 return
|
||||
end subroutine qra64c
|
||||
|
|
@ -0,0 +1,27 @@
|
|||
subroutine shell(n,a)
|
||||
integer n
|
||||
real a(n)
|
||||
integer i,j,inc
|
||||
real v
|
||||
|
||||
inc=1
|
||||
1 inc=3*inc+1
|
||||
if(inc.le.n) go to 1
|
||||
2 inc=inc/3
|
||||
|
||||
do i=inc+1,n
|
||||
v=a(i)
|
||||
j=i
|
||||
3 if(a(j-inc).gt.v) then
|
||||
a(j)=a(j-inc)
|
||||
j=j-inc
|
||||
if(j.le.inc) go to 4
|
||||
go to 3
|
||||
endif
|
||||
4 a(j)=v
|
||||
enddo
|
||||
|
||||
if(inc.gt.1) go to 2
|
||||
|
||||
return
|
||||
end subroutine shell
|
||||
|
|
@ -0,0 +1,19 @@
|
|||
subroutine smo(x,npts,y,nadd)
|
||||
|
||||
real x(npts)
|
||||
real y(npts)
|
||||
|
||||
nh=nadd/2
|
||||
do i=1+nh,npts-nh
|
||||
sum=0.
|
||||
do j=-nh,nh
|
||||
sum=sum + x(i+j)
|
||||
enddo
|
||||
y(i)=sum
|
||||
enddo
|
||||
x=y
|
||||
x(:nh)=0.
|
||||
x(npts-nh+1:)=0.
|
||||
|
||||
return
|
||||
end subroutine smo
|
||||
|
|
@ -0,0 +1,42 @@
|
|||
subroutine spec64(c0,npts2,mode64,jpk,s3,LL,NN)
|
||||
|
||||
parameter (NSPS=3456) !Samples per symbol at 6000 Hz
|
||||
complex c0(0:360000) !Complex spectrum of dd()
|
||||
complex cs(0:NSPS-1) !Complex symbol spectrum
|
||||
real s3(LL,NN) !Synchronized symbol spectra
|
||||
real xbase0(LL),xbase(LL)
|
||||
|
||||
nfft=nsps
|
||||
fac=1.0/nfft
|
||||
do j=1,NN
|
||||
jj=j+7 !Skip first Costas array
|
||||
if(j.ge.33) jj=j+14 !Skip middle Costas array
|
||||
ja=jpk + (jj-1)*nfft
|
||||
jb=ja+nfft-1
|
||||
cs(0:nfft-1)=fac*c0(ja:jb)
|
||||
call four2a(cs,nfft,1,-1,1)
|
||||
do ii=1,LL
|
||||
i=ii-65
|
||||
if(i.lt.0) i=i+nfft
|
||||
s3(ii,j)=real(cs(i))**2 + aimag(cs(i))**2
|
||||
enddo
|
||||
enddo
|
||||
|
||||
df=6000.0/nfft
|
||||
do i=1,LL
|
||||
call pctile2(s3(i,1:NN),NN,45,xbase0(i)) !Get baseline for passband shape
|
||||
enddo
|
||||
|
||||
nh=25
|
||||
xbase(1:nh-1)=sum(xbase0(1:nh-1))/(nh-1.0)
|
||||
xbase(LL-nh+1:LL)=sum(xbase0(LL-nh+1:LL))/(nh-1.0)
|
||||
do i=nh,LL-nh
|
||||
xbase(i)=sum(xbase0(i-nh+1:i+nh))/(2*nh+1) !Smoothed passband shape
|
||||
enddo
|
||||
|
||||
do i=1,LL
|
||||
s3(i,1:NN)=s3(i,1:NN)/(xbase(i)+0.001) !Apply frequency equalization
|
||||
enddo
|
||||
|
||||
return
|
||||
end subroutine spec64
|
||||
|
|
@ -0,0 +1,154 @@
|
|||
subroutine sync64(c0,nf1,nf2,nfqso,ntol,mode64,emedelay,dtx,f0,jpk,sync, &
|
||||
sync2,width)
|
||||
|
||||
! use timer_module, only: timer
|
||||
|
||||
parameter (NMAX=60*12000) !Max size of raw data at 12000 Hz
|
||||
parameter (NSPS=3456) !Samples per symbol at 6000 Hz
|
||||
parameter (NSPC=7*NSPS) !Samples per Costas array
|
||||
real s1(0:NSPC-1) !Power spectrum of Costas 1
|
||||
real s2(0:NSPC-1) !Power spectrum of Costas 2
|
||||
real s3(0:NSPC-1) !Power spectrum of Costas 3
|
||||
real s0(0:NSPC-1) !Sum of s1+s2+s3
|
||||
real s0a(0:NSPC-1) !Best synchromized spectrum (saved)
|
||||
real s0b(0:NSPC-1) !tmp
|
||||
real a(5)
|
||||
integer icos7(0:6) !Costas 7x7 tones
|
||||
integer ipk0(1)
|
||||
complex cc(0:NSPC-1) !Costas waveform
|
||||
complex c0(0:720000) !Complex spectrum of dd()
|
||||
complex c1(0:NSPC-1) !Complex spectrum of Costas 1
|
||||
complex c2(0:NSPC-1) !Complex spectrum of Costas 2
|
||||
complex c3(0:NSPC-1) !Complex spectrum of Costas 3
|
||||
data icos7/2,5,6,0,4,1,3/ !Costas 7x7 tone pattern
|
||||
data mode64z/-1/
|
||||
save
|
||||
|
||||
if(mode64.ne.mode64z) then
|
||||
twopi=8.0*atan(1.0)
|
||||
dfgen=mode64*12000.0/6912.0
|
||||
k=-1
|
||||
phi=0.
|
||||
do j=0,6 !Compute complex Costas waveform
|
||||
dphi=twopi*10.0*icos7(j)*dfgen/6000.0
|
||||
do i=1,NSPS
|
||||
phi=phi + dphi
|
||||
if(phi.gt.twopi) phi=phi-twopi
|
||||
k=k+1
|
||||
cc(k)=cmplx(cos(phi),sin(phi))
|
||||
enddo
|
||||
enddo
|
||||
mode64z=mode64
|
||||
endif
|
||||
|
||||
nfft3=NSPC
|
||||
nh3=nfft3/2
|
||||
df3=6000.0/nfft3
|
||||
|
||||
fa=max(nf1,nfqso-ntol)
|
||||
fb=min(nf2,nfqso+ntol)
|
||||
iaa=max(0,nint(fa/df3))
|
||||
ibb=min(NSPC-1,nint(fb/df3))
|
||||
|
||||
maxtol=max(ntol,500)
|
||||
fa=max(nf1,nfqso-maxtol)
|
||||
fb=min(nf2,nfqso+maxtol)
|
||||
ia=max(0,nint(fa/df3))
|
||||
ib=min(NSPC-1,nint(fb/df3))
|
||||
id=0.1*(ib-ia)
|
||||
iz=ib-ia+1
|
||||
sync=-1.e30
|
||||
smaxall=0.
|
||||
jpk=0
|
||||
ja=0
|
||||
jb=(5.0+emedelay)*6000
|
||||
jstep=100
|
||||
ipk=0
|
||||
kpk=0
|
||||
nadd=10*mode64
|
||||
if(mod(nadd,2).eq.0) nadd=nadd+1 !Make nadd odd
|
||||
nskip=max(49,nadd)
|
||||
|
||||
do j1=ja,jb,jstep
|
||||
call timer('sync64_1',0)
|
||||
j2=j1 + 39*NSPS
|
||||
j3=j1 + 77*NSPS
|
||||
c1=1.e-4*c0(j1:j1+NSPC-1) * conjg(cc)
|
||||
c2=1.e-4*c0(j2:j2+NSPC-1) * conjg(cc)
|
||||
c3=1.e-4*c0(j3:j3+NSPC-1) * conjg(cc)
|
||||
call four2a(c1,nfft3,1,-1,1)
|
||||
call four2a(c2,nfft3,1,-1,1)
|
||||
call four2a(c3,nfft3,1,-1,1)
|
||||
s1=0.
|
||||
s2=0.
|
||||
s3=0.
|
||||
s0b=0.
|
||||
do i=ia,ib
|
||||
freq=i*df3
|
||||
s1(i)=real(c1(i))**2 + aimag(c1(i))**2
|
||||
s2(i)=real(c2(i))**2 + aimag(c2(i))**2
|
||||
s3(i)=real(c3(i))**2 + aimag(c3(i))**2
|
||||
enddo
|
||||
call timer('sync64_1',1)
|
||||
|
||||
call timer('sync64_2',0)
|
||||
s0(ia:ib)=s1(ia:ib) + s2(ia:ib) + s3(ia:ib)
|
||||
s0(:ia-1)=0.
|
||||
s0(ib+1:)=0.
|
||||
if(nadd.ge.3) then
|
||||
do ii=1,3
|
||||
s0b(ia:ib)=s0(ia:ib)
|
||||
call smo(s0b(ia:ib),iz,s0(ia:ib),nadd)
|
||||
enddo
|
||||
endif
|
||||
call averms(s0(ia+id:ib-id),iz-2*id,nskip,ave,rms)
|
||||
s=(maxval(s0(iaa:ibb))-ave)/rms
|
||||
ipk0=maxloc(s0(iaa:ibb))
|
||||
ip=ipk0(1) + iaa - 1
|
||||
if(s.gt.sync) then
|
||||
jpk=j1
|
||||
s0a=(s0-ave)/rms
|
||||
sync=s
|
||||
dtx=jpk/6000.0 - 1.0
|
||||
ipk=ip
|
||||
f0=ip*df3
|
||||
endif
|
||||
call timer('sync64_2',1)
|
||||
enddo
|
||||
|
||||
s0a=s0a+2.0
|
||||
write(17) ia,ib,s0a(ia:ib) !Save data for red curve
|
||||
close(17)
|
||||
|
||||
nskip=50
|
||||
call lorentzian(s0a(ia+nskip:ib-nskip),iz-2*nskip,a)
|
||||
f0a=(a(3)+ia+49)*df3
|
||||
w1=df3*a(4)
|
||||
w2=2*nadd*df3
|
||||
width=w1
|
||||
if(w1.gt.1.2*w2) width=sqrt(w1**2 - w2**2)
|
||||
|
||||
sq=0.
|
||||
do i=1,20
|
||||
j=ia+nskip+1
|
||||
k=ib-nskip-21+i
|
||||
sq=sq + (s0a(j)-a(1))**2 + (s0a(k)-a(1))**2
|
||||
enddo
|
||||
rms2=sqrt(sq/40.0)
|
||||
sync2=10.0*log10(a(2)/rms2)
|
||||
|
||||
! do i=1,iz-2*nskip
|
||||
! x=i
|
||||
! z=(x-a(3))/(0.5*a(4))
|
||||
! yfit=a(1)
|
||||
! if(abs(z).lt.3.0) then
|
||||
! d=1.0 + z*z
|
||||
! yfit=a(1) + a(2)*(1.0/d - 0.1)
|
||||
! endif
|
||||
! j=i+ia+49
|
||||
! write(76,1110) j*df3,s0a(j),yfit
|
||||
!1110 format(3f10.3)
|
||||
! enddo
|
||||
|
||||
return
|
||||
end subroutine sync64
|
||||
|
|
@ -0,0 +1,26 @@
|
|||
subroutine twkfreq2(c3,c4,npts,fsample,a)
|
||||
|
||||
complex c3(npts)
|
||||
complex c4(npts)
|
||||
complex w,wstep
|
||||
real a(3)
|
||||
data twopi/6.283185307/
|
||||
|
||||
! Mix the complex signal
|
||||
w=1.0
|
||||
wstep=1.0
|
||||
x0=0.5*(npts+1)
|
||||
s=2.0/npts
|
||||
do i=1,npts
|
||||
x=s*(i-x0)
|
||||
p2=1.5*x*x - 0.5
|
||||
! p3=2.5*(x**3) - 1.5*x
|
||||
! p4=4.375*(x**4) - 3.75*(x**2) + 0.375
|
||||
dphi=(a(1) + x*a(2) + p2*a(3)) * (twopi/fsample)
|
||||
wstep=cmplx(cos(dphi),sin(dphi))
|
||||
w=w*wstep
|
||||
c4(i)=w*c3(i)
|
||||
enddo
|
||||
|
||||
return
|
||||
end subroutine twkfreq2
|
||||
Loading…
Reference in New Issue