diff --git a/CMakeLists.txt b/CMakeLists.txt index e33ffe408..1b1ba38ce 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -462,7 +462,7 @@ set (wsjt_FSRCS lib/ft8/genft8.f90 lib/genmsk_128_90.f90 lib/genmsk40.f90 - lib/fsk4hf/genft2.f90 + lib/fsk4hf/genft4.f90 lib/genqra64.f90 lib/ft8/genft8refsig.f90 lib/genwspr.f90 @@ -508,8 +508,8 @@ set (wsjt_FSRCS lib/msk144signalquality.f90 lib/msk144sim.f90 lib/mskrtd.f90 - lib/fsk4hf/ft2sim.f90 - lib/fsk4hf/ft2d.f90 + lib/fsk4hf/ft4sim.f90 + lib/fsk4hf/ft4d.f90 lib/77bit/my_hash.f90 lib/wsprd/osdwspr.f90 lib/ft8/osd174_91.f90 @@ -1261,11 +1261,11 @@ 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 (ft4sim lib/fsk4hf/ft4sim.f90 wsjtx.rc) +target_link_libraries (ft4sim wsjt_fort wsjt_cxx) -add_executable (ft2d lib/fsk4hf/ft2d.f90 wsjtx.rc) -target_link_libraries (ft2d wsjt_fort wsjt_cxx) +add_executable (ft4d lib/fsk4hf/ft4d.f90 wsjtx.rc) +target_link_libraries (ft4d wsjt_fort wsjt_cxx) endif(WSJT_BUILD_UTILS) diff --git a/lib/fsk4hf/ft2d.f90 b/lib/fsk4hf/ft2d.f90 index e924366c6..fda1f1826 100644 --- a/lib/fsk4hf/ft2d.f90 +++ b/lib/fsk4hf/ft2d.f90 @@ -13,6 +13,7 @@ program ft2d complex c1(0:9),c0(0:9) complex ccor(0:1,144) complex csum,cterm,cc0,cc1,csync1,csync2 + complex csync(16),csl(0:159) real*8 fMHz real a(5) @@ -26,9 +27,10 @@ program ft2d integer*2 iwave(NMAX) !Generated full-length waveform integer*1 message77(77),apmask(128),cw(128) integer*1 hbits(144),hbits1(144),hbits3(144) - integer*1 s16(16) + integer*1 s16(16),s45(45) logical unpk77_success data s16/0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0/ + data s45/0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,1,0,0,1,1,1,1,0,0,1,0,0,0,1,1,0,1,0,0,0,1,1,1,0,0/ fs=12000.0/NDOWN !Sample rate dt=1/fs !Sample interval after downsample (s) @@ -36,7 +38,7 @@ program ft2d 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) + h=0.800 !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 @@ -52,6 +54,18 @@ program ft2d the=twopi*h/2.0 cc1=cmplx(cos(the),-sin(the)) cc0=cmplx(cos(the),sin(the)) + + k=0 + do j=1,16 + dphi1=(2*s16(j)-1)*dphi + phi1=0.0 + do i=0,9 + csl(k)=cmplx(cos(phi1),sin(phi1)) + phi1=mod(phi1+dphi1,twopi) + k=k+1 + enddo + enddo + nargs=iargc() if(nargs.lt.1) then print*,'Usage: ft2d [-a ] [-f fMHz] file1 [file2 ...]' @@ -89,6 +103,22 @@ program ft2d xsnr=1.0 if( f0.le.375.0 .or. f0.ge.(5000.0-375.0) ) cycle call ft2_downsample(iwave,f0,c2) ! downsample from 160s/Symbol to 10s/Symbol + +!c2=c2/sqrt(sum(abs(c2(0:NMAX/16-1)))) +!ishift=-1 +!rccbest=-99. +!do is=0,435 +!rcc=0.0 +! do id=10,10 +! rcc=rcc+abs(sum(conjg(c2(is:is+159-id))*c2(is+id:is+159)*csl(0:159-id)*conjg(csl(id:159)))) +! enddo +! if(rcc.gt.rccbest) then +! rccbest=rcc +! ishift=is +! endif +!write(21,*) is,rcc +!enddo + ! 750 samples/second here ibest=-1 sybest=-99. @@ -102,9 +132,10 @@ program ft2d csync1=0. cterm=1 do ib=1,16 +! do ib=1,45 i1=(ib-1)*10+is - i2=i1+136*10 if(s16(ib).eq.1) then +! if(s45(ib).eq.1) then csync1=csync1+sum(cb(i1:i1+9)*conjg(c1(0:9)))*cterm cterm=cterm*cc1 else @@ -121,13 +152,18 @@ program ft2d enddo a=0. +!dfbest=1500.0-f0 a(1)=-dfbest + call twkfreq1(c2,NMAX/16,fs,a,cb) + +!ibest=197 ib=ibest + cd=cb(ib:ib+144*10-1) s2=sum(cd*conjg(cd))/(10*144) cd=cd/sqrt(s2) - do nseq=1,5 + do nseq=1,4 if( nseq.eq.1 ) then ! noncoherent single-symbol detection sbits1=0.0 do ibit=1,144 @@ -178,7 +214,7 @@ program ft2d hbits=hbits3 endif nsync_qual=count(hbits(1:16).eq.s16) - if(nsync_qual.lt.10) exit +! if(nsync_qual.lt.10) exit rxdata=sbits(17:144) rxav=sum(rxdata(1:128))/128.0 rx2av=sum(rxdata(1:128)*rxdata(1:128))/128.0 @@ -186,7 +222,11 @@ program ft2d rxdata=rxdata/rxsig sigma=0.80 llr(1:128)=2*rxdata/(sigma*sigma) +!xllrmax=maxval(abs(llr)) +!write(*,*) ifile,icand,nseq,nsync_qual apmask=0 +!apmask(1:29)=1 +!llr(1:29)=xllrmax*(2*s45(17:45)-1) max_iterations=40 do ibias=0,0 llr2=llr @@ -209,8 +249,8 @@ program ft2d nsnr=nint(xsnr) freq=f0+dfbest 1210 format(a11,2i4,f6.2,f12.7,2x,a22,i3) - write(*,1212) datetime(8:11),nsnr,ibest/750.0,freq,message,'*',nseq,nharderror -1212 format(a4,i4,f5.1,f11.1,2x,a22,a1,i5,i5) + write(*,1212) datetime(8:11),nsnr,ibest/750.0,freq,message,'*',nseq,nharderror,nsync_qual +1212 format(a4,i4,2x,f5.3,f11.1,2x,a22,a1,i5,i5,i5) goto 888 endif enddo ! nseq diff --git a/lib/fsk4hf/ft2sim.f90 b/lib/fsk4hf/ft2sim.f90 index 5c0e07431..ea0518a52 100644 --- a/lib/fsk4hf/ft2sim.f90 +++ b/lib/fsk4hf/ft2sim.f90 @@ -48,7 +48,7 @@ program ft2sim 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) + hmod=0.800 !Modulation index (0.5 is MSK, 1.0 is FSK) tt=NSPS*dt !Duration of symbols (s) baud=1.0/tt !Keying rate (baud) txt=NZ*dt !Transmission length (s) diff --git a/lib/fsk4hf/ft4_params.f90 b/lib/fsk4hf/ft4_params.f90 new file mode 100644 index 000000000..3f3e9769d --- /dev/null +++ b/lib/fsk4hf/ft4_params.f90 @@ -0,0 +1,14 @@ +! FT4 37.5 baud - 26.67 ms symbol duration +! LDPC (128,90) code + +parameter (KK=90) !Information bits (77 + CRC13) +parameter (ND=64) !Data symbols +parameter (NS=12) !Sync symbols (12) +parameter (NN=NS+ND) !Total channel symbols (76) +parameter (NSPS=320) !Samples per symbol at 12000 S/s +parameter (NZ=NSPS*NN) !Samples in full 1.92 s waveform (23040) +parameter (NMAX=2.5*12000) !Samples in iwave (36,000) +parameter (NFFT1=640, 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=16) !Downsample factor diff --git a/lib/fsk4hf/ft4d.f90 b/lib/fsk4hf/ft4d.f90 new file mode 100644 index 000000000..22c81cec8 --- /dev/null +++ b/lib/fsk4hf/ft4d.f90 @@ -0,0 +1,317 @@ +program ft4d + + use crc + use packjt77 + include 'ft4_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 cd2(0:NMAX/16-1) !Complex waveform + complex cb(0:NMAX/16-1) + complex cd(0:76*10-1) !Complex waveform + complex c3(0:19),c2(0:19),c1(0:19),c0(0:19) + complex ccor(0:3,76) + complex csum,cterm,cc0,cc1,csync1,csync2 + complex csync(12) + real*8 fMHz + + real a(5) + real rxdata(128),llr(128) !Soft symbols + real llr2(128) + real sbits(152),sbits1(152),sbits3(152) + real ps(0:8191),psbest(0:8191) + real candidates(100,2) + real savg(NH1),sbase(NH1) + integer ihdr(11) + integer*2 iwave(NMAX) !Generated full-length waveform + integer*1 message77(77),apmask(128),cw(128) + integer*1 hbits(152),hbits1(152),hbits3(152) + integer*1 s12(12) + logical unpk77_success + data s12/0,0,0,1,1,1,1,1,1,0,0,0/ + + 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=1.000 !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.5*dphi + dphi1=-0.5*dphi + dphi2=+0.5*dphi + dphi3=+1.5*dphi + phi0=0.0 + phi1=0.0 + phi2=0.0 + phi3=0.0 + do i=0,9 + c3(i)=cmplx(cos(phi3),sin(phi3)) + c2(i)=cmplx(cos(phi2),sin(phi2)) + c1(i)=cmplx(cos(phi1),sin(phi1)) + c0(i)=cmplx(cos(phi0),sin(phi0)) + phi3=mod(phi1+dphi3,twopi) + phi2=mod(phi1+dphi2,twopi) + phi1=mod(phi1+dphi1,twopi) + phi0=mod(phi0+dphi0,twopi) + enddo + the=twopi*h/2.0 + cc3=cmplx(cos(3*the),-sin(3*the)) + cc2=cmplx(cos(the),-sin(the)) + cc1=cmplx(cos(the),sin(the)) + cc0=cmplx(cos(3*the),-sin(3*the)) + + nargs=iargc() + if(nargs.lt.1) then + print*,'Usage: ft4d [-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 + ncoh=1 + + 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) + candidates=0.0 + ncand=0 + call getcandidates2(iwave,375.0,3000.0,0.2,2200.0,100,savg,candidates,ncand,sbase) + ndecodes=0 + do icand=1,ncand + f0=candidates(icand,1) + xsnr=1.0 + if( f0.le.375.0 .or. f0.ge.(5000.0-375.0) ) cycle + call ft4_downsample(iwave,f0,cd2) ! downsample from 160s/Symbol to 10s/Symbol + +! 750 samples/second here + ibest=-1 + sybest=-99. + dfbest=-1. + do if=-60,+60 + df=if + a=0. + a(1)=-df + call twkfreq1(cd2,NMAX/16,fs,a,cb) + do is=0,380 + csync1=0. + cterm=1 + do ib=1,12 + i1=(ib-1)*10+is + if(s12(ib).eq.0) then + csync1=csync1+sum(cb(i1:i1+19)*conjg(c0(0:19)))*cterm + cterm=cterm*cc0 + else + csync1=csync1+sum(cb(i1:i1+19)*conjg(c3(0:19)))*cterm + cterm=cterm*cc3 + endif + enddo + if(abs(csync1).gt.sybest) then + ibest=is + sybest=abs(csync1) + dfbest=df + endif + enddo + enddo + + a=0. +!dfbest=1500.0-f0 + a(1)=-dfbest + + call twkfreq1(c2,NMAX/16,fs,a,cb) + +!ibest=197 + ib=ibest +write(*,*) f0,f0+dfbest,ibest +goto 888 + cd=cb(ib:ib+144*10-1) + s2=sum(cd*conjg(cd))/(10*144) + cd=cd/sqrt(s2) + do nseq=1,1 + if( nseq.eq.1 ) then ! noncoherent single-symbol detection + sbits1=0.0 + do ibit=1,76 + 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 + nsync_qual=count(hbits(1:16).eq.s16) +! if(nsync_qual.lt.10) exit + rxdata=sbits(17: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.80 + llr(1:128)=2*rxdata/(sigma*sigma) +!xllrmax=maxval(abs(llr)) +!write(*,*) ifile,icand,nseq,nsync_qual + apmask=0 +!apmask(1:29)=1 +!llr(1:29)=xllrmax*(2*s45(17:45)-1) + max_iterations=40 + 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 + if(sum(message77).eq.0) cycle + if( nharderror.ge.0 ) then + write(c77,'(77i1)') message77(1:77) + call unpack77(c77,1,message,unpk77_success) + idupe=0 + 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=f0+dfbest +1210 format(a11,2i4,f6.2,f12.7,2x,a22,i3) + write(*,1212) datetime(8:11),nsnr,ibest/750.0,freq,message,'*',nseq,nharderror,nsync_qual +1212 format(a4,i4,2x,f5.3,f11.1,2x,a22,a1,i5,i5,i5) + goto 888 + endif + enddo ! nseq +888 continue + enddo !candidate list + enddo !files + + write(*,1120) +1120 format("") + +999 end program ft4d + +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 ft4_downsample(iwave,f0,c) + +! Input: i*2 data in iwave() at sample rate 12000 Hz +! Output: Complex data in c(), sampled at 1200 Hz + + include 'ft4_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) + + BW=4.0*75 + df=12000.0/NMAX + x=iwave + call four2a(x,NMAX,1,-1,0) !r2c FFT to freq domain + ibw=nint(BW/df) + i0=nint(f0/df) + c1=0. + c1(0)=cx(i0) + do i=1,NFFT2/2 + arg=(i-1)*df/bw + win=exp(-arg*arg) + c1(i)=cx(i0+i)*win + c1(NFFT2-i)=cx(i0-i)*win + 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 ft4_downsample + diff --git a/lib/fsk4hf/ft4sim.f90 b/lib/fsk4hf/ft4sim.f90 new file mode 100644 index 000000000..c5fb3e572 --- /dev/null +++ b/lib/fsk4hf/ft4sim.f90 @@ -0,0 +1,158 @@ +program ft4sim + +! Generate simulated signals for experimental "FT4" mode + + use wavhdr + use packjt77 + include 'ft4_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) + real dphi(0:NMAX-1) + real pulse(960) + 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: ft4sim "message" f0 DT fdop del width nfiles snr' + print*,'Examples: ft4sim "K1ABC W9XYZ EN37" 1500.0 0.0 0.1 1.0 0 10 -18' + print*,' ft4sim "WA9XYZ/R KA1ABC/R FN42" 1500.0 0.0 0.1 1.0 0 10 -18' + print*,' ft4sim "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 + + nfiles=abs(nfiles) + twopi=8.0*atan(1.0) + fs=12000.0 !Sample rate (Hz) + dt=1.0/fs !Sample interval (s) + hmod=1.000 !Modulation index (0.5 is MSK, 1.0 is FSK) + tt=NSPS*dt !Duration of symbols (s) + baud=1.0/tt !Keying rate (baud) + 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 genft4(msg37,0,msgsent37,itone) + write(*,*) + write(*,'(a9,a37,3x,a7,i1,a1,i1)') 'Message: ',msgsent37,'i3.n3: ',i3,'.',n3 + write(*,1000) f0,xdt,txt,snrdb +1000 format('f0:',f9.3,' DT:',f6.2,' TxT:',f6.1,' SNR:',f6.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(*,'(76i1)') itone + write(*,*) + + call sgran() + +! The filtered frequency pulse + do i=1,960 + tt=(i-480.5)/320.0 + pulse(i)=gfsk_pulse(1.0,tt) + enddo + +! Define the instantaneous frequency waveform + dphi_peak=twopi*(hmod/2.0)/real(NSPS) + dphi=0.0 + do j=1,NN + ib=(j-1)*320 + ie=ib+960-1 + dphi(ib:ie)=dphi(ib:ie)+dphi_peak*pulse*(2*itone(j)-3) + enddo + + phi=0.0 + c0=0.0 + dphi=dphi+twopi*f0*dt + do j=0,NMAX-1 + c0(j)=cmplx(cos(phi),sin(phi)) + phi=mod(phi+dphi(j),twopi) + enddo + + c0(0:319)=c0(0:319)*(1.0-cos(twopi*(/(i,i=0,319)/)/640.0) )/2.0 + c0(77*320:77*320+319)=c0(77*320:77*320+319)*(1.0+cos(twopi*(/(i,i=0,319)/)/640.0 ))/2.0 + c0(78*320:)=0. + + k=nint((xdt+0.25)/dt) + c0=cshift(c0,-k) + ia=k + +do i=0,NMAX-1 +write(21,*) i,real(c0(i)),imag(c0(i)),dphi(i) +enddo + + do ifile=1,nfiles + c=c0 + if(fspread.ne.0.0 .or. delay.ne.0.0) call watterson(c,NMAX,NWAVE,fs,delay,fspread) + c=sig*c + + 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 ft4sim diff --git a/lib/fsk4hf/genft4.f90 b/lib/fsk4hf/genft4.f90 new file mode 100644 index 000000000..da7ff40c3 --- /dev/null +++ b/lib/fsk4hf/genft4.f90 @@ -0,0 +1,68 @@ +subroutine genft4(msg0,ichk,msgsent,i4tone) +! s12 + 64symbols = 76 channel symbols 2.027s message duration +! +! Encode an FT4 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,1,2,3} + + 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(76) + integer*1 codeword(128) + integer*1 msgbits(77) + integer*1 s12(12) + real*8 xi(864),xq(864),pi,twopi + data s12/0,0,0,3,3,3,3,3,3,0,0,0/ + logical unpk77_success + + twopi=8.*atan(1.0) + pi=twopi/2.0 + + 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 + + i3=-1 + n3=-1 + call pack77(message,i3,n3,c77) + call unpack77(c77,0,msgsent,unpk77_success) !Unpack to get msgsent + + if(ichk.eq.1) go to 999 + read(c77,"(77i1)") msgbits + call encode_128_90(msgbits,codeword) + +! Grayscale mapping: +! bits tone +! 00 0 +! 01 1 +! 11 2 +! 10 3 + +!Create 144-bit channel vector: + i4tone(1:12)=s12 + do i=1,64 + is=codeword(2*i-1)+2*codeword(2*i) + if(is.le.1) i4tone(12+i)=is + if(is.eq.2) i4tone(12+i)=3 + if(is.eq.3) i4tone(12+i)=2 + enddo + +999 return +end subroutine genft4 diff --git a/lib/fsk4hf/getcandidates4.f90 b/lib/fsk4hf/getcandidates4.f90 new file mode 100644 index 000000000..428b8faf5 --- /dev/null +++ b/lib/fsk4hf/getcandidates4.f90 @@ -0,0 +1,63 @@ +subroutine getcandidates4(id,fa,fb,syncmin,nfqso,maxcand,savg,candidate, & + ncand,sbase) + +! For now, hardwired to find the largest peak in the average spectrum + + include 'ft2_params.f90' + real s(NH1,NHSYM) + real savg(NH1),savsm(NH1) + real sbase(NH1) + real x(NFFT1) + complex cx(0:NH1) + real candidate(3,maxcand) + integer*2 id(NMAX) + integer*1 s8(8) + integer indx(NH1) + data s8/0,1,1,1,0,0,1,0/ + equivalence (x,cx) + +! Compute symbol spectra, stepping by NSTEP steps. + savg=0. + tstep=NSTEP/12000.0 + df=12000.0/NFFT1 !3.125 Hz + fac=1.0/300.0 + do j=1,NHSYM + ia=(j-1)*NSTEP + 1 + ib=ia+NSPS-1 + x(1:NSPS)=fac*id(ia:ib) + x(NSPS+1:)=0. + call four2a(x,NFFT1,1,-1,0) !r2c FFT + do i=1,NH1 + s(i,j)=real(cx(i))**2 + aimag(cx(i))**2 + enddo + savg=savg + s(1:NH1,j) !Average spectrum + enddo + savsm=0. + do i=2,NH1-1 + savsm(i)=sum(savg(i-1:i+1))/3. + enddo + + nfa=fa/df + nfb=fb/df + np=nfb-nfa+1 + indx=0 + call indexx(savsm(nfa:nfb),np,indx) + xn=savsm(nfa+indx(nint(0.3*np))) + savsm=savsm/xn + imax=-1 + xmax=-99. + do i=2,NH1-1 + if(savsm(i).gt.savsm(i-1).and. & + savsm(i).gt.savsm(i+1).and. & + savsm(i).gt.xmax) then + xmax=savsm(i) + imax=i + endif + enddo + f0=imax*df + if(xmax.gt.1.2) then + ncand=ncand+1 + candidate(1,ncand)=f0 + endif +return +end subroutine getcandidates4