From c5502cda05fcb37c67ad9886725f182087447021 Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Thu, 8 Oct 2020 16:48:11 -0400 Subject: [PATCH] QRA65 now decodes using qra_loops() -- the same inner loops as QRA64. Very effective! --- CMakeLists.txt | 5 +- lib/ana64.f90 | 5 +- lib/qra64a.f90 | 2 - lib/qra65_decode.f90 | 139 +++++++++++++++---------------------------- lib/qra_loops.f90 | 31 +++------- lib/spec64.f90 | 67 ++++++++++++++++----- 6 files changed, 112 insertions(+), 137 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 4d1246489..8ff8f0ebf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -570,7 +570,6 @@ set (wsjt_FSRCS lib/softsym9w.f90 lib/shell.f90 lib/spec64.f90 - lib/spec_qra65.f90 lib/spec9f.f90 lib/stdmsg.f90 lib/subtract65.f90 @@ -1317,7 +1316,7 @@ if (${OPENMP_FOUND} OR APPLE) endif (APPLE) if (WIN32) set_target_properties (jt9 PROPERTIES - LINK_FLAGS -Wl,--stack,16777216 + LINK_FLAGS -Wl,--stack,28000000 ) endif () target_link_libraries (jt9 wsjt_fort_omp wsjt_cxx fort_qt) @@ -1454,7 +1453,7 @@ else () ) if (WIN32) set_target_properties (wsjtx PROPERTIES - LINK_FLAGS -Wl,--stack,16777216 + LINK_FLAGS -Wl,--stack,28000000 ) endif () endif () diff --git a/lib/ana64.f90 b/lib/ana64.f90 index 5681ee176..faebb84e5 100644 --- a/lib/ana64.f90 +++ b/lib/ana64.f90 @@ -3,13 +3,12 @@ subroutine ana64(dd,npts,c0) 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 dd(NMAX) !Raw data complex c0(0:720000) !Complex spectrum of dd() save - nfft1=672000 +! nfft1=672000 + nfft1=720000 nfft2=nfft1/2 df1=12000.0/nfft1 fac=2.0/nfft1 diff --git a/lib/qra64a.f90 b/lib/qra64a.f90 index 953a4789b..e5da55afd 100644 --- a/lib/qra64a.f90 +++ b/lib/qra64a.f90 @@ -81,8 +81,6 @@ subroutine qra64a(dd,npts,nf1,nf2,nfqso,ntol,mode64,minsync,ndepth, & snr2=0. endif nfreq=nint(f0) - write(71,3071) idf,idt,ncall,irc,nsnr,dtx,nfreq,decoded -3071 format(5i5,f7.2,i6,2x,a22) 900 if(irc.lt.0) then sy=max(1.0,sync) diff --git a/lib/qra65_decode.f90 b/lib/qra65_decode.f90 index 0f04565f7..d377a6d11 100644 --- a/lib/qra65_decode.f90 +++ b/lib/qra65_decode.f90 @@ -41,18 +41,20 @@ contains character(len=6) :: hisgrid character*37 decoded integer*2 iwave(NMAX) !Raw data + real dd(NMAX) !Raw data integer dat4(12) logical lapdx,ltext + complex, allocatable :: c00(:) !Analytic signal, 6000 S/s complex, allocatable :: c0(:) !Analytic signal, 6000 S/s real, allocatable, save :: s3(:,:) !Symbol spectra real, allocatable, save :: s3a(:,:) !Symbol spectra for avg messages - real a(5) data nc1z/-1/,nc2z/-1/,ng2z/-1/,maxaptypez/-1/,nsubmodez/-1/ save nc1z,nc2z,ng2z,maxaptypez,nsave,nsubmodez mode65=2**nsubmode nfft1=ntrperiod*12000 nfft2=ntrperiod*6000 + allocate (c00(0:nfft1-1)) allocate (c0(0:nfft1-1)) if(nsubmode.ne.nsubmodez) then @@ -74,106 +76,61 @@ contains nsps=41472 else stop 'Invalid TR period' - endif - baud=12000.0/nsps - df1=12000.0/nfft1 - -! do i=1,12000*ntrperiod -! write(61,3061) i/12000.0,iwave(i)/32767.0 -!3061 format(2f12.6) -! enddo - - this%callback => callback - - if(nutc.eq.-999) print*,lapdx,nfa,nfb,nfqso !Silence warning + endif + npts=ntrperiod*12000 + baud=12000.0/nsps + df1=12000.0/nfft1 + this%callback => callback + if(nutc.eq.-999) print*,lapdx,nfa,nfb,nfqso !Silence warning ! Prime the QRA decoder for possible use of AP - call packcall(mycall(1:6),nc1,ltext) - call packcall(hiscall(1:6),nc2,ltext) - call packgrid(hisgrid(1:4),ng2,ltext) - b90=20.0 !8 to 25 is OK; not very critical - nFadingModel=1 + call packcall(mycall(1:6),nc1,ltext) + call packcall(hiscall(1:6),nc2,ltext) + call packgrid(hisgrid(1:4),ng2,ltext) + b90=20.0 !8 to 25 is OK; not very critical + nFadingModel=1 ! AP control could be done differently, but this works well: - maxaptype=0 - if(ndepth.eq.2) maxaptype=3 - if(ndepth.eq.3) maxaptype=5 + maxaptype=0 + if(ndepth.eq.2) maxaptype=3 + if(ndepth.eq.3) 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 - s3a=0. - nsave=0 - endif - naptype=maxaptype + 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 + s3a=0. + nsave=0 + endif + naptype=maxaptype - call timer('sync_q65',0) - call sync_qra65(iwave,ntrperiod*12000,mode65,nsps,nfqso,ntol,xdt,f0,snr1) - call timer('sync_q65',1) + call timer('sync_q65',0) + call sync_qra65(iwave,ntrperiod*12000,mode65,nsps,nfqso,ntol,xdt,f0,snr1) + call timer('sync_q65',1) -! Downsample to give complex data at 6000 S/s - call timer('down_q65',0) - fac=2.0/nfft1 - c0=fac*iwave(1:nfft1) - call four2a(c0,nfft1,1,-1,1) !Forward c2c FFT - c0(nfft2/2+1:nfft2)=0. !Zero the top half - c0(0)=0.5*c0(0) - call four2a(c0,nfft2,1,1,1) !Inverse c2c FFT - a=0. - a(1)=-(f0 + mode65*baud) !Data tones start mode65 bins higher - call twkfreq(c0,c0,ntrperiod*6000,6000.0,a) - call timer('down_q65',1) + jpk0=(xdt+1.0)*6000 !### + if(jpk0.lt.0) jpk0=0 - jpk=(xdt+0.5)*6000 - 384 !### Empirical ### - if(ntrperiod.ge.60) jpk=(xdt+1.0)*6000 - 384 !### TBD ??? ### - if(jpk.lt.0) jpk=0 - xdt=jpk/6000.0 - 0.5 - LL=64*(mode65+2) - NN=63 - call timer('spec_q65',0) - call spec_qra65(c0(jpk:),nsps/2,s3,LL,NN) !Compute synced symbol spectra - call timer('spec_q65',1) + fac=1.0/32767.0 + dd=fac*iwave +! npts=648000 + minsync=-2 + nmode=65 - do j=1,63 !Normalize to symbol baseline - call pctile(s3(:,j),LL,40,base) - s3(:,j)=s3(:,j)/base - enddo + call ana64(dd,npts,c00) - LL2=64*(mode65+1)-1 - s3max=20.0 - do j=1,63 !Apply AGC to suppress pings - xx=maxval(s3(-64:LL2,j)) - if(xx.gt.s3max) s3(-64:LL2,j)=s3(-64:LL2,j)*s3max/xx - enddo - -! Call Nico's QRA64 decoder - call timer('qra64_de',0) - call qra64_dec(s3,nc1,nc2,ng2,naptype,0,nSubmode,b90, & - nFadingModel,dat4,snr2,irc) - call timer('qra64_de',1) - - if(irc.lt.0) then -! No luck so far. Try for an average decode. - call timer('qra64_av',0) - s3a=s3a+s3 - nsave=nsave+1 - if(nsave.ge.2) then - call qra64_dec(s3a,nc1,nc2,ng2,naptype,0,nSubmode,b90, & - nFadingModel,dat4,snr2,irc) - if(irc.ge.0) irc=100*nsave + irc - endif - call timer('qra64_av',1) - endif - snr2=snr2 + db(6912.0/nsps) - if(irc.gt.0) call badmsg(irc,dat4,nc1,nc2,ng2) + call timer('qraloops',0) + call qra_loops(c00,npts/2,nmode,mode65,nsubmode,nFadingModel,minsync, & + ndepth,nc1,nc2,ng2,naptype,jpk0,xdt,f0,width,snr2,s3,irc,dat4) + if(nmode.eq.65) xdt=xdt+0.4 !### Empirical WHY ??? ### + call timer('qraloops',1) decoded=' ' if(irc.ge.0) then diff --git a/lib/qra_loops.f90 b/lib/qra_loops.f90 index e82b708dd..0ea261047 100644 --- a/lib/qra_loops.f90 +++ b/lib/qra_loops.f90 @@ -1,5 +1,5 @@ subroutine qra_loops(c00,npts2,mode,mode64,nsubmode,nFadingModel,minsync, & - ndepth,nc1,nc2,ng2,naptype,jpk0,dtx,f0,width,snr2,s3,irc,dat4) + ndepth,nc1,nc2,ng2,naptype,jpk0,xdt,f0,width,snr2,s3,irc,dat4) use timer_module, only: timer parameter (LN=1152*63) @@ -10,7 +10,6 @@ subroutine qra_loops(c00,npts2,mode,mode64,nsubmode,nFadingModel,minsync, & integer dat4(12),dat4x(12) !Decoded message (as 12 integers) integer nap(0:11) !AP return codes data nap/0,2,3,2,3,4,2,3,6,4,6,6/ -! save irc=-99 s3lim=20. @@ -26,6 +25,8 @@ subroutine qra_loops(c00,npts2,mode,mode64,nsubmode,nFadingModel,minsync, & NN=63 napmin=99 ncall=0 + nsps=3456 !QRA64 + if(mode.eq.65) nsps=3840 !QRA65 do idf0=1,11 idf=idf0/2 @@ -36,25 +37,9 @@ subroutine qra_loops(c00,npts2,mode,mode64,nsubmode,nFadingModel,minsync, & do idt0=1,idtmax idt=idt0/2 if(mod(idt0,2).eq.0) idt=-idt - jpk=jpk0 + 750*idt - if(mode.eq.64) then - call spec64(c0,jpk,s3a,LL,NN) - else - if(jpk.lt.0) jpk=0 - call timer('spec_q65',0) - call spec_qra65(c0(jpk:),nsps2,s3,LL,NN) !Get synced symbol spectra - call timer('spec_q65',1) -! do j=1,63 !Normalize to symbol baseline -! call pctile(s3(:,j),LL,40,base) -! s3(:,j)=s3(:,j)/base -! enddo -! LL2=64*(mode65+1)-1 -! s3max=20.0 -! do j=1,63 !Apply AGC to suppress pings -! xx=maxval(s3(-64:LL2,j)) -! if(xx.gt.s3max) s3(-64:LL2,j)=s3(-64:LL2,j)*s3max/xx -! enddo - endif + jpk=jpk0 + 750*idt + if(jpk.lt.0) jpk=0 + call spec64(c0,nsps,mode,jpk,s3a,LL,NN) call pctile(s3a,LL*NN,40,base) s3a=s3a/base where(s3a(1:LL*NN)>s3lim) s3a(1:LL*NN)=s3lim @@ -77,7 +62,7 @@ subroutine qra_loops(c00,npts2,mode,mode64,nsubmode,nFadingModel,minsync, & snr2x=snr2 napmin=nap(iirc) irckeep=irc - dtxkeep=jpk/6000.0 - 1.0 + xdtkeep=jpk/6000.0 - 1.0 f0keep=-a(1) idfkeep=idf idtkeep=idt @@ -94,7 +79,7 @@ subroutine qra_loops(c00,npts2,mode,mode64,nsubmode,nFadingModel,minsync, & b90=b90x snr2=snr2x irc=irckeep - dtx=dtxkeep + xdt=xdtkeep f0=f0keep idt=idtkeep idf=idfkeep diff --git a/lib/spec64.f90 b/lib/spec64.f90 index 3404243c4..377e8dfcf 100644 --- a/lib/spec64.f90 +++ b/lib/spec64.f90 @@ -1,26 +1,50 @@ -subroutine spec64(c0,jpk,s3,LL,NN) +subroutine spec64(c0,nsps,mode,jpk,s3,LL,NN) - parameter (NSPS=3456) !Samples per symbol at 6000 Hz + parameter (MAXFFT=3840) complex c0(0:360000) !Complex spectrum of dd() - complex cs(0:NSPS-1) !Complex symbol spectrum + complex cs(0:MAXFFT-1) !Complex symbol spectrum real s3(LL,NN) !Synchronized symbol spectra real xbase0(LL),xbase(LL) + integer isync(22) !Indices of sync symbols + data isync/1,9,12,13,15,22,23,26,27,33,35,38,46,50,55,60,62,66,69,74,76,85/ 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 + + if(mode.eq.64) then + 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 - enddo + else + j=0 + n=1 + do k=1,84 + if(k.eq.isync(n)) then + n=n+1 + cycle + endif + j=j+1 + ja=(k-1)*nsps + jpk + jb=ja+nsps-1 + cs(0:nfft-1)=fac*c0(ja:jb) + call four2a(cs,nsps,1,-1,1) !c2c FFT to frequency + do ii=1,LL + i=ii-65 + if(i.lt.0) i=i+nsps + s3(ii,j)=real(cs(i))**2 + aimag(cs(i))**2 + enddo + enddo + endif df=6000.0/nfft do i=1,LL @@ -38,5 +62,18 @@ subroutine spec64(c0,jpk,s3,LL,NN) s3(i,1:NN)=s3(i,1:NN)/(xbase(i)+0.001) !Apply frequency equalization enddo +! print*,'a',LL,NN,jpk +! df=6000.0/nfft +! do i=1,LL +! write(71,3071) i,i-65,i*df,(s3(i,j),j=1,4) +!3071 format(2i8,f10.3,4e12.3) +! enddo +! +! do j=1,NN +! write(72,3072) j,maxloc(s3(1:LL,j)),maxloc(s3(1:LL,j))-65 +!3072 format(3i8) +! enddo +! if(nfft.ne.-999) stop + return end subroutine spec64