diff --git a/lib/fsk4hf/Makefile.win b/lib/fsk4hf/Makefile.win
new file mode 100644
index 000000000..e2687eac5
--- /dev/null
+++ b/lib/fsk4hf/Makefile.win
@@ -0,0 +1,47 @@
+# Compilers
+CC = gcc
+CXX = g++
+FC = gfortran
+
+FFLAGS = -O2 -fbounds-check -Wall -Wno-conversion
+CFLAGS = -O2 -I. 
+
+# Default rules
+%.o: %.c
+	${CC} ${CFLAGS} -c $<
+%.o: %.f
+	${FC} ${FFLAGS} -c $<
+%.o: %.F
+	${FC} ${FFLAGS} -c $<
+%.o: %.f90
+	${FC} ${FFLAGS} -c $<
+%.o: %.F90
+	${FC} ${FFLAGS} -c $<
+
+all:	fsk4sim.exe
+
+OBJS0 = testpsk.o four2a.o bpfilter.o nonlinear.o tweak1.o spectrum.o smo.o
+testpsk: $(OBJS0)
+	$(FC) -o testpsk $(OBJS0) C:\JTSDK\fftw3f\libfftw3f-3.dll
+
+OBJS1 = gmsk8.o four2a.o gaussfilt.o
+gmsk8: $(OBJS1)
+	$(FC) -o gmsk8 $(OBJS1) C:\JTSDK\fftw3f\libfftw3f-3.dll
+
+OBJS2 = testfsk.o four2a.o smo.o
+testfsk: $(OBJS2)
+	$(FC) -o testfsk $(OBJS2) C:\JTSDK\fftw3f\libfftw3f-3.dll
+
+OBJS3 = fsk2sim.o four2a.o smo.o wavhdr.o gran.o 
+fsk2sim: $(OBJS3)
+	$(FC) -o fsk2sim $(OBJS3) C:\JTSDK\fftw3f\libfftw3f-3.dll
+
+OBJS4 = fsk4sim.o four2a.o wavhdr.o gran.o tweak1.o genfsk4.o \
+      dopspread.o smo.o smo121.o
+fsk4sim.exe: $(OBJS4)
+	$(FC) -o fsk4sim.exe $(OBJS4) C:\JTSDK\fftw3f\libfftw3f-3.dll
+
+.PHONY : clean
+
+clean:
+	$(RM) *.o testpsk.exe testfsk.exe fsk2sim.exe fsk4sim.exe
diff --git a/lib/fsk4hf/dopspread.f90 b/lib/fsk4hf/dopspread.f90
new file mode 100644
index 000000000..cadfbecc6
--- /dev/null
+++ b/lib/fsk4hf/dopspread.f90
@@ -0,0 +1,62 @@
+subroutine dopspread(c,fspread)
+
+  parameter (NFFT=268800,NH=NFFT/2)
+  complex c(0:NFFT-1)
+  complex cspread(0:NFFT-1)
+
+  df=12000.0/nfft
+  twopi=8*atan(1.0)
+  cspread(0)=1.0
+  cspread(NH)=0.
+  b=6.0                                     !Lorenzian 3/28 onward
+  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
+        call random_number(r1)
+        phi1=twopi*r1
+        z=a*cmplx(cos(phi1),sin(phi1))
+     endif
+     cspread(i)=z
+     z=0.
+     if(x.lt.50.0) then
+        call random_number(r2)
+        phi2=twopi*r2
+        z=a*cmplx(cos(phi2),sin(phi2))
+     endif
+     cspread(NFFT-i)=z
+  enddo
+
+  izh=fspread/df
+  do i=-izh,izh
+     f=i*df
+     j=i
+     if(j.lt.0) j=j+nfft
+     s=real(cspread(j))**2 + aimag(cspread(j))**2
+!          write(23,3000) f,s,cspread(j)
+!3000      format(f10.3,3f12.6)
+  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
+  c=cspread*c                           !Apply Rayleigh fading to c()
+
+  do i=0,NFFT-1
+     p=real(cspread(i))**2 + aimag(cspread(i))**2
+!     write(24,3010) i,p,cspread(i)
+!3010 format(i8,3f12.6)
+  enddo
+
+  return
+end subroutine dopspread
diff --git a/lib/fsk4hf/genfsk4.f90 b/lib/fsk4hf/genfsk4.f90
new file mode 100644
index 000000000..4d6805cc1
--- /dev/null
+++ b/lib/fsk4hf/genfsk4.f90
@@ -0,0 +1,63 @@
+subroutine genfsk4(id,f0,c)
+
+  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 (NFFT=512*1024)
+  parameter (NSYNC=NS*NSPS)
+  parameter (NDOWN=168)
+  parameter (NFFT2=NZ/NDOWN,NH2=NFFT2/2) !3200
+  parameter (NSPSD=NFFT2/NN)
+
+  complex c(0:NFFT-1)                   !Complex waveform
+  complex cf(0:NFFT-1)
+  real*8 twopi,dt,fs,baud,f0,dphi,phi
+  integer id(NN)                        !Encoded 2-bit data (values 0-3)
+  integer icos4(4)                      !4x4 Costas array
+  data icos4/0,1,3,2/
+
+  twopi=8.d0*atan(1.d0)
+  fs=12000.d0
+  dt=1.0/fs
+  baud=1.d0/(NSPS*dt)
+
+! Generate the 4-FSK waveform
+  x=0.
+  c=0.
+  phi=0.d0
+  k=-1
+  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
+
+  nh=NFFT/2
+  df=12000.0/NFFT
+  cf=c
+  call four2a(cf,NFFT,1,-1,1)           !Transform to frequency domain
+
+  if(sum(id).ne.0) then
+     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.
+     enddo
+  endif
+
+  c=cf
+  call four2a(c,NFFT,1,1,1)            !Transform back to time domain
+  c=c/nfft
+  
+  return
+end subroutine genfsk4