w4kek
/
WSJT-X
mirror of https://github.com/saitohirga/WSJT-X.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Save fsk4sim.f90 with 3 @ 4x4 sync arrays. 

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7625 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Joe Taylor 2017-03-30 19:24:20 +00:00
parent 1298b13957
commit 9ebd755c37
2 changed files with 1794 additions and 2242 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

423
lib/fsk4hf/fsk4sim.f90
View File

@ -1,14 +1,17 @@
program fsk4sim
use wavhdr
parameter (NR=4) !Ramp up, ramp down
parameter (NS=12) !Sync symbols (2 @ Costas 4x4)
parameter (ND=84) !Data symbols: LDPC (168,84), r=1/2
parameter (NN=NR+NS+ND) !Total symbols (100)
parameter (NSPS=2688) !Samples per symbol at 12000 sps
parameter (NZ=NSPS*NN) !Samples in waveform (258048)
parameter (NR=4) !Ramp up, ramp down
parameter (NS=12) !Sync symbols (2 @ Costas 4x4)
parameter (ND=84) !Data symbols: LDPC (168,84), r=1/2
parameter (NN=NR+NS+ND) !Total symbols (100)
parameter (NSPS=2688) !Samples per symbol at 12000 sps
parameter (NZ=NSPS*NN) !Samples in waveform (268800)
parameter (NSYNC=NS*NSPS) !Samples in sync waveform (32256)
parameter (NFFT=512*1024)
parameter (NSYNC=NS*NSPS)
parameter (NDOWN=84) !Downsample factor
parameter (NFFT2=NZ/NDOWN,NH2=NFFT2/2) !3200
parameter (NSPSD=NFFT2/NN) !Samples per symbol after downsample
type(hdr) header !Header for .wav file
character*8 arg
@ -17,30 +20,42 @@ program fsk4sim
complex cs(0:NSYNC-1)
complex ct(0:NSPS-1)
complex csync(0:NSYNC-1)
complex z
complex c2(0:NFFT2-1)
complex c2a(0:NSPSD-1)
complex cf2(0:NFFT2-1)
complex cx(0:3,NN)
complex z,zpk
logical snrtest
real*8 twopi,dt,fs,baud,f0,dphi,phi
real tmp(NN) !For generating random data
real s(0:NFFT-1)
! real s2(0:NFFT-1)
real xnoise(NZ) !Generated random noise
real xnoise(0:NZ-1) !Generated random noise
real s(NSYNC/2)
real ps(0:3)
integer*2 iwave(NZ) !Generated waveform
! integer*2 iwave(NZ) !Generated waveform
integer id(NN) !Encoded 2-bit data (values 0-3)
integer id2(NN) !Decoded after downsampling
integer icos4(4) !4x4 Costas array
data icos4/0,1,3,2/
data icos4/0,1,3,2/,eps/1.e-8/
nargs=iargc()
if(nargs.ne.3) then
print*,'Usage: fsk8sim f0 iters snr'
if(nargs.ne.4) then
print*,'Usage: fsk8sim f0 fspread iters snr'
go to 999
endif
call getarg(1,arg)
read(arg,*) f0 !Low tone frequency
call getarg(2,arg)
read(arg,*) iters
read(arg,*) fspread
call getarg(3,arg)
read(arg,*) iters
call getarg(4,arg)
read(arg,*) snrdb
snrtest=.false.
if(iters.lt.0) then
snrtest=.true.
iters=abs(iters)
endif
twopi=8.d0*atan(1.d0)
fs=12000.d0
@ -48,43 +63,13 @@ program fsk4sim
ts=NSPS*dt
baud=1.d0/ts
txt=NZ*dt
isna=-20
isnb=-30
if(snrdb.ne.0.0) then
isna=nint(snrdb)
isnb=isna
endif
do isnr=isna,isnb,-1
snrdb=isnr
bandwidth_ratio=2500.0/6000.0
sig=sqrt(2*bandwidth_ratio) * 10.0**(0.05*snrdb)
if(snrdb.gt.90.0) sig=1.0
header=default_header(12000,NZ)
open(10,file='000000_0001.wav',access='stream',status='unknown')
nsyncerr=0
nharderr=0
nbiterr=0
do iter=1,iters
id=0
call random_number(tmp)
where(tmp.ge.0.25 .and. tmp.lt.0.50) id=1
where(tmp.ge.0.50 .and. tmp.lt.0.75) id=2
where(tmp.ge.0.75) id=3
id(1:2)=icos4(3:4) !Ramp up
id(45:48)=icos4 !Costas sync
id(49:52)=icos4 !Costas sync
id(53:56)=icos4 !Costas sync
id(NN-1:NN)=icos4(1:2) !Ramp down
! Generate sync waveform
phi=0.d0
k=-1
do j=45,56
dphi=twopi*(id(j)*baud)*dt
do j=1,12
n=mod(j-1,4) + 1
dphi=twopi*(icos4(n)*baud)*dt
do i=1,NSPS
k=k+1
phi=phi+dphi
@ -93,141 +78,229 @@ program fsk4sim
csync(k)=cmplx(cos(xphi),-sin(xphi))
enddo
enddo
bandwidth_ratio=2500.0/6000.0
header=default_header(12000,NZ)
! Generate the 4-FSK waveform
x=0.
c=0.
phi=0.d0
k=-1
u=0.5
do j=1,NN
dphi=twopi*(f0 + id(j)*baud)*dt
do i=1,NSPS
k=k+1
phi=phi+dphi
if(phi.gt.twopi) phi=phi-twopi
xphi=phi
c(k)=cmplx(cos(xphi),sin(xphi))
enddo
enddo
if(sig.ne.1.0) c=sig*c
nh=NFFT/2
df=12000.0/NFFT
s=0.
cf=c
call four2a(cf,NFFT,1,-1,1) !Transform to frequency domain
flo=f0-baud
fhi=f0+4*baud
do i=0,NFFT-1 !Remove spectral sidelobes
f=i*df
if(i.gt.nh) f=(i-nfft)*df
if(f.le.flo .or. f.ge.fhi) cf(i)=0.
s(i)=s(i) + real(cf(i))**2 + aimag(cf(i))**2
enddo
! s2=cshift(s,nh)
! s2=s2/maxval(s2)
! do i=0,NFFT-1
! f=(i-nh)*df
! write(13,1000) f,s2(i),10.0*log10(s2(i)+1.e-12)
!1000 format(3f12.3)
! enddo
c=cf
call four2a(c,NFFT,1,1,1) !Transform back to time domain
c=c/nfft
xnoise=0.
if(snrdb.lt.90) then
a=1.0/sqrt(2.0)
do i=0,NZ-1
xx=a*gran()
yy=a*gran()
c(i)=c(i) + cmplx(xx,yy) !Scale signal and add noise
enddo
write(*,1000) 2*ND,ND,NS,NN,NSPS,baud,txt,fspread
1000 format('LDPC('i3,',',i2,') SyncSym:',i2,' ChanSym:',i3,' NSPS:',i4, &
' Baud:',f6.3,' TxT:',f5.1,' fDop:',f5.2/)
if(snrtest) then
write(*,1002)
1002 format(5x,'SNR test'/'Requested Measured Difference')
else
write(*,1004)
1004 format(' SNR Sync Sym1 Sym2 Bits SyncErr Sym1Err BER'/ &
60('-'))
endif
fac=32767.0
rms=100.0
if(snrdb.ge.90.0) iwave(1:NZ)=nint(fac*aimag(c(0:NZ-1)))
if(snrdb.lt.90.0) iwave(1:NZ)=nint(rms*aimag(c(0:NZ-1)))
call set_wsjtx_wav_params(14.0,'JT65 ',1,30,iwave)
write(10) header,iwave !Save the .wav file
isna=-15
isnb=-27
if(snrdb.ne.0.0) then
isna=nint(snrdb)
isnb=isna
endif
do isnr=isna,isnb,-1
snrdb=isnr
sig=sqrt(2*bandwidth_ratio) * 10.0**(0.05*snrdb)
if(snrdb.gt.90.0) sig=1.0
! open(10,file='000000_0001.wav',access='stream',status='unknown')
! do i=0,NZ-1
! a=abs(c(i))
! j=mod(i,NSPS)
! write(14,1010) i*dt/ts,c(i),a
!1010 format(4f12.6)
! enddo
ppmax=0.
fpk=-99.
xdt=-99.
do j4=-40,40
ia=(44+0.25*j4)*NSPS
ib=ia+NSYNC-1
cs=csync*c(ia:ib)
call four2a(cs,NSYNC,1,-1,1) !Transform to frequency domain
df1=12000.0/NSYNC
fac=1.e-6
do i=0,NSYNC/2
pp=fac*(real(cs(i))**2 + aimag(cs(i))**2)
if(pp.gt.ppmax) then
fpk=i*df1
xdt=0.25*j4*ts
ppmax=pp
nsyncerr=0
nharderr=0
nherr=0
nbiterr=0
do iter=1,iters
id=0
if(.not.snrtest) then
! Generate random data
call random_number(tmp)
where(tmp.ge.0.25 .and. tmp.lt.0.50) id=1
where(tmp.ge.0.50 .and. tmp.lt.0.75) id=2
where(tmp.ge.0.75) id=3
id(1:2)=icos4(3:4) !Ramp up
id(45:48)=icos4 !Costas sync
id(49:52)=icos4 !Costas sync
id(53:56)=icos4 !Costas sync
id(NN-1:NN)=icos4(1:2) !Ramp down
endif
! if(j4.eq.0) then
! f=i*df1
! write(16,1030) f,pp,10.0*log10(pp)
!1030 format(3f15.3)
! endif
enddo
enddo
if(xdt.ne.0.0 .or. fpk.ne.1500.0) nsyncerr=nsyncerr+1
ipk=0
do j=1,NN
ia=(j-1)*NSPS + 1
ib=ia+NSPS
pmax=0.
do i=0,3
f=fpk + i*baud
call tweak1(c(ia:ib),NSPS,-f,ct)
z=sum(ct)
ps(i)=1.e-3*(real(z)**2 + aimag(z)**2)
if(ps(i).gt.pmax) then
ipk=i
pmax=ps(i)
call genfsk4(id,f0,c) !Generate the 4-FSK waveform
if(sig.ne.1.0) c=sig*c !Scale to requested SNR
if(snrdb.lt.90) then
do i=0,NZ-1 !Generate gaussian noise
xnoise(i)=gran()
enddo
endif
if(fspread.gt.0.0) call dopspread(c,fspread)
c(0:NZ-1)=real(c(0:NZ-1)) + xnoise !Add noise to signal
! fac=32767.0
! rms=100.0
! if(snrdb.ge.90.0) iwave(1:NZ)=nint(fac*aimag(c(0:NZ-1)))
! if(snrdb.lt.90.0) iwave(1:NZ)=nint(rms*aimag(c(0:NZ-1)))
! call set_wsjtx_wav_params(14.0,'JT65 ',1,30,iwave)
! write(10) header,iwave !Save the .wav file
ppmax=0.
fpk=-99.
xdt=-99.
df1=12000.0/NSYNC
iaa=nint(250.0/df1)
ibb=nint(2750.0/df1)
if(.not.snrtest) then
do j4=-40,40
ia=(44+0.25*j4)*NSPS
ib=ia+NSYNC-1
cs=csync*c(ia:ib)
call four2a(cs,NSYNC,1,-1,1) !Transform to frequency domain
s=0.
do i=iaa,ibb
s(i)=1.e-6*(real(cs(i))**2 + aimag(cs(i))**2)
enddo
if(j4.eq.0) then
do i=iaa,ibb
write(66,3301) i*df1,s(i)
3301 format(f10.3,2f12.6)
enddo
endif
call smo121(s,NSYNC/2)
if(j4.eq.0) then
do i=iaa,ibb
write(67,3301) i*df1,s(i)
enddo
endif
do i=iaa,ibb
if(s(i).gt.ppmax) then
fpk=i*df1
xdt=0.25*j4*ts
ppmax=s(i)
endif
enddo
enddo
endif
if(xdt.ne.0.0 .or. fpk.ne.1500.0) nsyncerr=nsyncerr+1
! Compute spectrum again
cf=c
df2=12000.0/NZ
call four2a(cf,NZ,1,-1,1) !Transform to frequency domain
if(snrtest) then
width=5.0*df2 + fspread
iz=nint(2500.0/df2) + 2
if(iter.eq.1) then
pnoise=0.
psig=0.
n=0
endif
do i=0,iz !Remove spectral sidelobes
f=i*df2
if(i.gt.NZ/2) f=(i-NZ)*df2
p=1.e-6*(real(cf(i))**2 + aimag(cf(i))**2)
if(abs(f-f0).lt.width) then
psig=psig+p
n=n+1
else
pnoise=pnoise + p
endif
enddo
if(iter.eq.iters) then
db=10.0*log10(psig/pnoise)
write(*,1010) snrdb,db,db-snrdb
1010 format(f7.1,2f9.1)
endif
go to 40
endif
! Select a small frequency slice around fpk.
cf=cf/NZ
ib=nint(fpk/df2)+NH2
ia=ib-NFFT2+1
cf2=cshift(cf(ia:ib),NH2-1)
flo=-baud
fhi=4*baud
do i=0,NFFT2-1
f=i*df2
if(i.gt.NH2) f=(i-NFFT2)*df2
if(f.le.flo .or. f.ge.fhi) cf2(i)=0.
s2=real(cf2(i))**2 + aimag(cf2(i))**2
write(15,3001) f,s2,10*log10(s2+eps)
3001 format(f10.3,2f15.6)
enddo
c2=cf2
call four2a(c2,NFFT2,1,1,1) !Back to time domain
fshift=NSPS*baud/NSPSD
dt2=dt*NDOWN
do j=1,NN
ia=(j-1)*NSPSD
ib=ia+NSPSD-1
c2a=c2(ia:ib)
call four2a(c2a,NSPSD,1,-1,1) !To freq domain
cx(0:3,j)=c2a(0:3)
ipk=-1
zpk=0.
pmax=0.0
do i=0,3
if(abs(cx(i,j)).gt.pmax) then
ipk=i
zpk=cx(i,j)
pmax=abs(zpk)
endif
enddo
id2(j)=ipk
if(ipk.ne.id(j)) nherr=nherr+1
write(16,3003) j,id(j),ipk,ipk-id(j),abs(zpk), &
atan2(aimag(zpk),real(zpk)),abs(cx(0:3,j))
3003 format(3i3,i4,6f9.3)
enddo
ipk=0
do j=1,NN
ia=(j-1)*NSPS + 1
ib=ia+NSPS
pmax=0.
do i=0,3
f=fpk + i*baud
call tweak1(c(ia:ib),NSPS,-f,ct)
z=sum(ct)
ps(i)=1.e-3*(real(z)**2 + aimag(z)**2)
if(ps(i).gt.pmax) then
ipk=i
pmax=ps(i)
endif
enddo
nlo=0
nhi=0
if(max(ps(1),ps(3)).ge.max(ps(0),ps(2))) nlo=1
if(max(ps(2),ps(3)).ge.max(ps(0),ps(1))) nhi=1
if(nlo.ne.iand(id(j),1)) nbiterr=nbiterr+1
if(nhi.ne.iand(id(j)/2,1)) nbiterr=nbiterr+1
if(ipk.ne.id(j)) nharderr=nharderr+1
write(17,1040) j,ps,ipk,id(j),id2(j),2*nhi+nlo,nhi,nlo,nbiterr
1040 format(i3,4f12.1,7i4)
enddo
40 continue
enddo
nlo=0
nhi=0
if(max(ps(1),ps(3)).ge.max(ps(0),ps(2))) nlo=1
if(max(ps(2),ps(3)).ge.max(ps(0),ps(1))) nhi=1
! if(ps(1)+ps(3).ge.ps(0)+ps(2)) nlo=1
! if(ps(2)+ps(3).ge.ps(0)+ps(1)) nhi=1
if(nlo.ne.iand(id(j),1)) nbiterr=nbiterr+1
if(nhi.ne.iand(id(j)/2,1)) nbiterr=nbiterr+1
if(ipk.ne.id(j)) nharderr=nharderr+1
write(17,1040) j,ps,ipk,id(j),2*nhi+nlo,nhi,nlo,nbiterr
1040 format(i3,4f12.1,6i4)
enddo
enddo
fsyncerr=float(nsyncerr)/iters
ser=float(nharderr)/(NN*iters)
ber=float(nbiterr)/(2*NN*iters)
write(*,1050) snrdb,nsyncerr,nharderr,nbiterr,fsyncerr,ser,ber
write(18,1050) snrdb,nsyncerr,nharderr,nbiterr,fsyncerr,ser,ber
1050 format(f6.1,3i6,3f10.6)
if(.not.snrtest) then
fsyncerr=float(nsyncerr)/iters
ser=float(nharderr)/(NN*iters)
ber=float(nbiterr)/(2*NN*iters)
write(*,1050) snrdb,nsyncerr,nharderr,nherr,nbiterr,fsyncerr,ser,ber
write(18,1050) snrdb,nsyncerr,nharderr,nherr,nbiterr,fsyncerr,ser,ber
1050 format(f6.1,4i6,3f10.6)
endif
enddo
if(.not.snrtest) write(*,1060) iters,100*iters,100*iters,200*iters
1060 format(60('-')/'Max: ',4i6)
999 end program fsk4sim

3613
lib/nfft.out
View File

File diff suppressed because it is too large Load Diff