w4kek/WSJT-X
1
0
Fork 0
mirror of https://github.com/saitohirga/WSJT-X.git synced 2026-06-07 08:24:53 -04:00
Code Issues Packages Projects Releases Wiki Activity

Merging code for v1.1 back into the main wsjtx branch.

Browse Source 
git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@3462 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Joe Taylor
2013-07-08 13:17:22 +00:00
parent ba61cfde6a
commit 3dc6abc0c1
75 changed files with 4807 additions and 3514 deletions
Download Patch File Download Diff File Expand all files Collapse all files
lib/Makefile.MinGW
+22 -16
View File
@@ -20,7 +20,7 @@ CFLAGS = -I. -fbounds-check -mno-stack-arg-probe
%.o: %.F90
${FC} ${FFLAGS} -c $<
all: libjt9.a jt9sim.exe jt9.exe jt9code.exe test9.exe
all: libjt9.a jt9sim.exe jt9.exe jt9code.exe jt65.exe
OBJS1 = pctile.o graycode.o sort.o ssort.o \
unpackmsg.o igray.o unpackcall.o unpackgrid.o \
@@ -30,10 +30,15 @@ OBJS1 = pctile.o graycode.o sort.o ssort.o \
symspec.o analytic.o db.o genjt9.o \
packbits.o unpackbits.o encode232.o interleave9.o \
entail.o fano232.o gran.o sync9.o decode9.o \
fil3.o redsync.o decoder.o grid2n.o n2grid.o timer.o \
fil3.o decoder.o grid2n.o n2grid.o timer.o \
softsym.o getlags.o afc9.o fchisq.o twkfreq.o downsam9.o \
peakdt9.o symspec2.o stdmsg.o morse.o azdist.o geodist.o \
fillcom.o chkss2.o zplot9.o
fillcom.o chkss2.o zplot9.o flat2.o \
jt65a.o symspec65.o flat65.o ccf65.o decode65a.o \
filbig.o fil6521.o afc65b.o decode65b.o setup65.o \
extract.o fchisq65.o demod64a.o chkhist.o interleave63.o ccf2.o \
move.o indexx.o graycode65.o twkfreq65.o smo121.o \
wrapkarn.o init_rs.o encode_rs.o decode_rs.o gen65.o
libjt9.a: $(OBJS1)
ar cr libjt9.a $(OBJS1)
@@ -54,15 +59,9 @@ OBJS4 = jt9code.o
jt9code.exe: $(OBJS4) libjt9.a
$(FC) -o jt9code.exe $(OBJS4) libjt9.a
OBJS5 = test9.o
test9.exe: $(OBJS5) libjt9.a
$(FC) -o test9.exe $(OBJS5) libjt9.a ../libfftw3f_win.a
OBJS6 = wsjt24d.o wsjt24.o sync24.o decode24.o ps24.o flat1.o \
xcor24.o slope.o peakup.o interleave24.o getmet24.o smo.o \
deep24.o encode4.o chkmsg.o avemsg4.o extract4.o
wsjt24d.exe: $(OBJS6) libjt9.a
$(FC) -o wsjt24d.exe $(OBJS6) libjt9.a ../libfftw3f_win.a
OBJS5 = jt65.o
jt65.exe: $(OBJS5) libjt9.a
$(FC) -o jt65.exe $(OBJS5) libjt9.a ../libfftw3f_win.a
INCPATH = -I'c:/QtSDK/Desktop/Qt/4.7.4/mingw/include/QtCore' \
-I'c:/QtSDK/Desktop/Qt/4.7.4/mingw/include' \
@@ -75,9 +74,6 @@ sync9.o: sync9.f90 jt9sync.f90
spec9.o: spec9.f90 jt9sync.f90
$(FC) $(FFLAGS) -c spec9.f90
peakdf9.o: peakdf9.f90 jt9sync.f90
$(FC) $(FFLAGS) -c peakdf9.f90
peakdt9.o: peakdt9.f90 jt9sync.f90
$(FC) $(FFLAGS) -c peakdt9.f90
@@ -105,7 +101,17 @@ sec_midn.o: sec_midn.f90
tstrig.o: tstrig.c
$(CC) -c -Wall -I..\..\..\hamlib-1.2.15.3\include tstrig.c
init_rs.o: init_rs.c
$(CC) -c -DBIGSYM=1 -o init_rs.o init_rs.c
encode_rs.o: encode_rs.c
$(CC) -c -DBIGSYM=1 -o encode_rs.o encode_rs.c
decode_rs.o: decode_rs.c
$(CC) -c -DBIGSYM=1 -o decode_rs.o decode_rs.c
.PHONY : clean
clean:
rm -f *.o libjt9.a wsjtx.exe jt9sim.exe jt9.exe jt9test.exe
rm -f *.o libjt9.a wsjtx.exe jt9sim.exe jt9.exe jt65.exe
lib/Makefile.linux
+18 -4
View File
@@ -23,14 +23,19 @@ OBJS1 = pctile.o graycode.o sort.o ssort.o \
grid2k.o unpacktext.o getpfx2.o packmsg.o deg2grid.o \
packtext.o getpfx1.o packcall.o k2grid.o packgrid.o \
nchar.o four2a.o grid2deg.o pfxdump.o f77_wisdom.o \
symspec.o analytic.o db.o genjt9.o ptt_unix.o \
symspec.o analytic.o db.o genjt9.o \
packbits.o unpackbits.o encode232.o interleave9.o \
entail.o fano232.o gran.o sync9.o decode9.o \
fil3.o redsync.o decoder.o grid2n.o n2grid.o timer.o \
fil3.o decoder.o grid2n.o n2grid.o timer.o \
softsym.o peakdt9.o getlags.o afc9.o fchisq.o \
twkfreq.o downsam9.o symspec2.o ipcomm.o sleep_msec.o \
stdmsg.o sec_midn.o cutil.o azdist.o geodist.o morse.o \
fillcom.o chkss2.o
fillcom.o chkss2.o zplot9.o flat2.o \
jt65a.o symspec65.o flat65.o ccf65.o decode65a.o \
filbig.o fil6521.o afc65b.o decode65b.o setup65.o \
extract.o fchisq65.o demod64a.o chkhist.o interleave63.o ccf2.o \
move.o indexx.o graycode65.o twkfreq65.o smo121.o \
wrapkarn.o init_rs.o encode_rs.o decode_rs.o gen65.o
libjt9.a: $(OBJS1)
ar cr libjt9.a $(OBJS1)
@@ -38,7 +43,7 @@ libjt9.a: $(OBJS1)
OBJS2 = jt9.o jt9a.o jt9b.o jt9c.o
jt9: $(OBJS2) libjt9.a
jt9: $(OBJS2) libjt9.a
g++ -o jt9 $(OBJS2) libjt9.a -lfftw3f -lgfortran -lQtCore
cp jt9 ../../wsjtx_install
@@ -84,6 +89,15 @@ ipcomm.o: ipcomm.cpp
sec_midn.o: sec_midn.f90
$(FC) -c -fno-second-underscore sec_midn.f90
init_rs.o: init_rs.c
$(CC) -c -DBIGSYM=1 -o init_rs.o init_rs.c
encode_rs.o: encode_rs.c
$(CC) -c -DBIGSYM=1 -o encode_rs.o encode_rs.c
decode_rs.o: decode_rs.c
$(CC) -c -DBIGSYM=1 -o decode_rs.o decode_rs.c
.PHONY : clean
clean:
lib/afc65b.f90
+59
View File
@@ -0,0 +1,59 @@
subroutine afc65b(cx,npts,fsample,nflip,a,ccfbest,dtbest)
! Find delta f, f1, f2 ==> a(1:3)
complex cx(npts)
real a(5),deltaa(5)
a(1)=0.
a(2)=0.
a(3)=0.
a(4)=0.
deltaa(1)=2.0
deltaa(2)=2.0
deltaa(3)=2.0
deltaa(4)=0.05
nterms=3 !Maybe 2 is enough?
! Start the iteration
chisqr=0.
chisqr0=1.e6
do iter=1,3 !One iteration is enough?
do j=1,nterms
chisq1=fchisq65(cx,npts,fsample,nflip,a,ccfmax,dtmax)
fn=0.
delta=deltaa(j)
10 a(j)=a(j)+delta
chisq2=fchisq65(cx,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=fchisq65(cx,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=fchisq65(cx,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
return
end subroutine afc65b
lib/ccf2.f90
+45
View File
@@ -0,0 +1,45 @@
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
lib/ccf65.f90
+117
View File
@@ -0,0 +1,117 @@
subroutine ccf65(ss,nhsym,ssmax,sync1,dt1,flipk,syncshort,snr2,dt2)
parameter (NFFT=512,NH=NFFT/2)
real ss(322) !Input: half-symbol normalized powers
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)
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
! Look for JT65 sync pattern and shorthand square-wave pattern.
ccfbest=0.
ccfbest2=0.
do i=1,nhsym-1
s(i)=min(ssmax,ss(i)+ss(i+1))
! s(i)=ss(i)+ss(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)=s(j)
! if(abs(ccf(lag)).gt.ccfbest) then
if(ccf(lag).gt.ccfbest) then !No inverted sync for use at HF
! ccfbest=abs(ccf(lag))
ccfbest=ccf(lag)
lagpk=lag
flipk=1.0
! if(ccf(lag).lt.0.0) flipk=-1.0
endif
enddo
! do lag=-11,54 !Check for best shorthand
! ccf2=s2(lag+28)
! if(ccf2.gt.ccfbest2) then
! ccfbest2=ccf2
! lagpk2=lag
! endif
! enddo
! Find rms level on baseline of "ccfblue", for normalization.
sum=0.
do lag=-11,54
if(abs(lag-lagpk).gt.1) sum=sum + ccf(lag)
enddo
base=sum/50.0
sq=0.
do lag=-11,54
if(abs(lag-lagpk).gt.1) sq=sq + (ccf(lag)-base)**2
enddo
rms=sqrt(sq/49.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(i)
enddo
call pctile(tmp1,nhsym,40,base)
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))
dt2=0.
return
end subroutine ccf65
lib/char.h
+57
View File
@@ -0,0 +1,57 @@
/* Include file to configure the RS codec for character symbols
*
* Copyright 2002, Phil Karn, KA9Q
* May be used under the terms of the GNU General Public License (GPL)
*/
#define DTYPE unsigned char
/* Reed-Solomon codec control block */
struct rs {
int mm; /* Bits per symbol */
int nn; /* Symbols per block (= (1<<mm)-1) */
DTYPE *alpha_to; /* log lookup table */
DTYPE *index_of; /* Antilog lookup table */
DTYPE *genpoly; /* Generator polynomial */
int nroots; /* Number of generator roots = number of parity symbols */
int fcr; /* First consecutive root, index form */
int prim; /* Primitive element, index form */
int iprim; /* prim-th root of 1, index form */
int pad; /* Padding bytes in shortened block */
};
static inline int modnn(struct rs *rs,int x){
while (x >= rs->nn) {
x -= rs->nn;
x = (x >> rs->mm) + (x & rs->nn);
}
return x;
}
#define MODNN(x) modnn(rs,x)
#define MM (rs->mm)
#define NN (rs->nn)
#define ALPHA_TO (rs->alpha_to)
#define INDEX_OF (rs->index_of)
#define GENPOLY (rs->genpoly)
#define NROOTS (rs->nroots)
#define FCR (rs->fcr)
#define PRIM (rs->prim)
#define IPRIM (rs->iprim)
#define PAD (rs->pad)
#define A0 (NN)
#define ENCODE_RS encode_rs_char
#define DECODE_RS decode_rs_char
#define INIT_RS init_rs_char
#define FREE_RS free_rs_char
void ENCODE_RS(void *p,DTYPE *data,DTYPE *parity);
int DECODE_RS(void *p,DTYPE *data,int *eras_pos,int no_eras);
void *INIT_RS(int symsize,int gfpoly,int fcr,
int prim,int nroots,int pad);
void FREE_RS(void *p);
lib/chkhist.f90
+21
View File
@@ -0,0 +1,21 @@
subroutine chkhist(mrsym,nmax,ipk)
integer mrsym(63)
integer hist(0:63)
hist=0
do j=1,63
i=mrsym(j)
hist(i)=hist(i)+1
enddo
nmax=0
do i=0,63
if(hist(i).gt.nmax) then
nmax=hist(i)
ipk=i+1
endif
enddo
return
end subroutine chkhist
lib/constants.f90
+5
View File
@@ -0,0 +1,5 @@
parameter (NTMAX=120)
parameter (NMAX=NTMAX*12000) !Total sample intervals per 30 minutes
parameter (NDMAX=NTMAX*1500) !Sample intervals at 1500 Hz rate
parameter (NSMAX=6827) !Max length of saved spectra
parameter (MAXFFT3=16384)
lib/decode24.f90
-176
View File
@@ -1,176 +0,0 @@
subroutine decode24(dat,npts,dtx,dfx,flip,mode,mode4,decoded,ncount, &
deepbest,qbest,submode)
! Decodes JT65 data, assuming that DT and DF have already been determined.
parameter (MAXAVE=120)
real dat(npts) !Raw data
character decoded*22,deepmsg*22,deepbest*22
character*12 mycall,hiscall
character*6 hisgrid
character*72 c72
character submode*1
real*8 dt,df,phi,f0,dphi,twopi,phi1,dphi1
complex*16 cz,cz1,c0,c1
integer*1 symbol(207)
real*4 rsymbol(207,7)
real*4 sym(207)
integer nsum(7)
integer*1 data1(13) !Decoded data (8-bit bytes)
integer data4a(9) !Decoded data (8-bit bytes)
integer data4(12) !Decoded data (6-bit bytes)
integer amp,delta
integer mettab(0:255,0:1) !Metric table
integer nch(7)
integer npr2(207)
common/ave/ppsave(207,7,MAXAVE),nflag(MAXAVE),nsave,iseg(MAXAVE)
data mode0/-999/
data nsum/7*0/,rsymbol/1449*0.0/
data npr2/ &
0,0,0,0,1,1,0,0,0,1,1,0,1,1,0,0,1,0,1,0,0,0,0,0,0,0,1,1,0,0, &
0,0,0,0,0,0,0,0,0,0,1,0,1,1,0,1,1,0,1,0,1,1,1,1,1,0,1,0,0,0, &
1,0,0,1,0,0,1,1,1,1,1,0,0,0,1,0,1,0,0,0,1,1,1,1,0,1,1,0,0,1, &
0,0,0,1,1,0,1,0,1,0,1,0,1,0,1,1,1,1,1,0,1,0,1,0,1,1,0,1,0,1, &
0,1,1,1,0,0,1,0,1,1,0,1,1,1,1,0,0,0,0,1,1,0,1,1,0,0,0,1,1,1, &
0,1,1,1,0,1,1,1,0,0,1,0,0,0,1,1,0,1,1,0,0,1,0,0,0,1,1,1,1,1, &
1,0,0,1,1,0,0,0,0,1,1,0,0,0,1,0,1,1,0,1,1,1,1,0,1,0,1/
data nch/1,2,4,9,18,36,72/
save mettab,mode0,nsum,rsymbol
if(mode.ne.mode0) call getmet24(mode,mettab)
mode0=mode
twopi=8*atan(1.d0)
dt=2.d0/11025 !Sample interval (2x downsampled data)
df=11025.d0/2520.d0
nsym=206
amp=15
istart=nint(dtx/dt) !Start index for synced FFTs
if(istart.lt.0) istart=0
nchips=0
ich=0
qbest=0.
deepmsg=' '
ichbest=-1
! Should amp be adjusted according to signal strength?
! Compute soft symbols using differential BPSK demodulation
c0=0. !### C0=amp ???
k=istart
phi=0.d0
phi1=0.d0
40 ich=ich+1
nchips=nch(ich)
nspchip=1260/nchips
k=istart
phi=0.d0
phi1=0.d0
fac2=1.e-8 * sqrt(float(mode4))
do j=1,nsym+1
if(flip.gt.0.0) then
f0=1270.46 + dfx + (npr2(j)-1.5)*mode4*df
f1=1270.46 + dfx + (2+npr2(j)-1.5)*mode4*df
else
f0=1270.46 + dfx + (1-npr2(j)-1.5)*mode4*df
f1=1270.46 + dfx + (3-npr2(j)-1.5)*mode4*df
endif
dphi=twopi*dt*f0
dphi1=twopi*dt*f1
sq0=0.
sq1=0.
do nc=1,nchips
phi=0.d0
phi1=0.d0
c0=0.
c1=0.
do i=1,nspchip
k=k+1
phi=phi+dphi
phi1=phi1+dphi1
cz=dcmplx(cos(phi),-sin(phi))
cz1=dcmplx(cos(phi1),-sin(phi1))
if(k.le.npts) then
c0=c0 + dat(k)*cz
c1=c1 + dat(k)*cz1
endif
enddo
sq0=sq0 + real(c0)**2 + aimag(c0)**2
sq1=sq1 + real(c1)**2 + aimag(c1)**2
enddo
sq0=fac2*sq0
sq1=fac2*sq1
rsym=amp*(sq1-sq0)
r=rsym+128.
if(r.gt.255.0) r=255.0
if(r.lt.0.0) r=0.0
i4=nint(r)
if(i4.gt.127) i4=i4-256
if(j.ge.1) then
symbol(j)=i4
! rsymbol(j,ich)=rsymbol(j,ich) + rsym
rsymbol(j,ich)=rsym
sym(j)=rsym
endif
enddo
!### The following does simple message averaging:
! nsum(ich)=nsum(ich)+1
! do j=1,207
! sym(j)=rsymbol(j,ich)/nsum(ich)
! r=sym(j) + 128.
! if(r.gt.255.0) r=255.0
! if(r.lt.0.0) r=0.0
! i4=nint(r)
! if(i4.gt.127) i4=i4-256
! symbol(j)=i4
! enddo
!###
call extract4(sym,nadd,ncount,decoded) !Do the KV decode
qual=0. !Now try deep search
neme=1
mycall='VK7MO'
hiscall='W5LUA'
hisgrid='EM13'
call deep24(sym(2),neme,flip,mycall,hiscall,hisgrid,deepmsg,qual)
if(qual.gt.qbest) then
qbest=qual
deepbest=deepmsg
ichbest=ich
endif
if(ncount.ge.0) go to 100
if(mode.eq.7 .and. nchips.lt.mode4) go to 40
100 continue
!100 do i=1,9
! i4=data1(i)
! if(i4.lt.0) i4=i4+256
! data4a(i)=i4
! enddo
! write(c72,1100) (data4a(i),i=1,9)
!1100 format(9b8.8)
! read(c72,1102) data4
!1102 format(12b6)
! decoded=' '
! submode=' '
if(ncount.lt.0) then
decoded=deepbest
submode=char(ichar('A')+ichbest-1)
qual=qbest
endif
! if(decoded(1:6).eq.'000AAA') then
! decoded='***WRONG MODE?***'
! ncount=-1
! endif
! Save symbol spectra for possible decoding of average.
ppsave(1:207,1:7,nsave)=rsymbol(1:207,1:7)
return
end subroutine decode24
lib/decode65a.f90
+95
View File
@@ -0,0 +1,95 @@
subroutine decode65a(dd,npts,newdat,f0,nflip,mode65,sync2,a,dt, &
nbmkv,nhist,decoded)
! Apply AFC corrections to a candidate JT65 signal, then decode it.
parameter (NMAX=60*12000) !Samples per 60 s
real*4 dd(NMAX) !92 MB: raw data from Linrad timf2
complex cx(NMAX/8) !Data at 1378.125 samples/s
complex c5x(NMAX/32) !Data at 344.53125 Hz
complex c5a(512)
real s2(66,126)
real a(5)
logical first
character decoded*22
data first/.true./,jjjmin/1000/,jjjmax/-1000/
data 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,npts,f0,newdat,cx,n5,sq0)
call timer('filbig ',1)
! NB: cx has sample rate 12000*77125/672000 = 1378.125 Hz
! Find best DF, f1, f2, and DT. Start by downsampling to 344.53125 Hz
call timer('fil6521 ',0)
! 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)
n6=n6+nadd
call timer('fil6521 ',1)
fsample=1378.125/4.
a(5)=dt00
i0=nint((a(5)+0.5)*fsample) - 2 + nadd
if(i0.lt.1) then
! write(23,*) 'i0 too small in decode1a:',i0,f0,a(5),fsample,nadd
! flush(23)
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.
call timer('afc65b ',0)
! Best fit for DF, f1, and f2
call afc65b(c5x(i0),nz,fsample,nflip,a,ccfbest,dtbest)
call timer('afc65b ',1)
sync2=3.7e-4*ccfbest/sq0 !Constant is empirical
! 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 timer('twkfreq ',0)
call twkfreq65(cx,n5,a)
call timer('twkfreq ',1)
! Compute spectrum 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
call timer('sh_ffts ',0)
do k=1,nsym
do i=1,nfft
j=j+1
c5a(i)=cx(j)
enddo
call four2a(c5a,nfft,1,1,1)
do i=1,66
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
enddo
call timer('sh_ffts ',1)
call timer('dec65b ',0)
call decode65b(s2,nflip,mode65,nbmkv,nhist,decoded)
dt=dt00 + dtbest + 1.7
call timer('dec65b ',1)
return
end subroutine decode65a
lib/decode65b.f90
+36
View File
@@ -0,0 +1,36 @@
subroutine decode65b(s2,nflip,mode65,nbmkv,nhist,decoded)
real s2(66,126)
real s3(64,63)
logical first,ltext
character decoded*22
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(nflip.lt.0) k=mdat2(j)
do i=1,64
s3(i,j)=s2(i+2,k)
enddo
k=mdat2(j) !Points to data symbol
if(nflip.lt.0) k=mdat(j)
enddo
nadd=mode65
call extract(s3,nadd,ncount,nhist,decoded,ltext,nbmkv) !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(nflip.lt.0 .and. ltext) ncount=-1
if(ncount.lt.0) then
nbmkv=0
decoded=' '
endif
return
end subroutine decode65b
lib/decode9.f90
+2 -3
View File
@@ -16,7 +16,7 @@ subroutine decode9(i1SoftSymbols,limit,nlim,msg)
logical first
integer*4 mettab(0:255,0:1)
data first/.true./
data xx0/ &
data xx0/ & !Metric table
1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
@@ -53,9 +53,8 @@ subroutine decode9(i1SoftSymbols,limit,nlim,msg)
if(first) then
! Get the metric table
! bias=0.37 !To be optimized, in decoder program
bias=0.5
scale=10 ! ... ditto ...
scale=10
do i=0,255
mettab(i,0)=nint(scale*(xx0(i)-bias))
if(i.ge.1) mettab(256-i,1)=mettab(i,0)
lib/decode_rs.c
+1 -1
View File
@@ -9,7 +9,7 @@
#include <string.h>
#define NULL ((void *)0)
//#define NULL ((void *)0)
#define min(a,b) ((a) < (b) ? (a) : (b))
#ifdef FIXED
lib/decoder.f90
+48 -68
View File
@@ -1,23 +1,21 @@
subroutine decoder(ss,c0,nstandalone)
subroutine decoder(ss,id2)
! Decoder for JT9.
parameter (NTMAX=120)
parameter (NMAX=NTMAX*12000) !Total sample intervals per 30 minutes
parameter (NDMAX=NTMAX*1500) !Sample intervals at 1500 Hz rate
parameter (NSMAX=1365) !Max length of saved spectra
include 'constants.f90'
real ss(184,NSMAX)
character*22 msg
character*80 fmt
character*20 datetime
real*4 ccfred(NSMAX)
real*4 red2(NSMAX)
logical ccfok(NSMAX)
logical done(NSMAX)
logical done65
integer*2 id2(NTMAX*12000)
real*4 dd(NTMAX*12000)
integer*1 i1SoftSymbols(207)
complex c0(NDMAX)
common/npar/nutc,ndiskdat,ntrperiod,nfqso,newdat,npts8,nfa,nfb,ntol, &
kin,nzhsym,nsave,nagain,ndepth,nrxlog,nfsample,datetime
kin,nzhsym,nsave,nagain,ndepth,ntxmode,nmode,datetime
common/tracer/limtrace,lu
save
@@ -27,48 +25,40 @@ subroutine decoder(ss,c0,nstandalone)
ndecodes0=0
ndecodes1=0
call timer('decoder ',0)
open(13,file='decoded.txt',status='unknown')
ntrMinutes=ntrperiod/60
newdat=1
open(22,file='kvasd.dat',access='direct',recl=1024,status='unknown')
npts65=52*12000
ntol65=20
done65=.false.
if(nmode.ge.65 .and. ntxmode.eq.65) then
if(newdat.ne.0) dd(1:npts65)=id2(1:npts65)
call jt65a(dd,npts65,newdat,nutc,nfa,nfqso,ntol65,nagain,ndecoded)
done65=.true.
endif
if(nmode.eq.65) go to 800
nsynced=0
ndecoded=0
nsps=0
if(ntrMinutes.eq.1) then
nsps=6912
df3=1500.0/2048.0
fmt='(i4.4,i4,i5,f6.1,f8.0,i4,3x,a22)'
else if(ntrMinutes.eq.2) then
nsps=15360
df3=1500.0/2048.0
fmt='(i4.4,i4,i5,f6.1,f8.0,i4,3x,a22)'
else if(ntrMinutes.eq.5) then
nsps=40960
df3=1500.0/6144.0
fmt='(i4.4,i4,i5,f6.1,f8.1,i4,3x,a22)'
else if(ntrMinutes.eq.10) then
nsps=82944
df3=1500.0/12288.0
fmt='(i4.4,i4,i5,f6.1,f8.2,i4,3x,a22)'
else if(ntrMinutes.eq.30) then
nsps=252000
df3=1500.0/32768.0
fmt='(i4.4,i4,i5,f6.1,f8.2,i4,3x,a22)'
endif
if(nsps.eq.0) stop 'Error: bad TRperiod' !Better: return an error code###
nsps=6912 !Params for JT9-1
df3=1500.0/2048.0
tstep=0.5*nsps/12000.0 !Half-symbol step (seconds)
done=.false.
ia=max(1,nint((nfa-1000)/df3))
ib=min(NSMAX,nint((nfb-1000)/df3))
nf0=0
ia=max(1,nint((nfa-nf0)/df3))
ib=min(NSMAX,nint((nfb-nf0)/df3))
lag1=-(2.5/tstep + 0.9999)
lag2=5.0/tstep + 0.9999
call timer('sync9 ',0)
call sync9(ss,nzhsym,lag1,lag2,ia,ib,ccfred,red2,ipk)
call timer('sync9 ',1)
if(newdat.ne.0) then
call timer('sync9 ',0)
call sync9(ss,nzhsym,lag1,lag2,ia,ib,ccfred,red2,ipk)
call timer('sync9 ',1)
endif
nsps8=nsps/8
df8=1500.0/nsps8
@@ -94,24 +84,22 @@ subroutine decoder(ss,c0,nstandalone)
ccfok(ia:ib)=.true.
nfa1=nfqso-ntol
nfb1=nfqso+ntol
ia=max(1,nint((nfa1-1000)/df3))
ib=min(NSMAX,nint((nfb1-1000)/df3))
ia=max(1,nint((nfa1-nf0)/df3))
ib=min(NSMAX,nint((nfb1-nf0)/df3))
ia1=ia
ib1=ib
else
nfa1=nfa
nfb1=nfb
ia=max(1,nint((nfa1-1000)/df3))
ib=min(NSMAX,nint((nfb1-1000)/df3))
ia=max(1,nint((nfa1-nf0)/df3))
ib=min(NSMAX,nint((nfb1-nf0)/df3))
do i=ia,ib
ccfok(i)=ccfred(i).gt.ccflim .and. red2(i).gt.red2lim
enddo
ccfok(ia1:ib1)=.false.
endif
nRxLog=0
fgood=0.
do i=ia,ib
f=(i-1)*df3
if(done(i) .or. (.not.ccfok(i)) .or. (ccfred(i).lt.ccflim-1.0)) cycle
@@ -122,9 +110,9 @@ subroutine decoder(ss,c0,nstandalone)
if(nqd.eq.1) nfreqs1=nfreqs1+1
call timer('softsym ',0)
fpk=1000.0 + df3*(i-1)
call softsym(c0,npts8,nsps8,newdat,fpk,syncpk,snrdb,xdt,freq, &
drift,schk,i1SoftSymbols)
fpk=nf0 + df3*(i-1)
call softsym(id2,npts8,nsps8,newdat,fpk,syncpk,snrdb,xdt, &
freq,drift,schk,i1SoftSymbols)
call timer('softsym ',1)
if(schk.ge.schklim) then
@@ -140,16 +128,14 @@ subroutine decoder(ss,c0,nstandalone)
nsnr=nint(snrdb)
ndrift=nint(drift/df3)
! write(38,3002) nutc,nqd,nsnr,i,freq,ndrift,ccfred(i), &
! red2(i),schk,nlim,msg
!3002 format(i4.4,i2,i4,i5,f7.1,i4,f5.1,f6.1,f5.1,i8,1x,a22)
if(msg.ne.' ') then
if(nqd.eq.0) ndecodes0=ndecodes0+1
if(nqd.eq.1) ndecodes1=ndecodes1+1
write(*,fmt) nutc,nsync,nsnr,xdt,freq,ndrift,msg
write(13,fmt) nutc,nsync,nsnr,xdt,freq,ndrift,msg
write(*,1000) nutc,nsnr,xdt,nint(freq),msg
1000 format(i4.4,i4,f5.1,i5,1x,'@',1x,a22)
write(13,1002) nutc,nsync,nsnr,xdt,freq,ndrift,msg
1002 format(i4.4,i4,i5,f6.1,f8.0,i4,3x,a22,' JT9')
iaa=max(1,i-1)
ibb=min(NSMAX,i+22)
@@ -160,8 +146,6 @@ subroutine decoder(ss,c0,nstandalone)
done(iaa:ibb)=.true.
call flush(6)
endif
else
! write(38,3002) nutc,nqd,-99,i,freq,ndrift,ccfred(i),red2(i),schk,0
endif
endif
enddo
@@ -169,21 +153,17 @@ subroutine decoder(ss,c0,nstandalone)
if(nagain.ne.0) exit
enddo
write(*,1010) nsynced,ndecoded
if(nmode.ge.65 .and. (.not.done65)) then
if(newdat.ne.0) dd(1:npts65)=id2(1:npts65)
call jt65a(dd,npts65,newdat,nutc,nfa,nfqso,ntol65,nagain,ndecoded)
endif
!### JT65 is not yet producing info for nsynced, ndecoded.
800 write(*,1010) nsynced,ndecoded
1010 format('<DecodeFinished>',2i4)
call flush(6)
close(13)
! call flush(14)
call timer('decoder ',1)
if(nstandalone.eq.0) call timer('decoder ',101)
call system_clock(iclock,iclock_rate,iclock_max)
! write(39,3001) nutc,nfreqs1,nfreqs0,ndecodes1,ndecodes0, &
! float(iclock-iclock0)/iclock_rate
!3001 format(5i8,f10.3)
! call flush(38)
! call flush(39)
close(22)
return
end subroutine decoder
lib/demod64a.f90
+73
View File
@@ -0,0 +1,73 @@
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)
if(nadd.eq.-999) return
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.
do i=1,64
x=min(afac*s3(i,j)/ave,50.d0)
fs(i)=exp(x)
fsum=fsum+fs(i)
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
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
return
end subroutine demod64a
lib/downsam9.f90
+32 -30
View File
@@ -1,58 +1,60 @@
subroutine downsam9(c0,npts8,nsps8,newdat,nspsd,fpk,c2,nz2)
subroutine downsam9(id2,npts8,nsps8,newdat,nspsd,fpk,c2,nz2)
!Downsample to nspsd samples per symbol, info centered at fpk
!Downsample from id2() into C2() so as to yield nspsd samples per symbol,
!mixing from fpk down to zero frequency.
parameter (NMAX=128*31500)
complex c0(0:npts8-1)
complex c1(0:NMAX-1)
include 'constants.f90'
parameter (NMAX1=1024*1920)
integer*2 id2(0:8*npts8-1)
real*4 x1(0:NMAX1-1)
complex c1(0:NMAX1/2)
complex c2(0:4096-1)
real s(1000)
real s(5000)
equivalence (c1,x1)
save
nfft1=128*nsps8 !Forward FFT length
nh1=nfft1/2
df1=1500.0/nfft1
nfft1=1024*nsps8 !Forward FFT length
df1=12000.0/nfft1
npts=8*npts8
if(newdat.eq.1) then
fac=1.e-4
do i=0,npts8-1,2
c1(i)=fac*conjg(c0(i))
c1(i+1)=-fac*conjg(c0(i+1))
fac=6.963e-6 !Why this weird constant?
do i=0,npts-1
x1(i)=fac*id2(i)
enddo
c1(npts8:)=0. !Zero the rest of c1
call four2a(c1,nfft1,1,-1,1) !Forward FFT
x1(npts:nfft1-1)=0. !Zero the rest of x1
call four2a(c1,nfft1,1,-1,0) !Forward FFT, r2c
nadd=1.0/df1
j=250/df1
s=0.
do i=1,1000
do i=1,5000
j=(i-1)/df1
do n=1,nadd
j=j+1
s(i)=s(i)+real(c1(j))**2 + aimag(c1(j))**2
enddo
! write(37,3001) i+1000,s(i),db(s(i)),nadd
!3001 format(i5,2f12.3,i8)
enddo
call pctile(s,1000,40,avenoise)
endif
ndown=nsps8/16 !Downsample factor
ndown=8*nsps8/nspsd !Downsample factor
nfft2=nfft1/ndown !Backward FFT length
nh2=nfft2/2
fshift=fpk-1500.0
i0=nh1 + fshift/df1
nf=nint(fpk)
i0=fpk/df1
nw=100
ia=max(1,nf-nw)
ib=min(5000,nf+nw)
call pctile(s(ia),ib-ia+1,40,avenoise)
fac=sqrt(1.0/avenoise)
do i=0,nfft2-1
j=i0+i
if(i.gt.nh2) j=j-nfft2
c2(i)=fac*c1(j)
enddo
call four2a(c2,nfft2,1,1,1) !Backward FFT
nspsd=nsps8/ndown
nz2=npts8/ndown
call four2a(c2,nfft2,1,1,1) !FFT back to time domain
nz2=8*npts8/ndown
return
end subroutine downsam9
lib/extract.F90
+102
View File
@@ -0,0 +1,102 @@
subroutine extract(s3,nadd,ncount,nhist,decoded,ltext,nbmkv)
real s3(64,63)
character decoded*22
integer era(51),dat4(12),indx(64)
integer mrsym(63),mr2sym(63),mrprob(63),mr2prob(63)
logical nokv,ltext
data nokv/.false./,nsec1/0/
save
nbmkv=0
nfail=0
1 continue
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
do j=1,63
s3(ipk,j)=base
enddo
if(nfail.gt.30) then
decoded=' '
ncount=-1
go to 900
endif
go to 1
endif
call graycode65(mrsym,63,-1)
call interleave63(mrsym,-1)
call interleave63(mrprob,-1)
! Decode using Berlekamp-Massey algorithm
nemax=30 !Max BM erasures
call indexx(63,mrprob,indx)
do i=1,nemax
j=indx(i)
if(mrprob(j).gt.120) then
ne2=i-1
go to 2
endif
era(i)=j-1
enddo
ne2=nemax
2 decoded=' '
do nerase=0,ne2,2
call rs_decode(mrsym,era,nerase,dat4,ncount)
if(ncount.ge.0) then
call unpackmsg(dat4,decoded)
if(iand(dat4(10),8).ne.0) ltext=.true.
nbmkv=1
go to 900
endif
enddo
! Berlekamp-Massey algorithm failed, try Koetter-Vardy
if(nokv) go to 900
maxe=8 !Max KV errors in 12 most reliable symbols
xlambda=10.0
call graycode65(mr2sym,63,-1)
call interleave63(mr2sym,-1)
call interleave63(mr2prob,-1)
nsec1=nsec1+1
write(22,rec=1) nsec1,xlambda,maxe,200,mrsym,mrprob,mr2sym,mr2prob
call flush(22)
call timer('kvasd ',0)
#ifdef UNIX
iret=system('./kvasd -q > dev_null')
#else
iret=system('kvasd -q > dev_null')
#endif
call timer('kvasd ',1)
if(iret.ne.0) then
if(.not.nokv) write(*,1000)
1000 format('Error in KV decoder, or no KV decoder present.')
nokv=.true.
go to 900
endif
read(22,rec=2,err=900) nsec2,ncount,dat4
j=nsec2 !Silence compiler warning
decoded=' '
ltext=.false.
if(ncount.ge.0) then
call unpackmsg(dat4,decoded) !Unpack the user message
if(iand(dat4(10),8).ne.0) ltext=.true.
nbmkv=2
endif
900 continue
return
end subroutine extract
lib/fchisq65.f90
+68
View File
@@ -0,0 +1,68 @@
real function fchisq65(cx,npts,fsample,nflip,a,ccfmax,dtmax)
parameter (NMAX=60*12000) !Samples per 60 s
complex cx(npts)
real a(5)
complex w,wstep,z
real ss(3000)
complex csx(0:NMAX/8)
data twopi/6.283185307/a1,a2,a3/99.,99.,99./
save
call timer('fchisq65',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 signal
csx(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
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)
enddo
endif
! Compute 1/2-symbol powers at 1/16-symbol steps.
fac=1.e-4
do i=1,nout
j=i*nsps/ndiv
k=j-nsph
ss(i)=0.
if(k.ge.1) then
z=csx(j)-csx(k)
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
fchisq65=-ccfmax
call timer('fchisq65',1)
return
end function fchisq65
lib/fil6521.f90
+44
View File
@@ -0,0 +1,44 @@
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/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
lib/filbig.f90
+127
View File
@@ -0,0 +1,127 @@
subroutine filbig(dd,npts,f0,newdat,c4a,n4,sq0)
! Filter and downsample the real data in array dd(npts), sampled at 12000 Hz.
! Output is complex, sampled at 1378.125 Hz.
parameter (NSZ=3413)
parameter (NFFT1=672000,NFFT2=77175)
parameter (NZ2=1000)
real*4 dd(npts) !Input data
complex ca(NFFT1) !FFT of input
complex c4a(NFFT2) !Output data
real*4 s(NZ2)
real*8 df
real halfpulse(8) !Impulse response of filter (one sided)
complex cfilt(NFFT2) !Filter (complex; imag = 0)
real rfilt(NFFT2) !Filter (real)
integer*8 plan1,plan2,plan3
logical first
include 'fftw3.f90'
equivalence (rfilt,cfilt)
data first/.true./,npatience/0/
data halfpulse/114.97547150,36.57879257,-20.93789101, &
5.89886379,1.59355187,-2.49138308,0.60910773,-0.04248129/
common/refspec/dfref,ref(NSZ)
save
if(npts.lt.0) go to 900 !Clean up at end of program
if(first) then
nflags=FFTW_ESTIMATE
if(npatience.eq.1) nflags=FFTW_ESTIMATE_PATIENT
if(npatience.eq.2) nflags=FFTW_MEASURE
if(npatience.eq.3) nflags=FFTW_PATIENT
if(npatience.eq.4) nflags=FFTW_EXHAUSTIVE
! Plan the FFTs just once
call timer('FFTplans ',0)
call sfftw_plan_dft_1d(plan1,nfft1,ca,ca,FFTW_BACKWARD,nflags)
call sfftw_plan_dft_1d(plan2,nfft2,c4a,c4a,FFTW_FORWARD,nflags)
call sfftw_plan_dft_1d(plan3,nfft2,cfilt,cfilt,FFTW_BACKWARD,nflags)
call timer('FFTplans ',1)
! 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 timer('FFTfilt ',0)
call sfftw_execute(plan3)
call timer('FFTfilt ',1)
base=cfilt(nfft2/2+1)
do i=1,nfft2
rfilt(i)=real(cfilt(i))-base
enddo
df=12000.d0/nfft1
first=.false.
endif
! When new data comes along, we need to compute a new "big FFT"
! If we just have a new f0, continue with the existing data in ca.
if(newdat.ne.0) then
nz=min(npts,nfft1)
ca(1:nz)=dd(1:nz)
ca(nz+1:)=0. !### Should change this to r2c FFT ###
call timer('FFTbig ',0)
call sfftw_execute(plan1)
call timer('FFTbig ',1)
do i=1,NFFT1/2 !Flatten the spectrum
j=nint(i*df/dfref)
if(j.lt.1) j=1
if(j.gt.NSZ) j=NSZ
fac=sqrt(min(30.0,1.0/ref(j)))
ca(i)=fac * ca(i)
enddo
endif
! NB: f0 is the frequency at which we want our filter centered.
! i0 is the bin number in ca closest to f0.
i0=nint(f0/df) + 1
nh=nfft2/2
do i=1,nh !Copy data into c4a and apply
j=i0+i-1 !the filter function
if(j.ge.1 .and. j.le.nfft1) then
c4a(i)=rfilt(i)*ca(j)
else
c4a(i)=0.
endif
enddo
do i=nh+1,nfft2
j=i0+i-1-nfft2
if(j.lt.1) j=j+nfft1 !nfft1 was nfft2
c4a(i)=rfilt(i)*ca(j)
enddo
nadd=nfft2/NZ2
i=0
do j=1,NZ2
s(j)=0.
do n=1,nadd
i=i+1
s(j)=s(j) + real(c4a(i))**2 + aimag(c4a(i))**2
enddo
enddo
call pctile(s,NZ2,30,sq0)
! Do the short reverse transform, to go back to time domain.
call timer('FFTsmall',0)
call sfftw_execute(plan2)
call timer('FFTsmall',1)
n4=min(npts/8,nfft2)
return
900 call sfftw_destroy_plan(plan1)
call sfftw_destroy_plan(plan2)
call sfftw_destroy_plan(plan3)
return
end subroutine filbig
lib/fillcom.f90
+6 -6
View File
@@ -1,24 +1,24 @@
subroutine fillcom(nutc0,ndepth0)
character*20 datetime
common/npar/nutc,ndiskdat,ntrperiod,nfqso,newdat,npts8,nfa,nfb,ntol, &
kin,nzhsym,nsave,nagain,ndepth,nrxlog,nfsample,datetime
kin,nzhsym,nsave,nagain,ndepth,ntxmode,nmode,datetime
save
nutc=nutc0
ndiskdat=1
ntrperiod=60
nfqso=1500
nfqso=1197
newdat=1
npts8=74736
nfa=1000
nfb=2000
nfa=2700
nfb=4007
ntol=3
kin=1024
nzhsym=173
nsave=0
ndepth=ndepth0
nrxlog=1
nfsample=12000
ntxmode=9
nmode=9+65
datetime="2013-Apr-16 15:13"
return
lib/flat2.f90
+18
View File
@@ -0,0 +1,18 @@
subroutine flat2(s,nz,ref)
parameter (NSMAX=6827)
real s(NSMAX)
real ref(NSMAX)
nsmo=10
ia=nsmo+1
ib=nz-nsmo-1
do i=ia,ib
call pctile(s(i-nsmo),2*nsmo+1,5,ref(i))
enddo
ref(:ia-1)=ref(ia)
ref(ib+1:)=ref(ib)
return
end subroutine flat2
lib/flat65.f90
+25
View File
@@ -0,0 +1,25 @@
subroutine flat65(ss,nhsym,maxhsym,nsz,ref)
real stmp(nsz)
real ss(maxhsym,nsz)
real ref(nsz)
npct=28 !Somewhat arbitrary
do i=1,nsz
call pctile(ss(1,i),nhsym,npct,stmp(i))
enddo
nsmo=33
ia=nsmo/2 + 1
ib=nsz - nsmo/2 - 1
do i=ia,ib
call pctile(stmp(i-nsmo/2),nsmo,npct,ref(i))
enddo
ref(:ia-1)=ref(ia)
ref(ib+1:)=ref(ib)
ref=4.0*ref
return
end subroutine flat65
lib/gen65.f90
+65
View File
@@ -0,0 +1,65 @@
subroutine gen65(msg0,ichk,msgsent,itone,itext)
! Encodes a JT65 message to yieild itone(1:126)
! Temporarily, does not implement EME shorthands
character*22 msg0
character*22 message !Message to be generated
character*22 msgsent !Message as it will be received
integer itone(126)
! character*3 cok !' ' or 'OOO'
integer dgen(13)
integer sent(63)
logical text
integer nprc(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/
save
message=msg0
do i=1,22
if(ichar(message(i:i)).eq.0) then
message(i:)=' '
exit
endif
enddo
do i=1,22 !Strip leading blanks
if(message(1:1).ne.' ') exit
message=message(i+1:)
enddo
nspecial=0
! call chkmsg(message,cok,nspecial,flip)
if(nspecial.eq.0) then
call packmsg(message,dgen,text) !Pack message into 72 bits
itext=0
if(text) itext=1
call unpackmsg(dgen,msgsent) !Unpack to get message sent
if(ichk.ne.0) go to 999 !Return if checking only
call rs_encode(dgen,sent) !Apply Reed-Solomon code
call interleave63(sent,1) !Apply interleaving
call graycode65(sent,63,1) !Apply Gray code
nsym=126 !Symbols per transmission
else
nsym=32
endif
k=0
do j=1,nsym
if(nprc(j).eq.0) then
k=k+1
itone(j)=sent(k)+2
else
itone(j)=0
endif
enddo
999 return
end subroutine gen65
lib/genjt9.f90
+1 -1
View File
@@ -26,7 +26,7 @@ subroutine genjt9(msg0,ichk,msgsent,i4tone,itext)
endif
enddo
do i=1,22 !Omit leading blanks
do i=1,22 !Strip leading blanks
if(message(1:1).ne.' ') exit
message=message(i+1:)
enddo
lib/graycode65.f90
+9
View File
@@ -0,0 +1,9 @@
subroutine graycode65(dat,n,idir)
integer dat(n)
do i=1,n
dat(i)=igray(dat(i),idir)
enddo
return
end subroutine graycode65
lib/init_rs.c
+1 -1
View File
@@ -13,7 +13,7 @@
#include "char.h"
#endif
#define NULL ((void *)0)
//#define NULL ((void *)0)
void FREE_RS(void *p){
struct rs *rs = (struct rs *)p;
lib/interleave63.f90
+25
View File
@@ -0,0 +1,25 @@
subroutine interleave63(d1,idir)
! Interleave (idir=1) or de-interleave (idir=-1) the array d1.
integer d1(0:6,0:8)
integer d2(0:8,0:6)
if(idir.ge.0) then
do i=0,6
do j=0,8
d2(j,i)=d1(i,j)
enddo
enddo
call move(d2,d1,63)
else
call move(d1,d2,63)
do i=0,6
do j=0,8
d1(i,j)=d2(j,i)
enddo
enddo
endif
return
end subroutine interleave63
lib/jt65.f90
+68
View File
@@ -0,0 +1,68 @@
program jt65
! Test the JT65 decoder for WSJT-X
parameter (NZMAX=60*12000)
integer*4 ihdr(11)
integer*2 id2(NZMAX)
real*4 dd(NZMAX)
character*80 infile
integer*2 nfmt2,nchan2,nbitsam2,nbytesam2
character*4 ariff,awave,afmt,adata
common/hdr/ariff,lenfile,awave,afmt,lenfmt,nfmt2,nchan2, &
nsamrate,nbytesec,nbytesam2,nbitsam2,adata,ndata
common/tracer/limtrace,lu
equivalence (ariff,ihdr)
nargs=iargc()
if(nargs.lt.1) then
print*,'Usage: jt65 file1 [file2 ...]'
go to 999
endif
limtrace=0
lu=12
newdat=1
ntol=50
nfa=2700
! nfb=4000
nfqso=933
nagain=0
open(12,file='timer.out',status='unknown')
open(22,file='kvasd.dat',access='direct',recl=1024,status='unknown')
call timer('jt65 ',0)
do ifile=1,nargs
call getarg(ifile,infile)
open(10,file=infile,access='stream',status='old',err=998)
call timer('read ',0)
read(10) ihdr
nutc=ihdr(1) !Silence compiler warning
i1=index(infile,'.wav')
read(infile(i1-4:i1-1),*,err=10) nutc
go to 20
10 nutc=0
20 npts=52*12000
read(10) id2(1:npts)
call timer('read ',1)
dd(1:npts)=id2(1:npts)
dd(npts+1:)=0.
call timer('jt65a ',0)
call jt65a(dd,npts,newdat,nutc,ntol,nfa,nfqso,nagain,ndecoded)
call timer('jt65a ',1)
enddo
call timer('jt65 ',1)
call timer('jt65 ',101)
call four2a(a,-1,1,1,1) !Free the memory used for plans
call filbig(a,-1,1,0.0,0,0,0,0,0) ! (ditto)
go to 999
998 print*,'Cannot open file:'
print*,infile
999 end program jt65
lib/jt65a.f90
+91
View File
@@ -0,0 +1,91 @@
subroutine jt65a(dd,npts,newdat,nutc,nfa,nfqso,ntol,nagain,ndecoded)
! Process dd() data to find and decode JT65 signals.
parameter (NSZ=3413)
parameter (NZMAX=60*12000)
parameter (NFFT=8192)
real dd(NZMAX)
real*4 ss(322,NSZ)
real*4 savg(NSZ)
logical done(NSZ)
real a(5)
character decoded*22
save
if(newdat.ne.0) then
call timer('symsp65 ',0)
call symspec65(dd,npts,ss,nhsym,savg) !Get normalized symbol spectra
call timer('symsp65 ',1)
endif
df=12000.0/NFFT !df = 12000.0/16384 = 0.732 Hz
ftol=15.0 !Frequency tolerance (Hz)
mode65=1 !Decoding JT65A only, for now.
done=.false.
do nqd=1,0,-1
if(nqd.eq.1) then !Quick decode, at fQSO
fa=nfqso - ntol
fb=nfqso + ntol
else !Wideband decode at all freqs
fa=200
fb=nfa
endif
ia=max(51,nint(fa/df))
ib=min(NSZ-51,nint(fb/df))
freq0=-999.
thresh0=1.5
do i=ia,ib !Search over freq range
if(savg(i).lt.thresh0 .or. done(i)) cycle
freq=i*df
call timer('ccf65 ',0)
call ccf65(ss(1,i),nhsym,savg(i),sync1,dt,flipk,syncshort,snr2,dt2)
call timer('ccf65 ',1)
! ########################### Search for Shorthand Messages #################
! include 'shorthand1.f90'
! ########################### Search for Normal Messages ###########
thresh1=1.0
! Use lower thresh1 at fQSO
if(nqd.eq.1 .and. ntol.le.100) thresh1=0.
! Is sync1 above threshold?
if(sync1.lt.thresh1) cycle
! Keep only the best candidate within ftol.
if(freq-freq0.lt.ftol) cycle
nflip=nint(flipk)
call timer('decod65a',0)
call decode65a(dd,npts,newdat,freq,nflip,mode65,sync2,a,dt, &
nbmkv,nhist,decoded)
call timer('decod65a',1)
if(decoded.ne.' ') then
ndecoded=1
nfreq=nint(freq+a(1))
ndrift=nint(2.0*a(2))
s2db=10.0*log10(sync2) - 32 !### empirical (was 40) ###
nsnr=nint(s2db)
if(nsnr.lt.-30) nsnr=-30
if(nsnr.gt.-1) nsnr=-1
write(*,1010) nutc,nsnr,dt,nfreq,decoded
1010 format(i4.4,i4,f5.1,i5,1x,'#',1x,a22)
write(13,1012) nutc,nint(sync1),nsnr,dt,float(nfreq),ndrift, &
decoded,nbmkv
1012 format(i4.4,i4,i5,f6.1,f8.0,i4,3x,a22,' JT65',i4)
freq0=freq
i2=min(NSZ,i+15) !### ??? ###
done(i:i2)=.true.
endif
enddo
if(nagain.eq.1) exit
enddo
return
end subroutine jt65a
lib/jt9.f90
+4 -20
View File
@@ -3,18 +3,12 @@ program jt9
! Decoder for JT9. Can run stand-alone, reading data from *.wav files;
! or as the back end of wsjt-x, with data placed in a shared memory region.
parameter (NTMAX=120)
parameter (NMAX=NTMAX*12000) !Total sample intervals per 30 minutes
parameter (NDMAX=NTMAX*1500) !Sample intervals at 1500 Hz rate
parameter (NSMAX=1365) !Max length of saved spectra
include 'constants.f90'
integer*4 ihdr(11)
real*4 s(NSMAX)
real*4 ccfred(NSMAX)
logical*1 lstrong(0:1023)
integer*2 id2
complex c0
character*80 arg,infile
common/jt9com/ss(184,NSMAX),savg(NSMAX),c0(NDMAX),id2(NMAX),nutc,ndiskdat, &
common/jt9com/ss(184,NSMAX),savg(NSMAX),id2(NMAX),nutc,ndiskdat, &
ntr,mousefqso,newdat,nfa,nfb,ntol,kin,nzhsym,nsynced,ndecoded
common/tracer/limtrace,lu
@@ -39,11 +33,6 @@ program jt9
limtrace=0
lu=12
nfa=1000
nfb=2000
mousefqso=1500
newdat=1
ndiskdat=1
do ifile=ifile1,nargs
call getarg(ifile,infile)
@@ -82,10 +71,6 @@ program jt9
call timer('jt9 ',0)
endif
! do i=1,npts
! id2(i)=100.0*sin(6.283185307*1600.0*i/12000.0)
! enddo
id2=0 !??? Why is this necessary ???
do iblk=1,npts/kstep
@@ -99,8 +84,7 @@ program jt9
! Emit signal readyForFFT
ingain=0
call timer('symspec ',0)
call symspec(k,ntrperiod,nsps,ingain,pxdb,s,ccfred,df3, &
ihsym,nzap,slimit,lstrong,npts8)
call symspec(k,ntrperiod,nsps,ingain,slope,pxdb,s,df3,ihsym,npts8)
call timer('symspec ',1)
nhsym0=nhsym
if(ihsym.ge.173) go to 10
@@ -109,7 +93,7 @@ program jt9
10 close(10)
call fillcom(nutc0,ndepth)
call decoder(ss,c0,1)
call decoder(ss,id2)
enddo
call timer('jt9 ',1)
lib/jt9a.f90
+6 -1
View File
@@ -20,6 +20,7 @@ subroutine jt9a
open(12,file='timer.out',status='unknown')
limtrace=0
! limtrace=-1 !Disable all calls to timer()
lu=12
i1=attach_jt9()
@@ -44,12 +45,16 @@ subroutine jt9a
go to 999
endif
p_jt9=>address_jt9()
call timer('jt9b ',0)
call jt9b(p_jt9,nbytes)
call timer('jt9b ',1)
100 inquire(file=trim(cwd)//'/.lock',exist=fileExists)
if(fileExists) go to 10
call sleep_msec(100)
go to 100
999 return
999 call timer('jt9b ',101)
return
end subroutine jt9a
lib/jt9b.f90
+3 -5
View File
@@ -1,14 +1,12 @@
subroutine jt9b(jt9com,nbytes)
parameter (NTMAX=120)
parameter (NSMAX=1365)
include 'constants.f90'
integer*1 jt9com(0:nbytes-1)
kss=0
ksavg=kss + 4*184*NSMAX
kc0=ksavg + 4*NSMAX
kid2=kc0 + 2*4*NTMAX*1500
kid2=ksavg + 4*NSMAX
knutc=kid2 + 2*NTMAX*12000
call jt9c(jt9com(kss),jt9com(ksavg),jt9com(kc0),jt9com(kid2),jt9com(knutc))
call jt9c(jt9com(kss),jt9com(ksavg),jt9com(kid2),jt9com(knutc))
return
end subroutine jt9b
lib/jt9c.f90
+6 -6
View File
@@ -1,21 +1,21 @@
subroutine jt9c(ss,savg,c0,id2,nparams0)
subroutine jt9c(ss,savg,id2,nparams0)
parameter (NSMAX=22000)
include 'constants.f90'
real*4 ss(184*NSMAX),savg(NSMAX)
complex c0(1800*1500)
integer*2 id2(1800*12000)
integer*2 id2(NTMAX*12000)
integer nparams0(21),nparams(21)
character*20 datetime
common/npar/nutc,ndiskdat,ntrperiod,nfqso,newdat,npts8,nfa,nfb,ntol, &
kin,nzhsym,nsave,nagain,ndepth,nrxlog,nfsample,datetime
kin,nzhsym,nsave,nagain,ndepth,ntxmode,nmode,datetime
equivalence (nparams,nutc)
nutc=id2(1)+int(savg(1)) !Silence compiler warning
nparams=nparams0 !Copy parameters into common/npar/
call flush(6)
if(sum(nparams).ne.0) call decoder(ss,c0,0)
! if(sum(nparams).ne.0) call decoder(ss,id2)
call decoder(ss,id2)
return
end subroutine jt9c
lib/jt9test.f90
-183
View File
@@ -1,183 +0,0 @@
program jt9test
! Decoder for JT9. Can run stand-alone, reading data from *.wav files;
! or as the back end of wsjt-x, with data placed in a shared memory region.
! NB: For unknown reason, ***MUST*** be compiled by g95 with -O0 !!!
parameter (NTMAX=120)
parameter (NMAX=NTMAX*12000) !Total sample intervals per 30 minutes
parameter (NDMAX=NTMAX*1500) !Sample intervals at 1500 Hz rate
parameter (NSMAX=1365) !Max length of saved spectra
integer*4 ihdr(11)
real*4 s(NSMAX)
real*4 ccfred(NSMAX)
logical*1 lstrong(0:1023)
integer*1 i1SoftSymbols(207)
character*22 msg
character*33 line
character*80 arg,infile
integer*2 id2
complex c0
complex c1(0:2700000)
common/jt9com/ss(184,NSMAX),savg(NSMAX),c0(NDMAX),id2(NMAX),nutc,ndiskdat, &
ntr,mousefqso,newdat,nfa,nfb,ntol,kin,nzhsym,nsynced,ndecoded
common/tracer/limtrace,lu
nargs=iargc()
if(nargs.lt.1) then
print*,'Usage: jt9 TRperiod file1 [file2 ...]'
print*,' Reads data from *.wav files.'
print*,''
print*,' jt9 -s'
print*,' Gets data from shared memory region.'
go to 999
endif
call getarg(1,arg)
if(arg(1:2).eq.'-s') then
! call jt9a
! call ftnquit
go to 999
endif
read(arg,*) ntrperiod
ifile1=2
limtrace=0
lu=12
nfa=1000
nfb=2000
ntol=500
nfqso=1500
newdat=1
nb=0
nbslider=100
limit=20000
ndiskdat=1
do ifile=ifile1,nargs
call getarg(ifile,infile)
open(10,file=infile,access='stream',status='old',err=998)
read(10) ihdr
i1=index(infile,'.wav')
read(infile(i1-4:i1-1),*,err=1) nutc0
go to 2
1 nutc0=0
2 nsps=0
if(ntrperiod.eq.1) then
nsps=6912
nzhsym=181
else if(ntrperiod.eq.2) then
nsps=15360
nzhsym=178
else if(ntrperiod.eq.5) then
nsps=40960
nzhsym=172
else if(ntrperiod.eq.10) then
nsps=82944
nzhsym=171
else if(ntrperiod.eq.30) then
nsps=252000
nzhsym=167
endif
if(nsps.eq.0) stop 'Error: bad TRperiod'
kstep=nsps/2
tstep=kstep/12000.0
k=0
nhsym0=-999
npts=(60*ntrperiod-6)*12000
if(ifile.eq.ifile1) then
open(12,file='timer.out',status='unknown')
call timer('jt9 ',0)
endif
! do i=1,npts
! id2(i)=100.0*sin(6.283185307*1600.0*i/12000.0)
! enddo
do iblk=1,npts/kstep
k=iblk*kstep
call timer('read_wav',0)
read(10,end=10) id2(k-kstep+1:k)
call timer('read_wav',1)
nhsym=(k-2048)/kstep
if(nhsym.ge.1 .and. nhsym.ne.nhsym0) then
! Emit signal readyForFFT
ingain=0
call timer('symspec ',0)
call symspec(k,ntrperiod,nsps,ingain,nb,nbslider,pxdb, &
s,ccfred,df3,ihsym,nzap,slimit,lstrong,npts8)
call timer('symspec ',1)
nhsym0=nhsym
if(ihsym.ge.184) go to 10
endif
enddo
10 close(10)
nsps8=nsps/8
iz=1000.0/df3
nutc=nutc0
call timer('sync9 ',0)
call sync9(ss,nzhsym,tstep,df3,ccfred,ia,ib,ipk) !Get sync, freq
call timer('sync9 ',1)
fgood=0.
df8=1500.0/(nsps/8)
sbest=0.
do i=ia,ib
f=(i-1)*df3
if((i.eq.ipk .or. ccfred(i).ge.3.0) .and. f.gt.fgood+10.0*df8) then
call timer('test9 ',0)
fpk=1000.0 + df3*(i-1)
c1(0:npts8-1)=conjg(c0(1:npts8))
call test9(c1,npts8,nsps8,fpk,syncpk,snrdb,xdt,freq,drift, &
i1SoftSymbols)
call timer('test9 ',1)
call timer('decode9 ',0)
call decode9(i1SoftSymbols,limit,nlim,msg)
call timer('decode9 ',1)
snr=snrdb
sync=syncpk - 2.0
if(sync.lt.0.0) sync=0.0
nsync=sync
if(nsync.gt.10) nsync=10
nsnr=nint(snr)
width=0.0
if(sync.gt.sbest .and. fgood.eq.0.0) then
sbest=sync
write(line,1010) nutc,nsync,nsnr,xdt,1000.0+fpk,width
if(nsync.gt.0) nsynced=1
endif
if(msg.ne.' ') then
write(*,1010) nutc,nsync,nsnr,xdt,freq,drift,msg
1010 format(i4.4,i4,i5,f6.1,f8.2,f6.2,3x,a22)
fgood=f
nsynced=1
ndecoded=1
endif
endif
enddo
if(fgood.eq.0.0) then
write(*,1020) line
1020 format(a33)
endif
enddo
call timer('jt9 ',1)
call timer('jt9 ',101)
! call ftnquit
go to 999
998 print*,'Cannot open file:'
print*,infile
999 end program jt9test
lib/move.f90
+7
View File
@@ -0,0 +1,7 @@
subroutine move(x,y,n)
real x(n),y(n)
do i=1,n
y(i)=x(i)
enddo
return
end subroutine move
lib/peakdf9.f90
-32
View File
@@ -1,32 +0,0 @@
subroutine peakdf9(c0,npts8,nsps8,istart,foffset,idfpk)
complex c0(0:npts8-1)
complex zsum
include 'jt9sync.f90'
twopi=8.0*atan(1.0)
df=1500.0/nsps8
smax=0.
do idf=-5,5
f0=foffset + 0.1*df*idf
dphi=twopi*f0/1500.0
sum=0.
do j=1,16
i1=(ii(j)-1)*nsps8 + istart
phi=0.
zsum=0.
do i=i1,i1+nsps8-1
if(i.lt.0 .or. i.gt.npts8-1) cycle
phi=phi + dphi
zsum=zsum + c0(i) * cmplx(cos(phi),-sin(phi))
enddo
sum=sum + real(zsum)**2 + aimag(zsum)**2
enddo
if(sum.gt.smax) then
idfpk=idf
smax=sum
endif
enddo
return
end subroutine peakdf9
lib/peakdt9.f90
+2 -7
View File
@@ -9,11 +9,8 @@ subroutine peakdt9(c2,nz2,nsps8,nspsd,c3,nz3,xdt)
p=0.
i0=5*nspsd
do i=0,nz2-1
z=1.e-3*sum(c2(max(i-(nspsd-1),0):i)) !Integrate
p(i0+i)=real(z)**2 + aimag(z)**2 !Symbol power at freq=0
! Option here for coherent processing ?
! write(53,3301) i,z,p(i0+i),atan2(aimag(z),real(z))
!3301 format(i6,4e12.3)
z=1.e-3*sum(c2(max(i-(nspsd-1),0):i))
p(i0+i)=real(z)**2 + aimag(z)**2 !Integrated symbol power at freq=0
enddo
call getlags(nsps8,lag0,lag1,lag2)
@@ -35,8 +32,6 @@ subroutine peakdt9(c2,nz2,nsps8,nspsd,c3,nz3,xdt)
enddo
ss=(sum1/16.0)/(sum0/69.0) - 1.0
xdt=(lag-lag0)*dtlag
! write(52,3001) lag,xdt,ss
!3001 format(i5,2f12.3)
if(ss.gt.smax) then
smax=ss
lagpk=lag
lib/ptt_unix.c
+19 -27
View File
@@ -68,7 +68,7 @@ int lp_ptt (int fd, int onoff);
/* parport functions */
int dev_is_parport(int fd);
int ptt_parallel(int fd, int *ntx, int *iptt);
int ptt_parallel(int fd, int ntx, int *iptt);
int ptt_serial(int fd, int *ntx, int *iptt);
int fd=-1; /* Used for both serial and parallel */
@@ -91,53 +91,42 @@ int fd=-1; /* Used for both serial and parallel */
#define STATE_PORT_OPEN_PARALLEL 1
#define STATE_PORT_OPEN_SERIAL 2
//int ptt_(int *unused, char *ptt_port, int *ntx, int *iptt)
int ptt_(int *unused, int *ntx, int *iptt)
//int ptt_(int *unused, int *ntx, int *iptt)
int ptt_(int nport, int ntx, int *iptt, int *nopen)
{
static int state=0;
char *p;
char ptt_port[]="/dev/ttyUSB0";
fflush(stdout);
// ### Temporary:
char* ptt_port;
if(*unused != -99) {
*iptt=*ntx;
return 0;
}
// ###
/* In the very unlikely event of a NULL pointer, just return.
* Yes, I realise this should not be possible in WSJT.
*/
// In the very unlikely event of a NULL pointer, just return.
if (ptt_port == NULL) {
*iptt = *ntx;
*iptt = ntx;
return (0);
}
switch (state) {
case STATE_PORT_CLOSED:
/* Remove trailing ' ' */
// Remove trailing ' '
if ((p = strchr(ptt_port, ' ')) != NULL)
*p = '\0';
/* If all that is left is a '\0' then also just return */
// If all that is left is a '\0' then also just return
if (*ptt_port == '\0') {
*iptt = *ntx;
*iptt = ntx;
return(0);
}
if ((fd = open(ptt_port, O_RDWR|O_NONBLOCK)) < 0) {
fprintf(stderr, "Can't open %s.\n", ptt_port);
return (1);
}
if (dev_is_parport(fd)) {
state = STATE_PORT_OPEN_PARALLEL;
lp_reset(fd);
ptt_parallel(fd, ntx, iptt);
} else {
state = STATE_PORT_OPEN_SERIAL;
ptt_serial(fd, ntx, iptt);
ptt_serial(fd, &ntx, iptt);
}
break;
@@ -146,7 +135,7 @@ int ptt_(int *unused, int *ntx, int *iptt)
break;
case STATE_PORT_OPEN_SERIAL:
ptt_serial(fd, ntx, iptt);
ptt_serial(fd, &ntx, iptt);
break;
default:
@@ -155,9 +144,11 @@ int ptt_(int *unused, int *ntx, int *iptt)
state = STATE_PORT_CLOSED;
break;
}
*iptt=ntx;
return(0);
}
/*
* ptt_serial
*
@@ -180,6 +171,8 @@ ptt_serial(int fd, int *ntx, int *iptt)
ioctl(fd, TIOCMBIC, &control);
*iptt = 0;
}
printf("ptt_serial: %d %d",*ntx,*iptt);
fflush(stdout);
return(0);
}
@@ -379,10 +372,9 @@ lp_ptt (int fd, int onoff)
* iptt - pointer to fortran command status on or off
*/
int
ptt_parallel(int fd, int *ntx, int *iptt)
int ptt_parallel(int fd, int ntx, int *iptt)
{
if(*ntx) {
if(ntx) {
lp_ptt(fd, 1);
*iptt=1;
} else {
lib/redsync.f90
-44
View File
@@ -1,44 +0,0 @@
subroutine redsync(ss,ntrperiod,ihsym,iz,red)
! Compute the red curve (approx JT9 sync amplitude).
! NB: red() is used for real-time display only. A better ccfred() is
! computed during the decode procedure.
Parameter (NSMAX=1365)
real*4 ss(184,NSMAX)
real*4 red(NSMAX)
include 'jt9sync.f90'
lagmax=9
if(ntrperiod.eq.2) lagmax=5
if(ntrperiod.eq.5) lagmax=2
if(ntrperiod.eq.10) lagmax=1
if(ntrperiod.eq.30) lagmax=1
do i=1,iz !Loop over frequency range
smax=0.
do lag=-lagmax,lagmax !Loop over DT lags
sig=0.
do j=1,16
k=ii2(j)+lag
if(k.ge.5 .and. k.le.ihsym) then
sig=sig + ss(k,i) - 0.5*(ss(k-2,i)+ss(k-4,i))
endif
enddo
if(sig.gt.smax) smax=sig
enddo
red(i)=smax
enddo
call pctile(red,iz,40,xmed)
if(xmed.le.0.0) xmed=1.0
red=red/xmed
smax=0.
do i=1,iz
red(i)=0.2*db(red(i))
smax=max(smax,red(i))
enddo
h=10.
if(smax.gt.h) red=red*(h/smax)
return
end subroutine redsync
lib/s3avg.f90
+60
View File
@@ -0,0 +1,60 @@
subroutine s3avg(nsave,mode65,nutc,nhz,xdt,npol,ntol,s3,nsum,nkv,decoded)
! Save the current synchronized spectra, s3(64,63), for possible
! decoding of average.
real s3(64,63) !Synchronized spectra for 63 symbols
real s3a(64,63,64) !Saved spectra
real s3b(64,63) !Average spectra
integer iutc(64),ihz(64),ipol(64)
real dt(64)
character*22 decoded
logical ltext,first
data first/.true./
save
if(first) then
iutc=-1
ihz=0
ipol=0
first=.false.
ihzdiff=min(100,ntol)
dtdiff=0.2
endif
do i=1,64
if(nutc.eq.iutc(i) .and. abs(nhz-ihz(i)).lt.ihzdiff) then
nsave=mod(nsave-1+64,64)+1
go to 10
endif
enddo
iutc(nsave)=nutc !Save UTC
ihz(nsave)=nhz !Save freq in Hz
ipol(nsave)=npol !Save pol
dt(nsave)=xdt !Save DT
s3a(1:64,1:63,nsave)=s3 !Save the spectra
10 s3b=0.
do i=1,64 !Accumulate avg spectra
if(iutc(i).lt.0) cycle
if(mod(iutc(i),2).ne.mod(nutc,2)) cycle !Use only same sequence
if(abs(nhz-ihz(i)).gt.ihzdiff) cycle !Freq must match
if(abs(xdt-dt(i)).gt.dtdiff) cycle !DT must match
s3b=s3b + s3a(1:64,1:63,i)
nsum=nsum+1
enddo
decoded=' '
if(nsum.ge.2) then !Try decoding the sverage
nadd=mode65*nsum
call extract(s3b,nadd,ncount,nhist,decoded,ltext) !Extract the message
nkv=nsum
if(ncount.lt.0) then
nkv=0
decoded=' '
endif
endif
return
end subroutine s3avg
lib/setup65.f90
+96
View File
@@ -0,0 +1,96 @@
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
cycle
20 mref(k,2)=m+n
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
cycle
120 mref2(k,2)=m+n
enddo
return
end subroutine setup65
lib/smo121.f90
+13
View File
@@ -0,0 +1,13 @@
subroutine smo121(x,nz)
real x(nz)
x0=x(1)
do i=2,nz-1
x1=x(i)
x(i)=0.5*x(i) + 0.25*(x0+x(i+1))
x0=x1
enddo
return
end subroutine smo121
lib/softsym.f90
+4 -5
View File
@@ -1,9 +1,8 @@
subroutine softsym(c0,npts8,nsps8,newdat,fpk,syncpk,snrdb,xdt,freq,drift, &
schk,i1SoftSymbols)
subroutine softsym(id2,npts8,nsps8,newdat,fpk,syncpk,snrdb,xdt, &
freq,drift,schk,i1SoftSymbols)
! Compute the soft symbols
complex c0(0:npts8-1)
complex c2(0:4096-1)
complex c3(0:4096-1)
complex c5(0:4096-1)
@@ -16,7 +15,7 @@ subroutine softsym(c0,npts8,nsps8,newdat,fpk,syncpk,snrdb,xdt,freq,drift, &
ndown=nsps8/nspsd
! Mix, low-pass filter, and downsample to 16 samples per symbol
call downsam9(c0,npts8,nsps8,newdat,nspsd,fpk,c2,nz2)
call downsam9(id2,npts8,nsps8,newdat,nspsd,fpk,c2,nz2)
call peakdt9(c2,nz2,nsps8,nspsd,c3,nz3,xdt) !Find DT
@@ -26,7 +25,7 @@ subroutine softsym(c0,npts8,nsps8,newdat,fpk,syncpk,snrdb,xdt,freq,drift, &
freq=fpk - a(1)
drift=-2.0*a(2)
call twkfreq(c3,c5,nz3,fsample,a) !Correct for deltaF, fDot, fDDot
call twkfreq(c3,c5,nz3,fsample,a) !Correct for delta f, f1, f2 ==> a(1:3)
! Compute soft symbols (in scrambled order)
call symspec2(c5,nz3,nsps8,nspsd,fsample,freq,drift,snrdb,schk, &
lib/symspec.f90
+56 -85
View File
@@ -1,4 +1,4 @@
subroutine symspec(k,ntrperiod,nsps,ingain,pxdb,s,red,df3,ihsym,npts8)
subroutine symspec(k,ntrperiod,nsps,ingain,slope,pxdb,s,df3,ihsym,npts8)
! Input:
! k pointer to the most recent new data
@@ -10,111 +10,89 @@ subroutine symspec(k,ntrperiod,nsps,ingain,pxdb,s,red,df3,ihsym,npts8)
! Output:
! pxdb power (0-60 dB)
! s() spectrum for waterfall display
! red() first cut at JT9 sync amplitude
! s() current spectrum for waterfall display
! ihsym index number of this half-symbol (1-184)
parameter (NTMAX=120)
parameter (NMAX=NTMAX*12000) !Total sample intervals per 30 minutes
parameter (NDMAX=NTMAX*1500) !Sample intervals at 1500 Hz rate
parameter (NSMAX=1365) !Max length of saved spectra
parameter (NFFT1=1024)
parameter (NFFT2=1024,NFFT2A=NFFT2/8)
parameter (MAXFFT3=32768)
real*4 s(NSMAX),w3(MAXFFT3)
real*4 x1(NFFT1)
real*4 x2(NFFT1+105)
! jt9com
! ss() JT9 symbol spectra at half-symbol steps
! savg() average spectra for waterfall display
include 'constants.f90'
real*4 w3(MAXFFT3)
real*4 s(NSMAX)
real*4 scale(NSMAX)
real*4 ssum(NSMAX)
real*4 red(NSMAX)
complex cx(0:MAXFFT3-1)
real*4 xc(0:MAXFFT3-1)
complex cx(0:MAXFFT3/2)
integer*2 id2
complex c0
common/jt9com/ss(184,NSMAX),savg(NSMAX),c0(NDMAX),id2(NMAX),nutc,ndiskdat, &
common/jt9com/ss(184,NSMAX),savg(NSMAX),id2(NMAX),nutc,ndiskdat, &
ntr,mousefqso,newdat,nfa,nfb,ntol,kin,nzhsym,nsynced,ndecoded
data rms/999.0/,k0/99999999/,ntrperiod0/0/,nfft3z/0/
data rms/999.0/,k0/99999999/,nfft3z/0/,slope0/0.0/
equivalence (xc,cx)
save
if(ntrperiod.eq.1) nfft3=2048
if(ntrperiod.eq.2) nfft3=2048
if(ntrperiod.eq.5) nfft3=6144
if(ntrperiod.eq.10) nfft3=12288
if(ntrperiod.eq.30) nfft3=32768
jstep=nsps/16 !Step size = half-symbol in c0()
if(k.gt.NMAX) go to 999
if(k.lt.nfft3) then
if(ntrperiod.eq.-999) stop !Silence compiler warning
nfft3=16384 !df=12000.0/16384 = 0.732422
jstep=nsps/2 !Step size = half-symbol in id2()
if(k.gt.NMAX) go to 900
if(k.lt.2048) then !(2048 was nfft3) (Any need for this ???)
ihsym=0
go to 999 !Wait for enough samples to start
go to 900 !Wait for enough samples to start
endif
if(nfft3.ne.nfft3z) then !New nfft3, compute window
if(nfft3.ne.nfft3z .or. slope.ne.slope0) then !New nfft3, compute window
pi=4.0*atan(1.0)
if(ntrperiod.eq.1) then !Compute window for nfft3 spectrun
do i=1,nfft3
xx=float(i-1)/(nfft3-1)
w3(i)=0.40897 -0.5*cos(2.0*pi*xx) + 0.09103*cos(4.0*pi*xx)
! w3(i)=0.355768 - 0.487306*cos(2.0*pi*xx) + 0.144232*cos(4.0*pi*xx) &
! - 0.012604*cos(6.0*pi*xx)
enddo
else
do i=1,nfft3
w3(i)=2.0*(sin(i*pi/nfft3))**2 !Window for nfft3 spectrum
enddo
endif
do i=1,nfft3
w3(i)=2.0*(sin(i*pi/nfft3))**2 !Window for nfft3 spectrum
enddo
nfft3z=nfft3
nh=NSMAX/2
do i=1,NSMAX
x=slope*float(i)/nh - 1.0 + 2.6
scale(i)=10.0**x
enddo
slope0=slope
endif
if(k.lt.k0) then !Start a new data block
ja=0
ssum=0.
ihsym=0
k1=0
k8=0
x2=0.
if(ndiskdat.eq.0) then
id2(k+1:)=0
c0=0. !This is necessary to prevent "ghosts". Not sure why.
endif
if(ndiskdat.eq.0) id2(k+1:)=0 !Needed to prevent "ghosts". Not sure why.
endif
k0=k
kstep1=NFFT1
fac=2.0/NFFT1
nblks=(k-k1)/kstep1
gain=10.0**(0.05*ingain)
sq=0.
do nblk=1,nblks
do i=1,NFFT1
x1(i)=gain*id2(k1+i)
enddo
sq=sq + dot_product(x1,x1)
! Mix at 1500 Hz, lowpass at +/-750 Hz, and downsample to 1500 Hz complex.
x2(106:105+kstep1)=x1(1:kstep1)
call fil3(x2,kstep1+105,c0(k8+1),n2)
x2(1:105)=x1(kstep1-104:kstep1) !Save 105 trailing samples
k1=k1+kstep1
k8=k8+kstep1/8
do i=k0+1,k
x1=id2(i)
sq=sq + x1*x1
enddo
npts8=k8
ja=ja+jstep !Index of first sample
rms=sqrt(sq/(nblks*NFFT1))
sq=sq * gain**2
rms=sqrt(sq/(k-k0))
pxdb=0.
if(rms.gt.0.0) pxdb=20.0*log10(rms)
if(pxdb.gt.60.0) pxdb=60.0
k0=k
ja=ja+jstep !Index of first sample
fac0=0.1
do i=0,nfft3-1 !Copy data into cx
j=ja+i-(nfft3-1)
cx(i)=0.
if(j.ge.1 .and. j.le.NDMAX) cx(i)=c0(j)
xc(i)=0.
if(j.ge.1) xc(i)=fac0*id2(j)
enddo
if(ihsym.lt.184) ihsym=ihsym+1
cx(0:nfft3-1)=w3(1:nfft3)*cx(0:nfft3-1) !Apply window w3
call four2a(cx,nfft3,1,1,1) !Third FFT (forward)
xc(0:nfft3-1)=w3(1:nfft3)*xc(0:nfft3-1) !Apply window w3
call four2a(xc,nfft3,1,-1,0) !Real-to-complex FFT
n=min(184,ihsym)
df3=1500.0/nfft3 !JT9-a: 0.732 Hz = 0.42 * tone spacing
i0=nint(-500.0/df3)
iz=min(NSMAX,nint(1000.0/df3))
df3=12000.0/nfft3 !JT9-1: 0.732 Hz = 0.42 * tone spacing
! i0=nint(1000.0/df3)
i0=0
iz=min(NSMAX,nint(5000.0/df3))
fac=(1.0/nfft3)**2
do i=1,iz
j=i0+i-1
@@ -125,18 +103,11 @@ subroutine symspec(k,ntrperiod,nsps,ingain,pxdb,s,red,df3,ihsym,npts8)
s(i)=sx
enddo
999 continue
! s=0.05*s/ref
s=scale*s
savg=scale*ssum/ihsym
fac00=0.35
npct=20
call pctile(s,iz,npct,xmed0)
fac0=fac00/max(xmed0,0.006)
s(1:iz)=fac0*s(1:iz)
call pctile(ssum,iz,npct,xmed1)
fac1=fac00/max(xmed1,0.006*ihsym)
savg(1:iz)=fac1*ssum(1:iz)
! savg(iz+1:iz+20)=savg(iz)
call redsync(ss,ntrperiod,ihsym,iz,red)
900 npts8=k/8
return
end subroutine symspec
lib/symspec65.f90
+47
View File
@@ -0,0 +1,47 @@
subroutine symspec65(dd,npts,ss,nhsym,savg)
! Compute JT65 symbol spectra at half-symbol steps
parameter (NFFT=8192)
parameter (NSZ=3413) !NFFT*5000/12000
parameter (MAXHSYM=322)
real*8 hstep
real*4 dd(npts)
real*4 ss(MAXHSYM,NSZ)
real*4 savg(NSZ)
real*4 x(NFFT)
complex c(0:NFFT/2)
common/refspec/dfref,ref(NSZ)
equivalence (x,c)
save /refspec/
hstep=2048.d0*12000.d0/11025.d0 !half-symbol = 2229.116 samples
nsps=nint(2*hstep)
df=12000.0/NFFT
nhsym=npts/hstep - 1.0
savg=0.
fac1=1.e-3
do j=1,nhsym
i0=(j-1)*hstep
x(1:nsps)=fac1*dd(i0+1:i0+nsps)
x(nsps+1:)=0.
call four2a(c,NFFT,1,-1,0) !r2c forward FFT
do i=1,NSZ
s=real(c(i))**2 + aimag(c(i))**2
ss(j,i)=s
savg(i)=savg(i)+s
enddo
enddo
savg=savg/nhsym
call flat65(ss,nhsym,MAXHSYM,NSZ,ref) !Flatten the 2d spectrum, saving
dfref=df ! the reference spectrum ref()
savg=savg/ref
do j=1,nhsym
ss(j,1:NSZ)=ss(j,1:NSZ)/ref
enddo
return
end subroutine symspec65
lib/sync9.f90
+2 -1
View File
@@ -1,6 +1,7 @@
subroutine sync9(ss,nzhsym,lag1,lag2,ia,ib,ccfred,red2,ipkbest)
parameter (NSMAX=1365) !Max length of saved spectra
include 'constants.f90'
! parameter (NSMAX=1365) !Max length of saved spectra
real ss(184,NSMAX)
real ss1(184)
real ccfred(NSMAX)
lib/timer.f90
+5 -3
View File
@@ -7,7 +7,7 @@ subroutine timer(dname,k)
character*8 dname,name(50),space,ename
character*16 sname
logical on(50)
real ut(50),ut0(50),dut(50),tt(2)
real ut(50),ut0(50),dut(50)
integer ncall(50),nlevel(50),nparent(50)
integer onlevel(0:10)
common/tracer/ limtrace,lu
@@ -36,7 +36,8 @@ subroutine timer(dname,k)
if(on(n)) print*,'Error in timer: ',dname,' already on.'
level=level+1 !Increment the level
on(n)=.true.
ut0(n)=etime(tt)
call system_clock(icount,irate)
ut0(n)=float(icount)/irate
ncall(n)=ncall(n)+1
if(ncall(n).gt.1.and.nlevel(n).ne.level) then
nlevel(n)=-1
@@ -49,7 +50,8 @@ subroutine timer(dname,k)
else if(k.eq.1) then !Get stop times and accumulate sums. (k=1)
if(on(n)) then
on(n)=.false.
ut1=etime(tt)
call system_clock(icount,irate)
ut1=float(icount)/irate
ut(n)=ut(n)+ut1-ut0(n)
endif
level=level-1
lib/twkfreq.f90
+2 -2
View File
@@ -13,13 +13,13 @@ subroutine twkfreq(c3,c4,npts,fsample,a)
s=2.0/npts
do i=1,npts
x=s*(i-x0)
! if(mod(i,100).eq.1) then
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
endif
w=w*wstep
c4(i)=w*c3(i)
enddo
lib/twkfreq65.f90
+25
View File
@@ -0,0 +1,25 @@
subroutine twkfreq65(c4aa,n5,a)
complex c4aa(n5)
real a(5)
complex w,wstep
data twopi/6.283185307/
! Apply AFC corrections to the c4aa data
w=1.0
wstep=1.0
x0=0.5*(n5+1)
s=2.0/n5
do i=1,n5
x=s*(i-x0)
if(mod(i,100).eq.1) then
p2=1.5*x*x - 0.5
dphi=(a(1) + x*a(2) + p2*a(3)) * (twopi/1378.125)
wstep=cmplx(cos(dphi),sin(dphi))
endif
w=w*wstep
c4aa(i)=w*c4aa(i)
enddo
return
end subroutine twkfreq65
lib/wrapkarn.c
+70
View File
@@ -0,0 +1,70 @@
#include <math.h>
#include <stdio.h>
#include <float.h>
#include <limits.h>
#include <stdlib.h>
#include "rs.h"
static void *rs;
static int first=1;
void rs_encode_(int *dgen, int *sent)
// Encode JT65 data dgen[12], producing sent[63].
{
int dat1[12];
int b[51];
int i;
if(first) {
// Initialize the JT65 codec
rs=init_rs_int(6,0x43,3,1,51,0);
first=0;
}
// Reverse data order for the Karn codec.
for(i=0; i<12; i++) {
dat1[i]=dgen[11-i];
}
// Compute the parity symbols
encode_rs_int(rs,dat1,b);
// Move parity symbols and data into sent[] array, in reverse order.
for (i = 0; i < 51; i++) sent[50-i] = b[i];
for (i = 0; i < 12; i++) sent[i+51] = dat1[11-i];
}
void rs_decode_(int *recd0, int *era0, int *numera0, int *decoded, int *nerr)
// Decode JT65 received data recd0[63], producing decoded[12].
// Erasures are indicated in era0[numera]. The number of corrected
// errors is *nerr. If the data are uncorrectable, *nerr=-1 is returned.
{
int numera;
int i;
int era_pos[50];
int recd[63];
if(first) {
rs=init_rs_int(6,0x43,3,1,51,0);
first=0;
}
numera=*numera0;
for(i=0; i<12; i++) recd[i]=recd0[62-i];
for(i=0; i<51; i++) recd[12+i]=recd0[50-i];
if(numera)
for(i=0; i<numera; i++) era_pos[i]=era0[i];
*nerr=decode_rs_int(rs,recd,era_pos,numera);
for(i=0; i<12; i++) decoded[i]=recd[11-i];
}
void rs_encode__(int *dgen, int *sent)
{
rs_encode_(dgen, sent);
}
void rs_decode__(int *recd0, int *era0, int *numera0, int *decoded, int *nerr)
{
rs_decode_(recd0, era0, numera0, decoded, nerr);
}