First good compile with map65a running inside map65.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/map65@334 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Joe Taylor 2007-01-11 21:25:52 +00:00
parent bc06bc211e
commit c43bfde2ed
94 changed files with 6836 additions and 6104 deletions

124
decode1a.f Normal file
View File

@ -0,0 +1,124 @@
subroutine decode1a(id,newdat,nfilt,freq,nflip,ipol,sync2,a,dt,
+ pol,nkv,nhist,qual,decoded)
C Apply AFC corrections to a candidate JT65 signal, and then try
C to decode it.
parameter (NFFT1=77760,NFFT2=2430)
parameter (NMAX=60*96000) !Samples per 60 s
integer*2 id(4,NMAX) !46 MB: raw data from Linrad timf2
complex c2x(NMAX/4), c2y(NMAX/4) !After 1/4 filter and downsample
complex c3x(NMAX/16),c3y(NMAX/16) !After 1/16 filter and downsample
complex c4x(NMAX/64),c4y(NMAX/64) !After 1/64 filter and downsample
complex cx(NMAX/64), cy(NMAX/64) !Data at 1378.125 samples/s
complex c5x(NMAX/256),c5y(NMAX/256)
complex c5a(256), c5b(256)
real s2(256,126)
real a(5)
real*8 samratio
integer resample
logical first
character decoded*22
data first/.true./,jjjmin/1000/,jjjmax/-1000/
save
C Mix sync tone to baseband, low-pass filter, and decimate by 64
dt00=dt
C If freq=125.0 kHz, f0=48000 Hz.
f0=1000*(freq-77.0) !Freq of sync tone (0-96000 Hz)
if(nfilt.eq.1) then
call filbig(id,NMAX,f0,newdat,cx,cy,n5)
joff=0
else
call fil659(id,NMAX,f0,c2x,c2y,n2) !Pass 1: mix and filter both pol'ns
call fil658(c2x,n2,c3x,n3) !Pass 2
call fil658(c2y,n2,c3y,n3)
call fil658(c3x,n3,c4x,n4) !Pass 3
call fil658(c3y,n3,c4y,n4)
joff=-8
C Resample from 96000/64 = 1500 Hz to 1378.125 Hz
C Converter type: 0=Best quality sinc (band limited), BW=97%
C 1=medium quality sinc, BW=90%
C 2=fastest sinc, BW=80%
C 3=stepwise (very fast)
C 4=linear (very fast)
nconv_type=2 !### test! ###
nchans=2
samratio=1378.125d0/1500.d0
i1=resample(c4x,n4,nconv_type,nchans,samratio,cx,n5)
i2=resample(c4y,n4,nconv_type,nchans,samratio,cy,n5)
endif
sqa=0.
sqb=0.
do i=1,n5
sqa=sqa + real(cx(i))**2 + aimag(cx(i))**2
sqb=sqb + real(cy(i))**2 + aimag(cy(i))**2
enddo
sqa=sqa/n5
sqb=sqb/n5
C Find best DF, f1, f2, DT, and pol
! a(5)=dt00
! fsample=1378.125
! i0=nint((a(5)+0.5)*fsample)
! if(i0.lt.1) i0=1
! nz=n5+1-i0
! call afc65b(cx(i0),cy(i0),nz,fsample,nflip,ipol,a,dt,
! + ccfbest,dtbest)
call fil6521(cx,n5,c5x,n6)
call fil6521(cy,n5,c5y,n6)
fsample=1378.125/4.
a(5)=dt00
i0=nint((a(5)+0.5)*fsample) - 2
if(i0.lt.1) i0=1
nz=n6+1-i0
call afc65b(c5x(i0),c5y(i0),nz,fsample,nflip,ipol,a,dt,
+ ccfbest,dtbest)
pol=a(4)/57.2957795
aa=cos(pol)
bb=sin(pol)
sq0=aa*aa*sqa + bb*bb*sqb
sync2=3.7*ccfbest/sq0
C Apply AFC corrections to the time-domain signal
call twkfreq(cx,cy,n5,a)
C Compute spectrum at best polarization for each half symbol.
C Adding or subtracting a small number (e.g., 5) to j may make it decode.
nsym=126
nfft=256
j=(dt00+dtbest+2.685)*1378.125 + joff
if(j.lt.0) j=0
j0=j
do k=1,nsym
do i=1,nfft
j=j+1
c5a(i)=aa*cx(j) + bb*cy(j)
enddo
call four2a(c5a,nfft,1,1,1)
do i=1,nfft
j=j+1
c5b(i)=aa*cx(j) + bb*cy(j)
enddo
call four2a(c5b,nfft,1,1,1)
do i=1,256
s2(i,k)=real(c5a(i))**2 + aimag(c5a(i))**2 +
+ real(c5b(i))**2 + aimag(c5b(i))**2
enddo
enddo
flip=nflip
call decode65b(s2,flip,nkv,nhist,qual,decoded)
dt=dt00 + dtbest
return
end

59
decode65b.f Normal file
View File

@ -0,0 +1,59 @@
subroutine decode65b(s2,flip,nkv,nhist,qual,decoded)
real s2(256,126)
real s3(64,63)
logical first
character decoded*22,deepmsg*22
character mycall*12,hiscall*12,hisgrid*6
! include 'avecom.h'
include 'prcom.h'
data first/.true./
save
if(first) call setup65
first=.false.
call setup65
do j=1,63
k=mdat(j) !Points to data symbol
if(flip.lt.0.0) k=mdat2(j)
do i=1,64
s3(i,j)=s2(i+2,k) !### Check the "i+2" ###
enddo
enddo
mode65=2
nadd=mode65
call extract(s3,nadd,ncount,nhist,decoded) !Extract the message
C Suppress "birdie messages":
if(decoded(1:7).eq.'000AAA ') ncount=-1
if(decoded(1:7).eq.'0L6MWK ') ncount=-1
nkv=1
if(ncount.lt.0) then
nkv=0
decoded=' '
endif
qual=0.
if(nkv.eq.0) then
mycall='K1JT'
hiscall='W1ABC'
hisgrid='EM79'
neme=0
nsked=0
ndepth=5
if(ndepth.ge.1) call deep65(s3,mode65,neme,
+ nsked,flip,mycall,hiscall,hisgrid,deepmsg,qual)
C Save symbol spectra for possible decoding of average.
! do j=1,63
! k=mdat(j)
! if(flip.lt.0.0) k=mdat2(j)
! call move(s2(8,k),ppsave(1,j,nsave),64)
! enddo
endif
if(nkv.eq.0 .and. qual.ge.1.0) decoded=deepmsg
return
end

View File

@ -5,18 +5,17 @@
real s3(64,63)
character callsign*12,grid*4,message*22,hisgrid*6,c*1,ceme*3
character*12 mycall,hiscall
character mycall0*12,hiscall0*12,hisgrid0*6
character*22 decoded
character*22 testmsg(2*MAXCALLS + 2 + MAXRPT)
character*15 callgrid(MAXCALLS)
character*180 line
character*4 rpt(MAXRPT)
integer ncode(63,2*MAXCALLS + 2 + MAXRPT)
integer nflip(2*MAXCALLS + 2 + MAXRPT)
integer istat23(13)
real pp(2*MAXCALLS + 2 + MAXRPT)
common/tmp9/ mrs(63),mrs2(63)
#ifdef Win32
C This prevents some optimizations that break this subroutine.
volatile p1,p2,bias
#endif
common/mrscom/ mrs(63),mrs2(63)
data neme0/-99/
data rpt/'-01','-02','-03','-04','-05',
@ -32,7 +31,13 @@ C This prevents some optimizations that break this subroutine.
+ 'R-21','R-22','R-23','R-24','R-25',
+ 'R-26','R-27','R-28','R-29','R-30',
+ 'RO','RRR','73'/
save
! call fstatqqq(23,istat23,ierr) !@@@
! modified=istat23(10) !@@@
modified=0 !@@@
if(mycall.eq.mycall0 .and. hiscall.eq.hiscall0 .and.
+ hisgrid.eq.hisgrid0 .and. modified.eq.modified0) go to 30
rewind 23
k=0
icall=0
@ -77,7 +82,7 @@ C This prevents some optimizations that break this subroutine.
mz=1
if(n.eq.1 .and. j3.lt.1 .and. j4.lt.1 .and.
+ flip.gt.0.0 .and. callsign(1:6).ne.' ') mz=MAXRPT+1
+ callsign(1:6).ne.' ') mz=MAXRPT+1
C Test for messages with MyCall + HisCall + report
do m=1,mz
if(m.gt.1) grid=rpt(m-1)
@ -87,12 +92,14 @@ C Test for messages with MyCall + HisCall + report
k=k+1
testmsg(k)=message
call encode65(message,ncode(1,k))
C Insert CQ message unless sync=OOO (flip=-1).
nflip(k)=flip
C Insert CQ message
if(m.eq.1 .and. flip.gt.0.0) then
message='CQ '//callgrid(icall)
k=k+1
testmsg(k)=message
call encode65(message,ncode(1,k))
nflip(k)=flip
endif
enddo
if(nsked.eq.1) go to 20
@ -101,28 +108,33 @@ C Insert CQ message unless sync=OOO (flip=-1).
20 ntot=k
neme0=neme
30 mycall0=mycall
hiscall0=hiscall
hisgrid0=hisgrid
modified0=modified
ref0=0.
do j=1,63
ref0=ref0 + s3(mrs(j),j)
enddo
p1=-1.e30
p2=-1.e30
do k=1,ntot
sum=0.
ref=ref0
do j=1,63
i=ncode(j,k)+1
sum=sum + s3(i,j)
if(i.eq.mrs(j)) then
ref=ref - s3(i,j) + s3(mrs2(j),j)
if(flip.gt.0.0 .or. nflip(k).lt.0) then !Skip CQ msg if flip=-1
sum=0.
ref=ref0
do j=1,63
i=ncode(j,k)+1
sum=sum + s3(i,j)
if(i.eq.mrs(j)) then
ref=ref - s3(i,j) + s3(mrs2(j),j)
endif
enddo
p=sum/ref
pp(k)=p
if(p.gt.p1) then
p1=p
ip1=k
endif
enddo
p=sum/ref
pp(k)=p
if(p.gt.p1) then
p1=p
ip1=k
endif
enddo
@ -131,10 +143,18 @@ C Insert CQ message unless sync=OOO (flip=-1).
if(pp(i).gt.p2 .and. pp(i).ne.p1) p2=pp(i)
enddo
C ### Find out why this needs to be here ###
C ### It's OK without it, in Linux, if compiled without optimization.
! rewind 77
! write(77,*) p1,p2
if(mode65.eq.1) bias=max(1.12*p2,0.335)
if(mode65.eq.2) bias=max(1.08*p2,0.405)
if(mode65.ge.4) bias=max(1.04*p2,0.505)
if(p2.eq.p1) stop 'Error in deep65'
qual=100.0*(p1-bias)
decoded=' '
c=' '
@ -145,6 +165,7 @@ C Insert CQ message unless sync=OOO (flip=-1).
qual=0.
endif
decoded(22:22)=c
C Make sure everything is upper case.
do i=1,22
if(decoded(i:i).ge.'a' .and. decoded(i:i).le.'z')

View File

@ -14,7 +14,7 @@ C mr2prob probability that mr2sym was the transmitted value
real*4 signal(64,63)
real*8 fs(64)
integer mrsym(63),mrprob(63),mr2sym(63),mr2prob(63)
common/tmp9/ mrs(63),mrs2(63)
common/mrscom/ mrs(63),mrs2(63)
afac=1.1 * float(nadd)**0.64
scale=255.999

50
display.f Normal file
View File

@ -0,0 +1,50 @@
subroutine display(nutc)
parameter (MAXLINES=500)
integer indx(MAXLINES)
character*80 line(MAXLINES)
real freqkHz(MAXLINES)
integer utc(MAXLINES)
real*8 f0
ftol=0.02
rewind 26
do i=1,MAXLINES
read(26,1010,end=10) line(i)
1010 format(a80)
read(line(i),1020) f0,ndf,utc(i)
1020 format(f7.3,i5,26x,i5)
freqkHz(i)=1000.d0*(f0-144.d0) + 0.001d0*ndf
enddo
10 nz=i-1
call indexx(nz,freqkHz,indx)
nstart=1
rewind 24
write(24,3101) line(indx(1))
3101 format(a80)
do i=2,nz
j0=indx(i-1)
j=indx(i)
if(freqkHz(j)-freqkHz(j0).gt.ftol) then
if(nstart.eq.0) write(24,3101)
endfile 24
if(nstart.eq.1) then
!@@@ call sysqqq('sort -k 1.40 fort.24 | uniq > fort.13')
nstart=0
else
!@@@ call sysqqq('sort -k 1.40 fort.24 | uniq >> fort.13')
endif
rewind 24
endif
if(i.eq.nz) write(24,3101)
write(24,3101) line(j)
j0=j
enddo
endfile 24
!@@@ call sysqqq('sort -k 1.40 fort.24 | uniq >> fort.13')
return
end

View File

@ -1,28 +1,39 @@
subroutine extract(s3,nadd,ncount,decoded)
subroutine extract(s3,nadd,ncount,nhist,decoded)
real s3(64,63)
real tmp(4032)
character decoded*22
integer era(51),dat4(12),indx(63)
integer era(51),dat4(12),indx(64)
integer mrsym(63),mr2sym(63),mrprob(63),mr2prob(63)
logical first
data first/.true./,nsec1/0/
save
call demod64a(s3,nadd,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow)
nfail=0
1 call demod64a(s3,nadd,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow)
if(ntest.lt.50 .or. nlow.gt.20) then
ncount=-999 !Flag bad data
go to 900
endif
call chkhist(mrsym,nhist,ipk)
if(nhist.ge.20) then
nfail=nfail+1
call pctile(s3,tmp,4032,50,base) ! ### or, use ave from demod64a
do j=1,63
s3(ipk,j)=base
enddo
go to 1
endif
call graycode(mrsym,63,-1)
call interleave63(mrsym,-1)
call interleave63(mrprob,-1)
ndec=1
nemax=30
nemax=30 !Was 200 (30)
maxe=8
xlambda=15.0
xlambda=12.0 !Was 15 (12)
if(ndec.eq.1) then
call graycode(mr2sym,63,-1)
@ -35,9 +46,9 @@
call flushqqq(22)
call runqqq('kvasd.exe','-q',iret)
if(iret.ne.0) then
if(first) write(*,1000)
if(first) write(*,1000) iret
1000 format('Error in KV decoder, or no KV decoder present.'/
+ 'Using BM algorithm.')
+ 'Return code:',i8,'. Will use BM algorithm.')
ndec=0
first=.false.
go to 20

104
filbig.f Normal file
View File

@ -0,0 +1,104 @@
subroutine filbig(id,nmax,f0,newdat,c4a,c4b,n4)
C Filter and downsample complex data for X and Y polarizations,
C stored in array id(4,nmax). Output is downsampled from 96000 Hz
C to 1500 Hz, and the low-pass filter has f_cutoff = 375 Hz and
C f_stop = 750 Hz.
parameter (NFFT1=5376000,NFFT2=77175)
integer*2 id(4,nmax) !Input data
complex c4a(NFFT2),c4b(NFFT2) !Output data
complex ca(NFFT1),cb(NFFT1) !FFTs of input
real*8 df
real halfpulse(8) !Impulse response of filter (one side)
complex cfilt(NFFT2) !Filter (complex; imag = 0)
real rfilt(NFFT2) !Filter (real)
integer*8 plan1,plan2,plan3,plan4,plan5
logical first
include 'fftw3.f'
equivalence (rfilt,cfilt)
data first/.true./
data halfpulse/114.97547150,36.57879257,-20.93789101,
+ 5.89886379,1.59355187,-2.49138308,0.60910773,-0.04248129/
save
if(first) then
npatience=FFTW_ESTIMATE
! npatience=FFTW_MEASURE
C Plan the FFTs just once
call sfftw_plan_dft_1d_(plan1,NFFT1,ca,ca,
+ FFTW_BACKWARD,npatience)
call sfftw_plan_dft_1d_(plan2,NFFT1,cb,cb,
+ FFTW_BACKWARD,npatience)
call sfftw_plan_dft_1d_(plan3,NFFT2,c4a,c4a,
+ FFTW_FORWARD,npatience)
call sfftw_plan_dft_1d_(plan4,NFFT2,c4b,c4b,
+ FFTW_FORWARD,npatience)
call sfftw_plan_dft_1d_(plan5,NFFT2,cfilt,cfilt,
+ FFTW_BACKWARD,npatience)
C Convert impulse response to filter function
do i=1,NFFT2
cfilt(i)=0.
enddo
fac=0.00625/NFFT1
cfilt(1)=fac*halfpulse(1)
do i=2,8
cfilt(i)=fac*halfpulse(i)
cfilt(NFFT2+2-i)=fac*halfpulse(i)
enddo
call sfftw_execute_(plan5)
base=cfilt(NFFT2/2+1)
do i=1,NFFT2
rfilt(i)=real(cfilt(i))-base
enddo
df=96000.d0/NFFT1
first=.false.
endif
C When new data comes along, we need to compute a new "big FFT"
C If we just have a new f0, continue with the existing ca and cb.
if(newdat.ne.0) then
nz=min(nmax,NFFT1)
do i=1,nz
ca(i)=cmplx(float(id(1,i)),float(id(2,i)))
cb(i)=cmplx(float(id(3,i)),float(id(4,i)))
enddo
if(nmax.lt.NFFT1) then
do i=nmax+1,NFFT1
ca(i)=0.
cb(i)=0.
enddo
endif
call sfftw_execute_(plan1)
call sfftw_execute_(plan2)
newdat=0
endif
C NB: f0 is the frequency at which we want our filter centered.
C i0 is the bin number in ca and cb closest to f0.
i0=nint(f0/df) + 1
nh=NFFT2/2
do i=1,nh !Copy data into c4a and c4b,
j=i0+i-1 !and apply the filter function
c4a(i)=rfilt(i)*ca(j)
c4b(i)=rfilt(i)*cb(j)
enddo
do i=nh+1,NFFT2
j=i0+i-1-NFFT2
if(j.lt.1) j=j+NFFT2
c4a(i)=rfilt(i)*ca(j)
c4b(i)=rfilt(i)*cb(j)
enddo
C Do the short reverse transform, to go back to time domain.
call sfftw_execute_(plan3)
call sfftw_execute_(plan4)
n4=min(nmax/64,NFFT2)
return
end

View File

@ -19,8 +19,9 @@ C The transform will be real and returned to the input array.
parameter (NPMAX=100)
complex a(nfft)
complex aa(32768)
integer nn(NPMAX),ns(NPMAX),nf(NPMAX),nl(NPMAX)
integer plan(NPMAX)
integer*8 plan(NPMAX)
data nplan/0/
include 'fftw3.f'
save
@ -43,7 +44,12 @@ C The transform will be real and returned to the input array.
C Planning: FFTW_ESTIMATE, FFTW_MEASURE, FFTW_PATIENT, FFTW_EXHAUSTIVE
nspeed=FFTW_ESTIMATE
if(nfft.le.16384) nspeed=FFTW_MEASURE
nspeed=FFTW_MEASURE
if(nfft.le.32768) then
do j=1,nfft
aa(j)=a(j)
enddo
endif
if(isign.eq.-1 .and. iform.eq.1) then
call sfftw_plan_dft_1d_(plan(i),nfft,a,a,
+ FFTW_FORWARD,nspeed)
@ -57,19 +63,19 @@ C Planning: FFTW_ESTIMATE, FFTW_MEASURE, FFTW_PATIENT, FFTW_EXHAUSTIVE
else
stop 'Unsupported request in four2a'
endif
i=nplan
! write(*,3001) i,nn(i),ns(i),nf(i),nl(i),plan(i)
! 3001 format(6i10)
if(nfft.le.32768) then
do j=1,nfft
a(j)=aa(j)
enddo
endif
10 call sfftw_execute_(plan(i))
return
999 do i=1,nplan
! print*,i,nn(i),ns(i),nf(i),nl(i),plan(i)
call sfftw_destroy_plan_(plan(i))
enddo
! print*,'FFTW plans destroyed:',nplan
return
end

View File

@ -69,6 +69,15 @@ subroutine ftn_init
err=940)
#endif
#ifdef Win32
open(19,file=appdir(:iz)//'/bandmap.txt',status='unknown', &
share='denynone',err=910)
#else
open(19,file=appdir(:iz)//'/bandmap.txt',status='unknown', &
err=910)
#endif
endfile 19
#ifdef Win32
open(21,file=appdir(:iz)//'/ALL.TXT',status='unknown', &
access='append',share='denynone',err=950)

View File

@ -164,6 +164,22 @@ def testmsgs():
tx5.insert(0,"@1000")
tx6.insert(0,"@2000")
#------------------------------------------------------ bandmap
def bandmap(event=NONE):
global Version,bm,bm_geom,bmtext
bm=Toplevel(root)
bm.geometry(bm_geom)
if g.Win32: bm.iconbitmap("wsjt.ico")
iframe_bm1 = Frame(bm, bd=1, relief=SUNKEN)
bmtext=Text(iframe_bm1, height=35, width=45, bg="Navy", fg="yellow")
bmtext.pack(side=LEFT, fill=X, padx=1)
bmsb = Scrollbar(iframe_bm1, orient=VERTICAL, command=bmtext.yview)
bmsb.pack(side=RIGHT, fill=Y)
bmtext.configure(yscrollcommand=bmsb.set)
# bmtext.insert(END,'144.103 CQ EA3DXU JN11\n')
# bmtext.insert(END,'144.118 OH6KTL RA3AQ KO85 OOO')
iframe_bm1.pack(expand=1, fill=X, padx=4)
#------------------------------------------------------ logqso
def logqso(event=NONE):
t=time.strftime("%Y-%b-%d,%H:%M",time.gmtime())
@ -1070,22 +1086,6 @@ def plot_yellow():
xy2.append(n)
graph1.create_line(xy2,fill="yellow")
#------------------------------------------------------ bandmap
def bandmap(event=NONE):
global Version,bm,bm_geom
bm=Toplevel(root)
bm.geometry(bm_geom)
if g.Win32: bm.iconbitmap("wsjt.ico")
iframe_bm1 = Frame(bm, bd=1, relief=SUNKEN)
text=Text(iframe_bm1, height=35, width=32, bg="Navy", fg="yellow")
text.pack(side=LEFT, fill=X, padx=1)
sb = Scrollbar(iframe_bm1, orient=VERTICAL, command=text.yview)
sb.pack(side=RIGHT, fill=Y)
text.configure(yscrollcommand=sb.set)
text.insert(END,'144.103 CQ EA3DXU JN11\n')
text.insert(END,'144.118 OH6KTL RA3AQ KO85 OOO')
iframe_bm1.pack(expand=1, fill=X, padx=4)
#------------------------------------------------------ update
def update():
global root_geom,isec0,naz,nel,ndmiles,ndkm,nopen, \
@ -1179,6 +1179,10 @@ def update():
bdecode.configure(bg='gray85',activebackground='gray95')
if Audio.gcom2.ndecoding: #Set button bg=light_blue while decoding
bdecode.configure(bg='#66FFFF',activebackground='#66FFFF')
# print 'A'
Audio.map65a0() # @@@ Temporary @@@
# print 'B'
tx1.configure(bg='white')
tx2.configure(bg='white')
tx3.configure(bg='white')
@ -1251,6 +1255,21 @@ def update():
avetext.insert(END,lines[0])
avetext.insert(END,lines[1])
# avetext.configure(state=DISABLED)
try:
f=open(appdir+'/bandmap.txt',mode='r')
lines=f.readlines()
f.close()
except:
lines=""
bmtext.configure(state=NORMAL)
bmtext.insert(END,' Freq DF Pol UTC\n')
bmtext.insert(END,'--------------------------------------------\n')
for i in range(len(lines)):
bmtext.insert(END,lines[i])
bmtext.see(END)
Audio.gcom2.ndecdone=2
if g.cmap != cmap0:
@ -1744,7 +1763,6 @@ msg7=Message(iframe6, text=' ', width=300,relief=SUNKEN)
msg7.pack(side=RIGHT, fill=X, padx=1)
iframe6.pack(expand=1, fill=X, padx=4)
frame.pack()
ldate.after(100,update)
lauto=0
isync=1

288
map65a.f Normal file
View File

@ -0,0 +1,288 @@
subroutine map65a
C Processes timf2 data from Linrad to find and decode JT65 signals.
parameter (NMAX=60*96000) !Samples per 60 s file
parameter (MAXMSG=1000) !Size of decoded message list
integer*2 id(4,NMAX) !46 MB: raw data from Linrad timf2
parameter (NFFT=32768) !Half symbol = 17833 samples;
real ss(4,322,NFFT) !169 MB: half-symbol spectra
real savg(4,NFFT)
real tavg(-50:50) !Temp for finding local base level
real base(4) !Local basel level at 4 pol'ns
real tmp (200) !Temp storage for pctile sorting
real short(3,NFFT) !SNR dt ipol for potential shorthands
real sig(MAXMSG,30) !Parameters of detected signals
real a(5)
character*22 msg(MAXMSG)
character*3 shmsg0(4),shmsg
character arg*12,infile*11,outfile*11
integer indx(MAXMSG),nsiz(MAXMSG)
logical done(MAXMSG)
integer rfile3
character decoded*22,blank*22,cbad*1
data blank/' '/
data shmsg0/'ATT','RO ','RRR','73 '/
tskip=0.
fselect=0.
! fselect=103.0
nmin=1
infile='061111.0745'
C Initialize some constants
! open(22,file='kvasd.dat',access='direct',recl=1024,
! + status='unknown')
open(23,file='CALL3.TXT',status='unknown')
! nbytes=8*(4*96000+9000) !Empirical, for 061111_0744.dat.48
! nskip=8*nint(96000*(tskip+4.09375))
! n=rfile3(infile,id,nskip) !Skip to start of minute
! if(n.ne.nskip) go to 9999
df=96000.0/NFFT !df = 96000/NFFT = 2.930 Hz
fa=0.0
fb=60000.0
ia=nint((fa+23000.0)/df + 1.0) ! 23000 = 48000 - 25000
ib=nint((fb+23000.0)/df + 1.0)
ftol=0.020 !Frequency tolerance (kHz)
kk=0
nkk=1
do nfile=1,nmin
! n=rfile3(infile,id,8*NMAX) !Read 60 s of data (approx 46 MB)
n=8*NMAX
call rfile3a(infile,id,n,ierr)
newdat=1
nz=n/8
read(infile(8:11),*) nutc
if(fselect.gt.0.0) then
nfilt=1 !nfilt=2 is faster for selected freq
freq=fselect+1.600
nflip=-1 !May need to try both +/- 1
ipol=4 !Try all four?
dt=2.314240 !Not needed?
call decode1a(id,newdat,nfilt,freq,nflip,ipol,sync2,
+ a,dt,pol,nkv,nhist,qual,decoded)
write(11,1010) nutc,nsync1,nsync2,dt,ndf,decoded,
+ nkv,nqual
1010 format(i4.4,i5,i4,f5.1,i5,2x,a22,2i3)
if(nfile.eq.1) go to 999
endif
nfilt=1
do i=1,NFFT
short(1,i)=0.
short(2,i)=0.
short(3,i)=0.
enddo
call symspec(id,nz,ss,savg)
freq0=-999.
sync10=-999.
fshort0=-999.
sync20=-999.
ntry=0
do i=ia,ib !Search over freq range
freq=0.001*((i-1)*df - 23000) + 100.0
C Find the local base level for each polarization; update every 10 bins.
if(mod(i-ia,10).eq.0) then
do jp=1,4
do ii=-50,50
tavg(ii)=savg(jp,i+ii)
enddo
call pctile(tavg,tmp,101,50,base(jp))
enddo
endif
C Find max signal at this frequency
smax=0.
do jp=1,4
if(savg(jp,i)/base(jp).gt.smax) smax=savg(jp,i)/base(jp)
enddo
if(smax.gt.1.1) then
ntry=ntry+1
C Look for JT65 sync patterns and shorthand square-wave patterns.
call ccf65(ss(1,1,i),sync1,ipol,dt,flipk,
+ syncshort,snr2,ipol2,dt2)
shmsg=' '
C Is there a shorthand tone above threshold?
if(syncshort.gt.1.0) then
C### Do shorthand AFC here (or maybe after finding a pair?) ###
short(1,i)=syncshort
short(2,i)=dt2
short(3,i)=ipol2
C Check to see if lower tone of shorthand pair was found.
do j=2,4
i0=i-nint(j*53.8330078/df)
C Should this be i0 +/- 1, or just i0?
C Should we also insist that difference in DT be either 1.5 or -1.5 s?
if(short(1,i0).gt.1.0) then
fshort=0.001*((i0-1)*df - 23000) + 100.0
C Keep only the best candidate within ftol.
if(fshort-fshort0.le.ftol .and. sync2.gt.sync20
+ .and. nkk.eq.2) kk=kk-1
if(fshort-fshort0.gt.ftol .or.
+ sync2.gt.sync20) then
kk=kk+1
sig(kk,1)=nfile
sig(kk,2)=nutc
sig(kk,3)=fshort
sig(kk,4)=syncshort
sig(kk,5)=dt2
sig(kk,6)=45*(ipol2-1)/57.2957795
sig(kk,7)=0
sig(kk,8)=snr2
sig(kk,9)=0
sig(kk,10)=0
! sig(kk,11)=rms0
sig(kk,12)=savg(ipol2,i)
sig(kk,13)=0
sig(kk,14)=0
sig(kk,15)=0
sig(kk,16)=0
! sig(kk,17)=0
sig(kk,18)=0
msg(kk)=shmsg0(j)
fshort0=fshort
sync20=sync2
nkk=2
endif
endif
enddo
endif
C Is sync1 above threshold?
if(sync1.gt.1.0) then
C Keep only the best candidate within ftol.
C (Am I deleting any good decodes by doing this? Any harm in omitting
C these statements??)
if(freq-freq0.le.ftol .and. sync1.gt.sync10 .and.
+ nkk.eq.1) kk=kk-1
if(freq-freq0.gt.ftol .or. sync1.gt.sync10) then
nflip=nint(flipk)
call decode1a(id,newdat,nfilt,freq,nflip,ipol,
+ sync2,a,dt,pol,nkv,nhist,qual,decoded)
kk=kk+1
sig(kk,1)=nfile
sig(kk,2)=nutc
sig(kk,3)=freq
sig(kk,4)=sync1
sig(kk,5)=dt
sig(kk,6)=pol
sig(kk,7)=flipk
sig(kk,8)=sync2
sig(kk,9)=nkv
sig(kk,10)=qual
! sig(kk,11)=rms0
sig(kk,12)=savg(ipol,i)
sig(kk,13)=a(1)
sig(kk,14)=a(2)
sig(kk,15)=a(3)
sig(kk,16)=a(4)
! sig(kk,17)=a(5)
sig(kk,18)=nhist
msg(kk)=decoded
freq0=freq
sync10=sync1
nkk=1
endif
endif
endif
enddo
! write(*,1010)
C Trim the list and produce a sorted index and sizes of groups.
C (Should trimlist remove all but best SNR for given UTC and message content?)
call trimlist(sig,kk,indx,nsiz,nz)
do i=1,kk
done(i)=.false.
enddo
j=0
ilatest=-1
do n=1,nz
ifile0=0
do m=1,nsiz(n)
i=indx(j+m)
ifile=sig(i,1)
if(ifile.gt.ifile0 .and.msg(i).ne.blank) then
ilatest=i
ifile0=ifile
endif
enddo
i=ilatest
if(i.ge.1) then
if(.not.done(i)) then
done(i)=.true.
nutc=sig(i,2)
freq=sig(i,3)
sync1=sig(i,4)
dt=sig(i,5)
npol=nint(57.2957795*sig(i,6))
flip=sig(i,7)
sync2=sig(i,8)
nkv=sig(i,9)
nqual=min(sig(i,10),10.0)
! rms0=sig(i,11)
nsavg=sig(i,12) !Was used for diagnostic ...
do k=1,5
a(k)=sig(i,12+k)
enddo
nhist=sig(i,18)
decoded=msg(i)
if(flip.lt.0.0) then
do i=22,1,-1
if(decoded(i:i).ne.' ') go to 10
enddo
stop 'Error in message format'
10 if(i.le.18) decoded(i+2:i+4)='OOO'
endif
nkHz=nint(freq-1.600)
f0=144.0+0.001*nkHz
ndf=nint(1000.0*(freq-1.600-nkHz))
ndf0=nint(a(1))
ndf1=nint(a(2))
ndf2=nint(a(3))
nsync1=sync1
nsync2=nint(10.0*log10(sync2)) - 40 !### empirical ###
cbad=' '
if(abs(f0-144.103).lt.0.001) then
write(11,1010) nutc,nsync1,nsync2,dt,ndf,decoded,
+ nkv,nqual
endif
write(19,1012) f0,ndf,npol,nutc,decoded
1012 format(f7.3,i5,i4,i5.4,2x,a22)
write(26,1014) f0,ndf,ndf0,ndf1,ndf2,dt,npol,nsync1,
+ nsync2,nutc,decoded,nkv,nqual,nhist
1014 format(f7.3,i5,3i3,f5.1,i5,i3,i4,i5.4,2x,a22,3i3)
endif
endif
j=j+nsiz(n)
enddo
call display(nutc)
! if(nfile.ge.1) go to 999
100 continue
enddo
999 call four2a(cx,-1,1,1,1) !Destroy the FFTW plans
9999 end

View File

@ -1,7 +1,8 @@
#include <stdio.h>
#include <samplerate.h>
int resample_( float din[], float dout[], double *samfac, int *jz)
int resample_(float din[], int *jzin, int *conv_type, int *channels,
double *samfac, float dout[], int *jzout)
{
SRC_DATA src_data;
int input_len;
@ -10,7 +11,7 @@ int resample_( float din[], float dout[], double *samfac, int *jz)
double src_ratio;
src_ratio=*samfac;
input_len=*jz;
input_len=*jzin;
output_len=(int) (input_len*src_ratio);
src_data.data_in=din;
@ -19,10 +20,7 @@ int resample_( float din[], float dout[], double *samfac, int *jz)
src_data.input_frames=input_len;
src_data.output_frames=output_len;
ierr=src_simple(&src_data,2,1);
*jz=output_len;
/* printf("%d %d %d %d %f\n",input_len,output_len,src_data.input_frames_used,
src_data.output_frames_gen,src_ratio);
*/
ierr=src_simple(&src_data,*conv_type,*channels);
*jzout=output_len;
return ierr;
}

View File

@ -15,26 +15,3 @@ subroutine runqqq(fname,cmnd,iret)
return
end subroutine runqqq
subroutine flushqqq(lu)
#ifdef Win32
use dfport
#endif
call flush(lu)
return
end subroutine flushqqq
subroutine sleepqqq(n)
#ifdef Win32
use dflib
call sleepqq(n)
#else
call usleep(n*1000)
#endif
return
end subroutine sleepqqq

67
symspec.f Normal file
View File

@ -0,0 +1,67 @@
subroutine symspec(id,nz,ss,savg)
C Compute spectra at four polarizations, using half-symbol steps.
parameter (NFFT=32768)
integer*2 id(4,nz)
real ss(4,322,NFFT)
real savg(4,NFFT)
complex cx(NFFT),cy(NFFT) ! pad to 32k with zeros
complex z
real*8 ts,hsym
fac=1.e-4
hsym=2048.d0*96000.d0/11025.d0 !Samples per half symbol
npts=hsym !Integral samples per half symbol
nsteps=322 !Half symbols per transmission
do ip=1,4
do i=1,NFFT
savg(ip,i)=0.
enddo
enddo
ts=1.d0 - hsym
do n=1,nsteps
ts=ts+hsym !Update exact sample pointer
i0=ts !Starting sample pointer
do i=1,npts !Copy data to FFT arrays
xr=fac*id(1,i0+i)
xi=fac*id(2,i0+i)
cx(i)=cmplx(xr,xi)
yr=fac*id(3,i0+i)
yi=fac*id(4,i0+i)
cy(i)=cmplx(yr,yi)
enddo
do i=npts+1,NFFT !Pad to 32k with zeros
cx(i)=0.
cy(i)=0.
enddo
call four2a(cx,NFFT,1,1,1) !Do the FFTs
call four2a(cy,NFFT,1,1,1)
do i=1,NFFT !Save and accumulate power spectra
s=real(cx(i))**2 + aimag(cx(i))**2
ss(1,n,i)=s ! Pol = 0
savg(1,i)=savg(1,i) + s
z=cx(i) + cy(i)
s=0.5*(real(z)**2 + aimag(z)**2)
ss(2,n,i)=s ! Pol = 45
savg(2,i)=savg(2,i) + s
s=real(cy(i))**2 + aimag(cy(i))**2
ss(3,n,i)=s ! Pol = 90
savg(3,i)=savg(3,i) + s
z=cx(i) - cy(i)
s=0.5*(real(z)**2 + aimag(z)**2)
ss(4,n,i)=s ! Pol = 135
savg(4,i)=savg(4,i) + s
enddo
enddo
return
end