diff --git a/CMakeLists.txt b/CMakeLists.txt index 8a39490a0..b97a9516f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -332,6 +332,7 @@ set (wsjt_FSRCS lib/chkss2.f90 lib/coord.f90 lib/db.f90 + lib/fsk4hf/dbpsksim.f90 lib/decode4.f90 lib/decode65a.f90 lib/decode65b.f90 @@ -381,6 +382,7 @@ set (wsjt_FSRCS lib/gen4.f90 lib/gen65.f90 lib/gen9.f90 + lib/fsk4hf/genbpsk.f90 lib/geniscat.f90 lib/genmsk144.f90 lib/genmsk40.f90 @@ -474,6 +476,7 @@ set (wsjt_FSRCS lib/unpackmsg144.f90 lib/update_recent_calls.f90 lib/update_hasharray.f90 + lib/fsk4hf/watterson.f90 lib/wav11.f90 lib/wav12.f90 lib/xcor.f90 @@ -1096,6 +1099,9 @@ target_link_libraries (ldpcsim40 wsjt_fort wsjt_cxx) add_executable (ldpcsim120 lib/fsk4hf/ldpcsim120.f90 wsjtx.rc) target_link_libraries (ldpcsim120 wsjt_fort wsjt_cxx) +add_executable (dbpsksim lib/fsk4hf/dbpsksim.f90 wsjtx.rc) +target_link_libraries (dbpsksim wsjt_fort wsjt_cxx) + add_executable (ldpcsim144 lib/ldpcsim144.f90 wsjtx.rc) target_link_libraries (ldpcsim144 wsjt_fort wsjt_cxx) diff --git a/lib/fsk4hf/Makefile.win b/lib/fsk4hf/Makefile.win index e2687eac5..611843354 100644 --- a/lib/fsk4hf/Makefile.win +++ b/lib/fsk4hf/Makefile.win @@ -18,7 +18,7 @@ CFLAGS = -O2 -I. %.o: %.F90 ${FC} ${FFLAGS} -c $< -all: fsk4sim.exe +all: dbpsksim.exe OBJS0 = testpsk.o four2a.o bpfilter.o nonlinear.o tweak1.o spectrum.o smo.o testpsk: $(OBJS0) @@ -36,12 +36,29 @@ OBJS3 = fsk2sim.o four2a.o smo.o wavhdr.o gran.o fsk2sim: $(OBJS3) $(FC) -o fsk2sim $(OBJS3) C:\JTSDK\fftw3f\libfftw3f-3.dll -OBJS4 = fsk4sim.o four2a.o wavhdr.o gran.o tweak1.o genfsk4.o \ - dopspread.o smo.o smo121.o +OBJS4 = fsk4sim.o four2a.o gran.o genfsk4.o smo.o getsnr.o spec4.o \ + watterson.o db.o snr2_wsprlf.o pctile.o shell.o snr_wsprlf.o fsk4sim.exe: $(OBJS4) $(FC) -o fsk4sim.exe $(OBJS4) C:\JTSDK\fftw3f\libfftw3f-3.dll +OBJS5 = wsprlf.o four2a.o downsample.o +wsprlf.exe: $(OBJS5) + $(FC) -o wsprlf.exe $(OBJS5) C:\JTSDK\fftw3f\libfftw3f-3.dll + +OBJS6 = wspr_gmsk.o four2a.o gaussfilt.o +wspr_gmsk.exe: $(OBJS6) + $(FC) -o wspr_gmsk.exe $(OBJS6) C:\JTSDK\fftw3f\libfftw3f-3.dll + +OBJS7 = wspr_msk.o four2a.o bpfilter.o +wspr_msk.exe: $(OBJS7) + $(FC) -o wspr_msk.exe $(OBJS7) C:\JTSDK\fftw3f\libfftw3f-3.dll + +OBJS8 = dbpsksim.o four2a.o gran.o genbpsk.o watterson.o db.o \ + encode120.o bpdecode120.o platanh.o +dbpsksim.exe: $(OBJS8) + $(FC) -o dbpsksim.exe $(OBJS8) C:\JTSDK\fftw3f\libfftw3f-3.dll + .PHONY : clean clean: - $(RM) *.o testpsk.exe testfsk.exe fsk2sim.exe fsk4sim.exe + $(RM) *.o testpsk.exe testfsk.exe fsk2sim.exe fsk4sim.exe wsprlf.exe diff --git a/lib/fsk4hf/bpdecode120.f90 b/lib/fsk4hf/bpdecode120.f90 index de6450deb..c9b11ac39 100644 --- a/lib/fsk4hf/bpdecode120.f90 +++ b/lib/fsk4hf/bpdecode120.f90 @@ -1,7 +1,7 @@ -subroutine bpdecode120(llr,apmask,maxiterations,decoded,niterations) -! +subroutine bpdecode120(llr,apmask,maxiterations,decoded,niterations,cw) + ! A log-domain belief propagation decoder for the (120,60) code. -! + integer, parameter:: N=120, K=60, M=N-K integer*1 codeword(N),cw(N),apmask(N) integer colorder(N) diff --git a/lib/fsk4hf/dbpsksim.f90 b/lib/fsk4hf/dbpsksim.f90 new file mode 100644 index 000000000..640cc8c61 --- /dev/null +++ b/lib/fsk4hf/dbpsksim.f90 @@ -0,0 +1,243 @@ +program dbpsksim + + parameter (ND=121) !Data symbols: LDPC (120,60), r=1/2 + parameter (NN=ND) !Total symbols (121) + parameter (NSPS=28800) !Samples per symbol at 12000 sps + parameter (NZ=NSPS*NN) !Samples in waveform (3484800) + parameter (NFFT1=65536,NH1=NFFT1/2) + + character*8 arg + complex c(0:NZ-1) !Complex waveform + complex c2(0:NFFT1-1) !Short spectra + complex cr(0:NZ-1) + complex ct(0:NZ-1) + complex z0,z,zp + real s(-NH1+1:NH1) + real xnoise(0:NZ-1) !Generated random noise + real ynoise(0:NZ-1) !Generated random noise + real rxdata(120),llr(120) + integer id(NN) !Encoded NRZ data (values +/-1) + integer id1(NN) !Recovered data (1st pass) + integer id2(NN) !Recovered data (2nd pass) +! integer icw(NN) + integer*1 msgbits(60),decoded(60),codeword(120),apmask(120),cw(120) + data msgbits/0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,& + 0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,1,0,0,1,0,1,1,0,1,0/ + + nnn=0 + nargs=iargc() + if(nargs.ne.5) then + print*,'Usage: dbpsksim f0(Hz) delay(ms) fspread(Hz) iters snr(dB)' + print*,'Example: dbpsksim 1500 0 0 10 -35' + print*,'Set snr=0 to cycle through a range' + go to 999 + endif + call getarg(1,arg) + read(arg,*) f0 !Low tone frequency + call getarg(2,arg) + read(arg,*) delay + call getarg(3,arg) + read(arg,*) fspread + call getarg(4,arg) + read(arg,*) iters + call getarg(5,arg) + read(arg,*) snrdb + + twopi=8.d0*atan(1.d0) + fs=12000.d0 + dt=1.0/fs + ts=NSPS*dt + baud=1.d0/ts + txt=NZ*dt + bandwidth_ratio=2500.0/6000.0 + ndiff=1 !Encode/decode differentially + write(*,1000) baud,5*baud,txt,delay,fspread +1000 format('Baud:',f6.3,' BW:',f4.1,' TxT:',f6.1,' Delay:',f5.2, & + ' fSpread:',f5.2/) + + write(*,1004) +1004 format(' SNR e1 e2 ber1 ber2 fer1 fer2 fsigma'/55('-')) + + call encode120(msgbits,codeword) !Encode the test message + isna=-28 + isnb=-40 + if(snrdb.ne.0.0) then + isna=nint(snrdb) + isnb=isna + endif + do isnr=isna,isnb,-1 + snrdb=isnr + sig=sqrt(bandwidth_ratio) * 10.0**(0.05*snrdb) + if(snrdb.gt.90.0) sig=1.0 + nhard1=0 + nhard2=0 + nfe1=0 + nfe2=0 + sqf=0. + do iter=1,iters + nnn=nnn+1 + id(1)=1 !First bit is always 1 + id(2:NN)=2*codeword-1 + call genbpsk(id,f0,ndiff,0,c) !Generate the 4-FSK waveform + if(delay.ne.0.0 .or. fspread.ne.0.0) call watterson(c,delay,fspread) + c=sig*c !Scale to requested SNR + if(snrdb.lt.90) then + do i=0,NZ-1 !Generate gaussian noise + xnoise(i)=gran() + ynoise(i)=gran() + enddo + c=c + cmplx(xnoise,ynoise) !Add noise to signal + endif + +! First attempt at finding carrier frequency fc + nspec=NZ/NFFT1 + df1=12000.0/NFFT1 + s=0. + do k=1,nspec + ia=(k-1)*NSPS + ib=ia+NSPS-1 + c2(0:NSPS-1)=c(ia:ib) + c2(NSPS:)=0. + call four2a(c2,NFFT1,1,-1,1) + do i=0,NFFT1-1 + j=i + if(j.gt.NH1) j=j-NFFT1 + s(j)=s(j) + real(c2(i))**2 + aimag(c2(i))**2 + enddo + enddo + s=1.e-6*s + smax=0. + ipk=0 + ia=(1400.0)/df1 + ib=(1600.0)/df1 + do i=ia,ib + f=i*df1 + if(s(i).gt.smax) then + smax=s(i) + ipk=i + fc=f + endif + enddo + a=(s(ipk+1)-s(ipk-1))/2.0 + b=(s(ipk+1)+s(ipk-1)-2.0*s(ipk))/2.0 + dx=-a/(2.0*b) + fc=fc + df1*dx + sqf=sqf + (fc-f0)**2 + +! The following is for testing SNR calibration: +! sp5n=(s(ipk-2)+s(ipk-1)+s(ipk)+s(ipk+1)+s(ipk+2)) !Sig + 5*noise +! base=(sum(s)-sp5n)/(NFFT1-5.0) !Noise per bin +! psig=sp5n-5*base !Sig only +! pnoise=(2500.0/df1)*base !Noise in 2500 Hz +! xsnrdb=db(psig/pnoise) + + call genbpsk(id,fc,ndiff,1,cr) !Generate reference carrier + c=c*conjg(cr) !Mix to baseband + + z0=1.0 + ie0=1 + do j=1,NN !Demodulate + ia=(j-1)*NSPS + ib=ia+NSPS-1 + z=sum(c(ia:ib)) + zp=z*conjg(z0) + p=1.e-4*real(zp) + id1(j)=-1 + if(p.ge.0.0) id1(j)=1 + if(j.ge.2) rxdata(j-1)=p + z0=z + +! For testing, treat every 3rd symbol as having a known value (i.e., as Sync): +! ie=id(j)*ie0 +! if(mod(j,3).eq.0) write(12,1010) j,ie,1.e-3*ie*z, & +! atan2(aimag(ie*z),real(ie*z)) +!1010 format(2i4,3f10.3) +! ie0=ie + enddo + nhard1=nhard1 + count(id1.ne.id) !Count hard errors + + rxav=sum(rxdata)/120 + rx2av=sum(rxdata*rxdata)/120 + rxsig=sqrt(rx2av-rxav*rxav) + rxdata=rxdata/rxsig + ss=0.84 + llr=2.0*rxdata/(ss*ss) + apmask=0 + max_iterations=10 + call bpdecode120(llr,apmask,max_iterations,decoded,niterations,cw) + +! Count the hard errors in id1() and icw() +! icw(1)=1 +! icw(2:NN)=2*cw-1 +! nid1=0 +! ncw=0 +! ie0=1 +! do j=2,NN +! ib=(id(j)+1)/2 +! ib1=(id1(j)+1)/2 +! if(ib1.ne.ib) nid1=nid1+1 +! if(cw(j-1).ne.ib) ncw=ncw+1 +! enddo +! print*,niterations,nid1,ncw + +! Count frame errors + if(niterations.lt.0 .or. count(msgbits.ne.decoded).gt.0) nfe1=nfe1+1 + +! Generate a new reference carrier, using first-pass hard bits + call genbpsk(id1,0.0,ndiff,0,cr) + ct=c*conjg(cr) + call four2a(ct,NZ,1,-1,1) + df2=12000.0/NZ + pmax=0. + do i=0,NZ-1 + f=i*df2 + if(i.gt.NZ/2) f=(i-NZ)*df2 + if(abs(f).lt.1.0) then + p=real(ct(i))**2 + aimag(ct(i))**2 + if(p.gt.pmax) then + pmax=p + fc2=f + ipk=i + endif + endif + enddo + + call genbpsk(id,fc2,ndiff,1,cr) !Generate new ref carrier at fc2 + c=c*conjg(cr) + z0=1.0 + do j=1,NN !Demodulate + ia=(j-1)*NSPS + ib=ia+NSPS-1 + z=sum(c(ia:ib)) + zp=z*conjg(z0) + p=1.e-4*real(zp) + id2(j)=-1 + if(p.ge.0.0) id2(j)=1 + if(j.ge.2) rxdata(j-1)=p + z0=z + enddo + nhard2=nhard2 + count(id2.ne.id) !Count hard errors + + rxav=sum(rxdata)/120 + rx2av=sum(rxdata*rxdata)/120 + rxsig=sqrt(rx2av-rxav*rxav) + rxdata=rxdata/rxsig + ss=0.84 + llr=2.0*rxdata/(ss*ss) !Soft symbols + apmask=0 + max_iterations=10 + call bpdecode120(llr,apmask,max_iterations,decoded,niterations,cw) + if(niterations.lt.0 .or. count(msgbits.ne.decoded).gt.0) nfe2=nfe2+1 + enddo + + fsigma=sqrt(sqf/iters) + ber1=float(nhard1)/(NN*iters) + ber2=float(nhard2)/(NN*iters) + fer1=float(nfe1)/iters + fer2=float(nfe2)/iters + write(*,1050) snrdb,nhard1,nhard2,ber1,ber2,fer1,fer2,fsigma + write(14,1050) snrdb,nhard1,nhard2,ber1,ber2,fer1,fer2,fsigma +1050 format(f6.1,2i5,2f8.4,2f7.3,f8.2,3i5) + enddo + +999 end program dbpsksim diff --git a/lib/fsk4hf/genbpsk.f90 b/lib/fsk4hf/genbpsk.f90 new file mode 100644 index 000000000..45e945286 --- /dev/null +++ b/lib/fsk4hf/genbpsk.f90 @@ -0,0 +1,42 @@ +subroutine genbpsk(id,f00,ndiff,nref,c) + + parameter (ND=121) !Data symbols: LDPC (120,60), r=1/2 + parameter (NN=ND) !Total symbols (121) + parameter (NSPS=28800) !Samples per symbol at 12000 sps + parameter (NZ=NSPS*NN) !Samples in waveform (3456000) + + complex c(0:NZ-1) !Complex waveform + real*8 twopi,dt,fs,baud,f0,dphi,phi + integer id(NN) !Encoded NRZ data (values +/-1) + integer ie(NN) !Differentially encoded data + + f0=f00 + twopi=8.d0*atan(1.d0) + fs=12000.d0 + dt=1.0/fs + baud=1.d0/(NSPS*dt) + + ie(1)=1 !First bit is always 1 + do i=2,NN !Differentially encode + ie(i)=id(i)*ie(i-1) + enddo + +! Generate the BPSK waveform + phi=0.d0 + k=-1 + do j=1,NN + dphi=twopi*f0*dt + x=id(j) + if(ndiff.ne.0) x=ie(j) !Differential + if(nref.ne.0) x=1.0 !Generate reference carrier + do i=1,NSPS + k=k+1 + phi=phi+dphi + if(phi.gt.twopi) phi=phi-twopi + xphi=phi + c(k)=x*cmplx(cos(xphi),sin(xphi)) + enddo + enddo + + return +end subroutine genbpsk diff --git a/lib/fsk4hf/ldpcsim120.f90 b/lib/fsk4hf/ldpcsim120.f90 index a4dd4026b..cb8f79f41 100644 --- a/lib/fsk4hf/ldpcsim120.f90 +++ b/lib/fsk4hf/ldpcsim120.f90 @@ -12,6 +12,7 @@ integer*1, target:: i1Msg8BitBytes(9) integer*1, target:: i1Dec8BitBytes(9) integer*1 msgbits(60) integer*1 apmask(120) +integer*1 cw(120) integer*2 checksum integer colorder(120) integer nerrtot(120),nerrdec(120),nmpcbad(60) @@ -156,7 +157,7 @@ do idb = -10, 24 apmask=0 ! max_iterations is max number of belief propagation iterations - call bpdecode120(llr, apmask, max_iterations, decoded, niterations) + call bpdecode120(llr, apmask, max_iterations, decoded, niterations, cw) ! If the decoder finds a valid codeword, niterations will be .ge. 0. diff --git a/lib/fsk4hf/watterson.f90 b/lib/fsk4hf/watterson.f90 new file mode 100644 index 000000000..aa7db106e --- /dev/null +++ b/lib/fsk4hf/watterson.f90 @@ -0,0 +1,59 @@ +subroutine watterson(c,delay,fspread) + + parameter (NZ=3456000) + complex c(0:NZ-1) + complex c2(0:NZ-1) + complex cs1(0:NZ-1) + complex cs2(0:NZ-1) + + df=12000.0/NZ + if(fspread.gt.0.0) then + do i=0,NZ-1 + xx=gran() + yy=gran() + cs1(i)=cmplx(xx,yy) + xx=gran() + yy=gran() + cs2(i)=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 + f=i*df + if(i.gt.NZ/2) f=(i-NZ)*df + x=(f/fspread)**2 + a=0. + if(x.le.50.0) then + a=exp(-x) + endif + cs1(i)=a*cs1(i) + cs2(i)=a*cs2(i) +! if(abs(f).lt.10.0) then +! p1=real(cs1(i))**2 + aimag(cs1(i))**2 +! p2=real(cs2(i))**2 + aimag(cs2(i))**2 +! write(62,3101) f,db(p1+1.e-12)-60,db(p2+1.e-12)-60 +!3101 format(3f10.3) +! endif + enddo + call four2a(cs1,NZ,1,1,1) !Back to time domain + call four2a(cs2,NZ,1,1,1) + cs1=cs1/NZ + cs2=cs2/NZ + endif + + nshift=0.001*delay*12000.0 + c2=cshift(c,nshift) + + sq=0. + do i=0,NZ-1 + if(fspread.eq.0.0) c(i)=0.5*(c(i) + c2(i)) + if(fspread.gt.0.0) c(i)=0.5*(cs1(i)*c(i) + cs2(i)*c2(i)) + sq=sq + real(c(i))**2 + aimag(c(i))**2 +! write(61,3001) i/12000.0,c(i) +!3001 format(3f12.6) + enddo + rms=sqrt(sq/NZ) + c=c/rms + + return +end subroutine watterson