From 3124648fbcd615c4a89e1acfaf4a1f37a9df15d8 Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Wed, 19 Apr 2017 16:06:42 +0000 Subject: [PATCH] Add/update some experimental routines. git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7636 ab8295b8-cf94-4d9e-aec4-7959e3be5d79 --- CMakeLists.txt | 10 +- lib/fsk4hf/Makefile.win | 20 +++ lib/fsk4hf/dbpsksim.f90 | 120 ++++++------- lib/fsk4hf/fsk4sim.f90 | 375 +++++++++++++-------------------------- lib/fsk4hf/genbpsk.f90 | 10 +- lib/fsk4hf/genfsk4.f90 | 47 ++--- lib/fsk4hf/genmskhf.f90 | 129 ++++++++++++++ lib/fsk4hf/mskhfsim.f90 | 347 ++++++++++++++++++++++++++++++++++++ lib/fsk4hf/polyfit4.f90 | 109 ++++++++++++ lib/fsk4hf/twkfreq1.f90 | 26 +++ lib/fsk4hf/watterson.f90 | 34 ++-- lib/fsk4hf/wsprlf.f90 | 105 ++++++----- 12 files changed, 911 insertions(+), 421 deletions(-) create mode 100644 lib/fsk4hf/genmskhf.f90 create mode 100644 lib/fsk4hf/mskhfsim.f90 create mode 100644 lib/fsk4hf/polyfit4.f90 create mode 100644 lib/fsk4hf/twkfreq1.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 876cfb952..2cab0c6dc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -329,12 +329,12 @@ set (wsjt_FSRCS lib/ccf2.f90 lib/ccf65.f90 lib/fsk4hf/chkcrc10.f90 + lib/fsk4hf/chkcrc12.f90 lib/chkhist.f90 lib/chkmsg.f90 lib/chkss2.f90 lib/coord.f90 lib/db.f90 - lib/fsk4hf/dbpsksim.f90 lib/decode4.f90 lib/decode65a.f90 lib/decode65b.f90 @@ -385,10 +385,10 @@ 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 + lib/fsk4hf/genmskhf.f90 lib/genqra64.f90 lib/genwspr.f90 lib/geodist.f90 @@ -435,11 +435,13 @@ set (wsjt_FSRCS lib/mskrtd.f90 lib/msk144signalquality.f90 lib/msk144sim.f90 + lib/fsk4hf/mskhfsim.f90 lib/mskrtd.f90 lib/pctile.f90 lib/peakdt9.f90 lib/peakup.f90 lib/polyfit.f90 + lib/fsk4hf/polyfit4.f90 lib/prog_args.f90 lib/ps4.f90 lib/qra64a.f90 @@ -476,6 +478,7 @@ set (wsjt_FSRCS lib/timf2.f90 lib/tweak1.f90 lib/twkfreq.f90 + lib/fsk4hf/twkfreq1.f90 lib/twkfreq65.f90 lib/unpackmsg144.f90 lib/update_recent_calls.f90 @@ -1115,6 +1118,9 @@ target_link_libraries (ldpcsim144 wsjt_fort wsjt_cxx) add_executable (ldpcsim168 lib/fsk4hf/ldpcsim168.f90 wsjtx.rc) target_link_libraries (ldpcsim168 wsjt_fort wsjt_cxx) +add_executable (mskhfsim lib/fsk4hf/mskhfsim.f90 wsjtx.rc) +target_link_libraries (mskhfsim wsjt_fort wsjt_cxx) + add_executable (msk144sim lib/msk144sim.f90 wsjtx.rc) target_link_libraries (msk144sim wsjt_fort wsjt_cxx) diff --git a/lib/fsk4hf/Makefile.win b/lib/fsk4hf/Makefile.win index 611843354..ab383557d 100644 --- a/lib/fsk4hf/Makefile.win +++ b/lib/fsk4hf/Makefile.win @@ -58,6 +58,26 @@ OBJS8 = dbpsksim.o four2a.o gran.o genbpsk.o watterson.o db.o \ dbpsksim.exe: $(OBJS8) $(FC) -o dbpsksim.exe $(OBJS8) C:\JTSDK\fftw3f\libfftw3f-3.dll +OBJS9 = fsk4a.o four2a.o gran.o genfsk4a.o spec4.o \ + watterson.o db.o +fsk4a.exe: $(OBJS9) + $(FC) -o fsk4a.exe $(OBJS9) C:\JTSDK\fftw3f\libfftw3f-3.dll + +OBJS10 = gmsk8.o gaussfilt.o four2a.o +gmsk8.exe: $(OBJS10) + $(FC) -o gmsk8.exe $(OBJS10) C:\JTSDK\fftw3f\libfftw3f-3.dll + +OBJS11 = gmsksim.o four2a.o gran.o gengmsk.o genbpsk.o watterson.o db.o \ + encode168.o bpdecode168.o platanh.o gaussfilt.o tweak1.o smo121.o +gmsksim.exe: $(OBJS11) + $(FC) -o gmsksim.exe $(OBJS11) C:\JTSDK\fftw3f\libfftw3f-3.dll + +OBJS12 = mskhfsim.o four2a.o gran.o genmskhf.o watterson.o db.o \ + encode168.o bpdecode168.o platanh.o twkfreq1.o smo121.o \ + polyfit4.o +mskhfsim.exe: $(OBJS12) + $(FC) -o mskhfsim.exe $(OBJS12) C:\JTSDK\fftw3f\libfftw3f-3.dll + .PHONY : clean clean: diff --git a/lib/fsk4hf/dbpsksim.f90 b/lib/fsk4hf/dbpsksim.f90 index ce460f081..8b9ea12c6 100644 --- a/lib/fsk4hf/dbpsksim.f90 +++ b/lib/fsk4hf/dbpsksim.f90 @@ -5,14 +5,17 @@ program dbpsksim parameter (NSPS=28800) !Samples per symbol at 12000 sps parameter (NZ=NSPS*NN) !Samples in waveform (3484800) parameter (NFFT1=65536,NH1=NFFT1/2) + parameter (NFFT2=128,NH2=NFFT2/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 cz(0:NFFT2-1) complex z0,z,zp real s(-NH1+1:NH1) + real s2(-NH2+1:NH2) real xnoise(0:NZ-1) !Generated random noise real ynoise(0:NZ-1) !Generated random noise real rxdata(120),llr(120) @@ -26,8 +29,8 @@ program dbpsksim nnn=0 nargs=iargc() - if(nargs.ne.5) then - print*,'Usage: dbpsksim f0(Hz) delay(ms) fspread(Hz) iters snr(dB)' + if(nargs.ne.6) then + print*,'Usage: dbpsksim f0(Hz) delay(ms) fspread(Hz) ndiff iters snr(dB)' print*,'Example: dbpsksim 1500 0 0 10 -35' print*,'Set snr=0 to cycle through a range' go to 999 @@ -39,8 +42,10 @@ program dbpsksim call getarg(3,arg) read(arg,*) fspread call getarg(4,arg) - read(arg,*) iters + read(arg,*) ndiff call getarg(5,arg) + read(arg,*) iters + call getarg(6,arg) read(arg,*) snrdb twopi=8.d0*atan(1.d0) @@ -50,13 +55,12 @@ program dbpsksim 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 + write(*,1000) baud,5*baud,txt,delay,fspread,ndiff 1000 format('Baud:',f6.3,' BW:',f4.1,' TxT:',f6.1,' Delay:',f5.2, & - ' fSpread:',f5.2/) + ' fSpread:',f5.2,' ndiff:',i2/) write(*,1004) -1004 format(' SNR e1 e2 ber1 ber2 fer1 fer2 fsigma'/55('-')) +1004 format(' SNR err ber fer fsigma'/35('-')) call encode120(msgbits,codeword) !Encode the test message isna=-28 @@ -69,8 +73,8 @@ program dbpsksim snrdb=isnr sig=sqrt(bandwidth_ratio) * 10.0**(0.05*snrdb) if(snrdb.gt.90.0) sig=1.0 - nhard1=0 - nhard2=0 + nhard=0 + nhardc=0 nfe1=0 nfe2=0 sqf=0. @@ -89,7 +93,7 @@ program dbpsksim c=c + cmplx(xnoise,ynoise) !Add noise to signal endif -! First attempt at finding carrier frequency fc +! First attempt at finding carrier frequency fc: 64k FFTs ==> avg power spectra nspec=NZ/NFFT1 df1=12000.0/NFFT1 s=0. @@ -121,7 +125,7 @@ program dbpsksim 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 + fc=fc + df1*dx !Estimated carrier frequency sqf=sqf + (fc-f0)**2 ! The following is for testing SNR calibration: @@ -132,29 +136,21 @@ program dbpsksim ! xsnrdb=db(psig/pnoise) call genbpsk(id,fc,ndiff,1,cr) !Generate reference carrier - c=c*conjg(cr) !Mix to baseband + c=c*conjg(cr) !Mix signal 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)) + cz(j-1)=z 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 @@ -166,42 +162,30 @@ program dbpsksim 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 +! Find carrier frequency from squared cz array. + cz(121:)=0. + cz=cz*cz + call four2a(cz,NFFT2,1,-1,1) + s2max=0. + do i=0,NFFT2-1 + j=i + if(i.gt.NH2) j=j-NFFT2 + s2(j)=real(cz(i))**2 + aimag(cz(i))**2 + if(s2(j).gt.s2max) then + s2max=s2(j) + jpk=j + endif +! write(16,1200) j*baud/NFFT2,1.e-12*s2(j) +!1200 format(2f12.3) + enddo + a=(s2(jpk+1)-s2(jpk-1))/2.0 + b=(s2(jpk+1)+s2(jpk-1)-2.0*s2(jpk))/2.0 + dx=-a/(2.0*b) + fc2=0.5*(jpk+dx)*baud/NFFT2 + call genbpsk(id,fc2,ndiff,1,cr) !Generate new ref carrier at fc2 c=c*conjg(cr) z0=1.0 @@ -209,14 +193,21 @@ program dbpsksim ia=(j-1)*NSPS ib=ia+NSPS-1 z=sum(c(ia:ib)) + if(j.eq.1) z0=z 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 + ierr=0 + if(id2(j).ne.id(j)) ierr=1 + id3=-1 + if(real(z).ge.0.0) id3=1 + if(j.ge.2 .and. id3.ne.id(j)) nhardc=nhardc+1 + if(j.ge.2 .and. ndiff.eq.0) rxdata(j-1)=real(z) z0=z enddo - nhard2=nhard2 + count(id2.ne.id) !Count hard errors + nhard=nhard + count(id2.ne.id) !Count hard errors rxav=sum(rxdata)/120 rx2av=sum(rxdata*rxdata)/120 @@ -226,20 +217,25 @@ program dbpsksim llr=2.0*rxdata/(ss*ss) !Soft symbols apmask=0 max_iterations=10 + decoded=0 call bpdecode120(llr,apmask,max_iterations,decoded,niterations,cw) +! if(niterations.lt.0) then +! llr=-llr +! call bpdecode120(llr,apmask,max_iterations,decoded,niterations,cw) +! if(niterations.ge.0) nhard=NN*iters-nhard +! endif if(niterations.ge.0) call chkcrc10(decoded,nbadcrc) if(niterations.lt.0 .or. count(msgbits.ne.decoded).gt.0 .or. & nbadcrc.ne.0) nfe2=nfe2+1 enddo + if(ndiff.eq.0) nhard=nhardc 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) + ber=float(nhard)/(NN*iters) + fer=float(nfe2)/iters + write(*,1050) snrdb,nhard,ber,fer,fsigma + write(14,1050) snrdb,nhard,ber,fer,fsigma +1050 format(f6.1,i5,f8.4,f7.3,f8.2) enddo 999 end program dbpsksim diff --git a/lib/fsk4hf/fsk4sim.f90 b/lib/fsk4hf/fsk4sim.f90 index 09f1c7f30..9800616e9 100644 --- a/lib/fsk4hf/fsk4sim.f90 +++ b/lib/fsk4hf/fsk4sim.f90 @@ -1,100 +1,62 @@ program fsk4sim - use wavhdr - parameter (NR=4) !Ramp up, ramp down - parameter (NS=12) !Sync symbols (2 @ Costas 4x4) - parameter (ND=84) !Data symbols: LDPC (168,84), r=1/2 - parameter (NN=NR+NS+ND) !Total symbols (100) - parameter (NSPS=2688) !Samples per symbol at 12000 sps - parameter (NZ=NSPS*NN) !Samples in waveform (268800) - parameter (NSYNC=NS*NSPS) !Samples in sync waveform (32256) - parameter (NFFT=512*1024) - parameter (NDOWN=84) !Downsample factor - parameter (NFFT2=NZ/NDOWN,NH2=NFFT2/2) !3200 - parameter (NSPSD=NFFT2/NN) !Samples per symbol after downsample + parameter (ND=60) !Data symbols: LDPC (120,60), r=1/2 + parameter (NN=ND) !Total symbols (60) + parameter (NSPS=57600) !Samples per symbol at 12000 sps + parameter (NZ=NSPS*NN) !Samples in waveform (3456000) - type(hdr) header !Header for .wav file character*8 arg - complex c(0:NFFT-1) !Complex waveform - complex cf(0:NFFT-1) - complex cs(0:NSYNC-1) - complex ct(0:NSPS-1) - complex csync(0:NSYNC-1) - complex c2(0:NFFT2-1) - complex c2a(0:NSPSD-1) - complex cf2(0:NFFT2-1) - complex cx(0:3,NN) - complex z,zpk - logical snrtest - real*8 twopi,dt,fs,baud,f0,dphi,phi + complex c(0:NZ-1) !Complex waveform + complex cr(0:NZ-1) + complex cs(NSPS,NN) + complex cps(0:3) + complex ct(0:2*NN-1) + complex z,w,zsum + real r(0:NZ-1) + real s(NSPS,NN) + real savg(NSPS) real tmp(NN) !For generating random data real xnoise(0:NZ-1) !Generated random noise - real s(NSYNC/2) real ps(0:3) -! integer*2 iwave(NZ) !Generated waveform - integer id(NN) !Encoded 2-bit data (values 0-3) - integer id2(NN) !Decoded after downsampling - integer icos4(4) !4x4 Costas array - data icos4/0,1,3,2/,eps/1.e-8/ + integer id(NN) !Encoded 2-bit data (values 0-3) + integer id2(NN) !Recovered data + equivalence (r,cr) + nnn=0 nargs=iargc() - if(nargs.ne.4) then - print*,'Usage: fsk8sim f0 fspread iters snr' + if(nargs.ne.6) then + print*,'Usage: fsk8sim f0 delay(ms) fspread(Hz) nts iters snr(dB)' go to 999 endif call getarg(1,arg) read(arg,*) f0 !Low tone frequency call getarg(2,arg) - read(arg,*) fspread + read(arg,*) delay call getarg(3,arg) - read(arg,*) iters + read(arg,*) fspread call getarg(4,arg) + read(arg,*) nts + call getarg(5,arg) + read(arg,*) iters + call getarg(6,arg) read(arg,*) snrdb - snrtest=.false. - if(iters.lt.0) then - snrtest=.true. - iters=abs(iters) - endif - twopi=8.d0*atan(1.d0) fs=12000.d0 dt=1.0/fs ts=NSPS*dt baud=1.d0/ts txt=NZ*dt - -! Generate sync waveform - phi=0.d0 - k=-1 - do j=1,12 - n=mod(j-1,4) + 1 - dphi=twopi*(icos4(n)*baud)*dt - do i=1,NSPS - k=k+1 - phi=phi+dphi - if(phi.gt.twopi) phi=phi-twopi - xphi=phi - csync(k)=cmplx(cos(xphi),-sin(xphi)) - enddo - enddo bandwidth_ratio=2500.0/6000.0 - header=default_header(12000,NZ) + write(*,1000) baud,5*baud,txt,delay,fspread,nts +1000 format('Baud:',f6.3,' BW:',f5.1,' TxT:',f5.1,' Delay:',f5.2, & + ' fSpread:',f5.2,' nts:',i3/) - write(*,1000) 2*ND,ND,NS,NN,NSPS,baud,txt,fspread -1000 format('LDPC('i3,',',i2,') SyncSym:',i2,' ChanSym:',i3,' NSPS:',i4, & - ' Baud:',f6.3,' TxT:',f5.1,' fDop:',f5.2/) - if(snrtest) then - write(*,1002) -1002 format(5x,'SNR test'/'Requested Measured Difference') - else - write(*,1004) -1004 format(' SNR Sync Sym1 Sym2 Bits SyncErr Sym1Err BER'/ & - 60('-')) - endif + write(*,1004) +1004 format(' SNR Sym Bit SER BER Sym Bit SER BER'/59('-')) - isna=-15 - isnb=-27 + isna=-25 + isnb=-40 if(snrdb.ne.0.0) then isna=nint(snrdb) isnb=isna @@ -103,204 +65,121 @@ program fsk4sim snrdb=isnr sig=sqrt(2*bandwidth_ratio) * 10.0**(0.05*snrdb) if(snrdb.gt.90.0) sig=1.0 -! open(10,file='000000_0001.wav',access='stream',status='unknown') - - nsyncerr=0 - nharderr=0 - nherr=0 - nbiterr=0 + nhard1=0 + nhard2=0 + nbit1=0 + nbit2=0 + nh2=0 + nb2=0 do iter=1,iters + nnn=nnn+1 id=0 - if(.not.snrtest) then - ! Generate random data - call random_number(tmp) - where(tmp.ge.0.25 .and. tmp.lt.0.50) id=1 - where(tmp.ge.0.50 .and. tmp.lt.0.75) id=2 - where(tmp.ge.0.75) id=3 - id(1:2)=icos4(3:4) !Ramp up - id(45:48)=icos4 !Costas sync - id(49:52)=icos4 !Costas sync - id(53:56)=icos4 !Costas sync - id(NN-1:NN)=icos4(1:2) !Ramp down - endif + call random_number(tmp) + where(tmp.ge.0.25 .and. tmp.lt.0.50) id=1 + where(tmp.ge.0.50 .and. tmp.lt.0.75) id=2 + where(tmp.ge.0.75) id=3 - call genfsk4(id,f0,c) !Generate the 4-FSK waveform + call genfsk4(id,f0,nts,c) !Generate the 4-FSK waveform + call watterson(c,delay,fspread) if(sig.ne.1.0) c=sig*c !Scale to requested SNR - if(snrdb.lt.90) then do i=0,NZ-1 !Generate gaussian noise xnoise(i)=gran() enddo endif - if(fspread.gt.0.0) call dopspread(c,fspread) - c(0:NZ-1)=real(c(0:NZ-1)) + xnoise !Add noise to signal - -! fac=32767.0 -! rms=100.0 -! if(snrdb.ge.90.0) iwave(1:NZ)=nint(fac*aimag(c(0:NZ-1))) -! if(snrdb.lt.90.0) iwave(1:NZ)=nint(rms*aimag(c(0:NZ-1))) -! call set_wsjtx_wav_params(14.0,'JT65 ',1,30,iwave) -! write(10) header,iwave !Save the .wav file - - ppmax=0. - fpk=-99. - xdt=-99. - df1=12000.0/NSYNC - iaa=nint(250.0/df1) - ibb=nint(2750.0/df1) - if(.not.snrtest) then - do j4=-40,40 - ia=(44+0.25*j4)*NSPS - ib=ia+NSYNC-1 - cs=csync*c(ia:ib) - call four2a(cs,NSYNC,1,-1,1) !Transform to frequency domain - s=0. - do i=iaa,ibb - s(i)=1.e-6*(real(cs(i))**2 + aimag(cs(i))**2) - enddo - - if(j4.eq.0) then - do i=iaa,ibb - write(66,3301) i*df1,s(i) -3301 format(f10.3,2f12.6) - enddo - endif - - call smo121(s,NSYNC/2) - - if(j4.eq.0) then - do i=iaa,ibb - write(67,3301) i*df1,s(i) - enddo - endif - - do i=iaa,ibb - if(s(i).gt.ppmax) then - fpk=i*df1 - xdt=0.25*j4*ts - ppmax=s(i) - endif - enddo - - enddo - endif - if(xdt.ne.0.0 .or. fpk.ne.1500.0) nsyncerr=nsyncerr+1 - -! Compute spectrum again - cf=c - df2=12000.0/NZ - call four2a(cf,NZ,1,-1,1) !Transform to frequency domain - - if(snrtest) then - width=5.0*df2 + fspread - iz=nint(2500.0/df2) + 2 - if(iter.eq.1) then - pnoise=0. - psig=0. - n=0 - endif - do i=0,iz !Remove spectral sidelobes - f=i*df2 - if(i.gt.NZ/2) f=(i-NZ)*df2 - p=1.e-6*(real(cf(i))**2 + aimag(cf(i))**2) - if(abs(f-f0).lt.width) then - psig=psig+p - n=n+1 - else - pnoise=pnoise + p - endif - enddo - if(iter.eq.iters) then - db=10.0*log10(psig/pnoise) - write(*,1010) snrdb,db,db-snrdb -1010 format(f7.1,2f9.1) - endif - go to 40 - endif - -! Select a small frequency slice around fpk. - cf=cf/NZ - ib=nint(fpk/df2)+NH2 - ia=ib-NFFT2+1 - cf2=cshift(cf(ia:ib),NH2-1) - flo=-baud - fhi=4*baud - do i=0,NFFT2-1 - f=i*df2 - if(i.gt.NH2) f=(i-NFFT2)*df2 - if(f.le.flo .or. f.ge.fhi) cf2(i)=0. - s2=real(cf2(i))**2 + aimag(cf2(i))**2 - write(15,3001) f,s2,10*log10(s2+eps) -3001 format(f10.3,2f15.6) - enddo + r(0:NZ-1)=real(c(0:NZ-1)) + xnoise !Add noise to signal - c2=cf2 - call four2a(c2,NFFT2,1,1,1) !Back to time domain + call snr2_wsprlf(r,freq,snr2500,width,1) + write(*,3001) freq,snr2500,width +3001 format(40x,3f10.3) + + df=12000.0/(2*NSPS) +! i0=nint(f0/df) +! i0=nint((1500.0+freq)/df) + i0=nint((f0+freq)/df) + call spec4(r,cs,s,savg) - fshift=NSPS*baud/NSPSD - dt2=dt*NDOWN do j=1,NN - ia=(j-1)*NSPSD - ib=ia+NSPSD-1 - c2a=c2(ia:ib) - call four2a(c2a,NSPSD,1,-1,1) !To freq domain - cx(0:3,j)=c2a(0:3) - ipk=-1 - zpk=0. - pmax=0.0 - do i=0,3 - if(abs(cx(i,j)).gt.pmax) then - ipk=i - zpk=cx(i,j) - pmax=abs(zpk) - endif - enddo - id2(j)=ipk - if(ipk.ne.id(j)) nherr=nherr+1 - write(16,3003) j,id(j),ipk,ipk-id(j),abs(zpk), & - atan2(aimag(zpk),real(zpk)),abs(cx(0:3,j)) -3003 format(3i3,i4,6f9.3) - enddo - - ipk=0 - do j=1,NN - ia=(j-1)*NSPS + 1 - ib=ia+NSPS - pmax=0. - do i=0,3 - f=fpk + i*baud - call tweak1(c(ia:ib),NSPS,-f,ct) - z=sum(ct) - ps(i)=1.e-3*(real(z)**2 + aimag(z)**2) - if(ps(i).gt.pmax) then - ipk=i - pmax=ps(i) - endif - enddo - nlo=0 nhi=0 + ps=s(i0:i0+6*nts:2*nts,j) + cps=cs(i0:i0+6*nts:2*nts,j) if(max(ps(1),ps(3)).ge.max(ps(0),ps(2))) nlo=1 if(max(ps(2),ps(3)).ge.max(ps(0),ps(1))) nhi=1 - if(nlo.ne.iand(id(j),1)) nbiterr=nbiterr+1 - if(nhi.ne.iand(id(j)/2,1)) nbiterr=nbiterr+1 - if(ipk.ne.id(j)) nharderr=nharderr+1 - write(17,1040) j,ps,ipk,id(j),id2(j),2*nhi+nlo,nhi,nlo,nbiterr -1040 format(i3,4f12.1,7i4) + id2(j)=2*nhi+nlo + z=cps(id2(j)) + ct(j-1)=z enddo + nh1=count(id.ne.id2) + nb1=count(iand(id,1).ne.iand(id2,1)) + count(iand(id,2).ne.iand(id2,2)) + + ct(NN:)=0. + call four2a(ct,2*NN,1,-1,1) + df2=baud/(2*NN) + ct=cshift(ct,NN) + ppmax=0. + dfpk=0. + do i=0,2*NN-1 + f=(i-NN)*df2 + pp=real(ct(i))**2 + aimag(ct(i))**2 + if(pp.gt.ppmax) then + ppmax=pp + dfpk=f + endif + enddo + + zsum=0. + do j=1,NN + phi=(j-1)*twopi*dfpk*ts + w=cmplx(cos(phi),sin(phi)) + cps=cs(i0:i0+6*nts:2*nts,j)*conjg(w) + z=cps(id2(j)) + ct(j)=z + zsum=zsum+z + write(12,1042) j,id(j),id2(j),20*ps,atan2(aimag(z),real(z)), & + atan2(aimag(zsum),real(zsum)),zsum +1042 format(3i2,6f8.3,2f8.1) + enddo + + phi0=atan2(aimag(zsum),real(zsum)) + zsum=0. + do j=1,NN + phi=(j-1)*twopi*dfpk*ts + phi0 + w=cmplx(cos(phi),sin(phi)) + nlo=0 + nhi=0 + cps=cs(i0:i0+6*nts:2*nts,j)*conjg(w) + ps=real(cps) + if(max(ps(1),ps(3)).ge.max(ps(0),ps(2))) nlo=1 + if(max(ps(2),ps(3)).ge.max(ps(0),ps(1))) nhi=1 + id2(j)=2*nhi+nlo + z=cps(id2(j)) + ct(j)=z + zsum=zsum+z + enddo + + nh2=count(id.ne.id2) + nb2=count(iand(id,1).ne.iand(id2,1)) + count(iand(id,2).ne.iand(id2,2)) + nhard1=nhard1+nh1 + nhard2=nhard2+nh2 + nbit1=nbit1+nb1 + nbit2=nbit2+nb2 + + fdiff=1500.0+freq - f0 + write(13,1040) snrdb,snr2500,f0,fdiff,width,nh1,nb1,nh2,nb2 +1040 format(2f7.1,f9.2,f7.2,f6.1,2(i8,i6)) 40 continue enddo - if(.not.snrtest) then - fsyncerr=float(nsyncerr)/iters - ser=float(nharderr)/(NN*iters) - ber=float(nbiterr)/(2*NN*iters) - write(*,1050) snrdb,nsyncerr,nharderr,nherr,nbiterr,fsyncerr,ser,ber - write(18,1050) snrdb,nsyncerr,nharderr,nherr,nbiterr,fsyncerr,ser,ber -1050 format(f6.1,4i6,3f10.6) - endif + ser1=float(nhard1)/(NN*iters) + ser2=float(nhard2)/(NN*iters) + ber1=float(nbit1)/(2*NN*iters) + ber2=float(nbit2)/(2*NN*iters) + write(*,1050) snrdb,nhard1,nbit1,ser1,ber1,nhard2,nbit2,ser2,ber2 + write(14,1050) snrdb,nhard1,nbit1,ser1,ber1,nhard2,nbit2,ser2,ber2 +1050 format(f6.1,2(2i5,2f8.4)) enddo - if(.not.snrtest) write(*,1060) iters,100*iters,100*iters,200*iters -1060 format(60('-')/'Max: ',4i6) + write(*,1060) NN*iters,2*NN*iters +1060 format(59('-')/'Max: ',2i5) 999 end program fsk4sim diff --git a/lib/fsk4hf/genbpsk.f90 b/lib/fsk4hf/genbpsk.f90 index 45e945286..6c283e72f 100644 --- a/lib/fsk4hf/genbpsk.f90 +++ b/lib/fsk4hf/genbpsk.f90 @@ -16,10 +16,12 @@ subroutine genbpsk(id,f00,ndiff,nref,c) 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 + if(ndiff.ne.0) then + ie(1)=1 !First bit is always 1 + do i=2,NN !Differentially encode + ie(i)=id(i)*ie(i-1) + enddo + endif ! Generate the BPSK waveform phi=0.d0 diff --git a/lib/fsk4hf/genfsk4.f90 b/lib/fsk4hf/genfsk4.f90 index 4d6805cc1..d8b0394b2 100644 --- a/lib/fsk4hf/genfsk4.f90 +++ b/lib/fsk4hf/genfsk4.f90 @@ -1,24 +1,16 @@ -subroutine genfsk4(id,f0,c) +subroutine genfsk4(id,f00,nts,c) - parameter (NR=4) !Ramp up, ramp down - parameter (NS=12) !Sync symbols (2 @ Costas 4x4) - parameter (ND=84) !Data symbols: LDPC (168,84), r=1/2 - parameter (NN=NR+NS+ND) !Total symbols (100) - parameter (NSPS=2688) !Samples per symbol at 12000 sps - parameter (NZ=NSPS*NN) !Samples in waveform (268800) - parameter (NFFT=512*1024) - parameter (NSYNC=NS*NSPS) - parameter (NDOWN=168) - parameter (NFFT2=NZ/NDOWN,NH2=NFFT2/2) !3200 - parameter (NSPSD=NFFT2/NN) + parameter (ND=60) !Data symbols: LDPC (120,60), r=1/2 + parameter (NN=ND) !Total symbols (60) + parameter (NSPS=57600) !Samples per symbol at 12000 sps + parameter (NZ=NSPS*NN) !Samples in waveform (3456000) + parameter (NFFT=NZ) !Full length FFT - complex c(0:NFFT-1) !Complex waveform - complex cf(0:NFFT-1) + complex c(0:NFFT-1) !Complex waveform real*8 twopi,dt,fs,baud,f0,dphi,phi - integer id(NN) !Encoded 2-bit data (values 0-3) - integer icos4(4) !4x4 Costas array - data icos4/0,1,3,2/ + integer id(NN) !Encoded 2-bit data (values 0-3) + f0=f00 twopi=8.d0*atan(1.d0) fs=12000.d0 dt=1.0/fs @@ -30,7 +22,7 @@ subroutine genfsk4(id,f0,c) phi=0.d0 k=-1 do j=1,NN - dphi=twopi*(f0 + id(j)*baud)*dt + dphi=twopi*(f0 + nts*id(j)*baud)*dt do i=1,NSPS k=k+1 phi=phi+dphi @@ -40,24 +32,5 @@ subroutine genfsk4(id,f0,c) enddo enddo - nh=NFFT/2 - df=12000.0/NFFT - cf=c - call four2a(cf,NFFT,1,-1,1) !Transform to frequency domain - - if(sum(id).ne.0) then - flo=f0-baud - fhi=f0+4*baud - do i=0,NFFT-1 !Remove spectral sidelobes - f=i*df - if(i.gt.nh) f=(i-nfft)*df - if(f.le.flo .or. f.ge.fhi) cf(i)=0. - enddo - endif - - c=cf - call four2a(c,NFFT,1,1,1) !Transform back to time domain - c=c/nfft - return end subroutine genfsk4 diff --git a/lib/fsk4hf/genmskhf.f90 b/lib/fsk4hf/genmskhf.f90 new file mode 100644 index 000000000..a9f2665c6 --- /dev/null +++ b/lib/fsk4hf/genmskhf.f90 @@ -0,0 +1,129 @@ +subroutine genmskhf(msgbits,id,icw,cbb,csync) + +!Encode an MSK-HF message, produce baseband waveform and sync vector. + + parameter (KK=84) !Information bits (72 + CRC12) + parameter (ND=168) !Data symbols: LDPC (168,84), r=1/2 + parameter (NS=65) !Sync symbols (2 x 26 + Barker 13) + parameter (NR=3) !Ramp up/down + parameter (NN=NR+NS+ND) !Total symbols (236) + parameter (NSPS=16) !Samples per MSK symbol (16) + parameter (N2=2*NSPS) !Samples per OQPSK symbol (32) + parameter (NZ=NSPS*NN) !Samples in baseband waveform (3760) + + complex cbb(0:NZ-1) + complex csync(0:NZ-1) + real x(0:NZ-1) + real y(0:NZ-1) + real pp(N2) + logical first + integer*1 msgbits(KK),codeword(ND) + integer icw(ND) + integer id(NS+ND) + integer isync(26) !Long sync vector + integer ib13(13) !Barker 13 code + data ib13/1,1,1,1,1,-1,-1,1,1,-1,1,-1,1/ + data first/.true./ + save first,isync,twopi,pp + + if(first) then + n=z'2c1aeb1' + do i=1,26 + isync(i)=-1 + if(iand(n,1).eq.1) isync(i)=1 + n=n/2 + enddo + + twopi=8.0*atan(1.0) + do i=1,N2 !Half-sine shaped pulse + pp(i)=sin(0.5*(i-1)*twopi/N2) + enddo + first=.false. + endif + + call random_number(x) + codeword=0 + where(x(1:ND).ge.0.5) codeword=1 + call encode168(msgbits,codeword) !Encode the test message + icw=2*codeword - 1 + +! Message structure: R1 26*(S1+D1) S13 26*(D1+S1) R1 +! Generate QPSK without any offset; then shift the y array to get OQPSK. + +! Do the I channel first: results in array x + n=0 + k=0 + ia=0 + ib=NSPS-1 + x(ia:ib)=0. !Ramp up (half-symbol; shape TBD) + do j=1,26 !Insert group of 26*(S1+D1) + ia=ib+1 + ib=ia+N2-1 + n=n+1 + id(n)=2*isync(j) + x(ia:ib)=isync(j)*pp !Insert Sync bit + ia=ib+1 + ib=ia+N2-1 + k=k+1 + n=n+1 + id(n)=icw(k) + x(ia:ib)=id(n)*pp !Insert data bit + enddo + + do j=1,13 !Insert Barker 13 code + ia=ib+1 + ib=ia+N2-1 + n=n+1 + id(n)=2*ib13(j) + x(ia:ib)=ib13(j)*pp + enddo + + do j=1,26 !Insert group of 26*(S1+D1) + ia=ib+1 + ib=ia+N2-1 + k=k+1 + n=n+1 + id(n)=icw(k) + x(ia:ib)=id(n)*pp !Insert data bit + ia=ib+1 + ib=ia+N2-1 + n=n+1 + id(n)=2*isync(j) + x(ia:ib)=isync(j)*pp !Insert Sync bit + enddo + ia=ib+1 + ib=ia+NSPS-1 + x(ia:ib)=0. !Ramp down (half-symbol; shape TBD) + +! Now do the Q channel: results in array y + ia=0 + ib=NSPS-1 + y(ia:ib)=0. !Ramp up (half-symbol; shape TBD) + do j=1,116 + ia=ib+1 + ib=ia+N2-1 + k=k+1 + n=n+1 + id(n)=icw(k) + y(ia:ib)=id(n)*pp + enddo + ia=ib+1 + ib=ia+NSPS-1 + y(ia:ib)=0. !Ramp down (half-symbol; shape TBD) + y=cshift(y,-NSPS) !Shift Q array to get OQPSK + cbb=cmplx(x,y) !Complex baseband waveform + + ib=NSPS-1 + ib2=NSPS-1+64*N2 + do j=1,26 !Zero all data symbols in x + ia=ib+1+N2 + ib=ia+N2-1 + x(ia:ib)=0. + ia2=ib2+1+N2 + ib2=ia2+N2-1 + x(ia2:ib2)=0. + enddo + csync=x + + return +end subroutine genmskhf diff --git a/lib/fsk4hf/mskhfsim.f90 b/lib/fsk4hf/mskhfsim.f90 new file mode 100644 index 000000000..915ef4983 --- /dev/null +++ b/lib/fsk4hf/mskhfsim.f90 @@ -0,0 +1,347 @@ +program msksim + +! Simulate characteristics of a potential "MSK10" mode using LDPC (168,84) +! code, OQPDK modulation, and 30 s 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 + + parameter (KK=84) !Information bits (72 + CRC12) + parameter (ND=168) !Data symbols: LDPC (168,84), r=1/2 + parameter (NS=65) !Sync symbols (2 x 26 + Barker 13) + parameter (NR=3) !Ramp up/down + parameter (NN=NR+NS+ND) !Total symbols (236) + 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 (3760) + parameter (NFFT1=4*NSPS,NH1=NFFT1/2) + + 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 cs(0:NZ-1) !For computing spectrum + complex c2(0:NFFT1-1) !Short spectra + complex zz(NS+ND) !Complex symbol values (intermediate) + complex z,z0 + real s(-NH1+1:NH1) !Coarse spectrum + real xnoise(0:NZ-1) !Generated random noise + real ynoise(0:NZ-1) !Generated random noise + real x(NS),yi(NS),yq(NS) !For complex polyfit + 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 icw(NN) + integer*1 msgbits(KK),decoded(KK),apmask(ND),cw(ND) +! 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,0,1,1,1,0,1,1,0,1,1, & + 1,1,0,1,0,1,1,0,0,0,0,0,1,0,0,0,0,0,1,0,1,0,1,0/ + + nargs=iargc() + if(nargs.ne.6) then + print*,'Usage: msksim f0(Hz) delay(ms) fspread(Hz) maxn iters snr(dB)' + print*,'Example: msksim 20 0 0 5 10 -20' + print*,'Set snr=0 to cycle through a range' + go to 999 + endif + call getarg(1,arg) + read(arg,*) f0 !Generated carrier frequency + call getarg(2,arg) + read(arg,*) delay !Delta_t (ms) for Watterson model + call getarg(3,arg) + read(arg,*) fspread !Fspread (Hz) for Watterson model + call getarg(4,arg) + read(arg,*) maxn !Max nterms for polyfit + call getarg(5,arg) + read(arg,*) iters !Iterations at each SNR + call getarg(6,arg) + read(arg,*) snrdb !Specified SNR_2500 + + twopi=8.0*atan(1.0) + fs=12000.0/72.0 !Sample rate = 166.6666667 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,1.5*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(*,1004) +1004 format(/' SNR err ber fer fsigma'/37('-')) + + do i=1,N2 !Half-sine pulse shape + pp(i)=sin(0.5*(i-1)*twopi/(2*NSPS)) + enddo + + 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 + + isna=-10 + isnb=-30 + if(snrdb.ne.0.0) then + isna=nint(snrdb) + isnb=isna + endif + do isnr=isna,isnb,-1 !Loop over SNR range + snrdb=isnr + sig=sqrt(bandwidth_ratio) * 10.0**(0.05*snrdb) + if(snrdb.gt.90.0) sig=1.0 + nhard=0 + nhardsync=0 + 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 + 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 AWGN noise + endif + +!----------------------------------------------------------------- fc1 +! First attempt at finding carrier frequency, fc1: low-resolution power spectra + nspec=NZ/NFFT1 + df1=fs/NFFT1 + s=0. + do k=1,nspec + ia=(k-1)*N2 + ib=ia+N2-1 + c2(0:N2-1)=c(ia:ib) + c2(N2:)=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 +! call smo121(s,NFFT1) + smax=0. + ipk=0 + fc1=0. + ia=nint(40.0/df1) + do i=-ia,ia + f=i*df1 + if(s(i).gt.smax) then + smax=s(i) + ipk=i + fc1=f + endif +! write(51,3001) f,s(i),db(s(i)) +! 3001 format(f10.3,e12.3,f10.3) + enddo + +! The following is for testing SNR calibration: +! sp3n=(s(ipk-1)+s(ipk)+s(ipk+1)) !Sig + 3*noise +! base=(sum(s)-sp3n)/(NFFT1-3.0) !Noise per bin +! psig=sp3n-3*base !Sig only +! pnoise=(2500.0/df1)*base !Noise in 2500 Hz +! xsnrdb=db(psig/pnoise) + + a(1)=-fc1 + a(2:5)=0. + call twkfreq1(c,NZ,fs,a,cs) !Mix down by 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 + cs=cs/NZ + cs=cs*cs !Square the data + call four2a(cs,NZ,1,-1,1) !Compute squared spectrum + +! Find two peaks separated by baud + pmax=0. + fc2=0. + ic=nint(baud/df) + ja=nint(0.5*baud/df) + do j=-ja,ja + f2=j*df + ia=nint((f2-0.5*baud)/df) + if(ia.lt.0) ia=ia+NZ + ib=nint((f2+0.5*baud)/df) + p=real(cs(ia))**2 + aimag(cs(ia))**2 + & + real(cs(ib))**2 + aimag(cs(ib))**2 + if(p.gt.pmax) then + pmax=p + fc2=0.5*f2 + endif +! write(52,1200) f2,p,db(p) +!1200 format(f10.3,2f15.3) + enddo + sqf=sqf + (fc1+fc2-f0)**2 + a(1)=-(fc1+fc2) + a(2:5)=0. + call twkfreq1(c,NZ,fs,a,c) !Mix c down by fc1+fc2 + +! z=sum(c(1680:2095)*cb13)/208.0 !Get phase from Barker 13 vector +! z0=z/abs(z) +! c=c*conjg(z0) + +!---------------------------------------------------------------- DT + amax=0. + jpk=0 + 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) + jpk=j + endif +! write(53,1220) j,j*dt,z +!1220 format(i6,3f10.4) + enddo + xdt=jpk/fs + +!------------------------------------------------------------------ cpolyfit + ib=NSPS-1 + ib2=N2-1 + n=0 + do j=1,117 !First-pass demodulation + ia=ib+1 + ib=ia+N2-1 + zz(j)=sum(pp*c(ia:ib))/NSPS + if(abs(id(j)).eq.2) then !Save all sync symbols + n=n+1 + 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) + endif + if(j.le.116) then + zz(j+117)=sum(pp*c(ia+NSPS:ib+NSPS))/NSPS + endif + enddo + + aa=0. + bb=0. + nterms=0 + if(maxn.gt.0) then +! Fit sync info with a complex polynomial + npts=n + mode=0 + chisqa0=1.e30 + chisqb0=1.e30 + 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 + +!-------------------------------------------------------------- Soft Symbols + n=0 + do j=1,117 + 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) + if(abs(id(j)).eq.2) then + if(real(z)*id(j).lt.0) then + nhardsync=nhardsync+1 + nhardsync0=nhardsync0+1 + endif +! write(55,2002) j,id(j)/2,xx,z*id(j)/2 !Sync bit +!2002 format(2i5,3f10.3) + else + p=real(z) !Data bit + n=n+1 + rxdata(n)=p + ierr=0 + if(id(j)*p.lt.0) ierr=1 + nhard0=nhard0+ierr + nhard=nhard+ierr +! write(56,2003) j,id(j),n,ierr,nhard,xx,p*id(j),z +!2003 format(5i6,4f10.3) + endif + enddo + + do j=118,233 + xx=(j-116.5)*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=aimag(z) + n=n+1 + rxdata(n)=p + ierr=0 + if(id(j)*p.lt.0) ierr=1 + nhard=nhard+ierr + enddo + + 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 + nfe=nfe+ifer + write(58,1045) snrdb,nhard0,nhardsync0,niterations,nbadcrc,ifer, & + nterms,fc1+fc2-f0,xdt + if(ifer.eq.1) write(59,1045) snrdb,nhard0,nhardsync0,niterations, & + nbadcrc,ifer,nterms,fc1+fc2-f0,xdt +1045 format(f6.1,6i5,2f8.3) + enddo + fsigma=sqrt(sqf/iters) + 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 +1050 format(f6.1,i7,f8.4,f7.3,f8.2) + enddo + +999 end program msksim diff --git a/lib/fsk4hf/polyfit4.f90 b/lib/fsk4hf/polyfit4.f90 new file mode 100644 index 000000000..892cd2cd6 --- /dev/null +++ b/lib/fsk4hf/polyfit4.f90 @@ -0,0 +1,109 @@ +subroutine polyfit4(x,y,sigmay,npts,nterms,mode,a,chisqr) + + parameter (MAXN=20) + implicit real*8 (a-h,o-z) + real x(npts), y(npts), sigmay(npts), a(nterms),chisqr + real*8 sumx(2*MAXN-1), sumy(MAXN), array(MAXN,MAXN) + +! Accumulate weighted sums + nmax = 2*nterms-1 + sumx=0. + sumy=0. + chisq=0. + do i=1,npts + xi=x(i) + yi=y(i) + if(mode.lt.0) then + weight=1./abs(yi) + else if(mode.eq.0) then + weight=1 + else + weight=1./sigmay(i)**2 + end if + xterm=weight + do n=1,nmax + sumx(n)=sumx(n)+xterm + xterm=xterm*xi + enddo + yterm=weight*yi + do n=1,nterms + sumy(n)=sumy(n)+yterm + yterm=yterm*xi + enddo + chisq=chisq+weight*yi**2 + enddo + +! Construct matrices and calculate coefficients + do j=1,nterms + do k=1,nterms + n=j+k-1 + array(j,k)=sumx(n) + enddo + enddo + + delta=determ4(array,nterms) + if(delta.eq.0) then + chisqr=0. + a=0. + else + do l=1,nterms + do j=1,nterms + do k=1,nterms + n=j+k-1 + array(j,k)=sumx(n) + enddo + array(j,l)=sumy(j) + enddo + a(l)=determ4(array,nterms)/delta + enddo + +! Calculate chi square + + do j=1,nterms + chisq=chisq-2*a(j)*sumy(j) + do k=1,nterms + n=j+k-1 + chisq=chisq+a(j)*a(k)*sumx(n) + enddo + enddo + free=npts-nterms + chisqr=chisq/free + end if + + return +end subroutine polyfit4 + +real*8 function determ4(array,norder) + + parameter (MAXN=20) + implicit real*8 (a-h,o-z) + real*8 array(MAXN,MAXN) + + determ4=1. + do k=1,norder + if (array(k,k).ne.0) go to 41 + do j=k,norder + if(array(k,j).ne.0) go to 31 + enddo + determ4=0. + go to 60 + +31 do i=k,norder + s8=array(i,j) + array(i,j)=array(i,k) + array(i,k)=s8 + enddo + determ4=-1.*determ4 +41 determ4=determ4*array(k,k) + if(k.lt.norder) then + k1=k+1 + do i=k1,norder + do j=k1,norder + array(i,j)=array(i,j)-array(i,k)*array(k,j)/array(k,k) + enddo + enddo + end if + enddo + +60 return +end function determ4 diff --git a/lib/fsk4hf/twkfreq1.f90 b/lib/fsk4hf/twkfreq1.f90 new file mode 100644 index 000000000..e6aa8f8bf --- /dev/null +++ b/lib/fsk4hf/twkfreq1.f90 @@ -0,0 +1,26 @@ +subroutine twkfreq1(ca,npts,fsample,a,cb) + + complex ca(npts) + complex cb(npts) + complex w,wstep + real a(5) + data twopi/6.283185307/ + +! Mix the complex signal + w=1.0 + wstep=1.0 + x0=0.5*(npts+1) + s=2.0/npts + do i=1,npts + x=s*(i-x0) + p2=1.5*x*x - 0.5 + p3=2.5*(x**3) - 1.5*x + p4=4.375*(x**4) - 3.75*(x**2) + 0.375 + dphi=(a(1) + x*a(2) + p2*a(3) + p3*a(4) + p4*a(5)) * (twopi/fsample) + wstep=cmplx(cos(dphi),sin(dphi)) + w=w*wstep + cb(i)=w*ca(i) + enddo + + return +end subroutine twkfreq1 diff --git a/lib/fsk4hf/watterson.f90 b/lib/fsk4hf/watterson.f90 index aa7db106e..4b6799cb2 100644 --- a/lib/fsk4hf/watterson.f90 +++ b/lib/fsk4hf/watterson.f90 @@ -1,39 +1,41 @@ -subroutine watterson(c,delay,fspread) +subroutine watterson(c,fs,delay,fspread) - parameter (NZ=3456000) + parameter (NZ=3840) complex c(0:NZ-1) complex c2(0:NZ-1) complex cs1(0:NZ-1) complex cs2(0:NZ-1) - df=12000.0/NZ + nonzero=0 + df=fs/NZ if(fspread.gt.0.0) then do i=0,NZ-1 xx=gran() yy=gran() - cs1(i)=cmplx(xx,yy) + cs1(i)=0.707*cmplx(xx,yy) xx=gran() yy=gran() - cs2(i)=cmplx(xx,yy) + 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 f=i*df if(i.gt.NZ/2) f=(i-NZ)*df - x=(f/fspread)**2 + x=(f/(0.5*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 + 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 + if(p1.gt.0.0) nonzero=nonzero+1 +! write(62,3101) f,p1,p2,db(p1+1.e-12)-60,db(p2+1.e-12)-60 +!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) @@ -43,11 +45,13 @@ subroutine watterson(c,delay,fspread) 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)) + if(nonzero.gt.1) then + c(i)=0.5*(cs1(i)*c(i) + cs2(i)*c2(i)) + else + c(i)=0.5*(c(i) + c2(i)) + endif sq=sq + real(c(i))**2 + aimag(c(i))**2 ! write(61,3001) i/12000.0,c(i) !3001 format(3f12.6) diff --git a/lib/fsk4hf/wsprlf.f90 b/lib/fsk4hf/wsprlf.f90 index 803f3c027..efc3e8463 100644 --- a/lib/fsk4hf/wsprlf.f90 +++ b/lib/fsk4hf/wsprlf.f90 @@ -1,41 +1,45 @@ program wsprlf parameter (NN=121) !Total symbols -! parameter (NSPS=28672) !Samples per symbol - parameter (NSPS=28800) !Samples per symbol + parameter (NSPS=28800) !Samples per symbol @ fs=12000 Hz parameter (NZ=NSPS*NN) !Samples in waveform - parameter (NFFT=11*NSPS) - + character*8 arg complex c(0:NZ-1) - complex ct(0:NFFT-1) - real*8 twopi,f0,dt,phi,dphi - real s(0:NZ-1) - real h0(0:NSPS/2) - real h1(0:NSPS/2) - real p(0:NFFT-1) + real*8 twopi,fs,f0,dt,phi,dphi + real x(0:NZ-1) + real p(0:NZ/2) + real h0(0:NSPS/2) !Pulse shape, rising edge + real h1(0:NSPS/2) !Pulse shape, trailing edge real tmp(NN) - integer id(NN) + integer id(NN) !Generated data + integer ie(NN) !Differentially encoded data + data fs/12000.d0/ nargs=iargc() - if(nargs.ne.2) then - print*,'Usage: wsprlf f0 t1' + if(nargs.ne.3) then + print*,'Usage: wsprlf f0 t1 snr' goto 999 endif call getarg(1,arg) read(arg,*) f0 call getarg(2,arg) read(arg,*) t1 + call getarg(3,arg) + read(arg,*) snrdb - call random_number(tmp) !Generate random data - id=0 - where(tmp.ge.0.5) id=1 - id(1)=0 + call random_number(tmp) !Generate random bipolar data + id=1 + where(tmp.lt.0.5) id=-1 + ie(1)=1 + do i=2,NN !Differentially encode + ie(i)=id(i)*ie(i-1) + enddo n1=nint(t1*NSPS) twopi=8.d0*atan(1.d0) - do i=0,2*n1-1 + do i=0,2*n1-1 !Define the shape functions if(i.le.n1-1) then h0(i)=0.5*(1.0-cos(0.5*i*twopi/n1)) else @@ -45,67 +49,62 @@ program wsprlf if(t1.eq.0.0) h0=1 if(t1.eq.0.0) h1=1 - s=1. - s(0:n1-1)=h0(0:n1-1) !Leading edge of 1st pulse +! Shape the channel pulses + x=1. + x(0:n1-1)=h0(0:n1-1) !Leading edge of 1st pulse do j=2,NN !Leading edges - if(id(j).ne.id(j-1)) then + if(ie(j).ne.ie(j-1)) then ia=(j-1)*NSPS + 1 ib=ia+n1-1 - s(ia:ib)=h0(0:n1-1) + x(ia:ib)=h0(0:n1-1) endif enddo do j=1,NN-1 !Trailing edges - if(id(j+1).ne.id(j)) then + if(ie(j+1).ne.ie(j)) then ib=j*NSPS ia=ib-n1+1 - s(ia:ib)=h1(0:n1-1) + x(ia:ib)=h1(0:n1-1) endif enddo ib=NN*NSPS-1 ia=ib-n1+1 - s(ia:ib)=h1(0:n1-1) !Trailing edge of last pulse + x(ia:ib)=h1(0:n1-1) !Trailing edge of last pulse - dt=1.d0/12000.d0 + dt=1.d0/fs ts=dt*NSPS - baud=12000.0/NSPS + baud=fs/NSPS write(*,1000) baud,ts 1000 format('Baud:',f6.3,' Tsym:',f6.3) - phi=0. + dphi=twopi*f0*dt + phi=0.d0 i=-1 - do j=1,NN - x=1. - if(id(j).eq.1) x=-1. + do j=1,NN !Generate the baseband waveform + a=ie(j) do k=1,NSPS i=i+1 + x(i)=a*x(i) phi=phi+dphi if(phi.gt.twopi) phi=phi-twopi - c(i)=x*s(i)*cmplx(cos(phi),sin(phi)) - t=i*dt - sym=t/ts - if(sym.ge.10.0 .and. sym.le.20.0) write(13,3001) t, & - sym,s(i),c(i) -3001 format(5f12.6,i10) + xphi=phi + c(i)=x(i)*cmplx(cos(xphi),sin(xphi)) + sym=i*dt/ts + if(j.le.20) write(13,1010) sym,x(i),c(i) +1010 format(4f12.6) enddo enddo - p=0. - do iblk=1,11 - ia=(iblk-1)*NFFT - ib=ia+NFFT-1 - ct=c(ia:ib) - call four2a(ct,NFFT,1,-1,1) - do i=0,NFFT-1 - p(i)=p(i) + real(ct(i))**2 + aimag(ct(i))**2 - enddo + call four2a(c,NZ,1,-1,1) !To freq domain + df=fs/NZ + nh=NZ/2 + do i=0,nh + f=i*df + p(i)=real(c(i))**2 + aimag(c(i))**2 enddo - - p=cshift(p,NFFT/2)/maxval(p) - df=12000.0/NFFT - do i=0,NFFT-1 - f=i*df - 6000.0 - write(14,1020) f,p(i),10.0*log10(p(i)+1.e-12) -1020 format(f12.4,2e12.3) + p=p/maxval(p) + do i=0,nh !Save spectrum for plotting + write(14,1020) i*df,p(i),10.0*log10(p(i)+1.e-8) +1020 format(f10.3,2e12.3) enddo 999 end program wsprlf