w4kek/WSJT-X
1
0
Fork 0
mirror of https://github.com/saitohirga/WSJT-X.git synced 2025-03-23 20:48:33 -04:00
Code Issues Packages Projects Releases Wiki Activity

Major prund files from Q65W source tree.

Browse Source
This commit is contained in:
Joe Taylor 2022-12-12 11:14:58 -05:00
parent 97b58a387b
commit ecaa0b8861
28 changed files with 0 additions and 2999 deletions

38
q65w/libm65/CMakeLists.txt
View File

@ -3,54 +3,33 @@ set (libm65_FSRCS
wideband_sync.f90
# Non-module Fortran routines:
afc65b.f90
astro.f90
astro0.f90
astrosub.f90
averms.f90
badmsg.f90
ccf2.f90
ccf65.f90
cgen65.f90
chkhist.f90
chkmsg.f90
coord.f90
dcoord.f90
decode0.f90
decode65b.f90
deep65.f90
deg2grid.f90
demod64a.f90
display.f90
dot.f90
dpol.f90
encode65.f90
extract.f90
fchisq.f90
fchisq0.f90
fil6521.f90
filbig.f90
fmtmsg.f90
four2a.f90
ftninit.f90
ftnquit.f90
q65b.f90
gen65.f90
gen_q65_cwave.f90
gen_q65_wave.f90
geocentric.f90
getdphi.f90
getpfx1.f90
getpfx2.f90
graycode.f90
graycode65.f90
grid2deg.f90
grid2k.f90
indexx.f90
interleave63.f90
iqcal.f90
iqfix.f90
jt65code.f90
k2grid.f90
lorentzian.f90
m65c.f90
@ -58,15 +37,12 @@ set (libm65_FSRCS
moon2.f90
moondop.f90
nchar.f90
noisegen.f90
packjt.f90
pfxdump.f90
recvpkt.f90
rfile3a.f90
s3avg.f90
sec_midn.f90
set.f90
setup65.f90
shell.f90
sleep_msec.f90
smo.f90
@ -78,7 +54,6 @@ set (libm65_FSRCS
trimlist.f90
twkfreq.f90
twkfreq_xy.f90
txpol.f90
wavhdr.f90
f77_wisdom.f
@ -128,16 +103,3 @@ set_property (SOURCE ${libm65_C_and_CXXSRCS} APPEND PROPERTY OBJECT_DEPENDS ${CM
add_library (m65impl STATIC ${libm65_FSRCS} ${libm65_CSRCS} ${libm65_CXXSRCS})
target_link_libraries (m65impl wsjt_fort wsjt_cxx Qt5::Core)
add_executable (mapsim mapsim.f90)
target_link_libraries (mapsim m65impl ${FFTW3_LIBRARIES})
#add_executable (synctest synctest.f90)
#target_link_libraries (synctest m65impl ${FFTW3_LIBRARIES})
if (WIN32)
install (
TARGETS mapsim
RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} COMPONENT runtime
BUNDLE DESTINATION . COMPONENT runtime
)
endif ()

70
q65w/libm65/afc65b.f90
View File

@ -1,70 +0,0 @@
subroutine afc65b(cx,cy,npts,fsample,nflip,ipol,xpol,ndphi,a,ccfbest,dtbest)
logical xpol
complex cx(npts)
complex cy(npts)
real a(5),deltaa(5)
a(1)=0.
a(2)=0.
a(3)=0.
a(4)=45.0*(ipol-1.0)
deltaa(1)=2.0
deltaa(2)=2.0
deltaa(3)=2.0
deltaa(4)=22.5
deltaa(5)=0.05
nterms=3
if(xpol) nterms=4
! Don't fit polarization when solving for dphi
if(ndphi.ne.0) nterms=3
! Start the iteration
chisqr=0.
chisqr0=1.e6
do iter=1,3 !One iteration is enough?
do j=1,nterms
chisq1=fchisq(cx,cy,npts,fsample,nflip,a,ccfmax,dtmax)
fn=0.
delta=deltaa(j)
10 a(j)=a(j)+delta
chisq2=fchisq(cx,cy,npts,fsample,nflip,a,ccfmax,dtmax)
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=fchisq(cx,cy,npts,fsample,nflip,a,ccfmax,dtmax)
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.
enddo
chisqr=fchisq(cx,cy,npts,fsample,nflip,a,ccfmax,dtmax)
if(chisqr/chisqr0.gt.0.9999) go to 30
chisqr0=chisqr
enddo
30 ccfbest=ccfmax * (1378.125/fsample)**2
dtbest=dtmax
if(a(4).lt.0.0) a(4)=a(4)+180.0
if(a(4).ge.180.0) a(4)=a(4)-180.0
if(nint(a(4)).eq.180) a(4)=0.
ipol=nint(a(4)/45.0) + 1
if(ipol.gt.4) ipol=ipol-4
return
end subroutine afc65b

20
q65w/libm65/averms.f90
View File

@ -1,20 +0,0 @@
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

46
q65w/libm65/badmsg.f90
View File

@ -1,46 +0,0 @@
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

45
q65w/libm65/ccf2.f90
View File

@ -1,45 +0,0 @@
subroutine ccf2(ss,nz,nflip,ccfbest,lagpk)
! parameter (LAGMAX=60)
parameter (LAGMAX=200)
real ss(nz)
real ccf(-LAGMAX:LAGMAX)
integer npr(126)
! The JT65 pseudo-random sync pattern:
data npr/ &
1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0, &
0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1, &
0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1, &
0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1, &
1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1, &
0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1, &
1,1,1,1,1,1/
save
ccfbest=0.
lag1=-LAGMAX
lag2=LAGMAX
do lag=lag1,lag2
s0=0.
s1=0.
do i=1,126
j=2*(8*i + 43) + lag
if(j.ge.1 .and. j.le.nz-8) then
x=ss(j)+ss(j+8) !Add two half-symbol contributions
if(npr(i).eq.0) then
s0=s0 + x
else
s1=s1 + x
endif
endif
enddo
ccf(lag)=nflip*(s1-s0)
if(ccf(lag).gt.ccfbest) then
ccfbest=ccf(lag)
lagpk=lag
endif
enddo
return
end subroutine ccf2

128
q65w/libm65/ccf65.f90
View File

@ -1,128 +0,0 @@
subroutine ccf65(ss,nhsym,ssmax,sync1,ipol1,jpz,dt1,flipk, &
syncshort,snr2,ipol2,dt2)
parameter (NFFT=512,NH=NFFT/2)
real ss(4,322) !Input: half-symbol powers, 4 pol'ns
real s(NFFT) !CCF = ss*pr
complex cs(0:NH) !Complex FT of s
real s2(NFFT) !CCF = ss*pr2
complex cs2(0:NH) !Complex FT of s2
real pr(NFFT) !JT65 pseudo-random sync pattern
complex cpr(0:NH) !Complex FT of pr
real pr2(NFFT) !JT65 shorthand pattern
complex cpr2(0:NH) !Complex FT of pr2
real tmp1(322)
real ccf(-11:54,4)
logical first
integer npr(126)
data first/.true./
equivalence (s,cs),(pr,cpr),(s2,cs2),(pr2,cpr2)
save
! The JT65 pseudo-random sync pattern:
data npr/ &
1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0, &
0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1, &
0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1, &
0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1, &
1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1, &
0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1, &
1,1,1,1,1,1/
if(first) then
! Initialize pr, pr2; compute cpr, cpr2.
fac=1.0/NFFT
do i=1,NFFT
pr(i)=0.
pr2(i)=0.
k=2*mod((i-1)/8,2)-1
if(i.le.NH) pr2(i)=fac*k
enddo
do i=1,126
j=2*i
pr(j)=fac*(2*npr(i)-1)
! Not sure why, but it works significantly better without the following line:
! pr(j-1)=pr(j)
enddo
call four2a(cpr,NFFT,1,-1,0)
call four2a(cpr2,NFFT,1,-1,0)
first=.false.
endif
syncshort=0.
snr2=0.
! Look for JT65 sync pattern and shorthand square-wave pattern.
ccfbest=0.
ccfbest2=0.
ipol1=1
ipol2=1
do ip=1,jpz !Do jpz polarizations
do i=1,nhsym-1
! s(i)=ss(ip,i)+ss(ip,i+1)
s(i)=min(ssmax,ss(ip,i)+ss(ip,i+1))
enddo
call pctile(s,nhsym-1,50,base)
s(1:nhsym-1)=s(1:nhsym-1)-base
s(nhsym:NFFT)=0.
call four2a(cs,NFFT,1,-1,0) !Real-to-complex FFT
do i=0,NH
cs2(i)=cs(i)*conjg(cpr2(i)) !Mult by complex FFT of pr2
cs(i)=cs(i)*conjg(cpr(i)) !Mult by complex FFT of pr
enddo
call four2a(cs,NFFT,1,1,-1) !Complex-to-real inv-FFT
call four2a(cs2,NFFT,1,1,-1) !Complex-to-real inv-FFT
do lag=-11,54 !Check for best JT65 sync
j=lag
if(j.lt.1) j=j+NFFT
ccf(lag,ip)=s(j)
if(abs(ccf(lag,ip)).gt.ccfbest) then
ccfbest=abs(ccf(lag,ip))
lagpk=lag
ipol1=ip
flipk=1.0
if(ccf(lag,ip).lt.0.0) flipk=-1.0
endif
enddo
!### Not sure why this is ever true???
if(sum(ccf).eq.0.0) return
!###
do lag=-11,54 !Check for best shorthand
ccf2=s2(lag+28)
if(ccf2.gt.ccfbest2) then
ccfbest2=ccf2
lagpk2=lag
ipol2=ip
endif
enddo
enddo
! Find rms level on baseline of "ccfblue", for normalization.
sumccf=0.
do lag=-11,54
if(abs(lag-lagpk).gt.1) sumccf=sumccf + ccf(lag,ipol1)
enddo
base=sumccf/50.0
sq=0.
do lag=-11,54
if(abs(lag-lagpk).gt.1) sq=sq + (ccf(lag,ipol1)-base)**2
enddo
rms=sqrt(sq/49.0)
sync1=-4.0
if(rms.gt.0.0) sync1=ccfbest/rms - 4.0
dt1=lagpk*(2048.0/11025.0) - 2.5
! Find base level for normalizing snr2.
do i=1,nhsym
tmp1(i)=ss(ipol2,i)
enddo
call pctile(tmp1,nhsym,40,base)
snr2=0.01
if(base.gt.0.0) snr2=0.398107*ccfbest2/base !### empirical
syncshort=0.5*ccfbest2/rms - 4.0 !### better normalizer than rms?
dt2=2.5 + lagpk2*(2048.0/11025.0)
return
end subroutine ccf65

99
q65w/libm65/cgen65.f90
View File

@ -1,99 +0,0 @@
subroutine cgen65(message,mode65,samfac,nsendingsh,msgsent,cwave,nwave)
! Encodes a JT65 message into a wavefile.
! Executes in 17 ms on opti-745.
use packjt
parameter (NMAX=60*96000) !Max length of wave file
character*22 message !Message to be generated
character*22 msgsent !Message as it will be received
character*3 cok !' ' or 'OOO'
real*8 t,dt,phi,f,f0,dfgen,dphi,twopi,samfac,tsymbol
complex cwave(NMAX) !Generated complex wave file
integer dgen(12)
integer sent(63)
logical first
integer nprc(126)
real pr(126)
data nprc/1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0, &
0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1, &
0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1, &
0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1, &
1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1, &
0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1, &
1,1,1,1,1,1/
data twopi/6.283185307179586476d0/,first/.true./
save
if(first) then
do i=1,126
pr(i)=2*nprc(i)-1
enddo
first=.false.
endif
call chkmsg(message,cok,nspecial,flip) !See if it's a shorthand
if(nspecial.eq.0) then
call packmsg(message,dgen,itype) !Pack message into 72 bits
nsendingsh=0
if(iand(dgen(10),8).ne.0) nsendingsh=-1 !Plain text flag
call rs_encode(dgen,sent)
call interleave63(sent,1) !Apply interleaving
call graycode(sent,63,1) !Apply Gray code
nsym=126 !Symbols per transmission
tsymbol=4096.d0/11025.d0 !Time per symbol
else
nsendingsh=1 !Flag for shorthand message
nsym=32
tsymbol=16384.d0/11025.d0
endif
! Set up necessary constants
dt=1.d0/(samfac*96000.d0)
f0=118*11025.d0/1024
dfgen=mode65*11025.d0/4096.d0
t=0.d0
phi=0.d0
k=0
j0=0
ndata=nsym*96000.d0*samfac*tsymbol
do i=1,ndata
t=t+dt
j=int(t/tsymbol) + 1 !Symbol number, 1-126
if(j.ne.j0) then
f=f0
if(nspecial.ne.0 .and. mod(j,2).eq.0) f=f0+10*nspecial*dfgen
if(nspecial.eq.0 .and. flip*pr(j).lt.0.0) then
k=k+1
f=f0+(sent(k)+2)*dfgen
endif
dphi=twopi*dt*f
j0=j
endif
phi=phi+dphi
if(phi.gt.twopi) phi=phi-twopi
xphi=phi
cwave(i)=cmplx(cos(xphi),-sin(xphi))
enddo
cwave(ndata+1:)=0
nwave=ndata + 48000
call unpackmsg(dgen,msgsent)
if(flip.lt.0.0) then
do i=22,1,-1
if(msgsent(i:i).ne.' ') goto 10
enddo
10 msgsent=msgsent(1:i)//' OOO'
endif
if(nsendingsh.eq.1) then
if(nspecial.eq.2) msgsent='RO'
if(nspecial.eq.3) msgsent='RRR'
if(nspecial.eq.4) msgsent='73'
endif
return
end subroutine cgen65

145
q65w/libm65/decode1a.f90
View File

@ -1,145 +0,0 @@
subroutine decode1a(dd,newdat,f0,nflip,mode65,nfsample,xpol, &
mycall,hiscall,hisgrid,neme,ndepth,nqd,dphi,ndphi, &
nutc,nkhz,ndf,ipol,ntol,sync2,a,dt,pol,nkv,nhist,nsum,nsave, &
qual,decoded)
! Apply AFC corrections to a candidate JT65 signal, then decode it.
use timer_module, only: timer
parameter (NMAX=60*96000) !Samples per 60 s
real*4 dd(4,NMAX) !92 MB: raw data from Linrad timf2
complex cx(NMAX/64), cy(NMAX/64) !Data at 1378.125 samples/s
complex c5x(NMAX/256),c5y(NMAX/256) !Data at 344.53125 Hz
complex c5a(512)
complex z
real s2(66,126)
real s3(64,63),sy(63)
real a(5)
logical first,xpol
character decoded*22
character mycall*12,hiscall*12,hisgrid*6
data first/.true./,jjjmin/1000/,jjjmax/-1000/
data nutc0/-999/,nhz0/-9999999/
save
! Mix sync tone to baseband, low-pass filter, downsample to 1378.125 Hz
dt00=dt
call timer('filbig ',0)
call filbig(dd,NMAX,f0,newdat,nfsample,xpol,cx,cy,n5)
! NB: cx, cy have sample rate 96000*77125/5376000 = 1378.125 Hz
call timer('filbig ',1)
if(mode65.eq.0) goto 900
sqa=0.
sqb=0.
do i=1,n5
sqa=sqa + real(cx(i))**2 + aimag(cx(i))**2
if(xpol) sqb=sqb + real(cy(i))**2 + aimag(cy(i))**2
enddo
sqa=sqa/n5
sqb=sqb/n5
! Find best DF, f1, f2, DT, and pol. Start by downsampling to 344.53125 Hz
if(xpol) then
z=cmplx(cos(dphi),sin(dphi))
cy(:n5)=z*cy(:n5) !Adjust for cable length difference
endif
! Add some zeros at start of c5 arrays -- empirical fix for negative DT's
nadd=1089
c5x(:nadd)=0.
call fil6521(cx,n5,c5x(nadd+1),n6)
if(xpol) then
c5y(:nadd)=0.
call fil6521(cy,n5,c5y(nadd+1),n6)
endif
n6=n6+nadd
fsample=1378.125/4.
a(5)=dt00
i0=nint((a(5)+0.5)*fsample) - 2 + nadd
if(i0.lt.1) then
write(13,*) 'i0 too small in decode1a:',i0,f0
flush(13)
i0=1
endif
nz=n6+1-i0
! We're looking only at sync tone here... so why not downsample by another
! factor of 1/8, say? Should be a significant execution speed-up.
! Best fit for DF, f1, f2, pol
call afc65b(c5x(i0),c5y(i0),nz,fsample,nflip,ipol,xpol,ndphi,a,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
! Apply AFC corrections to the time-domain signal
! Now we are back to using the 1378.125 Hz sample rate, enough to
! accommodate the full JT65C bandwidth.
call twkfreq_xy(cx,cy,n5,a)
! Compute spectrum at best polarization for each half symbol.
! Adding or subtracting a small number (e.g., 5) to j may make it decode.\
! NB: might want to try computing full-symbol spectra (nfft=512, even for
! submodes B and C).
nsym=126
nfft=512
j=(dt00+dtbest+2.685)*1378.125
if(j.lt.0) j=0
! Perhaps should try full-symbol-length FFTs even in B, C sub-modes?
! (Tried this, found no significant difference in decodes.)
do k=1,nsym
! do n=1,mode65
do n=1,1
do i=1,nfft
j=min(j+1,NMAX/64)
c5a(i)=aa*cx(j) + bb*cy(j)
enddo
call four2a(c5a,nfft,1,1,1)
if(n.eq.1) then
do i=1,66
! s2(i,k)=real(c5a(i))**2 + aimag(c5a(i))**2
jj=i
if(mode65.eq.2) jj=2*i-1
if(mode65.eq.4) jj=4*i-3
s2(i,k)=real(c5a(jj))**2 + aimag(c5a(jj))**2
enddo
else
do i=1,66
s2(i,k)=s2(i,k) + real(c5a(i))**2 + aimag(c5a(i))**2
enddo
endif
enddo
enddo
flip=nflip
call timer('dec65b ',0)
call decode65b(s2,flip,mycall,hiscall,hisgrid,mode65,neme,ndepth, &
nqd,nkv,nhist,qual,decoded,s3,sy)
dt=dt00 + dtbest + 1.7
call timer('dec65b ',1)
if(nqd.eq.1 .and. decoded.eq.' ') then
nhz=1000*nkhz + ndf
ihzdiff=min(500,ntol)
if(nutc.ne.nutc0 .or. abs(nhz-nhz0).ge.ihzdiff) syncbest=0.
if(sync2.gt.0.99999*syncbest) then
nsave=nsave+1
nsave=mod(nsave-1,64)+1
npol=nint(57.296*pol)
call s3avg(nsave,mode65,nutc,nhz,xdt,npol,ntol,s3,nsum,nkv,decoded)
syncbest=sync2
nhz0=nhz
endif
nutc0=nutc
endif
900 return
end subroutine decode1a

48
q65w/libm65/decode65b.f90
View File

@ -1,48 +0,0 @@
subroutine decode65b(s2,flip,mycall,hiscall,hisgrid,mode65,neme,ndepth, &
nqd,nkv,nhist,qual,decoded,s3,sy)
real s2(66,126)
real s3(64,63),sy(63)
logical first,ltext
character decoded*22,deepmsg*22
character mycall*12,hiscall*12,hisgrid*6
common/prcom/pr(126),mdat(126),mref(126,2),mdat2(126),mref2(126,2)
data first/.true./
save
if(first) call setup65
first=.false.
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)
enddo
k=mdat2(j) !Points to data symbol
if(flip.lt.0.0) k=mdat(j)
sy(j)=s2(1,k)
enddo
nadd=mode65
call extract(s3,nadd,ncount,nhist,decoded,ltext) !Extract the message
! Suppress "birdie messages" and other garbage decodes:
if(decoded(1:7).eq.'000AAA ') ncount=-1
if(decoded(1:7).eq.'0L6MWK ') ncount=-1
if(flip.lt.0.0 .and. ltext) ncount=-1
nkv=1
if(ncount.lt.0) then
nkv=0
decoded=' '
endif
qual=0.
if(ndepth.ge.1 .and. (nqd.eq.1 .or. flip.eq.1.0)) then
call deep65(s3,mode65,neme,flip,mycall,hiscall,hisgrid,deepmsg,qual)
if(nqd.ne.1 .and. qual.lt.10.0) qual=0.0
if(ndepth.lt.2 .and. qual.lt.6.0) qual=0.0
endif
if(nkv.eq.0 .and. qual.ge.1.0) decoded=deepmsg
return
end subroutine decode65b

170
q65w/libm65/deep65.f90
View File

@ -1,170 +0,0 @@
subroutine deep65(s3,mode65,neme,flip,mycall,hiscall,hisgrid,decoded,qual)
use timer_module, only: timer
parameter (MAXCALLS=10000,MAXRPT=63)
real s3(64,63)
character callsign*12,grid*4,message*22,hisgrid*6,c*1,ceme*3
character*12 mycall,hiscall
character*22 decoded,bestmsg
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)
real pp(2*MAXCALLS + 2 + MAXRPT)
common/mrscom/ mrs(63),mrs2(63)
common/c3com/ mcall3a
data rpt/'-01','-02','-03','-04','-05', &
'-06','-07','-08','-09','-10', &
'-11','-12','-13','-14','-15', &
'-16','-17','-18','-19','-20', &
'-21','-22','-23','-24','-25', &
'-26','-27','-28','-29','-30', &
'R-01','R-02','R-03','R-04','R-05', &
'R-06','R-07','R-08','R-09','R-10', &
'R-11','R-12','R-13','R-14','R-15', &
'R-16','R-17','R-18','R-19','R-20', &
'R-21','R-22','R-23','R-24','R-25', &
'R-26','R-27','R-28','R-29','R-30', &
'RO','RRR','73'/
save
if(mcall3a.eq.0) go to 30
call timer('deep65a ',0)
mcall3a=0
rewind 23
k=0
icall=0
do n=1,MAXCALLS
if(n.eq.1) then
callsign=hiscall
do i=4,12
if(ichar(callsign(i:i)).eq.0) callsign(i:i)=' '
enddo
grid=hisgrid(1:4)
if(ichar(grid(3:3)).eq.0) grid(3:3)=' '
if(ichar(grid(4:4)).eq.0) grid(4:4)=' '
else
read(23,1002,end=20) line
1002 format (A80)
if(line(1:4).eq.'ZZZZ') go to 20
if(line(1:2).eq.'//') go to 10
i1=index(line,',')
if(i1.lt.4) go to 10
i2=index(line(i1+1:),',')
if(i2.lt.5) go to 10
i2=i2+i1
i3=index(line(i2+1:),',')
if(i3.lt.1) i3=index(line(i2+1:),' ')
i3=i2+i3
callsign=line(1:i1-1)
grid=line(i1+1:i2-1)
ceme=line(i2+1:i3-1)
if(neme.eq.1 .and. ceme.ne.'EME') go to 10
endif
icall=icall+1
j1=index(mycall,' ') - 1
if(j1.le.-1) j1=12
if(j1.lt.3) j1=6
j2=index(callsign,' ') - 1
if(j2.le.-1) j2=12
if(j2.lt.3) j2=6
j3=index(mycall,'/') ! j3>0 means compound mycall
j4=index(callsign,'/') ! j4>0 means compound hiscall
callgrid(icall)=callsign(1:j2)
mz=1
! Allow MyCall + HisCall + rpt (?)
if(n.eq.1 .and. j3.lt.1 .and. j4.lt.1 .and. callsign(1:6).ne.' ') &
mz=MAXRPT+1
do m=1,mz
if(m.gt.1) grid=rpt(m-1)
if(j3.lt.1 .and.j4.lt.1) callgrid(icall)=callsign(1:j2)//' '//grid
message=mycall(1:j1)//' '//callgrid(icall)
k=k+1
testmsg(k)=message
call encode65(message,ncode(1,k))
! Insert CQ message
if(j4.lt.1) callgrid(icall)=callsign(1:j2)//' '//grid
message='CQ '//callgrid(icall)
k=k+1
testmsg(k)=message
call encode65(message,ncode(1,k))
enddo
10 continue
enddo
20 continue
ntot=k
call timer('deep65a ',1)
30 continue
call timer('deep65b ',0)
ref0=0.
do j=1,63
ref0=ref0 + s3(mrs(j),j)
enddo
p1=-1.e30
do k=1,ntot
pp(k)=0.
if(k.ge.2 .and. k.le.64 .and. flip.lt.0.0) cycle
! Test all messages if flip=+1; skip the CQ messages if flip=-1.
if(flip.gt.0.0 .or. testmsg(k)(1:3).ne.'CQ ') then
sum=0.
ref=ref0
do j=1,63
i=ncode(j,k)+1
sum=sum + s3(i,j)
if(i.eq.mrs(j)) ref=ref - s3(i,j) + s3(mrs2(j),j)
enddo
p=sum/ref
pp(k)=p
if(p.gt.p1) then
p1=p
ip1=k
bestmsg=testmsg(k)
endif
endif
enddo
p2=-1.e30
do i=1,ntot
if(pp(i).gt.p2 .and. testmsg(i).ne.bestmsg) p2=pp(i)
enddo
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 .and. p1.ne.-1.e30) then
open(77,file='error.log',status='unknown',access='append')
write(77,*) p1,p2,ip1,bestmsg
close(77)
endif
qual=100.0*(p1-bias)
decoded=' '
c=' '
if(qual.gt.1.0) then
if(qual.lt.6.0) c='?'
decoded=testmsg(ip1)
else
qual=0.
endif
decoded(22:22)=c
! Make sure everything is upper case.
do i=1,22
if(decoded(i:i).ge.'a' .and. decoded(i:i).le.'z') &
decoded(i:i)=char(ichar(decoded(i:i))-32)
enddo
call timer('deep65b ',1)
return
end subroutine deep65

77
q65w/libm65/demod64a.f90
View File

@ -1,77 +0,0 @@
subroutine demod64a(s3,nadd,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow)
! Demodulate the 64-bin spectra for each of 63 symbols in a frame.
! Parameters
! nadd number of spectra already summed
! mrsym most reliable symbol value
! mr2sym second most likely symbol value
! mrprob probability that mrsym was the transmitted value
! mr2prob probability that mr2sym was the transmitted value
implicit real*8 (a-h,o-z)
real*4 s3(64,63)
real*8 fs(64)
integer mrsym(63),mrprob(63),mr2sym(63),mr2prob(63)
common/mrscom/ mrs(63),mrs2(63)
afac=1.1 * float(nadd)**0.64
scale=255.999
! Compute average spectral value
sum=0.
do j=1,63
do i=1,64
sum=sum+s3(i,j)
enddo
enddo
ave=sum/(64.*63.)
i1=1 !Silence warning
i2=1
! Compute probabilities for most reliable symbol values
do j=1,63
s1=-1.e30
fsum=0.
psum=0.
do i=1,64
x=min(afac*s3(i,j)/ave,50.d0)
fs(i)=exp(x)
fsum=fsum+fs(i)
psum=psum + s3(i,j)
if(s3(i,j).gt.s1) then
s1=s3(i,j)
i1=i !Most reliable
endif
enddo
s2=-1.e30
do i=1,64
if(i.ne.i1 .and. s3(i,j).gt.s2) then
s2=s3(i,j)
i2=i !Second most reliable
endif
enddo
! p1=fs(i1)/fsum !Normalized probabilities
! p2=fs(i2)/fsum
p1=s1/psum
p2=s2/psum
mrsym(j)=i1-1
mr2sym(j)=i2-1
mrprob(j)=scale*p1
mr2prob(j)=scale*p2
mrs(j)=i1
mrs2(j)=i2
enddo
sum=0.
nlow=0
do j=1,63
sum=sum+mrprob(j)
if(mrprob(j).le.5) nlow=nlow+1
enddo
! ntest=sum/63
ntest=sum
return
end subroutine demod64a

183
q65w/libm65/display.f90
View File

@ -1,183 +0,0 @@
subroutine display(nkeep,ftol)
parameter (MAXLINES=400,MX=400,MAXCALLS=500)
integer indx(MAXLINES),indx2(MX)
character*83 line(MAXLINES),line2(MX),line3(MAXLINES)
character out*52,out0*52,cfreq0*3,livecq*58
character*6 callsign,callsign0
character*12 freqcall(MAXCALLS)
real freqkHz(MAXLINES)
integer utc(MAXLINES),utc2(MX),utcz
real*8 f0
save
out0=' '
rewind(26)
do i=1,MAXLINES
read(26,1010,end=10) line(i)
1010 format(a77)
read(line(i),1020) f0,ndf,nh,nm
1020 format(f8.3,i5,25x,i3,i2)
utc(i)=60*nh + nm
freqkHz(i)=1000.d0*(f0-144.d0) + 0.001d0*ndf
enddo
10 backspace(26)
nz=i-1
utcz=utc(nz)
nz=nz-1
if(nz.lt.1) go to 999
nquad=max(nkeep/4,3)
do i=1,nz
nage=utcz-utc(i)
if(nage.lt.0) nage=nage+1440
iage=nage/nquad
write(line(i)(73:74),1021) iage
1021 format(i2)
enddo
nage=utcz-utc(1)
if(nage.lt.0) nage=nage+1440
if(nage.gt.nkeep) then
do i=1,nz
nage=utcz-utc(i)
if(nage.lt.0) nage=nage+1440
if(nage.le.nkeep) go to 20
enddo
20 i0=i
nz=nz-i0+1
rewind(26)
if(nz.lt.1) go to 999
do i=1,nz
j=i+i0-1
line(i)=line(j)
utc(i)=utc(j)
freqkHz(i)=freqkHz(j)
write(26,1022) line(i)
1022 format(a77)
enddo
endif
call flush(26)
call indexx(freqkHz,nz,indx)
nstart=1
k3=0
k=1
m=indx(1)
if(m.lt.1 .or. m.gt.MAXLINES) then
print*,'Error in display.f90: ',nz,m
m=1
endif
line2(1)=line(m)
utc2(1)=utc(m)
do i=2,nz
j0=indx(i-1)
j=indx(i)
if(freqkHz(j)-freqkHz(j0).gt.2.0*ftol) then
if(nstart.eq.0) then
k=k+1
line2(k)=""
utc2(k)=-1
endif
kz=k
if(nstart.eq.1) then
call indexx(float(utc2(1:kz)),kz,indx2)
k3=0
do k=1,kz
k3=min(k3+1,400)
line3(k3)=line2(indx2(k))
enddo
nstart=0
else
call indexx(float(utc2(1:kz)),kz,indx2)
do k=1,kz
k3=min(k3+1,400)
line3(k3)=line2(indx2(k))
enddo
endif
k=0
endif
if(i.eq.nz) then
k=k+1
line2(k)=""
utc2(k)=-1
endif
k=k+1
line2(k)=line(j)
utc2(k)=utc(j)
j0=j
enddo
kz=k
call indexx(float(utc2(1:kz)),kz,indx2)
do k=1,kz
k3=min(k3+1,400)
line3(k3)=line2(indx2(k))
enddo
rewind 19
rewind 20
cfreq0=' '
nc=0
callsign0=' '
do k=1,k3
out=line3(k)(6:13)//line3(k)(28:31)//line3(k)(39:45)// &
line3(k)(35:38)//line3(k)(46:74)
if(out(1:3).ne.' ') then
cfreq0=out(1:3)
livecq=line3(k)(6:13)//line3(k)(28:31)//line3(k)(39:45)// &
line3(k)(23:27)//line3(k)(35:38)//line3(k)(46:70)// &
line3(k)(73:77)
if(livecq(56:56).eq.':') livecq(56:58)=' '//livecq(56:57)
if(index(livecq,' CQ ').gt.0 .or. index(livecq,' QRZ ').gt.0 .or. &
index(livecq,' QRT ').gt.0 .or. index(livecq,' CQV ').gt.0 .or. &
index(livecq,' CQH ').gt.0) write(19,1029) livecq
1029 format(a58)
! Suppress listing duplicate (same time, decoded message, and frequency)
if(out(14:17).ne.out0(14:17) .or. out(26:50).ne.out0(26:50) .or. &
out(1:3).ne.out0(1:3)) then
!###
! write(*,1030) out !Messages
!1030 format('@',a52)
!###
out0=out
endif
i1=index(out(26:),' ')
callsign=out(i1+26:)
i2=index(callsign,' ')
if(i2.gt.1) callsign(i2:)=' '
if(callsign.ne.' ' .and. callsign.ne.callsign0) then
len=i2-1
if(len.lt.0) len=6
if(len.ge.4) then !Omit short "callsigns"
if(nc.lt.MAXCALLS) nc=nc+1
freqcall(nc)=cfreq0//' '//callsign//line3(k)(73:74)
callsign0=callsign
endif
endif
if(callsign.ne.' ' .and. callsign.eq.callsign0) then
freqcall(nc)=cfreq0//' '//callsign//line3(k)(73:74)
endif
endif
enddo
flush(19)
if(nc.lt.MAXCALLS) nc=nc+1
freqcall(nc)=' '
if(nc.lt.MAXCALLS) nc=nc+1
freqcall(nc)=' '
freqcall(nc+1)=' '
freqcall(nc+2)=' '
!###
! do i=1,nc
! write(*,1042) freqcall(i) !Band Map
!1042 format('&',a12)
! enddo
!###
999 continue
return
end subroutine display

41
q65w/libm65/dpol.f90
View File

@ -1,41 +0,0 @@
real function dpol(mygrid,hisgrid)
! Compute spatial polartzation offset in degrees for the present
! time, between two specified grid locators.
character*6 MyGrid,HisGrid
real lat,lon,LST
character cdate*8,ctime2*10,czone*5
integer it(8)
data rad/57.2957795/
call date_and_time(cdate,ctime2,czone,it)
nyear=it(1)
month=it(2)
nday=it(3)
nh=it(5)-it(4)/60
nm=it(6)
ns=it(7)
uth=nh + nm/60.0 + ns/3600.0
call grid2deg(MyGrid,lon,lat)
call MoonDop(nyear,month,nday,uth,-lon,lat,RAMoon,DecMoon, &
LST,HA,AzMoon,ElMoon,vr,dist)
xx=sin(lat/rad)*cos(ElMoon/rad) - cos(lat/rad)* &
cos(AzMoon/rad)*sin(ElMoon/rad)
yy=cos(lat/rad)*sin(AzMoon/rad)
poloffset1=rad*atan2(yy,xx)
call grid2deg(hisGrid,lon,lat)
call MoonDop(nyear,month,nday,uth,-lon,lat,RAMoon,DecMoon, &
LST,HA,AzMoon,ElMoon,vr,dist)
xx=sin(lat/rad)*cos(ElMoon/rad) - cos(lat/rad)* &
cos(AzMoon/rad)*sin(ElMoon/rad)
yy=cos(lat/rad)*sin(AzMoon/rad)
poloffset2=rad*atan2(yy,xx)
dpol=mod(poloffset2-poloffset1+720.0,180.0)
if(dpol.gt.90.0) dpol=dpol-180.0
return
end function dpol

14
q65w/libm65/encode65.f90
View File

@ -1,14 +0,0 @@
subroutine encode65(message,sent)
use packjt
character message*22
integer dgen(12)
integer sent(63)
call packmsg(message,dgen,itype)
call rs_encode(dgen,sent)
call interleave63(sent,1)
call graycode(sent,63,1)
return
end subroutine encode65

136
q65w/libm65/extract.f90
View File

@ -1,136 +0,0 @@
subroutine extract(s3,nadd,ncount,nhist,decoded,ltext)
use packjt
use timer_module, only: timer
real s3(64,63)
character decoded*22
integer dat4(12)
integer mrsym(63),mr2sym(63),mrprob(63),mr2prob(63)
logical first,ltext
integer correct(63),itmp(63)
integer param(0:8)
integer h0(0:11),d0(0:11)
real r0(0:11)
common/test001/s3a(64,63),mrs(63),mrs2(63) !### TEST ONLY ###
! 0 1 2 3 4 5 6 7 8 9 10 11
data h0/41,42,43,43,44,45,46,47,48,48,49,49/
data d0/71,72,73,74,76,77,78,80,81,82,83,83/
! 0 1 2 3 4 5 6 7 8 9 10 11
data r0/0.70,0.72,0.74,0.76,0.78,0.80,0.82,0.84,0.86,0.88,0.90,0.90/
data first/.true./,nsec1/0/
save
nfail=0
call pctile(s3,4032,50,base) ! ### or, use ave from demod64a
s3=s3/base
s3a=s3
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,4032,50,base) ! ### or, use ave from demod64a
s3(ipk,1:63)=base
if(nfail.gt.30) then
decoded=' '
ncount=-1
go to 900
endif
go to 1
endif
mrs=mrsym
mrs2=mr2sym
call graycode(mrsym,63,-1)
call interleave63(mrsym,-1)
call interleave63(mrprob,-1)
call graycode(mr2sym,63,-1)
call interleave63(mr2sym,-1)
call interleave63(mr2prob,-1)
ntrials=10000
naggressive=10
ntry=0
param=0
call timer('ftrsd ',0)
call ftrsd2(mrsym,mrprob,mr2sym,mr2prob,ntrials,correct,param,ntry)
call timer('ftrsd ',1)
ncandidates=param(0)
nhard=param(1)
nsoft=param(2)
nerased=param(3)
rtt=0.001*param(4)
ntotal=param(5)
qual=0.001*param(7)
nd0=81
r00=0.87
if(naggressive.eq.10) then
nd0=83
r00=0.90
endif
if(ntotal.le.nd0 .and. rtt.le.r00) nft=1
n=naggressive
if(nhard.gt.50) nft=0
if(nhard.gt.h0(n)) nft=0
if(ntotal.gt.d0(n)) nft=0
if(rtt.gt.r0(n)) nft=0
ncount=-1
decoded=' '
ltext=.false.
if(nft.gt.0) then
! Turn the corrected symbol array into channel symbols for subtraction;
! pass it back to jt65a via common block "chansyms65".
do i=1,12
dat4(i)=correct(13-i)
enddo
do i=1,63
itmp(i)=correct(64-i)
enddo
correct(1:63)=itmp(1:63)
call interleave63(correct,1)
call graycode65(correct,63,1)
call unpackmsg(dat4,decoded) !Unpack the user message
ncount=0
if(iand(dat4(10),8).ne.0) ltext=.true.
endif
900 continue
if(nft.eq.1 .and. nhard.lt.0) decoded=' '
! write(81,3001) naggressive,ncandidates,nhard,ntotal,rtt,qual,decoded
!3001 format(i2,i6,i3,i4,2f8.2,2x,a22)
return
end subroutine extract
subroutine getpp(workdat,p)
integer workdat(63)
integer a(63)
common/test001/s3a(64,63),mrs(63),mrs2(63)
a(1:63)=workdat(63:1:-1)
call interleave63(a,1)
call graycode(a,63,1)
psum=0.
do j=1,63
i=a(j)+1
x=s3a(i,j)
s3a(i,j)=0.
psum=psum + x
s3a(i,j)=x
enddo
p=psum/63.0
return
end subroutine getpp

77
q65w/libm65/fchisq.f90
View File

@ -1,77 +0,0 @@
real function fchisq(cx,cy,npts,fsample,nflip,a,ccfmax,dtmax)
use timer_module, only: timer
parameter (NMAX=60*96000) !Samples per 60 s
complex cx(npts),cy(npts)
real a(5)
complex w,wstep,za,zb,z
real ss(3000)
complex csx(0:NMAX/64),csy(0:NMAX/64)
data twopi/6.283185307/a1,a2,a3/99.,99.,99./
save
call timer('fchisq ',0)
baud=11025.0/4096.0
nsps=nint(fsample/baud) !Samples per symbol
nsph=nsps/2 !Samples per half-symbol
ndiv=16 !Output ss() steps per symbol
nout=ndiv*npts/nsps
dtstep=1.0/(ndiv*baud) !Time per output step
if(a(1).ne.a1 .or. a(2).ne.a2 .or. a(3).ne.a3) then
a1=a(1)
a2=a(2)
a3=a(3)
! Mix and integrate the complex X and Y signals
csx(0)=0.
csy(0)=0.
w=1.0
x0=0.5*(npts+1)
s=2.0/npts
do i=1,npts
x=s*(i-x0)
if(mod(i,100).eq.1) then
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))
endif
w=w*wstep
csx(i)=csx(i-1) + w*cx(i)
csy(i)=csy(i-1) + w*cy(i)
enddo
endif
! Compute 1/2-symbol powers at 1/16-symbol steps.
fac=1.e-4
pol=a(4)/57.2957795
aa=cos(pol)
bb=sin(pol)
do i=1,nout
j=i*nsps/ndiv
k=j-nsph
ss(i)=0.
if(k.ge.1) then
za=csx(j)-csx(k)
zb=csy(j)-csy(k)
z=aa*za + bb*zb
ss(i)=fac*(real(z)**2 + aimag(z)**2)
endif
enddo
ccfmax=0.
call timer('ccf2 ',0)
call ccf2(ss,nout,nflip,ccf,lagpk)
call timer('ccf2 ',1)
if(ccf.gt.ccfmax) then
ccfmax=ccf
dtmax=lagpk*dtstep
endif
fchisq=-ccfmax
call timer('fchisq ',1)
return
end function fchisq

44
q65w/libm65/fil6521.f90
View File

@ -1,44 +0,0 @@
subroutine fil6521(c1,n1,c2,n2)
! FIR lowpass filter designed using ScopeFIR
! Pass #1 Pass #2
!-----------------------------------------------
! fsample (Hz) 1378.125 Input sample rate
! Ntaps 21 Number of filter taps
! fc (Hz) 40 Cutoff frequency
! fstop (Hz) 172.266 Lower limit of stopband
! Ripple (dB) 0.1 Ripple in passband
! Stop Atten (dB) 38 Stopband attenuation
! fout (Hz) 344.531 Output sample rate
parameter (NTAPS=21)
parameter (NH=(NTAPS-1)/2)
parameter (NDOWN=4) !Downsample ratio = 1/4
complex c1(n1)
complex c2(n1/NDOWN)
! Filter coefficients:
real a(-NH:NH)
data a/ &
-0.011958606980,-0.013888627387,-0.015601306443,-0.010602249570, &
0.003804023436, 0.028320058273, 0.060903935217, 0.096841904411, &
0.129639871228, 0.152644580853, 0.160917511283, 0.152644580853, &
0.129639871228, 0.096841904411, 0.060903935217, 0.028320058273, &
0.003804023436,-0.010602249570,-0.015601306443,-0.013888627387, &
-0.011958606980/
n2=(n1-NTAPS+NDOWN)/NDOWN
k0=NH-NDOWN+1
! Loop over all output samples
do i=1,n2
c2(i)=0.
k=k0 + NDOWN*i
do j=-NH,NH
c2(i)=c2(i) + c1(j+k)*a(j)
enddo
enddo
return
end subroutine fil6521

21
q65w/libm65/fmtmsg.f90
View File

@ -1,21 +0,0 @@
subroutine fmtmsg(msg,iz)
character*22 msg
! Convert all letters to upper case
iz=22
do i=1,22
if(msg(i:i).ge.'a' .and. msg(i:i).le.'z') &
msg(i:i)= char(ichar(msg(i:i))+ichar('A')-ichar('a'))
if(msg(i:i).ne.' ') iz=i
enddo
do iter=1,5 !Collapse multiple blanks into one
ib2=index(msg(1:iz),' ')
if(ib2.lt.1) go to 100
msg=msg(1:ib2)//msg(ib2+2:)
iz=iz-1
enddo
100 return
end subroutine fmtmsg

99
q65w/libm65/gen65.f90
View File

@ -1,99 +0,0 @@
subroutine gen65(message,mode65,samfac,nsendingsh,msgsent,iwave,nwave)
! Encodes a JT65 message into a wavefile.
! Executes in 17 ms on opti-745.
use packjt
parameter (NMAX=2*60*11025) !Max length of wave file
character*22 message !Message to be generated
character*22 msgsent !Message as it will be received
character*3 cok !' ' or 'OOO'
real*8 dt,phi,f,f0,dfgen,dphi,twopi,samfac
integer*2 iwave(NMAX) !Generated wave file
integer dgen(12)
integer sent(63)
logical first
integer nprc(126)
real pr(126)
data nprc/1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0, &
0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1, &
0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1, &
0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1, &
1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1, &
0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1, &
1,1,1,1,1,1/
data twopi/6.283185307179586476d0/,first/.true./
save
if(first) then
do i=1,126
pr(i)=2*nprc(i)-1
enddo
first=.false.
endif
call chkmsg(message,cok,nspecial,flip)
if(nspecial.eq.0) then
call packmsg(message,dgen,itype) !Pack message into 72 bits
nsendingsh=0
if(iand(dgen(10),8).ne.0) nsendingsh=-1 !Plain text flag
call rs_encode(dgen,sent)
call interleave63(sent,1) !Apply interleaving
call graycode(sent,63,1) !Apply Gray code
nsym=126 !Symbols per transmission
nsps=4096
else
nsym=32
nsps=16384
nsendingsh=1 !Flag for shorthand message
endif
if(mode65.eq.0) go to 900
! Set up necessary constants
dt=1.d0/(samfac*11025.d0)
f0=118*11025.d0/1024
dfgen=mode65*11025.d0/4096.d0
phi=0.d0
dphi=twopi*dt*f0
i=0
k=0
do j=1,nsym
if(message(1:5).ne.'@TUNE') then
f=f0
if(nspecial.ne.0 .and. mod(j,2).eq.0) f=f0+10*nspecial*dfgen
if(nspecial.eq.0 .and. flip*pr(j).lt.0.0) then
k=k+1
f=f0+(sent(k)+2)*dfgen
endif
dphi=twopi*dt*f
endif
do ii=1,nsps
phi=phi+dphi
if(phi.gt.twopi) phi=phi-twopi
xphi=phi
i=i+1
iwave(2*i-1)=32767.0*cos(xphi)
iwave(2*i)=32767.0*sin(xphi)
enddo
enddo
iwave(2*nsym*nsps+1:)=0
nwave=2*nsym*nsps + 5512
call unpackmsg(dgen,msgsent)
if(flip.lt.0.0) then
do i=22,1,-1
if(msgsent(i:i).ne.' ') goto 10
enddo
10 msgsent=msgsent(1:i)//' OOO'
endif
if(nsendingsh.eq.1) then
if(nspecial.eq.2) msgsent='RO'
if(nspecial.eq.3) msgsent='RRR'
if(nspecial.eq.4) msgsent='73'
endif
900 return
end subroutine gen65

52
q65w/libm65/gen_q65_cwave.f90
View File

@ -1,52 +0,0 @@
subroutine gen_q65_cwave(msg,ntxfreq,ntone_spacing,msgsent,cwave,nwave)
! Encodes a Q65 message to yield complex cwave() at fsample = 96000 Hz
use packjt
use q65_encoding
parameter (NMAX=60*96000)
character*22 msg
character*22 msgsent !Message as it will be received
character*37 msg37
real*8 t,dt,phi,f,f0,dfgen,dphi,twopi,tsym
complex cwave(NMAX)
integer codeword(65),itone(85)
integer icos7(0:6)
data icos7/2,5,6,0,4,1,3/ !Defines a 7x7 Costas array
data twopi/6.283185307179586476d0/
save
msgsent=msg
msg37=''
msg37(1:22)=msg
call get_q65_tones(msg37,codeword,itone)
! Set up necessary constants
nsym=85
tsym=7200.d0/12000.d0
dt=1.d0/96000.d0
f0=ntxfreq
dfgen=ntone_spacing*12000.d0/7200.d0
phi=0.d0
dphi=twopi*dt*f0
i=0
nwave=85*7200*96000.d0/12000.d0
t=0.d0
j0=0
do i=1,nwave
t=t+dt
j=t/tsym + 1
if(j.gt.85) exit
if(j.ne.j0) then
f=f0 + itone(j)*dfgen
dphi=twopi*dt*f
j0=j
endif
phi=phi+dphi
if(phi.gt.twopi) phi=phi-twopi
xphi=phi
cwave(i)=cmplx(cos(xphi),-sin(xphi))
enddo
999 return
end subroutine gen_q65_cwave

54
q65w/libm65/gen_q65_wave.f90
View File

@ -1,54 +0,0 @@
subroutine gen_q65_wave(msg,ntxfreq,mode65,msgsent,iwave,nwave)
! Encodes a Q65 message to yield complex iwave() at fsample = 11025 Hz
use packjt
use q65_encoding
parameter (NMAX=2*60*11025)
character*22 msg
character*22 msgsent !Message as it will be received
character*37 msg37
real*8 t,dt,phi,f,f0,dfgen,dphi,twopi,tsym
integer codeword(65),itone(85)
integer*2 iwave(NMAX)
integer icos7(0:6)
data icos7/2,5,6,0,4,1,3/ !Defines a 7x7 Costas array
data twopi/6.283185307179586476d0/
save
msgsent=msg
msg37=''
msg37(1:22)=msg
call get_q65_tones(msg37,codeword,itone)
! Set up necessary constants
nsym=85
tsym=7200.d0/12000.d0
dt=1.d0/11025.d0
f0=ntxfreq
ndf=2**(mode65-1)
dfgen=ndf*12000.d0/7200.d0
phi=0.d0
dphi=twopi*dt*f0
i=0
iz=85*7200*11025.d0/12000.d0
t=0.d0
j0=0
do i=1,iz
t=t+dt
j=t/tsym + 1.0
if(j.ne.j0) then
f=f0 + itone(j)*dfgen
dphi=twopi*dt*f
j0=j
endif
phi=phi+dphi
if(phi.gt.twopi) phi=phi-twopi
xphi=phi
iwave(2*i-1)=32767.0*cos(xphi)
iwave(2*i)=32767.0*sin(xphi)
enddo
nwave=2*iz
999 return
end subroutine gen_q65_wave

9
q65w/libm65/graycode65.f90
View File

@ -1,9 +0,0 @@
subroutine graycode65(dat,n,idir)
integer dat(n)
do i=1,n
dat(i)=igray(dat(i),idir)
enddo
return
end subroutine graycode65

47
q65w/libm65/jt65code.f90
View File

@ -1,47 +0,0 @@
program JT65code
! Provides examples of message packing, bit and symbol ordering,
! Reed Solomon encoding, and other necessary details of the JT65
! protocol.
character*22 msg0,msg,decoded,cok*3
integer dgen(12),sent(63),recd(12),era(51)
nargs=iargc()
if(nargs.ne.1) then
print*,'Usage: JT65code "message"'
go to 999
endif
call getarg(1,msg0) !Get message from command line
msg=msg0
call chkmsg(msg,cok,nspecial,flip) !See if it includes "OOO" report
if(nspecial.gt.0) then !or is a shorthand message
write(*,1010)
1010 format('Shorthand message.')
go to 999
endif
call packmsg(msg,dgen) !Pack message into 72 bits
write(*,1020) msg0
1020 format('Message: ',a22) !Echo input message
if(iand(dgen(10),8).ne.0) write(*,1030) !Is plain text bit set?
1030 format('Plain text.')
write(*,1040) dgen
1040 format('Packed message, 6-bit symbols: ',12i3) !Display packed symbols
call rs_encode(dgen,sent) !RS encode
call interleave63(sent,1) !Interleave channel symbols
call graycode(sent,63,1) !Apply Gray code
write(*,1050) sent
1050 format('Channel symbols, including FEC:'/(i5,20i3))
call graycode(sent,63,-1)
call interleave63(sent,-1)
call rs_decode(sent,era,0,recd,nerr)
call unpackmsg(recd,decoded) !Unpack the user message
write(*,1060) decoded,cok
1060 format('Decoded message: ',a22,2x,a3)
999 end program JT65code

2
q65w/libm65/map65a.f90
View File

@ -211,8 +211,6 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, &
if(npol.lt.0) npol=npol+180
endif
call txpol(xpol,decoded,mygrid,npol,nxant,ntxpol,cp)
cmode='#A'
if(mode65.eq.2) cmode='#B'
if(mode65.eq.4) cmode='#C'

229
q65w/libm65/mapsim.f90
View File

@ -1,229 +0,0 @@
program mapsim
! Generate simulated data for testing of MAP65
parameter (NMAX=60*96000)
real*4 d4(4,NMAX) !Floating-point data
integer*2 id4(4,NMAX) !i*2 data, dual polarization
integer*2 id2(2,NMAX) !i*2 data, single polarization
complex cwave(NMAX) !Generated complex waveform (no noise)
complex z,zx,zy
real*8 fcenter,fsample,samfac,f,dt,twopi,phi,dphi
logical bq65
character msg0*22,message*22,msgsent*22,arg*8,fname*11,mode*2
character*16 msg_list(60)
data msg_list/ &
'W1AAA K2BBB EM00','W2CCC K3DDD EM01','W3EEE K4FFF EM02', &
'W5GGG K6HHH EM03','W7III K8JJJ EM04','W9KKK K0LLL EM05', &
'G0MMM F1NNN JN06','G2OOO F3PPP JN07','G4QQQ F5RRR JN08', &
'G6SSS F7TTT JN09','W1XAA K2XBB EM10','W2XCC K3XDD EM11', &
'W3XEE K4XFF EM12','W5XGG K6XHH EM13','W7XII K8XJJ EM14', &
'W9XKK K0XLL EM15','G0XMM F1XNN JN16','G2XOO F3XPP JN17', &
'G4XQQ F5XRR JN18','G6XSS F7XTT JN19','W1YAA K2YBB EM20', &
'W2YCC K3YDD EM21','W3YEE K4YFF EM22','W5YGG K6YHH EM23', &
'W7YII K8YJJ EM24','W9YKK K0YLL EM25','G0YMM F1YNN JN26', &
'G2YOO F3YPP JN27','G4YQQ F5YRR JN28','G6YSS F7YTT JN29', &
'W1ZAA K2ZBB EM30','W2ZCC K3ZDD EM31','W3ZEE K4ZFF EM32', &
'W5ZGG K6ZHH EM33','W7ZII K8ZJJ EM34','W9ZKK K0ZLL EM35', &
'G0ZMM F1ZNN JN36','G2ZOO F3ZPP JN37','G4ZQQ F5ZRR JN38', &
'G6ZSS F7ZTT JN39','W1AXA K2BXB EM40','W2CXC K3DXD EM41', &
'W3EXE K4FXF EM42','W5GXG K6HXH EM43','W7IXI K8JXJ EM44', &
'W9KXK K0LXL EM45','G0MXM F1NXN JN46','G2OXO F3PXP JN47', &
'G4QXQ F5RXR JN48','G6SXS F7TXT JN49','W1AYA K2BYB EM50', &
'W2CYC K3DYD EM51','W3EYE K4FYF EM52','W5GYG K6HYH EM53', &
'W7IYI K8JYJ EM54','W9KYK K0LYL EM55','G0MYM F1NYN JN56', &
'G2OYO F3PYP JN57','G4QYQ F5RYR JN58','G6SYS F7TYT JN59'/
nargs=iargc()
if(nargs.ne.10) then
print*,'Usage: mapsim "message" mode DT fa fb nsigs pol fDop SNR nfiles'
print*,'Example: mapsim "CQ K1ABC FN42" B 2.5 -20 20 21 45 0.0 -20 1'
print*,' '
print*,' mode = A B C for JT65; QA-QE for Q65-60A'
print*,' fa = lowest freq in kHz, relative to center'
print*,' fb = highest freq in kHz, relative to center'
print*,' message = "list" to use callsigns from list'
print*,' pol = -1 to generate a range of polarization angles.'
print*,' SNR = 0 to generate a range of SNRs.'
go to 999
endif
call getarg(1,msg0)
call getarg(2,mode) !JT65 sub-mode (A B C QA-QE)
call getarg(3,arg)
read(arg,*) dt0 !Time delay
call getarg(4,arg)
read(arg,*) fa !Lowest freq (kHz, relative to fcenter)
call getarg(5,arg)
read(arg,*) fb !Highest freq
call getarg(6,arg)
read(arg,*) nsigs !Number of signals in each file
call getarg(7,arg)
read(arg,*) npol !Polarization in degrees
pol=npol
call getarg(8,arg)
read(arg,*) fdop !Doppler spread
call getarg(9,arg)
read(arg,*) snrdb !S/N
call getarg(10,arg)
read(arg,*) nfiles !Number of files
message=msg0 !Transmitted message
rmsdb=25.
rms=10.0**(0.05*rmsdb)
fcenter=144.125d0 !Center frequency (MHz)
fsample=96000.d0 !Sample rate (Hz)
dt=1.d0/fsample !Sample interval (s)
twopi=8.d0*atan(1.d0)
rad=360.0/twopi
samfac=1.d0
bq65=(mode(1:1).eq.'Q')
ntone_spacing=1
ntxfreq=1270
fac=1.0/32767.0
if(mode(1:1).eq.'B' .or. mode(2:2).eq.'B') ntone_spacing=2
if(mode(1:1).eq.'C' .or. mode(2:2).eq.'C') ntone_spacing=4
if(mode(2:2).eq.'D') ntone_spacing=8
if(mode(2:2).eq.'E') ntone_spacing=16
npts=NMAX
write(*,1000)
1000 format('File N Mode DT freq pol fDop SNR Message'/68('-'))
do ifile=1,nfiles
ilist=0
nmin=ifile-1
if(mode(2:2).eq.' ') nmin=2*nmin
write(fname,1002) nmin !Create the output filenames
1002 format('000000_',i4.4)
open(10,file=fname//'.iq',access='stream',status='unknown')
open(11,file=fname//'.tf2',access='stream',status='unknown')
call noisegen(d4,npts) !Generate Gaussuian noise
if(msg0(1:4).ne.'list') then
if(bq65) then
call gen_q65_cwave(message,ntxfreq,ntone_spacing,msgsent, &
cwave,nwave)
else
call cgen65(message,ntone_spacing,samfac,nsendingsh,msgsent, &
cwave,nwave)
endif
endif
if(fdop.gt.0.0) call dopspread(cwave,fdop)
do isig=1,nsigs
if(msg0(1:4).eq.'list') then
ilist=ilist+1
message=msg_list(ilist)
if(bq65) then
call gen_q65_cwave(message,ntxfreq,ntone_spacing,msgsent, &
cwave,nwave)
else
call cgen65(message,ntone_spacing,samfac,nsendingsh,msgsent, &
cwave,nwave)
endif
endif
if(npol.lt.0) pol=(isig-1)*180.0/nsigs
a=cos(pol/rad)
b=sin(pol/rad)
f=1000.0*(fa+fb)/2.0
if(nsigs.gt.1) f=1000.0*(fa + (isig-1)*(fb-fa)/(nsigs-1.0))
dphi=twopi*f*dt + 0.5*twopi
snrdbx=snrdb
if(snrdb.eq.0.0) snrdbx=-15.0 - 15.0*(isig-1.0)/nsigs
sig=sqrt(2.2*2500.0/96000.0) * 10.0**(0.05*snrdbx)
write(*,1020) ifile,isig,mode,dt0,0.001*f,nint(pol),fDop,snrdbx,msgsent
1020 format(i3,i3,2x,a2,f6.2,f8.3,i5,2f7.1,2x,a22)
phi=0.
! i0=fsample*(3.5d0+0.05d0*(isig-1))
i0=fsample*(1.d0 + dt0)
do i=1,nwave
phi=phi + dphi
if(phi.lt.-twopi) phi=phi+twopi
if(phi.gt.twopi) phi=phi-twopi
xphi=phi
z=sig*cwave(i)*cmplx(cos(xphi),-sin(xphi))
zx=a*z
zy=b*z
j=i+i0
d4(1,j)=d4(1,j) + real(zx)
d4(2,j)=d4(2,j) + aimag(zx)
d4(3,j)=d4(3,j) + real(zy)
d4(4,j)=d4(4,j) + aimag(zy)
enddo
enddo
do i=1,npts
id4(1,i)=nint(rms*d4(1,i))
id4(2,i)=nint(rms*d4(2,i))
id4(3,i)=nint(rms*d4(3,i))
id4(4,i)=nint(rms*d4(4,i))
id2(1,i)=id4(1,i)
id2(2,i)=id4(2,i)
enddo
write(10) fcenter,id2(1:2,1:npts)
write(11) fcenter,id4(1:4,1:npts)
close(10)
close(11)
enddo
999 end program mapsim
subroutine dopspread(cwave,fspread)
parameter (NMAX=60*96000)
parameter (NFFT=NMAX,NH=NFFT/2)
complex cwave(NMAX)
complex cspread(0:NMAX-1)
twopi=8.0*atan(1.0)
df=96000.0/nfft
cspread(0)=1.0
cspread(NH)=0.
b=6.0 !Use truncated Lorenzian shape for fspread
do i=1,NH
f=i*df
x=b*f/fspread
z=0.
a=0.
if(x.lt.3.0) then !Cutoff beyond x=3
a=sqrt(1.111/(1.0+x*x)-0.1) !Lorentzian amplitude
phi1=twopi*rran() !Random phase
z=a*cmplx(cos(phi1),sin(phi1))
endif
cspread(i)=z
z=0.
if(x.lt.3.0) then !Same thing for negative freqs
phi2=twopi*rran()
z=a*cmplx(cos(phi2),sin(phi2))
endif
cspread(nfft-i)=z
enddo
call four2a(cspread,nfft,1,1,1) !Transform to time domain
sum=0.
do i=0,nfft-1
p=real(cspread(i))**2 + aimag(cspread(i))**2
sum=sum+p
enddo
avep=sum/nfft
fac=sqrt(1.0/avep)
cspread=fac*cspread !Normalize to constant avg power
cwave=cspread*cwave !Apply Rayleigh fading
! do i=0,nfft-1
! p=real(cspread(i))**2 + aimag(cspread(i))**2
! write(14,3010) i,p,cspread(i)
!3010 format(i8,3f12.6)
! enddo
return
end subroutine dopspread

13
q65w/libm65/noisegen.f90
View File

@ -1,13 +0,0 @@
subroutine noisegen(d4,nmax)
real*4 d4(4,nmax)
do i=1,nmax
d4(1,i)=gran()
d4(2,i)=gran()
d4(3,i)=gran()
d4(4,i)=gran()
enddo
return
end subroutine noisegen

996
q65w/libm65/packjt.f90
View File

@ -1,996 +0,0 @@
module packjt
contains
subroutine packbits(dbits,nsymd,m0,sym)
! Pack 0s and 1s from dbits() into sym() with m0 bits per word.
! NB: nsymd is the number of packed output words.
integer sym(:)
integer*1 dbits(:)
k=0
do i=1,nsymd
n=0
do j=1,m0
k=k+1
m=dbits(k)
n=ior(ishft(n,1),m)
enddo
sym(i)=n
enddo
return
end subroutine packbits
subroutine unpackbits(sym,nsymd,m0,dbits)
! Unpack bits from sym() into dbits(), one bit per byte.
! NB: nsymd is the number of input words, and m0 their length.
! there will be m0*nsymd output bytes, each 0 or 1.
integer sym(:)
integer*1 dbits(:)
k=0
do i=1,nsymd
mask=ishft(1,m0-1)
do j=1,m0
k=k+1
dbits(k)=0
if(iand(mask,sym(i)).ne.0) dbits(k)=1
mask=ishft(mask,-1)
enddo
enddo
return
end subroutine unpackbits
subroutine packcall(callsign,ncall,text)
! Pack a valid callsign into a 28-bit integer.
parameter (NBASE=37*36*10*27*27*27)
character callsign*6,c*1,tmp*6
logical text
text=.false.
! Work-around for Swaziland prefix:
if(callsign(1:4).eq.'3DA0') callsign='3D0'//callsign(5:6)
if(callsign(1:3).eq.'CQ ') then
ncall=NBASE + 1
if(callsign(4:4).ge.'0' .and. callsign(4:4).le.'9' .and. &
callsign(5:5).ge.'0' .and. callsign(5:5).le.'9' .and. &
callsign(6:6).ge.'0' .and. callsign(6:6).le.'9') then
read(callsign(4:6),*) nfreq
ncall=NBASE + 3 + nfreq
endif
return
else if(callsign(1:4).eq.'QRZ ') then
ncall=NBASE + 2
return
else if(callsign(1:3).eq.'DE ') then
ncall=267796945
return
endif
tmp=' '
if(callsign(3:3).ge.'0' .and. callsign(3:3).le.'9') then
tmp=callsign
else if(callsign(2:2).ge.'0' .and. callsign(2:2).le.'9') then
if(callsign(6:6).ne.' ') then
text=.true.
return
endif
tmp=' '//callsign(:5)
else
text=.true.
return
endif
do i=1,6
c=tmp(i:i)
if(c.ge.'a' .and. c.le.'z') &
tmp(i:i)=char(ichar(c)-ichar('a')+ichar('A'))
enddo
n1=0
if((tmp(1:1).ge.'A'.and.tmp(1:1).le.'Z').or.tmp(1:1).eq.' ') n1=1
if(tmp(1:1).ge.'0' .and. tmp(1:1).le.'9') n1=1
n2=0
if(tmp(2:2).ge.'A' .and. tmp(2:2).le.'Z') n2=1
if(tmp(2:2).ge.'0' .and. tmp(2:2).le.'9') n2=1
n3=0
if(tmp(3:3).ge.'0' .and. tmp(3:3).le.'9') n3=1
n4=0
if((tmp(4:4).ge.'A'.and.tmp(4:4).le.'Z').or.tmp(4:4).eq.' ') n4=1
n5=0
if((tmp(5:5).ge.'A'.and.tmp(5:5).le.'Z').or.tmp(5:5).eq.' ') n5=1
n6=0
if((tmp(6:6).ge.'A'.and.tmp(6:6).le.'Z').or.tmp(6:6).eq.' ') n6=1
if(n1+n2+n3+n4+n5+n6 .ne. 6) then
text=.true.
return
endif
ncall=nchar(tmp(1:1))
ncall=36*ncall+nchar(tmp(2:2))
ncall=10*ncall+nchar(tmp(3:3))
ncall=27*ncall+nchar(tmp(4:4))-10
ncall=27*ncall+nchar(tmp(5:5))-10
ncall=27*ncall+nchar(tmp(6:6))-10
return
end subroutine packcall
subroutine unpackcall(ncall,word,iv2,psfx)
parameter (NBASE=37*36*10*27*27*27)
character word*12,c*37,psfx*4
data c/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ '/
word='......'
psfx=' '
n=ncall
iv2=0
if(n.ge.262177560) go to 20
word='......'
! if(n.ge.262177560) go to 999 !Plain text message ...
i=mod(n,27)+11
word(6:6)=c(i:i)
n=n/27
i=mod(n,27)+11
word(5:5)=c(i:i)
n=n/27
i=mod(n,27)+11
word(4:4)=c(i:i)
n=n/27
i=mod(n,10)+1
word(3:3)=c(i:i)
n=n/10
i=mod(n,36)+1
word(2:2)=c(i:i)
n=n/36
i=n+1
word(1:1)=c(i:i)
do i=1,4
if(word(i:i).ne.' ') go to 10
enddo
go to 999
10 word=word(i:)
go to 999
20 if(n.ge.267796946) go to 999
! We have a JT65v2 message
if((n.ge.262178563) .and. (n.le.264002071)) then
! CQ with prefix
iv2=1
n=n-262178563
i=mod(n,37)+1
psfx(4:4)=c(i:i)
n=n/37
i=mod(n,37)+1
psfx(3:3)=c(i:i)
n=n/37
i=mod(n,37)+1
psfx(2:2)=c(i:i)
n=n/37
i=n+1
psfx(1:1)=c(i:i)
else if((n.ge.264002072) .and. (n.le.265825580)) then
! QRZ with prefix
iv2=2
n=n-264002072
i=mod(n,37)+1
psfx(4:4)=c(i:i)
n=n/37
i=mod(n,37)+1
psfx(3:3)=c(i:i)
n=n/37
i=mod(n,37)+1
psfx(2:2)=c(i:i)
n=n/37
i=n+1
psfx(1:1)=c(i:i)
else if((n.ge.265825581) .and. (n.le.267649089)) then
! DE with prefix
iv2=3
n=n-265825581
i=mod(n,37)+1
psfx(4:4)=c(i:i)
n=n/37
i=mod(n,37)+1
psfx(3:3)=c(i:i)
n=n/37
i=mod(n,37)+1
psfx(2:2)=c(i:i)
n=n/37
i=n+1
psfx(1:1)=c(i:i)
else if((n.ge.267649090) .and. (n.le.267698374)) then
! CQ with suffix
iv2=4
n=n-267649090
i=mod(n,37)+1
psfx(3:3)=c(i:i)
n=n/37
i=mod(n,37)+1
psfx(2:2)=c(i:i)
n=n/37
i=n+1
psfx(1:1)=c(i:i)
else if((n.ge.267698375) .and. (n.le.267747659)) then
! QRZ with suffix
iv2=5
n=n-267698375
i=mod(n,37)+1
psfx(3:3)=c(i:i)
n=n/37
i=mod(n,37)+1
psfx(2:2)=c(i:i)
n=n/37
i=n+1
psfx(1:1)=c(i:i)
else if((n.ge.267747660) .and. (n.le.267796944)) then
! DE with suffix
iv2=6
n=n-267747660
i=mod(n,37)+1
psfx(3:3)=c(i:i)
n=n/37
i=mod(n,37)+1
psfx(2:2)=c(i:i)
n=n/37
i=n+1
psfx(1:1)=c(i:i)
else if(n.eq.267796945) then
! DE with no prefix or suffix
iv2=7
psfx = ' '
endif
999 if(word(1:3).eq.'3D0') word='3DA0'//word(4:)
return
end subroutine unpackcall
subroutine packgrid(grid,ng,text)
parameter (NGBASE=180*180)
character*4 grid
character*1 c1
logical text
text=.false.
if(grid.eq.' ') go to 90 !Blank grid is OK
! First, handle signal reports in the original range, -01 to -30 dB
if(grid(1:1).eq.'-') then
read(grid(2:3),*,err=800,end=800) n
if(n.ge.1 .and. n.le.30) then
ng=NGBASE+1+n
go to 900
endif
go to 10
else if(grid(1:2).eq.'R-') then
read(grid(3:4),*,err=800,end=800) n
if(n.ge.1 .and. n.le.30) then
ng=NGBASE+31+n
go to 900
endif
go to 10
! Now check for RO, RRR, or 73 in the message field normally used for grid
else if(grid(1:4).eq.'RO ') then
ng=NGBASE+62
go to 900
else if(grid(1:4).eq.'RRR ') then
ng=NGBASE+63
go to 900
else if(grid(1:4).eq.'73 ') then
ng=NGBASE+64
go to 900
endif
! Now check for extended-range signal reports: -50 to -31, and 0 to +49.
10 n=99
c1=grid(1:1)
read(grid,*,err=20,end=20) n
go to 30
20 read(grid(2:4),*,err=30,end=30) n
30 if(n.ge.-50 .and. n.le.49) then
if(c1.eq.'R') then
write(grid,1002) n+50
1002 format('LA',i2.2)
else
write(grid,1003) n+50
1003 format('KA',i2.2)
endif
go to 40
endif
! Maybe it's free text ?
if(grid(1:1).lt.'A' .or. grid(1:1).gt.'R') text=.true.
if(grid(2:2).lt.'A' .or. grid(2:2).gt.'R') text=.true.
if(grid(3:3).lt.'0' .or. grid(3:3).gt.'9') text=.true.
if(grid(4:4).lt.'0' .or. grid(4:4).gt.'9') text=.true.
if(text) go to 900
! OK, we have a properly formatted grid locator
40 call grid2deg(grid//'mm',dlong,dlat)
long=int(dlong)
lat=int(dlat+ 90.0)
ng=((long+180)/2)*180 + lat
go to 900
90 ng=NGBASE + 1
go to 900
800 text=.true.
900 continue
return
end subroutine packgrid
subroutine unpackgrid(ng,grid)
parameter (NGBASE=180*180)
character grid*4,grid6*6
grid=' '
if(ng.ge.32400) go to 10
dlat=mod(ng,180)-90
dlong=(ng/180)*2 - 180 + 2
call deg2grid(dlong,dlat,grid6)
grid=grid6(:4)
if(grid(1:2).eq.'KA') then
read(grid(3:4),*) n
n=n-50
write(grid,1001) n
1001 format(i3.2)
if(grid(1:1).eq.' ') grid(1:1)='+'
else if(grid(1:2).eq.'LA') then
read(grid(3:4),*) n
n=n-50
write(grid,1002) n
1002 format('R',i3.2)
if(grid(2:2).eq.' ') grid(2:2)='+'
endif
go to 900
10 n=ng-NGBASE-1
if(n.ge.1 .and.n.le.30) then
write(grid,1012) -n
1012 format(i3.2)
else if(n.ge.31 .and.n.le.60) then
n=n-30
write(grid,1022) -n
1022 format('R',i3.2)
else if(n.eq.61) then
grid='RO'
else if(n.eq.62) then
grid='RRR'
else if(n.eq.63) then
grid='73'
endif
900 return
end subroutine unpackgrid
subroutine packmsg(msg0,dat,itype)
! Packs a JT4/JT9/JT65 message into twelve 6-bit symbols
! itype Message Type
!--------------------
! 1 Standardd message
! 2 Type 1 prefix
! 3 Type 1 suffix
! 4 Type 2 prefix
! 5 Type 2 suffix
! 6 Free text
! -1 Does not decode correctly
parameter (NBASE=37*36*10*27*27*27)
parameter (NBASE2=262178562)
character*22 msg0,msg
integer dat(:)
character*12 c1,c2
character*4 c3
character*6 grid6
logical text1,text2,text3
msg=msg0
itype=1
call fmtmsg(msg,iz)
if(msg(1:6).eq.'CQ DX ') msg(3:3)='9'
if(msg(1:3).eq."CQ " .and. &
msg(4:4).ge.'A' .and. msg(4:4).le.'Z' .and. &
msg(5:5).ge.'A' .and. msg(5:5).le.'Z' .and. &
msg(6:6).eq.' ') msg='E9'//msg(4:)
! See if it's a CQ message
if(msg(1:3).eq.'CQ ') then
i=3
! ... and if so, does it have a reply frequency?
if(msg(4:4).ge.'0' .and. msg(4:4).le.'9' .and. &
msg(5:5).ge.'0' .and. msg(5:5).le.'9' .and. &
msg(6:6).ge.'0' .and. msg(6:6).le.'9') i=7
go to 1
endif
do i=1,22
if(msg(i:i).eq.' ') go to 1 !Get 1st blank
enddo
go to 10 !Consider msg as plain text
1 ia=i
c1=msg(1:ia-1)
do i=ia+1,22
if(msg(i:i).eq.' ') go to 2 !Get 2nd blank
enddo
go to 10 !Consider msg as plain text
2 ib=i
c2=msg(ia+1:ib-1)
do i=ib+1,22
if(msg(i:i).eq.' ') go to 3 !Get 3rd blank
enddo
go to 10 !Consider msg as plain text
3 ic=i
c3=' '
if(ic.ge.ib+1) c3=msg(ib+1:ic)
if(c3.eq.'OOO ') c3=' ' !Strip out the OOO flag
call getpfx1(c1,k1,nv2a)
if(nv2a.ge.4) go to 10
call packcall(c1,nc1,text1)
if(text1) go to 10
call getpfx1(c2,k2,nv2b)
call packcall(c2,nc2,text2)
if(text2) go to 10
if(nv2a.eq.2 .or. nv2a.eq.3 .or. nv2b.eq.2 .or. nv2b.eq.3) then
if(k1.lt.0 .or. k2.lt.0 .or. k1*k2.ne.0) go to 10
if(k2.gt.0) k2=k2+450
k=max(k1,k2)
if(k.gt.0) then
call k2grid(k,grid6)
c3=grid6(:4)
endif
endif
call packgrid(c3,ng,text3)
if(nv2a.lt.4 .and. nv2b.lt.4 .and. (.not.text1) .and. (.not.text2) .and. &
(.not.text3)) go to 20
nc1=0
if(nv2b.eq.4) then
if(c1(1:3).eq.'CQ ') nc1=262178563 + k2
if(c1(1:4).eq.'QRZ ') nc1=264002072 + k2
if(c1(1:3).eq.'DE ') nc1=265825581 + k2
else if(nv2b.eq.5) then
if(c1(1:3).eq.'CQ ') nc1=267649090 + k2
if(c1(1:4).eq.'QRZ ') nc1=267698375 + k2
if(c1(1:3).eq.'DE ') nc1=267747660 + k2
endif
if(nc1.ne.0) go to 20
! The message will be treated as plain text.
10 itype=6
call packtext(msg,nc1,nc2,ng)
ng=ng+32768
! Encode data into 6-bit words
20 continue
if(itype.ne.6) itype=max(nv2a,nv2b)
dat(1)=iand(ishft(nc1,-22),63) !6 bits
dat(2)=iand(ishft(nc1,-16),63) !6 bits
dat(3)=iand(ishft(nc1,-10),63) !6 bits
dat(4)=iand(ishft(nc1, -4),63) !6 bits
dat(5)=4*iand(nc1,15)+iand(ishft(nc2,-26),3) !4+2 bits
dat(6)=iand(ishft(nc2,-20),63) !6 bits
dat(7)=iand(ishft(nc2,-14),63) !6 bits
dat(8)=iand(ishft(nc2, -8),63) !6 bits
dat(9)=iand(ishft(nc2, -2),63) !6 bits
dat(10)=16*iand(nc2,3)+iand(ishft(ng,-12),15) !2+4 bits
dat(11)=iand(ishft(ng,-6),63)
dat(12)=iand(ng,63)
return
end subroutine packmsg
subroutine unpackmsg(dat,msg)
parameter (NBASE=37*36*10*27*27*27)
parameter (NGBASE=180*180)
integer dat(:)
character c1*12,c2*12,grid*4,msg*22,grid6*6,psfx*4,junk2*4
logical cqnnn
cqnnn=.false.
nc1=ishft(dat(1),22) + ishft(dat(2),16) + ishft(dat(3),10)+ &
ishft(dat(4),4) + iand(ishft(dat(5),-2),15)
nc2=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)
ng=ishft(iand(dat(10),15),12) + ishft(dat(11),6) + dat(12)
if(ng.ge.32768) then
call unpacktext(nc1,nc2,ng,msg)
go to 100
endif
call unpackcall(nc1,c1,iv2,psfx)
if(iv2.eq.0) then
! This is an "original JT65" message
if(nc1.eq.NBASE+1) c1='CQ '
if(nc1.eq.NBASE+2) c1='QRZ '
nfreq=nc1-NBASE-3
if(nfreq.ge.0 .and. nfreq.le.999) then
write(c1,1002) nfreq
1002 format('CQ ',i3.3)
cqnnn=.true.
endif
endif
call unpackcall(nc2,c2,junk1,junk2)
call unpackgrid(ng,grid)
if(iv2.gt.0) then
! This is a JT65v2 message
do i=1,4
if(ichar(psfx(i:i)).eq.0) psfx(i:i)=' '
enddo
n1=len_trim(psfx)
n2=len_trim(c2)
if(iv2.eq.1) msg='CQ '//psfx(:n1)//'/'//c2(:n2)//' '//grid
if(iv2.eq.2) msg='QRZ '//psfx(:n1)//'/'//c2(:n2)//' '//grid
if(iv2.eq.3) msg='DE '//psfx(:n1)//'/'//c2(:n2)//' '//grid
if(iv2.eq.4) msg='CQ '//c2(:n2)//'/'//psfx(:n1)//' '//grid
if(iv2.eq.5) msg='QRZ '//c2(:n2)//'/'//psfx(:n1)//' '//grid
if(iv2.eq.6) msg='DE '//c2(:n2)//'/'//psfx(:n1)//' '//grid
if(iv2.eq.7) msg='DE '//c2(:n2)//' '//grid
if(iv2.eq.8) msg=' '
go to 100
else
endif
grid6=grid//'ma'
call grid2k(grid6,k)
if(k.ge.1 .and. k.le.450) call getpfx2(k,c1)
if(k.ge.451 .and. k.le.900) call getpfx2(k,c2)
i=index(c1,char(0))
if(i.ge.3) c1=c1(1:i-1)//' '
i=index(c2,char(0))
if(i.ge.3) c2=c2(1:i-1)//' '
msg=' '
j=0
if(cqnnn) then
msg=c1//' '
j=7 !### ??? ###
go to 10
endif
do i=1,12
j=j+1
msg(j:j)=c1(i:i)
if(c1(i:i).eq.' ') go to 10
enddo
j=j+1
msg(j:j)=' '
10 do i=1,12
if(j.le.21) j=j+1
msg(j:j)=c2(i:i)
if(c2(i:i).eq.' ') go to 20
enddo
if(j.le.21) j=j+1
msg(j:j)=' '
20 if(k.eq.0) then
do i=1,4
if(j.le.21) j=j+1
msg(j:j)=grid(i:i)
enddo
if(j.le.21) j=j+1
msg(j:j)=' '
endif
100 continue
if(msg(1:6).eq.'CQ9DX ') msg(3:3)=' '
if(msg(1:2).eq.'E9' .and. &
msg(3:3).ge.'A' .and. msg(3:3).le.'Z' .and. &
msg(4:4).ge.'A' .and. msg(4:4).le.'Z' .and. &
msg(5:5).eq.' ') msg='CQ '//msg(3:)
return
end subroutine unpackmsg
subroutine packtext(msg,nc1,nc2,nc3)
parameter (MASK28=2**28 - 1)
character*13 msg
character*42 c
data c/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ +-./?'/
nc1=0
nc2=0
nc3=0
do i=1,5 !First 5 characters in nc1
do j=1,42 !Get character code
if(msg(i:i).eq.c(j:j)) go to 10
enddo
j=37
10 j=j-1 !Codes should start at zero
nc1=42*nc1 + j
enddo
do i=6,10 !Characters 6-10 in nc2
do j=1,42 !Get character code
if(msg(i:i).eq.c(j:j)) go to 20
enddo
j=37
20 j=j-1 !Codes should start at zero
nc2=42*nc2 + j
enddo
do i=11,13 !Characters 11-13 in nc3
do j=1,42 !Get character code
if(msg(i:i).eq.c(j:j)) go to 30
enddo
j=37
30 j=j-1 !Codes should start at zero
nc3=42*nc3 + j
enddo
! We now have used 17 bits in nc3. Must move one each to nc1 and nc2.
nc1=nc1+nc1
if(iand(nc3,32768).ne.0) nc1=nc1+1
nc2=nc2+nc2
if(iand(nc3,65536).ne.0) nc2=nc2+1
nc3=iand(nc3,32767)
return
end subroutine packtext
subroutine unpacktext(nc1,nc2,nc3,msg)
character*22 msg
character*44 c
data c/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ +-./?'/
nc3=iand(nc3,32767) !Remove the "plain text" bit
if(iand(nc1,1).ne.0) nc3=nc3+32768
nc1=nc1/2
if(iand(nc2,1).ne.0) nc3=nc3+65536
nc2=nc2/2
do i=5,1,-1
j=mod(nc1,42)+1
msg(i:i)=c(j:j)
nc1=nc1/42
enddo
do i=10,6,-1
j=mod(nc2,42)+1
msg(i:i)=c(j:j)
nc2=nc2/42
enddo
do i=13,11,-1
j=mod(nc3,42)+1
msg(i:i)=c(j:j)
nc3=nc3/42
enddo
msg(14:22) = ' '
return
end subroutine unpacktext
subroutine getpfx1(callsign,k,nv2)
character*12 callsign0,callsign,lof,rof
character*8 c
character addpfx*8,tpfx*4,tsfx*3
logical ispfx,issfx,invalid
common/pfxcom/addpfx
include 'pfx.f90'
callsign0=callsign
nv2=1
iz=index(callsign,' ') - 1
if(iz.lt.0) iz=12
islash=index(callsign(1:iz),'/')
k=0
! if(k.eq.0) go to 10 !Tnx to DL9RDZ for reminder:this was for tests only!
c=' '
if(islash.gt.0 .and. islash.le.(iz-4)) then
! Add-on prefix
c=callsign(1:islash-1)
callsign=callsign(islash+1:iz)
do i=1,NZ
if(pfx(i)(1:4).eq.c) then
k=i
nv2=2
go to 10
endif
enddo
if(addpfx.eq.c) then
k=449
nv2=2
go to 10
endif
else if(islash.eq.(iz-1)) then
! Add-on suffix
c=callsign(islash+1:iz)
callsign=callsign(1:islash-1)
do i=1,NZ2
if(sfx(i).eq.c(1:1)) then
k=400+i
nv2=3
go to 10
endif
enddo
endif
10 if(islash.ne.0 .and.k.eq.0) then
! Original JT65 would force this compound callsign to be treated as
! plain text. In JT65v2, we will encode the prefix or suffix into nc1.
! The task here is to compute the proper value of k.
lof=callsign0(:islash-1)
rof=callsign0(islash+1:)
llof=len_trim(lof)
lrof=len_trim(rof)
ispfx=(llof.gt.0 .and. llof.le.4)
issfx=(lrof.gt.0 .and. lrof.le.3)
invalid=.not.(ispfx.or.issfx)
if(ispfx.and.issfx) then
if(llof.lt.3) issfx=.false.
if(lrof.lt.3) ispfx=.false.
if(ispfx.and.issfx) then
i=ichar(callsign0(islash-1:islash-1))
if(i.ge.ichar('0') .and. i.le.ichar('9')) then
issfx=.false.
else
ispfx=.false.
endif
endif
endif
if(invalid) then
k=-1
else
if(ispfx) then
tpfx=lof(1:4)
k=nchar(tpfx(1:1))
k=37*k + nchar(tpfx(2:2))
k=37*k + nchar(tpfx(3:3))
k=37*k + nchar(tpfx(4:4))
nv2=4
i=index(callsign0,'/')
callsign=callsign0(:i-1)
callsign=callsign0(i+1:)
endif
if(issfx) then
tsfx=rof(1:3)
k=nchar(tsfx(1:1))
k=37*k + nchar(tsfx(2:2))
k=37*k + nchar(tsfx(3:3))
nv2=5
i=index(callsign0,'/')
callsign=callsign0(:i-1)
endif
endif
endif
return
end subroutine getpfx1
subroutine getpfx2(k0,callsign)
character callsign*12
include 'pfx.f90'
character addpfx*8
common/pfxcom/addpfx
k=k0
if(k.gt.450) k=k-450
if(k.ge.1 .and. k.le.NZ) then
iz=index(pfx(k),' ') - 1
callsign=pfx(k)(1:iz)//'/'//callsign
else if(k.ge.401 .and. k.le.400+NZ2) then
iz=index(callsign,' ') - 1
callsign=callsign(1:iz)//'/'//sfx(k-400)
else if(k.eq.449) then
iz=index(addpfx,' ') - 1
if(iz.lt.1) iz=8
callsign=addpfx(1:iz)//'/'//callsign
endif
return
end subroutine getpfx2
subroutine grid2k(grid,k)
character*6 grid
call grid2deg(grid,xlong,xlat)
nlong=nint(xlong)
nlat=nint(xlat)
k=0
if(nlat.ge.85) k=5*(nlong+179)/2 + nlat-84
return
end subroutine grid2k
subroutine k2grid(k,grid)
character grid*6
nlong=2*mod((k-1)/5,90)-179
if(k.gt.450) nlong=nlong+180
nlat=mod(k-1,5)+ 85
dlat=nlat
dlong=nlong
call deg2grid(dlong,dlat,grid)
return
end subroutine k2grid
subroutine grid2n(grid,n)
character*4 grid
i1=ichar(grid(1:1))-ichar('A')
i2=ichar(grid(3:3))-ichar('0')
i=10*i1 + i2
n=-i - 31
return
end subroutine grid2n
subroutine n2grid(n,grid)
character*4 grid
if(n.gt.-31 .or. n.lt.-70) stop 'Error in n2grid'
i=-(n+31) !NB: 0 <= i <= 39
i1=i/10
i2=mod(i,10)
grid(1:1)=char(ichar('A')+i1)
grid(2:2)='A'
grid(3:3)=char(ichar('0')+i2)
grid(4:4)='0'
return
end subroutine n2grid
function nchar(c)
! Convert ascii number, letter, or space to 0-36 for callsign packing.
character c*1
n=0 !Silence compiler warning
if(c.ge.'0' .and. c.le.'9') then
n=ichar(c)-ichar('0')
else if(c.ge.'A' .and. c.le.'Z') then
n=ichar(c)-ichar('A') + 10
else if(c.ge.'a' .and. c.le.'z') then
n=ichar(c)-ichar('a') + 10
else if(c.ge.' ') then
n=36
else
Print*,'Invalid character in callsign ',c,' ',ichar(c)
stop
endif
nchar=n
return
end function nchar
subroutine pack50(n1,n2,dat)
integer*1 dat(:),i1
i1=iand(ishft(n1,-20),255) !8 bits
dat(1)=i1
i1=iand(ishft(n1,-12),255) !8 bits
dat(2)=i1
i1=iand(ishft(n1, -4),255) !8 bits
dat(3)=i1
i1=16*iand(n1,15)+iand(ishft(n2,-18),15) !4+4 bits
dat(4)=i1
i1=iand(ishft(n2,-10),255) !8 bits
dat(5)=i1
i1=iand(ishft(n2, -2),255) !8 bits
dat(6)=i1
i1=64*iand(n2,3) !2 bits
dat(7)=i1
dat(8)=0
dat(9)=0
dat(10)=0
dat(11)=0
return
end subroutine pack50
subroutine packpfx(call1,n1,ng,nadd)
character*12 call1,call0
character*3 pfx
logical text
i1=index(call1,'/')
if(call1(i1+2:i1+2).eq.' ') then
! Single-character add-on suffix (maybe also fourth suffix letter?)
call0=call1(:i1-1)
call packcall(call0,n1,text)
nadd=1
nc=ichar(call1(i1+1:i1+1))
if(nc.ge.48 .and. nc.le.57) then
n=nc-48
else if(nc.ge.65 .and. nc.le.90) then
n=nc-65+10
else
n=38
endif
nadd=1
ng=60000-32768+n
else if(call1(i1+3:i1+3).eq.' ') then
! Two-character numerical suffix, /10 to /99
call0=call1(:i1-1)
call packcall(call0,n1,text)
nadd=1
n=10*(ichar(call1(i1+1:i1+1))-48) + ichar(call1(i1+2:i1+2)) - 48
nadd=1
ng=60000 + 26 + n
else
! Prefix of 1 to 3 characters
pfx=call1(:i1-1)
if(pfx(3:3).eq.' ') pfx=' '//pfx(1:2)
if(pfx(3:3).eq.' ') pfx=' '//pfx(1:2)
call0=call1(i1+1:)
call packcall(call0,n1,text)
ng=0
do i=1,3
nc=ichar(pfx(i:i))
if(nc.ge.48 .and. nc.le.57) then
n=nc-48
else if(nc.ge.65 .and. nc.le.90) then
n=nc-65+10
else
n=36
endif
ng=37*ng + n
enddo
nadd=0
if(ng.ge.32768) then
ng=ng-32768
nadd=1
endif
endif
return
end subroutine packpfx
end module packjt

96
q65w/libm65/setup65.f90
View File

@ -1,96 +0,0 @@
subroutine setup65
! Defines arrays related to the JT65 pseudo-random synchronizing pattern.
! Executed at program start.
integer nprc(126)
common/prcom/pr(126),mdat(126),mref(126,2),mdat2(126),mref2(126,2)
! JT65
data nprc/ &
1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0, &
0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1, &
0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1, &
0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1, &
1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1, &
0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1, &
1,1,1,1,1,1/
data mr2/0/ !Silence compiler warning
! Put the appropriate pseudo-random sequence into pr
nsym=126
do i=1,nsym
pr(i)=2*nprc(i)-1
enddo
! Determine locations of data and reference symbols
k=0
mr1=0
do i=1,nsym
if(pr(i).lt.0.0) then
k=k+1
mdat(k)=i
else
mr2=i
if(mr1.eq.0) mr1=i
endif
enddo
nsig=k
! Determine the reference symbols for each data symbol.
do k=1,nsig
m=mdat(k)
mref(k,1)=mr1
do n=1,10 !Get ref symbol before data
if((m-n).gt.0) then
if (pr(m-n).gt.0.0) go to 10
endif
enddo
go to 12
10 mref(k,1)=m-n
12 mref(k,2)=mr2
do n=1,10 !Get ref symbol after data
if((m+n).le.nsym) then
if (pr(m+n).gt.0.0) go to 20
endif
enddo
go to 22
20 mref(k,2)=m+n
22 enddo
! Now do it all again, using opposite logic on pr(i)
k=0
mr1=0
do i=1,nsym
if(pr(i).gt.0.0) then
k=k+1
mdat2(k)=i
else
mr2=i
if(mr1.eq.0) mr1=i
endif
enddo
nsig=k
do k=1,nsig
m=mdat2(k)
mref2(k,1)=mr1
do n=1,10
if((m-n).gt.0) then
if (pr(m-n).lt.0.0) go to 110
endif
enddo
go to 112
110 mref2(k,1)=m-n
112 mref2(k,2)=mr2
do n=1,10
if((m+n).le.nsym) then
if (pr(m+n).lt.0.0) go to 120
endif
enddo
go to 122
120 mref2(k,2)=m+n
122 enddo
return
end subroutine setup65