Updated to some test routines.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7646 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Joe Taylor 2017-04-20 19:21:54 +00:00
parent 55fc0c7ced
commit d76a51b86f
5 changed files with 98 additions and 54 deletions

View File

@ -1,4 +1,4 @@
subroutine cpolyfit(c,pp,id,maxn,aa,bb,zz)
subroutine cpolyfit(c,pp,id,maxn,aa,bb,zz,nhardsync)
parameter (KK=84) !Information bits (72 + CRC12)
parameter (ND=168) !Data symbols: LDPC (168,84), r=1/2
@ -13,6 +13,7 @@ subroutine cpolyfit(c,pp,id,maxn,aa,bb,zz)
complex c(0:NZ-1) !Complex waveform
complex zz(NS+ND) !Complex symbol values (intermediate)
complex z,z0
real x(NS),yi(NS),yq(NS) !For complex polyfit
real pp(2*NSPS) !Shaped pulse for OQPSK
real aa(20),bb(20) !Fitted polyco's
@ -30,8 +31,8 @@ subroutine cpolyfit(c,pp,id,maxn,aa,bb,zz)
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.116) then
zz(j+117)=sum(pp*c(ia+NSPS:ib+NSPS))/NSPS
@ -41,21 +42,35 @@ subroutine cpolyfit(c,pp,id,maxn,aa,bb,zz)
aa=0.
bb=0.
nterms=0
chisqa=0.
chisqb=0.
if(maxn.gt.0) then
! Fit sync info with a complex polynomial
npts=n
mode=0
chisqa0=1.e30
chisqb0=1.e30
nterms=maxn
! do nterms=1,maxn
call polyfit4(x,yi,yi,npts,nterms,mode,aa,chisqa)
call polyfit4(x,yq,yq,npts,nterms,mode,bb,chisqb)
! if(chisqa/chisqa0.ge.0.98 .and. chisqb/chisqb0.ge.0.98) exit
chisqa0=chisqa
chisqb0=chisqb
! enddo
endif
nhardsync=0
do j=1,117
if(abs(id(j)).ne.2) cycle
xx=j*2.0/117.0 - 1.0
yii=1.
yqq=0.
if(nterms.gt.0) then
yii=aa(1)
yqq=bb(1)
do i=2,nterms
yii=yii + aa(i)*xx**(i-1)
yqq=yqq + bb(i)*xx**(i-1)
enddo
endif
z0=cmplx(yii,yqq)
z=zz(j)*conjg(z0)
p=real(z)
if(p*id(j).lt.0) nhardsync=nhardsync+1
enddo
return
end subroutine cpolyfit

View File

@ -43,7 +43,7 @@ subroutine getfc1(c,fc1)
ipk=i
fc1=f
endif
! write(51,3001) f,s(i),db(s(i))
! write(51,3001) f,s(i),db(s(i))
! 3001 format(f10.3,e12.3,f10.3)
enddo

View File

@ -1,4 +1,4 @@
subroutine getfc2(c,fc1,fc2)
subroutine getfc2(c,csync,fc1,fc2,fc3)
parameter (KK=84) !Information bits (72 + CRC12)
parameter (ND=168) !Data symbols: LDPC (168,84), r=1/2
@ -13,9 +13,11 @@ subroutine getfc2(c,fc1,fc2)
complex c(0:NZ-1) !Complex waveform
complex cs(0:NZ-1) !For computing spectrum
complex csync(0:NZ-1) !Sync symbols only, from cbb
real a(5)
fs=12000.0/72.0
df=fs/NZ
baud=fs/NSPS
a(1)=-fc1
a(2:5)=0.
@ -23,7 +25,6 @@ subroutine getfc2(c,fc1,fc2)
! Filter, square, then FFT to get refined carrier frequency fc2.
call four2a(cs,NZ,1,-1,1) !To freq domain
df=fs/NZ
ia=nint(0.75*baud/df)
cs(ia:NZ-1-ia)=0. !Save only freqs around fc1
call four2a(cs,NZ,1,1,1) !Back to time domain
@ -51,5 +52,23 @@ subroutine getfc2(c,fc1,fc2)
!1200 format(f10.3,2f15.3)
enddo
a(1)=-fc1
a(2:5)=0.
call twkfreq1(c,NZ,fs,a,cs) !Mix down by fc1
cs=cs*conjg(csync)
call four2a(cs,NZ,1,-1,1) !To freq domain
pmax=0.
do i=0,NZ-1
f=i*df
if(i.gt.NZ/2) f=(i-NZ)*df
p=real(cs(i))**2 + aimag(cs(i))**2
! write(51,3001) f,p,db(p)
!3001 format(f10.3,e12.3,f10.3)
if(p.gt.pmax) then
pmax=p
fc3=f
endif
enddo
return
end subroutine getfc2

View File

@ -28,6 +28,7 @@ program msksim
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 zz(NS+ND) !Complex symbol values (intermediate)
complex z
real xnoise(0:NZ-1) !Generated random noise
@ -53,7 +54,7 @@ program msksim
go to 999
endif
call getarg(1,arg)
read(arg,*) f0 !Generated carrier frequency
read(arg,*) f0 !Generated carrier frequency
call getarg(2,arg)
read(arg,*) delay !Delta_t (ms) for Watterson model
call getarg(3,arg)
@ -83,9 +84,8 @@ program msksim
pp(i)=sin(0.5*(i-1)*twopi/(2*NSPS))
enddo
call genmskhf(msgbits,id,icw,cbb,csync) !Generate baseband waveform and csync
call genmskhf(msgbits,id,icw,cbb,csync)!Generate baseband waveform and csync
cb13=csync(1680:2095) !Copy the Barker 13 waveform
a=0.
a(1)=f0
call twkfreq1(cbb,NZ,fs,a,cbb) !Mix to specified frequency
@ -105,25 +105,25 @@ program msksim
nfe=0
sqf=0.
do iter=1,iters !Loop over requested iterations
nhard0=0
nhardsync0=0
c=cbb
if(delay.ne.0.0 .or. fspread.ne.0.0) call watterson(c,fs,delay,fspread)
c=sig*c !Scale to requested SNR
c=sig*c !Scale to requested SNR
if(snrdb.lt.90) then
do i=0,NZ-1 !Generate gaussian noise
do i=0,NZ-1 !Generate gaussian noise
xnoise(i)=gran()
ynoise(i)=gran()
enddo
c=c + cmplx(xnoise,ynoise) !Add AWGN noise
c=c + cmplx(xnoise,ynoise) !Add AWGN noise
endif
call getfc1(c,fc1) !First approx for freq
call getfc2(c,fc1,fc2) !Refined freq
call getfc1(c,fc1) !First approx for freq
call getfc2(c,csync,fc1,fc2,fc3) !Refined freq
sqf=sqf + (fc1+fc2-f0)**2
!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 twkfreq1(c,NZ,fs,a,c) !Mix c down by fc1+fc2
! The following may not be necessary?
! z=sum(c(1680:2095)*cb13)/208.0 !Get phase from Barker 13 vector
@ -134,7 +134,7 @@ program msksim
! Not presently used:
amax=0.
jpk=0
do j=-20*NSPS,20*NSPS !Get jpk
do j=-20*NSPS,20*NSPS !Get jpk
z=sum(c(1680+j:2095+j)*cb13)/208.0
if(abs(z).gt.amax) then
amax=abs(z)
@ -142,28 +142,42 @@ program msksim
endif
enddo
xdt=jpk/fs
call cpolyfit(c,pp,id,maxn,aa,bb,zz)
nterms=maxn
call msksoftsym(zz,aa,bb,id,nterms,ierror,rxdata,nhard0,nhard, &
nhardsync0,nhardsync)
!---------------------------------------------------------------- Decode
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
call bpdecode168(llr,apmask,max_iterations,decoded,niterations,cw)
nbadcrc=0
ifer=0
if(niterations.ge.0) call chkcrc12(decoded,nbadcrc)
if(niterations.lt.0 .or. count(msgbits.ne.decoded).gt.0 .or. &
nbadcrc.ne.0) ifer=1
nterms=maxn
c0=c
do itry=1,10
idf=itry/2
if(mod(itry,2).eq.0) idf=-idf
nhard0=0
nhardsync0=0
ifer=1
a(1)=idf*0.01
a(2:5)=0.
call twkfreq1(c0,NZ,fs,a,c) !Mix c0 into c
call cpolyfit(c,pp,id,maxn,aa,bb,zz,nhs)
call msksoftsym(zz,aa,bb,id,nterms,ierror,rxdata,nhard0,nhardsync0)
if(nhardsync0.gt.12) 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 bpdecode168(llr,apmask,max_iterations,decoded,niterations,cw)
nbadcrc=0
if(niterations.ge.0) call chkcrc12(decoded,nbadcrc)
if(niterations.lt.0 .or. count(msgbits.ne.decoded).gt.0 .or. &
nbadcrc.ne.0) ifer=1
! if(ifer.eq.0) write(67,1301) snrdb,itry,idf,niterations, &
! nhardsync0,nhard0
!1301 format(f6.1,5i6)
if(ifer.eq.0) exit
enddo !Freq dither loop
nhard=nhard+nhard0
nhardsync=nharsdync+nhardsync0
nfe=nfe+ifer
enddo
@ -171,7 +185,7 @@ program msksim
ber=float(nhard)/((NS+ND)*iters)
fer=float(nfe)/iters
write(*,1050) snrdb,nhard,ber,fer,fsigma
write(60,1050) snrdb,nhard,ber,fer,fsigma
! write(60,1050) snrdb,nhard,ber,fer,fsigma
1050 format(f6.1,i7,f8.4,f7.3,f8.2)
enddo

View File

@ -1,5 +1,4 @@
subroutine msksoftsym(zz,aa,bb,id,nterms,ierror,rxdata,nhard0,nhard, &
nhardsync0,nhardsync)
subroutine msksoftsym(zz,aa,bb,id,nterms,ierror,rxdata,nhard0,nhardsync0)
parameter (KK=84) !Information bits (72 + CRC12)
parameter (ND=168) !Data symbols: LDPC (168,84), r=1/2
@ -38,7 +37,6 @@ subroutine msksoftsym(zz,aa,bb,id,nterms,ierror,rxdata,nhard0,nhard, &
p=real(z)
if(abs(id(j)).eq.2) then
if(real(z)*id(j).lt.0) then !Sync bit
nhardsync=nhardsync+1
nhardsync0=nhardsync0+1
ierror(j)=2
endif
@ -51,7 +49,6 @@ subroutine msksoftsym(zz,aa,bb,id,nterms,ierror,rxdata,nhard0,nhard, &
ierror(j)=1
endif
nhard0=nhard0+ierr
nhard=nhard+ierr
endif
enddo
@ -78,7 +75,6 @@ subroutine msksoftsym(zz,aa,bb,id,nterms,ierror,rxdata,nhard0,nhard, &
ierror(j)=1
endif
nhard0=nhard0+ierr
nhard=nhard+ierr
enddo
return