diff --git a/CMakeLists.txt b/CMakeLists.txt
index 4e5cabaa4..722ab1b3f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -385,6 +385,7 @@ set (wsjt_FSRCS
   lib/demod64a.f90
   lib/determ.f90
   lib/downsam9.f90
+  lib/echosim.f90
   lib/echo_snr.f90
   lib/encode232.f90
   lib/encode4.f90
@@ -1184,6 +1185,9 @@ target_link_libraries (ft8code wsjt_fort wsjt_cxx)
 add_executable (ft4code lib/ft4/ft4code.f90)
 target_link_libraries (ft4code wsjt_fort wsjt_cxx)
 
+add_executable (echosim lib/echosim.f90)
+target_link_libraries (echosim wsjt_fort wsjt_cxx)
+
 add_executable (ft8sim lib/ft8/ft8sim.f90)
 target_link_libraries (ft8sim wsjt_fort wsjt_cxx)
 
diff --git a/lib/echosim.f90 b/lib/echosim.f90
new file mode 100644
index 000000000..b90a5f541
--- /dev/null
+++ b/lib/echosim.f90
@@ -0,0 +1,101 @@
+program echosim
+
+! Generate simulated echo-mode files -- self-echo or "measure" 
+
+  use wavhdr
+  parameter (NWAVE=27648,NMAX=32768,NZ=36000)
+  type(hdr) h                            !Header for .wav file
+  character arg*12,fname*17
+  complex c0(0:NMAX-1)
+  complex c(0:NMAX-1)
+!  complex cwave(0:NWAVE-1)
+  real*8 f0,dt,twopi,phi,dphi
+  real wave(NZ)
+  integer*2 iwave(NZ)                  !Generated full-length waveform
+  equivalence (nDop0,iwave(1))
+  equivalence (nDopAudio0,iwave(3))
+  equivalence (nfrit0,iwave(5))
+  equivalence (f10,iwave(7))
+  equivalence (fspread0,iwave(9))
+
+! 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'
+     go to 999
+  endif
+
+  call getarg(1,arg)
+  read(arg,*) f0                         !Tone frequency
+  call getarg(2,arg)
+  read(arg,*) fdop                       !Doppler shift (Hz)
+  call getarg(3,arg)
+  read(arg,*) fspread                    !Watterson frequency spread (Hz)
+  call getarg(4,arg)
+  read(arg,*) delay                      !Watterson delay (ms)
+  call getarg(5,arg)
+  read(arg,*) nfiles                     !Number of files
+  call getarg(6,arg)
+  read(arg,*) snrdb                      !SNR_2500
+
+  twopi=8.d0*atan(1.d0)
+  fs=12000.0                             !Sample rate (Hz)
+  dt=1.d0/fs                              !Sample interval (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
+  dphi=twopi*(f0+fdop)*dt
+
+  do ifile=1,nfiles
+     phi=0.d0
+     do i=0,NWAVE-1
+        phi=phi + dphi
+        if(phi.gt.twopi) phi=phi-twopi
+        xphi=phi
+        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)
+     c=sig*c0
+     wave(1:NWAVE)=imag(c(1:NWAVE))
+     peak=maxval(abs(wave))
+     if(snrdb.lt.90) then
+        do i=1,NWAVE                   !Add gaussian noise at specified SNR
+           xnoise=gran()
+           wave(i)=wave(i) + xnoise
+        enddo
+        do i=NWAVE+1,NZ
+           xnoise=gran()
+           wave(i)=xnoise
+        enddo
+     endif
+
+     gain=100.0
+     if(snrdb.lt.90.0) then
+       wave=gain*wave
+     else
+       datpk=maxval(abs(wave))
+       fac=32766.9/datpk
+       wave=fac*wave
+     endif
+     if(any(abs(wave).gt.32767.0)) print*,"Warning - data will be clipped."
+     iwave=nint(wave)
+
+     nDop0=nint(fdop)
+     nDopAudio0=0
+     nfrit0=0
+     f10=f0 + fdop
+     fspread0=fspread
+     
+     h=default_header(12000,NMAX)
+     write(fname,1102) 3*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 echosim