From b4145e9b38162e11de785200f201b8fa45dd9685 Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Fri, 13 Jan 2017 15:12:27 +0000 Subject: [PATCH] QRA64 decoding basically alive within MAP65. git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/map65@7487 ab8295b8-cf94-4d9e-aec4-7959e3be5d79 --- libm65/CMakeLists.txt | 20 +++++ libm65/averms.f90 | 20 +++++ libm65/badmsg.f90 | 46 +++++++++++ libm65/fchisq0.f90 | 23 ++++++ libm65/lorentzian.f90 | 102 ++++++++++++++++++++++++ libm65/map65a.f90 | 6 +- libm65/pctile2.f90 | 22 +++++ libm65/qra64b.f90 | 14 ++-- libm65/qra64c.f90 | 181 ++++++++++++++++++++++++++++++++++++++++++ libm65/shell.f90 | 27 +++++++ libm65/smo.f90 | 19 +++++ libm65/spec64.f90 | 42 ++++++++++ libm65/sync64.f90 | 154 +++++++++++++++++++++++++++++++++++ libm65/twkfreq2.f90 | 26 ++++++ 14 files changed, 694 insertions(+), 8 deletions(-) create mode 100644 libm65/averms.f90 create mode 100644 libm65/badmsg.f90 create mode 100644 libm65/fchisq0.f90 create mode 100644 libm65/lorentzian.f90 create mode 100644 libm65/pctile2.f90 create mode 100644 libm65/qra64c.f90 create mode 100644 libm65/shell.f90 create mode 100644 libm65/smo.f90 create mode 100644 libm65/spec64.f90 create mode 100644 libm65/sync64.f90 create mode 100644 libm65/twkfreq2.f90 diff --git a/libm65/CMakeLists.txt b/libm65/CMakeLists.txt index 124c2fc42..ebb0fb5b4 100644 --- a/libm65/CMakeLists.txt +++ b/libm65/CMakeLists.txt @@ -60,6 +60,8 @@ set (FSRCS astro.f90 astro0.f90 astrosub.f90 + averms.f90 + badmsg.f90 ccf2.f90 ccf65.f90 cgen65.f90 @@ -79,6 +81,7 @@ set (FSRCS encode65.f90 extract.f90 fchisq.f90 + fchisq0.f90 fil6521.f90 filbig.f90 fmtmsg.f90 @@ -100,6 +103,7 @@ set (FSRCS iqfix.f90 jt65code.f90 k2grid.f90 + lorentzian.f90 map65a.f90 moon2.f90 moondop.f90 @@ -107,24 +111,31 @@ set (FSRCS noisegen.f90 packjt.f90 pctile.f90 + pctile2.f90 pfxdump.f90 qra64b.f90 + qra64c.f90 recvpkt.f90 rfile3a.f90 s3avg.f90 sec_midn.f90 set.f90 setup65.f90 + shell.f90 sleep_msec.f90 + smo.f90 sort.f90 + spec64.f90 sun.f90 symspec.f90 + sync64.f90 timer.f90 timf2.f90 tm2.f90 toxyz.f90 trimlist.f90 twkfreq.f90 + twkfreq2.f90 zplot.f90 f77_wisdom.f @@ -141,6 +152,15 @@ set (CSRCS tmoonsub.c usleep.c wrapkarn.c + + ../../wsjtx/lib/qra/qra64/qra64.c + ../../wsjtx/lib/qra/qra64/qra64_subs.c + ../../wsjtx/lib/qra/qracodes/npfwht.c + ../../wsjtx/lib/qra/qracodes/pdmath.c + ../../wsjtx/lib/qra/qracodes/qra12_63_64_irr_b.c + ../../wsjtx/lib/qra/qracodes/qra13_64_64_irr_e.c + ../../wsjtx/lib/qra/qracodes/qracodes.c + ../../wsjtx/lib/qra/qracodes/normrnd.c ) if (WIN32) diff --git a/libm65/averms.f90 b/libm65/averms.f90 new file mode 100644 index 000000000..d4f41846c --- /dev/null +++ b/libm65/averms.f90 @@ -0,0 +1,20 @@ +subroutine averms(x,n,nskip,ave,rms) + real x(n) + integer ipk(1) + + ns=0 + s=0. + sq=0. + ipk=maxloc(x) + do i=1,n + if(abs(i-ipk(1)).gt.nskip) then + s=s + x(i) + sq=sq + x(i)**2 + ns=ns+1 + endif + enddo + ave=s/ns + rms=sqrt(sq/ns - ave*ave) + + return +end subroutine averms diff --git a/libm65/badmsg.f90 b/libm65/badmsg.f90 new file mode 100644 index 000000000..1ebc12510 --- /dev/null +++ b/libm65/badmsg.f90 @@ -0,0 +1,46 @@ +subroutine badmsg(irc,dat,nc1,nc2,ng2) + +! Get rid of a few QRA64 false decodes that cannot be correct messages. + + integer dat(12) !Decoded message (as 12 integers) + + ic1=ishft(dat(1),22) + ishft(dat(2),16) + ishft(dat(3),10)+ & + ishft(dat(4),4) + iand(ishft(dat(5),-2),15) + +! Test for "......" or "CQ 000" + if(ic1.eq.262177560 .or. ic1.eq.262177563) then + irc=-1 + return + endif + + ic2=ishft(iand(dat(5),3),26) + ishft(dat(6),20) + & + ishft(dat(7),14) + ishft(dat(8),8) + ishft(dat(9),2) + & + iand(ishft(dat(10),-4),3) + + ig=ishft(iand(dat(10),15),12) + ishft(dat(11),6) + dat(12) + +! Test for blank, -01 to -30, R-01 to R-30, RO, RRR, 73 + if(ig.ge.32401 .and. ig.le.32464) return + + if(ig.ge.14220 .and. ig.le.14229) return !-41 to -50 + if(ig.ge.14040 .and. ig.le.14049) return !-31 to -40 + + if(ig.ge.13320 .and. ig.le.13329) return !+00 to +09 + if(ig.ge.13140 .and. ig.le.13149) return !+10 to +19 + if(ig.ge.12960 .and. ig.le.12969) return !+20 to +29 + if(ig.ge.12780 .and. ig.le.12789) return !+30 to +39 + if(ig.ge.12600 .and. ig.le.12609) return !+40 to +49 + + if(ig.ge.12420 .and. ig.le.12429) return !R-41 to R-50 + if(ig.ge.12240 .and. ig.le.12249) return !R-31 to R-40 + + if(ig.ge.11520 .and. ig.le.11529) return !R+00 to R+09 + if(ig.ge.11340 .and. ig.le.11349) return !R+10 to R+19 + if(ig.ge.11160 .and. ig.le.11169) return !R+20 to R+29 + if(ig.ge.10980 .and. ig.le.10989) return !R+30 to R+39 + if(ig.ge.10800 .and. ig.le.10809) return !R+40 to R+49 + + if(ic1.eq.nc1 .and. ic2.eq.nc2 .and. ng2.ne.32401 .and. ig.ne.ng2) irc=-1 + + return +end subroutine badmsg diff --git a/libm65/fchisq0.f90 b/libm65/fchisq0.f90 new file mode 100644 index 000000000..2c62f2757 --- /dev/null +++ b/libm65/fchisq0.f90 @@ -0,0 +1,23 @@ +real function fchisq0(y,npts,a) + + real y(npts),a(4) + +! rewind 51 + chisq = 0. + do i=1,npts + x=i + z=(x-a(3))/(0.5*a(4)) + yfit=a(1) + if(abs(z).lt.3.0) then + d=1.0 + z*z + yfit=a(1) + a(2) * (1.0/d - 0.1) + endif + chisq=chisq + (y(i) - yfit)**2 +! write(51,3001) i,y(i),yfit,y(i)-yfit +!3001 format(i5,3f10.4) + enddo + fchisq0=chisq + + return +end function fchisq0 + diff --git a/libm65/lorentzian.f90 b/libm65/lorentzian.f90 new file mode 100644 index 000000000..cd2257a75 --- /dev/null +++ b/libm65/lorentzian.f90 @@ -0,0 +1,102 @@ +subroutine lorentzian(y,npts,a) + +! Input: y(npts); assume x(i)=i, i=1,npts +! Output: a(5) +! a(1) = baseline +! a(2) = amplitude +! a(3) = x0 +! a(4) = width +! a(5) = chisqr + + real y(npts) + real a(5) + real deltaa(4) + + a=0. + df=12000.0/8192.0 !df = 1.465 Hz + width=0. + ipk=0 + ymax=-1.e30 + do i=1,npts + if(y(i).gt.ymax) then + ymax=y(i) + ipk=i + endif +! write(50,3001) i,i*df,y(i) +!3001 format(i6,2f12.3) + enddo +! base=(sum(y(ipk-149:ipk-50)) + sum(y(ipk+51:ipk+150)))/200.0 + base=(sum(y(1:20)) + sum(y(npts-19:npts)))/40.0 + stest=ymax - 0.5*(ymax-base) + ssum=y(ipk) + do i=1,50 + if(ipk+i.gt.npts) exit + if(y(ipk+i).lt.stest) exit + ssum=ssum + y(ipk+i) + enddo + do i=1,50 + if(ipk-i.lt.1) exit + if(y(ipk-i).lt.stest) exit + ssum=ssum + y(ipk-i) + enddo + ww=ssum/y(ipk) + width=2 + t=ww*ww - 5.67 + if(t.gt.0.0) width=sqrt(t) + a(1)=base + a(2)=ymax-base + a(3)=ipk + a(4)=width + +! Now find Lorentzian parameters + + deltaa(1)=0.1 + deltaa(2)=0.1 + deltaa(3)=1.0 + deltaa(4)=1.0 + nterms=4 + +! Start the iteration + chisqr=0. + chisqr0=1.e6 + do iter=1,5 + do j=1,nterms + chisq1=fchisq0(y,npts,a) + fn=0. + delta=deltaa(j) +10 a(j)=a(j)+delta + chisq2=fchisq0(y,npts,a) + if(chisq2.eq.chisq1) go to 10 + if(chisq2.gt.chisq1) then + delta=-delta !Reverse direction + a(j)=a(j)+delta + tmp=chisq1 + chisq1=chisq2 + chisq2=tmp + endif +20 fn=fn+1.0 + a(j)=a(j)+delta + chisq3=fchisq0(y,npts,a) + if(chisq3.lt.chisq2) then + chisq1=chisq2 + chisq2=chisq3 + go to 20 + endif + +! Find minimum of parabola defined by last three points + delta=delta*(1./(1.+(chisq1-chisq2)/(chisq3-chisq2))+0.5) + a(j)=a(j)-delta + deltaa(j)=deltaa(j)*fn/3. +! write(*,4000) iter,j,a,chisq2 +!4000 format(i1,i2,4f10.4,f11.3) + enddo + chisqr=fchisq0(y,npts,a) +! write(*,4000) 0,0,a,chisqr + if(chisqr/chisqr0.gt.0.99) exit + chisqr0=chisqr + enddo + a(5)=chisqr + + return +end subroutine lorentzian + diff --git a/libm65/map65a.f90 b/libm65/map65a.f90 index 06488bd0a..ba74837bc 100644 --- a/libm65/map65a.f90 +++ b/libm65/map65a.f90 @@ -33,7 +33,7 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, & data nfile/0/,nutc0/-999/,nid/0/,ip000/1/,ip001/1/,mousefqso0/-999/ save - bqra64=nfast.ge.100 + bqra64=nfast.ge.100 nfast=mod(nfast,100) mcall3a=mcall3b mousefqso0=mousefqso @@ -218,7 +218,8 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, & a,dt,pol,nkv,nhist,nsum,nsave,qual,decoded) call timer('decode1a',1) if(nqd.eq.2) then - call qra64b(nutc,ikhz) + call qra64b(nutc,nqd,ikhz,mousedf,ntol,xpol,mycall, & + hiscall,hisgrid) cycle endif @@ -348,6 +349,7 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, & 1013 format('',2i4) flush(6) endif + if(nqd.eq.2) exit !### TESTING: do only QRA64 if(nagain.eq.1 .and. nqd.eq.1) go to 999 enddo diff --git a/libm65/pctile2.f90 b/libm65/pctile2.f90 new file mode 100644 index 000000000..07856fdef --- /dev/null +++ b/libm65/pctile2.f90 @@ -0,0 +1,22 @@ +subroutine pctile2(x,npts,npct,xpct) + + parameter (NMAX=100000) + real*4 x(npts) + real*4 tmp(NMAX) + + if(npts.le.0) then + xpct=1.0 + go to 900 + endif + if(npts.gt.NMAX) stop + + tmp(1:npts)=x + call shell(npts,tmp) + j=nint(npts*0.01*npct) + if(j.lt.1) j=1 + if(j.gt.npts) j=npts + xpct=tmp(j) + +900 continue + return +end subroutine pctile2 diff --git a/libm65/qra64b.f90 b/libm65/qra64b.f90 index 9723f2cdd..5f6d22836 100644 --- a/libm65/qra64b.f90 +++ b/libm65/qra64b.f90 @@ -1,14 +1,14 @@ -subroutine qra64b(nutc,ikhz) +subroutine qra64b(nutc,nqd,ikhz,idf,ntol,xpol,mycall_12,hiscall_12, & + hisgrid_6) parameter (NFFT1=5376000) !56*96000 parameter (NFFT2=336000) !56*6000 (downsampled by 1/16) complex cx(0:NFFT2-1),cy(0:NFFT2-1) complex ca(NFFT1),cb(NFFT1) !FFTs of raw x,y data + logical xpol + character*12 mycall_12,hiscall_12 + character*6 hisgrid_6 common/cacb/ca,cb - -!### - if(nutc.ne.-999) return -!### df=96000.0/NFFT1 k0=(ikhz-75.7)*1000.0/df @@ -25,7 +25,9 @@ subroutine qra64b(nutc,ikhz) call four2a(cx,NFFT2,1,-1,1) call four2a(cy,NFFT2,1,-1,1) - write(67) nutc,cx,cy +! write(67) nutc,cx,cy + call qra64c(cx,cy,nutc,nqd,ikhz,idf,ntol,xplo,mycall_12, & + hiscall_12,hisgrid_6) return end subroutine qra64b diff --git a/libm65/qra64c.f90 b/libm65/qra64c.f90 new file mode 100644 index 000000000..f0f089f5a --- /dev/null +++ b/libm65/qra64c.f90 @@ -0,0 +1,181 @@ +subroutine qra64c(cx,cy,nutc,nqd,ikhz,nfqso,ntol,xpol,mycall_12, & + hiscall_12,hisgrid_6) + + use packjt + parameter (NFFT2=336000) !56*6000 (downsampled by 1/16) + parameter (NMAX=60*12000,LN=1152*63) + +! Required input data: +! nutc,cx,cy,nf1,nf2,nfqso,ntol,mode64,minsync,ndepth,emedelay, +! mycall_12,hiscall_12,hisgrid_6 + + character decoded*22 + character*12 mycall_12,hiscall_12 + character*6 mycall,hiscall,hisgrid_6 + character*4 hisgrid + character*1 cp + logical xpol,ltext + complex cx(0:NFFT2-1),cy(0:NFFT2-1) + complex c00(0:720000) !Complex spectrum of dd() + complex c0(0:720000) !Complex data for dd() + real a(3) + real s3(LN) !Symbol spectra + real s3a(LN) !Symbol spectra + integer dat4(12) !Decoded message (as 12 integers) + integer dat4x(12) + integer nap(0:11) + data nap/0,2,3,2,3,4,2,3,6,4,6,6/ + data nc1z/-1/,nc2z/-1/,ng2z/-1/,maxaptypez/-1/ + save + +! For now: + nf1=200 + nf2=2200 +! nfqso=808 +! ntol=1000 + mode64=1 + minsync=-1 + ndepth=3 + emedelay=2.5 + + irc=-1 + decoded=' ' + nft=99 + mycall=mycall_12(1:6) + hiscall=hiscall_12(1:6) + hisgrid=hisgrid_6(1:4) + call packcall(mycall,nc1,ltext) + call packcall(hiscall,nc2,ltext) + call packgrid(hisgrid,ng2,ltext) + nSubmode=0 + if(mode64.eq.2) nSubmode=1 + if(mode64.eq.4) nSubmode=2 + if(mode64.eq.8) nSubmode=3 + if(mode64.eq.16) nSubmode=4 + b90=1.0 + nFadingModel=1 + maxaptype=4 + if(iand(ndepth,64).ne.0) maxaptype=5 + if(nc1.ne.nc1z .or. nc2.ne.nc2z .or. ng2.ne.ng2z .or. & + maxaptype.ne.maxaptypez) then + do naptype=0,maxaptype + if(naptype.eq.2 .and. maxaptype.eq.4) cycle + call qra64_dec(s3,nc1,nc2,ng2,naptype,1,nSubmode,b90, & + nFadingModel,dat4,snr2,irc) + enddo + nc1z=nc1 + nc2z=nc2 + ng2z=ng2 + maxaptypez=maxaptype + endif + naptype=maxaptype + npts2=NFFT2 + +!1 read(67,end=999) nutc,cx,cy + c00(0:NFFT2-1)=conjg(cy) + call sync64(c00,nf1,nf2,nfqso,ntol,mode64,emedelay,dtx,f0,jpk0,sync, & + sync2,width) + + nfreq=nint(f0) + if(mode64.eq.1 .and. minsync.ge.0 .and. (sync-7.0).lt.minsync) go to 900 + a=0. + a(1)=-f0 + call twkfreq2(c00,c0,npts2,6000.0,a) + + irc=-99 + s3lim=20. + itz=11 + if(mode64.eq.4) itz=9 + if(mode64.eq.2) itz=7 + if(mode64.eq.1) itz=5 + + LL=64*(mode64+2) + NN=63 + napmin=99 + do itry0=1,5 + idt=itry0/2 + if(mod(itry0,2).eq.0) idt=-idt + jpk=jpk0 + 750*idt + call spec64(c0,npts2,mode64,jpk,s3a,LL,NN) + call pctile2(s3a,LL*NN,40,base) + s3a=s3a/base + where(s3a(1:LL*NN)>s3lim) s3a(1:LL*NN)=s3lim + do iter=itz,0,-2 + b90=1.728**iter + if(b90.gt.230.0) cycle + if(b90.lt.0.15*width) exit + s3(1:LL*NN)=s3a(1:LL*NN) + call qra64_dec(s3,nc1,nc2,ng2,naptype,0,nSubmode,b90, & + nFadingModel,dat4,snr2,irc) + if(irc.eq.0) go to 10 + if(irc.gt.0) call badmsg(irc,dat4,nc1,nc2,ng2) + iirc=max(0,min(irc,11)) + if(irc.gt.0 .and. nap(iirc).lt.napmin) then + dat4x=dat4 + b90x=b90 + snr2x=snr2 + napmin=nap(iirc) + irckeep=irc + dtxkeep=jpk/6000.0 - 1.0 + itry0keep=itry0 + iterkeep=iter + endif + enddo + if(irc.eq.0) exit + enddo + + if(napmin.ne.99) then + dat4=dat4x + b90=b90x + snr2=snr2x + irc=irckeep + dtx=dtxkeep + itry0=itry0keep + iter=iterkeep + endif +10 decoded=' ' + + if(irc.ge.0) then + call unpackmsg(dat4,decoded) !Unpack the user message + call fmtmsg(decoded,iz) + if(index(decoded,"000AAA ").ge.1) then + ! Suppress a certain type of garbage decode. + decoded=' ' + irc=-1 + endif + nft=100 + irc + nsnr=nint(snr2) + else + snr2=0. + endif + +900 if(irc.lt.0) then + sy=max(1.0,sync) + if(nSubmode.eq.0) nsnr=nint(10.0*log10(sy)-35.0) !A + if(nSubmode.eq.1) nsnr=nint(10.0*log10(sy)-34.0) !B + if(nSubmode.eq.2) nsnr=nint(10.0*log10(sy)-29.0) !C + if(nSubmode.eq.3) nsnr=nint(10.0*log10(sy)-29.0) !D + if(nSubmode.eq.4) nsnr=nint(10.0*log10(sy)-24.0) !E + endif + +! write(*,1011) nutc/100,nsnr,dtx,nfreq,decoded +!1011 format(i4.4,i4,f5.1,i5,1x,2x,1x,a22) + + nkhz=108 + npol=0 + cp='H' + nsync=sync + ntxpol=0 + if(irc.ge.0) then + write(*,1010) nkHz,nfreq,npol,nutc/100,dtx,nsnr,decoded,irc,ntxpol,cp +!1010 format('!',i3,i5,i4,i7.6,f5.1,i4,2x,a22,i2,i5,i5,1x,a1) +!1010 format(i3,i5,i4,i5.4,f5.1,i5,2x,a22,i2,i5,1x,a1) +1010 format('!',i3,i5,i4,i6.4,f5.1,i5,2x,a22,i2,i5,1x,a1) + else + write(*,1010) nkHz,nfreq,npol,nutc/100,dtx,nsync + endif + +! goto 1 + +999 return +end subroutine qra64c diff --git a/libm65/shell.f90 b/libm65/shell.f90 new file mode 100644 index 000000000..1e9b44889 --- /dev/null +++ b/libm65/shell.f90 @@ -0,0 +1,27 @@ +subroutine shell(n,a) + integer n + real a(n) + integer i,j,inc + real v + + inc=1 +1 inc=3*inc+1 + if(inc.le.n) go to 1 +2 inc=inc/3 + + do i=inc+1,n + v=a(i) + j=i +3 if(a(j-inc).gt.v) then + a(j)=a(j-inc) + j=j-inc + if(j.le.inc) go to 4 + go to 3 + endif +4 a(j)=v + enddo + + if(inc.gt.1) go to 2 + + return +end subroutine shell diff --git a/libm65/smo.f90 b/libm65/smo.f90 new file mode 100644 index 000000000..c42de7c60 --- /dev/null +++ b/libm65/smo.f90 @@ -0,0 +1,19 @@ +subroutine smo(x,npts,y,nadd) + + real x(npts) + real y(npts) + + nh=nadd/2 + do i=1+nh,npts-nh + sum=0. + do j=-nh,nh + sum=sum + x(i+j) + enddo + y(i)=sum + enddo + x=y + x(:nh)=0. + x(npts-nh+1:)=0. + + return +end subroutine smo diff --git a/libm65/spec64.f90 b/libm65/spec64.f90 new file mode 100644 index 000000000..fc303ae02 --- /dev/null +++ b/libm65/spec64.f90 @@ -0,0 +1,42 @@ +subroutine spec64(c0,npts2,mode64,jpk,s3,LL,NN) + + parameter (NSPS=3456) !Samples per symbol at 6000 Hz + complex c0(0:360000) !Complex spectrum of dd() + complex cs(0:NSPS-1) !Complex symbol spectrum + real s3(LL,NN) !Synchronized symbol spectra + real xbase0(LL),xbase(LL) + + nfft=nsps + fac=1.0/nfft + do j=1,NN + jj=j+7 !Skip first Costas array + if(j.ge.33) jj=j+14 !Skip middle Costas array + ja=jpk + (jj-1)*nfft + jb=ja+nfft-1 + cs(0:nfft-1)=fac*c0(ja:jb) + call four2a(cs,nfft,1,-1,1) + do ii=1,LL + i=ii-65 + if(i.lt.0) i=i+nfft + s3(ii,j)=real(cs(i))**2 + aimag(cs(i))**2 + enddo + enddo + + df=6000.0/nfft + do i=1,LL + call pctile2(s3(i,1:NN),NN,45,xbase0(i)) !Get baseline for passband shape + enddo + + nh=25 + xbase(1:nh-1)=sum(xbase0(1:nh-1))/(nh-1.0) + xbase(LL-nh+1:LL)=sum(xbase0(LL-nh+1:LL))/(nh-1.0) + do i=nh,LL-nh + xbase(i)=sum(xbase0(i-nh+1:i+nh))/(2*nh+1) !Smoothed passband shape + enddo + + do i=1,LL + s3(i,1:NN)=s3(i,1:NN)/(xbase(i)+0.001) !Apply frequency equalization + enddo + + return +end subroutine spec64 diff --git a/libm65/sync64.f90 b/libm65/sync64.f90 new file mode 100644 index 000000000..d52ad23d2 --- /dev/null +++ b/libm65/sync64.f90 @@ -0,0 +1,154 @@ +subroutine sync64(c0,nf1,nf2,nfqso,ntol,mode64,emedelay,dtx,f0,jpk,sync, & + sync2,width) + +! use timer_module, only: timer + + parameter (NMAX=60*12000) !Max size of raw data at 12000 Hz + parameter (NSPS=3456) !Samples per symbol at 6000 Hz + parameter (NSPC=7*NSPS) !Samples per Costas array + real s1(0:NSPC-1) !Power spectrum of Costas 1 + real s2(0:NSPC-1) !Power spectrum of Costas 2 + real s3(0:NSPC-1) !Power spectrum of Costas 3 + real s0(0:NSPC-1) !Sum of s1+s2+s3 + real s0a(0:NSPC-1) !Best synchromized spectrum (saved) + real s0b(0:NSPC-1) !tmp + real a(5) + integer icos7(0:6) !Costas 7x7 tones + integer ipk0(1) + complex cc(0:NSPC-1) !Costas waveform + complex c0(0:720000) !Complex spectrum of dd() + complex c1(0:NSPC-1) !Complex spectrum of Costas 1 + complex c2(0:NSPC-1) !Complex spectrum of Costas 2 + complex c3(0:NSPC-1) !Complex spectrum of Costas 3 + data icos7/2,5,6,0,4,1,3/ !Costas 7x7 tone pattern + data mode64z/-1/ + save + + if(mode64.ne.mode64z) then + twopi=8.0*atan(1.0) + dfgen=mode64*12000.0/6912.0 + k=-1 + phi=0. + do j=0,6 !Compute complex Costas waveform + dphi=twopi*10.0*icos7(j)*dfgen/6000.0 + do i=1,NSPS + phi=phi + dphi + if(phi.gt.twopi) phi=phi-twopi + k=k+1 + cc(k)=cmplx(cos(phi),sin(phi)) + enddo + enddo + mode64z=mode64 + endif + + nfft3=NSPC + nh3=nfft3/2 + df3=6000.0/nfft3 + + fa=max(nf1,nfqso-ntol) + fb=min(nf2,nfqso+ntol) + iaa=max(0,nint(fa/df3)) + ibb=min(NSPC-1,nint(fb/df3)) + + maxtol=max(ntol,500) + fa=max(nf1,nfqso-maxtol) + fb=min(nf2,nfqso+maxtol) + ia=max(0,nint(fa/df3)) + ib=min(NSPC-1,nint(fb/df3)) + id=0.1*(ib-ia) + iz=ib-ia+1 + sync=-1.e30 + smaxall=0. + jpk=0 + ja=0 + jb=(5.0+emedelay)*6000 + jstep=100 + ipk=0 + kpk=0 + nadd=10*mode64 + if(mod(nadd,2).eq.0) nadd=nadd+1 !Make nadd odd + nskip=max(49,nadd) + + do j1=ja,jb,jstep + call timer('sync64_1',0) + j2=j1 + 39*NSPS + j3=j1 + 77*NSPS + c1=1.e-4*c0(j1:j1+NSPC-1) * conjg(cc) + c2=1.e-4*c0(j2:j2+NSPC-1) * conjg(cc) + c3=1.e-4*c0(j3:j3+NSPC-1) * conjg(cc) + call four2a(c1,nfft3,1,-1,1) + call four2a(c2,nfft3,1,-1,1) + call four2a(c3,nfft3,1,-1,1) + s1=0. + s2=0. + s3=0. + s0b=0. + do i=ia,ib + freq=i*df3 + s1(i)=real(c1(i))**2 + aimag(c1(i))**2 + s2(i)=real(c2(i))**2 + aimag(c2(i))**2 + s3(i)=real(c3(i))**2 + aimag(c3(i))**2 + enddo + call timer('sync64_1',1) + + call timer('sync64_2',0) + s0(ia:ib)=s1(ia:ib) + s2(ia:ib) + s3(ia:ib) + s0(:ia-1)=0. + s0(ib+1:)=0. + if(nadd.ge.3) then + do ii=1,3 + s0b(ia:ib)=s0(ia:ib) + call smo(s0b(ia:ib),iz,s0(ia:ib),nadd) + enddo + endif + call averms(s0(ia+id:ib-id),iz-2*id,nskip,ave,rms) + s=(maxval(s0(iaa:ibb))-ave)/rms + ipk0=maxloc(s0(iaa:ibb)) + ip=ipk0(1) + iaa - 1 + if(s.gt.sync) then + jpk=j1 + s0a=(s0-ave)/rms + sync=s + dtx=jpk/6000.0 - 1.0 + ipk=ip + f0=ip*df3 + endif + call timer('sync64_2',1) + enddo + + s0a=s0a+2.0 + write(17) ia,ib,s0a(ia:ib) !Save data for red curve + close(17) + + nskip=50 + call lorentzian(s0a(ia+nskip:ib-nskip),iz-2*nskip,a) + f0a=(a(3)+ia+49)*df3 + w1=df3*a(4) + w2=2*nadd*df3 + width=w1 + if(w1.gt.1.2*w2) width=sqrt(w1**2 - w2**2) + + sq=0. + do i=1,20 + j=ia+nskip+1 + k=ib-nskip-21+i + sq=sq + (s0a(j)-a(1))**2 + (s0a(k)-a(1))**2 + enddo + rms2=sqrt(sq/40.0) + sync2=10.0*log10(a(2)/rms2) + +! do i=1,iz-2*nskip +! x=i +! z=(x-a(3))/(0.5*a(4)) +! yfit=a(1) +! if(abs(z).lt.3.0) then +! d=1.0 + z*z +! yfit=a(1) + a(2)*(1.0/d - 0.1) +! endif +! j=i+ia+49 +! write(76,1110) j*df3,s0a(j),yfit +!1110 format(3f10.3) +! enddo + + return +end subroutine sync64 diff --git a/libm65/twkfreq2.f90 b/libm65/twkfreq2.f90 new file mode 100644 index 000000000..2333d9fb2 --- /dev/null +++ b/libm65/twkfreq2.f90 @@ -0,0 +1,26 @@ +subroutine twkfreq2(c3,c4,npts,fsample,a) + + complex c3(npts) + complex c4(npts) + complex w,wstep + real a(3) + data twopi/6.283185307/ + +! Mix the complex signal + w=1.0 + wstep=1.0 + x0=0.5*(npts+1) + s=2.0/npts + do i=1,npts + x=s*(i-x0) + p2=1.5*x*x - 0.5 +! p3=2.5*(x**3) - 1.5*x +! p4=4.375*(x**4) - 3.75*(x**2) + 0.375 + dphi=(a(1) + x*a(2) + p2*a(3)) * (twopi/fsample) + wstep=cmplx(cos(dphi),sin(dphi)) + w=w*wstep + c4(i)=w*c3(i) + enddo + + return +end subroutine twkfreq2