From c81bcfa3ef9d5d45356b93177445f4e25784b8a4 Mon Sep 17 00:00:00 2001
From: Joe Taylor <joe@princeton.edu>
Date: Tue, 30 Aug 2022 14:57:44 -0400
Subject: [PATCH] Modify echosim to use the same Lorentzian model for fspread
as used in q65sim.
---
CMakeLists.txt | 1 +
lib/avecho.f90 | 8 +++++--
lib/echosim.f90 | 14 ++++++------
lib/fspread_lorentz.f90 | 47 +++++++++++++++++++++++++++++++++++++++++
4 files changed, 60 insertions(+), 10 deletions(-)
create mode 100644 lib/fspread_lorentz.f90
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 722ab1b3f..2341257d3 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -420,6 +420,7 @@ set (wsjt_FSRCS
lib/fmtmsg.f90
lib/foldspec9f.f90
lib/four2a.f90
+ lib/fspread_lorentz.f90
lib/ft8/foxfilt.f90
lib/ft8/foxgen.f90
lib/ft8/foxgen_wrap.f90
diff --git a/lib/avecho.f90 b/lib/avecho.f90
index 69b006e41..476e4321e 100644
--- a/lib/avecho.f90
+++ b/lib/avecho.f90
@@ -87,8 +87,12 @@ subroutine avecho(id2,ndop,nfrit,nauto,nqual,f1,xlevel,snrdb,db_err, &
call smo121(blue,NZ)
enddo
-! write(*,3001) snrdb,db_err,dfreq,snr_detect,redmax,nqual,nsmo,nclearave,nsum
-!3001 format('A',5f10.1,4i4)
+ ia=200.0/df
+ ib=400.0/df
+ call pctile(red(ia:ib),ib-ia+1,50,bred)
+ red=red-bred
+ call pctile(blue(ia:ib),ib-ia+1,50,bblue)
+ blue=blue-bblue
900 return
end subroutine avecho
diff --git a/lib/echosim.f90 b/lib/echosim.f90
index b90a5f541..0cf193de0 100644
--- a/lib/echosim.f90
+++ b/lib/echosim.f90
@@ -20,9 +20,9 @@ program echosim
! Get command-line argument(s)
nargs=iargc()
- if(nargs.ne.6) then
- print*,'Usage: echosim f0 fdop fspread delay nfiles snr'
- print*,'Examples: echosim 1500 0.0 4.0 5.0 10 -22'
+ if(nargs.ne.5) then
+ print*,'Usage: echosim f0 fdop fspread nfiles snr'
+ print*,'Examples: echosim 1500 0.0 4.0 10 -22'
go to 999
endif
@@ -31,12 +31,10 @@ program echosim
call getarg(2,arg)
read(arg,*) fdop !Doppler shift (Hz)
call getarg(3,arg)
- read(arg,*) fspread !Watterson frequency spread (Hz)
+ read(arg,*) fspread !Frequency spread (Hz) (JHT Lorentzian model)
call getarg(4,arg)
- read(arg,*) delay !Watterson delay (ms)
- call getarg(5,arg)
read(arg,*) nfiles !Number of files
- call getarg(6,arg)
+ call getarg(5,arg)
read(arg,*) snrdb !SNR_2500
twopi=8.d0*atan(1.d0)
@@ -56,7 +54,7 @@ program echosim
c0(i)=cmplx(cos(xphi),sin(xphi))
enddo
c0(NWAVE:)=0.
- if(fspread.ne.0.0 .or. delay.ne.0.0) call watterson(c0,NMAX,NWAVE,fs,delay,fspread)
+ if(fspread.gt.0.0) call fspread_lorentz(c0,fspread)
c=sig*c0
wave(1:NWAVE)=imag(c(1:NWAVE))
peak=maxval(abs(wave))
diff --git a/lib/fspread_lorentz.f90 b/lib/fspread_lorentz.f90
new file mode 100644
index 000000000..697d308c5
--- /dev/null
+++ b/lib/fspread_lorentz.f90
@@ -0,0 +1,47 @@
+subroutine fspread_lorentz(cdat,fspread)
+
+ parameter (NZ=3*12000)
+ complex cdat(0:NZ-1)
+ complex cspread(0:NZ-1)
+ complex z
+
+ twopi=8.0*atan(1.0)
+ nfft=NZ
+ nh=nfft/2
+ df=12000.0/nfft
+ cspread(0)=1.0
+ cspread(nh)=0.
+ b=6.0 !Use truncated Lorenzian shape for fspread
+ do i=1,nh
+ f=i*df
+ x=b*f/fspread
+ z=0.
+ a=0.
+ if(x.lt.3.0) then !Cutoff beyond x=3
+ a=sqrt(1.111/(1.0+x*x)-0.1) !Lorentzian amplitude
+ phi1=twopi*rran() !Random phase
+ z=a*cmplx(cos(phi1),sin(phi1))
+ endif
+ cspread(i)=z
+ z=0.
+ if(x.lt.3.0) then !Same thing for negative freqs
+ phi2=twopi*rran()
+ z=a*cmplx(cos(phi2),sin(phi2))
+ endif
+ cspread(nfft-i)=z
+ enddo
+
+ call four2a(cspread,nfft,1,1,1) !Transform to time domain
+
+ sum=0.
+ do i=0,nfft-1
+ p=real(cspread(i))**2 + aimag(cspread(i))**2
+ sum=sum+p
+ enddo
+ avep=sum/nfft
+ fac=sqrt(1.0/avep)
+ cspread=fac*cspread !Normalize to constant avg power
+ cdat=cspread*cdat !Apply Rayleigh fading
+
+ return
+end subroutine fspread_lorentz