# Last edited on 2012-09-25 22:42:09 by stolfilocal

GOAL

  Develop methods for generating orthonormal bases for discrete signals,
  other than the Hartley and Fourier basis.  Specifically we want a basis
  that will capture the notion of "detail scale" or "frequency" but 
  can be used for non-periodic signals without the need for a 
  windowing function. 

METHODS

  For this note, a discrete signal is a vector of {R^n} whose coordinates are assumed
  to be equally spaced samples of a signal {z} along a time or space axis {x}.
  Unlike in Fourier theory, we do not assume that the signal is periodic;
  that is, sample {z[n-1]} is not followed by {z[0]}.
  
  An orthonormal basis for such space is a matrix {phi[i,j]} whose rows 
  have unit norm and are pairwise orthogonal under the inner product
  
    { <y,z> = (1/n)*SUM{ y[j]*z[j] : j IN {0..n-1} } }
  
  In this note we assume that the row index {i} of the matrix {phi[i,j]}
  is cyclic, i.e. {phi[i+n] == phi[i]} for any integer {i}. We do NOT
  assume that the column index {j} is cyclic.

  We consider two methods for generating an orthogonal basis {phi[i,j]} for such signals:
  
    (1) The two-term recurrence of orthogonal polynomials, 
        modified for {n}-vectors.
        
    (2) Rotating the first basis vector by a basis permutation matrix.
    
TWO-TERM RECURRENCE

  We use the two-term recurrence 
  
    { phi[n] = phi[n-1]*(a[n-1]*I + b[n-1]*D) +  phi[n-2]*c[n-1] }
  
  where {I} is the identity operator of {R^n}, {D} is a fixed linear operator of
  {R^n} (the /degree-raising operator/), and {a[k],b[k],c[k]} are 
  coefficients that depend on {k}.
  
  The recurrence starts with
  
    { phi[0] = (1,1,1,...,1) }
    { phi[1] = phi[0](a[0]*I + b[0]*D) }
    
  The coefficients {a[k],b[k],c[k]} are chosen so as to ensure that
  <phi[n],phi[n]> = 1, <phi[n],phi[n-1]> = 0, and <phi[n],phi[n-2]> = 0.
  Namely,

    cDD1 = < D(phi[k-1]),D(phi[k-1]) >;
    cD1 =  < D(phi[k-1]),  phi[k-1]  >;
    cD2 =  < D(phi[k-1]),  phi[k-2]  > (ou 0 se k = 1);
    h = cDD1 - cD1*cD1 - cD2*cD2;
    b[k] = 1/sqrt(h);
    a[k] = -b[k]*cD1;
    c[k] = -b[k]*cD2;
  
  In the theory of (continuous) orthogonal polynomials, the {D} operator
  is pointwise multiplication by the argument {x}. In the discrete
  theory, it can be a fixed diagonal operator, i.e, coordinatewise
  multiplication by some signal, that turns the first all-ones vector {phi[0]}
  into the second basis vector {phi[1]}.
  
  Without loss of generality we cn assume that the coordinate axes have
  been permuted so that $phi[1]$ is non-decreasing. This implies that
  {phi[1]} is the lowest-frequency (or lowest-degree) basis signal ---
  some ramp-like signal, not necessarily real.
  
  The script "orthopol.gawk" implements this method for some choices 
  of the operator {D}.  The choices "hcos", "hqua" give trig functions.
  The choice "poly" gives polynomials. 
  
  The "hqua" basis is particularly nice.  The first element {phi[0]} 
  is all ones (frequency 0).  For {i != 0}, the generic element is 
  
     { phi[i,j] = sqrt(2)*cos(pi*i*(2*n-2*j-1)/(2*n)) }
     
  Note that this element is a trig wave with frequency {i/2}. Unlike the
  Hartley and Fourier bases, the *apparent* frequency of {phi[i]}
  increases smoothly as {i} ranges from 1 to {n-1}.
  
  For {i == 1}, for instance, it gives a sinsoudal ramp that rises from
  almost {-sqrt(2)} to almost {+sqrt(2)} as {j} ranges from {0} to
  {n-1}.
  
CYCLIC RECURRENCE
    
  For every orthonormal basis {phi} of {R^n} there is a rigid rotation
  operator {M} of {R^n} that cycles over the basis elements, that is
  
    { M(phi[i]) = phi[i+1] }
  
  The matrix of the operator {M} is 
  
    { M[i,j] = (1/n)*SUM{ phi[k,i]*phi[k+1,j] : k IN {0..n-1} } }
    
  Note that the summation is NOT the inner product of two rows of {phi}.
  
  For the canonical orthonormal basis of {R^n}, with {phi[i,j] =
  sqrt(n)*(i == j)}, the matrix {M} is the cyclic forward shift matrix
  
    { M[i,j] = ((i+1)%n == j) }
    
  For the (complex) Fourier basis, {M} is the (complex) diagonal matrix
  
    { M[i,j] = (w^j)*(i == j) }
    
  where {w} is the first non-trivial {n}th root of unity,
  {w = exp(2*pi*im/n)} where {im = sqrt(-1)}.
  
  For the Hartley basis, 
  
    { M[i,j] = (i == j ? cos(2*pi*j/n) : 0) + ((i+j)%n == 0 ? sin(2*pi*j/n) : 0) }
    
  Not clear yet what it is for the "hqua" basis defined above.
  
INTERESTING QUESTIONS FOR DISCRETE ORTHONORMAL SIGNAL BASES

  For some operators {D} the recurrence fails before generating {n}
  basis elements, because the equations for {a[k],b[k],c[k]} cannot be solved.
  What are the conditions on {D} to avoid this problem?

  How do we get from the {D} frequency-raising operator to the rotation matrix {M}?

  The FFT algorithm is based on a factorization of the {phi} matrix into 
  the product of {log(n)} sparse matrices, each with {~2n} nonzero entries.
  Which orthonormal bases of {R^n} have a similar decomposition?

  If the rotation matrix {M} is sparse, can we find such a sparse 
  factorization of {phi}?
  
INTERESTING QUESTIONS FOR CONTINUOUS ORTHONORMAL POLYNOMIALS

  What is the rotation operator {M} for the standard families of 
  continuous orthonormal polynomials?
  
  What do we get if the use {D} operators that are not just multiplication
  by {x}; say, indefinite integration from the midpoint of the domain?