From octave-sources-request at bevo dot che dot wisc dot edu Thu Aug 30 15:27:42 2001
Subject: Re: "expsinesweep" -- generates swept sine signal 
From: "J. Milgram" <milgramjh at nswccd dot navy dot mil>
To: octave-sources at bevo dot che dot wisc dot edu
Date: Thu, 30 Aug 2001 16:30:59 -0400



OK, I see it's up on the archive already, and no one will be able to
see the source without the hassle of de-mimifying it. So here it is in
plaintext. I've put the  writeup at

http://www.tux.org/~milgram/expsinesweep

- Judah


#
# expsinesweep - exponential sine sweep.
#
# 8/23/00, Judah Milgram, Carderock Division, Naval Surface
# Warfare Center, milgramjh at nswccd dot navy dot mil dot  (C) 2000 J.M.
# "Permission to use, copy, modify, and distribute this software
# for any purpose and without fee is hereby granted, provided that
# this notice appear in all copies and relevant documentation.
# This is free software and comes with no warranty."



function b=l2norm(v)
#
# Input: v, a vector
#
   b=0.;
   for x=v'   # transpose required for octave
       b=b+x^2;
   end
   b=sqrt(b);
end

function [ tsw, sw ] = expsinesweep(dt,frq0,frq1,alpha)
#
#  expsinesweep(dt,frq0,frq1,alpha) - exponential time sweep (chirp)
#
#
#
# Inputs:
#   	dt = time step (e.g. s)
# 	frq0 = starting frequency (Hz -- not 1/s)
#	frq1 = ending frequency   (Hz)
#      Note: If frq1 < frq0, you get a reverse sweep. No need to adjust alpha.
#	alpha = sweep speed parameter:
#	        forward sweep: n'th period = alpha * (n-1)'th period
#              reverse sweep: n'th period = alpha * (n+1)'th period
#              Try  something like  alpha = .95
#              alpha always  < 1.0
#
# Returns:
#      tsw = vector of times
#      sw = vector of sweep signal
# 

    twopi = 2*4*atan(1);

    if frq0 < 0
       error("expsinesweep: initial frequency must be > 0");
    elseif frq1 < 0
       error("expsinesweep: final frequency must be > 0");
    elseif frq0 > 1/2./dt
       error("expsinesweep: initial frequency doesn't satisfy Nyquist")
    elseif frq1 > 1/2./dt
       error("expsinesweep: final frequency doesn't satisfy Nyquist")
    elseif frq1 == frq0
       # someday we'll jazz it up by allowing  frq0==frq1,
       # then return a dwell, duration defined somehow by alpha.
       error("expsinesweep: doesn't do dwells yet.")
    elseif alpha >= 1.0
       error("expsinesweep: alpha must be < 1")
    elseif frq1 < frq0
       temp=frq1
       frq1=frq0
       frq0=temp
       backwards = 1;
    else
       backwards = 0;
    end

    # tau0, tauf in sec.
    tau0 = 1./frq0;     # initial period, nominal
    tauf = 1./frq1;     # final period,  nominal

#  Calculate duration of sweep:
    n = ceil(log(tauf/tau0)/log(alpha));  # how many cycles
    p=0; t=0.;
    while p<=n
       t = t + tau0*alpha^p;              # total up sweep time required
       p = p+1;
    end    
    nsw = ceil(t/dt);                     # how many points  in sweep.
    tsw = dt*(0:nsw-1)';                  # vector of  times

#  Solve nonlinear system for [ k a ]
#  Interative Newton-Raphson method
# 
#  k,a will have units of 1/s, not Hz.
#
    tol=.00001;   #  you can play with  these parameters. 
    damp=1;       #   "     "       "     "       "
    x0=[ twopi*frq0; 0.];   #  initial guess
    err=1.;
    i=0; maxiter=1000;
    while abs(err) > tol
      i=i+1;
      if i > maxiter
         error("expsinesweep: failed")
      elseif i > 1
         x0=x0-inv(fp)*f/damp;
      end
      k=x0(1);
      a=x0(2);
#     I really should document this.
#     Basically, get first and last periods to match tau0 and tau1
      f = [ tau0*k*exp(a*tau0)-twopi; ...
            k*t-k*exp(-a*tauf)*(t-tauf)-twopi*exp(-a*t) ];
      fp=[ tau0*exp(a*tau0) k*tau0^2*exp(a*tau0); ...
           t-exp(-a*tauf)*(t-tauf) ...
           k*tauf*(t-tauf)*exp(-a*tauf)+twopi*t*exp(-a*t) ];
      err=l2norm(f);
    end

#  Now we know the sweep parameters, so generate the  signal:

    frq = exp(a*tsw)*k;        # vector of frequencies, 1/s
    sw =  sin(frq.*tsw);       # generate the signal
    if backwards == 1          # if this is reverse sweep, flip the signal
       temp = sw;
       p=1;
       while p<=nsw
         sw(p)=temp(nsw-p+1);
         p=p+1;
       end
    end

end


#
# flattop - modified  flat top weighting function
#
# JM 10/11/00

function w = flattop(n,d)
# n = number of points
# d = fraction of interval to be constant "flat top"
# w = vector of weights

pi=4*atan(1);

# How many points in each segment:


n2=round(d*n)         #  preliminary
n1=round((n-n2)/2)
n3=n1
n2=n-n1-n3            # final


# check to see if n1, n2, n3 are sane

#  First segment:

w1 = .5-cos((0:(n1-1))*pi/(n1-1))/2;
w3 = .5+cos((0:(n3-1))*pi/(n3-1))/2;
w2 = ones(1,n2);

w = [ w1 w2 w3 ]';

end