mirror of https://github.com/saitohirga/WSJT-X.git
242 lines
7.6 KiB
Fortran
242 lines
7.6 KiB
Fortran
program dbpsksim
|
|
|
|
parameter (ND=121) !Data symbols: LDPC (120,60), r=1/2
|
|
parameter (NN=ND) !Total symbols (121)
|
|
parameter (NSPS=28800) !Samples per symbol at 12000 sps
|
|
parameter (NZ=NSPS*NN) !Samples in waveform (3484800)
|
|
parameter (NFFT1=65536,NH1=NFFT1/2)
|
|
parameter (NFFT2=128,NH2=NFFT2/2)
|
|
|
|
character*8 arg
|
|
complex c(0:NZ-1) !Complex waveform
|
|
complex c2(0:NFFT1-1) !Short spectra
|
|
complex cr(0:NZ-1)
|
|
complex ct(0:NZ-1)
|
|
complex cz(0:NFFT2-1)
|
|
complex z0,z,zp
|
|
real s(-NH1+1:NH1)
|
|
real s2(-NH2+1:NH2)
|
|
real xnoise(0:NZ-1) !Generated random noise
|
|
real ynoise(0:NZ-1) !Generated random noise
|
|
real rxdata(120),llr(120)
|
|
integer id(NN) !Encoded NRZ data (values +/-1)
|
|
integer id1(NN) !Recovered data (1st pass)
|
|
integer id2(NN) !Recovered data (2nd pass)
|
|
! integer icw(NN)
|
|
integer*1 msgbits(60),decoded(60),codeword(120),apmask(120),cw(120)
|
|
data msgbits/0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,&
|
|
0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,1,0,0,1,0,1,1,0,1,0/
|
|
|
|
nnn=0
|
|
nargs=iargc()
|
|
if(nargs.ne.6) then
|
|
print*,'Usage: dbpsksim f0(Hz) delay(ms) fspread(Hz) ndiff iters snr(dB)'
|
|
print*,'Example: dbpsksim 1500 0 0 10 -35'
|
|
print*,'Set snr=0 to cycle through a range'
|
|
go to 999
|
|
endif
|
|
call getarg(1,arg)
|
|
read(arg,*) f0 !Low tone frequency
|
|
call getarg(2,arg)
|
|
read(arg,*) delay
|
|
call getarg(3,arg)
|
|
read(arg,*) fspread
|
|
call getarg(4,arg)
|
|
read(arg,*) ndiff
|
|
call getarg(5,arg)
|
|
read(arg,*) iters
|
|
call getarg(6,arg)
|
|
read(arg,*) snrdb
|
|
|
|
twopi=8.d0*atan(1.d0)
|
|
fs=12000.d0
|
|
dt=1.0/fs
|
|
ts=NSPS*dt
|
|
baud=1.d0/ts
|
|
txt=NZ*dt
|
|
bandwidth_ratio=2500.0/6000.0
|
|
write(*,1000) baud,5*baud,txt,delay,fspread,ndiff
|
|
1000 format('Baud:',f6.3,' BW:',f4.1,' TxT:',f6.1,' Delay:',f5.2, &
|
|
' fSpread:',f5.2,' ndiff:',i2/)
|
|
|
|
write(*,1004)
|
|
1004 format(' SNR err ber fer fsigma'/35('-'))
|
|
|
|
call encode120(msgbits,codeword) !Encode the test message
|
|
isna=-28
|
|
isnb=-40
|
|
if(snrdb.ne.0.0) then
|
|
isna=nint(snrdb)
|
|
isnb=isna
|
|
endif
|
|
do isnr=isna,isnb,-1
|
|
snrdb=isnr
|
|
sig=sqrt(bandwidth_ratio) * 10.0**(0.05*snrdb)
|
|
if(snrdb.gt.90.0) sig=1.0
|
|
nhard=0
|
|
nhardc=0
|
|
nfe1=0
|
|
nfe2=0
|
|
sqf=0.
|
|
do iter=1,iters
|
|
nnn=nnn+1
|
|
id(1)=1 !First bit is always 1
|
|
id(2:NN)=2*codeword-1
|
|
call genbpsk(id,f0,ndiff,0,c) !Generate the 4-FSK waveform
|
|
if(delay.ne.0.0 .or. fspread.ne.0.0) call watterson(c,delay,fspread)
|
|
c=sig*c !Scale to requested SNR
|
|
if(snrdb.lt.90) then
|
|
do i=0,NZ-1 !Generate gaussian noise
|
|
xnoise(i)=gran()
|
|
ynoise(i)=gran()
|
|
enddo
|
|
c=c + cmplx(xnoise,ynoise) !Add noise to signal
|
|
endif
|
|
|
|
! First attempt at finding carrier frequency fc: 64k FFTs ==> avg power spectra
|
|
nspec=NZ/NFFT1
|
|
df1=12000.0/NFFT1
|
|
s=0.
|
|
do k=1,nspec
|
|
ia=(k-1)*NSPS
|
|
ib=ia+NSPS-1
|
|
c2(0:NSPS-1)=c(ia:ib)
|
|
c2(NSPS:)=0.
|
|
call four2a(c2,NFFT1,1,-1,1)
|
|
do i=0,NFFT1-1
|
|
j=i
|
|
if(j.gt.NH1) j=j-NFFT1
|
|
s(j)=s(j) + real(c2(i))**2 + aimag(c2(i))**2
|
|
enddo
|
|
enddo
|
|
s=1.e-6*s
|
|
smax=0.
|
|
ipk=0
|
|
ia=(1400.0)/df1
|
|
ib=(1600.0)/df1
|
|
do i=ia,ib
|
|
f=i*df1
|
|
if(s(i).gt.smax) then
|
|
smax=s(i)
|
|
ipk=i
|
|
fc=f
|
|
endif
|
|
enddo
|
|
a=(s(ipk+1)-s(ipk-1))/2.0
|
|
b=(s(ipk+1)+s(ipk-1)-2.0*s(ipk))/2.0
|
|
dx=-a/(2.0*b)
|
|
fc=fc + df1*dx !Estimated carrier frequency
|
|
sqf=sqf + (fc-f0)**2
|
|
|
|
! The following is for testing SNR calibration:
|
|
! sp5n=(s(ipk-2)+s(ipk-1)+s(ipk)+s(ipk+1)+s(ipk+2)) !Sig + 5*noise
|
|
! base=(sum(s)-sp5n)/(NFFT1-5.0) !Noise per bin
|
|
! psig=sp5n-5*base !Sig only
|
|
! pnoise=(2500.0/df1)*base !Noise in 2500 Hz
|
|
! xsnrdb=db(psig/pnoise)
|
|
|
|
call genbpsk(id,fc,ndiff,1,cr) !Generate reference carrier
|
|
c=c*conjg(cr) !Mix signal to baseband
|
|
|
|
z0=1.0
|
|
do j=1,NN !Demodulate
|
|
ia=(j-1)*NSPS
|
|
ib=ia+NSPS-1
|
|
z=sum(c(ia:ib))
|
|
cz(j-1)=z
|
|
zp=z*conjg(z0)
|
|
p=1.e-4*real(zp)
|
|
id1(j)=-1
|
|
if(p.ge.0.0) id1(j)=1
|
|
if(j.ge.2) rxdata(j-1)=p
|
|
z0=z
|
|
enddo
|
|
|
|
rxav=sum(rxdata)/120
|
|
rx2av=sum(rxdata*rxdata)/120
|
|
rxsig=sqrt(rx2av-rxav*rxav)
|
|
rxdata=rxdata/rxsig
|
|
ss=0.84
|
|
llr=2.0*rxdata/(ss*ss)
|
|
apmask=0
|
|
max_iterations=10
|
|
call bpdecode120(llr,apmask,max_iterations,decoded,niterations,cw)
|
|
|
|
! Count frame errors
|
|
if(niterations.lt.0 .or. count(msgbits.ne.decoded).gt.0) nfe1=nfe1+1
|
|
|
|
! Find carrier frequency from squared cz array.
|
|
cz(121:)=0.
|
|
cz=cz*cz
|
|
call four2a(cz,NFFT2,1,-1,1)
|
|
s2max=0.
|
|
do i=0,NFFT2-1
|
|
j=i
|
|
if(i.gt.NH2) j=j-NFFT2
|
|
s2(j)=real(cz(i))**2 + aimag(cz(i))**2
|
|
if(s2(j).gt.s2max) then
|
|
s2max=s2(j)
|
|
jpk=j
|
|
endif
|
|
! write(16,1200) j*baud/NFFT2,1.e-12*s2(j)
|
|
!1200 format(2f12.3)
|
|
enddo
|
|
a=(s2(jpk+1)-s2(jpk-1))/2.0
|
|
b=(s2(jpk+1)+s2(jpk-1)-2.0*s2(jpk))/2.0
|
|
dx=-a/(2.0*b)
|
|
fc2=0.5*(jpk+dx)*baud/NFFT2
|
|
|
|
call genbpsk(id,fc2,ndiff,1,cr) !Generate new ref carrier at fc2
|
|
c=c*conjg(cr)
|
|
z0=1.0
|
|
do j=1,NN !Demodulate
|
|
ia=(j-1)*NSPS
|
|
ib=ia+NSPS-1
|
|
z=sum(c(ia:ib))
|
|
if(j.eq.1) z0=z
|
|
zp=z*conjg(z0)
|
|
p=1.e-4*real(zp)
|
|
id2(j)=-1
|
|
if(p.ge.0.0) id2(j)=1
|
|
if(j.ge.2) rxdata(j-1)=p
|
|
ierr=0
|
|
if(id2(j).ne.id(j)) ierr=1
|
|
id3=-1
|
|
if(real(z).ge.0.0) id3=1
|
|
if(j.ge.2 .and. id3.ne.id(j)) nhardc=nhardc+1
|
|
if(j.ge.2 .and. ndiff.eq.0) rxdata(j-1)=real(z)
|
|
z0=z
|
|
enddo
|
|
nhard=nhard + count(id2.ne.id) !Count hard errors
|
|
|
|
rxav=sum(rxdata)/120
|
|
rx2av=sum(rxdata*rxdata)/120
|
|
rxsig=sqrt(rx2av-rxav*rxav)
|
|
rxdata=rxdata/rxsig
|
|
ss=0.84
|
|
llr=2.0*rxdata/(ss*ss) !Soft symbols
|
|
apmask=0
|
|
max_iterations=10
|
|
decoded=0
|
|
call bpdecode120(llr,apmask,max_iterations,decoded,niterations,cw)
|
|
! if(niterations.lt.0) then
|
|
! llr=-llr
|
|
! call bpdecode120(llr,apmask,max_iterations,decoded,niterations,cw)
|
|
! if(niterations.ge.0) nhard=NN*iters-nhard
|
|
! endif
|
|
if(niterations.ge.0) call chkcrc10(decoded,nbadcrc)
|
|
if(niterations.lt.0 .or. count(msgbits.ne.decoded).gt.0 .or. &
|
|
nbadcrc.ne.0) nfe2=nfe2+1
|
|
enddo
|
|
|
|
if(ndiff.eq.0) nhard=nhardc
|
|
fsigma=sqrt(sqf/iters)
|
|
ber=float(nhard)/(NN*iters)
|
|
fer=float(nfe2)/iters
|
|
write(*,1050) snrdb,nhard,ber,fer,fsigma
|
|
write(14,1050) snrdb,nhard,ber,fer,fsigma
|
|
1050 format(f6.1,i5,f8.4,f7.3,f8.2)
|
|
enddo
|
|
|
|
999 end program dbpsksim
|