chris/js8call
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Updated to r8541

Browse Source
This commit is contained in:
Jordan Sherer
2018-03-05 14:49:51 -05:00
parent a4fa5b9988
commit a32fe6a4dc
200 changed files with 20394 additions and 4957 deletions
Download Patch File Download Diff File Expand all files Collapse all files
lib/fsk4hf/baseline.f90
-48
View File
@@ -1,48 +0,0 @@
subroutine baseline(s,nfa,nfb,sbase)
! Fit baseline to spectrum (for FT8)
! Input: s(npts) Linear scale in power
! Output: sbase(npts) Baseline
implicit real*8 (a-h,o-z)
real*4 s(1920)
real*4 sbase(1920)
real*4 base
real*8 x(1000),y(1000),a(5)
data nseg/10/,npct/10/
df=12000.0/3840.0 !3.125 Hz
ia=max(1,nint(nfa/df))
ib=nint(nfb/df)
do i=ia,ib
s(i)=10.0*log10(s(i)) !Convert to dB scale
enddo
nterms=5
nlen=(ib-ia+1)/nseg !Length of test segment
i0=(ib-ia+1)/2 !Midpoint
k=0
do n=1,nseg !Loop over all segments
ja=ia + (n-1)*nlen
jb=ja+nlen-1
call pctile(s(ja),nlen,npct,base) !Find lowest npct of points
do i=ja,jb
if(s(i).le.base) then
if (k.lt.1000) k=k+1 !Save all "lower envelope" points
x(k)=i-i0
y(k)=s(i)
endif
enddo
enddo
kz=k
a=0.
call polyfit(x,y,y,kz,nterms,0,a,chisqr) !Fit a low-order polynomial
do i=ia,ib
t=i-i0
sbase(i)=a(1)+t*(a(2)+t*(a(3)+t*(a(4)+t*(a(5))))) + 0.65
! write(51,3051) i*df,s(i),sbase(i)
!3051 format(3f12.3)
enddo
return
end subroutine baseline
lib/fsk4hf/bpdecode174.f90
-401
View File
@@ -1,401 +0,0 @@
subroutine bpdecode174(llr,apmask,maxiterations,decoded,cw,nharderror,iter)
!
! A log-domain belief propagation decoder for the (174,87) code.
!
integer, parameter:: N=174, K=87, M=N-K
integer*1 codeword(N),cw(N),apmask(N)
integer colorder(N)
integer*1 decoded(K)
integer Nm(7,M) ! 5, 6, or 7 bits per check
integer Mn(3,N) ! 3 checks per bit
integer synd(M)
real tov(3,N)
real toc(7,M)
real tanhtoc(7,M)
real zn(N)
real llr(N)
real Tmn
integer nrw(M)
data colorder/ &
0, 1, 2, 3, 30, 4, 5, 6, 7, 8, 9, 10, 11, 32, 12, 40, 13, 14, 15, 16,&
17, 18, 37, 45, 29, 19, 20, 21, 41, 22, 42, 31, 33, 34, 44, 35, 47, 51, 50, 43,&
36, 52, 63, 46, 25, 55, 27, 24, 23, 53, 39, 49, 59, 38, 48, 61, 60, 57, 28, 62,&
56, 58, 65, 66, 26, 70, 64, 69, 68, 67, 74, 71, 54, 76, 72, 75, 78, 77, 80, 79,&
73, 83, 84, 81, 82, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99,&
100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,&
120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,&
140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,&
160,161,162,163,164,165,166,167,168,169,170,171,172,173/
data Mn/ &
1, 25, 69, &
2, 5, 73, &
3, 32, 68, &
4, 51, 61, &
6, 63, 70, &
7, 33, 79, &
8, 50, 86, &
9, 37, 43, &
10, 41, 65, &
11, 14, 64, &
12, 75, 77, &
13, 23, 81, &
15, 16, 82, &
17, 56, 66, &
18, 53, 60, &
19, 31, 52, &
20, 67, 84, &
21, 29, 72, &
22, 24, 44, &
26, 35, 76, &
27, 36, 38, &
28, 40, 42, &
30, 54, 55, &
34, 49, 87, &
39, 57, 58, &
45, 74, 83, &
46, 62, 80, &
47, 48, 85, &
59, 71, 78, &
1, 50, 53, &
2, 47, 84, &
3, 25, 79, &
4, 6, 14, &
5, 7, 80, &
8, 34, 55, &
9, 36, 69, &
10, 43, 83, &
11, 23, 74, &
12, 17, 44, &
13, 57, 76, &
15, 27, 56, &
16, 28, 29, &
18, 19, 59, &
20, 40, 63, &
21, 35, 52, &
22, 54, 64, &
24, 62, 78, &
26, 32, 77, &
30, 72, 85, &
31, 65, 87, &
33, 39, 51, &
37, 48, 75, &
38, 70, 71, &
41, 42, 68, &
45, 67, 86, &
46, 81, 82, &
49, 66, 73, &
58, 60, 66, &
61, 65, 85, &
1, 14, 21, &
2, 13, 59, &
3, 67, 82, &
4, 32, 73, &
5, 36, 54, &
6, 43, 46, &
7, 28, 75, &
8, 33, 71, &
9, 49, 76, &
10, 58, 64, &
11, 48, 68, &
12, 19, 45, &
15, 50, 61, &
16, 22, 26, &
17, 72, 80, &
18, 40, 55, &
20, 35, 51, &
23, 25, 34, &
24, 63, 87, &
27, 39, 74, &
29, 78, 83, &
30, 70, 77, &
31, 69, 84, &
22, 37, 86, &
38, 41, 81, &
42, 44, 57, &
47, 53, 62, &
52, 56, 79, &
60, 75, 81, &
1, 39, 77, &
2, 16, 41, &
3, 31, 54, &
4, 36, 78, &
5, 45, 65, &
6, 57, 85, &
7, 14, 49, &
8, 21, 46, &
9, 15, 72, &
10, 20, 62, &
11, 17, 71, &
12, 34, 47, &
13, 68, 86, &
18, 23, 43, &
19, 64, 73, &
24, 48, 79, &
25, 70, 83, &
26, 80, 87, &
27, 32, 40, &
28, 56, 69, &
29, 63, 66, &
30, 42, 50, &
33, 37, 82, &
35, 60, 74, &
38, 55, 84, &
44, 52, 61, &
51, 53, 72, &
58, 59, 67, &
47, 56, 76, &
1, 19, 37, &
2, 61, 75, &
3, 8, 66, &
4, 60, 84, &
5, 34, 39, &
6, 26, 53, &
7, 32, 57, &
9, 52, 67, &
10, 12, 15, &
11, 51, 69, &
13, 14, 65, &
16, 31, 43, &
17, 20, 36, &
18, 80, 86, &
21, 48, 59, &
22, 40, 46, &
23, 33, 62, &
24, 30, 74, &
25, 42, 64, &
27, 49, 85, &
28, 38, 73, &
29, 44, 81, &
35, 68, 70, &
41, 63, 76, &
45, 49, 71, &
50, 58, 87, &
48, 54, 83, &
13, 55, 79, &
77, 78, 82, &
1, 2, 24, &
3, 6, 75, &
4, 56, 87, &
5, 44, 53, &
7, 50, 83, &
8, 10, 28, &
9, 55, 62, &
11, 29, 67, &
12, 33, 40, &
14, 16, 20, &
15, 35, 73, &
17, 31, 39, &
18, 36, 57, &
19, 46, 76, &
21, 42, 84, &
22, 34, 59, &
23, 26, 61, &
25, 60, 65, &
27, 64, 80, &
30, 37, 66, &
32, 45, 72, &
38, 51, 86, &
41, 77, 79, &
43, 56, 68, &
47, 74, 82, &
40, 52, 78, &
54, 61, 71, &
46, 58, 69/
data Nm/ &
1, 30, 60, 89, 118, 147, 0, &
2, 31, 61, 90, 119, 147, 0, &
3, 32, 62, 91, 120, 148, 0, &
4, 33, 63, 92, 121, 149, 0, &
2, 34, 64, 93, 122, 150, 0, &
5, 33, 65, 94, 123, 148, 0, &
6, 34, 66, 95, 124, 151, 0, &
7, 35, 67, 96, 120, 152, 0, &
8, 36, 68, 97, 125, 153, 0, &
9, 37, 69, 98, 126, 152, 0, &
10, 38, 70, 99, 127, 154, 0, &
11, 39, 71, 100, 126, 155, 0, &
12, 40, 61, 101, 128, 145, 0, &
10, 33, 60, 95, 128, 156, 0, &
13, 41, 72, 97, 126, 157, 0, &
13, 42, 73, 90, 129, 156, 0, &
14, 39, 74, 99, 130, 158, 0, &
15, 43, 75, 102, 131, 159, 0, &
16, 43, 71, 103, 118, 160, 0, &
17, 44, 76, 98, 130, 156, 0, &
18, 45, 60, 96, 132, 161, 0, &
19, 46, 73, 83, 133, 162, 0, &
12, 38, 77, 102, 134, 163, 0, &
19, 47, 78, 104, 135, 147, 0, &
1, 32, 77, 105, 136, 164, 0, &
20, 48, 73, 106, 123, 163, 0, &
21, 41, 79, 107, 137, 165, 0, &
22, 42, 66, 108, 138, 152, 0, &
18, 42, 80, 109, 139, 154, 0, &
23, 49, 81, 110, 135, 166, 0, &
16, 50, 82, 91, 129, 158, 0, &
3, 48, 63, 107, 124, 167, 0, &
6, 51, 67, 111, 134, 155, 0, &
24, 35, 77, 100, 122, 162, 0, &
20, 45, 76, 112, 140, 157, 0, &
21, 36, 64, 92, 130, 159, 0, &
8, 52, 83, 111, 118, 166, 0, &
21, 53, 84, 113, 138, 168, 0, &
25, 51, 79, 89, 122, 158, 0, &
22, 44, 75, 107, 133, 155, 172, &
9, 54, 84, 90, 141, 169, 0, &
22, 54, 85, 110, 136, 161, 0, &
8, 37, 65, 102, 129, 170, 0, &
19, 39, 85, 114, 139, 150, 0, &
26, 55, 71, 93, 142, 167, 0, &
27, 56, 65, 96, 133, 160, 174, &
28, 31, 86, 100, 117, 171, 0, &
28, 52, 70, 104, 132, 144, 0, &
24, 57, 68, 95, 137, 142, 0, &
7, 30, 72, 110, 143, 151, 0, &
4, 51, 76, 115, 127, 168, 0, &
16, 45, 87, 114, 125, 172, 0, &
15, 30, 86, 115, 123, 150, 0, &
23, 46, 64, 91, 144, 173, 0, &
23, 35, 75, 113, 145, 153, 0, &
14, 41, 87, 108, 117, 149, 170, &
25, 40, 85, 94, 124, 159, 0, &
25, 58, 69, 116, 143, 174, 0, &
29, 43, 61, 116, 132, 162, 0, &
15, 58, 88, 112, 121, 164, 0, &
4, 59, 72, 114, 119, 163, 173, &
27, 47, 86, 98, 134, 153, 0, &
5, 44, 78, 109, 141, 0, 0, &
10, 46, 69, 103, 136, 165, 0, &
9, 50, 59, 93, 128, 164, 0, &
14, 57, 58, 109, 120, 166, 0, &
17, 55, 62, 116, 125, 154, 0, &
3, 54, 70, 101, 140, 170, 0, &
1, 36, 82, 108, 127, 174, 0, &
5, 53, 81, 105, 140, 0, 0, &
29, 53, 67, 99, 142, 173, 0, &
18, 49, 74, 97, 115, 167, 0, &
2, 57, 63, 103, 138, 157, 0, &
26, 38, 79, 112, 135, 171, 0, &
11, 52, 66, 88, 119, 148, 0, &
20, 40, 68, 117, 141, 160, 0, &
11, 48, 81, 89, 146, 169, 0, &
29, 47, 80, 92, 146, 172, 0, &
6, 32, 87, 104, 145, 169, 0, &
27, 34, 74, 106, 131, 165, 0, &
12, 56, 84, 88, 139, 0, 0, &
13, 56, 62, 111, 146, 171, 0, &
26, 37, 80, 105, 144, 151, 0, &
17, 31, 82, 113, 121, 161, 0, &
28, 49, 59, 94, 137, 0, 0, &
7, 55, 83, 101, 131, 168, 0, &
24, 50, 78, 106, 143, 149, 0/
data nrw/ &
6,6,6,6,6,6,6,6,6,6, &
6,6,6,6,6,6,6,6,6,6, &
6,6,6,6,6,6,6,6,6,6, &
6,6,6,6,6,6,6,6,6,7, &
6,6,6,6,6,7,6,6,6,6, &
6,6,6,6,6,7,6,6,6,6, &
7,6,5,6,6,6,6,6,6,5, &
6,6,6,6,6,6,6,6,6,6, &
5,6,6,6,5,6,6/
ncw=3
decoded=0
toc=0
tov=0
tanhtoc=0
! initialize messages to checks
do j=1,M
do i=1,nrw(j)
toc(i,j)=llr((Nm(i,j)))
enddo
enddo
ncnt=0
do iter=0,maxiterations
! Update bit log likelihood ratios (tov=0 in iteration 0).
do i=1,N
if( apmask(i) .ne. 1 ) then
zn(i)=llr(i)+sum(tov(1:ncw,i))
else
zn(i)=llr(i)
endif
enddo
! Check to see if we have a codeword (check before we do any iteration).
cw=0
where( zn .gt. 0. ) cw=1
ncheck=0
do i=1,M
synd(i)=sum(cw(Nm(1:nrw(i),i)))
if( mod(synd(i),2) .ne. 0 ) ncheck=ncheck+1
! if( mod(synd(i),2) .ne. 0 ) write(*,*) 'check ',i,' unsatisfied'
enddo
! write(*,*) 'number of unsatisfied parity checks ',ncheck
if( ncheck .eq. 0 ) then ! we have a codeword - reorder the columns and return it
codeword=cw(colorder+1)
decoded=codeword(M+1:N)
nerr=0
do i=1,N
if( (2*cw(i)-1)*llr(i) .lt. 0.0 ) nerr=nerr+1
enddo
nharderror=nerr
return
endif
if( iter.gt.0 ) then ! this code block implements an early stopping criterion
! if( iter.gt.10000 ) then ! this code block implements an early stopping criterion
nd=ncheck-nclast
if( nd .lt. 0 ) then ! # of unsatisfied parity checks decreased
ncnt=0 ! reset counter
else
ncnt=ncnt+1
endif
! write(*,*) iter,ncheck,nd,ncnt
if( ncnt .ge. 5 .and. iter .ge. 10 .and. ncheck .gt. 15) then
nharderror=-1
return
endif
endif
nclast=ncheck
! Send messages from bits to check nodes
do j=1,M
do i=1,nrw(j)
ibj=Nm(i,j)
toc(i,j)=zn(ibj)
do kk=1,ncw ! subtract off what the bit had received from the check
if( Mn(kk,ibj) .eq. j ) then
toc(i,j)=toc(i,j)-tov(kk,ibj)
endif
enddo
enddo
enddo
! send messages from check nodes to variable nodes
do i=1,M
tanhtoc(1:7,i)=tanh(-toc(1:7,i)/2)
enddo
do j=1,N
do i=1,ncw
ichk=Mn(i,j) ! Mn(:,j) are the checks that include bit j
Tmn=product(tanhtoc(1:nrw(ichk),ichk),mask=Nm(1:nrw(ichk),ichk).ne.j)
call platanh(-Tmn,y)
! y=atanh(-Tmn)
tov(i,j)=2*y
enddo
enddo
enddo
nharderror=-1
return
end subroutine bpdecode174
lib/fsk4hf/chkcrc12a.f90
-24
View File
@@ -1,24 +0,0 @@
subroutine chkcrc12a(decoded,nbadcrc)
use crc
integer*1 decoded(87)
integer*1, target:: i1Dec8BitBytes(11)
character*87 cbits
! Write decoded bits into cbits: 75-bit message plus 12-bit CRC
write(cbits,1000) decoded
1000 format(87i1)
read(cbits,1001) i1Dec8BitBytes
1001 format(11b8)
read(cbits,1002) ncrc12 !Received CRC12
1002 format(75x,b12)
i1Dec8BitBytes(10)=iand(i1Dec8BitBytes(10),128+64+32)
i1Dec8BitBytes(11)=0
icrc12=crc12(c_loc(i1Dec8BitBytes),11) !CRC12 computed from 75 msg bits
nbadcrc=1
if(ncrc12.eq.icrc12) nbadcrc=0
return
end subroutine chkcrc12a
lib/fsk4hf/encode174.f90
-50
View File
@@ -1,50 +0,0 @@
subroutine encode174(message,codeword)
! Encode an 101-bit message and return a 174-bit codeword.
! The generator matrix has dimensions (73,101).
! The code is a (174,101) regular ldpc code with column weight 3.
! The code was generated using the PEG algorithm.
! After creating the codeword, the columns are re-ordered according to
! "colorder" to make the codeword compatible with the parity-check matrix
!
include "ldpc_174_87_params.f90"
integer*1 codeword(N)
integer*1 gen(M,K)
integer*1 itmp(N)
integer*1 message(K)
integer*1 pchecks(M)
logical first
data first/.true./
save first,gen
if( first ) then ! fill the generator matrix
gen=0
do i=1,M
do j=1,11
read(g(i)( (j-1)*2+1:(j-1)*2+2 ),"(Z2)") istr
do jj=1, 8
icol=(j-1)*8+jj
if( icol .le. 87 ) then
if( btest(istr,8-jj) ) gen(i,icol)=1
endif
enddo
enddo
enddo
first=.false.
endif
do i=1,M
nsum=0
do j=1,K
nsum=nsum+message(j)*gen(i,j)
enddo
pchecks(i)=mod(nsum,2)
enddo
itmp(1:M)=pchecks
itmp(M+1:N)=message(1:K)
codeword(colorder+1)=itmp(1:N)
return
end subroutine encode174
lib/fsk4hf/extractmessage174.f90
-48
View File
@@ -1,48 +0,0 @@
subroutine extractmessage174(decoded,msgreceived,ncrcflag,recent_calls,nrecent)
use iso_c_binding, only: c_loc,c_size_t
use crc
use packjt
character*22 msgreceived
character*12 call1,call2
character*12 recent_calls(nrecent)
character*87 cbits
integer*1 decoded(87)
integer*1, target:: i1Dec8BitBytes(11)
integer*4 i4Dec6BitWords(12)
! Write decoded bits into cbits: 75-bit message plus 12-bit CRC
write(cbits,1000) decoded
1000 format(87i1)
read(cbits,1001) i1Dec8BitBytes
1001 format(11b8)
read(cbits,1002) ncrc12 !Received CRC12
1002 format(75x,b12)
i1Dec8BitBytes(10)=iand(i1Dec8BitBytes(10),128+64+32)
i1Dec8BitBytes(11)=0
icrc12=crc12(c_loc(i1Dec8BitBytes),11) !CRC12 computed from 75 msg bits
if(ncrc12.eq.icrc12) then
! CRC12 checks out --- unpack 72-bit message
do ibyte=1,12
itmp=0
do ibit=1,6
itmp=ishft(itmp,1)+iand(1,decoded((ibyte-1)*6+ibit))
enddo
i4Dec6BitWords(ibyte)=itmp
enddo
call unpackmsg144(i4Dec6BitWords,msgreceived,call1,call2)
ncrcflag=1
if( call1(1:2) .ne. 'CQ' .and. call1(1:2) .ne. ' ' ) then
call update_recent_calls(call1,recent_calls,nrecent)
endif
if( call2(1:2) .ne. ' ' ) then
call update_recent_calls(call2,recent_calls,nrecent)
endif
else
msgreceived=' '
ncrcflag=-1
endif
return
end subroutine extractmessage174
lib/fsk4hf/ft8_downsample.f90
-42
View File
@@ -1,42 +0,0 @@
subroutine ft8_downsample(dd,newdat,f0,c1)
! Downconvert to complex data sampled at 200 Hz ==> 32 samples/symbol
parameter (NMAX=15*12000,NSPS=1920)
parameter (NFFT1=192000,NFFT2=3200) !192000/60 = 3200
logical newdat
complex c1(0:NFFT2-1)
complex cx(0:NFFT1/2)
real dd(NMAX),x(NFFT1)
equivalence (x,cx)
save cx
if(newdat) then
! Data in dd have changed, recompute the long FFT
x(1:NMAX)=dd
x(NMAX+1:NFFT1)=0. !Zero-pad the x array
call four2a(cx,NFFT1,1,-1,0) !r2c FFT to freq domain
newdat=.false.
endif
df=12000.0/NFFT1
baud=12000.0/NSPS
i0=nint(f0/df)
ft=f0+8.0*baud
it=min(nint(ft/df),NFFT1/2)
fb=f0-1.0*baud
ib=max(1,nint(fb/df))
k=0
c1=0.
do i=ib,it
c1(k)=cx(i)
k=k+1
enddo
c1=cshift(c1,i0-ib)
call four2a(c1,NFFT2,1,1,1) !c2c FFT back to time domain
fac=1.0/sqrt(float(NFFT1)*NFFT2)
c1=fac*c1
return
end subroutine ft8_downsample
lib/fsk4hf/ft8_params.f90
-12
View File
@@ -1,12 +0,0 @@
! LDPC (174,87) code
parameter (KK=87) !Information bits (75 + CRC12)
parameter (ND=58) !Data symbols
parameter (NS=21) !Sync symbols (3 @ Costas 7x7)
parameter (NN=NS+ND) !Total channel symbols (79)
parameter (NSPS=1920) !Samples per symbol at 12000 S/s
parameter (NZ=NSPS*NN) !Samples in full 15 s waveform (151,680)
parameter (NMAX=15*12000) !Samples in iwave (180,000)
parameter (NFFT1=2*NSPS, NH1=NFFT1/2) !Length of FFTs for symbol spectra
parameter (NSTEP=NSPS/4) !Rough time-sync step size
parameter (NHSYM=NMAX/NSTEP-3) !Number of symbol spectra (1/4-sym steps)
parameter (NDOWN=60) !Downsample factor
lib/fsk4hf/ft8apset.f90
-31
View File
@@ -1,31 +0,0 @@
subroutine ft8apset(mycall12,mygrid6,hiscall12,hisgrid6,bcontest,apsym,iaptype)
parameter(NAPM=4,KK=87)
character*12 mycall12,hiscall12
character*22 msg,msgsent
character*6 mycall,hiscall
character*6 mygrid6,hisgrid6
character*4 hisgrid
logical bcontest
integer apsym(KK)
integer*1 msgbits(KK)
integer itone(KK)
mycall=mycall12(1:6)
hiscall=hiscall12(1:6)
hisgrid=hisgrid6(1:4)
if(len_trim(hiscall).eq.0) then
iaptype=1
hiscall="K9AN"
else
iaptype=2
endif
hisgrid=hisgrid6(1:4)
! if(len_trim(hisgrid).eq.0) hisgrid="EN50"
if(index(hisgrid," ").eq.0) hisgrid="EN50"
msg=mycall//' '//hiscall//' '//hisgrid
i3bit=0 ! ### TEMPORARY ??? ###
call genft8(msg,mygrid6,bcontest,i3bit,msgsent,msgbits,itone)
apsym=2*msgbits-1
return
end subroutine ft8apset
lib/fsk4hf/ft8b.f90
-437
View File
@@ -1,437 +0,0 @@
subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid, &
lsubtract,nagain,iaptype,mygrid6,bcontest,sync0,f1,xdt,xbase,apsym, &
nharderrors,dmin,nbadcrc,ipass,iera,message,xsnr)
use timer_module, only: timer
include 'ft8_params.f90'
parameter(NRECENT=10,NP2=2812)
character message*22,msgsent*22
character*12 recent_calls(NRECENT)
character*6 mygrid6
logical bcontest
real a(5)
real s1(0:7,ND),s2(0:7,NN),s1sort(8*ND)
real ps(0:7),psl(0:7)
real bmeta(3*ND),bmetb(3*ND),bmetap(3*ND)
real llr(3*ND),llra(3*ND),llr0(3*ND),llr1(3*ND),llrap(3*ND) !Soft symbols
real dd0(15*12000)
integer*1 decoded(KK),apmask(3*ND),cw(3*ND)
integer*1 msgbits(KK)
integer apsym(KK)
integer mcq(28),mde(28),mrrr(16),m73(16),mrr73(16)
integer itone(NN)
integer indxs1(8*ND)
integer icos7(0:6),ip(1)
integer nappasses(0:5) ! the number of decoding passes to use for each QSO state
integer naptypes(0:5,4) ! (nQSOProgress, decoding pass) maximum of 4 passes for now
complex cd0(3200)
complex ctwk(32)
complex csymb(32)
logical first,newdat,lsubtract,lapon,nagain
equivalence (s1,s1sort)
data icos7/2,5,6,0,4,1,3/
data mcq/1,1,1,1,1,0,1,0,0,0,0,0,1,0,0,0,0,0,1,1,0,0,0,1,1,0,0,1/
data mrrr/0,1,1,1,1,1,1,0,1,1,0,0,1,1,1,1/
data m73/0,1,1,1,1,1,1,0,1,1,0,1,0,0,0,0/
data mde/1,1,1,1,1,1,1,1,0,1,1,0,0,1,0,0,0,0,0,1,1,1,0,1,0,0,0,1/
data mrr73/0,0,0,0,0,0,1,0,0,0,0,1,0,1,0,1/
data first/.true./
save nappasses,naptypes
if(first) then
mcq=2*mcq-1
mde=2*mde-1
mrrr=2*mrrr-1
m73=2*m73-1
mrr73=2*mrr73-1
nappasses(0)=2
nappasses(1)=2
nappasses(2)=2
nappasses(3)=4
nappasses(4)=4
nappasses(5)=3
! iaptype
!------------------------
! 1 CQ ??? ???
! 2 MyCall ??? ???
! 3 MyCall DxCall ???
! 4 MyCall DxCall RRR
! 5 MyCall DxCall 73
! 6 MyCall DxCall RR73
! 7 ??? DxCall ???
naptypes(0,1:4)=(/1,2,0,0/)
naptypes(1,1:4)=(/2,3,0,0/)
naptypes(2,1:4)=(/2,3,0,0/)
naptypes(3,1:4)=(/3,4,5,6/)
naptypes(4,1:4)=(/3,4,5,6/)
naptypes(5,1:4)=(/3,1,2,0/) !?
first=.false.
endif
max_iterations=30
nharderrors=-1
fs2=12000.0/NDOWN
dt2=1.0/fs2
twopi=8.0*atan(1.0)
delfbest=0.
ibest=0
call timer('ft8_down',0)
call ft8_downsample(dd0,newdat,f1,cd0) !Mix f1 to baseband and downsample
call timer('ft8_down',1)
i0=nint((xdt+0.5)*fs2) !Initial guess for start of signal
smax=0.0
do idt=i0-8,i0+8 !Search over +/- one quarter symbol
call sync8d(cd0,idt,ctwk,0,sync)
if(sync.gt.smax) then
smax=sync
ibest=idt
endif
enddo
xdt2=ibest*dt2 !Improved estimate for DT
! Now peak up in frequency
i0=nint(xdt2*fs2)
smax=0.0
do ifr=-5,5 !Search over +/- 2.5 Hz
delf=ifr*0.5
dphi=twopi*delf*dt2
phi=0.0
do i=1,32
ctwk(i)=cmplx(cos(phi),sin(phi))
phi=mod(phi+dphi,twopi)
enddo
call sync8d(cd0,i0,ctwk,1,sync)
if( sync .gt. smax ) then
smax=sync
delfbest=delf
endif
enddo
a=0.0
a(1)=-delfbest
call twkfreq1(cd0,NP2,fs2,a,cd0)
xdt=xdt2
f1=f1+delfbest !Improved estimate of DF
call sync8d(cd0,i0,ctwk,2,sync)
j=0
do k=1,NN
i1=ibest+(k-1)*32
csymb=cmplx(0.0,0.0)
if( i1.ge.1 .and. i1+31 .le. NP2 ) csymb=cd0(i1:i1+31)
call four2a(csymb,32,1,-1,1)
s2(0:7,k)=abs(csymb(1:8))/1e3
enddo
! sync quality check
is1=0
is2=0
is3=0
do k=1,7
ip=maxloc(s2(:,k))
if(icos7(k-1).eq.(ip(1)-1)) is1=is1+1
ip=maxloc(s2(:,k+36))
if(icos7(k-1).eq.(ip(1)-1)) is2=is2+1
ip=maxloc(s2(:,k+72))
if(icos7(k-1).eq.(ip(1)-1)) is3=is3+1
enddo
! hard sync sum - max is 21
nsync=is1+is2+is3
if(nsync .le. 6) then ! bail out
nbadcrc=1
return
endif
j=0
do k=1,NN
if(k.le.7) cycle
if(k.ge.37 .and. k.le.43) cycle
if(k.gt.72) cycle
j=j+1
s1(0:7,j)=s2(0:7,k)
enddo
call indexx(s1sort,8*ND,indxs1)
xmeds1=s1sort(indxs1(nint(0.5*8*ND)))
s1=s1/xmeds1
do j=1,ND
i4=3*j-2
i2=3*j-1
i1=3*j
! Max amplitude
ps=s1(0:7,j)
r1=max(ps(1),ps(3),ps(5),ps(7))-max(ps(0),ps(2),ps(4),ps(6))
r2=max(ps(2),ps(3),ps(6),ps(7))-max(ps(0),ps(1),ps(4),ps(5))
r4=max(ps(4),ps(5),ps(6),ps(7))-max(ps(0),ps(1),ps(2),ps(3))
bmeta(i4)=r4
bmeta(i2)=r2
bmeta(i1)=r1
bmetap(i4)=r4
bmetap(i2)=r2
bmetap(i1)=r1
! Max log metric
psl=log(ps)
r1=max(psl(1),psl(3),psl(5),psl(7))-max(psl(0),psl(2),psl(4),psl(6))
r2=max(psl(2),psl(3),psl(6),psl(7))-max(psl(0),psl(1),psl(4),psl(5))
r4=max(psl(4),psl(5),psl(6),psl(7))-max(psl(0),psl(1),psl(2),psl(3))
bmetb(i4)=r4
bmetb(i2)=r2
bmetb(i1)=r1
! Metric for Cauchy noise
! r1=log(ps(1)**3+ps(3)**3+ps(5)**3+ps(7)**3)- &
! log(ps(0)**3+ps(2)**3+ps(4)**3+ps(6)**3)
! r2=log(ps(2)**3+ps(3)**3+ps(6)**3+ps(7)**3)- &
! log(ps(0)**3+ps(1)**3+ps(4)**3+ps(5)**3)
! r4=log(ps(4)**3+ps(5)**3+ps(6)**3+ps(7)**3)- &
! log(ps(0)**3+ps(1)**3+ps(2)**3+ps(3)**3)
! Metric for AWGN, no fading
! bscale=2.5
! b0=bessi0(bscale*ps(0))
! b1=bessi0(bscale*ps(1))
! b2=bessi0(bscale*ps(2))
! b3=bessi0(bscale*ps(3))
! b4=bessi0(bscale*ps(4))
! b5=bessi0(bscale*ps(5))
! b6=bessi0(bscale*ps(6))
! b7=bessi0(bscale*ps(7))
! r1=log(b1+b3+b5+b7)-log(b0+b2+b4+b6)
! r2=log(b2+b3+b6+b7)-log(b0+b1+b4+b5)
! r4=log(b4+b5+b6+b7)-log(b0+b1+b2+b3)
if(nQSOProgress .eq. 0 .or. nQSOProgress .eq. 5) then
! When bits 88:115 are set as ap bits, bit 115 lives in symbol 39 along
! with no-ap bits 116 and 117. Take care of metrics for bits 116 and 117.
if(j.eq.39) then ! take care of bits that live in symbol 39
if(apsym(28).lt.0) then
bmetap(i2)=max(ps(2),ps(3))-max(ps(0),ps(1))
bmetap(i1)=max(ps(1),ps(3))-max(ps(0),ps(2))
else
bmetap(i2)=max(ps(6),ps(7))-max(ps(4),ps(5))
bmetap(i1)=max(ps(5),ps(7))-max(ps(4),ps(6))
endif
endif
endif
! When bits 116:143 are set as ap bits, bit 115 lives in symbol 39 along
! with ap bits 116 and 117. Take care of metric for bit 115.
! if(j.eq.39) then ! take care of bit 115
! iii=2*(apsym(29)+1)/2 + (apsym(30)+1)/2 ! known values of bits 116 & 117
! if(iii.eq.0) bmetap(i4)=ps(4)-ps(0)
! if(iii.eq.1) bmetap(i4)=ps(5)-ps(1)
! if(iii.eq.2) bmetap(i4)=ps(6)-ps(2)
! if(iii.eq.3) bmetap(i4)=ps(7)-ps(3)
! endif
! bit 144 lives in symbol 48 and will be 1 if it is set as an ap bit.
! take care of metrics for bits 142 and 143
if(j.eq.48) then ! bit 144 is always 1
bmetap(i4)=max(ps(5),ps(7))-max(ps(1),ps(3))
bmetap(i2)=max(ps(3),ps(7))-max(ps(1),ps(5))
endif
! bit 154 lives in symbol 52 and will be 0 if it is set as an ap bit
! take care of metrics for bits 155 and 156
if(j.eq.52) then ! bit 154 will be 0 if it is set as an ap bit.
bmetap(i2)=max(ps(2),ps(3))-max(ps(0),ps(1))
bmetap(i1)=max(ps(1),ps(3))-max(ps(0),ps(2))
endif
enddo
call normalizebmet(bmeta,3*ND)
call normalizebmet(bmetb,3*ND)
call normalizebmet(bmetap,3*ND)
scalefac=2.83
llr0=scalefac*bmeta
llr1=scalefac*bmetb
llra=scalefac*bmetap ! llr's for use with ap
apmag=scalefac*(maxval(abs(bmetap))*1.01)
! pass #
!------------------------------
! 1 regular decoding
! 2 erase 24
! 3 erase 48
! 4 ap pass 1
! 5 ap pass 2
! 6 ap pass 3
! 7 ap pass 4, etc.
if(lapon) then
npasses=4+nappasses(nQSOProgress)
else
npasses=4
endif
do ipass=1,npasses
llr=llr0
if(ipass.eq.2) llr=llr1
if(ipass.eq.3) llr(1:24)=0.
if(ipass.eq.4) llr(1:48)=0.
if(ipass.le.4) then
apmask=0
llrap=llr
iaptype=0
endif
if(ipass .gt. 4) then
iaptype=naptypes(nQSOProgress,ipass-4)
if(iaptype.ge.3 .and. (abs(f1-nfqso).gt.napwid .and. abs(f1-nftx).gt.napwid) ) cycle
if(iaptype.eq.1 .or. iaptype.eq.2 ) then ! AP,???,???
apmask=0
apmask(88:115)=1 ! first 28 bits are AP
apmask(144)=1 ! not free text
llrap=llr
if(iaptype.eq.1) llrap(88:115)=apmag*mcq
if(iaptype.eq.2) llrap(88:115)=apmag*apsym(1:28)
llrap(116:117)=llra(116:117)
llrap(142:143)=llra(142:143)
llrap(144)=-apmag
endif
if(iaptype.eq.3) then ! mycall, dxcall, ???
apmask=0
apmask(88:115)=1 ! mycall
apmask(116:143)=1 ! hiscall
apmask(144)=1 ! not free text
llrap=llr
llrap(88:143)=apmag*apsym(1:56)
llrap(144)=-apmag
endif
if(iaptype.eq.4 .or. iaptype.eq.5 .or. iaptype.eq.6) then
apmask=0
apmask(88:115)=1 ! mycall
apmask(116:143)=1 ! hiscall
apmask(144:159)=1 ! RRR or 73 or RR73
llrap=llr
llrap(88:143)=apmag*apsym(1:56)
if(iaptype.eq.4) llrap(144:159)=apmag*mrrr
if(iaptype.eq.5) llrap(144:159)=apmag*m73
if(iaptype.eq.6) llrap(144:159)=apmag*mrr73
endif
if(iaptype.eq.7) then ! ???, dxcall, ???
apmask=0
apmask(116:143)=1 ! hiscall
apmask(144)=1 ! not free text
llrap=llr
llrap(115)=llra(115)
llrap(116:143)=apmag*apsym(29:56)
llrap(144)=-apmag
endif
endif
cw=0
call timer('bpd174 ',0)
call bpdecode174(llrap,apmask,max_iterations,decoded,cw,nharderrors, &
niterations)
call timer('bpd174 ',1)
dmin=0.0
if(ndepth.eq.3 .and. nharderrors.lt.0) then
ndeep=3
if(abs(nfqso-f1).le.napwid .or. abs(nftx-f1).le.napwid) then
if((ipass.eq.3 .or. ipass.eq.4) .and. .not.nagain) then
ndeep=3
else
ndeep=4
endif
endif
if(nagain) ndeep=5
call timer('osd174 ',0)
call osd174(llrap,apmask,ndeep,decoded,cw,nharderrors,dmin)
call timer('osd174 ',1)
endif
nbadcrc=1
message=' '
xsnr=-99.0
if(count(cw.eq.0).eq.174) cycle !Reject the all-zero codeword
!### if(any(decoded(73:75).ne.0)) cycle !Reject if any of the 3 extra bits is nonzero
if(nharderrors.ge.0 .and. nharderrors+dmin.lt.60.0 .and. &
.not.(sync.lt.2.0 .and. nharderrors.gt.35) .and. &
.not.(ipass.gt.2 .and. nharderrors.gt.39) .and. &
.not.(ipass.eq.4 .and. nharderrors.gt.30) &
) then
call chkcrc12a(decoded,nbadcrc)
else
nharderrors=-1
cycle
endif
!###
i3bit=4*decoded(73) + 2*decoded(74) + decoded(75)
iFreeText=decoded(57)
! if(nbadcrc.eq.0) write(*,3001) nharderrors,nbadcrc,i3bit
!3001 format('A',3i5)
!###
if(nbadcrc.eq.0) then
call extractmessage174(decoded,message,ncrcflag,recent_calls,nrecent)
call genft8(message,mygrid6,bcontest,i3bit,msgsent,msgbits,itone)
if(i3bit.eq.1 .and. iFreeText.eq.0) message(21:21)='1'
if(i3bit.eq.2 .and. iFreeText.eq.0) message(21:21)='2'
if(lsubtract) call subtractft8(dd0,itone,f1,xdt2)
xsig=0.0
xnoi=0.0
do i=1,79
xsig=xsig+s2(itone(i),i)**2
ios=mod(itone(i)+4,7)
xnoi=xnoi+s2(ios,i)**2
enddo
xsnr=0.001
if(xnoi.gt.0 .and. xnoi.lt.xsig) xsnr=xsig/xnoi-1.0
xsnr=10.0*log10(xsnr)-27.0
xsnr2=db(xsig/xbase - 1.0) - 32.0
! write(52,3052) f1,xdt,xsig,xnoi,xbase,xsnr,xsnr2
!3052 format(7f10.2)
if(.not.nagain) xsnr=xsnr2
if(xsnr .lt. -24.0) xsnr=-24.0
return
endif
enddo
return
end subroutine ft8b
subroutine normalizebmet(bmet,n)
real bmet(n)
bmetav=sum(bmet)/real(n)
bmet2av=sum(bmet*bmet)/real(n)
var=bmet2av-bmetav*bmetav
if( var .gt. 0.0 ) then
bmetsig=sqrt(var)
else
bmetsig=sqrt(bmet2av)
endif
bmet=bmet/bmetsig
return
end subroutine normalizebmet
function bessi0(x)
! From Numerical Recipes
real bessi0,x
double precision p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y
save p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9
data p1,p2,p3,p4,p5,p6,p7/1.0d0,3.5156229d0,3.0899424d0,1.2067492d0, &
0.2659732d0,0.360768d-1,0.45813d-2/
data q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,0.1328592d-1, &
0.225319d-2,-0.157565d-2,0.916281d-2,-0.2057706d-1, &
0.2635537d-1,-0.1647633d-1,0.392377d-2/
if (abs(x).lt.3.75) then
y=(x/3.75)**2
bessi0=p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))
else
ax=abs(x)
y=3.75/ax
bessi0=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4 &
+y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))))
endif
return
end function bessi0
lib/fsk4hf/ft8d.f90
-63
View File
@@ -1,63 +0,0 @@
program ft8d
! Decode FT8 data read from *.wav files.
include 'ft8_params.f90'
character*12 arg
character infile*80,datetime*13,message*22
real s(NH1,NHSYM)
real candidate(3,100)
integer ihdr(11)
integer*2 iwave(NMAX) !Generated full-length waveform
real dd(NMAX)
nargs=iargc()
if(nargs.lt.3) then
print*,'Usage: ft8d MaxIt Norder file1 [file2 ...]'
print*,'Example ft8d 40 2 *.wav'
go to 999
endif
call getarg(1,arg)
read(arg,*) max_iterations
call getarg(2,arg)
read(arg,*) norder
nfiles=nargs-2
twopi=8.0*atan(1.0)
fs=12000.0 !Sample rate
dt=1.0/fs !Sample interval (s)
tt=NSPS*dt !Duration of "itone" symbols (s)
ts=2*NSPS*dt !Duration of OQPSK symbols (s)
baud=1.0/tt !Keying rate (baud)
txt=NZ*dt !Transmission length (s)
nfa=100.0
nfb=3000.0
nfqso=1500.0
do ifile=1,nfiles
call getarg(ifile+2,infile)
open(10,file=infile,status='old',access='stream')
read(10,end=999) ihdr,iwave
close(10)
j2=index(infile,'.wav')
read(infile(j2-6:j2-1),*) nutc
datetime=infile(j2-13:j2-1)
call sync8(iwave,nfa,nfb,nfqso,s,candidate,ncand)
syncmin=2.0
dd=iwave
do icand=1,ncand
sync=candidate(3,icand)
if( sync.lt.syncmin) cycle
f1=candidate(1,icand)
xdt=candidate(2,icand)
nsnr=min(99,nint(10.0*log10(sync)-25.5))
call ft8b(dd,nfqso,f1,xdt,nharderrors,dmin,nbadcrc,message,xsnr)
nsnr=xsnr
xdt=xdt-0.6
write(*,1110) datetime,0,nsnr,xdt,f1,message,nharderrors,dmin
1110 format(a13,2i4,f6.2,f7.1,' ~ ',a22,i6,f7.1)
enddo
enddo ! ifile loop
999 end program ft8d
lib/fsk4hf/ft8sim.f90
-126
View File
@@ -1,126 +0,0 @@
program ft8sim
! Generate simulated data for a 15-second HF/6m mode using 8-FSK.
! Output is saved to a *.wav file.
use wavhdr
include 'ft8_params.f90' !Set various constants
type(hdr) h !Header for .wav file
character arg*12,fname*17,sorm*1
character msg*22,msgsent*22
character*6 mygrid6
logical bcontest
complex c0(0:NMAX-1)
complex c(0:NMAX-1)
integer itone(NN)
integer*1 msgbits(KK)
integer*2 iwave(NMAX) !Generated full-length waveform
data mygrid6/'EM48 '/
! Get command-line argument(s)
nargs=iargc()
if(nargs.ne.8) then
print*,'Usage: ft8sim "message" s|m f0 DT fdop del nfiles snr'
print*,'Example: ft8sim "K1ABC W9XYZ EN37" m 1500.0 0.0 0.1 1.0 10 -18'
print*,'s|m: "s" for single signal at 1500 Hz, "m" for 25 signals'
print*,'f0 is ignored when sorm = m'
print*,'Make nfiles negative to invoke 72-bit contest mode.'
go to 999
endif
call getarg(1,msg) !Message to be transmitted
call getarg(2,sorm) !s for single signal, m for multiple sigs
if(sorm.eq."s") then
print*,"Generating single signal at 1500 Hz."
nsig=1
elseif( sorm.eq."m") then
print*,"Generating 25 signals per file."
nsig=25
else
print*,"sorm parameter must be s (single) or m (multiple)."
goto 999
endif
call getarg(3,arg)
read(arg,*) f0 !Frequency (only used for single-signal)
call getarg(4,arg)
read(arg,*) xdt !Time offset from nominal (s)
call getarg(5,arg)
read(arg,*) fspread !Watterson frequency spread (Hz)
call getarg(6,arg)
read(arg,*) delay !Watterson delay (ms)
call getarg(7,arg)
read(arg,*) nfiles !Number of files
call getarg(8,arg)
read(arg,*) snrdb !SNR_2500
bcontest=nfiles.lt.0
nfiles=abs(nfiles)
twopi=8.0*atan(1.0)
fs=12000.0 !Sample rate (Hz)
dt=1.0/fs !Sample interval (s)
tt=NSPS*dt !Duration of symbols (s)
baud=1.0/tt !Keying rate (baud)
bw=8*baud !Occupied bandwidth (Hz)
txt=NZ*dt !Transmission length (s)
bandwidth_ratio=2500.0/(fs/2.0)
sig=sqrt(2*bandwidth_ratio) * 10.0**(0.05*snrdb)
if(snrdb.gt.90.0) sig=1.0
txt=NN*NSPS/12000.0
i3bit=0 ! ### TEMPORARY ??? ###
! Source-encode, then get itone()
call genft8(msg,mygrid6,bcontest,i3bit,msgsent,msgbits,itone)
write(*,1000) f0,xdt,txt,snrdb,bw,msgsent
1000 format('f0:',f9.3,' DT:',f6.2,' TxT:',f6.1,' SNR:',f6.1, &
' BW:',f4.1,2x,a22)
write(*,'(28i1,1x,28i1)') msgbits(1:56)
write(*,'(16i1)') msgbits(57:72)
write(*,'(3i1)') msgbits(73:75)
write(*,'(12i1)') msgbits(76:87)
! call sgran()
do ifile=1,nfiles
c=0.
do isig=1,nsig
c0=0.
if(nsig.eq.25) then
f0=(isig+2)*100.0
endif
k=-1 + nint((xdt+0.5+0.01*gran())/dt)
! k=-1 + nint((xdt+0.5)/dt)
phi=0.0
do j=1,NN !Generate complex waveform
dphi=twopi*(f0+itone(j)*baud)*dt
do i=1,NSPS
k=k+1
phi=mod(phi+dphi,twopi)
if(k.ge.0 .and. k.lt.NMAX) c0(k)=cmplx(cos(phi),sin(phi))
enddo
enddo
if(fspread.ne.0.0 .or. delay.ne.0.0) call watterson(c0,NMAX,fs,delay,fspread)
c=c+sig*c0
enddo
if(snrdb.lt.90) then
do i=0,NMAX-1 !Add gaussian noise at specified SNR
xnoise=gran()
ynoise=gran()
c(i)=c(i) + cmplx(xnoise,ynoise)
enddo
endif
fac=32767.0
rms=100.0
if(snrdb.ge.90.0) iwave(1:NMAX)=nint(fac*real(c))
if(snrdb.lt.90.0) iwave(1:NMAX)=nint(rms*real(c))
h=default_header(12000,NMAX)
write(fname,1102) ifile
1102 format('000000_',i6.6,'.wav')
open(10,file=fname,status='unknown',access='stream')
write(10) h,iwave !Save to *.wav file
close(10)
write(*,1110) ifile,xdt,f0,snrdb,fname
1110 format(i4,f7.2,f8.2,f7.1,2x,a17)
enddo
999 end program ft8sim
lib/fsk4hf/genft8.f90
-56
View File
@@ -1,56 +0,0 @@
subroutine genft8(msg,mygrid,bcontest,i3bit,msgsent,msgbits,itone)
! Encode an FT8 message, producing array itone().
use crc
use packjt
include 'ft8_params.f90'
character*22 msg,msgsent
character*6 mygrid
character*87 cbits
logical bcontest
integer*4 i4Msg6BitWords(12) !72-bit message as 6-bit words
integer*1 msgbits(KK),codeword(3*ND)
integer*1, target:: i1Msg8BitBytes(11)
integer itone(NN)
integer icos7(0:6)
data icos7/2,5,6,0,4,1,3/ !Costas 7x7 tone pattern
call packmsg(msg,i4Msg6BitWords,itype,bcontest) !Pack into 12 6-bit bytes
call unpackmsg(i4Msg6BitWords,msgsent,bcontest,mygrid) !Unpack to get msgsent
write(cbits,1000) i4Msg6BitWords,32*i3bit
1000 format(12b6.6,b8.8)
read(cbits,1001) i1Msg8BitBytes(1:10)
1001 format(10b8)
i1Msg8BitBytes(10)=iand(i1Msg8BitBytes(10),128+64+32)
i1Msg8BitBytes(11)=0
icrc12=crc12(c_loc(i1Msg8BitBytes),11)
! For reference, here's how to check the CRC
! i1Msg8BitBytes(10)=icrc12/256
! i1Msg8BitBytes(11)=iand (icrc12,255)
! checksumok = crc12_check(c_loc (i1Msg8BitBytes), 11)
! if( checksumok ) write(*,*) 'Good checksum'
write(cbits,1003) i4Msg6BitWords,i3bit,icrc12
1003 format(12b6.6,b3.3,b12.12)
read(cbits,1004) msgbits
1004 format(87i1)
call encode174(msgbits,codeword) !Encode the test message
! Message structure: S7 D29 S7 D29 S7
itone(1:7)=icos7
itone(36+1:36+7)=icos7
itone(NN-6:NN)=icos7
k=7
do j=1,ND
i=3*j -2
k=k+1
if(j.eq.30) k=k+7
itone(k)=codeword(i)*4 + codeword(i+1)*2 + codeword(i+2)
enddo
return
end subroutine genft8
lib/fsk4hf/genft8refsig.f90
-22
View File
@@ -1,22 +0,0 @@
subroutine genft8refsig(itone,cref,f0)
complex cref(79*1920)
integer itone(79)
real*8 twopi,phi,dphi,dt,xnsps
data twopi/0.d0/
save twopi
if( twopi .lt. 0.1 ) twopi=8.d0*atan(1.d0)
xnsps=1920.d0
dt=1.d0/12000.d0
phi=0.d0
k=1
do i=1,79
dphi=twopi*(f0*dt+itone(i)/xnsps)
do is=1,1920
cref(k)=cmplx(cos(phi),sin(phi))
phi=mod(phi+dphi,twopi)
k=k+1
enddo
enddo
return
end subroutine genft8refsig
lib/fsk4hf/ldpc_174_87_params.f90
-102
View File
@@ -1,102 +0,0 @@
integer, parameter:: N=174, K=87, M=N-K
character*22 g(87)
integer colorder(N)
data g/ & !parity generator matrix for (174,87) code
"23bba830e23b6b6f50982e", &
"1f8e55da218c5df3309052", &
"ca7b3217cd92bd59a5ae20", &
"56f78313537d0f4382964e", &
"29c29dba9c545e267762fe", &
"6be396b5e2e819e373340c", &
"293548a138858328af4210", &
"cb6c6afcdc28bb3f7c6e86", &
"3f2a86f5c5bd225c961150", &
"849dd2d63673481860f62c", &
"56cdaec6e7ae14b43feeee", &
"04ef5cfa3766ba778f45a4", &
"c525ae4bd4f627320a3974", &
"fe37802941d66dde02b99c", &
"41fd9520b2e4abeb2f989c", &
"40907b01280f03c0323946", &
"7fb36c24085a34d8c1dbc4", &
"40fc3e44bb7d2bb2756e44", &
"d38ab0a1d2e52a8ec3bc76", &
"3d0f929ef3949bd84d4734", &
"45d3814f504064f80549ae", &
"f14dbf263825d0bd04b05e", &
"f08a91fb2e1f78290619a8", &
"7a8dec79a51e8ac5388022", &
"ca4186dd44c3121565cf5c", &
"db714f8f64e8ac7af1a76e", &
"8d0274de71e7c1a8055eb0", &
"51f81573dd4049b082de14", &
"d037db825175d851f3af00", &
"d8f937f31822e57c562370", &
"1bf1490607c54032660ede", &
"1616d78018d0b4745ca0f2", &
"a9fa8e50bcb032c85e3304", &
"83f640f1a48a8ebc0443ea", &
"eca9afa0f6b01d92305edc", &
"3776af54ccfbae916afde6", &
"6abb212d9739dfc02580f2", &
"05209a0abb530b9e7e34b0", &
"612f63acc025b6ab476f7c", &
"0af7723161ec223080be86", &
"a8fc906976c35669e79ce0", &
"45b7ab6242b77474d9f11a", &
"b274db8abd3c6f396ea356", &
"9059dfa2bb20ef7ef73ad4", &
"3d188ea477f6fa41317a4e", &
"8d9071b7e7a6a2eed6965e", &
"a377253773ea678367c3f6", &
"ecbd7c73b9cd34c3720c8a", &
"b6537f417e61d1a7085336", &
"6c280d2a0523d9c4bc5946", &
"d36d662a69ae24b74dcbd8", &
"d747bfc5fd65ef70fbd9bc", &
"a9fa2eefa6f8796a355772", &
"cc9da55fe046d0cb3a770c", &
"f6ad4824b87c80ebfce466", &
"cc6de59755420925f90ed2", &
"164cc861bdd803c547f2ac", &
"c0fc3ec4fb7d2bb2756644", &
"0dbd816fba1543f721dc72", &
"a0c0033a52ab6299802fd2", &
"bf4f56e073271f6ab4bf80", &
"57da6d13cb96a7689b2790", &
"81cfc6f18c35b1e1f17114", &
"481a2a0df8a23583f82d6c", &
"1ac4672b549cd6dba79bcc", &
"c87af9a5d5206abca532a8", &
"97d4169cb33e7435718d90", &
"a6573f3dc8b16c9d19f746", &
"2c4142bf42b01e71076acc", &
"081c29a10d468ccdbcecb6", &
"5b0f7742bca86b8012609a", &
"012dee2198eba82b19a1da", &
"f1627701a2d692fd9449e6", &
"35ad3fb0faeb5f1b0c30dc", &
"b1ca4ea2e3d173bad4379c", &
"37d8e0af9258b9e8c5f9b2", &
"cd921fdf59e882683763f6", &
"6114e08483043fd3f38a8a", &
"2e547dd7a05f6597aac516", &
"95e45ecd0135aca9d6e6ae", &
"b33ec97be83ce413f9acc8", &
"c8b5dffc335095dcdcaf2a", &
"3dd01a59d86310743ec752", &
"14cd0f642fc0c5fe3a65ca", &
"3a0a1dfd7eee29c2e827e0", &
"8abdb889efbe39a510a118", &
"3f231f212055371cf3e2a2"/
data colorder/ &
0, 1, 2, 3, 30, 4, 5, 6, 7, 8, 9, 10, 11, 32, 12, 40, 13, 14, 15, 16,&
17, 18, 37, 45, 29, 19, 20, 21, 41, 22, 42, 31, 33, 34, 44, 35, 47, 51, 50, 43,&
36, 52, 63, 46, 25, 55, 27, 24, 23, 53, 39, 49, 59, 38, 48, 61, 60, 57, 28, 62,&
56, 58, 65, 66, 26, 70, 64, 69, 68, 67, 74, 71, 54, 76, 72, 75, 78, 77, 80, 79,&
73, 83, 84, 81, 82, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99,&
100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,&
120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,&
140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,&
160,161,162,163,164,165,166,167,168,169,170,171,172,173/
lib/fsk4hf/ldpcsim174.f90
-238
View File
@@ -1,238 +0,0 @@
program ldpcsim174
! End to end test of the (174,75)/crc12 encoder and decoder.
use crc
use packjt
parameter(NRECENT=10)
character*12 recent_calls(NRECENT)
character*22 msg,msgsent,msgreceived
character*8 arg
integer*1, allocatable :: codeword(:), decoded(:), message(:)
integer*1, target:: i1Msg8BitBytes(11)
integer*1 msgbits(87)
integer*1 apmask(174), cw(174)
integer*2 checksum
integer*4 i4Msg6BitWords(13)
integer colorder(174)
integer nerrtot(174),nerrdec(174),nmpcbad(87)
logical checksumok,fsk,bpsk
real*8, allocatable :: rxdata(:)
real, allocatable :: llr(:)
data colorder/ &
0, 1, 2, 3, 30, 4, 5, 6, 7, 8, 9, 10, 11, 32, 12, 40, 13, 14, 15, 16,&
17, 18, 37, 45, 29, 19, 20, 21, 41, 22, 42, 31, 33, 34, 44, 35, 47, 51, 50, 43,&
36, 52, 63, 46, 25, 55, 27, 24, 23, 53, 39, 49, 59, 38, 48, 61, 60, 57, 28, 62,&
56, 58, 65, 66, 26, 70, 64, 69, 68, 67, 74, 71, 54, 76, 72, 75, 78, 77, 80, 79,&
73, 83, 84, 81, 82, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99,&
100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,&
120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,&
140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,&
160,161,162,163,164,165,166,167,168,169,170,171,172,173/
do i=1,NRECENT
recent_calls(i)=' '
enddo
nerrtot=0
nerrdec=0
nmpcbad=0 ! Used to collect the number of errors in the message+crc part of the codeword
nargs=iargc()
if(nargs.ne.4) then
print*,'Usage: ldpcsim niter ndepth #trials s '
print*,'eg: ldpcsim 10 2 1000 0.84'
print*,'belief propagation iterations: niter, ordered-statistics depth: ndepth'
print*,'If s is negative, then value is ignored and sigma is calculated from SNR.'
return
endif
call getarg(1,arg)
read(arg,*) max_iterations
call getarg(2,arg)
read(arg,*) ndepth
call getarg(3,arg)
read(arg,*) ntrials
call getarg(4,arg)
read(arg,*) s
fsk=.false.
bpsk=.true.
! don't count crc bits as data bits
N=174
K=87
! scale Eb/No for a (174,87) code
rate=real(K)/real(N)
write(*,*) "rate: ",rate
write(*,*) "niter= ",max_iterations," s= ",s
allocate ( codeword(N), decoded(K), message(K) )
allocate ( rxdata(N), llr(N) )
msg="K1JT K9AN EN50"
! msg="G4WJS K9AN EN50"
call packmsg(msg,i4Msg6BitWords,itype,.false.) !Pack into 12 6-bit bytes
call unpackmsg(i4Msg6BitWords,msgsent,.false.,'') !Unpack to get msgsent
write(*,*) "message sent ",msgsent
i4=0
ik=0
im=0
do i=1,12
nn=i4Msg6BitWords(i)
do j=1, 6
ik=ik+1
i4=i4+i4+iand(1,ishft(nn,j-6))
i4=iand(i4,255)
if(ik.eq.8) then
im=im+1
! if(i4.gt.127) i4=i4-256
i1Msg8BitBytes(im)=i4
ik=0
endif
enddo
enddo
i1Msg8BitBytes(10:11)=0
checksum = crc12 (c_loc (i1Msg8BitBytes), 11)
! For reference, the next 3 lines show how to check the CRC
i1Msg8BitBytes(10)=checksum/256
i1Msg8BitBytes(11)=iand (checksum,255)
checksumok = crc12_check(c_loc (i1Msg8BitBytes), 11)
if( checksumok ) write(*,*) 'Good checksum'
! K=87, For now:
! msgbits(1:72) JT message bits
! msgbits(73:75) 3 free message bits (set to 0)
! msgbits(76:87) CRC12
mbit=0
do i=1, 9
i1=i1Msg8BitBytes(i)
do ibit=1,8
mbit=mbit+1
msgbits(mbit)=iand(1,ishft(i1,ibit-8))
enddo
enddo
msgbits(73:75)=0 ! the three extra message bits go here
i1=i1Msg8BitBytes(10) ! First 4 bits of crc12 are LSB of this byte
do ibit=1,4
msgbits(75+ibit)=iand(1,ishft(i1,ibit-4))
enddo
i1=i1Msg8BitBytes(11) ! Now shift in last 8 bits of the CRC
do ibit=1,8
msgbits(79+ibit)=iand(1,ishft(i1,ibit-8))
enddo
write(*,*) 'message'
write(*,'(11(8i1,1x))') msgbits
call encode174(msgbits,codeword)
call init_random_seed()
! call sgran()
write(*,*) 'codeword'
write(*,'(22(8i1,1x))') codeword
write(*,*) "Es/N0 SNR2500 ngood nundetected nbadcrc sigma"
do idb = 20,-10,-1
!do idb = -3,-3,-1
db=idb/2.0-1.0
sigma=1/sqrt( 2*(10**(db/10.0)) )
ngood=0
nue=0
nbadcrc=0
nberr=0
do itrial=1, ntrials
! Create a realization of a noisy received word
do i=1,N
if( bpsk ) then
rxdata(i) = 2.0*codeword(i)-1.0 + sigma*gran()
elseif( fsk ) then
if( codeword(i) .eq. 1 ) then
r1=(1.0 + sigma*gran())**2 + (sigma*gran())**2
r2=(sigma*gran())**2 + (sigma*gran())**2
elseif( codeword(i) .eq. 0 ) then
r2=(1.0 + sigma*gran())**2 + (sigma*gran())**2
r1=(sigma*gran())**2 + (sigma*gran())**2
endif
! rxdata(i)=0.35*(sqrt(r1)-sqrt(r2))
! rxdata(i)=0.35*(exp(r1)-exp(r2))
rxdata(i)=0.12*(log(r1)-log(r2))
endif
enddo
nerr=0
do i=1,N
if( rxdata(i)*(2*codeword(i)-1.0) .lt. 0 ) nerr=nerr+1
enddo
nerrtot(nerr)=nerrtot(nerr)+1
nberr=nberr+nerr
! Correct signal normalization is important for this decoder.
rxav=sum(rxdata)/N
rx2av=sum(rxdata*rxdata)/N
rxsig=sqrt(rx2av-rxav*rxav)
rxdata=rxdata/rxsig
! To match the metric to the channel, s should be set to the noise standard deviation.
! For now, set s to the value that optimizes decode probability near threshold.
! The s parameter can be tuned to trade a few tenth's dB of threshold for an order of
! magnitude in UER
if( s .lt. 0 ) then
ss=sigma
else
ss=s
endif
llr=2.0*rxdata/(ss*ss)
nap=0 ! number of AP bits
llr(colorder(174-87+1:174-87+nap)+1)=5*(2.0*msgbits(1:nap)-1.0)
apmask=0
apmask(colorder(174-87+1:174-87+nap)+1)=1
! max_iterations is max number of belief propagation iterations
call bpdecode174(llr, apmask, max_iterations, decoded, cw, nharderrors,niterations)
if( ndepth .ge. 0 .and. nharderrors .lt. 0 ) call osd174(llr, apmask, ndepth, decoded, cw, nharderrors, dmin)
! If the decoder finds a valid codeword, nharderrors will be .ge. 0.
if( nharderrors .ge. 0 ) then
call extractmessage174(decoded,msgreceived,ncrcflag,recent_calls,nrecent)
if( ncrcflag .ne. 1 ) then
nbadcrc=nbadcrc+1
endif
nueflag=0
nerrmpc=0
do i=1,K ! find number of errors in message+crc part of codeword
if( msgbits(i) .ne. decoded(i) ) then
nueflag=1
nerrmpc=nerrmpc+1
endif
enddo
nmpcbad(nerrmpc)=nmpcbad(nerrmpc)+1
if( ncrcflag .eq. 1 ) then
if( nueflag .eq. 0 ) then
ngood=ngood+1
nerrdec(nerr)=nerrdec(nerr)+1
else if( nueflag .eq. 1 ) then
nue=nue+1;
endif
endif
endif
enddo
baud=12000/1920
snr2500=db+10.0*log10((baud/2500.0))
pberr=real(nberr)/(real(ntrials*N))
write(*,"(f4.1,4x,f5.1,1x,i8,1x,i8,1x,i8,8x,f5.2,8x,e10.3)") db,snr2500,ngood,nue,nbadcrc,ss,pberr
enddo
open(unit=23,file='nerrhisto.dat',status='unknown')
do i=1,174
write(23,'(i4,2x,i10,i10,f10.2)') i,nerrdec(i),nerrtot(i),real(nerrdec(i))/real(nerrtot(i)+1e-10)
enddo
close(23)
open(unit=25,file='nmpcbad.dat',status='unknown')
do i=1,87
write(25,'(i4,2x,i10)') i,nmpcbad(i)
enddo
close(25)
end program ldpcsim174
lib/fsk4hf/ldpcsim300.f90
+16 -14
View File
@@ -47,17 +47,19 @@ nerrdec=0
nmpcbad=0 ! Used to collect the number of errors in the message+crc part of the codeword
nargs=iargc()
if(nargs.ne.3) then
print*,'Usage: ldpcsim niter #trials s '
print*,'eg: ldpcsim 100 1000 0.84'
if(nargs.ne.4) then
print*,'Usage: ldpcsim niter ndeep #trials s '
print*,'eg: ldpcsim 100 4 1000 0.84'
print*,'If s is negative, then value is ignored and sigma is calculated from SNR.'
return
endif
call getarg(1,arg)
read(arg,*) max_iterations
call getarg(2,arg)
read(arg,*) ntrials
read(arg,*) ndeep
call getarg(3,arg)
read(arg,*) ntrials
call getarg(4,arg)
read(arg,*) s
fsk=.false.
@@ -147,14 +149,14 @@ do idb = 20,-16,-1
do i=1,N
if( rxdata(i)*(2*codeword(i)-1.0) .lt. 0 ) nerr=nerr+1
enddo
nerrtot(nerr)=nerrtot(nerr)+1
if(nerr.ge.1) nerrtot(nerr)=nerrtot(nerr)+1
nberr=nberr+nerr
! Correct signal normalization is important for this decoder.
! rxav=sum(rxdata)/N
! rx2av=sum(rxdata*rxdata)/N
! rxsig=sqrt(rx2av-rxav*rxav)
! rxdata=rxdata/rxsig
rxav=sum(rxdata)/N
rx2av=sum(rxdata*rxdata)/N
rxsig=sqrt(rx2av-rxav*rxav)
rxdata=rxdata/rxsig
! To match the metric to the channel, s should be set to the noise standard deviation.
! For now, set s to the value that optimizes decode probability near threshold.
! The s parameter can be tuned to trade a few tenth's dB of threshold for an order of
@@ -169,9 +171,9 @@ do idb = 20,-16,-1
apmask=0
! max_iterations is max number of belief propagation iterations
call bpdecode300(llr, apmask, max_iterations, decoded, niterations, cw)
if( niterations .lt. 0 ) then
norder=3
call osd300(llr, norder, decoded, niterations, cw)
if( (niterations .lt. 0) .and. (ndeep .ge. 0) ) then
call osd300(llr, apmask, ndeep, decoded, cw, nhardmin, dmin)
niterations=nhardmin
endif
n2err=0
do i=1,N
@@ -221,10 +223,10 @@ do idb = 20,-16,-1
nerrmpc=nerrmpc+1
endif
enddo
nmpcbad(nerrmpc)=nmpcbad(nerrmpc)+1 ! This histogram should inform our selection of CRC poly
if(nerrmpc.ge.1) nmpcbad(nerrmpc)=nmpcbad(nerrmpc)+1 ! This histogram should inform our selection of CRC poly
if( ncrcflag .eq. 1 .and. nueflag .eq. 0 ) then
ngood=ngood+1
nerrdec(nerr)=nerrdec(nerr)+1
if(nerr.ge.1) nerrdec(nerr)=nerrdec(nerr)+1
else if( ncrcflag .eq. 1 .and. nueflag .eq. 1 ) then
nue=nue+1;
endif
lib/fsk4hf/osd174.f90
-377
View File
@@ -1,377 +0,0 @@
subroutine osd174(llr,apmask,ndeep,decoded,cw,nhardmin,dmin)
!
! An ordered-statistics decoder for the (174,87) code.
!
include "ldpc_174_87_params.f90"
integer*1 apmask(N),apmaskr(N)
integer*1 gen(K,N)
integer*1 genmrb(K,N),g2(N,K)
integer*1 temp(K),m0(K),me(K),mi(K),misub(K),e2sub(N-K),e2(N-K),ui(N-K)
integer*1 r2pat(N-K)
integer indices(N),nxor(N)
integer*1 cw(N),ce(N),c0(N),hdec(N)
integer*1 decoded(K)
integer indx(N)
real llr(N),rx(N),absrx(N)
logical first,reset
data first/.true./
save first,gen
if( first ) then ! fill the generator matrix
gen=0
do i=1,M
do j=1,22
read(g(i)(j:j),"(Z1)") istr
do jj=1, 4
irow=(j-1)*4+jj
if( btest(istr,4-jj) ) gen(irow,i)=1
enddo
enddo
enddo
do irow=1,K
gen(irow,M+irow)=1
enddo
first=.false.
endif
! Re-order received vector to place systematic msg bits at the end.
rx=llr(colorder+1)
apmaskr=apmask(colorder+1)
! Hard decisions on the received word.
hdec=0
where(rx .ge. 0) hdec=1
! Use magnitude of received symbols as a measure of reliability.
absrx=abs(rx)
call indexx(absrx,N,indx)
! Re-order the columns of the generator matrix in order of decreasing reliability.
do i=1,N
genmrb(1:K,i)=gen(1:K,indx(N+1-i))
indices(i)=indx(N+1-i)
enddo
! Do gaussian elimination to create a generator matrix with the most reliable
! received bits in positions 1:K in order of decreasing reliability (more or less).
do id=1,K ! diagonal element indices
do icol=id,K+20 ! The 20 is ad hoc - beware
iflag=0
if( genmrb(id,icol) .eq. 1 ) then
iflag=1
if( icol .ne. id ) then ! reorder column
temp(1:K)=genmrb(1:K,id)
genmrb(1:K,id)=genmrb(1:K,icol)
genmrb(1:K,icol)=temp(1:K)
itmp=indices(id)
indices(id)=indices(icol)
indices(icol)=itmp
endif
do ii=1,K
if( ii .ne. id .and. genmrb(ii,id) .eq. 1 ) then
genmrb(ii,1:N)=ieor(genmrb(ii,1:N),genmrb(id,1:N))
endif
enddo
exit
endif
enddo
enddo
g2=transpose(genmrb)
! The hard decisions for the K MRB bits define the order 0 message, m0.
! Encode m0 using the modified generator matrix to find the "order 0" codeword.
! Flip various combinations of bits in m0 and re-encode to generate a list of
! codewords. Return the member of the list that has the smallest Euclidean
! distance to the received word.
hdec=hdec(indices) ! hard decisions from received symbols
m0=hdec(1:K) ! zero'th order message
absrx=absrx(indices)
rx=rx(indices)
apmaskr=apmaskr(indices)
call mrbencode(m0,c0,g2,N,K)
nxor=ieor(c0,hdec)
nhardmin=sum(nxor)
dmin=sum(nxor*absrx)
cw=c0
ntotal=0
nrejected=0
if(ndeep.eq.0) goto 998 ! norder=0
if(ndeep.gt.5) ndeep=5
if( ndeep.eq. 1) then
nord=1
npre1=0
npre2=0
nt=40
ntheta=12
elseif(ndeep.eq.2) then
nord=1
npre1=1
npre2=0
nt=40
ntheta=12
elseif(ndeep.eq.3) then
nord=1
npre1=1
npre2=1
nt=40
ntheta=12
ntau=14
elseif(ndeep.eq.4) then
nord=2
npre1=1
npre2=0
nt=40
ntheta=12
ntau=19
elseif(ndeep.eq.5) then
nord=2
npre1=1
npre2=1
nt=40
ntheta=12
ntau=19
endif
do iorder=1,nord
if( iorder.eq. 1 ) then
misub(1:K-1)=0
misub(K)=1
iflag=K
elseif( iorder.eq. 2 ) then
misub(1:K-2)=0
misub(K-1:K)=1
iflag=K-1
endif
do while(iflag .ge.0)
if(iorder.eq.nord .and. npre1.eq.0) then
iend=iflag
else
iend=1
endif
do n1=iflag,iend,-1
mi=misub
mi(n1)=1
if(any(iand(apmaskr(1:K),mi).eq.1)) cycle
ntotal=ntotal+1
me=ieor(m0,mi)
if(n1.eq.iflag) then
call mrbencode(me,ce,g2,N,K)
e2sub=ieor(ce(K+1:N),hdec(K+1:N))
e2=e2sub
nd1Kpt=sum(e2sub(1:nt))+1
d1=sum(ieor(me(1:K),hdec(1:K))*absrx(1:K))
else
e2=ieor(e2sub,g2(K+1:N,n1))
nd1Kpt=sum(e2(1:nt))+2
endif
if(nd1Kpt .le. ntheta) then
call mrbencode(me,ce,g2,N,K)
nxor=ieor(ce,hdec)
if(n1.eq.iflag) then
dd=d1+sum(e2sub*absrx(K+1:N))
else
dd=d1+ieor(ce(n1),hdec(n1))*absrx(n1)+sum(e2*absrx(K+1:N))
endif
if( dd .lt. dmin ) then
dmin=dd
cw=ce
nhardmin=sum(nxor)
nd1Kptbest=nd1Kpt
endif
else
nrejected=nrejected+1
endif
enddo
! Get the next test error pattern, iflag will go negative
! when the last pattern with weight iorder has been generated.
call nextpat(misub,k,iorder,iflag)
enddo
enddo
if(npre2.eq.1) then
reset=.true.
ntotal=0
do i1=K,1,-1
do i2=i1-1,1,-1
ntotal=ntotal+1
mi(1:ntau)=ieor(g2(K+1:K+ntau,i1),g2(K+1:K+ntau,i2))
call boxit(reset,mi(1:ntau),ntau,ntotal,i1,i2)
enddo
enddo
ncount2=0
ntotal2=0
reset=.true.
! Now run through again and do the second pre-processing rule
if(nord.eq.1) then
misub(1:K-1)=0
misub(K)=1
iflag=K
elseif(nord.eq.2) then
misub(1:K-1)=0
misub(K-1:K)=1
iflag=K-1
endif
do while(iflag .ge.0)
me=ieor(m0,misub)
call mrbencode(me,ce,g2,N,K)
e2sub=ieor(ce(K+1:N),hdec(K+1:N))
do i2=0,ntau
ntotal2=ntotal2+1
ui=0
if(i2.gt.0) ui(i2)=1
r2pat=ieor(e2sub,ui)
778 continue
call fetchit(reset,r2pat(1:ntau),ntau,in1,in2)
if(in1.gt.0.and.in2.gt.0) then
ncount2=ncount2+1
mi=misub
mi(in1)=1
mi(in2)=1
if(sum(mi).lt.nord+npre1+npre2.or.any(iand(apmaskr(1:K),mi).eq.1)) cycle
me=ieor(m0,mi)
call mrbencode(me,ce,g2,N,K)
nxor=ieor(ce,hdec)
dd=sum(nxor*absrx)
if( dd .lt. dmin ) then
dmin=dd
cw=ce
nhardmin=sum(nxor)
endif
goto 778
endif
enddo
call nextpat(misub,K,nord,iflag)
enddo
endif
998 continue
! Re-order the codeword to place message bits at the end.
cw(indices)=cw
hdec(indices)=hdec
decoded=cw(K+1:N)
cw(colorder+1)=cw ! put the codeword back into received-word order
return
end subroutine osd174
subroutine mrbencode(me,codeword,g2,N,K)
integer*1 me(K),codeword(N),g2(N,K)
! fast encoding for low-weight test patterns
codeword=0
do i=1,K
if( me(i) .eq. 1 ) then
codeword=ieor(codeword,g2(1:N,i))
endif
enddo
return
end subroutine mrbencode
subroutine nextpat(mi,k,iorder,iflag)
integer*1 mi(k),ms(k)
! generate the next test error pattern
ind=-1
do i=1,k-1
if( mi(i).eq.0 .and. mi(i+1).eq.1) ind=i
enddo
if( ind .lt. 0 ) then ! no more patterns of this order
iflag=ind
return
endif
ms=0
ms(1:ind-1)=mi(1:ind-1)
ms(ind)=1
ms(ind+1)=0
if( ind+1 .lt. k ) then
nz=iorder-sum(ms)
ms(k-nz+1:k)=1
endif
mi=ms
do i=1,k ! iflag will point to the lowest-index 1 in mi
if(mi(i).eq.1) then
iflag=i
exit
endif
enddo
return
end subroutine nextpat
subroutine boxit(reset,e2,ntau,npindex,i1,i2)
integer*1 e2(1:ntau)
integer indexes(4000,2),fp(0:525000),np(4000)
logical reset
common/boxes/indexes,fp,np
if(reset) then
patterns=-1
fp=-1
np=-1
sc=-1
indexes=-1
reset=.false.
endif
indexes(npindex,1)=i1
indexes(npindex,2)=i2
ipat=0
do i=1,ntau
if(e2(i).eq.1) then
ipat=ipat+ishft(1,ntau-i)
endif
enddo
ip=fp(ipat) ! see what's currently stored in fp(ipat)
if(ip.eq.-1) then
fp(ipat)=npindex
else
do while (np(ip).ne.-1)
ip=np(ip)
enddo
np(ip)=npindex
endif
return
end subroutine boxit
subroutine fetchit(reset,e2,ntau,i1,i2)
integer indexes(4000,2),fp(0:525000),np(4000)
integer lastpat
integer*1 e2(ntau)
logical reset
common/boxes/indexes,fp,np
save lastpat,inext
if(reset) then
lastpat=-1
reset=.false.
endif
ipat=0
do i=1,ntau
if(e2(i).eq.1) then
ipat=ipat+ishft(1,ntau-i)
endif
enddo
index=fp(ipat)
if(lastpat.ne.ipat .and. index.gt.0) then ! return first set of indices
i1=indexes(index,1)
i2=indexes(index,2)
inext=np(index)
elseif(lastpat.eq.ipat .and. inext.gt.0) then
i1=indexes(inext,1)
i2=indexes(inext,2)
inext=np(inext)
else
i1=-1
i2=-1
inext=-1
endif
lastpat=ipat
return
end subroutine fetchit
lib/fsk4hf/osd300.f90
+239 -58
View File
@@ -1,30 +1,31 @@
subroutine osd300(llr,norder,decoded,niterations,cw)
subroutine osd300(llr,apmask,ndeep,decoded,cw,nhardmin,dmin)
!
! An ordered-statistics decoder for the (300,60) code.
!
!
include "ldpc_300_60_params.f90"
integer*1 apmask(N),apmaskr(N)
integer*1 gen(K,N)
integer*1 genmrb(K,N),g2(N,K)
integer*1 temp(K),m0(K),me(K),mi(K)
integer*1 temp(K),m0(K),me(K),mi(K),misub(K),e2sub(N-K),e2(N-K),ui(N-K)
integer*1 r2pat(N-K)
integer indices(N),nxor(N)
integer*1 cw(N),ce(N),c0(N),hdec(N)
integer*1 decoded(K)
integer indx(N)
real llr(N),rx(N),absrx(N)
logical first
logical first,reset
data first/.true./
save first,gen
if( first ) then ! fill the generator matrix
gen=0
do i=1,M
do j=1,15
do j=1, 15
read(g(i)(j:j),"(Z1)") istr
do jj=1, 4
irow=(j-1)*4+jj
if( btest(istr,4-jj) ) gen(irow,i)=1
if( btest(istr,4-jj) ) gen(irow,i)=1
enddo
enddo
enddo
@@ -34,28 +35,26 @@ if( first ) then ! fill the generator matrix
first=.false.
endif
! re-order received vector to place systematic msg bits at the end
! Re-order received vector to place systematic msg bits at the end.
rx=llr(colorder+1)
apmaskr=apmask(colorder+1)
! hard decode the received word
! Hard decisions on the received word.
hdec=0
where(rx .ge. 0) hdec=1
! use magnitude of received symbols as a measure of reliability.
! Use magnitude of received symbols as a measure of reliability.
absrx=abs(rx)
call indexx(absrx,N,indx)
! re-order the columns of the generator matrix in order of decreasing reliability.
! Re-order the columns of the generator matrix in order of decreasing reliability.
do i=1,N
genmrb(1:K,i)=gen(1:K,indx(N+1-i))
indices(i)=indx(N+1-i)
enddo
! do gaussian elimination to create a generator matrix with the most reliable
! Do gaussian elimination to create a generator matrix with the most reliable
! received bits in positions 1:K in order of decreasing reliability (more or less).
! reliability will not be strictly decreasing because column re-ordering is needed
! to put the generator matrix in systematic form. the "indices" array tracks
! column permutations caused by reliability sorting and gaussian elimination.
do id=1,K ! diagonal element indices
do icol=id,K+20 ! The 20 is ad hoc - beware
iflag=0
@@ -71,7 +70,7 @@ do id=1,K ! diagonal element indices
endif
do ii=1,K
if( ii .ne. id .and. genmrb(ii,id) .eq. 1 ) then
genmrb(ii,1:N)=mod(genmrb(ii,1:N)+genmrb(id,1:N),2)
genmrb(ii,1:N)=ieor(genmrb(ii,1:N),genmrb(id,1:N))
endif
enddo
exit
@@ -84,66 +83,168 @@ g2=transpose(genmrb)
! The hard decisions for the K MRB bits define the order 0 message, m0.
! Encode m0 using the modified generator matrix to find the "order 0" codeword.
! Flip various combinations of bits in m0 and re-encode to generate a list of
! codewords. Test all such codewords against the received word to decide which
! codeword is most likely to be correct.
! codewords. Return the member of the list that has the smallest Euclidean
! distance to the received word.
hdec=hdec(indices) ! hard decisions from received symbols
m0=hdec(1:K) ! zero'th order message
absrx=absrx(indices)
rx=rx(indices)
apmaskr=apmaskr(indices)
s1=sum(absrx(1:K))
s2=sum(absrx(K+1:N))
xlam=5.0
rho=s1/(s1+xlam*s2)
call mrbencode(m0,c0,g2,N,K)
nxor=ieor(c0,hdec)
nhardmin=sum(nxor)
dmin=sum(nxor*absrx)
thresh=rho*dmin
cw=c0
nt=0
ntotal=0
nrejected=0
do iorder=1,norder
mi(1:K-iorder)=0
mi(K-iorder+1:K)=1
iflag=0
do while(iflag .ge. 0 )
dpat=sum(mi*absrx(1:K))
nt=nt+1
if( dpat .lt. thresh ) then ! reject unlikely error patterns
me=ieor(m0,mi)
call mrbencode(me,ce,g2,N,K)
nxor=ieor(ce,hdec)
dd=sum(nxor*absrx)
if( dd .lt. dmin ) then
dmin=dd
cw=ce
nhardmin=sum(nxor)
thresh=rho*dmin
if(ndeep.eq.0) goto 998 ! norder=0
if(ndeep.gt.5) ndeep=5
if( ndeep.eq. 1) then
nord=1
npre1=0
npre2=0
nt=120
ntheta=12
elseif(ndeep.eq.2) then
nord=1
npre1=1
npre2=0
nt=120
ntheta=12
elseif(ndeep.eq.3) then
nord=1
npre1=1
npre2=1
nt=120
ntheta=12
ntau=15
elseif(ndeep.eq.4) then
nord=2
npre1=1
npre2=0
nt=120
ntheta=12
ntau=15
elseif(ndeep.eq.5) then
nord=4
npre1=1
npre2=1
nt=120
ntheta=20
ntau=15
endif
do iorder=1,nord
misub(1:K-iorder)=0
misub(K-iorder+1:K)=1
iflag=K-iorder+1
do while(iflag .ge.0)
if(iorder.eq.nord .and. npre1.eq.0) then
iend=iflag
else
iend=1
endif
else
nrejected=nrejected+1
endif
! get the next test error pattern, iflag will go negative
! when the last pattern with weight iorder has been generated
call nextpat(mi,k,iorder,iflag)
enddo
do n1=iflag,iend,-1
mi=misub
mi(n1)=1
if(any(iand(apmaskr(1:K),mi).eq.1)) cycle
ntotal=ntotal+1
me=ieor(m0,mi)
if(n1.eq.iflag) then
call mrbencode(me,ce,g2,N,K)
e2sub=ieor(ce(K+1:N),hdec(K+1:N))
e2=e2sub
nd1Kpt=sum(e2sub(1:nt))+1
d1=sum(ieor(me(1:K),hdec(1:K))*absrx(1:K))
else
e2=ieor(e2sub,g2(K+1:N,n1))
nd1Kpt=sum(e2(1:nt))+2
endif
if(nd1Kpt .le. ntheta) then
call mrbencode(me,ce,g2,N,K)
nxor=ieor(ce,hdec)
if(n1.eq.iflag) then
dd=d1+sum(e2sub*absrx(K+1:N))
else
dd=d1+ieor(ce(n1),hdec(n1))*absrx(n1)+sum(e2*absrx(K+1:N))
endif
if( dd .lt. dmin ) then
dmin=dd
cw=ce
nhardmin=sum(nxor)
nd1Kptbest=nd1Kpt
endif
else
nrejected=nrejected+1
endif
enddo
! Get the next test error pattern, iflag will go negative
! when the last pattern with weight iorder has been generated.
call nextpat(misub,k,iorder,iflag)
enddo
enddo
!write(*,*) 'nhardmin ',nhardmin
!write(*,*) 'total patterns ',nt,' number rejected ',nrejected
if(npre2.eq.1) then
reset=.true.
ntotal=0
do i1=K,1,-1
do i2=i1-1,1,-1
ntotal=ntotal+1
mi(1:ntau)=ieor(g2(K+1:K+ntau,i1),g2(K+1:K+ntau,i2))
call boxit(reset,mi(1:ntau),ntau,ntotal,i1,i2)
enddo
enddo
! re-order the codeword to place message bits at the end
ncount2=0
ntotal2=0
reset=.true.
! Now run through again and do the second pre-processing rule
misub(1:K-nord)=0
misub(K-nord+1:K)=1
iflag=K-nord+1
do while(iflag .ge.0)
me=ieor(m0,misub)
call mrbencode(me,ce,g2,N,K)
e2sub=ieor(ce(K+1:N),hdec(K+1:N))
do i2=0,ntau
ntotal2=ntotal2+1
ui=0
if(i2.gt.0) ui(i2)=1
r2pat=ieor(e2sub,ui)
778 continue
call fetchit(reset,r2pat(1:ntau),ntau,in1,in2)
if(in1.gt.0.and.in2.gt.0) then
ncount2=ncount2+1
mi=misub
mi(in1)=1
mi(in2)=1
if(sum(mi).lt.nord+npre1+npre2.or.any(iand(apmaskr(1:K),mi).eq.1)) cycle
me=ieor(m0,mi)
call mrbencode(me,ce,g2,N,K)
nxor=ieor(ce,hdec)
dd=sum(nxor*absrx)
if( dd .lt. dmin ) then
dmin=dd
cw=ce
nhardmin=sum(nxor)
endif
goto 778
endif
enddo
call nextpat(misub,K,nord,iflag)
enddo
endif
998 continue
! Re-order the codeword to place message bits at the end.
cw(indices)=cw
hdec(indices)=hdec
decoded=cw(M+1:N)
nerr=0
do i=1,N
if( hdec(i) .ne. cw(i) ) nerr=nerr+1
enddo
niterations=nerr
decoded=cw(M+1:N)
cw(colorder+1)=cw ! put the codeword back into received-word order
return
end subroutine osd300
@@ -179,6 +280,86 @@ subroutine nextpat(mi,k,iorder,iflag)
ms(k-nz+1:k)=1
endif
mi=ms
iflag=ind
do i=1,k ! iflag will point to the lowest-index 1 in mi
if(mi(i).eq.1) then
iflag=i
exit
endif
enddo
return
end subroutine nextpat
subroutine boxit(reset,e2,ntau,npindex,i1,i2)
integer*1 e2(1:ntau)
integer indexes(4000,2),fp(0:525000),np(4000)
logical reset
common/boxes/indexes,fp,np
if(reset) then
patterns=-1
fp=-1
np=-1
sc=-1
indexes=-1
reset=.false.
endif
indexes(npindex,1)=i1
indexes(npindex,2)=i2
ipat=0
do i=1,ntau
if(e2(i).eq.1) then
ipat=ipat+ishft(1,ntau-i)
endif
enddo
ip=fp(ipat) ! see what's currently stored in fp(ipat)
if(ip.eq.-1) then
fp(ipat)=npindex
else
do while (np(ip).ne.-1)
ip=np(ip)
enddo
np(ip)=npindex
endif
return
end subroutine boxit
subroutine fetchit(reset,e2,ntau,i1,i2)
integer indexes(4000,2),fp(0:525000),np(4000)
integer lastpat
integer*1 e2(ntau)
logical reset
common/boxes/indexes,fp,np
save lastpat,inext
if(reset) then
lastpat=-1
reset=.false.
endif
ipat=0
do i=1,ntau
if(e2(i).eq.1) then
ipat=ipat+ishft(1,ntau-i)
endif
enddo
index=fp(ipat)
if(lastpat.ne.ipat .and. index.gt.0) then ! return first set of indices
i1=indexes(index,1)
i2=indexes(index,2)
inext=np(index)
elseif(lastpat.eq.ipat .and. inext.gt.0) then
i1=indexes(inext,1)
i2=indexes(inext,2)
inext=np(inext)
else
i1=-1
i2=-1
inext=-1
endif
lastpat=ipat
return
end subroutine fetchit
lib/fsk4hf/spb.m
+89
View File
@@ -0,0 +1,89 @@
clear all;
global N
global R
global A
#-------------------------------------------------------------------------------
function retval = f1(theta)
global N;
global R;
retval=0.0;
gterm = gammaln(N/2) - gammaln((N+1)/2) - log(2*sqrt(pi));
rhs = -N*R*log(2);
lhs=gterm + (N-1)*log(sin(theta)) + log(1-(tan(theta).^2)/N) - log(cos(theta));
retval = rhs-real(lhs);
endfunction
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
function retval = d(N,i,x)
t1=(x.^2)/2;
t2=gammaln(N/2);
t3=-gammaln(i/2+1);
t4=-gammaln(N-i);
t5=(N-1-i)*log(sqrt(2)*x);
t6=-log(2)/2;
t7arg=1+(-1)^i * gammainc((x.^2)/2.0,(i+1)/2);
t7=log(t7arg);
retval=t1+t2+t3+t4+t5+t6+t7;
endfunction
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
function retval = maxstar(x1,x2)
retval = max(x1,x2)+log(1+exp(-abs(x1-x2)));
endfunction
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
function retval = spb_integrand(x)
global N;
global A;
t1=log(N-1);
t2=-N*(A^2)/2;
t3=-0.5*log(2*pi);
t4=(N-2)*log(sin(x));
arg=sqrt(N)*A*cos(x);
t5=maxstar(d(N,0,arg),d(N,1,arg));
for i=2:N-1
t5=maxstar(t5,d(N,i,arg));
endfor
retval=exp(t1+t2+t3+t4+t5);
endfunction
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
function retval = qfunc(x)
retval = 0.5 * erfc(x/sqrt(2));
endfunction
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# Calculate sphere packing lower bound on the probability of word error
# given block length (N), code rate (R), and Eb/No.
#
# Ref:
# "Log-Domain Calculation of the 1959 Sphere-Packing Bound with Application to
# M-ary PSK Block Coded Modulation," Igal Sason and Gil Weichman,
# doi: 10.1109/EEEI.2006.321097
#-------------------------------------------------------------------------------
N=174
K=75
R=K/N
delta=0.01;
[ths,fval,info,output]=fzero(@f1,[delta,pi/2-delta], optimset ("jacobian", "off"));
for ebnodb=-6:0.5:4
ebno=10^(ebnodb/10.0);
esno=ebno*R;
A=sqrt(2*esno);
term1=quadcc(@spb_integrand,ths,pi/2);
term2=qfunc(sqrt(N)*A);
pe=term1+term2;
ps=1-pe;
printf("%f %f\n",ebnodb,ps);
endfor
lib/fsk4hf/subtractft8.f90
-61
View File
@@ -1,61 +0,0 @@
subroutine subtractft8(dd,itone,f0,dt)
! Subtract an ft8 signal
!
! Measured signal : dd(t) = a(t)cos(2*pi*f0*t+theta(t))
! Reference signal : cref(t) = exp( j*(2*pi*f0*t+phi(t)) )
! Complex amp : cfilt(t) = LPF[ dd(t)*CONJG(cref(t)) ]
! Subtract : dd(t) = dd(t) - 2*REAL{cref*cfilt}
use timer_module, only: timer
parameter (NMAX=15*12000,NFRAME=1920*79)
parameter (NFFT=NMAX,NFILT=1400)
real*4 dd(NMAX), window(-NFILT/2:NFILT/2)
complex cref,camp,cfilt,cw
integer itone(79)
logical first
data first/.true./
common/heap8/cref(NFRAME),camp(NMAX),cfilt(NMAX),cw(NMAX)
save first
nstart=dt*12000+1
call genft8refsig(itone,cref,f0)
camp=0.
do i=1,nframe
id=nstart-1+i
if(id.ge.1.and.id.le.NMAX) camp(i)=dd(id)*conjg(cref(i))
enddo
if(first) then
! Create and normalize the filter
pi=4.0*atan(1.0)
fac=1.0/float(nfft)
sum=0.0
do j=-NFILT/2,NFILT/2
window(j)=cos(pi*j/NFILT)**2
sum=sum+window(j)
enddo
cw=0.
cw(1:NFILT+1)=window/sum
cw=cshift(cw,NFILT/2+1)
call four2a(cw,nfft,1,-1,1)
cw=cw*fac
first=.false.
endif
cfilt=0.0
cfilt(1:nframe)=camp(1:nframe)
call four2a(cfilt,nfft,1,-1,1)
cfilt(1:nfft)=cfilt(1:nfft)*cw(1:nfft)
call four2a(cfilt,nfft,1,1,1)
! Subtract the reconstructed signal
do i=1,nframe
j=nstart+i-1
if(j.ge.1 .and. j.le.NMAX) dd(j)=dd(j)-2*REAL(cfilt(i)*cref(i))
enddo
return
end subroutine subtractft8
lib/fsk4hf/sync8.f90
-151
View File
@@ -1,151 +0,0 @@
subroutine sync8(dd,nfa,nfb,syncmin,nfqso,s,candidate,ncand,sbase)
include 'ft8_params.f90'
! Search over +/- 2.5s relative to 0.5s TX start time.
parameter (JZ=62)
complex cx(0:NH1)
real s(NH1,NHSYM)
real savg(NH1)
real sbase(NH1)
real x(NFFT1)
real sync2d(NH1,-JZ:JZ)
real red(NH1)
real candidate0(3,200)
real candidate(3,200)
real dd(NMAX)
integer jpeak(NH1)
integer indx(NH1)
integer ii(1)
integer icos7(0:6)
data icos7/2,5,6,0,4,1,3/ !Costas 7x7 tone pattern
equivalence (x,cx)
! Compute symbol spectra, stepping by NSTEP steps.
savg=0.
tstep=NSTEP/12000.0
df=12000.0/NFFT1 !3.125 Hz
fac=1.0/300.0
do j=1,NHSYM
ia=(j-1)*NSTEP + 1
ib=ia+NSPS-1
x(1:NSPS)=fac*dd(ia:ib)
x(NSPS+1:)=0.
call four2a(x,NFFT1,1,-1,0) !r2c FFT
do i=1,NH1
s(i,j)=real(cx(i))**2 + aimag(cx(i))**2
enddo
savg=savg + s(1:NH1,j) !Average spectrum
enddo
call baseline(savg,nfa,nfb,sbase)
! savg=savg/NHSYM
! do i=1,NH1
! write(51,3051) i*df,savg(i),db(savg(i))
!3051 format(f10.3,e12.3,f12.3)
! enddo
ia=max(1,nint(nfa/df))
ib=nint(nfb/df)
nssy=NSPS/NSTEP ! # steps per symbol
nfos=NFFT1/NSPS ! # frequency bin oversampling factor
jstrt=0.5/tstep
do i=ia,ib
do j=-JZ,+JZ
ta=0.
tb=0.
tc=0.
t0a=0.
t0b=0.
t0c=0.
do n=0,6
k=j+jstrt+nssy*n
if(k.ge.1.and.k.le.NHSYM) then
ta=ta + s(i+nfos*icos7(n),k)
t0a=t0a + sum(s(i:i+nfos*6:nfos,k))
endif
tb=tb + s(i+nfos*icos7(n),k+nssy*36)
t0b=t0b + sum(s(i:i+nfos*6:nfos,k+nssy*36))
if(k+nssy*72.le.NHSYM) then
tc=tc + s(i+nfos*icos7(n),k+nssy*72)
t0c=t0c + sum(s(i:i+nfos*6:nfos,k+nssy*72))
endif
enddo
t=ta+tb+tc
t0=t0a+t0b+t0c
t0=(t0-t)/6.0
sync_abc=t/t0
t=tb+tc
t0=t0b+t0c
t0=(t0-t)/6.0
sync_bc=t/t0
sync2d(i,j)=max(sync_abc,sync_bc)
enddo
enddo
red=0.
do i=ia,ib
ii=maxloc(sync2d(i,-JZ:JZ)) - 1 - JZ
j0=ii(1)
jpeak(i)=j0
red(i)=sync2d(i,j0)
! write(52,3052) i*df,red(i),db(red(i))
!3052 format(3f12.3)
enddo
iz=ib-ia+1
call indexx(red(ia:ib),iz,indx)
ibase=indx(nint(0.40*iz)) - 1 + ia
base=red(ibase)
red=red/base
candidate0=0.
k=0
do i=1,200
n=ia + indx(iz+1-i) - 1
if(red(n).lt.syncmin) exit
if(k.lt.200) k=k+1
candidate0(1,k)=n*df
candidate0(2,k)=(jpeak(n)-1)*tstep
candidate0(3,k)=red(n)
enddo
ncand=k
! Put nfqso at top of list, and save only the best of near-dupe freqs.
do i=1,ncand
if(abs(candidate0(1,i)-nfqso).lt.10.0) candidate0(1,i)=-candidate0(1,i)
if(i.ge.2) then
do j=1,i-1
fdiff=abs(candidate0(1,i))-abs(candidate0(1,j))
if(abs(fdiff).lt.4.0) then
if(candidate0(3,i).ge.candidate0(3,j)) candidate0(3,j)=0.
if(candidate0(3,i).lt.candidate0(3,j)) candidate0(3,i)=0.
endif
enddo
! write(*,3001) i,candidate0(1,i-1),candidate0(1,i),candidate0(3,i-1), &
! candidate0(3,i)
!3001 format(i2,4f8.1)
endif
enddo
fac=20.0/maxval(s)
s=fac*s
! Sort by sync
! call indexx(candidate0(3,1:ncand),ncand,indx)
! Sort by frequency
call indexx(candidate0(1,1:ncand),ncand,indx)
k=1
! do i=ncand,1,-1
do i=1,ncand
j=indx(i)
! if( candidate0(3,j) .ge. syncmin .and. candidate0(2,j).ge.-1.5 ) then
if( candidate0(3,j) .ge. syncmin ) then
candidate(1,k)=abs(candidate0(1,j))
candidate(2,k)=candidate0(2,j)
candidate(3,k)=candidate0(3,j)
k=k+1
endif
enddo
ncand=k-1
return
end subroutine sync8
lib/fsk4hf/sync8d.f90
-54
View File
@@ -1,54 +0,0 @@
subroutine sync8d(cd0,i0,ctwk,itwk,sync)
! Compute sync power for a complex, downsampled FT8 signal.
parameter(NP2=2812,NDOWN=60)
complex cd0(3125)
complex csync(0:6,32)
complex csync2(32)
complex ctwk(32)
complex z1,z2,z3
logical first
integer icos7(0:6)
data icos7/2,5,6,0,4,1,3/
data first/.true./
save first,twopi,fs2,dt2,taus,baud,csync
p(z1)=real(z1)**2 + aimag(z1)**2 !Statement function for power
! Set some constants and compute the csync array.
if( first ) then
twopi=8.0*atan(1.0)
fs2=12000.0/NDOWN !Sample rate after downsampling
dt2=1/fs2 !Corresponding sample interval
taus=32*dt2 !Symbol duration
baud=1.0/taus !Keying rate
do i=0,6
phi=0.0
dphi=twopi*icos7(i)*baud*dt2
do j=1,32
csync(i,j)=cmplx(cos(phi),sin(phi)) !Waveform for 7x7 Costas array
phi=mod(phi+dphi,twopi)
enddo
enddo
first=.false.
endif
sync=0
do i=0,6 !Sum over 7 Costas frequencies and
i1=i0+i*32 !three Costas arrays
i2=i1+36*32
i3=i1+72*32
csync2=csync(i,1:32)
if(itwk.eq.1) csync2=ctwk*csync2 !Tweak the frequency
z1=0.
z2=0.
z3=0.
if(i1.ge.1 .and. i1+31.le.NP2) z1=sum(cd0(i1:i1+31)*conjg(csync2))
if(i2.ge.1 .and. i2+31.le.NP2) z2=sum(cd0(i2:i2+31)*conjg(csync2))
if(i3.ge.1 .and. i3+31.le.NP2) z3=sum(cd0(i3:i3+31)*conjg(csync2))
sync = sync + p(z1) + p(z2) + p(z3)
enddo
return
end subroutine sync8d
lib/fsk4hf/twkfreq1.f90
-26
View File
@@ -1,26 +0,0 @@
subroutine twkfreq1(ca,npts,fsample,a,cb)
complex ca(npts)
complex cb(npts)
complex w,wstep
real a(5)
data twopi/6.283185307/
! Mix the complex signal
w=1.0
wstep=1.0
x0=0.5*(npts+1)
s=2.0/npts
do i=1,npts
x=s*(i-x0)
p2=1.5*x*x - 0.5
p3=2.5*(x**3) - 1.5*x
p4=4.375*(x**4) - 3.75*(x**2) + 0.375
dphi=(a(1) + x*a(2) + p2*a(3) + p3*a(4) + p4*a(5)) * (twopi/fsample)
wstep=cmplx(cos(dphi),sin(dphi))
w=w*wstep
cb(i)=w*ca(i)
enddo
return
end subroutine twkfreq1
lib/fsk4hf/watterson.f90
-62
View File
@@ -1,62 +0,0 @@
subroutine watterson(c,npts,fs,delay,fspread)
complex c(0:npts-1)
complex c2(0:npts-1)
complex cs1(0:npts-1)
complex cs2(0:npts-1)
nonzero=0
df=fs/npts
if(fspread.gt.0.0) then
do i=0,npts-1
xx=gran()
yy=gran()
cs1(i)=0.707*cmplx(xx,yy)
xx=gran()
yy=gran()
cs2(i)=0.707*cmplx(xx,yy)
enddo
call four2a(cs1,npts,1,-1,1) !To freq domain
call four2a(cs2,npts,1,-1,1)
do i=0,npts-1
f=i*df
if(i.gt.npts/2) f=(i-npts)*df
x=(f/(0.5*fspread))**2
a=0.
if(x.le.50.0) then
a=exp(-x)
endif
cs1(i)=a*cs1(i)
cs2(i)=a*cs2(i)
if(abs(f).lt.10.0) then
p1=real(cs1(i))**2 + aimag(cs1(i))**2
p2=real(cs2(i))**2 + aimag(cs2(i))**2
if(p1.gt.0.0) nonzero=nonzero+1
! write(62,3101) f,p1,p2,db(p1+1.e-12)-60,db(p2+1.e-12)-60
!3101 format(f10.3,2f12.3,2f10.3)
endif
enddo
call four2a(cs1,npts,1,1,1) !Back to time domain
call four2a(cs2,npts,1,1,1)
cs1(0:npts-1)=cs1(0:npts-1)/npts
cs2(0:npts-1)=cs2(0:npts-1)/npts
endif
nshift=nint(0.001*delay*fs)
c2(0:npts-1)=cshift(c(0:npts-1),nshift)
sq=0.
do i=0,npts-1
if(nonzero.gt.1) then
c(i)=0.5*(cs1(i)*c(i) + cs2(i)*c2(i))
else
c(i)=0.5*(c(i) + c2(i))
endif
sq=sq + real(c(i))**2 + aimag(c(i))**2
! write(61,3001) i/12000.0,c(i)
!3001 format(3f12.6)
enddo
rms=sqrt(sq/npts)
c=c/rms
return
end subroutine watterson
lib/fsk4hf/wspr5d.f90
+6 -3
View File
@@ -63,7 +63,8 @@ program wspr5d
open(13,file=trim(data_dir)//'/ALL_WSPR.TXT',status='unknown', &
position='append')
maxn=8 !Default value
! maxn=8 !Default value
maxn=20
twopi=8.0*atan(1.0)
fs=NSPS*12000.0/NSPS0 !Sample rate
dt=1.0/fs !Sample interval (s)
@@ -203,12 +204,14 @@ jpk=fs*xdt
call wqdecode(idat,message,itype)
nsnr=nint(xsnr)
! freq=fMHz + 1.d-6*(fc1+fc2)
freq=fMHz + 1.d-6*(fc1+fpks(itry))
! freq=fMHz + 1.d-6*(fc1+fpks(itry))
freq=fc1+fpks(itry)
nfdot=0
write(13,1110) datetime,0,nsnr,xdt,freq,message,nfdot
1110 format(a11,2i4,f6.2,f12.7,2x,a22,i3)
write(*,1112) datetime(8:11),nsnr,xdt,freq,nfdot,message,itry
1112 format(a4,i4,f5.1,f11.6,i3,2x,a22,i4)
!1112 format(a4,i4,f5.1,f11.6,i3,2x,a22,i4)
1112 format(a4,i4,f8.3,f8.3,i3,2x,a22,i4)
endif
enddo ! ifile loop
write(*,1120)
lib/fsk4hf/wspr5d_exp.f90
+4 -2
View File
@@ -208,12 +208,14 @@ program wspr5d
idat(7)=ishft(idat(7),6)
call wqdecode(idat,message,itype)
nsnr=nint(xsnr)
freq=fMHz + 1.d-6*(fc1+fc2)
! freq=fMHz + 1.d-6*(fc1+fc2)
freq=fc1+fc2
nfdot=0
write(13,1210) datetime,0,nsnr,xdt,freq,message,nfdot
1210 format(a11,2i4,f6.2,f12.7,2x,a22,i3)
write(*,1212) datetime(8:11),nsnr,xdt,freq,nfdot,message,'*',idf,nseq,is,niterations
1212 format(a4,i4,f5.1,f11.6,i3,2x,a22,a1,i3,i3,i3,i4)
!1212 format(a4,i4,f5.1,f11.6,i3,2x,a22,a1,i3,i3,i3,i4)
1212 format(a4,i4,f8.3,f8.3,i3,2x,a22,a1,i3,i3,i3,i4)
goto 888
endif
enddo !iseq
lib/fsk4hf/wspr_fsk8d.f90
+3 -3
View File
@@ -147,8 +147,8 @@ program wspr_fsk8d
do j=1,ND
k=j+7
ps=s(0:7,k)
! ps=sqrt(ps) !### ??? ###
ps=log(ps)
ps=sqrt(ps) !### ??? ###
! ps=log(ps)
r1=max(ps(1),ps(3),ps(5),ps(7))-max(ps(0),ps(2),ps(4),ps(6))
r2=max(ps(2),ps(3),ps(6),ps(7))-max(ps(0),ps(1),ps(4),ps(5))
r4=max(ps(4),ps(5),ps(6),ps(7))-max(ps(0),ps(1),ps(2),ps(3))
@@ -161,7 +161,7 @@ program wspr_fsk8d
rx2av=sum(rxdata*rxdata)/ND
rxsig=sqrt(rx2av-rxav*rxav)
rxdata=rxdata/rxsig
s0=0.84
s0=1.1
llr=2.0*rxdata/(s0*s0)
apmask=0
max_iterations=40
lib/fsk4hf/wspr_params.f90
+23
View File
@@ -0,0 +1,23 @@
parameter (NN=162)
parameter (NSPS0=8192) !Samples per symbol at 12000 S/s
parameter (NDOWN=32)
parameter (NSPS=NSPS0/NDOWN)
parameter (NZ=NSPS*NN) !Samples in waveform at 12000 S/s
parameter (NZ0=NSPS0*NN) !Samples in waveform at 375 S/s
parameter (NMAX=120*12000) !Samples in waveform at 375 S/s
! Define the sync vector:
integer*1 sync(162)
data sync/ &
1,1,0,0,0,0,0,0,1,0,0,0,1,1,1,0,0,0,1,0, &
0,1,0,1,1,1,1,0,0,0,0,0,0,0,1,0,0,1,0,1, &
0,0,0,0,0,0,1,0,1,1,0,0,1,1,0,1,0,0,0,1, &
1,0,1,0,0,0,0,1,1,0,1,0,1,0,1,0,1,0,0,1, &
0,0,1,0,1,1,0,0,0,1,1,0,1,0,1,0,0,0,1,0, &
0,0,0,0,1,0,0,1,0,0,1,1,1,0,1,1,0,0,1,1, &
0,1,0,0,0,1,1,1,0,0,0,0,0,1,0,1,0,0,1,1, &
0,0,0,0,0,0,0,1,1,0,1,0,1,1,0,0,0,1,1,0, &
0,0/
lib/fsk4hf/wspr_wav.f90
+49
View File
@@ -0,0 +1,49 @@
subroutine wspr_wav(baud,xdt,h,f0,itone,snrdb,iwave)
! Generate iwave() from itone().
include 'wspr_params.f90'
integer itone(NN)
integer*2 iwave(NMAX)
real*8 twopi,dt,dphi0,dphi1,dphi,phi
real dat(NMAX)
twopi=8.d0*atan(1.d0)
dt=1.d0/12000.d0
baud=375.0/256.0
dat=0.
if(snrdb.lt.90) then
do i=1,NMAX
dat(i)=gran() !Generate gaussian noise
enddo
bandwidth_ratio=2500.0/6000.0
sig=sqrt(2*bandwidth_ratio)*10.0**(0.05*snrdb)
else
sig=1.0
endif
phi=0.d0
k=nint(xdt/dt)
do j=1,NN
dphi=twopi*(f0+h*(itone(j)-1.5)*baud)*dt
do i=1,NSPS0
k=k+1
phi=mod(phi+dphi,twopi)
if(k.gt.0 .and. k.le.NMAX) dat(k)=dat(k) + sig*sin(phi)
enddo
enddo
rms=100.0
if(snrdb.lt.90.0) then
dat=rms*dat;
if(maxval(abs(dat)).gt.32767.0) print*,"Warning - data will be clipped."
else
datpk=maxval(abs(dat))
fac=32767.9/datpk
dat=fac*dat
endif
iwave=nint(dat)
return
end subroutine wspr_wav
lib/fsk4hf/wsprlfsim.f90
+3 -2
View File
@@ -124,7 +124,7 @@ program wsprlfsim
call system_clock(count2,clkfreq)
do iter=1,iters !Loop over requested iterations
c=c0
write(*,*) 'iter ',iter
call system_clock(count0,clkfreq)
if(delay.ne.0.0 .or. fspread.ne.0.0) then
call watterson(c,NZ,fs,delay,fspread)
@@ -151,9 +151,10 @@ program wsprlfsim
call system_clock(count1,clkfreq)
t(5)=t(5)+float(count1-count0)/float(clkfreq)
nc(5)=nc(5)+1
write(*,*) 'fc1 ',fc1
call system_clock(count0,clkfreq)
call getfc2w(c,csync,fs,fc1,fc2,fc3) !Refined freq
write(*,*) 'fc1,fc2,fc3 ',fc1,fc2,fc3
call system_clock(count1,clkfreq)
t(6)=t(6)+float(count1-count0)/float(clkfreq)
nc(6)=nc(6)+1
lib/fsk4hf/wsprsimf.f90
+113
View File
@@ -0,0 +1,113 @@
!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
!-------------------------------------------------------------------------------
program wsprsim
use wavhdr
include 'wspr_params.f90'
type(hdr) hwav
character arg*12,fname14*14,fname15*15
character*22 msg,msgsent
complex c0(0:NMAX/NDOWN-1)
complex c(0:NMAX/NDOWN-1)
integer itone(NN)
integer*2 iwave(NMAX)
real*8 fMHz
! Get command-line argument(s)
nargs=iargc()
if(nargs.ne.8) then
print*,'Usage: wsprsim "message" f0 DT fsp del nwav nfiles snr'
print*,'Example: wsprsim "K1ABC FN42 30" 50 0.0 0.1 1.0 1 10 -33'
go to 999
endif
call getarg(1,msg) !Message to be transmitted
call getarg(2,arg)
read(arg,*) f0 !Freq relative to WSPR-band center (Hz)
call getarg(3,arg)
read(arg,*) xdt !Time offset from nominal (s)
call getarg(4,arg)
read(arg,*) fspread !Watterson frequency spread (Hz)
call getarg(5,arg)
read(arg,*) delay !Watterson delay (ms)
call getarg(6,arg)
read(arg,*) nwav !1 for *.wav file, 0 for *.c2 file
call getarg(7,arg)
read(arg,*) nfiles !Number of files
call getarg(8,arg)
read(arg,*) snrdb !SNR_2500
twopi=8.0*atan(1.0)
fs=12000.0/NDOWN
dt=1.0/fs
tt=NSPS*dt
baud=12000.0/8192.0
txt=NZ*dt !Transmission length (s)
bandwidth_ratio=2500.0/(fs/2.0)
sig=sqrt(bandwidth_ratio) * 10.0**(0.05*snrdb)
if(snrdb.gt.90.0) sig=1.0
txt=NN*NSPS0/12000.0
call genwspr(msg,msgsent,itone) !Encode the message, get itone
write(*,1000) f0,xdt,txt,snrdb,fspread,delay,nfiles,msgsent
1000 format('f0:',f9.3,' DT:',f6.2,' txt:',f6.1,' SNR:',f6.1, &
' fspread:',f6.1,' delay:',f6.1,' nfiles:',i3,2x,a22)
! write(*,*) "Channel symbols: "
! write(*,'(162i2)') itone
h=1.0
phi=0.0
c0=0.
k=-1 + nint(xdt/dt)
do j=1,NN
dphi=-twopi*(f0+h*(itone(j)-1.5)*baud)*dt
do i=1,NSPS
k=k+1
phi=mod(phi+dphi,twopi)
if(k.ge.0 .and. k.lt.NMAX/NDOWN) c0(k)=cmplx(cos(phi),sin(phi))
enddo
enddo
call sgran()
do ifile=1,nfiles
c=c0
if(nwav.eq.0) then
if( fspread .ne. 0.0 .or. delay .ne. 0.0 ) then
call watterson(c,NMAX/NDOWN,fs,delay,fspread)
endif
c=c*sig
if(snrdb.lt.90) then
do i=0,NMAX/NDOWN-1 !Add gaussian noise at specified SNR
xnoise=gran()
ynoise=gran()
c(i)=c(i) + cmplx(xnoise,ynoise)
enddo
endif
write(fname14,1100) ifile
1100 format('000000_',i4.4,'.c2')
open(10,file=fname14,status='unknown',access='stream')
fMHz=10.1387d0
nmin=2
write(10) fname14,nmin,fMHz,c !Save to *.c2 file
close(10)
write(*,1108) ifile,xdt,f0,snrdb,fname14
1108 format(i4,f7.2,f8.2,f7.1,2x,a14)
else
freq=1500.0+f0
call wspr_wav(baud,xdt,h,freq,itone,snrdb,iwave)
hwav=default_header(12000,NMAX)
write(fname15,1102) ifile
1102 format('000000_',i4.4,'.wav')
open(10,file=fname15,status='unknown',access='stream')
write(10) hwav,iwave !Save to *.wav file
close(10)
write(*,1110) ifile,xdt,f0,snrdb,fname15
1110 format(i4,f7.2,f8.2,f7.1,2x,a15)
endif
enddo
999 end program wsprsim