1module fst4_decode 2 3 type :: fst4_decoder 4 procedure(fst4_decode_callback), pointer :: callback 5 contains 6 procedure :: decode 7 end type fst4_decoder 8 9 abstract interface 10 subroutine fst4_decode_callback (this,nutc,sync,nsnr,dt,freq, & 11 decoded,nap,qual,ntrperiod,lwspr,fmid,w50) 12 import fst4_decoder 13 implicit none 14 class(fst4_decoder), intent(inout) :: this 15 integer, intent(in) :: nutc 16 real, intent(in) :: sync 17 integer, intent(in) :: nsnr 18 real, intent(in) :: dt 19 real, intent(in) :: freq 20 character(len=37), intent(in) :: decoded 21 integer, intent(in) :: nap 22 real, intent(in) :: qual 23 integer, intent(in) :: ntrperiod 24 logical, intent(in) :: lwspr 25 real, intent(in) :: fmid 26 real, intent(in) :: w50 27 end subroutine fst4_decode_callback 28 end interface 29 30contains 31 32 subroutine decode(this,callback,iwave,nutc,nQSOProgress,nfa,nfb,nfqso, & 33 ndepth,ntrperiod,nexp_decode,ntol,emedelay,lagain,lapcqonly,mycall, & 34 hiscall,iwspr) 35 36 use prog_args 37 use timer_module, only: timer 38 use packjt77 39 use, intrinsic :: iso_c_binding 40 include 'fst4/fst4_params.f90' 41 parameter (MAXCAND=100,MAXWCALLS=100) 42 class(fst4_decoder), intent(inout) :: this 43 procedure(fst4_decode_callback) :: callback 44 character*37 decodes(100) 45 character*37 msg,msgsent 46 character*20 wcalls(MAXWCALLS), wpart 47 character*77 c77 48 character*12 mycall,hiscall 49 character*12 mycall0,hiscall0 50 complex, allocatable :: c2(:) 51 complex, allocatable :: cframe(:) 52 complex, allocatable :: c_bigfft(:) !Complex waveform 53 real llr(240),llrs(240,4) 54 real candidates0(200,5),candidates(200,5) 55 real bitmetrics(320,4) 56 real s4(0:3,NN) 57 real minsync 58 logical lagain,lapcqonly 59 integer itone(NN) 60 integer hmod 61 integer*1 apmask(240),cw(240),hdec(240) 62 integer*1 message101(101),message74(74),message77(77) 63 integer*1 rvec(77) 64 integer apbits(240) 65 integer nappasses(0:5) ! # of decoding passes for QSO states 0-5 66 integer naptypes(0:5,4) ! (nQSOProgress,decoding pass) 67 integer mcq(29),mrrr(19),m73(19),mrr73(19) 68 69 logical badsync,unpk77_success,single_decode 70 logical first,nohiscall,lwspr 71 logical new_callsign,plotspec_exists,wcalls_exists,do_k50_decode 72 logical decdata_exists 73 74 integer*2 iwave(30*60*12000) 75 76 data mcq/0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0/ 77 data mrrr/0,1,1,1,1,1,1,0,1,0,0,1,0,0,1,0,0,0,1/ 78 data m73/0,1,1,1,1,1,1,0,1,0,0,1,0,1,0,0,0,0,1/ 79 data mrr73/0,1,1,1,1,1,1,0,0,1,1,1,0,1,0,1,0,0,1/ 80 data rvec/0,1,0,0,1,0,1,0,0,1,0,1,1,1,1,0,1,0,0,0,1,0,0,1,1,0,1,1,0, & 81 1,0,0,1,0,1,1,0,0,0,0,1,0,0,0,1,0,1,0,0,1,1,1,1,0,0,1,0,1, & 82 0,1,0,1,0,1,1,0,1,1,1,1,1,0,0,0,1,0,1/ 83 data first/.true./,hmod/1/ 84 save first,apbits,nappasses,naptypes,mycall0,hiscall0 85 save wcalls,nwcalls 86 87 this%callback => callback 88 dxcall13=hiscall ! initialize for use in packjt77 89 mycall13=mycall 90 91 if(iwspr.ne.0.and.iwspr.ne.1) return 92 if(lagain) continue ! use lagain to keep compiler happy 93 94 if(first) then 95! read the fst4_calls.txt file 96 inquire(file=trim(data_dir)//'/fst4w_calls.txt',exist=wcalls_exists) 97 if( wcalls_exists ) then 98 open(42,file=trim(data_dir)//'/fst4w_calls.txt',status='unknown') 99 do i=1,MAXWCALLS 100 wcalls(i)='' 101 read(42,fmt='(a)',end=2867) wcalls(i) 102 wcalls(i)=adjustl(wcalls(i)) 103 if(len(trim(wcalls(i))).eq.0) exit 104 enddo 1052867 nwcalls=i-1 106 close(42) 107 endif 108 109 mcq=2*mod(mcq+rvec(1:29),2)-1 110 mrrr=2*mod(mrrr+rvec(59:77),2)-1 111 m73=2*mod(m73+rvec(59:77),2)-1 112 mrr73=2*mod(mrr73+rvec(59:77),2)-1 113 114 nappasses(0)=2 115 nappasses(1)=2 116 nappasses(2)=2 117 nappasses(3)=2 118 nappasses(4)=2 119 nappasses(5)=3 120 121! iaptype 122!------------------------ 123! 1 CQ ??? ??? (29 ap bits) 124! 2 MyCall ??? ??? (29 ap bits) 125! 3 MyCall DxCall ??? (58 ap bits) 126! 4 MyCall DxCall RRR (77 ap bits) 127! 5 MyCall DxCall 73 (77 ap bits) 128! 6 MyCall DxCall RR73 (77 ap bits) 129!******** 130 131 naptypes(0,1:4)=(/1,2,0,0/) ! Tx6 selected (CQ) 132 naptypes(1,1:4)=(/2,3,0,0/) ! Tx1 133 naptypes(2,1:4)=(/2,3,0,0/) ! Tx2 134 naptypes(3,1:4)=(/3,6,0,0/) ! Tx3 135 naptypes(4,1:4)=(/3,6,0,0/) ! Tx4 136 naptypes(5,1:4)=(/3,1,2,0/) ! Tx5 137 138 mycall0='' 139 hiscall0='' 140 first=.false. 141 endif 142 143 l1=index(mycall,char(0)) 144 if(l1.ne.0) mycall(l1:)=" " 145 l1=index(hiscall,char(0)) 146 if(l1.ne.0) hiscall(l1:)=" " 147 if(mycall.ne.mycall0 .or. hiscall.ne.hiscall0) then 148 apbits=0 149 apbits(1)=99 150 apbits(30)=99 151 152 if(len(trim(mycall)) .lt. 3) go to 10 153 154 nohiscall=.false. 155 hiscall0=hiscall 156 if(len(trim(hiscall0)).lt.3) then 157 hiscall0=mycall ! use mycall for dummy hiscall - mycall won't be hashed. 158 nohiscall=.true. 159 endif 160 msg=trim(mycall)//' '//trim(hiscall0)//' RR73' 161 i3=-1 162 n3=-1 163 call pack77(msg,i3,n3,c77) 164 call unpack77(c77,1,msgsent,unpk77_success) 165 if(i3.ne.1 .or. (msg.ne.msgsent) .or. .not.unpk77_success) go to 10 166 read(c77,'(77i1)') message77 167 message77=mod(message77+rvec,2) 168 apbits(1:77)=2*message77-1 169 if(nohiscall) apbits(30)=99 170 17110 continue 172 mycall0=mycall 173 hiscall0=hiscall 174 endif 175!************************************ 176 177 if(nfqso+nqsoprogress.eq.-999) return 178 Keff=91 179 nmax=15*12000 180 if(ntrperiod.eq.15) then 181 nsps=720 182 nmax=15*12000 183 ndown=18 !nss=40,80,160,400 184 nfft1=int(nmax/ndown)*ndown 185 else if(ntrperiod.eq.30) then 186 nsps=1680 187 nmax=30*12000 188 ndown=42 !nss=40,80,168,336 189 nfft1=359856 !nfft2=8568=2^3*3^2*7*17 190 else if(ntrperiod.eq.60) then 191 nsps=3888 192 nmax=60*12000 193 ndown=108 194 nfft1=7500*96 ! nfft2=7500=2^2*3*5^4 195 else if(ntrperiod.eq.120) then 196 nsps=8200 197 nmax=120*12000 198 ndown=205 !nss=40,82,164,328 199 nfft1=7200*200 ! nfft2=7200=2^5*3^2*5^2 200 else if(ntrperiod.eq.300) then 201 nsps=21504 202 nmax=300*12000 203 ndown=512 !nss=42,84,168,336 204 nfft1=7020*512 ! nfft2=7020=2^2*3^3*5*13 205 else if(ntrperiod.eq.900) then 206 nsps=66560 207 nmax=900*12000 208 ndown=1664 !nss=40,80,160,320 209 nfft1=6480*1664 ! nfft2=6480=2^4*3^4*5 210 else if(ntrperiod.eq.1800) then 211 nsps=134400 212 nmax=1800*12000 213 ndown=3360 !nss=40,80,160,320 214 nfft1=6426*3360 ! nfft2=6426=2*3^3*7*17 215 end if 216 nss=nsps/ndown 217 fs=12000.0 !Sample rate 218 fs2=fs/ndown 219 nspsec=nint(fs2) 220 dt=1.0/fs !Sample interval (s) 221 dt2=1.0/fs2 222 tt=nsps*dt !Duration of "itone" symbols (s) 223 baud=1.0/tt 224 sigbw=4.0*baud 225 nfft2=nfft1/ndown !make sure that nfft1 is exactly nfft2*ndown 226 nfft1=nfft2*ndown 227 nh1=nfft1/2 228 229 allocate( c_bigfft(0:nfft1/2) ) 230 allocate( c2(0:nfft2-1) ) 231 allocate( cframe(0:160*nss-1) ) 232 233 jittermax=2 234 do_k50_decode=.false. 235 if(ndepth.eq.3) then 236 nblock=4 237 jittermax=2 238 do_k50_decode=.true. 239 elseif(ndepth.eq.2) then 240 nblock=4 241 jittermax=2 242 do_k50_decode=.false. 243 elseif(ndepth.eq.1) then 244 nblock=4 245 jittermax=0 246 do_k50_decode=.false. 247 endif 248 249! Noise blanker setup 250 ndropmax=1 251 single_decode=iand(nexp_decode,32).ne.0 252 npct=0 253 nb=nexp_decode/256 - 3 254 if(nb.ge.0) npct=nb 255 inb1=20 256 inb2=5 257 if(nb.eq.-1) then 258 inb2=5 !Try NB = 0, 5, 10, 15, 20% 259 else if(nb.eq.-2) then 260 inb2=2 !Try NB = 0, 2, 4,... 20% 261 else if(nb.eq.-3) then 262 inb2=1 !Try NB = 0, 1, 2,... 20% 263 else 264 inb1=0 !Fixed NB value, 0 to 25% 265 endif 266 267 268! nfa,nfb: define the noise-baseline analysis window 269! fa, fb: define the signal search window 270! We usually make nfa<fa and nfb>fb so that noise baseline analysis 271! window extends outside of the [fa,fb] window where we think the signals are. 272! 273 if(iwspr.eq.1) then !FST4W 274 nfa=max(100,nfqso-ntol-100) 275 nfb=min(4800,nfqso+ntol+100) 276 fa=max(100,nint(nfqso+1.5*baud-ntol)) ! signal search window 277 fb=min(4800,nint(nfqso+1.5*baud+ntol)) 278 else if(iwspr.eq.0) then 279 if(single_decode) then 280 fa=max(100,nint(nfa+1.5*baud)) 281 fb=min(4800,nint(nfb+1.5*baud)) 282 ! extend noise fit 100 Hz outside of search window 283 nfa=max(100,nfa-100) 284 nfb=min(4800,nfb+100) 285 else 286 fa=max(100,nint(nfa+1.5*baud)) 287 fb=min(4800,nint(nfb+1.5*baud)) 288 ! extend noise fit 100 Hz outside of search window 289 nfa=max(100,nfa-100) 290 nfb=min(4800,nfb+100) 291 endif 292 endif 293 294 ndecodes=0 295 decodes=' ' 296 new_callsign=.false. 297 do inb=0,inb1,inb2 298 if(nb.lt.0) npct=inb ! we are looping over blanker settings 299 call blanker(iwave,nfft1,ndropmax,npct,c_bigfft) 300 301! The big fft is done once and is used for calculating the smoothed spectrum 302! and also for downconverting/downsampling each candidate. 303 call four2a(c_bigfft,nfft1,1,-1,0) !r2c 304 nhicoh=1 305 nsyncoh=8 306 minsync=1.20 307 if(ntrperiod.eq.15) minsync=1.15 308 309! Get first approximation of candidate frequencies 310 call get_candidates_fst4(c_bigfft,nfft1,nsps,hmod,fs,fa,fb,nfa,nfb, & 311 minsync,ncand,candidates0) 312 isbest=0 313 fc2=0. 314 do icand=1,ncand 315 fc0=candidates0(icand,1) 316 if(iwspr.eq.0 .and. nb.lt.0 .and. npct.ne.0 .and. & 317 abs(fc0-(nfqso+1.5*baud)).gt.ntol) cycle ! blanker loop only near nfqso 318 detmet=candidates0(icand,2) 319 320! Downconvert and downsample a slice of the spectrum centered on the 321! rough estimate of the candidates frequency. 322! Output array c2 is complex baseband sampled at 12000/ndown Sa/sec. 323! The size of the downsampled c2 array is nfft2=nfft1/ndown 324 call timer('dwnsmpl ',0) 325 call fst4_downsample(c_bigfft,nfft1,ndown,fc0,sigbw,c2) 326 call timer('dwnsmpl ',1) 327 328 call timer('sync240 ',0) 329 call fst4_sync_search(c2,nfft2,hmod,fs2,nss,ntrperiod,nsyncoh,emedelay,sbest,fcbest,isbest) 330 call timer('sync240 ',1) 331 332 fc_synced = fc0 + fcbest 333 dt_synced = (isbest-fs2)*dt2 !nominal dt is 1 second so frame starts at sample fs2 334 candidates0(icand,3)=fc_synced 335 candidates0(icand,4)=isbest 336 enddo 337 338! remove duplicate candidates 339 do icand=1,ncand 340 fc=candidates0(icand,3) 341 isbest=nint(candidates0(icand,4)) 342 do ic2=icand+1,ncand 343 fc2=candidates0(ic2,3) 344 isbest2=nint(candidates0(ic2,4)) 345 if(fc2.gt.0.0) then 346 if(abs(fc2-fc).lt.0.10*baud) then ! same frequency 347 if(abs(isbest2-isbest).le.2) then 348 candidates0(ic2,3)=-1 349 endif 350 endif 351 endif 352 enddo 353 enddo 354 ic=0 355 do icand=1,ncand 356 if(candidates0(icand,3).gt.0) then 357 ic=ic+1 358 candidates0(ic,:)=candidates0(icand,:) 359 endif 360 enddo 361 ncand=ic 362 363! If FST4 mode and Single Decode is not checked, then find candidates 364! within 20 Hz of nfqso and put them at the top of the list 365 if(iwspr.eq.0 .and. .not.single_decode) then 366 nclose=count(abs(candidates0(:,3)-(nfqso+1.5*baud)).le.20) 367 k=0 368 do i=1,ncand 369 if(abs(candidates0(i,3)-(nfqso+1.5*baud)).le.20) then 370 k=k+1 371 candidates(k,:)=candidates0(i,:) 372 endif 373 enddo 374 do i=1,ncand 375 if(abs(candidates0(i,3)-(nfqso+1.5*baud)).gt.20) then 376 k=k+1 377 candidates(k,:)=candidates0(i,:) 378 endif 379 enddo 380 else 381 candidates=candidates0 382 endif 383 384 xsnr=0. 385 do icand=1,ncand 386 sync=candidates(icand,2) 387 fc_synced=candidates(icand,3) 388 isbest=nint(candidates(icand,4)) 389 xdt=(isbest-nspsec)/fs2 390 if(ntrperiod.eq.15) xdt=(isbest-real(nspsec)/2.0)/fs2 391 call timer('dwnsmpl ',0) 392 call fst4_downsample(c_bigfft,nfft1,ndown,fc_synced,sigbw,c2) 393 call timer('dwnsmpl ',1) 394 395 do ijitter=0,jittermax 396 if(ijitter.eq.0) ioffset=0 397 if(ijitter.eq.1) ioffset=1 398 if(ijitter.eq.2) ioffset=-1 399 is0=isbest+ioffset 400 iend=is0+160*nss-1 401 if( is0.lt.0 .or. iend.gt.(nfft2-1) ) cycle 402 cframe=c2(is0:iend) 403 bitmetrics=0 404 call timer('bitmetrc',0) 405 call get_fst4_bitmetrics(cframe,nss,bitmetrics, & 406 s4,nsync_qual,badsync) 407 call timer('bitmetrc',1) 408 if(badsync) cycle 409 410 do il=1,4 411 llrs( 1: 60,il)=bitmetrics( 17: 76, il) 412 llrs( 61:120,il)=bitmetrics( 93:152, il) 413 llrs(121:180,il)=bitmetrics(169:228, il) 414 llrs(181:240,il)=bitmetrics(245:304, il) 415 enddo 416 417 apmag=maxval(abs(llrs(:,4)))*1.1 418 ntmax=nblock+nappasses(nQSOProgress) 419 if(lapcqonly) ntmax=nblock+1 420 if(ndepth.eq.1) ntmax=nblock ! no ap for ndepth=1 421 apmask=0 422 423 if(iwspr.eq.1) then ! 50-bit msgs, no ap decoding 424 nblock=4 425 ntmax=nblock 426 endif 427 428 do itry=1,ntmax 429 if(itry.eq.1) llr=llrs(:,1) 430 if(itry.eq.2.and.itry.le.nblock) llr=llrs(:,2) 431 if(itry.eq.3.and.itry.le.nblock) llr=llrs(:,3) 432 if(itry.eq.4.and.itry.le.nblock) llr=llrs(:,4) 433 if(itry.le.nblock) then 434 apmask=0 435 iaptype=0 436 endif 437 438 if(itry.gt.nblock .and. iwspr.eq.0) then ! do ap passes 439 llr=llrs(:,nblock) ! Use largest blocksize as the basis for AP passes 440 iaptype=naptypes(nQSOProgress,itry-nblock) 441 if(lapcqonly) iaptype=1 442 if(iaptype.ge.2 .and. apbits(1).gt.1) cycle ! No, or nonstandard, mycall 443 if(iaptype.ge.3 .and. apbits(30).gt.1) cycle ! No, or nonstandard, dxcall 444 if(iaptype.eq.1) then ! CQ 445 apmask=0 446 apmask(1:29)=1 447 llr(1:29)=apmag*mcq(1:29) 448 endif 449 450 if(iaptype.eq.2) then ! MyCall ??? ??? 451 apmask=0 452 apmask(1:29)=1 453 llr(1:29)=apmag*apbits(1:29) 454 endif 455 456 if(iaptype.eq.3) then ! MyCall DxCall ??? 457 apmask=0 458 apmask(1:58)=1 459 llr(1:58)=apmag*apbits(1:58) 460 endif 461 462 if(iaptype.eq.4 .or. iaptype.eq.5 .or. iaptype .eq.6) then 463 apmask=0 464 apmask(1:77)=1 465 llr(1:58)=apmag*apbits(1:58) 466 if(iaptype.eq.4) llr(59:77)=apmag*mrrr(1:19) 467 if(iaptype.eq.5) llr(59:77)=apmag*m73(1:19) 468 if(iaptype.eq.6) llr(59:77)=apmag*mrr73(1:19) 469 endif 470 endif 471 472 dmin=0.0 473 nharderrors=-1 474 unpk77_success=.false. 475 if(iwspr.eq.0) then 476 maxosd=2 477 Keff=91 478 norder=3 479 call timer('d240_101',0) 480 call decode240_101(llr,Keff,maxosd,norder,apmask,message101, & 481 cw,ntype,nharderrors,dmin) 482 call timer('d240_101',1) 483 if(count(cw.eq.1).eq.0) then 484 nharderrors=-nharderrors 485 cycle 486 endif 487 write(c77,'(77i1)') mod(message101(1:77)+rvec,2) 488 call unpack77(c77,1,msg,unpk77_success) 489 elseif(iwspr.eq.1) then 490! Try decoding with Keff=66 491 maxosd=2 492 call timer('d240_74 ',0) 493 Keff=66 494 norder=3 495 call decode240_74(llr,Keff,maxosd,norder,apmask,message74,cw, & 496 ntype,nharderrors,dmin) 497 call timer('d240_74 ',1) 498 if(nharderrors.lt.0) goto 3465 499 if(count(cw.eq.1).eq.0) then 500 nharderrors=-nharderrors 501 cycle 502 endif 503 write(c77,'(50i1)') message74(1:50) 504 c77(51:77)='000000000000000000000110000' 505 call unpack77(c77,1,msg,unpk77_success) 506 if(unpk77_success .and. do_k50_decode) then 507! If decode was obtained with Keff=66, save call/grid in fst4w_calls.txt if not there already. 508 i1=index(msg,' ') 509 i2=i1+index(msg(i1+1:),' ') 510 wpart=trim(msg(1:i2)) 511! Only save callsigns/grids from type 1 messages 512 if(index(wpart,'/').eq.0 .and. index(wpart,'<').eq.0) then 513 ifound=0 514 do i=1,nwcalls 515 if(index(wcalls(i),wpart).ne.0) ifound=1 516 enddo 517 518 if(ifound.eq.0) then ! This is a new callsign 519 new_callsign=.true. 520 if(nwcalls.lt.MAXWCALLS) then 521 nwcalls=nwcalls+1 522 wcalls(nwcalls)=wpart 523 else 524 wcalls(1:nwcalls-1)=wcalls(2:nwcalls) 525 wcalls(nwcalls)=wpart 526 endif 527 endif 528 endif 529 endif 5303465 continue 531 532! If no decode then try Keff=50 533 iaptype=0 534 if( .not. unpk77_success .and. do_k50_decode ) then 535 maxosd=1 536 call timer('d240_74 ',0) 537 Keff=50 538 norder=4 539 call decode240_74(llr,Keff,maxosd,norder,apmask,message74,cw, & 540 ntype,nharderrors,dmin) 541 call timer('d240_74 ',1) 542 if(count(cw.eq.1).eq.0) then 543 nharderrors=-nharderrors 544 cycle 545 endif 546 write(c77,'(50i1)') message74(1:50) 547 c77(51:77)='000000000000000000000110000' 548 call unpack77(c77,1,msg,unpk77_success) 549! No CRC in this mode, so only accept the decode if call/grid have been seen before 550 if(unpk77_success) then 551 unpk77_success=.false. 552 do i=1,nwcalls 553 if(index(msg,trim(wcalls(i))).gt.0) then 554 unpk77_success=.true. 555 endif 556 enddo 557 endif 558 endif 559 560 endif 561 562 if(nharderrors .ge.0 .and. unpk77_success) then 563 idupe=0 564 do i=1,ndecodes 565 if(decodes(i).eq.msg) idupe=1 566 enddo 567 if(idupe.eq.1) goto 800 568 ndecodes=ndecodes+1 569 decodes(ndecodes)=msg 570 571 if(iwspr.eq.0) then 572 call get_fst4_tones_from_bits(message101,itone,0) 573 else 574 call get_fst4_tones_from_bits(message74,itone,1) 575 endif 576 inquire(file='plotspec',exist=plotspec_exists) 577 fmid=-999.0 578 call timer('dopsprd ',0) 579 if(plotspec_exists) then 580 call dopspread(itone,iwave,nsps,nmax,ndown,hmod, & 581 isbest,fc_synced,fmid,w50) 582 endif 583 call timer('dopsprd ',1) 584 xsig=0 585 do i=1,NN 586 xsig=xsig+s4(itone(i),i) 587 enddo 588 base=candidates(icand,5) 589 arg=600.0*(xsig/base)-1.0 590 if(arg.gt.0.0) then 591 xsnr=10*log10(arg)-35.5-12.5*log10(nsps/8200.0) 592 if(ntrperiod.eq. 15) xsnr=xsnr+2 593 if(ntrperiod.eq. 30) xsnr=xsnr+1 594 if(ntrperiod.eq. 900) xsnr=xsnr+1 595 if(ntrperiod.eq.1800) xsnr=xsnr+2 596 else 597 xsnr=-99.9 598 endif 599 nsnr=nint(xsnr) 600 qual=0.0 601 fsig=fc_synced - 1.5*baud 602 inquire(file=trim(data_dir)//'/decdata',exist=decdata_exists) 603 if(decdata_exists) then 604 hdec=0 605 where(llrs(:,1).ge.0.0) hdec=1 606 nhp=count(hdec.ne.cw) ! # hard errors wrt N=1 soft symbols 607 hd=sum(ieor(hdec,cw)*abs(llrs(:,1))) ! weighted distance wrt N=1 symbols 608 open(21,file=trim(data_dir)//'/fst4_decodes.dat',status='unknown',position='append') 609 write(21,3021) nutc,icand,itry,nsyncoh,iaptype, & 610 ijitter,npct,ntype,Keff,nsync_qual,nharderrors,dmin,nhp,hd, & 611 sync,xsnr,xdt,fsig,w50,trim(msg) 6123021 format(i6.6,i4,6i3,3i4,f6.1,i4,f6.1,f9.2,f6.1,f6.2,f7.1,f7.3,1x,a) 613 close(21) 614 endif 615 call this%callback(nutc,smax1,nsnr,xdt,fsig,msg, & 616 iaptype,qual,ntrperiod,lwspr,fmid,w50) 617! if(iwspr.eq.0 .and. nb.lt.0) go to 900 618 goto 800 619 endif 620 enddo ! metrics 621 enddo ! istart jitter 622800 enddo !candidate list 623 enddo ! noise blanker loop 624 625 if(new_callsign .and. do_k50_decode) then ! re-write the fst4w_calls.txt file 626 open(42,file=trim(data_dir)//'/fst4w_calls.txt',status='unknown') 627 do i=1,nwcalls 628 write(42,'(a20)') trim(wcalls(i)) 629 enddo 630 close(42) 631 endif 632 633 return 634 end subroutine decode 635 636 subroutine sync_fst4(cd0,i0,f0,hmod,ncoh,np,nss,ntr,fs,sync) 637 638! Compute sync power for a complex, downsampled FST4 signal. 639 640 use timer_module, only: timer 641 include 'fst4/fst4_params.f90' 642 complex cd0(0:np-1) 643 complex csync1,csync2,csynct1,csynct2 644 complex ctwk(3200) 645 complex z1,z2,z3,z4,z5 646 integer hmod,isyncword1(0:7),isyncword2(0:7) 647 real f0save 648 common/sync240com/csync1(3200),csync2(3200),csynct1(3200),csynct2(3200) 649 data isyncword1/0,1,3,2,1,0,2,3/ 650 data isyncword2/2,3,1,0,3,2,0,1/ 651 data f0save/-99.9/,nss0/-1/,ntr0/-1/ 652 save twopi,dt,fac,f0save,nss0,ntr0 653 654 p(z1)=(real(z1*fac)**2 + aimag(z1*fac)**2)**0.5 !Compute power 655 656 nz=8*nss 657 call timer('sync240a',0) 658 if(nss.ne.nss0 .or. ntr.ne.ntr0) then 659 twopi=8.0*atan(1.0) 660 dt=1/fs 661 k=1 662 phi1=0.0 663 phi2=0.0 664 do i=0,7 665 dphi1=twopi*hmod*(isyncword1(i)-1.5)/real(nss) 666 dphi2=twopi*hmod*(isyncword2(i)-1.5)/real(nss) 667 do j=1,nss 668 csync1(k)=cmplx(cos(phi1),sin(phi1)) 669 csync2(k)=cmplx(cos(phi2),sin(phi2)) 670 phi1=mod(phi1+dphi1,twopi) 671 phi2=mod(phi2+dphi2,twopi) 672 k=k+1 673 enddo 674 enddo 675 fac=1.0/(8.0*nss) 676 nss0=nss 677 ntr0=ntr 678 f0save=-1.e30 679 endif 680 681 if(f0.ne.f0save) then 682 dphi=twopi*f0*dt 683 phi=0.0 684 do i=1,nz 685 ctwk(i)=cmplx(cos(phi),sin(phi)) 686 phi=mod(phi+dphi,twopi) 687 enddo 688 csynct1(1:nz)=ctwk(1:nz)*csync1(1:nz) 689 csynct2(1:nz)=ctwk(1:nz)*csync2(1:nz) 690 f0save=f0 691 nss0=nss 692 endif 693 call timer('sync240a',1) 694 695 i1=i0 !Costas arrays 696 i2=i0+38*nss 697 i3=i0+76*nss 698 i4=i0+114*nss 699 i5=i0+152*nss 700 701 s1=0.0 702 s2=0.0 703 s3=0.0 704 s4=0.0 705 s5=0.0 706 707 if(ncoh.gt.0) then 708 nsec=8/ncoh 709 do i=1,nsec 710 is=(i-1)*ncoh*nss 711 z1=0 712 if(i1+is.ge.1) then 713 z1=sum(cd0(i1+is:i1+is+ncoh*nss-1)*conjg(csynct1(is+1:is+ncoh*nss))) 714 endif 715 z2=sum(cd0(i2+is:i2+is+ncoh*nss-1)*conjg(csynct2(is+1:is+ncoh*nss))) 716 z3=sum(cd0(i3+is:i3+is+ncoh*nss-1)*conjg(csynct1(is+1:is+ncoh*nss))) 717 z4=sum(cd0(i4+is:i4+is+ncoh*nss-1)*conjg(csynct2(is+1:is+ncoh*nss))) 718 z5=0 719 if(i5+is+ncoh*nss-1.le.np) then 720 z5=sum(cd0(i5+is:i5+is+ncoh*nss-1)*conjg(csynct1(is+1:is+ncoh*nss))) 721 endif 722 s1=s1+abs(z1)/nz 723 s2=s2+abs(z2)/nz 724 s3=s3+abs(z3)/nz 725 s4=s4+abs(z4)/nz 726 s5=s5+abs(z5)/nz 727 enddo 728 else 729 nsub=-ncoh 730 nps=nss/nsub 731 do i=1,8 732 do isub=1,nsub 733 is=(i-1)*nss+(isub-1)*nps 734 z1=0.0 735 if(i1+is.ge.1) then 736 z1=sum(cd0(i1+is:i1+is+nps-1)*conjg(csynct1(is+1:is+nps))) 737 endif 738 z2=sum(cd0(i2+is:i2+is+nps-1)*conjg(csynct2(is+1:is+nps))) 739 z3=sum(cd0(i3+is:i3+is+nps-1)*conjg(csynct1(is+1:is+nps))) 740 z4=sum(cd0(i4+is:i4+is+nps-1)*conjg(csynct2(is+1:is+nps))) 741 z5=0.0 742 if(i5+is+ncoh*nss-1.le.np) then 743 z5=sum(cd0(i5+is:i5+is+nps-1)*conjg(csynct1(is+1:is+nps))) 744 endif 745 s1=s1+abs(z1)/(8*nss) 746 s2=s2+abs(z2)/(8*nss) 747 s3=s3+abs(z3)/(8*nss) 748 s4=s4+abs(z4)/(8*nss) 749 s5=s5+abs(z5)/(8*nss) 750 enddo 751 enddo 752 endif 753 sync = s1+s2+s3+s4+s5 754 return 755 end subroutine sync_fst4 756 757 subroutine fst4_downsample(c_bigfft,nfft1,ndown,f0,sigbw,c1) 758 759! Output: Complex data in c(), sampled at 12000/ndown Hz 760 761 complex c_bigfft(0:nfft1/2) 762 complex c1(0:nfft1/ndown-1) 763 764 df=12000.0/nfft1 765 i0=nint(f0/df) 766 ih=nint( ( f0 + 1.3*sigbw/2.0 )/df) 767 nbw=ih-i0+1 768 c1=0. 769 c1(0)=c_bigfft(i0) 770 nfft2=nfft1/ndown 771 do i=1,nbw 772 if(i0+i.le.nfft1/2) c1(i)=c_bigfft(i0+i) 773 if(i0-i.ge.0) c1(nfft2-i)=c_bigfft(i0-i) 774 enddo 775 c1=c1/nfft2 776 call four2a(c1,nfft2,1,1,1) !c2c FFT back to time domain 777 return 778 779 end subroutine fst4_downsample 780 781 subroutine get_candidates_fst4(c_bigfft,nfft1,nsps,hmod,fs,fa,fb,nfa,nfb, & 782 minsync,ncand,candidates) 783 784 complex c_bigfft(0:nfft1/2) !Full length FFT of raw data 785 integer hmod !Modulation index (submode) 786 integer im(1) !For maxloc 787 real candidates(200,5) !Candidate list 788 real, allocatable :: s(:) !Low resolution power spectrum 789 real, allocatable :: s2(:) !CCF of s() with 4 tones 790 real, allocatable :: sbase(:) !noise baseline estimate 791 real xdb(-3:3) !Model 4-tone CCF peaks 792 real minsync 793 data xdb/0.25,0.50,0.75,1.0,0.75,0.50,0.25/ 794 795 nh1=nfft1/2 796 df1=fs/nfft1 797 baud=fs/nsps !Keying rate 798 df2=baud/2.0 799 nd=df2/df1 !s() sums this many bins of big FFT 800 ndh=nd/2 801 ia=nint(max(100.0,fa)/df2) !Low frequency search limit 802 ib=nint(min(4800.0,fb)/df2) !High frequency limit 803 ina=nint(max(100.0,real(nfa))/df2) !Low freq limit for noise baseline fit 804 inb=nint(min(4800.0,real(nfb))/df2) !High freq limit for noise fit 805 if(ia.lt.ina) ia=ina 806 if(ib.gt.inb) ib=inb 807 808 nnw=nint(48000.*nsps*2./fs) 809 allocate (s(nnw)) 810 s=0. !Compute low-resolution power spectrum 811 do i=ina,inb ! noise analysis window includes signal analysis window 812 j0=nint(i*df2/df1) 813 do j=j0-ndh,j0+ndh 814 s(i)=s(i) + real(c_bigfft(j))**2 + aimag(c_bigfft(j))**2 815 enddo 816 enddo 817 818 ina=max(ina,1+3*hmod) !Don't run off the ends 819 inb=min(inb,nnw-3*hmod) 820 allocate (s2(nnw)) 821 allocate (sbase(nnw)) 822 s2=0. 823 do i=ina,inb !Compute CCF of s() and 4 tones 824 s2(i)=s(i-hmod*3) + s(i-hmod) +s(i+hmod) +s(i+hmod*3) 825 enddo 826 npctile=30 827 call fst4_baseline(s2,nnw,ina+hmod*3,inb-hmod*3,npctile,sbase) 828 if(any(sbase(ina:inb).le.0.0)) return 829 s2(ina:inb)=s2(ina:inb)/sbase(ina:inb) !Normalize wrt noise level 830 831 ncand=0 832 candidates=0 833 if(ia.lt.3) ia=3 834 if(ib.gt.nnw-2) ib=nnw-2 835 836! Find candidates, using the CLEAN algorithm to remove a model of each one 837! from s2() after it has been found. 838 pval=99.99 839 do while(ncand.lt.200) 840 im=maxloc(s2(ia:ib)) 841 iploc=ia+im(1)-1 !Index of CCF peak 842 pval=s2(iploc) !Peak value 843 if(pval.lt.minsync) exit 844 do i=-3,+3 !Remove 0.9 of a model CCF at 845 k=iploc+2*hmod*i !this frequency from s2() 846 if(k.ge.ia .and. k.le.ib) then 847 s2(k)=max(0.,s2(k)-0.9*pval*xdb(i)) 848 endif 849 enddo 850 ncand=ncand+1 851 candidates(ncand,1)=df2*iploc !Candidate frequency 852 candidates(ncand,2)=pval !Rough estimate of SNR 853 candidates(ncand,5)=sbase(iploc) 854 enddo 855 return 856 end subroutine get_candidates_fst4 857 858 subroutine fst4_sync_search(c2,nfft2,hmod,fs2,nss,ntrperiod,nsyncoh,emedelay,sbest,fcbest,isbest) 859 complex c2(0:nfft2-1) 860 integer hmod 861 nspsec=int(fs2) 862 baud=fs2/real(nss) 863 fc1=0.0 864 if(emedelay.lt.0.1) then ! search offsets from 0 s to 2 s 865 is0=1.5*nspsec 866 ishw=1.5*nspsec 867 else ! search plus or minus 1.5 s centered on emedelay 868 is0=nint((emedelay+1.0)*nspsec) 869 ishw=1.5*nspsec 870 endif 871 872 sbest=-1.e30 873 do if=-12,12 874 fc=fc1 + 0.1*baud*if 875 do istart=max(1,is0-ishw),is0+ishw,4*hmod 876 call sync_fst4(c2,istart,fc,hmod,nsyncoh,nfft2,nss, & 877 ntrperiod,fs2,sync) 878 if(sync.gt.sbest) then 879 fcbest=fc 880 isbest=istart 881 sbest=sync 882 endif 883 enddo 884 enddo 885 886 fc1=fcbest 887 is0=isbest 888 ishw=4*hmod 889 isst=1*hmod 890 891 sbest=0.0 892 do if=-7,7 893 fc=fc1 + 0.02*baud*if 894 do istart=max(1,is0-ishw),is0+ishw,isst 895 call sync_fst4(c2,istart,fc,hmod,nsyncoh,nfft2,nss, & 896 ntrperiod,fs2,sync) 897 if(sync.gt.sbest) then 898 fcbest=fc 899 isbest=istart 900 sbest=sync 901 endif 902 enddo 903 enddo 904 end subroutine fst4_sync_search 905 906 subroutine dopspread(itone,iwave,nsps,nmax,ndown,hmod,i0,fc,fmid,w50) 907 908! On "plotspec" special request, compute Doppler spread for a decoded signal 909 910 include 'fst4/fst4_params.f90' 911 complex, allocatable :: cwave(:) !Reconstructed complex signal 912 complex, allocatable :: g(:) !Channel gain, g(t) in QEX paper 913 real,allocatable :: ss(:) !Computed power spectrum of g(t) 914 integer itone(160) !Tones for this message 915 integer*2 iwave(nmax) !Raw Rx data 916 integer hmod !Modulation index 917 data ncall/0/ 918 save ncall 919 920 ncall=ncall+1 921 nfft=2*nmax 922 nwave=max(nmax,(NN+2)*nsps) 923 allocate(cwave(0:nwave-1)) 924 allocate(g(0:nfft-1)) 925 wave=0 926 fsample=12000.0 927 call gen_fst4wave(itone,NN,nsps,nwave,fsample,hmod,fc,1,cwave,wave) 928 cwave=cshift(cwave,-i0*ndown) 929 fac=1.0/32768 930 g(0:nmax-1)=fac*float(iwave)*conjg(cwave(:nmax-1)) 931 g(nmax:)=0. 932 call four2a(g,nfft,1,-1,1) !Forward c2c FFT 933 934 df=12000.0/nfft 935 ia=1.0/df 936 smax=0. 937 do i=-ia,ia !Find smax in +/- 1 Hz around 0. 938 j=i 939 if(j.lt.0) j=i+nfft 940 s=real(g(j))**2 + aimag(g(j))**2 941 smax=max(s,smax) 942 enddo 943 944 ia=10.1/df 945 allocate(ss(-ia:ia)) !Allocate space for +/- 10 Hz 946 sum1=0. 947 sum2=0. 948 nns=0 949 do i=-ia,ia 950 j=i 951 if(j.lt.0) j=i+nfft 952 ss(i)=(real(g(j))**2 + aimag(g(j))**2)/smax 953 f=i*df 954 if(f.ge.-4.0 .and. f.le.-2.0) then 955 sum1=sum1 + ss(i) !Power between -2 and -4 Hz 956 nns=nns+1 957 else if(f.ge.2.0 .and. f.le.4.0) then 958 sum2=sum2 + ss(i) !Power between +2 and +4 Hz 959 endif 960 enddo 961 avg=min(sum1/nns,sum2/nns) !Compute avg from smaller sum 962 963 sum1=0. 964 do i=-ia,ia 965 f=i*df 966 if(abs(f).le.1.0) sum1=sum1 + ss(i)-avg !Power in abs(f) < 1 Hz 967 enddo 968 969 ia=nint(1.0/df) + 1 970 sum2=0.0 971 xi1=-999 972 xi2=-999 973 xi3=-999 974 sum2z=0. 975 do i=-ia,ia !Find freq range that has 50% of signal power 976 sum2=sum2 + ss(i)-avg 977 if(sum2.ge.0.25*sum1 .and. xi1.eq.-999.0) then 978 xi1=i - 1 + (sum2-0.25*sum1)/(sum2-sum2z) 979 endif 980 if(sum2.ge.0.50*sum1 .and. xi2.eq.-999.0) then 981 xi2=i - 1 + (sum2-0.50*sum1)/(sum2-sum2z) 982 endif 983 if(sum2.ge.0.75*sum1) then 984 xi3=i - 1 + (sum2-0.75*sum1)/(sum2-sum2z) 985 exit 986 endif 987 sum2z=sum2 988 enddo 989 xdiff=sqrt(1.0+(xi3-xi1)**2) !Keep small values from fluctuating too widely 990 w50=xdiff*df !Compute Doppler spread 991 fmid=xi2*df !Frequency midpoint of signal powere 992 993 do i=-ia,ia !Save the spectrum for plotting 994 y=ncall-1 995 j=i+nint(xi2) 996 if(abs(j*df).lt.10.0) y=0.99*ss(i+nint(xi2)) + ncall-1 997 write(52,1010) i*df,y 9981010 format(f12.6,f12.6) 999 enddo 1000 1001 return 1002 end subroutine dopspread 1003 1004end module fst4_decode 1005