From b220643237291cfe6c4b07814cfb0b086ddb5e81 Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Mon, 11 Mar 2024 19:02:38 -0400 Subject: [PATCH] Add routines that will be needed to decode SuperFox. (Presently, that's done in sf.exe.) --- CMakeLists.txt | 2 +- lib/decoder.f90 | 11 +- lib/superfox/decode_sf.f90 | 70 ++++++++++++ lib/superfox/ftrsd3.f90 | 186 +++++++++++++++++++++++++++++++ lib/superfox/getpp3.f90 | 22 ++++ lib/superfox/ran1.f90 | 28 +++++ lib/superfox/rs_sf.a | Bin 0 -> 8296 bytes lib/superfox/sfox_ana.f90 | 21 ++++ lib/superfox/sfox_demod.f90 | 33 ++++++ lib/superfox/sfox_prob.f90 | 55 +++++++++ lib/superfox/sfox_sync.f90 | 153 +++++++++++++++++++++++++ lib/superfox/sfox_unpack.f90 | 51 +++++++++ lib/{ => superfox}/sfox_wave.f90 | 0 13 files changed, 628 insertions(+), 4 deletions(-) create mode 100644 lib/superfox/decode_sf.f90 create mode 100644 lib/superfox/ftrsd3.f90 create mode 100644 lib/superfox/getpp3.f90 create mode 100644 lib/superfox/ran1.f90 create mode 100644 lib/superfox/rs_sf.a create mode 100644 lib/superfox/sfox_ana.f90 create mode 100644 lib/superfox/sfox_demod.f90 create mode 100644 lib/superfox/sfox_prob.f90 create mode 100644 lib/superfox/sfox_sync.f90 create mode 100644 lib/superfox/sfox_unpack.f90 rename lib/{ => superfox}/sfox_wave.f90 (100%) diff --git a/CMakeLists.txt b/CMakeLists.txt index 172ce1d6b..dfdb227a8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -521,7 +521,6 @@ set (wsjt_FSRCS lib/sec0.f90 lib/sec_midn.f90 lib/setup65.f90 - lib/sfox_wave.f90 lib/sh65.f90 lib/sh65snr.f90 lib/slasubs.f @@ -587,6 +586,7 @@ set (wsjt_FSRCS lib/fst4/get_crc24.f90 lib/fst4/fst4_baseline.f90 lib/77bit/hash22calc.f90 + lib/superfox/sfox_wave.f90 ) # temporary workaround for a gfortran v7.3 ICE on Fedora 27 64-bit diff --git a/lib/decoder.f90 b/lib/decoder.f90 index 3386d5c41..1a051e7b6 100644 --- a/lib/decoder.f90 +++ b/lib/decoder.f90 @@ -139,10 +139,15 @@ subroutine multimode_decoder(ss,id2,params,nfsample) open(19,file=trim(temp_dir)//'/houndcallers.txt',status='unknown') endif - if(ncontest.eq.7 .and. params%b_superfox .and. params%b_even_seq) then + if(((ncontest.eq.7 .and. params%b_superfox .and. params%b_even_seq) & + .or. params%ndiskdat) .and. params%nzhsym.ge.50) then ! Call the superFox decoder - print*,'Calling SuperFox decoder',params%nzhsym,params%b_superfox, & - params%b_even_seq + open(47,file='fort.47',status='unknown',access='stream') + write(47) params%nutc/15,id2(1:20),id2(1:180000) + close(47) +! print*,'Calling SuperFox decoder',params%nzhsym,params%b_superfox, & +! params%b_even_seq,params%nutc + call execute_command_line('.\sf fort.47') else call timer('decft8 ',0) newdat=params%newdat diff --git a/lib/superfox/decode_sf.f90 b/lib/superfox/decode_sf.f90 new file mode 100644 index 000000000..853b0b88f --- /dev/null +++ b/lib/superfox/decode_sf.f90 @@ -0,0 +1,70 @@ +subroutine decode_sf(iwave) + + use sfox_mod + integer*2 iwave(NMAX) + integer msg1(0:47) + integer, allocatable :: rxdat(:) + integer, allocatable :: rxprob(:) + integer, allocatable :: rxdat2(:) + integer, allocatable :: rxprob2(:) + integer, allocatable :: correct(:) + real a(3) + real, allocatable :: s3(:,:) !Symbol spectra: will be s3(NQ,NN) + complex crcvd(NMAX) + integer isync(24) !Symbol numbers for sync tones + data isync/ 1, 2, 4, 7,11,16,22,29,37,39, & + 42,43,45,48,52,57,63,70,78,80, & + 83,84,86,89/ + +! Temporary, for initial tests: + data msg1/ 5, 126, 55, 29, 5, 127, 86, 117, 6, 0, & + 118, 77, 6, 2, 22, 37, 6, 3, 53, 125, & + 1, 27, 124, 110, 54, 12, 9, 43, 43, 64, & + 96, 94, 85, 92, 6, 7, 21, 5, 104, 48, & + 67, 37, 110, 67, 4, 106, 26, 64/ + + mm0=7 !Symbol size (bits) + nn0=127 !Number of information + parity symbols + kk0=48 !Number of information symbols + fspread=0.0 + delay=0.0 + fsample=12000.0 !Sample rate (Hz) + ns0=24 !Number of sync symbols + call sfox_init(mm0,nn0,kk0,'no',fspread,delay,fsample,ns0) + +! Allocate storage for arrays that depend on code parameters. + allocate(s3(0:NQ-1,0:NN-1)) + allocate(rxdat(0:NN-1)) + allocate(rxprob(0:NN-1)) + allocate(rxdat2(0:NN-1)) + allocate(rxprob2(0:NN-1)) + allocate(correct(0:NN-1)) + + call rs_init_sf(MM,NQ,NN,KK,NFZ) !Initialize the Karn codec + + call sfox_ana(iwave,NMAX,crcvd,NMAX) + + call sfox_sync(iwave,fsample,isync,f,t,fwidth) !Find freq, DT, width + + a=0. + a(1)=1500.0-f - baud !Shift frequencies down by one bin + call twkfreq(crcvd,crcvd,NMAX,fsample,a) + f=1500.0 + call sfox_demod(crcvd,f,t,isync,s3) !Get s3(0:NQ-1,0:127) + call sfox_prob(s3,rxdat,rxprob,rxdat2,rxprob2) + + do i=0,KK-1 + write(60,3060) i,msg1(i),rxdat(i),rxprob(i),rxdat2(i),rxprob2(i) +3060 format(6i8) + enddo + + ntrials=1000 + call ftrsd3(s3,rxdat,rxprob,rxdat2,rxprob2,ntrials, & + correct,param,ntry) + if(ntry.lt.ntrials) then + print*,'A',ntry,count(rxdat(0:KK-1).ne.msg1),count(correct(0:KK-1).ne.msg1) + call sfox_unpack(correct(0:KK-1)) + endif + + return +end subroutine decode_sf diff --git a/lib/superfox/ftrsd3.f90 b/lib/superfox/ftrsd3.f90 new file mode 100644 index 000000000..2db1b82c7 --- /dev/null +++ b/lib/superfox/ftrsd3.f90 @@ -0,0 +1,186 @@ +subroutine ftrsd3(s3,rxdat,rxprob,rxdat2,rxprob2,ntrials0,correct,param,ntry) + +! Soft-decision decoder for Reed-Solomon codes. + +! This decoding scheme is built around Phil Karn's Berlekamp-Massey +! errors and erasures decoder. The approach is inspired by a number of +! publications, including the stochastic Chase decoder described +! in "Stochastic Chase Decoding of Reed-Solomon Codes", by Leroux et al., +! IEEE Communications Letters, Vol. 14, No. 9, September 2010 and +! "Soft-Decision Decoding of Reed-Solomon Codes Using Successive Error- +! and-Erasure Decoding," by Soo-Woong Lee and B. V. K. Vijaya Kumar. + +! Steve Franke K9AN and Joe Taylor K1JT + + use sfox_mod + + real s3(0:NQ-1,0:NN-1) !Symbol spectra + integer rxdat(0:NN-1) !Hard-decision symbol values + integer rxprob(0:NN-1) !Probabilities that rxdat values are correct + integer rxdat2(0:NN-1) !Second most probable symbol values + integer rxprob2(0:NN-1) !Probabilities that rxdat2 values are correct + integer workdat(0:NN-1) !Work array + integer correct(0:NN-1) !Corrected codeword + integer indexes(0:NN-1) !For sorting probabilities + integer probs(0:NN-1) !Temp array for sorting probabilities + integer thresh0(0:NN-1) !Temp array for thresholds + integer era_pos(0:NN-KK-1) !Index values for erasures + integer param(0:8) + integer*8 nseed,ir !No unsigned int in Fortran + integer pass,tmp,thresh + + integer perr(0:7,0:7) + data perr/ 4, 9,11,13,14,14,15,15, & + 2,20,20,30,40,50,50,50, & + 7,24,27,40,50,50,50,50, & + 13,25,35,46,52,70,50,50, & + 17,30,42,54,55,64,71,70, & + 25,39,48,57,64,66,77,77, & + 32,45,54,63,66,75,78,83, & + 51,58,57,66,72,77,82,86/ + + ntrials=ntrials0 + nhard=0 + nhard_min=32768 + nsoft=0 + nsoft_min=32768 + ntotal=0 + ntotal_min=32768 + nera_best=0 + nsym=nn + + do i=0,NN-1 + indexes(i)=i + probs(i)=rxprob(i) + enddo + + do pass=1,nsym-1 + do k=0,nsym-pass-1 + if(probs(k).lt.probs(k+1)) then + tmp=probs(k) + probs(k)=probs(k+1) + probs(k+1)=tmp + tmp=indexes(k) + indexes(k)=indexes(k+1) + indexes(k+1)=tmp + endif + enddo + enddo + + correct=-1 + era_pos=0 + numera=0 + workdat=rxdat + call rs_decode_sf(workdat,era_pos,numera,nerr) !Call the decoder + nerr=-1 + + if(nerr.ge.0) then +! Hard-decision decoding succeeded. Save codeword and some parameters. + nhard=count(workdat.ne.rxdat) + correct=workdat + param(0)=0 + param(1)=nhard + param(2)=0 + param(3)=0 + param(4)=0 + param(5)=0 + param(7)=1000*1000 !??? + ntry=0 + go to 900 + endif + +! Hard-decision decoding failed. Try the FT soft-decision method. +! Generate random erasure-locator vectors and see if any of them +! decode. This will generate a list of "candidate" codewords. The +! soft distance between each candidate codeword and the received +! word is estimated by finding the largest (pp1) and second-largest +! (pp2) outputs from a synchronized filter-bank operating on the +! symbol spectra, and using these to decide which candidate +! codeword is "best". + + nseed=1 !Seed for random numbers + ncandidates=0 + nsum=0 + do i=0,NN-1 + nsum=nsum+rxprob(i) + j=indexes(NN-1-i) + ratio=float(rxprob2(j))/(float(rxprob(j))+0.01) + ii=7.999*ratio + jj=int((7.999/NN)*(NN-1-i)) + thresh0(i)=0.90*perr(jj,ii) + enddo + if(nsum.le.0) return + + pp1=0. + pp2=0. + do k=1,ntrials + era_pos=0 + workdat=rxdat + +! Mark a subset of the symbols as erasures. +! Run through the ranked symbols, starting with the worst, i=0. +! NB: j is the symbol-vector index of the symbol with rank i. + + numera=0 + do i=0,NN-1 + j=indexes(NN-1-i) + thresh=thresh0(i) +! Generate a random number ir, 0 <= ir <= 100 (see POSIX.1-2001 example). + ir=100.0*ran1(nseed) + if((ir.lt.thresh) .and. numera.lt. 0.69*(NN-KK)) then + era_pos(numera)=j + numera=numera+1 + endif + enddo + call rs_decode_sf(workdat,era_pos,numera,nerr) !Call the decoder + if( nerr.ge.0) then + ! We have a candidate codeword. Find its hard and soft distance from + ! the received word. Also find pp1 and pp2 from the full array + ! s3(NQ,NN) of synchronized symbol spectra. + ncandidates=ncandidates+1 + nhard=0 + nsoft=0 + do i=0,NN-1 + if(workdat(i).ne. rxdat(i)) then + nhard=nhard+1; + if(workdat(i) .ne. rxdat2(i)) nsoft=nsoft+rxprob(i) + endif + enddo + nsoft=NN*nsoft/nsum + ntotal=nsoft+nhard + + pp=0. + call getpp3(s3,workdat,pp) +! write(*,5001) ncandidates,nhard,nsoft,ntotal,pp,pp1,pp2 +!5001 format(4i8,3f7.3) + if(pp.gt.pp1) then + pp2=pp1 + pp1=pp + nsoft_min=nsoft + nhard_min=nhard + ntotal_min=ntotal + correct=workdat + nera_best=numera + ntry=k + else + if(pp.gt.pp2 .and. pp.ne.pp1) pp2=pp + endif + if(nhard_min.le.60 .and. ntotal_min.le.90) exit !### Needs tuning + endif + if(k.eq.ntrials) ntry=k + enddo + + param(0)=ncandidates + param(1)=nhard_min + param(2)=nsoft_min + param(3)=nera_best + param(4)=1000 + if(pp1.gt.0.0) param(4)=1000.0*pp2/pp1 + param(5)=ntotal_min + param(6)=ntry + param(7)=1000.0*pp2 + param(8)=1000.0*pp1 + if(param(0).eq.0) param(2)=-1 + +900 return +end subroutine ftrsd3 diff --git a/lib/superfox/getpp3.f90 b/lib/superfox/getpp3.f90 new file mode 100644 index 000000000..38d96357d --- /dev/null +++ b/lib/superfox/getpp3.f90 @@ -0,0 +1,22 @@ +subroutine getpp3(s3,workdat,p) + + use sfox_mod + real s3(NQ,NN) + integer workdat(NN) + integer a(NN) + +! a(1:NN)=workdat(NN:1:-1) + a=workdat + + psum=0. + do j=1,NN + i=a(j)+1 + x=s3(i,j) + s3(i,j)=0. + psum=psum + x + s3(i,j)=x + enddo + p=psum/NN + + return +end subroutine getpp3 diff --git a/lib/superfox/ran1.f90 b/lib/superfox/ran1.f90 new file mode 100644 index 000000000..b09b9baba --- /dev/null +++ b/lib/superfox/ran1.f90 @@ -0,0 +1,28 @@ +FUNCTION ran1(idum) + INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV + REAL ran1,AM,EPS,RNMX + PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836, & + NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS) + INTEGER j,k,iv(NTAB),iy + SAVE iv,iy + DATA iv /NTAB*0/, iy /0/ + if (idum.le.0.or.iy.eq.0) then + idum=max(-idum,1) + do j=NTAB+8,1,-1 + k=idum/IQ + idum=IA*(idum-k*IQ)-IR*k + if (idum.lt.0) idum=idum+IM + if (j.le.NTAB) iv(j)=idum + enddo + iy=iv(1) + endif + k=idum/IQ + idum=IA*(idum-k*IQ)-IR*k + if (idum.lt.0) idum=idum+IM + j=1+iy/NDIV + iy=iv(j) + iv(j)=idum + ran1=min(AM*iy,RNMX) + + return +END FUNCTION ran1 diff --git a/lib/superfox/rs_sf.a b/lib/superfox/rs_sf.a new file mode 100644 index 0000000000000000000000000000000000000000..81b0c949ecfa03ff2e68e9ebb8c9253e4f62f040 GIT binary patch literal 8296 zcmcgxe{dAneczQ9a{^f_jGPeL8BUXHDiy5g3vnq|v)|qB<~f2O%*OK*`*EESg!?xL!o$0BPwTFhmZn`DO&zVm*6mw^S+Zq& zV8@nR{g#%z7qD7<9j#3Q=qzLQEz3!7Q^IO?kp7O&-a?(dwYAHZEknDq0>4z6=WIE6 zAJi3uKfqswMM*+Uuw_rsqV1oFVnL8ZLD)9Se!m}oZlR&MX3OT_<~+Wc$W=l^wC354 zj+xr>>-H1hou!U+VQ`;oj2aC7WO}`z= z=OK_Z70rGVk z&{AE)c>D2U1cII_gPqmb9!X>6Snrzs8G^y5 zFQO$`<2p3Z*b48EHFtGDdwddF2=YQ`5*MRtILg*W1i5EEP0>3M>KC=j zw5q(bZ`2!Ae@yDj3wx-CptFr(K^I5#a8mAh3oW|$mqsWm>Xjd>%D}!+HT1D4hd)6B zCUpWGrnKDXsfrk)5g0ULSBtXxqRyIO3^QCXrSDIt`m=9&G#030^=^$NUr}#PhKJZC zx#tWxU6jOtna;;bg{gPzq>g4F?)HWDWAYmHqB zjfqBV6y}G`eSx#Bx%u+7VL`mDDYsc-Zw2h7Ri{}$jC1+6n9!7T*GI;2O?=-}u%0T{ z!h*^~xn~q&>jt!CE|TvHU67uNNXy2{Nnf3n`$A)qsayyezJQ?fhRJX=nBrzw;9r+a z{cXX{`!)Nkj@%NPUmTZcP#J)L_SyJ-;z>K zWHK3(r7Qy*?27f@T9ccPcf?>xY@ZcRKhv~H-KcI-pHZ96PPI?1rTkL)rd zIIBpQ^_W;ImX`?S&z9ECZJx8Uq`?)vFCmQ6RxE2; zw)D}RJ36-Rd9pHh?3%Z|u-tFJIea-qu*DqV@-nnbXpSJ{p)$SEF z-Wmub+?$gJRK*-@=j~Uo2~h6k+T* zna}xq*%xXl7XK)}$)YLVPgTMiLSXZ@Z95!uQBC_DHDs5XT@>LwFTC)AF!P}yk0Jm3 zK=J=GqrdktoC~N-MgL#Jq~v?#o71CzGb~M@b&kxALVv%D=ue)vB|9dKOJ%0$PhPtR ze90^?t^D1he@~W&WKp5$pBKribnkn97Bd4XjsyQ4KomRx?24tB`Jsh8oWk8cOxJkumU1U4V@Dq zHY-w7uF8Es1X~?1{X2C@iOKbHL@<><(yLq!o(QGI?(dm=lgmC-n7na;Rvk5ttGx6@ z@d|Li6)LI35p&nmqPVk>3|Q|iE%Y0D);GtNT@RYngM|r(giui z4=KrgWz}cdaOGJ9pUXSBTz^AC(E2@4OG6=6EVKaVhWCUR$oQ-qVKDQ8@j_Ae0ok-^Fp#m)X zEF{a{DKwd4k>*tBV!u7yqWw`=?!+zH8?+w^%Oly>yWv)QOM6F4FxSUIsh0Nk9lLX2 zMm4{W=DUDUe*K*Mgr)$UrsY$)++@k)k1w5$POC8ud$xC?Q!SR1(CPaWq)^K6{4UUm zY`8Q#CXGwwZs^nizA9P@l>ZH#Hr^qxVO3);*`QqG!1Z2;Fd~2N*^hsXZ*{8?T{KwZ z27~>~2qc)EG-B6^Or|#*ktbx0`P+=p)d_>eE&MX7&lwTT^CcrThC`=F7_V9h-sgbdHZdBSB4- zy{j$iuGPaIlENe?A@_Y7$+f}t4Qj-n3Ee8vxF8?A=AzU72d%2SJ#Mnm*at zaUB0WYA7SheSbp+VqIy2cWqGriMF{OY0Nn23KhFsR`b$-^bU>B$9Os`JqCM$K?@G&k92=O7Ve;MO{)l!&f(5*TMxY!smP#(Eyy9in zFzHjGo#RL_dp3d{z#cUsYb5WG!I1yxC+&shdVD~#e1Rs&^)}=!ee*yY&=&!rnLkzK zWBrl(jLKq*ikLqwpNd)9msKZOe`S(V6R3W}IBv3`P*YDPK*7@aMulgWc{R{-SSrtkKa$Zi!YRt{L01PSCO1Rhl5Hoc%4++`*|#x zCJ1EgC>cAgiBa3w;8MLRg0R)I^MWsQ7B2%3^Y5J@jkm`Q4kdY4)Z{HEOzxl1SjrnU zBT|HBfmCERSz2WkBH7(CH!KK8Ms$TsF+D9G>4hy(iOr?53vl8IOmo7B7&uUzGnhZ_ z_!IZ1y+gWp(B$jNwa7Oz;X#un_B}WbhW!`~&fG6X>N2Wwq1!Y9rQd~bU56o|Q?fc? zMfO9O@)da3VUu^eraWxgq)Cpxr1X$2Yu}Ijfx#yu(qWA+P?Zt6rx8=*vn@~ak20~l z&18*9-HUsHi{lSye7k!MGhc+2kc2DeyI0n8sb@7e58;Fpyb3R-p_L>r!)i@)4;nMS z@4j;|6rmWLB?=BO{0Nk-#aL-Yv zKBT8b6&nRIu+Qbbe?Yp2SzzTRW6+DC$JgQt~6^Aq# zxTwSoAVcSW06<5vSf??#?#*DJlQa-5o3(8;zI$!Y}Jki!gOaJU8d%e@cTA^{KO8o*Kaj+ngjRZJZ# zqQ`w~ASd6?uKF1~XK4HZYyaUQLG*=U5^l)$tWmj0`xEEfL&}xldm?t)I@@o)kn2zT z^*rrY>-1`YV@`y&q~b-0dH!{OU|M-x(?r)H4G@W0B=vAc?tKo2GMEOA zgn*={ae&3lcDu~T=1Bm}0o12apQO7fu8(0h8n!6bJE`%Ny}V49PYsOEv-CGYA4qz2 z+Nb<$a5*o6_P?-jWbcCP{v>^4r0*>(q>V4xF8Y#34x(GEb0a=khwPSvzpt z8ajKI1APWaAL}kBXFsdUlEnwfe1m#;B=`l)I{~kS$yV3aap@@K1d$;NaN5%%TGMo; z`6YreL}(Tx3`o6&-DO3|mxj}s%z!=$e!$$Q5ovKD9mrkvY|6junRzd#>L)Xq>NvJ4 z)7#|UO1$|3=eJcdYJWB!w{}XV3n3rkQNs6L- z-@XgpkN*^ro+6V!Fyaq<=nGW5yoQ(lm0n!>ZH-+NhwB$r{1B^>8?Sm;WtS?4t@XT! zfb$W2Wdr3~fgwWMikFQL?BlI~4_PK+KZi|>P**b|-VKnj2HgR^Pqn-;L05j<-;Cg4 zP0NCj)I~I6VN)Fz6Ja?XY(Uz4saIeMZ#4DfO~=5gt(Z29jcuAq^`TFGz9l|9H~et% zd~u$Tz1ND)PRtrDo)8;oZ3MwXDGp414t36ahx`vbft}Ci`cOu=6fb!hp<`tF+?7|~)YKGcd#)q+TvK}o zebwq{wLe-}NuP6fb(icAU)vs^o*nkXbLc*md>}s~ryBA;$_^KUhYDr9_mUlw$4@W& zXZUhO+3JrlZho05JM@4d4!*XmesF2wZg*GoJte=;-4%U(Qn|oOaHzeR&=-ltMWI`> z&oJ->9|+xYov8<)Llo$U%|A zlYV5AVoN)_j(o$WsGr?R_r`V|?O&U6JFR!#L)rTw0FM0fS@%#@uWSdL%d9&Fd~*iX zj7Jn`cT%awvmcLLlfPGJJ&HxLHtmq9@|)oqS*Nh%f7vCA7K>z~!gA;Bl?MzZg-m!` zo~^tf1zbc*{d6nAG~Z@mZhV^SW64_EI?4iG+OxKse`r E0E8e+RR910 literal 0 HcmV?d00001 diff --git a/lib/superfox/sfox_ana.f90 b/lib/superfox/sfox_ana.f90 new file mode 100644 index 000000000..51c5a38c3 --- /dev/null +++ b/lib/superfox/sfox_ana.f90 @@ -0,0 +1,21 @@ +subroutine sfox_ana(iwave,npts,c0,npts2) + + use timer_module, only: timer + + integer*2 iwave(npts) !Raw data at 12000 Hz + complex c0(0:npts2-1) !Complex data at 6000 Hz + save + + nfft1=npts +! nfft2=nfft1/2 + nfft2=nfft1 +! df1=12000.0/nfft1 + fac=2.0/(32767.0*nfft1) + c0(0:npts-1)=fac*iwave(1:npts) + call four2a(c0,nfft1,1,-1,1) !Forward c2c FFT + c0(nfft2/2+1:nfft2-1)=0. !Remove negative frequencies + c0(0)=0.5*c0(0) !Scale the DC term to 1/2 + call four2a(c0,nfft2,1,1,1) !Inverse c2c FFT; c0 is analytic sig + + return +end subroutine sfox_ana diff --git a/lib/superfox/sfox_demod.f90 b/lib/superfox/sfox_demod.f90 new file mode 100644 index 000000000..98167611f --- /dev/null +++ b/lib/superfox/sfox_demod.f90 @@ -0,0 +1,33 @@ +subroutine sfox_demod(crcvd,f,t,isync,s3) + + use sfox_mod + complex crcvd(NMAX) !Signal as received + complex c(0:NSPS-1) !Work array, one symbol long + real s3(0:NQ-1,0:NN-1) !Synchronized symbol spectra + integer isync(44) +! integer ipk(1) + + j0=nint(12000.0*(t+0.5)) + df=12000.0/NSPS + i0=nint(f/df)-NQ/2 + k=-1 + do n=1,NDS !Loop over all symbols + if(any(isync(1:NS).eq.n)) cycle + jb=n*NSPS + j0 + ja=jb-NSPS+1 + if(ja.lt.1 .or. jb.gt.NMAX) cycle + k=k+1 + c=crcvd(ja:jb) + call four2a(c,NSPS,1,-1,1) !Compute symbol spectrum + do i=0,NQ-1 + s3(i,k)=real(c(i0+i))**2 + aimag(c(i0+i))**2 + enddo +! ipk=maxloc(s3(0:NQ-1,k)) +! if(k.lt.10) print*,'AAA',k,ipk(1)-1 + enddo + + call pctile(s3,NQ*NN,50,base) + s3=s3/base + + return +end subroutine sfox_demod diff --git a/lib/superfox/sfox_prob.f90 b/lib/superfox/sfox_prob.f90 new file mode 100644 index 000000000..c8e9f25c9 --- /dev/null +++ b/lib/superfox/sfox_prob.f90 @@ -0,0 +1,55 @@ +subroutine sfox_prob(s3,rxdat,rxprob,rxdat2,rxprob2) + +! Demodulate the 64-bin spectra for each of 63 symbols in a frame. + +! Parameters +! rxdat most reliable symbol value +! rxdat2 second most likely symbol value +! rxprob probability that rxdat was the transmitted value +! rxprob2 probability that rxdat2 was the transmitted value + + use sfox_mod + implicit real*8 (a-h,o-z) + real*4 s3(0:NQ-1,0:NN-1) + integer rxdat(0:NN-1),rxprob(0:NN-1),rxdat2(0:NN-1),rxprob2(0:NN-1) + + afac=1.1 +! scale=255.999 + scale=2047.999 + +! Compute average spectral value + ave=sum(s3)/(NQ*ND) + i1=1 !Silence warning + i2=1 + +! Compute probabilities for most reliable symbol values + do j=0,NN-1 !Loop over all symbols + s1=-1.e30 + psum=0. + do i=0,NQ-1 !Loop over frequency bins + x=min(afac*s3(i,j)/ave,50.d0) + psum=psum+s3(i,j) + if(s3(i,j).gt.s1) then + s1=s3(i,j) !Find max signal+noise power + i1=i !Find most reliable symbol value + endif + enddo + if(psum.eq.0.0) psum=1.e-6 !Guard against zero signal+noise + + s2=-1.e30 + do i=0,NQ-1 + if(i.ne.i1 .and. s3(i,j).gt.s2) then + s2=s3(i,j) !Second largest signal+noise power + i2=i !Bin number for second largest power + endif + enddo + p1=s1/psum !p1, p2 are symbol metrics for ftrsd + p2=s2/psum + rxdat(j)=i1 + rxdat2(j)=i2 + rxprob(j)=scale*p1 !Scaled probabilities, 0 - 255 + rxprob2(j)=scale*p2 + enddo + + return +end subroutine sfox_prob diff --git a/lib/superfox/sfox_sync.f90 b/lib/superfox/sfox_sync.f90 new file mode 100644 index 000000000..c080fa516 --- /dev/null +++ b/lib/superfox/sfox_sync.f90 @@ -0,0 +1,153 @@ +subroutine sfox_sync(iwave,fsample,isync,f,t,fwidth) + + use sfox_mod + parameter (NSTEP=8) + integer*2 iwave(0:NMAX-1) + integer isync(44) + integer ipeak(2) + integer ipeak2(1) + complex, allocatable :: c(:) !Work array + real, allocatable :: s(:,:) !Symbol spectra, stepped by NSTEP + real, allocatable :: savg(:) !Average spectrum + real, allocatable :: ccf(:,:) + real, allocatable :: s2(:) !Fine spectrum of sync tone + + nfft=nsps + nh=nfft/2 + istep=NSPS/NSTEP + jz=(13.5*fsample)/istep + df=fsample/nfft + dtstep=istep/fsample + fsync=1500.0-bw/2 + ftol=50.0 + ia=nint((fsync-ftol)/df) + ib=nint((fsync+ftol)/df) + lagmax=1.5/dtstep + lag1=-lagmax + lag2=lagmax + + allocate(s(0:nh/2,jz)) + allocate(savg(0:nh/2)) + allocate(c(0:nfft-1)) + allocate(ccf(ia:ib,lag1:lag2)) + + s=0. + savg=0. + fac=1.0/nfft + +! Compute symbol spectra with df=baud/2 and NSTEP steps per symbol. + do j=1,jz + i1=(j-1)*istep + i2=i1+nsps-1 + k=-1 + do i=i1,i2,2 !Load iwave data into complex array c0, for r2c FFT + xx=iwave(i) + yy=iwave(i+1) + k=k+1 + c(k)=fac*cmplx(xx,yy) + enddo + c(k+1:)=0. + call four2a(c,nfft,1,-1,0) !r2c FFT + do i=1,nh/2 + s(i,j)=real(c(i))**2 + aimag(c(i))**2 + savg(i)=savg(i) + s(i,j) + enddo + enddo + savg=savg/jz + + ccfbest=0. + ibest=0 + lagpk=0 + lagbest=0 + j0=0.5/dtstep !Nominal start-signal index + + do i=ia,ib + ccfmax=0. + do lag=lag1,lag2 + ccft=0. + do m=1,NS + k=isync(m) + n=NSTEP*(k-1) + 1 + j=n+lag+j0 + if(j.ge.1 .and. j.le.jz) ccft=ccft + s(i,j) + enddo ! m + ccft=ccft - NS*savg(i) + ccf(i,lag)=ccft + if(ccft.gt.ccfmax) then + ccfmax=ccft + lagpk=lag + endif + enddo ! lag + + if(ccfmax.gt.ccfbest) then + ccfbest=ccfmax + ibest=i + lagbest=lagpk + endif + enddo ! i + + ipeak=maxloc(ccf) + ipk=ipeak(1)-1+ia + jpk=ipeak(2)-1+lag1 + + dxj=0. + if(jpk.gt.lag1 .and. jpk.lt.lag2) then + call peakup(ccf(ipk,jpk-1),ccf(ipk,jpk),ccf(ipk,jpk+1),dxj) + endif + + f=ibest*df + bw/2 + dxi*df + t=(lagbest+dxj)*dtstep + t=t-0.01 !### Why is this needed? ### + + nfft2=4*NSPS + deallocate(c) + allocate(c(0:nfft2-1)) + allocate(s2(0:nfft2-1)) + + i0=(t+0.5)*fsample + s2=0. + df2=fsample/nfft2 + do m=1,NS + i1=i0+(isync(m)-1)*NSPS + i2=i1+NSPS-1 + k=-1 + do i=i1,i2,2 !Load iwave data into complex array c0, for r2c FFT + if(i.gt.0) then + xx=iwave(i) + yy=iwave(i+1) + else + xx=0. + yy=0. + endif + k=k+1 + c(k)=fac*cmplx(xx,yy) + enddo + c(k+1:)=0. + call four2a(c,nfft2,1,-1,0) !r2c FFT + do i=1,nfft2/4 + s2(i)=s2(i) + real(c(i))**2 + aimag(c(i))**2 + enddo + enddo + + ia=nint((fsync-ftol)/df2) + ib=nint((fsync+ftol)/df2) + ipeak2=maxloc(s2(ia:ib)) + ipk=ipeak2(1)-1+ia + + dxi=0. + if(ipk.gt.1 .and. ipk.lt.nfft/4) then + call peakup(s2(ipk-1),s2(ipk),s2(ipk+1),dxi) + endif + f=(ipk+dxi)*df2 + bw/2.0 + fwidth=0. + + if(ipk.gt.100 .and. ipk.lt.nfft2/4-100) then + call pctile(s2(ipk-100:ipk+100),201,48,base) + s2=s2-base + smax=maxval(s2(ipk-10:ipk+10)) + w=count(s2(ipk-10:ipk+10).gt.0.5*smax) + if(w.gt.4.0) fwidth=sqrt(w*w - 4*4)*df2 + endif + + return +end subroutine sfox_sync diff --git a/lib/superfox/sfox_unpack.f90 b/lib/superfox/sfox_unpack.f90 new file mode 100644 index 000000000..27244de87 --- /dev/null +++ b/lib/superfox/sfox_unpack.f90 @@ -0,0 +1,51 @@ +subroutine sfox_unpack(imsg) + + use packjt77 + integer imsg(48) + character*336 msgbits + character*22 msg(10) + character*13 foxcall,c13 + character*4 crpt(5) + logical success + + write(msgbits,1000) imsg +1000 format(48b7.7) + read(msgbits(331:336),'(b6)') ntype !Message type + + if(ntype.eq.1) then !Get the Fox callsign + read(msgbits(271:328),'(b58)') n58 !Compound Fox call + call unpack28(n58,foxcall,success) + else + read(msgbits(303:330),'(b28)') n28 !Standard Fox call + call unpack28(n28,foxcall,success) + endif + + j=171 + do i=1,5 !Extract the reports + read(msgbits(j:j+3),'(b4)') n + if(n.eq.15) then + crpt(i)='RR73' + else + write(crpt(i),1006) 2*n-18 +1006 format(i3.2) + if(crpt(i)(1:1).eq.' ') crpt(i)(1:1)='+' + endif + j=j+32 + enddo + +! Unpack and format user-level messages: + do i=1,10 + j=28*i - 27 + if(i.gt.5) j=143 + (i-5)*32 + read(msgbits(j:j+27),'(b28)') n28 + if(n28.eq.0) cycle + call unpack28(n28,c13,success) + msg(i)=trim(c13)//' '//trim(foxcall) + if(i.le.5) msg(i)=trim(msg(i))//' RR73' + if(i.gt.5) msg(i)=trim(msg(i))//' '//crpt(i-5) + write(*,3001) i,trim(msg(i)) +3001 format(i2,2x,a) + enddo + + return +end subroutine sfox_unpack diff --git a/lib/sfox_wave.f90 b/lib/superfox/sfox_wave.f90 similarity index 100% rename from lib/sfox_wave.f90 rename to lib/superfox/sfox_wave.f90