diff --git a/CMakeLists.txt b/CMakeLists.txt index 4a5b5c341..29abb9064 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -546,6 +546,7 @@ set (wsjt_FSRCS lib/sync4.f90 lib/sync64.f90 lib/sync65.f90 + lib/fsk4hf/getcandidates2.f90 lib/ft8/sync8.f90 lib/ft8/sync8d.f90 lib/sync9.f90 diff --git a/lib/fsk4hf/ft2_params.f90 b/lib/fsk4hf/ft2_params.f90 index 027d72bd8..351119b4b 100644 --- a/lib/fsk4hf/ft2_params.f90 +++ b/lib/fsk4hf/ft2_params.f90 @@ -5,8 +5,8 @@ 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 (NMAX=2.5*12000) !Samples in iwave (36,000) +parameter (NFFT1=400, 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 +parameter (NDOWN=16) !Downsample factor diff --git a/lib/fsk4hf/ft2d.f90 b/lib/fsk4hf/ft2d.f90 index be43b7760..68cdcde32 100644 --- a/lib/fsk4hf/ft2d.f90 +++ b/lib/fsk4hf/ft2d.f90 @@ -7,23 +7,28 @@ program ft2d character*37 decodes(100) character*120 data_dir character*90 dmsg - complex c2(0:3*1200-1) !Complex waveform + complex c2(0:NMAX/16-1) !Complex waveform + complex cb(0:NMAX/16-1) 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 + complex csum,cterm,cc0,cc1,csync1,csync2 real*8 fMHz + real a(5) 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) + 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(144),hbits1(144),hbits3(144) + integer*1 s16(16) logical unpk77_success + data s16/0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0/ fs=12000.0/NDOWN !Sample rate dt=1/fs !Sample interval after downsample (s) @@ -81,7 +86,8 @@ program ft2d fr2=0.0 nav=0 ngood=0 - +fsum=0.0 +fsum2=0.0 do ifile=iarg,nargs call getarg(ifile,infile) j2=index(infile,'.wav') @@ -91,16 +97,55 @@ program ft2d datetime=infile(j2-11:j2-1) close(10) + call getcandidates2(iwave,100.0,3000.0,0.2,2200.0,100,savg,candidates,ncand,sbase) + ndecodes=0 ncand=1 do icand=1,ncand - fc0=1500.0 + f0=candidates(icand,1) 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) + call ft2_downsample(iwave,f0,c2) ! downsample from 160s/Symbol to 10s/Symbol +ibest=-1 +sybest=-99. +dfbest=-1. +do if=-15,+15 + df=if + a=0. + a(1)=-df +call twkfreq1(c2,NMAX/16,fs,a,cb) +! 750 samples/second here + do is=0,374 + csync1=0. + cterm=1 + do ib=1,16 + i1=(ib-1)*10+is + i2=i1+136*10 + if(s16(ib).eq.1) then + csync1=csync1+sum(cb(i1:i1+9)*conjg(c1(0:9)))*cterm + cterm=cterm*cc1 + else + csync1=csync1+sum(cb(i1:i1+9)*conjg(c0(0:9)))*cterm + cterm=cterm*cc0 + endif + enddo + if(abs(csync1).gt.sybest) then + ibest=is + sybest=abs(csync1) + dfbest=df + endif + enddo +enddo +freq=f0+dfbest +fsum=fsum+freq +fsum2=fsum+freq*freq +a=0. +a(1)=-dfbest +!write(*,*) 'dfbest ',dfbest +!dfbest=0.0 +!ibest=187 +call twkfreq1(c2,NMAX/16,fs,a,cb) + ib=ibest + cd=cb(ib:ib+144*10-1) s2=sum(cd*conjg(cd))/(10*144) cd=cd/sqrt(s2) do nseq=1,4 @@ -153,8 +198,7 @@ program ft2d sbits=sbits3 hbits=hbits3 endif - rxdata(1:48)=sbits(9:56) - rxdata(49:128)=sbits(65:144) + 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) @@ -175,6 +219,7 @@ program ft2d if( nharderror.ge.0 ) then write(c77,'(77i1)') message77(1:77) call unpack77(c77,message,unpk77_success) + idupe=0 do i=1,ndecodes if(decodes(i).eq.message) idupe=1 enddo @@ -182,10 +227,9 @@ program ft2d 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) + write(*,1212) datetime(8:11),nsnr,ibest/750.0,freq,message,'*',idf,nseq,ijitter,nharderror,nhardmin +1212 format(a4,i4,f5.1,f11.1,2x,a22,a1,i5,i5,i5,i5,i5) goto 888 endif enddo ! nseq @@ -195,7 +239,9 @@ program ft2d write(*,1120) 1120 format("") - +favg=fsum/1000.0 +fstd=sqrt(fsum2/1000.0-favg*favg) +write(*,*) "Mean, std frequency: ",favg,fstd 999 end program ft2d subroutine getbitmetric(ib,ps,ns,xmet) @@ -234,66 +280,7 @@ subroutine downsample2(ci,f0,co) 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) +subroutine ft2_downsample(iwave,f0,c) ! Input: i*2 data in iwave() at sample rate 12000 Hz ! Output: Complex data in c(), sampled at 1200 Hz @@ -310,7 +297,7 @@ subroutine ft2_downsample(iwave,c) df=12000.0/NMAX x=iwave call four2a(x,NMAX,1,-1,0) !r2c FFT to freq domain - i0=nint(1500.0/df) + i0=nint(f0/df) c1(0)=cx(i0) do i=1,NFFT2/2 c1(i)=cx(i0+i) diff --git a/lib/fsk4hf/ft2sim.f90 b/lib/fsk4hf/ft2sim.f90 index 565b57474..6addbe8d8 100644 --- a/lib/fsk4hf/ft2sim.f90 +++ b/lib/fsk4hf/ft2sim.f90 @@ -89,7 +89,7 @@ program ft2sim call sgran() do ifile=1,nfiles - k=nint((xdt+0.5)/dt) + k=nint((xdt+0.25)/dt) ia=k phi=0.0 c0=0.0 diff --git a/lib/fsk4hf/genft2.f90 b/lib/fsk4hf/genft2.f90 index d6d874e99..f416c3616 100644 --- a/lib/fsk4hf/genft2.f90 +++ b/lib/fsk4hf/genft2.f90 @@ -9,13 +9,8 @@ subroutine genft2(msg0,ichk,msgsent,i4tone,itype) ! - 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" +! 1 = 77 bit message +! 7 = 16 bit message " Rpt" use iso_c_binding, only: c_loc,c_size_t use packjt77 @@ -27,24 +22,15 @@ subroutine genft2(msg0,ichk,msgsent,i4tone,itype) 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) + integer*1 s16(16) real*8 xi(864),xq(864),pi,twopi - data s8/0,1,1,1,0,0,1,0/ + data s16/0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0/ equivalence (ihash,i1hash) - logical first,unpk77_success - data first/.true./ - save + logical unpk77_success - 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 + nsym=128 + pi=4.0*atan(1.0) + twopi=8.*atan(1.0) message(1:37)=' ' itype=1 @@ -89,36 +75,12 @@ subroutine genft2(msg0,ichk,msgsent,i4tone,itype) 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) + bitseq(1:16)=s16 + bitseq(17:144)=codeword 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/getcandidates2.f90 b/lib/fsk4hf/getcandidates2.f90 new file mode 100644 index 000000000..2f61e8623 --- /dev/null +++ b/lib/fsk4hf/getcandidates2.f90 @@ -0,0 +1,51 @@ +subroutine getcandidates2(id,nfa,nfb,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) + 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 + + 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 + ncand=1 + candidate(1,1)=imax*df + +return +end subroutine getcandidates2