Updates and additions to test progams.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7654 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Joe Taylor 2017-04-27 17:43:21 +00:00
parent 7fdad7b78b
commit 34b30ed664
13 changed files with 832 additions and 117 deletions

View File

@ -393,6 +393,7 @@ set (wsjt_FSRCS
lib/genmsk144.f90
lib/genmsk40.f90
lib/fsk4hf/genmskhf.f90
lib/fsk4hf/genwspr5.f90
lib/fsk4hf/genwsprlf.f90
lib/fsk4hf/getfc1.f90
lib/fsk4hf/getfc2.f90
@ -500,7 +501,10 @@ set (wsjt_FSRCS
lib/wav12.f90
lib/xcor.f90
lib/xcor4.f90
lib/wqdecode.f90
lib/wqencode.f90
lib/fsk4hf/wspr5d.f90
lib/fsk4hf/wspr5sim.f90
lib/fsk4hf/wsprlfsim.f90
lib/wspr_downsample.f90
lib/zplot9.f90
@ -865,7 +869,7 @@ if (Fortran_COMPILER_NAME MATCHES "gfortran.*")
set (CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -isysroot ${CMAKE_OSX_SYSROOT}")
endif (CMAKE_OSX_SYSROOT)
set (CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE} -fbounds-check -funroll-all-loops -fno-f2c ${General_FFLAGS}")
set (CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE} -funroll-all-loops -fno-f2c ${General_FFLAGS}")
set (CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -fbounds-check -fno-f2c ${General_FFLAGS}")
elseif (Fortran_COMPILER_NAME MATCHES "ifort.*")
# ifort (untested)
@ -1131,6 +1135,12 @@ target_link_libraries (ldpcsim168 wsjt_fort wsjt_cxx)
add_executable (fsk4hf lib/fsk4hf/fsk4hf.f90 wsjtx.rc)
target_link_libraries (fsk4hf wsjt_fort wsjt_cxx)
add_executable (wspr5d lib/fsk4hf/wspr5d.f90 wsjtx.rc)
target_link_libraries (wspr5d wsjt_fort wsjt_cxx)
add_executable (wspr5sim lib/fsk4hf/wspr5sim.f90 wsjtx.rc)
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)

View File

@ -1,15 +1,6 @@
subroutine cpolyfitw(c,pp,id,maxn,aa,bb,zz,nhardsync)
parameter (KK=60) !Information bits (50 + CRC10)
parameter (ND=300) !Data symbols: LDPC (300,60), r=1/5
parameter (NS=109) !Sync symbols (2 x 48 + Barker 13)
parameter (NR=3) !Ramp up/down
parameter (NN=NR+NS+ND) !Total symbols (412)
parameter (NSPS=16) !Samples per MSK symbol (16)
parameter (N2=2*NSPS) !Samples per OQPSK symbol (32)
parameter (N13=13*N2) !Samples in central sync vector (416)
parameter (NZ=NSPS*NN) !Samples in baseband waveform (6592)
parameter (NFFT1=4*NSPS,NH1=NFFT1/2)
include 'wsprlf_params.f90'
complex c(0:NZ-1) !Complex waveform
complex zz(NS+ND) !Complex symbol values (intermediate)
@ -22,7 +13,8 @@ subroutine cpolyfitw(c,pp,id,maxn,aa,bb,zz,nhardsync)
ib=NSPS-1
ib2=N2-1
n=0
do j=1,205 !First-pass demodulation
jz=(NS+ND+1)/2
do j=1,jz !First-pass demodulation
ia=ib+1
ib=ia+N2-1
zz(j)=sum(pp*c(ia:ib))/NSPS
@ -31,11 +23,11 @@ subroutine cpolyfitw(c,pp,id,maxn,aa,bb,zz,nhardsync)
x(n)=float(ia+ib)/NZ - 1.0
yi(n)=real(zz(j))*0.5*id(j)
yq(n)=aimag(zz(j))*0.5*id(j)
! write(54,1225) n,x(n),yi(n),yq(n)
!1225 format(i5,3f12.4)
write(54,1225) n,x(n),yi(n),yq(n)
1225 format(i5,3f12.4)
endif
if(j.le.204) then
zz(j+205)=sum(pp*c(ia+NSPS:ib+NSPS))/NSPS
if(j.lt.jz) then
zz(j+jz)=sum(pp*c(ia+NSPS:ib+NSPS))/NSPS
endif
enddo

107
lib/fsk4hf/genwspr5.f90 Normal file
View File

@ -0,0 +1,107 @@
subroutine genwspr5(msg,ichk,msgsent,itone,itype)
!Encode a WSPR-LF message, produce itone() array.
use crc
include 'wsprlf_params.f90'
character*22 msg,msgsent
character*60 cbits
integer*1,target :: idat(9)
integer*1 msgbits(KK),codeword(ND)
logical first
integer icw(ND)
integer id(NS+ND)
integer jd(NS+ND)
integer isync(48) !Long sync vector
integer ib13(13) !Barker 13 code
integer itone(NN)
integer*8 n8
data ib13/1,1,1,1,1,-1,-1,1,1,-1,1,-1,1/
data first/.true./
save first,isync
if(first) then
n8=z'cbf089223a51'
do i=1,48
isync(i)=-1
if(iand(n8,1).eq.1) isync(i)=1
n8=n8/2
enddo
first=.false.
endif
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
icw=2*codeword - 1 !NRZ codeword
! Message structure:
! I channel: R1 48*(S1+D1) S13 48*(D1+S1) R1
! Q channel: R1 D109 R1
! Generate QPSK with no offset, then shift the y array to get OQPSK.
! I channel:
n=0
k=0
do j=1,48 !Insert group of 48*(S1+D1)
n=n+1
id(n)=2*isync(j)
k=k+1
n=n+1
id(n)=icw(k)
enddo
do j=1,13 !Insert Barker 13 code
n=n+1
id(n)=2*ib13(j)
enddo
do j=1,48 !Insert group of 48*(S1+D1)
k=k+1
n=n+1
id(n)=icw(k)
n=n+1
id(n)=2*isync(j)
enddo
! Q channel
do j=1,204
k=k+1
n=n+1
id(n)=icw(k)
enddo
! Map I and Q to tones.
n=0
jz=(NS+ND+1)/2
do j=1,jz-1
jd(2*j-1)=id(j)/abs(id(j))
jd(2*j)=id(j+jz)/abs(id(j+jz))
enddo
jd(NS+ND)=id(jz)/abs(id(jz))
itone=0
do j=1,jz-1
itone(2*j-1)=(jd(2*j)*jd(2*j-1)+1)/2;
itone(2*j)=-(jd(2*j)*jd(2*j+1)-1)/2;
enddo
itone(NS+ND)=jd(NS+ND) !### Is this correct ??? ###
return
end subroutine genwspr5

View File

@ -2,16 +2,7 @@ subroutine genwsprlf(msgbits,id,icw,cbb,csync,itone)
!Encode a WSPR-LF message, produce baseband waveform and sync vector.
parameter (KK=60) !Information bits (50 + CRC10)
parameter (ND=300) !Data symbols: LDPC (300,60), r=1/5
parameter (NS=109) !Sync symbols (2 x 48 + Barker 13)
parameter (NR=3) !Ramp up/down (2 x half-length symbols)
parameter (NN=NR+NS+ND) !Total symbols (412)
parameter (NSPS=16) !Samples per MSK symbol (16)
parameter (N2=2*NSPS) !Samples per OQPSK symbol (32)
parameter (N13=13*N2) !Samples in central sync vector (416)
parameter (NZ=NSPS*NN) !Samples in baseband waveform (6560)
parameter (NFFT1=4*NSPS,NH1=NFFT1/2)
include 'wsprlf_params.f90'
complex cbb(0:NZ-1)
complex csync(0:NZ-1)
@ -129,19 +120,18 @@ subroutine genwsprlf(msgbits,id,icw,cbb,csync,itone)
! Map I and Q to tones.
n=0
do i=1,204
n=n+1
jd(n)=id(i)/abs(id(i))
n=n+1
jd(n)=id(i+205)/abs(id(i+205))
jz=(NS+ND+1)/2
do j=1,jz-1
jd(2*j-1)=id(j)/abs(id(j))
jd(2*j)=id(j+jz)/abs(id(j+jz))
enddo
jd(409)=id(205)/abs(id(205))
jd(NS+ND)=id(jz)/abs(id(jz))
itone=0
do i=1,204
itone(2*i-1)=(jd(2*i)*jd(2*i-1)+1)/2;
itone(2*i)=-(jd(2*i)*jd(2*i+1)-1)/2;
do j=1,jz-1
itone(2*j-1)=(jd(2*j)*jd(2*j-1)+1)/2;
itone(2*j)=-(jd(2*j)*jd(2*j+1)-1)/2;
enddo
itone(409)=jd(409) !### Is this correct ??? ###
itone(NS+ND)=jd(NS+ND) !### Is this correct ??? ###
return
end subroutine genwsprlf

View File

@ -1,15 +1,6 @@
subroutine getfc1w(c,fs,fc1)
parameter (KK=60) !Information bits (50 + CRC10)
parameter (ND=300) !Data symbols: LDPC (300,60), r=1/5
parameter (NS=109) !Sync symbols (2 x 48 + Barker 13)
parameter (NR=3) !Ramp up/down
parameter (NN=NR+NS+ND) !Total symbols (412)
parameter (NSPS=16) !Samples per MSK symbol (16)
parameter (N2=2*NSPS) !Samples per OQPSK symbol (32)
parameter (N13=13*N2) !Samples in central sync vector (416)
parameter (NZ=NSPS*NN) !Samples in baseband waveform (6592)
parameter (NFFT1=4*NSPS,NH1=NFFT1/2)
include 'wsprlf_params.f90'
complex c(0:NZ-1) !Complex waveform
complex c2(0:NFFT1-1) !Short spectra

View File

@ -1,15 +1,6 @@
subroutine getfc2w(c,csync,fs,fc1,fc2,fc3)
parameter (KK=60) !Information bits (50 + CRC10)
parameter (ND=300) !Data symbols: LDPC (300,60), r=1/5
parameter (NS=109) !Sync symbols (2 x 48 + Barker 13)
parameter (NR=3) !Ramp up/down
parameter (NN=NR+NS+ND) !Total symbols (412)
parameter (NSPS=16) !Samples per MSK symbol (16)
parameter (N2=2*NSPS) !Samples per OQPSK symbol (32)
parameter (N13=13*N2) !Samples in central sync vector (416)
parameter (NZ=NSPS*NN) !Samples in baseband waveform (6592)
parameter (NFFT1=4*NSPS,NH1=NFFT1/2)
include 'wsprlf_params.f90'
complex c(0:NZ-1) !Complex waveform
complex cs(0:NZ-1) !For computing spectrum

View File

@ -1,15 +1,6 @@
subroutine msksoftsymw(zz,aa,bb,id,nterms,ierror,rxdata,nhard0,nhardsync0)
parameter (KK=60) !Information bits (50 + CRC10)
parameter (ND=300) !Data symbols: LDPC (300,60), r=1/5
parameter (NS=109) !Sync symbols (2 x 48 + Barker 13)
parameter (NR=3) !Ramp up/down
parameter (NN=NR+NS+ND) !Total symbols (412)
parameter (NSPS=16) !Samples per MSK symbol (16)
parameter (N2=2*NSPS) !Samples per OQPSK symbol (32)
parameter (N13=13*N2) !Samples in central sync vector (416)
parameter (NZ=NSPS*NN) !Samples in baseband waveform (6592)
parameter (NFFT1=4*NSPS,NH1=NFFT1/2)
include 'wsprlf_params.f90'
complex zz(NS+ND) !Complex symbol values (intermediate)
complex z,z0
@ -21,8 +12,9 @@ subroutine msksoftsymw(zz,aa,bb,id,nterms,ierror,rxdata,nhard0,nhardsync0)
n=0
ierror=0
ierr=0
do j=1,205
xx=j*2.0/205.0 - 1.0
jz=(NS+ND+1)/2
do j=1,jz
xx=j*2.0/jz - 1.0
yii=1.
yqq=0.
if(nterms.gt.0) then
@ -51,12 +43,12 @@ subroutine msksoftsymw(zz,aa,bb,id,nterms,ierror,rxdata,nhard0,nhardsync0)
endif
nhard0=nhard0+ierr
endif
! write(41,3301) j,id(j),ierror(j),ierr,p*id(j)
!3301 format(4i6,f10.3)
! write(41,3301) j,id(j),ierror(j),ierr,n,p,p*id(j)
!3301 format(5i6,2f10.3)
enddo
do j=206,409
xx=(j-204.5)*2.0/205.0 - 1.0
do j=jz+1,NS+ND
xx=(j-jz+0.5)*2.0/jz - 1.0
yii=1.
yqq=0.
if(nterms.gt.0) then
@ -78,7 +70,7 @@ subroutine msksoftsymw(zz,aa,bb,id,nterms,ierror,rxdata,nhard0,nhardsync0)
ierror(j)=1
endif
nhard0=nhard0+ierr
! write(41,3301) j,id(j),ierror(j),ierr,p*id(j)
! write(41,3301) j,id(j),ierror(j),ierr,n,p,p*id(j)
enddo
return

View File

@ -1,15 +1,15 @@
subroutine watterson(c,nz,fs,delay,fspread)
subroutine watterson(c,npts,fs,delay,fspread)
parameter (NZMAX=6592)
complex c(0:nz-1)
include 'wsprlf_params.f90'
complex c(0:npts-1)
complex c2(0:NZMAX-1)
complex cs1(0:NZMAX-1)
complex cs2(0:NZMAX-1)
nonzero=0
df=fs/nz
df=fs/npts
if(fspread.gt.0.0) then
do i=0,nz-1
do i=0,npts-1
xx=gran()
yy=gran()
cs1(i)=0.707*cmplx(xx,yy)
@ -17,11 +17,11 @@ subroutine watterson(c,nz,fs,delay,fspread)
yy=gran()
cs2(i)=0.707*cmplx(xx,yy)
enddo
call four2a(cs1,nz,1,-1,1) !To freq domain
call four2a(cs2,nz,1,-1,1)
do i=0,nz-1
call four2a(cs1,npts,1,-1,1) !To freq domain
call four2a(cs2,npts,1,-1,1)
do i=0,npts-1
f=i*df
if(i.gt.nz/2) f=(i-nz)*df
if(i.gt.npts/2) f=(i-npts)*df
x=(f/(0.5*fspread))**2
a=0.
if(x.le.50.0) then
@ -37,16 +37,16 @@ subroutine watterson(c,nz,fs,delay,fspread)
!3101 format(f10.3,2f12.3,2f10.3)
endif
enddo
call four2a(cs1,nz,1,1,1) !Back to time domain
call four2a(cs2,nz,1,1,1)
cs1(0:nz-1)=cs1(0:nz-1)/nz
cs2(0:nz-1)=cs2(0:nz-1)/nz
call four2a(cs1,npts,1,1,1) !Back to time domain
call four2a(cs2,npts,1,1,1)
cs1(0:npts-1)=cs1(0:npts-1)/npts
cs2(0:npts-1)=cs2(0:npts-1)/npts
endif
nshift=nint(0.001*delay*fs)
c2(0:nz-1)=cshift(c(0:nz-1),nshift)
c2(0:npts-1)=cshift(c(0:npts-1),nshift)
sq=0.
do i=0,nz-1
do i=0,npts-1
if(nonzero.gt.1) then
c(i)=0.5*(cs1(i)*c(i) + cs2(i)*c2(i))
else
@ -56,7 +56,7 @@ subroutine watterson(c,nz,fs,delay,fspread)
! write(61,3001) i/12000.0,c(i)
!3001 format(3f12.6)
enddo
rms=sqrt(sq/nz)
rms=sqrt(sq/npts)
c=c/rms
return

162
lib/fsk4hf/wspr5d.f90 Normal file
View File

@ -0,0 +1,162 @@
program wspr5d
! Simulate characteristics of a potential "WSPR-LF" mode using LDPC (300,60)
! code, OQPSK modulation, and 5 minute T/R sequences.
! 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. Square, FFT; find peaks near -R/2 and +R/2 to get fc2
! 4. Mix from fc2 to 0
! 5. Fit cb13 (central part of csync) to c -> lag, phase
! 6. Fit complex ploynomial for channel equalization
! 7. Get soft bits from equalized data
include 'wsprlf_params.f90'
! Q: Would it be better for central Sync array to use both I and Q channels?
character arg*8,message*22,cbits*50
complex csync(0:NZ-1) !Sync symbols only, from cbb
complex c(0:NZ-1) !Complex waveform
complex c1(0:NZ-1) !Complex waveform
complex zz(NS+ND) !Complex symbol values (intermediate)
complex z
real rxdata(ND),llr(ND) !Soft symbols
real pp(2*NSPS) !Shaped pulse for OQPSK
real a(5) !For twkfreq1
real aa(20),bb(20) !Fitted polyco's
integer id(NS+ND) !NRZ values (+/-1) for Sync and Data
integer ierror(NS+ND)
integer isync(48) !Long sync vector
integer ib13(13) !Barker 13 code
integer*8 n8
integer*1 idat(7)
integer*1 decoded(KK),apmask(ND),cw(ND)
data ib13/1,1,1,1,1,-1,-1,1,1,-1,1,-1,1/
nargs=iargc()
if(nargs.ne.1) then
print*,'Usage: wspr5d maxn'
! print*,'Example: wsprlfsim 0 0 0 5 10 -20'
go to 999
endif
call getarg(1,arg)
read(arg,*) maxn
twopi=8.0*atan(1.0)
fs=NSPS*12000.0/NSPS0 !Sample rate
dt=1.0/fs !Sample interval (s)
tt=NSPS*dt !Duration of "itone" symbols (s)
ts=2*NSPS*dt !Duration of OQPSK symbols (s)
baud=1.0/tt !Keying rate for "itone" symbols (baud)
txt=NZ*dt !Transmission length (s)
do i=1,N2 !Half-sine pulse shape
pp(i)=sin(0.5*(i-1)*twopi/(2*NSPS))
enddo
n8=z'cbf089223a51'
do i=1,48
isync(i)=-1
if(iand(n8,1).eq.1) isync(i)=1
n8=n8/2
enddo
! Defind array id() for sync symbols
id=0
do j=1,48 !First group of 48
id(2*j-1)=2*isync(j)
enddo
do j=1,13 !Barker 13 code
id(j+96)=2*ib13(j)
enddo
do j=1,48 !Second group of 48
id(2*j+109)=2*isync(j)
enddo
do ifile=1,9999
read(10,end=999) c
! do i=0,NZ-1
! write(40,4001) i,c(i)
!4001 format(i8,2f10.6)
! enddo
call getfc1w(c,fs,fc1) !First approx for freq
call getfc2w(c,csync,fs,fc1,fc2,fc3) !Refined freq
!NB: Measured performance is about equally good using fc2 or fc3 here:
a(1)=-(fc1+fc2)
a(2:5)=0.
call twkfreq1(c,NZ,fs,a,c) !Mix c down by fc1+fc2
!---------------------------------------------------------------- DT
! Not presently used:
amax=0.
jpk=0
iaa=0
ibb=NZ-1
do j=-20*NSPS,20*NSPS,NSPS/8
ia=j
ib=NZ-1+j
if(ia.lt.0) then
ia=0
iaa=-j
else
iaa=0
endif
if(ib.gt.NZ-1) then
ib=NZ-1
ibb=NZ-1-j
endif
z=sum(c(ia:ib)*conjg(csync(iaa:ibb)))
if(abs(z).gt.amax) then
amax=abs(z)
jpk=j
endif
enddo
xdt=jpk/fs
!-----------------------------------------------------------------
nterms=maxn
c1=c
! do itry=1,1000
do itry=1,20
idf=itry/2
if(mod(itry,2).eq.0) idf=-idf
nhard0=0
nhardsync0=0
ifer=1
a(1)=idf*0.00085
a(2:5)=0.
call twkfreq1(c1,NZ,fs,a,c) !Mix c1 into c
call cpolyfitw(c,pp,id,maxn,aa,bb,zz,nhs)
call msksoftsymw(zz,aa,bb,id,nterms,ierror,rxdata,nhard0,nhardsync0)
if(nhardsync0.gt.35) cycle
rxav=sum(rxdata)/ND
rx2av=sum(rxdata*rxdata)/ND
rxsig=sqrt(rx2av-rxav*rxav)
rxdata=rxdata/rxsig
ss=0.84
llr=2.0*rxdata/(ss*ss)
apmask=0
max_iterations=40
ifer=0
call bpdecode300(llr,apmask,max_iterations,decoded,niterations,cw)
nbadcrc=0
if(niterations.ge.0) call chkcrc10(decoded,nbadcrc)
if(niterations.lt.0 .or. nbadcrc.ne.0) ifer=1
if(ifer.eq.0) exit
enddo !Freq dither loop
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)
endif
write(*,1110) nsnr,xdt,fc1+fc2,message
1110 format(i4,f7.2,f7.2,2x,a22)
enddo
999 end program wspr5d

74
lib/fsk4hf/wspr5sim.f90 Normal file
View File

@ -0,0 +1,74 @@
program wspr5sim
include 'wsprlf_params.f90'
character*12 arg
character*22 msg,msgsent
complex c0(0:NZ-1)
complex c(0:NZ-1)
integer itone(NN)
! Get command-line argument(s)
nargs=iargc()
if(nargs.ne.5) then
print*,'Usage: wspr5sim "message" f0 DT nfiles snr'
print*,'Example: wspr5sim "K1ABC FN42 30" 1500 0.0 10 -33'
go to 999
endif
call getarg(1,msg) !Get message from command line
call getarg(2,arg)
read(arg,*) f0
call getarg(3,arg)
read(arg,*) xdt
call getarg(4,arg)
read(arg,*) nfiles
call getarg(5,arg)
read(arg,*) snrdb
call genwspr5(msg,ichk,msgsent,itone,itype)
write(*,1000) f0,xdt,snr,nfiles,msgsent
1000 format('f0:',f9.3,' DT:',f6.2,' SNR:',f6.1,' nfiles:',i3,2x,a22)
twopi=8.0*atan(1.0)
fs=NSPS*12000.0/NSPS0 !Sample rate
dt=1.0/fs !Sample interval (s)
tt=NSPS*dt !Duration of "itone" symbols (s)
ts=2*NSPS*dt !Duration of OQPSK symbols (s)
baud=1.0/tt !Keying rate for "itone" symbols (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
dphi0=twopi*(f0-0.25d0*baud)*dt
dphi1=twopi*(f0+0.25d0*baud)*dt
phi=0.d0
c0=0.
k=-1
k=k+N2
do j=1,NN
dphi=dphi0
if(itone(j).eq.1) dphi=dphi1
if(k.eq.0) phi=-dphi
do i=1,NSPS
k=k+1
phi=phi+dphi
if(phi.gt.twopi) phi=phi-twopi
xphi=phi
if(k.ge.NZ) exit !### FIX THIS ###
c0(k)=cmplx(cos(xphi),sin(xphi))
enddo
enddo
c0=sig*c0 !Scale to requested sig level
do ifile=1,nfiles
if(snrdb.lt.90) then
do i=0,NZ-1 !Add gaussian noise at specified SNR
xnoise=gran()
ynoise=gran()
c(i)=c0(i) + cmplx(xnoise,ynoise)
enddo
endif
write(10) c
enddo
999 end program wspr5sim

View File

@ -12,23 +12,13 @@ program wsprlfsim
! 6. Fit complex ploynomial for channel equalization
! 7. Get soft bits from equalized data
parameter (KK=60) !Information bits (50 + CRC10)
parameter (ND=300) !Data symbols: LDPC (300,60), r=1/5
parameter (NS=109) !Sync symbols (2 x 48 + Barker 13)
parameter (NR=3) !Ramp up/down
parameter (NN=NR+NS+ND) !Total symbols (412)
parameter (NSPS=16) !Samples per MSK symbol (16)
parameter (N2=2*NSPS) !Samples per OQPSK symbol (32)
parameter (N13=13*N2) !Samples in central sync vector (416)
parameter (NZ=NSPS*NN) !Samples in baseband waveform (6592)
parameter (NFFT1=4*NSPS,NH1=NFFT1/2)
include 'wsprlf_params.f90'
! Q: Would it be better for central Sync array to use both I and Q channels?
character*8 arg
complex cbb(0:NZ-1) !Complex baseband waveform
complex csync(0:NZ-1) !Sync symbols only, from cbb
complex cb13(0:N13-1) !Barker 13 waveform
complex c(0:NZ-1) !Complex waveform
complex c0(0:NZ-1) !Complex waveform
complex c1(0:NZ-1) !Complex waveform
@ -36,11 +26,14 @@ program wsprlfsim
complex z
real xnoise(0:NZ-1) !Generated random noise
real ynoise(0:NZ-1) !Generated random noise
real s(0:NZ-1)
real rxdata(ND),llr(ND) !Soft symbols
real pp(2*NSPS) !Shaped pulse for OQPSK
real a(5) !For twkfreq1
real aa(20),bb(20) !Fitted polyco's
real t(11)
character*12 label(11)
integer*8 count0,count1,count2,count3,clkfreq
integer nc(11)
integer id(NS+ND) !NRZ values (+/-1) for Sync and Data
integer ierror(NS+ND)
integer icw(NN)
@ -49,6 +42,10 @@ program wsprlfsim
! integer*1 codeword(ND)
data msgbits/0,0,1,0,0,1,1,1,1,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,1,1,0,0,0,1, &
1,1,1,0,1,1,1,1,1,1,1,0,0,1,0,0,1,1,0,1,1,0,1,0,1,1,0,0,1,1/
data label/'genwsprlf','twkfreq1 a','watterson','noise gen','getfc1w', &
'getfc2w','twkfreq1 b','xdt loop','cpolyfitw','msksoftsym', &
'bpdecode300'/
nargs=iargc()
if(nargs.ne.6) then
print*,'Usage: wsprlfsim f0(Hz) delay(ms) fspread(Hz) maxn iters snr(dB)'
@ -69,30 +66,45 @@ program wsprlfsim
call getarg(6,arg)
read(arg,*) snrdb !Specified SNR_2500
nc=0
twopi=8.0*atan(1.0)
fs=12000.0/540.0 !Sample rate = 22.2222... Hz
fs=NSPS*12000.0/NSPS0 !Sample rate = 22.2222... Hz
dt=1.0/fs !Sample interval (s)
tt=NSPS*dt !Duration of "itone" symbols (s)
ts=2*NSPS*dt !Duration of OQPSK symbols (s)
baud=1.0/tt !Keying rate for "itone" symbols (baud)
txt=NZ*dt !Transmission length (s)
bandwidth_ratio=2500.0/(fs/2.0)
write(*,1000) f0,delay,fspread,maxn,iters,baud,3*baud,txt
1000 format('f0:',f5.1,' Delay:',f4.1,' fSpread:',f5.2,' maxn:',i3, &
' Iters:',i6/'Baud:',f7.3,' BW:',f5.1,' TxT:',f5.1,f5.2/)
write(*,1000) fs,f0,delay,fspread,maxn,baud,3*baud,txt,iters
1000 format('fs:',f10.3,' f0:',f5.1,' Delay:',f4.1,' fSpread:',f5.2, &
' maxn:',i3,/'Baud:',f8.3,' BW:',f5.1,' TxT:',f6.1,' iters:',i4/)
write(*,1004)
1004 format(/' SNR sync data ser ber fer fsigma tsigma'/60('-'))
1004 format(/' SNR sync data ser ber fer fsigma tsigma', &
' tsec'/68('-'))
do i=1,N2 !Half-sine pulse shape
pp(i)=sin(0.5*(i-1)*twopi/(2*NSPS))
enddo
t=0.
call system_clock(count0,clkfreq)
call genwsprlf(msgbits,id,icw,cbb,csync,itone)!Generate baseband waveform
cb13=csync(3088:3503) !Copy the Barker 13 waveform
call system_clock(count1,clkfreq)
t(1)=float(count1-count0)/float(clkfreq)
nc(1)=nc(1)+1
do i=0,NZ-1
write(40,4001) i,cbb(i),csync(i)
4001 format(i8,4f12.6)
enddo
call system_clock(count0,clkfreq)
a=0.
a(1)=f0
call twkfreq1(cbb,NZ,fs,a,c0) !Mix baseband to specified frequency
call system_clock(count1,clkfreq)
t(2)=float(count1-count0)/float(clkfreq)
nc(2)=nc(2)+1
isna=-20
isnb=-40
if(snrdb.ne.0.0) then
@ -108,11 +120,20 @@ program wsprlfsim
nfe=0
sqf=0.
sqt=0.
call system_clock(count2,clkfreq)
do iter=1,iters !Loop over requested iterations
c=c0
c=c0
call system_clock(count0,clkfreq)
if(delay.ne.0.0 .or. fspread.ne.0.0) then
call watterson(c,NZ,fs,delay,fspread)
endif
call system_clock(count1,clkfreq)
t(3)=t(3)+float(count1-count0)/float(clkfreq)
nc(3)=nc(3)+1
call system_clock(count0,clkfreq)
c=sig*c !Scale to requested SNR
if(snrdb.lt.90) then
do i=0,NZ-1 !Generate gaussian noise
@ -121,29 +142,45 @@ program wsprlfsim
enddo
c=c + cmplx(xnoise,ynoise) !Add AWGN noise
endif
call system_clock(count1,clkfreq)
t(4)=t(4)+float(count1-count0)/float(clkfreq)
nc(4)=nc(4)+1
call system_clock(count0,clkfreq)
call getfc1w(c,fs,fc1) !First approx for freq
call system_clock(count1,clkfreq)
t(5)=t(5)+float(count1-count0)/float(clkfreq)
nc(5)=nc(5)+1
call system_clock(count0,clkfreq)
call getfc2w(c,csync,fs,fc1,fc2,fc3) !Refined freq
call system_clock(count1,clkfreq)
t(6)=t(6)+float(count1-count0)/float(clkfreq)
nc(6)=nc(6)+1
sqf=sqf + (fc1+fc2-f0)**2
call system_clock(count0,clkfreq)
!NB: Measured performance is about equally good using fc2 or fc3 here:
a(1)=-(fc1+fc2)
a(2:5)=0.
call twkfreq1(c,NZ,fs,a,c) !Mix c down by fc1+fc2
call system_clock(count1,clkfreq)
t(7)=t(7)+float(count1-count0)/float(clkfreq)
nc(7)=nc(7)+1
! The following may not be necessary?
! z=sum(c(3088:3503)*cb13)/208.0 !Get phase from Barker 13 vector
! z0=z/abs(z)
! c=c*conjg(z0)
call system_clock(count0,clkfreq)
!---------------------------------------------------------------- DT
! Not presently used:
amax=0.
jpk=0
iaa=0
ibb=NZ-1
do j=-20*NSPS,20*NSPS !Get jpk
! z=sum(c(3088+j:3503+j)*conjg(cb13))/208.0
do j=-20*NSPS,20*NSPS,NSPS/8
ia=j
ib=NZ-1+j
if(ia.lt.0) then
@ -164,6 +201,10 @@ program wsprlfsim
enddo
xdt=jpk/fs
sqt=sqt + xdt**2
call system_clock(count1,clkfreq)
t(8)=t(8)+float(count1-count0)/float(clkfreq)
nc(8)=nc(8)+1
!-----------------------------------------------------------------
nterms=maxn
@ -176,9 +217,19 @@ program wsprlfsim
ifer=1
a(1)=idf*0.00085
a(2:5)=0.
call system_clock(count0,clkfreq)
call twkfreq1(c1,NZ,fs,a,c) !Mix c1 into c
call cpolyfitw(c,pp,id,maxn,aa,bb,zz,nhs)
call system_clock(count1,clkfreq)
t(9)=t(9)+float(count1-count0)/float(clkfreq)
nc(9)=nc(9)+1
call system_clock(count0,clkfreq)
call msksoftsymw(zz,aa,bb,id,nterms,ierror,rxdata,nhard0,nhardsync0)
call system_clock(count1,clkfreq)
t(10)=t(10)+float(count1-count0)/float(clkfreq)
nc(10)=nc(10)+1
if(nhardsync0.gt.35) cycle
rxav=sum(rxdata)/ND
rx2av=sum(rxdata*rxdata)/ND
@ -189,7 +240,11 @@ program wsprlfsim
apmask=0
max_iterations=40
ifer=0
call system_clock(count0,clkfreq)
call bpdecode300(llr,apmask,max_iterations,decoded,niterations,cw)
call system_clock(count1,clkfreq)
t(11)=t(11)+float(count1-count0)/float(clkfreq)
nc(11)=nc(11)+1
nbadcrc=0
if(niterations.ge.0) call chkcrc10(decoded,nbadcrc)
if(niterations.lt.0 .or. count(msgbits.ne.decoded).gt.0 .or. &
@ -203,16 +258,28 @@ program wsprlfsim
nhardsync0,nhard0,niterations,ifer,xdt
1045 format(f6.1,4i6,f8.2)
enddo
call system_clock(count3,clkfreq)
tsec=float(count3-count2)/float(clkfreq)
fsigma=sqrt(sqf/iters)
tsigma=sqrt(sqt/iters)
ser=float(nhardsync)/(NS*iters)
ber=float(nhard)/(ND*iters)
fer=float(nfe)/iters
write(*,1050) snrdb,nhardsync,nhard,ser,ber,fer,fsigma,tsigma
1050 format(f6.1,2i7,2f8.4,f7.3,2f8.2)
write(*,1050) snrdb,nhardsync,nhard,ser,ber,fer,fsigma,tsigma,tsec
1050 format(f6.1,2i7,2f8.4,f7.3,2f8.2f8.3)
enddo
write(*,1060) NS*iters,ND*iters
1060 format(60('-')/6x,2i7)
1060 format(68('-')/6x,2i7)
write(*,1065)
1065 format(/'Timing sec frac calls'/39('-'))
do i=1,11
write(*,1070) label(i),t(i),t(i)/sum(t),nc(i)
1070 format(a12,2f9.3,i8)
enddo
write(*,1072) sum(t),1.0
1072 format(39('-')/12x,2f10.3)
999 end program wsprlfsim

338
lib/wqdecode.f90 Normal file
View File

@ -0,0 +1,338 @@
!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
! File Name: wqdecode.f90
! Description:
!
! Copyright (C) 2001-2014 Joseph Taylor, K1JT
! License: GPL-3
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 3 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!-------------------------------------------------------------------------------
subroutine wqdecode(data0,message,ntype)
parameter (N15=32768)
integer*1 data0(11)
character*22 message
character*12 callsign
character*3 cdbm
character grid4*4,grid6*6
logical first
character*12 dcall(0:N15-1)
data first/.true./
save first,dcall
! May want to have a timeout (say, one hour?) on calls fetched
! from the hash table.
if(first) then
dcall=' '
first=.false.
endif
message=' '
call unpack50(data0,n1,n2)
! print*,data0,n1,n2
call unpackcall(n1,callsign)
i1=index(callsign,' ')
call unpackgrid(n2/128,grid4)
ntype=iand(n2,127) -64
! Standard WSPR message (types 0 3 7 10 13 17 ... 60)
if(ntype.ge.0 .and. ntype.le.62) then
nu=mod(ntype,10)
if(nu.eq.0 .or. nu.eq.3 .or. nu.eq.7) then
write(cdbm,'(i3)'),ntype
if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
message=callsign(1:i1)//grid4//' '//cdbm
call hash(callsign,i1-1,ih)
dcall(ih)=callsign(:i1)
else
nadd=nu
if(nu.gt.3) nadd=nu-3
if(nu.gt.7) nadd=nu-7
ng=n2/128 + 32768*(nadd-1)
call unpackpfx(ng,callsign)
ndbm=ntype-nadd
write(cdbm,'(i3)'),ndbm
if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
i2=index(callsign,' ')
message=callsign(:i2)//cdbm
call hash(callsign,i2-1,ih)
dcall(ih)=callsign(:i2)
endif
else if(ntype.lt.0) then
ndbm=-(ntype+1)
grid6=callsign(6:6)//callsign(1:5)
ih=(n2-ntype-64)/128
callsign=dcall(ih)
write(cdbm,'(i3)'),ndbm
if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
i2=index(callsign,' ')
if(dcall(ih)(1:1).ne.' ') then
message='<'//callsign(:i2-1)//'> '//grid6//' '//cdbm
else
message='<...> '//grid6//' '//cdbm
endif
endif
return
end subroutine wqdecode
!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
! File Name: unpack50.f90
! Description:
!
! Copyright (C) 2001-2014 Joseph Taylor, K1JT
! License: GPL-3
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 3 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!-------------------------------------------------------------------------------
subroutine unpack50(dat,n1,n2)
integer*1 dat(11)
i=dat(1)
i4=iand(i,255)
n1=ishft(i4,20)
i=dat(2)
i4=iand(i,255)
n1=n1 + ishft(i4,12)
i=dat(3)
i4=iand(i,255)
n1=n1 + ishft(i4,4)
i=dat(4)
i4=iand(i,255)
n1=n1 + iand(ishft(i4,-4),15)
n2=ishft(iand(i4,15),18)
i=dat(5)
i4=iand(i,255)
n2=n2 + ishft(i4,10)
i=dat(6)
i4=iand(i,255)
n2=n2 + ishft(i4,2)
i=dat(7)
i4=iand(i,255)
n2=n2 + iand(ishft(i4,-6),3)
return
end subroutine unpack50
!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
! File Name: unpackcall.f90
! Description:
!
! Copyright (C) 2001-2014 Joseph Taylor, K1JT
! License: GPL-3
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 3 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!-------------------------------------------------------------------------------
subroutine unpackcall(ncall,word)
character word*12,c*37
data c/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ '/
n=ncall
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:)
999 if(word(1:3).eq.'3D0') word='3DA0'//word(4:)
return
end subroutine unpackcall
!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
! File Name: unpackgrid.f90
! Description:
!
! Copyright (C) 2001-2014 Joseph Taylor, K1JT
! License: GPL-3
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 3 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!-------------------------------------------------------------------------------
subroutine unpackgrid(ng,grid)
parameter (NGBASE=180*180)
character grid*4,grid6*6,digit*10
data digit/'0123456789'/
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(1:4) !XXX explicitly truncate this -db
go to 100
10 n=ng-NGBASE-1
if(n.ge.1 .and.n.le.30) then
grid(1:1)='-'
grid(2:2)=char(48+n/10)
grid(3:3)=char(48+mod(n,10))
else if(n.ge.31 .and.n.le.60) then
n=n-30
grid(1:2)='R-'
grid(3:3)=char(48+n/10)
grid(4:4)=char(48+mod(n,10))
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
100 return
end subroutine unpackgrid
!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
! File Name: unpackpfx.f90
! Description:
!
! Copyright (C) 2001-2014 Joseph Taylor, K1JT
! License: GPL-3
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 3 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!-------------------------------------------------------------------------------
subroutine unpackpfx(ng,call1)
character*12 call1
character*3 pfx
if(ng.lt.60000) then
! Add-on prefix of 1 to 3 characters
n=ng
do i=3,1,-1
nc=mod(n,37)
if(nc.ge.0 .and. nc.le.9) then
pfx(i:i)=char(nc+48)
else if(nc.ge.10 .and. nc.le.35) then
pfx(i:i)=char(nc+55)
else
pfx(i:i)=' '
endif
n=n/37
enddo
call1=pfx//'/'//call1
if(call1(1:1).eq.' ') call1=call1(2:)
if(call1(1:1).eq.' ') call1=call1(2:)
else
! Add-on suffix, one or teo characters
i1=index(call1,' ')
nc=ng-60000
if(nc.ge.0 .and. nc.le.9) then
call1=call1(:i1-1)//'/'//char(nc+48)
else if(nc.ge.10 .and. nc.le.35) then
call1=call1(:i1-1)//'/'//char(nc+55)
else if(nc.ge.36 .and. nc.le.125) then
nc1=(nc-26)/10
nc2=mod(nc-26,10)
call1=call1(:i1-1)//'/'//char(nc1+48)//char(nc2+48)
endif
endif
return
end subroutine unpackpfx

View File

@ -122,8 +122,9 @@ program wsprcode
data0=0
call wqencode(msg,ntype0,data0) !Source encoding
write(*,1002) data0
1002 format(/'Source-encoded message (50 bits, hex):',7z3.2)
write(*,1002) data0(1:7),data0(1:6),data0(7)/64
1002 format(/'Source-encoded message, 50 bits:'/'Hex: ',7z3.2/ &
'Binary: ',6b9.8,b3.2)
call encode232(data0,nbytes,dat,MAXSYM) !Convolutional encoding
call inter_mept(dat,1) !Interleaving