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' character arg*8,message*22,cbits*50,infile*80,fname*16,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:NMAXD-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 ss0(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.7) then print*,'Usage: wspr_fsk8d -d db -f fMHz -a data_dir file1 [file2 ...]' go to 999 endif call getarg(1,arg) if(arg(1:3).ne.'-d ') go to 999 call getarg(2,arg) read(arg,*) degrade rxbw=3000. call getarg(3,arg) if(arg(1:3).ne.'-f ') go to 999 call getarg(4,arg) read(arg,*) fMHz call getarg(5,arg) if(arg(1:3).ne.'-a ') go to 999 call getarg(6,data_dir) 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=7,nargs call getarg(ifile,infile) open(10,file=infile,status='old',access='stream') j1=index(infile,'.c4') j2=index(infile,'.wav') if(j1.gt.0) then read(10,end=999) fname,ntrmin,fMHz,c(0:NZ-1) read(fname(8:11),*) nutc write(datetime,'(i11)') nutc else if(j2.gt.0) then read(10,end=999) ihdr,iwave read(infile(j2-4:j2-1),*) nutc datetime=infile(j2-11:j2-1) if(degrade.gt.0.0) call degrade_snr(iwave,NMAX,degrade,rxbw) call wspr_fsk8_downsample(iwave,c) else print*,'Wrong file format?' go to 999 endif close(10) pmax=0. df1=fs/(2*N7) ia=nint(100.0/df1) ib=nint(150.0/df1) ipk=0 jpk=0 ! Get xdt and f0 from the sync arrays do j0=0,1000,50 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) do i=0,N7 p=1.e-9*(real(c1(i))**2 + aimag(c1(i))**2 + & real(c2(i))**2 + aimag(c2(i))**2) ss(i)=p enddo do i=ia,ib p=ss(i) if(p.gt.pmax) then pmax=p ipk=i jpk=j0 endif enddo if(jpk.eq.j0) ss0=ss enddo f0=ipk*df1 xdt=jpk*dt - 1.0 sp3n=(ss0(ipk-1)+ss0(ipk)+ss0(ipk+1)) !Sig + 3*noise call pctile(ss0,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 if(jpk.ge.0) c(0:NMAXD-1-jpk)=c(jpk:NMAXD-1) 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 do j=1,ND k=j+7 ps=s(0:7,k) ps=sqrt(ps) !### ??? ### ! ps=log(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=1.1 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 freq=1.d-6*(f0+1500) 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