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