From 14cbb51e479b5155d7befb2025ae4832dc285a06 Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Thu, 20 Apr 2017 19:21:54 +0000 Subject: [PATCH] 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 --- lib/fsk4hf/cpolyfit.f90 | 39 +++++++++++++------ lib/fsk4hf/getfc1.f90 | 2 +- lib/fsk4hf/getfc2.f90 | 23 ++++++++++- lib/fsk4hf/mskhfsim.f90 | 82 +++++++++++++++++++++++---------------- lib/fsk4hf/msksoftsym.f90 | 6 +-- 5 files changed, 98 insertions(+), 54 deletions(-) diff --git a/lib/fsk4hf/cpolyfit.f90 b/lib/fsk4hf/cpolyfit.f90 index ccc2bc317..fd0d5987f 100644 --- a/lib/fsk4hf/cpolyfit.f90 +++ b/lib/fsk4hf/cpolyfit.f90 @@ -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 diff --git a/lib/fsk4hf/getfc1.f90 b/lib/fsk4hf/getfc1.f90 index aef97b774..b82b0303b 100644 --- a/lib/fsk4hf/getfc1.f90 +++ b/lib/fsk4hf/getfc1.f90 @@ -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 diff --git a/lib/fsk4hf/getfc2.f90 b/lib/fsk4hf/getfc2.f90 index a261da239..9506168dd 100644 --- a/lib/fsk4hf/getfc2.f90 +++ b/lib/fsk4hf/getfc2.f90 @@ -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 diff --git a/lib/fsk4hf/mskhfsim.f90 b/lib/fsk4hf/mskhfsim.f90 index 5b441417d..7383b0438 100644 --- a/lib/fsk4hf/mskhfsim.f90 +++ b/lib/fsk4hf/mskhfsim.f90 @@ -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 diff --git a/lib/fsk4hf/msksoftsym.f90 b/lib/fsk4hf/msksoftsym.f90 index ac844d017..016b23e1e 100644 --- a/lib/fsk4hf/msksoftsym.f90 +++ b/lib/fsk4hf/msksoftsym.f90 @@ -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