Add some new test programs.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7680 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Joe Taylor 2017-05-16 19:20:14 +00:00
parent 4731e02cdb
commit 28c09d4a76
10 changed files with 428 additions and 14 deletions

View File

@ -395,6 +395,7 @@ set (wsjt_FSRCS
lib/fsk4hf/genmskhf.f90
lib/fsk4hf/genwspr5.f90
lib/fsk4hf/genwsprlf.f90
lib/fsk4hf/genwspr_fsk8.f90
lib/fsk4hf/getfc1.f90
lib/fsk4hf/getfc2.f90
lib/fsk4hf/getfc1w.f90
@ -475,6 +476,7 @@ set (wsjt_FSRCS
lib/shell.f90
lib/fsk4hf/spec4.f90
lib/spec64.f90
lib/fsk4hf/spec8.f90
lib/spec9f.f90
lib/stdmsg.f90
lib/subtract65.f90
@ -506,7 +508,10 @@ set (wsjt_FSRCS
lib/wqencode.f90
lib/fsk4hf/wspr5d.f90
lib/fsk4hf/wspr5sim.f90
lib/fsk4hf/wspr_fsk8d.f90
lib/fsk4hf/wspr_fsk8_sim.f90
lib/fsk4hf/wspr5_downsample.f90
lib/fsk4hf/wspr_fsk8_downsample.f90
lib/fsk4hf/wspr5_wav.f90
lib/fsk4hf/wsprlfsim.f90
lib/wspr_downsample.f90
@ -1147,6 +1152,12 @@ target_link_libraries (wspr5sim wsjt_fort wsjt_cxx)
add_executable (wsprlfsim lib/fsk4hf/wsprlfsim.f90 wsjtx.rc)
target_link_libraries (wsprlfsim wsjt_fort wsjt_cxx)
add_executable (wspr_fsk8d lib/fsk4hf/wspr_fsk8d.f90 wsjtx.rc)
target_link_libraries (wspr_fsk8d wsjt_fort wsjt_cxx)
add_executable (wspr_fsk8_sim lib/fsk4hf/wspr_fsk8_sim.f90 wsjtx.rc)
target_link_libraries (wspr_fsk8_sim wsjt_fort wsjt_cxx)
add_executable (mskhfsim lib/fsk4hf/mskhfsim.f90 wsjtx.rc)
target_link_libraries (mskhfsim wsjt_fort wsjt_cxx)

View File

@ -0,0 +1,45 @@
subroutine genwspr_fsk8(msg,msgsent,itone)
! Encode a WSPR-LF 8-FSK message, producing array itone().
use crc
include 'wspr_fsk8_params.f90'
character*22 msg,msgsent
character*60 cbits
integer*1,target :: idat(9)
integer*1 msgbits(KK),codeword(3*ND)
integer itone(NN)
integer icos7(0:6)
data icos7/2,5,6,0,4,1,3/ !Costas 7x7 tone pattern
idat=0
call wqencode(msg,ntype0,idat) !Source encoding
id7=idat(7)
if(id7.lt.0) id7=id7+256
id7=id7/64
icrc=crc10(c_loc(idat),9) !Compute the 10-bit CRC
idat(8)=icrc/256 !Insert CRC into idat(8:9)
idat(9)=iand(icrc,255)
call wqdecode(idat,msgsent,itype)
write(cbits,1004) idat(1:6),id7,icrc
1004 format(6b8.8,b2.2,b10.10)
read(cbits,1006) msgbits
1006 format(60i1)
! call chkcrc10(msgbits,nbadcrc)
! print*,msgsent,itype,crc10_check(c_loc(idat),9),nbadcrc
call encode300(msgbits,codeword) !Encode the test message
! Message structure: S7 D100 S7
itone(1:7)=icos7
itone(NN-6:NN)=icos7
do j=1,ND
i=3*j -2
itone(j+7)=codeword(i)*4 + codeword(i+1)*2 + codeword(i+2)
enddo
return
end subroutine genwspr_fsk8

View File

@ -158,21 +158,24 @@ return
end subroutine mrbencode
subroutine nextpat(mi,k,iorder,iflag)
integer*1 mi(k),ms(k)
integer*1 mi(k),ms(k)
! generate the next test error pattern
ind=-1
do i=1,k-1
if( mi(i).eq.0 .and. mi(i+1).eq.1) ind=i
if( mi(i).eq.0 .and. mi(i+1).eq.1) ind=i
enddo
ms=0
ms(1:ind-1)=mi(1:ind-1)
ms(ind)=1
ms(ind+1)=0
if( ind+1 .lt. k ) then
nz=iorder-sum(ms)
ms(k-nz+1:k)=1
endif
mi=ms
if(ind.ge.1 .and. ind.le.k) then !### JHT temporary ### Correct ???
ms=0
ms(1:ind-1)=mi(1:ind-1)
ms(ind)=1
ms(ind+1)=0
if( ind+1 .lt. k ) then
nz=iorder-sum(ms)
ms(k-nz+1:k)=1
endif
mi=ms
endif !### JHT temporary ###
iflag=ind
return
end subroutine nextpat

31
lib/fsk4hf/spec8.f90 Normal file
View File

@ -0,0 +1,31 @@
subroutine spec8(c,s,savg)
include 'wspr_fsk8_params.f90'
complex c(0:NZ-1)
complex c1(0:NSPS-1)
real s(0:NH2,NN)
real savg(0:NH2)
fs=12000.0/NDOWN
df=fs/NSPS
savg=0.
do j=1,NN
ia=(j-1)*NSPS
ib=ia + NSPS-1
c1(0:NSPS-1)=c(ia:ib)
c1(NSPS:)=0.
call four2a(c1,NSPS,1,-1,1)
do i=0,NH2
s(i,j)=real(c1(i))**2 + aimag(c1(i))**2
enddo
savg=savg+s(0:NH2,j)
enddo
s=s/NZ
savg=savg/(NN*NZ)
do i=0,NH2
write(31,3101) i*df,savg(i)
3101 format(f10.3,f12.3)
enddo
return
end subroutine spec8

View File

@ -120,7 +120,7 @@ program wspr5d
go to 999
endif
close(10)
fa=100.0
fa=102.0
fb=150.0
call getfc1w(c,fs,fa,fb,fc1,xsnr) !First approx for freq
call getfc2w(c,csync,fs,fc1,fc2,fc3) !Refined freq
@ -215,5 +215,3 @@ program wspr5d
1120 format("<DecodeFinished>")
999 end program wspr5d
include 'wspr5_downsample.f90'

View File

@ -0,0 +1,29 @@
subroutine wspr_fsk8_downsample(iwave,c)
! Input: i*2 data in iwave() at sample rate 12000 Hz
! Output: Complex data in c(), sampled at 4=500 Hz
include 'wspr_fsk8_params.f90'
parameter (NMAX=240*12000,NFFTD=NMAX/24)
integer*2 iwave(NMAX)
complex c(0:NZ-1)
complex c1(0:NFFTD-1)
complex cx(0:NMAX/2)
real x(NMAX)
equivalence (x,cx)
df=12000.0/NMAX
x=iwave
call four2a(x,NMAX,1,-1,0) !r2c FFT to freq domain
i0=nint(1500.0/df)
c1(0)=cx(i0)
do i=1,NFFTD/2
c1(i)=cx(i0+i)
c1(NFFTD-i)=cx(i0-i)
enddo
c1=c1/NFFTD
call four2a(c1,NFFTD,1,1,1) !c2c FFT back to time domain
c=c1(0:NZ-1)
return
end subroutine wspr_fsk8_downsample

View File

@ -0,0 +1,13 @@
! LDPC (300,60) code
parameter (NDOWN=24) !Downsample factor
parameter (KK=60) !Information bits (50 + CRC10)
parameter (ND=100) !Data symbols
parameter (NS=14) !Sync symbols (2 @ Costas 7x7)
parameter (NN=NS+ND) !Total symbols (114)
parameter (NSPS0=24576) !Samples per symbol at 12000 S/s
parameter (NSPS=NSPS0/NDOWN) !Sam/sym, downsampled (1024)
parameter (N7=7*NSPS) !Samples in Costas 7x7 array (7168)
parameter (NZ=NSPS*NN) !Samples in downsampled waveform (116,736)
parameter (NZMAX=NSPS0*NN) !Samples in *.wav (2,801,664)
parameter (NFFT1=4*NSPS,NH1=NFFT1/2)
parameter (NH2=NSPS/2)

View File

@ -0,0 +1,64 @@
program wspr_fsk8_sim
! Generate simulated data for a 4-minute "WSPR-LF" mode using 8-FSK.
! Output is saved to a *.wav file.
use wavhdr
include 'wspr_fsk8_params.f90' !Set various constants
parameter (NMAX=300*12000)
type(hdr) h !Header for .wav file
character arg*12,fname*16
character msg*22,msgsent*22
integer itone(NN)
integer*2 iwave(NMAX) !Generated full-length waveform
! Get command-line argument(s)
nargs=iargc()
if(nargs.ne.5) then
print*,'Usage: wspr_fsk8_sim "message" f0 DT nfiles snr'
print*,'Example: wspr_fsk8_sim "K1ABC FN42 30" 1640 0.0 10 -33'
go to 999
endif
call getarg(1,msg) !Message to be transmitted
call getarg(2,arg)
read(arg,*) f0 !Freq relative to WSPR-band center (Hz)
call getarg(3,arg)
read(arg,*) xdt !Time offset from nominal (s)
call getarg(4,arg)
read(arg,*) nfiles !Number of files
call getarg(5,arg)
read(arg,*) snrdb !SNR_2500
twopi=8.0*atan(1.0)
fs=12000.0/NDOWN !Sample rate after downsampling
dt=1.0/fs !Sample interval (s)
tt=NSPS*dt !Duration of symbols (s)
baud=1.0/tt !Keying rate
bw=8*baud
txt=NZ*dt !Transmission length (s)
bandwidth_ratio=2500.0/(fs/2.0)
sig=sqrt(bandwidth_ratio) * 10.0**(0.05*snrdb)
if(snrdb.gt.90.0) sig=1.0
txt=NN*NSPS0/12000.0
call genwspr_fsk8(msg,msgsent,itone) !Encode the message, get itone
write(*,1000) f0,xdt,txt,snrdb,bw,msgsent
1000 format('f0:',f9.3,' DT:',f6.2,' TxT:',f6.1,' SNR:',f6.1, &
' BW:',f4.1,2x,a22)
do ifile=1,nfiles
call wspr_fsk8_wav(baud,xdt,f0,itone,snrdb,iwave)
h=default_header(12000,NMAX)
write(fname,1102) ifile
1102 format('000000_',i4.4,'.wav')
open(10,file=fname,status='unknown',access='stream')
write(10) h,iwave !Save to *.wav file
close(10)
write(*,1110) ifile,xdt,f0,snrdb,fname
1110 format(i4,f7.2,f8.2,f7.1,2x,a16)
enddo
999 end program wspr_fsk8_sim
include 'wspr_fsk8_wav.f90'

View File

@ -0,0 +1,45 @@
subroutine wspr_fsk8_wav(baud,xdt,f0,itone,snrdb,iwave)
! Generate iwave() from itone().
include 'wspr_fsk8_params.f90'
parameter (NMAX=240*12000)
integer itone(NN)
integer*2 iwave(NMAX)
real*8 twopi,dt,dphi,phi
real dat(NMAX)
twopi=8.d0*atan(1.d0)
dt=1.d0/12000.d0
dat=0.
if(snrdb.lt.90) then
do i=1,NMAX
dat(i)=gran() !Generate gaussian noise
enddo
bandwidth_ratio=2500.0/6000.0
sig=sqrt(2*bandwidth_ratio)*10.0**(0.05*snrdb)
else
sig=1.0
endif
phi=0.d0
k=nint(xdt/dt)
do j=1,NN
dphi=twopi*(f0+ itone(j)*baud)*dt
if(k.eq.0) phi=-dphi
do i=1,NSPS0
k=k+1
phi=phi+dphi
if(phi.gt.twopi) phi=phi-twopi
xphi=phi
if(k.gt.0 .and. k.le.NMAX) dat(k)=dat(k) + sig*sin(xphi)
enddo
enddo
fac=32767.0
rms=100.0
if(snrdb.ge.90.0) iwave=nint(fac*dat)
if(snrdb.lt.90.0) iwave=nint(rms*dat)
return
end subroutine wspr_fsk8_wav

175
lib/fsk4hf/wspr_fsk8d.f90 Normal file
View File

@ -0,0 +1,175 @@
program wspr_fsk8d
! WSPR-LF is a potential WSPR-like mode intended for use at LF and MF.
! This version uses 4-minute T/R sequences, an LDPC (300,60) code,
! 8-FSK modulation, and noncoherent demodulation. This decoder reads
! data from *.wav files.
! Reception and Demodulation algorithm:
! 1. Compute coarse spectrum; find fc1 = approx carrier freq
! 2. Mix from fc1 to 0; LPF at +/- 0.75*R
! 3. Find two 7x7 Costas arrays to get xdt and fc2
! 4. Mix from fc2 to 0, compute aligned symbol spectra
! 5. Get soft bits from symbol spectra
! Still to do: find and decode more than one signal in the specified passband.
include 'wspr_fsk8_params.f90'
parameter (NMAX=240*12000)
character arg*8,message*22,cbits*50,infile*80,datetime*11
character*120 data_dir
complex csync(0:N7-1) !Sync symbols for Costas 7x7 array
complex c1(0:2*N7-1)
complex c2(0:2*N7-1)
complex c(0:NMAX-1) !Complex waveform
real*8 fMHz
real rxdata(3*ND),llr(3*ND) !Soft symbols
real a(5) !For twkfreq1
real s(0:NH2,NN)
real savg(0:NH2)
real ss(0:N7)
real ps(0:7)
integer ihdr(11)
integer*2 iwave(NMAX) !Generated full-length waveform
integer*1 idat(7)
integer*1 decoded(KK),apmask(3*ND),cw(3*ND)
integer icos7(0:6)
data icos7/2,5,6,0,4,1,3/ !Costas 7x7 tone pattern
nargs=iargc()
if(nargs.lt.2) then
print*,'Usage: wspr_fsk8d [-a <data_dir>] [-f fMHz] file1 [file2 ...]'
go to 999
endif
iarg=1
data_dir="."
call getarg(iarg,arg)
if(arg(1:2).eq.'-a') then
call getarg(iarg+1,data_dir)
iarg=iarg+2
endif
call getarg(iarg,arg)
if(arg(1:2).eq.'-f') then
call getarg(iarg+1,arg)
read(arg,*) fMHz
iarg=iarg+2
endif
open(13,file=trim(data_dir)//'/ALL_WSPR.TXT',status='unknown', &
position='append')
twopi=8.0*atan(1.0)
fs=NSPS*12000.0/NSPS0 !Sample rate after downsampling (Hz)
dt=1.0/fs !Sample interval (s)
ts=NSPS*dt !Symbol duration (s)
baud=1.0/ts !Keying rate (Hz)
txt=NZ*dt !Transmission length (s)
phi=0.
k=-1
do j=0,6
dphi=twopi*baud*icos7(j)*dt
do i=1,NSPS
phi=phi+dphi
if(phi.gt.twopi) phi=phi-twopi
k=k+1
csync(k)=cmplx(cos(phi),sin(phi))
enddo
enddo
do ifile=iarg,nargs
call getarg(ifile,infile)
open(10,file=infile,status='old',access='stream')
read(10,end=999) ihdr,iwave
close(10)
j2=index(infile,'.wav')
read(infile(j2-4:j2-1),*) nutc
datetime=infile(j2-11:j2-1)
call wspr_fsk8_downsample(iwave,c)
c(NZ:)=0.
j0=0
! Need to get jpk ==> xdt
j0b=j0+107*NSPS
c1(0:N7-1)=c(j0:j0+N7-1)*conjg(csync)
c1(N7:)=0.
c2(0:N7-1)=c(j0b:j0b+N7-1)*conjg(csync)
c2(N7:)=0.
call four2a(c1,2*N7,1,-1,1)
call four2a(c2,2*N7,1,-1,1)
pmax=0.
df1=fs/(2*N7)
do i=0,N7
f=1500.0 + i*df1
p=1.e-9*(real(c1(i))**2 + aimag(c1(i))**2 + &
real(c2(i))**2 + aimag(c2(i))**2)
ss(i)=p
if(p.gt.pmax) then
pmax=p
ipk=i
endif
write(32,3201) f,ss(i)
3201 format(f12.4,e12.4)
enddo
sp3n=(ss(ipk-1)+ss(ipk)+ss(ipk+1)) !Sig + 3*noise
call pctile(ss,N7,45,base)
psig=sp3n-3*base !Sig only
pnoise=(2500.0/df1)*base !Noise in 2500 Hz
xsnr=db(psig/pnoise) !SNR from sync tones
f0=ipk*df1
a(1)=-f0
a(2:5)=0.
call twkfreq1(c,NZ,fs,a,c) !Mix from f0 down to 0
call spec8(c,s,savg) !Get symbol spectra
j0=0
do j=1,ND
k=j+7
ps=s(0:7,k)
! ps=sqrt(ps) !### ??? ###
r1=max(ps(1),ps(3),ps(5),ps(7))-max(ps(0),ps(2),ps(4),ps(6))
r2=max(ps(2),ps(3),ps(6),ps(7))-max(ps(0),ps(1),ps(4),ps(5))
r4=max(ps(4),ps(5),ps(6),ps(7))-max(ps(0),ps(1),ps(2),ps(3))
rxdata(3*j-2)=r4
rxdata(3*j-1)=r2
rxdata(3*j)=r1
enddo
rxav=sum(rxdata)/ND
rx2av=sum(rxdata*rxdata)/ND
rxsig=sqrt(rx2av-rxav*rxav)
rxdata=rxdata/rxsig
s0=0.84
llr=2.0*rxdata/(s0*s0)
apmask=0
max_iterations=40
ifer=0
call bpdecode300(llr,apmask,max_iterations,decoded,niterations,cw)
if(niterations.lt.0) call osd300(llr,4,decoded,niterations,cw)
nbadcrc=0
if(niterations.ge.0) call chkcrc10(decoded,nbadcrc)
if(niterations.lt.0 .or. nbadcrc.ne.0) ifer=1
nsnr=nint(xsnr)
freq=fMHz + 1.d-6*f0
nfdot=0
message=' '
if(ifer.eq.0) then
write(cbits,1100) decoded(1:50)
1100 format(50i1)
read(cbits,1102) idat
1102 format(6b8,b2)
idat(7)=ishft(idat(7),6)
call wqdecode(idat,message,itype)
write(*,1112) datetime(8:11),nsnr,xdt,freq,nfdot,message
1112 format(a4,i4,f5.1,f11.6,i3,2x,a22)
endif
write(13,1110) datetime,0,nsnr,xdt,freq,message,nfdot
1110 format(a11,2i4,f6.2,f12.7,2x,a22,i3)
enddo ! ifile loop
write(*,1120)
1120 format("<DecodeFinished>")
999 end program wspr_fsk8d