From 28c09d4a760834339d0a03c2f0f9dd8b067f9b95 Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Tue, 16 May 2017 19:20:14 +0000 Subject: [PATCH] Add some new test programs. git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7680 ab8295b8-cf94-4d9e-aec4-7959e3be5d79 --- CMakeLists.txt | 11 ++ lib/fsk4hf/genwspr_fsk8.f90 | 45 +++++++ lib/fsk4hf/osd300.f90 | 25 ++-- lib/fsk4hf/spec8.f90 | 31 +++++ lib/fsk4hf/wspr5d.f90 | 4 +- lib/fsk4hf/wspr_fsk8_downsample.f90 | 29 +++++ lib/fsk4hf/wspr_fsk8_params.f90 | 13 +++ lib/fsk4hf/wspr_fsk8_sim.f90 | 64 ++++++++++ lib/fsk4hf/wspr_fsk8_wav.f90 | 45 +++++++ lib/fsk4hf/wspr_fsk8d.f90 | 175 ++++++++++++++++++++++++++++ 10 files changed, 428 insertions(+), 14 deletions(-) create mode 100644 lib/fsk4hf/genwspr_fsk8.f90 create mode 100644 lib/fsk4hf/spec8.f90 create mode 100644 lib/fsk4hf/wspr_fsk8_downsample.f90 create mode 100644 lib/fsk4hf/wspr_fsk8_params.f90 create mode 100644 lib/fsk4hf/wspr_fsk8_sim.f90 create mode 100644 lib/fsk4hf/wspr_fsk8_wav.f90 create mode 100644 lib/fsk4hf/wspr_fsk8d.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 1947f21bd..435920737 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -395,6 +395,7 @@ set (wsjt_FSRCS lib/fsk4hf/genmskhf.f90 lib/fsk4hf/genwspr5.f90 lib/fsk4hf/genwsprlf.f90 + lib/fsk4hf/genwspr_fsk8.f90 lib/fsk4hf/getfc1.f90 lib/fsk4hf/getfc2.f90 lib/fsk4hf/getfc1w.f90 @@ -475,6 +476,7 @@ set (wsjt_FSRCS lib/shell.f90 lib/fsk4hf/spec4.f90 lib/spec64.f90 + lib/fsk4hf/spec8.f90 lib/spec9f.f90 lib/stdmsg.f90 lib/subtract65.f90 @@ -506,7 +508,10 @@ set (wsjt_FSRCS lib/wqencode.f90 lib/fsk4hf/wspr5d.f90 lib/fsk4hf/wspr5sim.f90 + lib/fsk4hf/wspr_fsk8d.f90 + lib/fsk4hf/wspr_fsk8_sim.f90 lib/fsk4hf/wspr5_downsample.f90 + lib/fsk4hf/wspr_fsk8_downsample.f90 lib/fsk4hf/wspr5_wav.f90 lib/fsk4hf/wsprlfsim.f90 lib/wspr_downsample.f90 @@ -1147,6 +1152,12 @@ 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) +add_executable (wspr_fsk8d lib/fsk4hf/wspr_fsk8d.f90 wsjtx.rc) +target_link_libraries (wspr_fsk8d wsjt_fort wsjt_cxx) + +add_executable (wspr_fsk8_sim lib/fsk4hf/wspr_fsk8_sim.f90 wsjtx.rc) +target_link_libraries (wspr_fsk8_sim wsjt_fort wsjt_cxx) + add_executable (mskhfsim lib/fsk4hf/mskhfsim.f90 wsjtx.rc) target_link_libraries (mskhfsim wsjt_fort wsjt_cxx) diff --git a/lib/fsk4hf/genwspr_fsk8.f90 b/lib/fsk4hf/genwspr_fsk8.f90 new file mode 100644 index 000000000..2ef6ace30 --- /dev/null +++ b/lib/fsk4hf/genwspr_fsk8.f90 @@ -0,0 +1,45 @@ +subroutine genwspr_fsk8(msg,msgsent,itone) + +! Encode a WSPR-LF 8-FSK message, producing array itone(). + + use crc + include 'wspr_fsk8_params.f90' + + character*22 msg,msgsent + character*60 cbits + integer*1,target :: idat(9) + integer*1 msgbits(KK),codeword(3*ND) + integer itone(NN) + integer icos7(0:6) + data icos7/2,5,6,0,4,1,3/ !Costas 7x7 tone pattern + + 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 + +! Message structure: S7 D100 S7 + itone(1:7)=icos7 + itone(NN-6:NN)=icos7 + do j=1,ND + i=3*j -2 + itone(j+7)=codeword(i)*4 + codeword(i+1)*2 + codeword(i+2) + enddo + + return +end subroutine genwspr_fsk8 diff --git a/lib/fsk4hf/osd300.f90 b/lib/fsk4hf/osd300.f90 index 5a858938a..434447067 100644 --- a/lib/fsk4hf/osd300.f90 +++ b/lib/fsk4hf/osd300.f90 @@ -158,21 +158,24 @@ return end subroutine mrbencode subroutine nextpat(mi,k,iorder,iflag) -integer*1 mi(k),ms(k) + integer*1 mi(k),ms(k) ! generate the next test error pattern ind=-1 do i=1,k-1 - if( mi(i).eq.0 .and. mi(i+1).eq.1) ind=i + if( mi(i).eq.0 .and. mi(i+1).eq.1) ind=i enddo - ms=0 - ms(1:ind-1)=mi(1:ind-1) - ms(ind)=1 - ms(ind+1)=0 - if( ind+1 .lt. k ) then - nz=iorder-sum(ms) - ms(k-nz+1:k)=1 - endif - mi=ms + if(ind.ge.1 .and. ind.le.k) then !### JHT temporary ### Correct ??? + ms=0 + ms(1:ind-1)=mi(1:ind-1) + ms(ind)=1 + ms(ind+1)=0 + if( ind+1 .lt. k ) then + nz=iorder-sum(ms) + ms(k-nz+1:k)=1 + endif + mi=ms + endif !### JHT temporary ### iflag=ind + return end subroutine nextpat diff --git a/lib/fsk4hf/spec8.f90 b/lib/fsk4hf/spec8.f90 new file mode 100644 index 000000000..7efc993d5 --- /dev/null +++ b/lib/fsk4hf/spec8.f90 @@ -0,0 +1,31 @@ +subroutine spec8(c,s,savg) + + include 'wspr_fsk8_params.f90' + complex c(0:NZ-1) + complex c1(0:NSPS-1) + real s(0:NH2,NN) + real savg(0:NH2) + + fs=12000.0/NDOWN + df=fs/NSPS + savg=0. + do j=1,NN + ia=(j-1)*NSPS + ib=ia + NSPS-1 + c1(0:NSPS-1)=c(ia:ib) + c1(NSPS:)=0. + call four2a(c1,NSPS,1,-1,1) + do i=0,NH2 + s(i,j)=real(c1(i))**2 + aimag(c1(i))**2 + enddo + savg=savg+s(0:NH2,j) + enddo + s=s/NZ + savg=savg/(NN*NZ) + do i=0,NH2 + write(31,3101) i*df,savg(i) +3101 format(f10.3,f12.3) + enddo + + return +end subroutine spec8 diff --git a/lib/fsk4hf/wspr5d.f90 b/lib/fsk4hf/wspr5d.f90 index cdda78214..6597af58d 100644 --- a/lib/fsk4hf/wspr5d.f90 +++ b/lib/fsk4hf/wspr5d.f90 @@ -120,7 +120,7 @@ program wspr5d go to 999 endif close(10) - fa=100.0 + fa=102.0 fb=150.0 call getfc1w(c,fs,fa,fb,fc1,xsnr) !First approx for freq call getfc2w(c,csync,fs,fc1,fc2,fc3) !Refined freq @@ -215,5 +215,3 @@ program wspr5d 1120 format("") 999 end program wspr5d - - include 'wspr5_downsample.f90' diff --git a/lib/fsk4hf/wspr_fsk8_downsample.f90 b/lib/fsk4hf/wspr_fsk8_downsample.f90 new file mode 100644 index 000000000..ed173bc05 --- /dev/null +++ b/lib/fsk4hf/wspr_fsk8_downsample.f90 @@ -0,0 +1,29 @@ +subroutine wspr_fsk8_downsample(iwave,c) + +! Input: i*2 data in iwave() at sample rate 12000 Hz +! Output: Complex data in c(), sampled at 4=500 Hz + + include 'wspr_fsk8_params.f90' + parameter (NMAX=240*12000,NFFTD=NMAX/24) + integer*2 iwave(NMAX) + complex c(0:NZ-1) + complex c1(0:NFFTD-1) + complex cx(0:NMAX/2) + real x(NMAX) + equivalence (x,cx) + + df=12000.0/NMAX + x=iwave + call four2a(x,NMAX,1,-1,0) !r2c FFT to freq domain + i0=nint(1500.0/df) + c1(0)=cx(i0) + do i=1,NFFTD/2 + c1(i)=cx(i0+i) + c1(NFFTD-i)=cx(i0-i) + enddo + c1=c1/NFFTD + call four2a(c1,NFFTD,1,1,1) !c2c FFT back to time domain + c=c1(0:NZ-1) + + return +end subroutine wspr_fsk8_downsample diff --git a/lib/fsk4hf/wspr_fsk8_params.f90 b/lib/fsk4hf/wspr_fsk8_params.f90 new file mode 100644 index 000000000..3bc17fa98 --- /dev/null +++ b/lib/fsk4hf/wspr_fsk8_params.f90 @@ -0,0 +1,13 @@ +! LDPC (300,60) code +parameter (NDOWN=24) !Downsample factor +parameter (KK=60) !Information bits (50 + CRC10) +parameter (ND=100) !Data symbols +parameter (NS=14) !Sync symbols (2 @ Costas 7x7) +parameter (NN=NS+ND) !Total symbols (114) +parameter (NSPS0=24576) !Samples per symbol at 12000 S/s +parameter (NSPS=NSPS0/NDOWN) !Sam/sym, downsampled (1024) +parameter (N7=7*NSPS) !Samples in Costas 7x7 array (7168) +parameter (NZ=NSPS*NN) !Samples in downsampled waveform (116,736) +parameter (NZMAX=NSPS0*NN) !Samples in *.wav (2,801,664) +parameter (NFFT1=4*NSPS,NH1=NFFT1/2) +parameter (NH2=NSPS/2) diff --git a/lib/fsk4hf/wspr_fsk8_sim.f90 b/lib/fsk4hf/wspr_fsk8_sim.f90 new file mode 100644 index 000000000..4049cdf18 --- /dev/null +++ b/lib/fsk4hf/wspr_fsk8_sim.f90 @@ -0,0 +1,64 @@ +program wspr_fsk8_sim + +! Generate simulated data for a 4-minute "WSPR-LF" mode using 8-FSK. +! Output is saved to a *.wav file. + + use wavhdr + include 'wspr_fsk8_params.f90' !Set various constants + parameter (NMAX=300*12000) + type(hdr) h !Header for .wav file + character arg*12,fname*16 + character msg*22,msgsent*22 + integer itone(NN) + integer*2 iwave(NMAX) !Generated full-length waveform + +! Get command-line argument(s) + nargs=iargc() + if(nargs.ne.5) then + print*,'Usage: wspr_fsk8_sim "message" f0 DT nfiles snr' + print*,'Example: wspr_fsk8_sim "K1ABC FN42 30" 1640 0.0 10 -33' + go to 999 + endif + call getarg(1,msg) !Message to be transmitted + call getarg(2,arg) + read(arg,*) f0 !Freq relative to WSPR-band center (Hz) + call getarg(3,arg) + read(arg,*) xdt !Time offset from nominal (s) + call getarg(4,arg) + read(arg,*) nfiles !Number of files + call getarg(5,arg) + read(arg,*) snrdb !SNR_2500 + + twopi=8.0*atan(1.0) + fs=12000.0/NDOWN !Sample rate after downsampling + dt=1.0/fs !Sample interval (s) + tt=NSPS*dt !Duration of symbols (s) + baud=1.0/tt !Keying rate + bw=8*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 + txt=NN*NSPS0/12000.0 + + call genwspr_fsk8(msg,msgsent,itone) !Encode the message, get itone + write(*,1000) f0,xdt,txt,snrdb,bw,msgsent +1000 format('f0:',f9.3,' DT:',f6.2,' TxT:',f6.1,' SNR:',f6.1, & + ' BW:',f4.1,2x,a22) + + do ifile=1,nfiles + call wspr_fsk8_wav(baud,xdt,f0,itone,snrdb,iwave) + h=default_header(12000,NMAX) + write(fname,1102) ifile +1102 format('000000_',i4.4,'.wav') + open(10,file=fname,status='unknown',access='stream') + write(10) h,iwave !Save to *.wav file + close(10) + write(*,1110) ifile,xdt,f0,snrdb,fname +1110 format(i4,f7.2,f8.2,f7.1,2x,a16) + enddo + +999 end program wspr_fsk8_sim + + include 'wspr_fsk8_wav.f90' + diff --git a/lib/fsk4hf/wspr_fsk8_wav.f90 b/lib/fsk4hf/wspr_fsk8_wav.f90 new file mode 100644 index 000000000..00039ae9e --- /dev/null +++ b/lib/fsk4hf/wspr_fsk8_wav.f90 @@ -0,0 +1,45 @@ +subroutine wspr_fsk8_wav(baud,xdt,f0,itone,snrdb,iwave) + +! Generate iwave() from itone(). + + include 'wspr_fsk8_params.f90' + parameter (NMAX=240*12000) + integer itone(NN) + integer*2 iwave(NMAX) + real*8 twopi,dt,dphi,phi + real dat(NMAX) + + twopi=8.d0*atan(1.d0) + dt=1.d0/12000.d0 + + dat=0. + if(snrdb.lt.90) then + do i=1,NMAX + dat(i)=gran() !Generate gaussian noise + enddo + bandwidth_ratio=2500.0/6000.0 + sig=sqrt(2*bandwidth_ratio)*10.0**(0.05*snrdb) + else + sig=1.0 + endif + + phi=0.d0 + k=nint(xdt/dt) + do j=1,NN + dphi=twopi*(f0+ itone(j)*baud)*dt + if(k.eq.0) phi=-dphi + do i=1,NSPS0 + k=k+1 + phi=phi+dphi + if(phi.gt.twopi) phi=phi-twopi + xphi=phi + if(k.gt.0 .and. k.le.NMAX) dat(k)=dat(k) + sig*sin(xphi) + enddo + enddo + fac=32767.0 + rms=100.0 + if(snrdb.ge.90.0) iwave=nint(fac*dat) + if(snrdb.lt.90.0) iwave=nint(rms*dat) + + return +end subroutine wspr_fsk8_wav diff --git a/lib/fsk4hf/wspr_fsk8d.f90 b/lib/fsk4hf/wspr_fsk8d.f90 new file mode 100644 index 000000000..c706a6ad9 --- /dev/null +++ b/lib/fsk4hf/wspr_fsk8d.f90 @@ -0,0 +1,175 @@ +program wspr_fsk8d + +! WSPR-LF is a potential WSPR-like mode intended for use at LF and MF. +! This version uses 4-minute T/R sequences, an LDPC (300,60) code, +! 8-FSK modulation, and noncoherent demodulation. This decoder reads +! data from *.wav files. + +! 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. Find two 7x7 Costas arrays to get xdt and fc2 +! 4. Mix from fc2 to 0, compute aligned symbol spectra +! 5. Get soft bits from symbol spectra + +! Still to do: find and decode more than one signal in the specified passband. + + include 'wspr_fsk8_params.f90' + parameter (NMAX=240*12000) + character arg*8,message*22,cbits*50,infile*80,datetime*11 + character*120 data_dir + complex csync(0:N7-1) !Sync symbols for Costas 7x7 array + complex c1(0:2*N7-1) + complex c2(0:2*N7-1) + complex c(0:NMAX-1) !Complex waveform + real*8 fMHz + real rxdata(3*ND),llr(3*ND) !Soft symbols + real a(5) !For twkfreq1 + real s(0:NH2,NN) + real savg(0:NH2) + real ss(0:N7) + real ps(0:7) + integer ihdr(11) + integer*2 iwave(NMAX) !Generated full-length waveform + integer*1 idat(7) + integer*1 decoded(KK),apmask(3*ND),cw(3*ND) + integer icos7(0:6) + data icos7/2,5,6,0,4,1,3/ !Costas 7x7 tone pattern + + nargs=iargc() + if(nargs.lt.2) then + print*,'Usage: wspr_fsk8d [-a ] [-f fMHz] file1 [file2 ...]' + go to 999 + endif + iarg=1 + data_dir="." + call getarg(iarg,arg) + if(arg(1:2).eq.'-a') then + call getarg(iarg+1,data_dir) + iarg=iarg+2 + endif + call getarg(iarg,arg) + if(arg(1:2).eq.'-f') then + call getarg(iarg+1,arg) + read(arg,*) fMHz + iarg=iarg+2 + endif + + open(13,file=trim(data_dir)//'/ALL_WSPR.TXT',status='unknown', & + position='append') + + twopi=8.0*atan(1.0) + fs=NSPS*12000.0/NSPS0 !Sample rate after downsampling (Hz) + dt=1.0/fs !Sample interval (s) + ts=NSPS*dt !Symbol duration (s) + baud=1.0/ts !Keying rate (Hz) + txt=NZ*dt !Transmission length (s) + + phi=0. + k=-1 + do j=0,6 + dphi=twopi*baud*icos7(j)*dt + do i=1,NSPS + phi=phi+dphi + if(phi.gt.twopi) phi=phi-twopi + k=k+1 + csync(k)=cmplx(cos(phi),sin(phi)) + enddo + enddo + + do ifile=iarg,nargs + call getarg(ifile,infile) + open(10,file=infile,status='old',access='stream') + read(10,end=999) ihdr,iwave + close(10) + j2=index(infile,'.wav') + read(infile(j2-4:j2-1),*) nutc + datetime=infile(j2-11:j2-1) + call wspr_fsk8_downsample(iwave,c) + c(NZ:)=0. + + j0=0 +! Need to get jpk ==> xdt + j0b=j0+107*NSPS + c1(0:N7-1)=c(j0:j0+N7-1)*conjg(csync) + c1(N7:)=0. + c2(0:N7-1)=c(j0b:j0b+N7-1)*conjg(csync) + c2(N7:)=0. + call four2a(c1,2*N7,1,-1,1) + call four2a(c2,2*N7,1,-1,1) + pmax=0. + df1=fs/(2*N7) + do i=0,N7 + f=1500.0 + i*df1 + p=1.e-9*(real(c1(i))**2 + aimag(c1(i))**2 + & + real(c2(i))**2 + aimag(c2(i))**2) + ss(i)=p + if(p.gt.pmax) then + pmax=p + ipk=i + endif + write(32,3201) f,ss(i) +3201 format(f12.4,e12.4) + enddo + + sp3n=(ss(ipk-1)+ss(ipk)+ss(ipk+1)) !Sig + 3*noise + call pctile(ss,N7,45,base) + psig=sp3n-3*base !Sig only + pnoise=(2500.0/df1)*base !Noise in 2500 Hz + xsnr=db(psig/pnoise) !SNR from sync tones + + f0=ipk*df1 + a(1)=-f0 + a(2:5)=0. + call twkfreq1(c,NZ,fs,a,c) !Mix from f0 down to 0 + call spec8(c,s,savg) !Get symbol spectra + + j0=0 + do j=1,ND + k=j+7 + ps=s(0:7,k) +! ps=sqrt(ps) !### ??? ### + r1=max(ps(1),ps(3),ps(5),ps(7))-max(ps(0),ps(2),ps(4),ps(6)) + r2=max(ps(2),ps(3),ps(6),ps(7))-max(ps(0),ps(1),ps(4),ps(5)) + r4=max(ps(4),ps(5),ps(6),ps(7))-max(ps(0),ps(1),ps(2),ps(3)) + rxdata(3*j-2)=r4 + rxdata(3*j-1)=r2 + rxdata(3*j)=r1 + enddo + + rxav=sum(rxdata)/ND + rx2av=sum(rxdata*rxdata)/ND + rxsig=sqrt(rx2av-rxav*rxav) + rxdata=rxdata/rxsig + s0=0.84 + llr=2.0*rxdata/(s0*s0) + apmask=0 + max_iterations=40 + ifer=0 + call bpdecode300(llr,apmask,max_iterations,decoded,niterations,cw) + if(niterations.lt.0) call osd300(llr,4,decoded,niterations,cw) + nbadcrc=0 + if(niterations.ge.0) call chkcrc10(decoded,nbadcrc) + if(niterations.lt.0 .or. nbadcrc.ne.0) ifer=1 + nsnr=nint(xsnr) + freq=fMHz + 1.d-6*f0 + nfdot=0 + 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) + write(*,1112) datetime(8:11),nsnr,xdt,freq,nfdot,message +1112 format(a4,i4,f5.1,f11.6,i3,2x,a22) + endif + write(13,1110) datetime,0,nsnr,xdt,freq,message,nfdot +1110 format(a11,2i4,f6.2,f12.7,2x,a22,i3) + enddo ! ifile loop + write(*,1120) +1120 format("") + +999 end program wspr_fsk8d +