1subroutine ft2_decode(cdatetime0,nfqso,iwave,ndecodes,mycall,hiscall,nrx,line) 2 3 use crc 4 use packjt77 5 include 'ft2_params.f90' 6 character message*37,c77*77 7 character*61 line 8 character*37 decodes(100) 9 character*120 data_dir 10 character*17 cdatetime0,cdatetime 11 character*6 mycall,hiscall,hhmmss 12 complex c2(0:NMAX/16-1) !Complex waveform 13 complex cb(0:NMAX/16-1) 14 complex cd(0:144*10-1) !Complex waveform 15 complex c1(0:9),c0(0:9) 16 complex ccor(0:1,144) 17 complex csum,cterm,cc0,cc1,csync1 18 real*8 fMHz 19 20 real a(5) 21 real rxdata(128),llr(128) !Soft symbols 22 real llr2(128) 23 real sbits(144),sbits1(144),sbits3(144) 24 real ps(0:8191),psbest(0:8191) 25 real candidate(3,100) 26 real savg(NH1) 27 integer*2 iwave(NMAX) !Generated full-length waveform 28 integer*1 message77(77),apmask(128),cw(128) 29 integer*1 hbits(144),hbits1(144),hbits3(144) 30 integer*1 s16(16) 31 logical unpk77_success 32 data s16/0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0/ 33 34 hhmmss=cdatetime0(8:13) 35 fs=12000.0/NDOWN !Sample rate 36 dt=1/fs !Sample interval after downsample (s) 37 tt=NSPS*dt !Duration of "itone" symbols (s) 38 baud=1.0/tt !Keying rate for "itone" symbols (baud) 39 txt=NZ*dt !Transmission length (s) 40 twopi=8.0*atan(1.0) 41 h=0.8 !h=0.8 seems to be optimum for AWGN sensitivity (not for fading) 42 43 dphi=twopi/2*baud*h*dt*16 ! dt*16 is samp interval after downsample 44 dphi0=-1*dphi 45 dphi1=+1*dphi 46 phi0=0.0 47 phi1=0.0 48 do i=0,9 49 c1(i)=cmplx(cos(phi1),sin(phi1)) 50 c0(i)=cmplx(cos(phi0),sin(phi0)) 51 phi1=mod(phi1+dphi1,twopi) 52 phi0=mod(phi0+dphi0,twopi) 53 enddo 54 the=twopi*h/2.0 55 cc1=cmplx(cos(the),-sin(the)) 56 cc0=cmplx(cos(the),sin(the)) 57 58 data_dir="." 59 fMHz=7.074 60 ncoh=1 61 candidate=0.0 62 ncand=0 63 fa=375.0 64 fb=3000.0 65 syncmin=0.2 66 maxcand=100 67 nfqso=-1 68 call getcandidates2a(iwave,fa,fb,maxcand,savg,candidate,ncand) 69 ndecodes=0 70 do icand=1,ncand 71 f0=candidate(1,icand) 72 if( f0.le.375.0 .or. f0.ge.(5000.0-375.0) ) cycle 73 call ft2_downsample(iwave,f0,c2) ! downsample from 160s/Symbol to 10s/Symbol 74! 750 samples/second here 75 ibest=-1 76 sybest=-99. 77 dfbest=-1. 78!### do if=-15,+15 79 do if=-30,30 80 df=if 81 a=0. 82 a(1)=-df 83 call twkfreq1(c2,NMAX/16,fs,a,cb) 84 do is=0,374 !DT search range is 0 - 0.5 s 85 csync1=0. 86 cterm=1 87 do ib=1,16 88 i1=(ib-1)*10+is 89 i2=i1+136*10 90 if(s16(ib).eq.1) then 91 csync1=csync1+sum(cb(i1:i1+9)*conjg(c1(0:9)))*cterm 92 cterm=cterm*cc1 93 else 94 csync1=csync1+sum(cb(i1:i1+9)*conjg(c0(0:9)))*cterm 95 cterm=cterm*cc0 96 endif 97 enddo 98 if(abs(csync1).gt.sybest) then 99 ibest=is 100 sybest=abs(csync1) 101 dfbest=df 102 endif 103 enddo 104 enddo 105 106 a=0. 107 a(1)=-dfbest 108 call twkfreq1(c2,NMAX/16,fs,a,cb) 109 ib=ibest 110 cd=cb(ib:ib+144*10-1) 111 s2=sum(real(cd*conjg(cd)))/(10*144) 112 cd=cd/sqrt(s2) 113 do nseq=1,5 114 if( nseq.eq.1 ) then ! noncoherent single-symbol detection 115 sbits1=0.0 116 do ibit=1,144 117 ib=(ibit-1)*10 118 ccor(1,ibit)=sum(cd(ib:ib+9)*conjg(c1(0:9))) 119 ccor(0,ibit)=sum(cd(ib:ib+9)*conjg(c0(0:9))) 120 sbits1(ibit)=abs(ccor(1,ibit))-abs(ccor(0,ibit)) 121 hbits1(ibit)=0 122 if(sbits1(ibit).gt.0) hbits1(ibit)=1 123 enddo 124 sbits=sbits1 125 hbits=hbits1 126 sbits3=sbits1 127 hbits3=hbits1 128 elseif( nseq.ge.2 ) then 129 nbit=2*nseq-1 130 numseq=2**(nbit) 131 ps=0 132 do ibit=nbit/2+1,144-nbit/2 133 ps=0.0 134 pmax=0.0 135 do iseq=0,numseq-1 136 csum=0.0 137 cterm=1.0 138 k=1 139 do i=nbit-1,0,-1 140 ibb=iand(iseq/(2**i),1) 141 csum=csum+ccor(ibb,ibit-(nbit/2+1)+k)*cterm 142 if(ibb.eq.0) cterm=cterm*cc0 143 if(ibb.eq.1) cterm=cterm*cc1 144 k=k+1 145 enddo 146 ps(iseq)=abs(csum) 147 if( ps(iseq) .gt. pmax ) then 148 pmax=ps(iseq) 149 ibflag=1 150 endif 151 enddo 152 if( ibflag .eq. 1 ) then 153 psbest=ps 154 ibflag=0 155 endif 156 call getbitmetric(2**(nbit/2),psbest,numseq,sbits3(ibit)) 157 hbits3(ibit)=0 158 if(sbits3(ibit).gt.0) hbits3(ibit)=1 159 enddo 160 sbits=sbits3 161 hbits=hbits3 162 endif 163 nsync_qual=count(hbits(1:16).eq.s16) 164 if(nsync_qual.lt.10) exit 165 rxdata=sbits(17:144) 166 rxav=sum(rxdata(1:128))/128.0 167 rx2av=sum(rxdata(1:128)*rxdata(1:128))/128.0 168 rxsig=sqrt(rx2av-rxav*rxav) 169 rxdata=rxdata/rxsig 170 sigma=0.80 171 llr(1:128)=2*rxdata/(sigma*sigma) 172 apmask=0 173 max_iterations=40 174 do ibias=0,0 175 llr2=llr 176 if(ibias.eq.1) llr2=llr+0.4 177 if(ibias.eq.2) llr2=llr-0.4 178 call bpdecode128_90(llr2,apmask,max_iterations,message77,cw,nharderror,niterations) 179 if(nharderror.ge.0) exit 180 enddo 181 nhardmin=-1 182 if(sum(message77).eq.0) cycle 183 if( nharderror.ge.0 ) then 184 write(c77,'(77i1)') message77(1:77) 185 call unpack77(c77,nrx,message,unpk77_success) 186 idupe=0 187 do i=1,ndecodes 188 if(decodes(i).eq.message) idupe=1 189 enddo 190 if(idupe.eq.1) exit 191 ndecodes=ndecodes+1 192 decodes(ndecodes)=message 193 xsnr=db(sybest*sybest) - 115.0 !### Rough estimate of S/N ### 194 nsnr=nint(xsnr) 195 freq=f0+dfbest 196 write(line,1000) hhmmss,nsnr,ibest/750.0,nint(freq),message 1971000 format(a6,i4,f5.2,i5,' + ',1x,a37) 198 open(24,file='all_ft2.txt',status='unknown',position='append') 199 write(24,1002) cdatetime0,nsnr,ibest/750.0,nint(freq),message, & 200 nseq,nharderror,nhardmin 201 if(hhmmss.eq.' ') write(*,1002) cdatetime0,nsnr, & 202 ibest/750.0,nint(freq),message,nseq,nharderror,nhardmin 2031002 format(a17,i4,f6.2,i5,' Rx ',a37,3i5) 204 close(24) 205 206!### Temporary: assume most recent decoded message conveys "hiscall". 207 i0=index(message,' ') 208 if(i0.ge.3 .and. i0.le.7) then 209 hiscall=message(i0+1:i0+6) 210 i1=index(hiscall,' ') 211 if(i1.gt.0) hiscall=hiscall(1:i1) 212 endif 213 nrx=-1 214 if(index(message,'CQ ').eq.1) nrx=1 215 if((index(message,trim(mycall)//' ').eq.1) .and. & 216 (index(message,' '//trim(hiscall)//' ').ge.4)) then 217 if(index(message,' 559 ').gt.8) nrx=2 218 if(index(message,' R 559 ').gt.8) nrx=3 219 if(index(message,' RR73 ').gt.8) nrx=4 220 endif 221!### 222 exit 223 endif 224 enddo ! nseq 225 enddo !candidate list 226 227 return 228end subroutine ft2_decode 229 230subroutine getbitmetric(ib,ps,ns,xmet) 231 real ps(0:ns-1) 232 xm1=0 233 xm0=0 234 do i=0,ns-1 235 if( iand(i/ib,1) .eq. 1 .and. ps(i) .gt. xm1 ) xm1=ps(i) 236 if( iand(i/ib,1) .eq. 0 .and. ps(i) .gt. xm0 ) xm0=ps(i) 237 enddo 238 xmet=xm1-xm0 239 return 240end subroutine getbitmetric 241 242subroutine downsample2(ci,f0,co) 243 parameter(NI=144*160,NH=NI/2,NO=NI/16) ! downsample from 200 samples per symbol to 10 244 complex ci(0:NI-1),ct(0:NI-1) 245 complex co(0:NO-1) 246 fs=12000.0 247 df=fs/NI 248 ct=ci 249 call four2a(ct,NI,1,-1,1) !c2c FFT to freq domain 250 i0=nint(f0/df) 251 ct=cshift(ct,i0) 252 co=0.0 253 co(0)=ct(0) 254 b=8.0 255 do i=1,NO/2 256 arg=(i*df/b)**2 257 filt=exp(-arg) 258 co(i)=ct(i)*filt 259 co(NO-i)=ct(NI-i)*filt 260 enddo 261 co=co/NO 262 call four2a(co,NO,1,1,1) !c2c FFT back to time domain 263 return 264end subroutine downsample2 265 266subroutine ft2_downsample(iwave,f0,c) 267 268! Input: i*2 data in iwave() at sample rate 12000 Hz 269! Output: Complex data in c(), sampled at 1200 Hz 270 271 include 'ft2_params.f90' 272 parameter (NFFT2=NMAX/16) 273 integer*2 iwave(NMAX) 274 complex c(0:NMAX/16-1) 275 complex c1(0:NFFT2-1) 276 complex cx(0:NMAX/2) 277 real x(NMAX) 278 equivalence (x,cx) 279 280 BW=4.0*75 281 df=12000.0/NMAX 282 x=iwave 283 call four2a(x,NMAX,1,-1,0) !r2c FFT to freq domain 284 ibw=nint(BW/df) 285 i0=nint(f0/df) 286 c1=0. 287 c1(0)=cx(i0) 288 do i=1,NFFT2/2 289 arg=(i-1)*df/bw 290 win=exp(-arg*arg) 291 c1(i)=cx(i0+i)*win 292 c1(NFFT2-i)=cx(i0-i)*win 293 enddo 294 c1=c1/NFFT2 295 call four2a(c1,NFFT2,1,1,1) !c2c FFT back to time domain 296 c=c1(0:NMAX/16-1) 297 return 298end subroutine ft2_downsample 299