Fri Mar 8 11:12:36 EST 2013

Exponentials: exp(spline(t))

Goal: compute the pole path of a 2nd order IIR filter as an
exponential of a complex spline curve, approximating a parameter curve
which is itself exponential or sinusoidal.

Simplification: exact curve tracking is not so important (say 1-2% is
ok), but instabilities should not accumulate.

Two subproblems:

- Kernel routine: exp ( a0 + a1 t + a2 t^2 ) for ai complex and t
  typically 0->63.

- Spline parameter update at block boundaries, compensating for any
  inaccuracies of the kernel routine.

a1 and a2 are probably << 1, so some approximation might be used.


* amplitude accuracy is probably more important than phase.  Use
  followed by normalization?

* use linear or multiplicative path.  problem with linear path == more

* are there any useful bounds on the delta over 64 samples?  e.g. more
  than pi/2 is in 64 samples is probably exceptional.  it seems delta
  is going to be small enough to allow for approximation of
  exponentials occuring in endpoint/increment operations.  -> what is
  the most natural thing to do here?  ignoring numerical errors, it's
  probably possible to get exact figures too.

Practically: limit to linear spline

This makes the update equations for the complex fiter coefficient
curve something like this (for maginitude 1, c=cos(dp), s=sin(dp)

   x+ = cx + sy
   y+ = cy - sx

For numerical stability, best ordered as increments:

   dx = (c-1)x + sy
   dy = (c-1)y - sy

here, (c + i s) comes from  (c' + i s')^{1/64}


- when c and s are taken from approximations, what does (c + i s)^64
  look like, ideally and numerically?

- can we make sure that errors do not introduce instability,
  i.e. poles should stay inside the unit circle, or, filters should be

- should phasor and amplitude be computed separately?

- subproblem: given current phase/amplitude, compute linear curve to

Subproblem: linear curve to setpoint

- In: phasor p(0), desired phase, log amplitude

- Out: p(n) = p(0)q^n

Given the polynomial used to compute sin/cos-1/exp, this is just a
nonlinear equation, which can probably be solved using newton-raphson.

- construct equation
- construct derivative expression (AI)
- construct NR expression

Find q   s.t.   p q ^n = exp (P') = p'

Here p' is computed using some exp() approximation that is good enough
for "humanly" accurate parameter mapping.  Trouble is the presence of
q, otherwise the ^n could be absorbed by the exp.

What about expressing this differently:

q = (exp(P') / p) ^ {1/n}

q is very close to 1, so it might be better to express

q'(P',p) = (exp(P') / p) ^ {1/n} - 1

To simplify, forget about the correction that takes into account
current p vs desired p (previous setpoint): apply it later.  The focus
should really be on the evaluation of the exponential.

      q01(P1) = exp((P1-P0) / n)

Given q01, we have p1 = p0 q01^n, which can be compared to the desired
value exp(P1).  If we don't make any correction here, eventually we'll
get off track: any systematic error in q_n,n+1 is going to be

Correction: approximate (exp(P1)/p1)^{1/n}, and use it to scale q01.
Up to approximation, this will cancel *all* the errors accumulated.

The approximation could use the log of p1, i.e. exp((P1-log(p1)/n)).
The weak point in all of this is log(p1), since p1 has the largest
deviation wrt. the base of any expansion = 1.

What about incrementally computing the estimate of exp(P1) by
successive squaring?  It could even be spread over time.


- Core routine is exp((P_k+1-P_k) / n), which is a small number ->
  gives fast convergence.

- There is systematic drift.  Some absolute value feedback is
  necessary to compensate.  Combine setpoint exp(P_k) with computed
  value p_k.  The exp(P_k) has larger error, but might be accurate
  enough to perform compensation.

- Stability is probably most important: approximation exp'(x) such
  that |exp'(x)| < |exp(x)|