Are there any analytic approaches to the design of quadratic polynomial filters? The following paper describes an application of the Teager-Kaiser energy operator to x-ray image enhancement:


  Reinhard Bernstein, Michael S. Moore and Sanjit K. Mitra, "Adjustable
  Quadratic Filters for Image Enhancement" Proc. IEEE International
  Conference on Image Processing (ICIP), Santa Barbara, CA, vol. 1, pp.
  287-290, Oct. 1997. http://vision.ece.ucsb.edu/publications/view_abstract.cgi?52


The authors develop intuition for the behavior of the filter through analogy with a similar linear operator (i.e. "Thus the output of a Teager filter is approximately equal to a highpass filter response weighted by the local mean."). For the sake of precision, by quadratic polynomial filters, I mean non-linear, non-recursive filters that can be completely characterized by a truncated Volterra Series, as follows (for the 1D case):

$$
y(n) = \sum_{m_1=0}^{N_1-1}{ h_1(m_1)x(n - m_1) } + \sum_{m_1=0}^{N_2-1}{\hphantom{.}\sum_{m_2=0}^{N_2-1}{ h_2(m_1,m_2)x(n - m_1)x(n - m_2) } } 
$$

It seems that most approaches to the design of low-order polynomial filters involve system identification frameworks, but without any deep understanding as to why the estimated filters work.  Are analytic approaches based on linear analogies currently the state-of-the-art, or are there any known mathematical tools that can be used?
Not really an answer (hence this is community wiki), but I thought we should capture @Mohammad's code:

%Mohammad Z

%Teager-Kaiser Non-Linear Energy Operator. 
function [out] = TKEO(x)
    N = length(x);
    x = x(:).';
    temp = x(2:N-1).^2 - x(3:N).*x(1:N-2);
    out = [temp(1) temp temp(end)];    
end