Tue Jan 4 11:09:54 EST 2011

Difference Calculus

Following up on the references from Tim Stilson's PhD[2] thesis.
Already have the book "Finite Difference Equations"[1] from 1961 by
Levy and Lessman.  ( I found the original in a used book store without
ISBN, [1] is the Dover reissue. )

The references Tim mentions are [3][4][5][6][7][8] (books and
paywalled papers).  I found one accessible paper on the subject[9].

For x[k] a real valued sequence |N - > |R, we define the operator z:

      z x[k] = x[k+1]

Then we can define d = z - 1 or

      d x[k] = x[k+1] - x[k]

These are normalized equations with time step equal to 1.  Naively
relating d to the s-plane requires s/f_s normalization.

Before getting into cargo cult DSP math, one needs to keep in mind that:

  * The d and z transforms of a discrete filter contain the same
    information: they are merely described in terms of different
    primitives: the unit integrator (accumulator) vs. the unit delay.

  * In the limit case of infinite sampling rate filter coefficients of
    d-plane and s-plane coincide.  However, naively converting between
    the discrete and continuous case for finite sample rates is still

    In fact, naive conversion is the forward difference transform wich
    maps s -> z - 1, and which has known stability issues for
    coefficients away from z = 1.

More general, all the caveats concerning traditional filter
discretizations still apply.  One can derive reformulated versions of
the well known transforms: forward/backward difference, bilinear,
impulse response invariance (pole mapping) and pole-zero mapping.
These then look like naive (f_s scaled) direct s -> d mappings + some
correction terms O(s / f_s).

The main use seems to be one of implementation: coefficient
sensitivity will be less in d-plane formulation for low-frequency
poles and zeros.

Of course one can start by simply substituting z = d + 1 in a fully
expanded z-plane formulation, but that is not very insightful.
Instead we do the following:

  * Factor the z-plane formulation first as a product of first order
    terms (z-z_k), including complementary complex ones.  Performing
    the substitution in this form yields terms (d + 1 - z_k).  We name
    d_k = 1 - z_k.  For low-frequency phenomena it happens that the
    z_k are all close to one meaning the d_k are small in magnitude.

    So essentially, a formulation in d-plane form eliminates a
    subtraction that is the cause of loss of precision in the z-plane

  * 1st order IIR sections in the form of b/(d-a) can be implemented
    directly using leaky integrators[11].  A state variable 2nd order
    section can be used to implement the complex conjugate 1st order

For IIR filter design (when everything is analytic) there is little
practical difference between z-plane and d-plane approach: simply
transform (analytically) when you're done.

For FIR filter design which is mostly numerical in nature, the d-plane
approach seems also useful for avoiding numerical instabilities.
I.e. to formulate the optimization problem in such a way that the FIR
filter appears in factored form, instead of multiplied out.

( From my own experience I've found that FIR filter design problems
  involving low frequencies are seriously plagued by convergence
  issues due to insufficient numerical precision.  I recall that using
  a spectral factorization method to design an FFT window filter I
  would obtain results that where clearly wrong after visual
  inspection, i.e. one coefficient value "sticking out" while there
  should be overall symmetry.  It would surprise me if this hasn't
  been realized long before.  Look into ladder FIR filters.  )

[1] isbn://0486672603
[2] https://ccrma.stanford.edu/~stilti/papers/Welcome.html
[3] http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=1104162
[4] http://md1.csa.com/partners/viewrecord.php?requester=gs&collection=TRD&recid=2270973CI&q=&uid=790228084&setcookie=yes
[5] isbn://9780817639341
[6] http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=0123294
[7] http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=73564
[8] http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=219661
[9] http://sites.google.com/site/jamesgibsonhomepage/projects/DeltaOperatorCaseStudy.pdf
[11] entry://20110104-124734