diff --git a/CMakeLists.txt b/CMakeLists.txt index 3e7e816bc..4a5b5c341 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -459,6 +459,7 @@ set (wsjt_FSRCS lib/ft8/genft8.f90 lib/genmsk_128_90.f90 lib/genmsk40.f90 + lib/fsk4hf/genft2.f90 lib/genqra64.f90 lib/ft8/genft8refsig.f90 lib/genwspr.f90 @@ -503,6 +504,8 @@ set (wsjt_FSRCS lib/msk144signalquality.f90 lib/msk144sim.f90 lib/mskrtd.f90 + lib/fsk4hf/ft2sim.f90 + lib/fsk4hf/ft2d.f90 lib/77bit/my_hash.f90 lib/wsprd/osdwspr.f90 lib/ft8/osd174_91.f90 @@ -1247,6 +1250,12 @@ target_link_libraries (ft8sim wsjt_fort wsjt_cxx) add_executable (msk144sim lib/msk144sim.f90 wsjtx.rc) target_link_libraries (msk144sim wsjt_fort wsjt_cxx) +add_executable (ft2sim lib/fsk4hf/ft2sim.f90 wsjtx.rc) +target_link_libraries (ft2sim wsjt_fort wsjt_cxx) + +add_executable (ft2d lib/fsk4hf/ft2d.f90 wsjtx.rc) +target_link_libraries (ft2d wsjt_fort wsjt_cxx) + endif(WSJT_BUILD_UTILS) # build the main application diff --git a/lib/fsk4hf/ft2_params.f90 b/lib/fsk4hf/ft2_params.f90 new file mode 100644 index 000000000..027d72bd8 --- /dev/null +++ b/lib/fsk4hf/ft2_params.f90 @@ -0,0 +1,12 @@ +! LDPC (128,90) code +parameter (KK=90) !Information bits (77 + CRC13) +parameter (ND=128) !Data symbols +parameter (NS=16) !Sync symbols (2x8) +parameter (NN=NS+ND) !Total channel symbols (144) +parameter (NSPS=160) !Samples per symbol at 12000 S/s +parameter (NZ=NSPS*NN) !Samples in full 1.92 s waveform (23040) +parameter (NMAX=3*12000) !Samples in iwave (36,000) +parameter (NFFT1=2*NSPS, NH1=NFFT1/2) !Length of FFTs for symbol spectra +parameter (NSTEP=NSPS/4) !Rough time-sync step size +parameter (NHSYM=NMAX/NSTEP-3) !Number of symbol spectra (1/4-sym steps) +parameter (NDOWN=10) !Downsample factor diff --git a/lib/fsk4hf/ft2d.f90 b/lib/fsk4hf/ft2d.f90 new file mode 100644 index 000000000..14e63c687 --- /dev/null +++ b/lib/fsk4hf/ft2d.f90 @@ -0,0 +1,325 @@ +program ft2d + + use crc + use packjt77 + include 'ft2_params.f90' + character arg*8,message*37,c77*77,infile*80,fname*16,datetime*11 + character*37 decodes(100) + character*120 data_dir + character*90 dmsg + complex c2(0:3*1200-1) !Complex waveform + complex cd(0:144*10-1) !Complex waveform + complex c1(0:9),c0(0:9) + complex ccor(0:1,144) + complex csum,cterm,cc0,cc1 + real*8 fMHz + + real rxdata(128),llr(128) !Soft symbols + real llr2(128) + real sbits(144),sbits1(144),sbits3(144) + real ps(0:8191),psbest(0:8191) + real candidates(100,2) + integer ihdr(11) + integer*2 iwave(NMAX) !Generated full-length waveform + integer*1 message77(77),apmask(128),cw(128) + integer*1 hbits(144),hbits1(144),hbits3(144) + logical unpk77_success + + fs=12000.0/NDOWN !Sample rate + dt=1/fs !Sample interval after downsample (s) + tt=NSPS*dt !Duration of "itone" symbols (s) + baud=1.0/tt !Keying rate for "itone" symbols (baud) + txt=NZ*dt !Transmission length (s) + twopi=8.0*atan(1.0) + h=0.8 !h=0.8 seems to be optimum for AWGN sensitivity (not for fading) + + dphi=twopi/2*baud*h*dt*16 ! dt*16 is samp interval after downsample + dphi0=-1*dphi + dphi1=+1*dphi + phi0=0.0 + phi1=0.0 + do i=0,9 + c1(i)=cmplx(cos(phi1),sin(phi1)) + c0(i)=cmplx(cos(phi0),sin(phi0)) + phi1=mod(phi1+dphi1,twopi) + phi0=mod(phi0+dphi0,twopi) + enddo + the=twopi*h/2.0 + cc1=cmplx(cos(the),-sin(the)) + cc0=cmplx(cos(the),sin(the)) + nargs=iargc() + if(nargs.lt.1) then + print*,'Usage: ft2d [-a ] [-f fMHz] [-c ncoh] 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 + ncoh=1 + npdi=16 + if(arg(1:2).eq.'-c') then + call getarg(iarg+1,arg) + read(arg,*) ncoh + iarg=iarg+2 + npdi=16/ncoh + endif +! write(*,*) 'ncoh: ',ncoh,' npdi: ',npdi + + xs1=0.0 + xs2=0.0 + fr1=0.0 + fr2=0.0 + nav=0 + ngood=0 + + do ifile=iarg,nargs + call getarg(ifile,infile) + j2=index(infile,'.wav') + open(10,file=infile,status='old',access='stream') + read(10,end=999) ihdr,iwave + read(infile(j2-4:j2-1),*) nutc + datetime=infile(j2-11:j2-1) + close(10) + + ndecodes=0 + ncand=1 + do icand=1,ncand + fc0=1500.0 + xsnr=1.0 + istart=6000+8 + call ft2_downsample(iwave,c2) ! downsample from 160s/Symbol to 10s/Symbol + + ib=istart/16 + cd=c2(ib:ib+144*10-1) + s2=sum(cd*conjg(cd))/(10*144) + cd=cd/sqrt(s2) + do nseq=1,7 + if( nseq.eq.1 ) then ! noncoherent single-symbol detection + sbits1=0.0 + do ibit=1,144 + ib=(ibit-1)*10 + ccor(1,ibit)=sum(cd(ib:ib+9)*conjg(c1(0:9))) + ccor(0,ibit)=sum(cd(ib:ib+9)*conjg(c0(0:9))) + sbits1(ibit)=abs(ccor(1,ibit))-abs(ccor(0,ibit)) + hbits1(ibit)=0 + if(sbits1(ibit).gt.0) hbits1(ibit)=1 + enddo + sbits=sbits1 + hbits=hbits1 + sbits3=sbits1 + hbits3=hbits1 + elseif( nseq.ge.2 ) then + nbit=2*nseq-1 + numseq=2**(nbit) + ps=0 + do ibit=nbit/2+1,144-nbit/2 + ps=0.0 + pmax=0.0 + do iseq=0,numseq-1 + csum=0.0 + cterm=1.0 + k=1 + do i=nbit-1,0,-1 + ibb=iand(iseq/(2**i),1) + csum=csum+ccor(ibb,ibit-(nbit/2+1)+k)*cterm + if(ibb.eq.0) cterm=cterm*cc0 + if(ibb.eq.1) cterm=cterm*cc1 + k=k+1 + enddo + ps(iseq)=abs(csum) + if( ps(iseq) .gt. pmax ) then + pmax=ps(iseq) + ibflag=1 + endif + enddo + if( ibflag .eq. 1 ) then + psbest=ps + ibflag=0 + endif + call getbitmetric(2**(nbit/2),psbest,numseq,sbits3(ibit)) + hbits3(ibit)=0 + if(sbits3(ibit).gt.0) hbits3(ibit)=1 + enddo + sbits=sbits3 + hbits=hbits3 + endif + rxdata(1:48)=sbits(9:56) + rxdata(49:128)=sbits(65:144) + rxav=sum(rxdata(1:128))/128.0 + rx2av=sum(rxdata(1:128)*rxdata(1:128))/128.0 + rxsig=sqrt(rx2av-rxav*rxav) + rxdata=rxdata/rxsig + sigma=0.90 + llr(1:128)=2*rxdata/(sigma*sigma) + apmask=0 + max_iterations=40 + ifer=0 + do ibias=0,0 + llr2=llr + if(ibias.eq.1) llr2=llr+0.4 + if(ibias.eq.2) llr2=llr-0.4 + call bpdecode128_90(llr2,apmask,max_iterations,message77,cw,nharderror,niterations) + if(nharderror.ge.0) exit + enddo + nhardmin=-1 + if(sum(message77).eq.0) cycle + if( nharderror.ge.0 ) then + write(c77,'(77i1)') message77(1:77) + call unpack77(c77,message,unpk77_success) + do i=1,ndecodes + if(decodes(i).eq.message) idupe=1 + enddo + if(idupe.eq.1) goto 888 + ndecodes=ndecodes+1 + decodes(ndecodes)=message + nsnr=nint(xsnr) + freq=fMHz + 1.d-6*(fc1+fbest) +1210 format(a11,2i4,f6.2,f12.7,2x,a22,i3) + write(*,1212) datetime(8:11),nsnr,xdt,freq,message,'*',idf,nseq,ijitter,nharderror,nhardmin +1212 format(a4,i4,f5.1,f11.6,2x,a22,a1,i5,i5,i5,i5,i5) + goto 888 + endif + enddo ! nseq +888 continue + enddo !candidate list + enddo !files + + write(*,1120) +1120 format("") + +999 end program ft2d + +subroutine getbitmetric(ib,ps,ns,xmet) + real ps(0:ns-1) + xm1=0 + xm0=0 + do i=0,ns-1 + if( iand(i/ib,1) .eq. 1 .and. ps(i) .gt. xm1 ) xm1=ps(i) + if( iand(i/ib,1) .eq. 0 .and. ps(i) .gt. xm0 ) xm0=ps(i) + enddo + xmet=xm1-xm0 + return +end subroutine getbitmetric + +subroutine downsample2(ci,f0,co) + parameter(NI=144*160,NH=NI/2,NO=NI/16) ! downsample from 200 samples per symbol to 10 + complex ci(0:NI-1),ct(0:NI-1) + complex co(0:NO-1) + fs=12000.0 + df=fs/NI + ct=ci + call four2a(ct,NI,1,-1,1) !c2c FFT to freq domain + i0=nint(f0/df) + ct=cshift(ct,i0) + co=0.0 + co(0)=ct(0) + b=8.0 + do i=1,NO/2 + arg=(i*df/b)**2 + filt=exp(-arg) + co(i)=ct(i)*filt + co(NO-i)=ct(NI-i)*filt + enddo + co=co/NO + call four2a(co,NO,1,1,1) !c2c FFT back to time domain + return +end subroutine downsample2 + +subroutine getcandidate2(c,npts,fs,fa,fb,ncand,candidates) + parameter(NDAT=200,NFFT1=120*12000/32,NH1=NFFT1/2,NFFT2=120*12000/320,NH2=NFFT2/2) + complex c(0:npts-1) !Complex waveform + complex cc(0:NFFT1-1) + complex csfil(0:NFFT2-1) + complex cwork(0:NFFT2-1) + real bigspec(0:NFFT2-1) + complex c2(0:NFFT1-1) !Short spectra + real s(-NH1+1:NH1) !Coarse spectrum + real ss(-NH1+1:NH1) !Smoothed coarse spectrum + real candidates(100,2) + integer indx(NFFT2-1) + logical first + data first/.true./ + save first,w,df,csfil + + if(first) then + df=10*fs/NFFT1 + csfil=cmplx(0.0,0.0) + do i=0,NFFT2-1 + csfil(i)=exp(-((i-NH2)/20.0)**2) + enddo + csfil=cshift(csfil,NH2) + call four2a(csfil,NFFT2,1,-1,1) + first=.false. + endif + + cc=cmplx(0.0,0.0) + cc(0:npts-1)=c; + call four2a(cc,NFFT1,1,-1,1) + cc=abs(cc)**2 + call four2a(cc,NFFT1,1,-1,1) + cwork(0:NH2)=cc(0:NH2)*conjg(csfil(0:NH2)) + cwork(NH2+1:NFFT2-1)=cc(NFFT1-NH2+1:NFFT1-1)*conjg(csfil(NH2+1:NFFT2-1)) + + call four2a(cwork,NFFT2,1,+1,1) + bigspec=cshift(real(cwork),-NH2) + il=NH2+fa/df + ih=NH2+fb/df + nnl=ih-il+1 + call indexx(bigspec(il:il+nnl-1),nnl,indx) + xn=bigspec(il-1+indx(nint(0.3*nnl))) + bigspec=bigspec/xn + ncand=0 + do i=il,ih + if((bigspec(i).gt.bigspec(i-1)).and. & + (bigspec(i).gt.bigspec(i+1)).and. & + (bigspec(i).gt.1.15).and.ncand.lt.100) then + ncand=ncand+1 + candidates(ncand,1)=df*(i-NH2) + candidates(ncand,2)=10*log10(bigspec(i))-30.0 + endif + enddo +! do i=1,ncand +! write(*,*) i,candidates(i,1),candidates(i,2) +! enddo + return +end subroutine getcandidate2 + +subroutine ft2_downsample(iwave,c) + +! Input: i*2 data in iwave() at sample rate 12000 Hz +! Output: Complex data in c(), sampled at 1200 Hz + + include 'ft2_params.f90' + parameter (NFFT2=NMAX/16) + integer*2 iwave(NMAX) + complex c(0:NMAX/16-1) + complex c1(0:NFFT2-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,NFFT2/2 + c1(i)=cx(i0+i) + c1(NFFT2-i)=cx(i0-i) + enddo + c1=c1/NFFT2 + call four2a(c1,NFFT2,1,1,1) !c2c FFT back to time domain + c=c1(0:NMAX/16-1) + return +end subroutine ft2_downsample + diff --git a/lib/fsk4hf/ft2sim.f90 b/lib/fsk4hf/ft2sim.f90 new file mode 100644 index 000000000..565b57474 --- /dev/null +++ b/lib/fsk4hf/ft2sim.f90 @@ -0,0 +1,139 @@ +program ft2sim + +! Generate simulated signals for experimental "FT2" mode + + use wavhdr + use packjt77 + include 'ft2_params.f90' !Set various constants + parameter (NWAVE=NN*NSPS) + type(hdr) h !Header for .wav file + character arg*12,fname*17 + character msg37*37,msgsent37*37 + character c77*77 + complex c0(0:NMAX-1) + complex c(0:NMAX-1) + real wave(NMAX) + integer itone(NN) + integer*1 msgbits(77) + integer*2 iwave(NMAX) !Generated full-length waveform + +! Get command-line argument(s) + nargs=iargc() + if(nargs.ne.8) then + print*,'Usage: ft2sim "message" f0 DT fdop del width nfiles snr' + print*,'Examples: ft2sim "K1ABC W9XYZ EN37" 1500.0 0.0 0.1 1.0 0 10 -18' + print*,' ft2sim "WA9XYZ/R KA1ABC/R FN42" 1500.0 0.0 0.1 1.0 0 10 -18' + print*,' ft2sim "K1ABC RR73; W9XYZ -11" 300 0 0 0 25 1 -10' + go to 999 + endif + call getarg(1,msg37) !Message to be transmitted + call getarg(2,arg) + read(arg,*) f0 !Frequency (only used for single-signal) + call getarg(3,arg) + read(arg,*) xdt !Time offset from nominal (s) + call getarg(4,arg) + read(arg,*) fspread !Watterson frequency spread (Hz) + call getarg(5,arg) + read(arg,*) delay !Watterson delay (ms) + call getarg(6,arg) + read(arg,*) width !Filter transition width (Hz) + call getarg(7,arg) + read(arg,*) nfiles !Number of files + call getarg(8,arg) + read(arg,*) snrdb !SNR_2500 + + nsig=1 + if(f0.lt.100.0) then + nsig=f0 + f0=1500 + endif + + nfiles=abs(nfiles) + twopi=8.0*atan(1.0) + fs=12000.0 !Sample rate (Hz) + dt=1.0/fs !Sample interval (s) + hmod=0.8 !Modulation index (0.5 is MSK, 1.0 is FSK) + tt=NSPS*dt !Duration of symbols (s) + baud=1.0/tt !Keying rate (baud) + bw=1.5*baud !Occupied bandwidth (Hz) + txt=NZ*dt !Transmission length (s) + bandwidth_ratio=2500.0/(fs/2.0) + sig=sqrt(2*bandwidth_ratio) * 10.0**(0.05*snrdb) + if(snrdb.gt.90.0) sig=1.0 + txt=NN*NSPS/12000.0 + + ! Source-encode, then get itone() + i3=-1 + n3=-1 + call pack77(msg37,i3,n3,c77) + read(c77,'(77i1)') msgbits + call genft2(msg37,0,msgsent37,itone,itype) + write(*,*) + write(*,'(a9,a37,3x,a7,i1,a1,i1)') 'Message: ',msgsent37,'i3.n3: ',i3,'.',n3 + write(*,1000) f0,xdt,txt,snrdb,bw +1000 format('f0:',f9.3,' DT:',f6.2,' TxT:',f6.1,' SNR:',f6.1, & + ' BW:',f5.1) + write(*,*) + if(i3.eq.1) then + write(*,*) ' mycall hiscall hisgrid' + write(*,'(28i1,1x,i1,1x,28i1,1x,i1,1x,i1,1x,15i1,1x,3i1)') msgbits(1:77) + else + write(*,'(a14)') 'Message bits: ' + write(*,'(77i1)') msgbits + endif + write(*,*) + write(*,'(a17)') 'Channel symbols: ' + write(*,'(79i1)') itone + write(*,*) + + call sgran() + + do ifile=1,nfiles + k=nint((xdt+0.5)/dt) + ia=k + phi=0.0 + c0=0.0 + do j=1,NN !Generate complex waveform + dphi=twopi*(f0*dt+(hmod/2.0)*(2*itone(j)-1)/real(NSPS)) + do i=1,NSPS + if(k.ge.0 .and. k.lt.NMAX) c0(k)=cmplx(cos(phi),sin(phi)) + k=k+1 + phi=mod(phi+dphi,twopi) + enddo + enddo + if(fspread.ne.0.0 .or. delay.ne.0.0) call watterson(c0,NMAX,NWAVE,fs,delay,fspread) + c=sig*c0 + + ib=k + wave=real(c) + peak=maxval(abs(wave(ia:ib))) + nslots=1 + if(width.gt.0.0) call filt8(f0,nslots,width,wave) + + if(snrdb.lt.90) then + do i=1,NMAX !Add gaussian noise at specified SNR + xnoise=gran() + wave(i)=wave(i) + xnoise + enddo + endif + + gain=100.0 + if(snrdb.lt.90.0) then + wave=gain*wave + else + datpk=maxval(abs(wave)) + fac=32766.9/datpk + wave=fac*wave + endif + if(any(abs(wave).gt.32767.0)) print*,"Warning - data will be clipped." + iwave=nint(wave) + h=default_header(12000,NMAX) + write(fname,1102) ifile +1102 format('000000_',i6.6,'.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,a17) + enddo +999 end program ft2sim diff --git a/lib/fsk4hf/genft2.f90 b/lib/fsk4hf/genft2.f90 new file mode 100644 index 000000000..d6d874e99 --- /dev/null +++ b/lib/fsk4hf/genft2.f90 @@ -0,0 +1,124 @@ +subroutine genft2(msg0,ichk,msgsent,i4tone,itype) +! s8 + 48bits + s8 + 80 bits = 144 bits (72ms message duration) +! +! Encode an MSK144 message +! Input: +! - msg0 requested message to be transmitted +! - ichk if ichk=1, return only msgsent +! if ichk.ge.10000, set imsg=ichk-10000 for short msg +! - msgsent message as it will be decoded +! - i4tone array of audio tone values, 0 or 1 +! - itype message type +! 1 = standard message "Call_1 Call_2 Grid/Rpt" +! 2 = type 1 prefix +! 3 = type 1 suffix +! 4 = type 2 prefix +! 5 = type 2 suffix +! 6 = free text (up to 13 characters) +! 7 = short message " Rpt" + + use iso_c_binding, only: c_loc,c_size_t + use packjt77 + character*37 msg0 + character*37 message !Message to be generated + character*37 msgsent !Message as it will be received + character*77 c77 + integer*4 i4tone(144) + integer*1 codeword(128) + integer*1 msgbits(77) + integer*1 bitseq(144) !Tone #s, data and sync (values 0-1) + integer*1 s8(8) + real*8 pp(12) + real*8 xi(864),xq(864),pi,twopi + data s8/0,1,1,1,0,0,1,0/ + equivalence (ihash,i1hash) + logical first,unpk77_success + data first/.true./ + save + + if(first) then + first=.false. + nsym=128 + pi=4.0*atan(1.0) + twopi=8.*atan(1.0) + do i=1,12 + pp(i)=sin((i-1)*pi/12) + enddo + endif + + message(1:37)=' ' + itype=1 + if(msg0(1:1).eq.'@') then !Generate a fixed tone + read(msg0(2:5),*,end=1,err=1) nfreq !at specified frequency + go to 2 +1 nfreq=1000 +2 i4tone(1)=nfreq + else + message=msg0 + + do i=1, 37 + if(ichar(message(i:i)).eq.0) then + message(i:37)=' ' + exit + endif + enddo + do i=1,37 !Strip leading blanks + if(message(1:1).ne.' ') exit + message=message(i+1:) + enddo + + if(message(1:1).eq.'<') then + i2=index(message,'>') + i1=0 + if(i2.gt.0) i1=index(message(1:i2),' ') + if(i1.gt.0) then + call genmsk40(message,msgsent,ichk,i4tone,itype) + if(itype.lt.0) go to 999 + i4tone(41)=-40 + go to 999 + endif + endif + + i3=-1 + n3=-1 + call pack77(message,i3,n3,c77) + call unpack77(c77,msgsent,unpk77_success) !Unpack to get msgsent + + if(ichk.eq.1) go to 999 + read(c77,"(77i1)") msgbits + call encode_128_90(msgbits,codeword) + +!Create 144-bit channel vector: +!8-bit sync word + 48 bits + 8-bit sync word + 80 bits + bitseq=0 + bitseq(1:8)=s8 + bitseq(9:56)=codeword(1:48) + bitseq(57:64)=s8 + bitseq(65:144)=codeword(49:128) + + i4tone=bitseq + +! bitseq=2*bitseq-1 +! xq(1:6)=bitseq(1)*pp(7:12) !first bit is mapped to 1st half-symbol on q +! do i=1,71 +! is=(i-1)*12+7 +! xq(is:is+11)=bitseq(2*i+1)*pp +! enddo +! xq(864-5:864)=bitseq(1)*pp(1:6) !last half symbol +! do i=1,72 +! is=(i-1)*12+1 +! xi(is:is+11)=bitseq(2*i)*pp +! enddo +! Map I and Q to tones. +! i4tone=0 +! do i=1,72 +! i4tone(2*i-1)=(bitseq(2*i)*bitseq(2*i-1)+1)/2; +! i4tone(2*i)=-(bitseq(2*i)*bitseq(mod(2*i,144)+1)-1)/2; +! enddo + endif + +! Flip polarity +! i4tone=-i4tone+1 + +999 return +end subroutine genft2 diff --git a/lib/fsk4hf/mskhfsim.f90 b/lib/fsk4hf/mskhfsim.f90 deleted file mode 100644 index 86140e1ce..000000000 --- a/lib/fsk4hf/mskhfsim.f90 +++ /dev/null @@ -1,194 +0,0 @@ -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=1152/72) !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 (3776) - 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 c0(0:NZ-1) !Complex waveform - complex zz(NS+ND) !Complex symbol values (intermediate) - complex z - real xnoise(0:NZ-1) !Generated random noise - real ynoise(0:NZ-1) !Generated random noise - 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 ierror(NS+ND) - 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: mskhfsim f0(Hz) delay(ms) fspread(Hz) maxn iters snr(dB)' - print*,'Example: mskhfsim 0 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,3*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 - c=cbb - if(delay.ne.0.0 .or. fspread.ne.0.0) then - call watterson(c,NZ,fs,delay,fspread) - endif - 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 - - call getfc1(c,fc1) !First approx for freq - call getfc2(c,csync,fc1,fc2,fc3) !Refined freq - sqf=sqf + (fc1+fc2-f0)**2 - -!NB: Measured performance is about equally good using fc2 or fc3 here: - a(1)=-(fc1+fc2) - a(2:5)=0. - call twkfreq1(c,NZ,fs,a,c) !Mix c down by fc1+fc2 - -! The following may not be necessary? -! z=sum(c(1680:2095)*cb13)/208.0 !Get phase from Barker 13 vector -! z0=z/abs(z) -! c=c*conjg(z0) - -!---------------------------------------------------------------- DT -! Not presently used: - 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 - enddo - xdt=jpk/fs - - nterms=maxn - c0=c - do itry=1,10 - idf=itry/2 - if(mod(itry,2).eq.0) idf=-idf - nhard0=0 - nhardsync0=0 - ifer=1 - a(1)=idf*0.01 - a(2:5)=0. - call twkfreq1(c0,NZ,fs,a,c) !Mix c0 into c - call cpolyfit(c,pp,id,maxn,aa,bb,zz,nhs) - call msksoftsym(zz,aa,bb,id,nterms,ierror,rxdata,nhard0,nhardsync0) - if(nhardsync0.gt.12) cycle - 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 - ifer=0 - call bpdecode168(llr,apmask,max_iterations,decoded,niterations,cw) - nbadcrc=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 -! if(ifer.eq.0) write(67,1301) snrdb,itry,idf,niterations, & -! nhardsync0,nhard0 -!1301 format(f6.1,5i6) - if(ifer.eq.0) exit - enddo !Freq dither loop - nhard=nhard+nhard0 - nhardsync=nharsdync+nhardsync0 - nfe=nfe+ifer - 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/spb.m b/lib/fsk4hf/spb.m index ed1761925..9bb164506 100644 --- a/lib/fsk4hf/spb.m +++ b/lib/fsk4hf/spb.m @@ -70,14 +70,14 @@ endfunction # M-ary PSK Block Coded Modulation," Igal Sason and Gil Weichman, # doi: 10.1109/EEEI.2006.321097 #------------------------------------------------------------------------------- -N=174 -K=75 +N=128 +K=90 R=K/N delta=0.01; [ths,fval,info,output]=fzero(@f1,[delta,pi/2-delta], optimset ("jacobian", "off")); -for ebnodb=-6:0.5:4 +for ebnodb=-3:0.5:4 ebno=10^(ebnodb/10.0); esno=ebno*R; A=sqrt(2*esno); diff --git a/lib/fsk4hf/spb_128_90.dat b/lib/fsk4hf/spb_128_90.dat new file mode 100644 index 000000000..9e32e28e9 --- /dev/null +++ b/lib/fsk4hf/spb_128_90.dat @@ -0,0 +1,19 @@ +N = 128 +K = 90 +R = 0.70312 +-3.000000 0.000341 +-2.500000 0.001513 +-2.000000 0.006049 +-1.500000 0.021280 +-1.000000 0.064283 +-0.500000 0.162755 +0.000000 0.338430 +0.500000 0.571867 +1.000000 0.791634 +1.500000 0.930284 +2.000000 0.985385 +2.500000 0.998258 +3.000000 0.999893 +3.500000 0.999997 +4.000000 1.000000 +