diff --git a/lib/mskber.f90 b/lib/mskber.f90 index 2d2f049e5..4b31127c1 100644 --- a/lib/mskber.f90 +++ b/lib/mskber.f90 @@ -3,80 +3,119 @@ program mskber ! Generate an MSK waveform, pass it through an AWGN channel, apply coherent ! MSK receiver, and count number of errors at each Eb/No. - parameter (NSYM=100000) !Number of symbols to test - parameter (NSPS=6) !Samples per symbol - real ct(-NSPS:NSPS*NSYM-1) !cos(pi*t/2T) - real st(-NSPS:NSPS*NSYM-1) !sin(pi*t/2T) - real r(NSYM) !Random numbers to determine test bits - real xsym(0:NSYM-1) !Soft Rx symbols + parameter (MAXSYM=1000*1000) + parameter (NSPS=5) !Samples per symbol + real ct(-NSPS:NSPS*MAXSYM-1) !cos(pi*t/2T) + real st(-NSPS:NSPS*MAXSYM-1) !sin(pi*t/2T) + real r(0:MAXSYM-1) !Random numbers to set test bits + real xsym(0:MAXSYM-1) !Soft Rx symbols - complex xt(0:NSPS*(NSYM+1)-1) !Complex baseband Tx waveform - complex nt(0:NSPS*(NSYM+1)-1) !Generated AWGN channel noise - complex yt(0:NSPS*(NSYM+1)-1) !Received signal, yt = xt + fac*nt + complex xt(-NSPS:NSPS*MAXSYM-1) !Complex baseband Tx waveform + complex nt(-NSPS:NSPS*MAXSYM-1) !Generated AWGN channel noise + complex yt(-NSPS:NSPS*MAXSYM-1) !Received signal, yt = xt + fac*nt + complex cwave(-NSPS:NSPS*MAXSYM-1) !Audio waveform, Tx real part + complex z - integer sym0(0:NSYM-1) !Generated test bits - integer sym(0:NSYM-1) !Hard-copy received bits + integer sym0(0:MAXSYM-1) !Generated test bits + integer sym(0:MAXSYM-1) !Hard-copy received bits + integer sym1(0:7) + + character*12 arg + data sym1/1,1,0,0,0,1,1,1/ + + nargs=iargc() + if(nargs.ne.2) then + print*,'Usage: mskber nsym EbNo' + go to 999 + endif + call getarg(1,arg) + read(arg,*) nsym + call getarg(2,arg) + read(arg,*) EbNo pi=4.0*atan(1.0) - iz=NSPS*(NSYM+1) - do i=-NSPS,NSPS*NSYM-1 !Define ct, st arrays + do i=-NSPS,NSPS*nsym-1 !Define ct, st arrays t=i*pi/(2.0*NSPS) ct(i)=cos(t) st(i)=sin(t) enddo fac=1.0/sqrt(float(NSPS)) - do iEbNo=0,10 !Loop over a range of Eb/No + do iEbNo=0,10 !Loop over a range of Eb/No sym0=0 call random_number(r) - where(r.gt.0.5) sym0=1 !Generate random data bits - call mskmod(sym0,NSYM,NSPS,ct,st,xt) !Generate Tx waveform at baseband -! NB: In WSJT-X, will mix xt upward from 0 to 1500 Hz. + where(r(0:nsym-1).gt.0.5) sym0(0:nsym-1)=1 !Generate random data bits + if(nsym.eq.8) sym0(0:nsym-1)=sym1 + call mskmod(sym0,nsym,NSPS,ct,st,xt,cwave) !Generate Tx waveform - do i=0,iz-1 !Generate Gaussian noise + do i=-NSPS,NSPS*nsym-1 !Generate Gaussian noise xx=0.707*gran() yy=0.707*gran() nt(i)=cmplx(xx,yy) enddo fac_noise=10.0**(-iEbNo/20.0) - yt=xt + fac_noise*nt !Rx signal, with noise + if(EbNo.ne.0.0) fac_noise=10.0**(-EbNo/20.0) + yt=xt + fac_noise*nt !Rx signal, with noise - call mskdemod(yt,NSYM,NSPS,ct,st,xsym) !MSK demodulator + call mskdemod(yt,nsym,NSPS,ct,st,xsym) !MSK demodulator sym=0 where(xsym.gt.0.0) sym=1 + if(nsym.le.160 .and. EbNo.ne.0.0) then + write(*,1012) sym0(0:nsym-1) + if(nsym.gt.50) write(*,1012) + write(*,1012) sym(0:nsym-1) +1012 format(50i1) + do i=-nsps,nsps*nsym-1 + phi=i*2.0*pi*1500/12000.0 + z=cwave(i)*cmplx(cos(phi),sin(phi)) !Mix back to baseband + write(51,1014) float(i)/nsps,xt(i),abs(xt(i)),cwave(i),z +1014 format(8f8.4) + enddo + endif + ! Count the hard errors - nerr=count(sym(0:NSYM-1).ne.sym0(0:NSYM-1)) + nerr=count(sym(0:nsym-1).ne.sym0(0:nsym-1)) thber=0.5*erfc(10.0**(iEbNo/20.0)) - write(*,1000) iEbNo,thber,float(nerr)/NSYM -1000 format(i3,2f10.6) + xEbNo=iEbNo + if(EbNo.ne.0.0) xEbNo=EbNo + write(*,1000) xEbNo,thber,float(nerr)/nsym +1000 format(f6.1,2f10.6) + if(EbNo.ne.0.0) exit enddo -end program mskber +999 end program mskber -subroutine mskmod(sym,nsym,nsps,ct,st,xt) +subroutine mskmod(sym,nsym,nsps,ct,st,xt,cwave) ! Generate MSK Tx waveform at baseband. integer sym(0:nsym-1) !Hard-copy received bits - complex xt(0:nsps*(nsym+1)-1) !Complex baseband Tx waveform + complex xt(-nsps:nsps*nsym-1) !Complex baseband Tx waveform + complex cwave(-nsps:nsps*nsym-1) !Audio waveform, fc=1500 Hz. real ct(-nsps:nsps*nsym-1) !cos(pi*t/2T) real st(-nsps:nsps*nsym-1) !sin(pi*t/2T) - real ai(0:nsps*(nsym+1)-1) !Rectangular pulses for even symbols - real aq(0:nsps*(nsym+1)-1) !Rectangular pulses for odd symbols + real ai(-nsps:nsps*nsym-1) !Rectangular pulses for even symbols + real aq(-nsps:nsps*nsym-1) !Rectangular pulses for odd symbols + ai=0. + aq=0. fac=1.0/sqrt(float(nsps)) do j=0,nsym-1,2 - ia=j*nsps + ia=(j-1)*nsps ib=ia+2*nsps-1 ai(ia:ib)=2*sym(j)-1 !Even bits as rectangular pulses aq(ia+nsps:ib+nsps)=2*sym(j+1)-1 !Odd bits as rectangular pulses enddo - ai(ib+1:)=0 !Pad ai with zeros at end - aq(0:nsps-1)=0 !Pad aq with zeros at start - xt=fac*cmplx(ai*ct,aq*st) !Baseband Tx waveform + xt=fac*cmplx(ai*ct,aq*st) !Baseband Tx waveform + + twopi=8.0*atan(1.0) + do i=-nsps,nsps*nsym-1 + phi=i*twopi*1500/12000.0 + cwave(i)=xt(i)*cmplx(cos(phi),-sin(phi)) + enddo return end subroutine mskmod @@ -86,21 +125,21 @@ subroutine mskdemod(yt,nsym,nsps,ct,st,xsym) ! MSK demodulator ! Rx phase must be known and stable; symbol sync must be established. - complex yt(0:nsps*(nsym+1)-1) !Received signal + complex yt(-nsps:nsps*nsym-1) !Received signal real ct(-nsps:nsps*nsym-1) !cos(pi*t/2T) real st(-nsps:nsps*nsym-1) !sin(pi*t/2T) - real xe(0:nsps*(nsym+3)-1) !Temp array for received even symbols - real xo(0:nsps*(nsym+3)-1) !Temp array for received odd symbols + real xe(-nsps:nsps*nsym-1) !Temp array for received even symbols + real xo(-nsps:nsps*nsym-1) !Temp array for received odd symbols real xsym(0:nsym-1) !Soft Rx symbols iz=nsps*(nsym+1) - xe(0:iz-1)=real(yt)*ct - xo(0:iz-1)=aimag(yt)*st + xe(-nsps:nsps*nsym-1)=real(yt)*ct !Apply matched filters + xo(-nsps:nsps*nsym-1)=aimag(yt)*st do j=0,nsym-1,2 - ia=j*nsps + ia=(j-1)*nsps ib=ia+2*nsps-1 - xsym(j)=sum(xe(ia:ib)) !Integrate over 2 successive symbols - xsym(j+1)=sum(xo(ia+6:ib+6)) + xsym(j)=sum(xe(ia:ib)) !Integrate over 2 symbol lengths + xsym(j+1)=sum(xo(ia+nsps:ib+nsps)) enddo return