1subroutine twkfreq(c3,c4,npts,fsample,a)
2
3  complex c3(npts)
4  complex c4(npts)
5  complex w,wstep
6  real a(3)
7  data twopi/6.283185307/
8
9! Mix the complex signal
10  w=1.0
11  wstep=1.0
12  x0=0.5*(npts+1)
13  s=2.0/npts
14  do i=1,npts
15     x=s*(i-x0)
16     p2=1.5*x*x - 0.5
17!          p3=2.5*(x**3) - 1.5*x
18!          p4=4.375*(x**4) - 3.75*(x**2) + 0.375
19     dphi=(a(1) + x*a(2) + p2*a(3)) * (twopi/fsample)
20     wstep=cmplx(cos(dphi),sin(dphi))
21     w=w*wstep
22     c4(i)=w*c3(i)
23  enddo
24
25  return
26end subroutine twkfreq
27