Legendere / "Optimum-L" filter designs SOLVED!

[Music Engineer]

http://www.crbond.com/papers/lopt.pdf

ah. finally some info available. great. the CRBond seems to be a nice guy. OK - so constructing the legendre polynomials themselves is straightforward by applying recurrence relations. and then there's this definite integral involving a weighted sum of such legendre polynomials. can't this be solved analytically? i have seen some hints in your code on numeric integration. can you explain how this works? after all, the result should be our set of polynomial coefficients, whereas i only know about numeric-integration algorithms that result in a value (the approximated value of the definite integral)

[Music Engineer]

or wait - no - it says:

Code: Select all

 
  //
  //  calculate definite integral
  //
  for (i=1;i<=n;i++) {
    if (i > 1) {
      c0 = -m_s[0];
      for (j=1;j<i+1;j++) {
        c1 = -m_s[j] + 2.0*m_s[j-1];
        m_s[j-1] = c0;
        c0 = c1;
      }
      c1 = 2.0*m_s[i];
      m_s[i] = c0;
      m_s[i+1] = c1;
    }
    for (j=i;j>0;j--) {
      m_w[j] += (m_v[i]*m_s[j]);
    }
  }

is this how the integral is evaluated analytically, resulting in a set of polynomial coefficients?

[thevinn]

is this how the integral is evaluated analytically, resulting in a set of polynomial coefficients?

I'm going to go with yes, and keep in mind I didn't write that part of the code, it was provided by C. R. Bond.

If there's an analytical solution, I'm all ears! I prefer analytical to numerical in most cases.

[Music Engineer]

I'm going to go with yes, and keep in mind I didn't write that part of the code, it was provided by C. R. Bond.

very nice of him.

If there's an analytical solution, I'm all ears! I prefer analytical to numerical in most cases.

so do i. i tend to think that this bit of code indeed implements the analytical solution. because numerical integration algorithms do generally just result in some number whereas this algorithm apparently computes the polynomial coefficients.

[lubomir.ivanov]

so do i. i tend to think that this bit of code indeed implements the analytical solution. because numerical integration algorithms do generally just result in some number whereas this algorithm apparently computes the polynomial coefficients.

robin,
yes this is an analytical solution, which comes from Papoulis originally i believe.

i haven't looked at the code, but in general the legendre ortogonal relation (integral) has the property to be simplified into a sum of integrals which should provide an analytical solution.

such as:
http://integrals.wolfram.com/index.jsp? ... ndom=false

also alternatively a closed form (a sloane's sequence) can be expressed as:
http://mathworld.wolfram.com/images/equ ... ation2.gif
(from: http://mathworld.wolfram.com/LegendrePolynomial.html)


lubomir

edit: oops polynomal -> integral wording

--

[Music Engineer]

robin,
yes this is an analytical solution, which comes from Papoulis originally i believe.

thanks for clarifying this.

i haven't looked at the code, but in general the legendre ortogonal relation (integral) has the property to be simplified into a sum of integrals which should provide an analytical solution.

i've been thinking about this a bit and to me, it seems like an analytical solution should be possible regardless of any orthogonality or whatever condition. basically, what we have to solve is something of the form (pseudo-latex):

int_{-1}^{2u-1} p(x) dx

where u = w^2, p(x) is this weighted sum of legendre polynomials squared (which is itself a polynomial, the coefficients of which we can compute). i assume, that we represent our polynomial as an array a[0],...,a[N] of coefficients such that p(x) = a[0] x^0 + ...+ a[N]x^N. now, to find the antiderivative of p(x), we would just shift the whole array one position to the right (resulting in an array of length N+1) and multiply each coeff by its new index and fill up the new a[0] with our integration constant (probably, we just choose 0). Lets denote the antiderivative of p(x) as P(x). to perform the definite integration, we would now need to plug in upper and lower integration limits and do: int = P(2u-1) - P(-1). P(-1) will just be a number, P(2u-1) will be another polynomial in u, the coefficients of which can possibly be calculated by application of the binomial theorem. take all this with a grain of salt - i didn't really work it out thoroughly. but this is the line of reasoning along which i would try to tackle this. i might well be missing something here.

[lubomir.ivanov]

i've been thinking about this a bit and to me, it seems like an analytical solution should be possible regardless of any orthogonality or whatever condition.

yes, in this case there is a given simplicity of the polynomials and the anti-derivatives. of course, i believe it can be made "broken" in therms of possible analytical solutions quite easily.

the ease considering anti-derivatives in this case:
int_{a}^{b} x^a dx = (x^(a + 1) / (a + 1)) + c
(c quite possibly set to 0 as you mention bellow)

in the case of p3(x) there will be:
int_(5*x^3 - 3*x)/2 dx = (x^2*(-6 + 5*x^2))/8 = P3(x)

basically, what we have to solve is something of the form (pseudo-latex):

int_{-1}^{2u-1} p(x) dx

where u = w^2, p(x) is this weighted sum of legendre polynomials squared (which is itself a polynomial, the coefficients of which we can compute). i assume, that we represent our polynomial as an array a[0],...,a[N] of coefficients such that p(x) = a[0] x^0 + ...+ a[N]x^N. now, to find the antiderivative of p(x), we would just shift the whole array one position to the right (resulting in an array of length N+1) and multiply each coeff by its new index and fill up the new a[0] with our integration constant (probably, we just choose 0). Lets denote the antiderivative of p(x) as P(x). to perform the definite integration, we would now need to plug in upper and lower integration limits and do: int = P(2u-1) - P(-1). P(-1) will just be a number, P(2u-1) will be another polynomial in u, the coefficients of which can possibly be calculated by application of the binomial theorem. take all this with a grain of salt - i didn't really work it out thoroughly. but this is the line of reasoning along which i would try to tackle this. i might well be missing something here.

i confirm that overall this is pretty much spot on...and even with the code considerations in mind. not exactly sure on the binomial theorem though since there will be forms such as: (5*(2*(w^2) - 1)^3 - 3*(2*(w^2) - 1))/2
but there are alternatives...

a quick maxima experiment:

Code: Select all

/*
    example of calculations for Optimal "L" filter coefficients
    lubomir i. ivanov (neolit123 [at] gmail)
*/

/*
    Legendre polynomial of the first kind from Rodrigues representation
    (with expansion)
*/

p_legendre(n) :=
(
    ratsimp(1/(2^n)*sum
    (
        ((-1)^k * (2*n - 2*k)!) /
        (k! * (n - k)! * (n - 2*k)!) *
        x^(n - 2*k),
        k,
        0,
        floor(n/2)
    ))
);

/*
    Papoulis optimal
    1. A.Papoulis,OnMonotonicResponseFilters,Proc.IRE,47,Feb.1959,332-333
    2. http://www.crbond.com/papers/lopt.pdf
*/

/*
    odd order
*/
p_optimal_odd(n, u) := 
(
    k: (n - 1)/2,
    ratsimp(integrate
    (
        sum
        (
            /*
                ai * p_legendre(i)
            */
            ((2*i + 1) / (sqrt(2)*(k + 1))) * p_legendre(i),
            i,
            0,
            k
        )^2,
        x,
        -1,
        2*u - 1
    ))
);

/*
    even order
*/
p_optimal_even(n, u) := 
(
    k: (n - 2)/2,
    if mod(k, 2) = 0
    then mi: mod(i + 1, 2)
    else mi: mod(i, 2),
    ratsimp(integrate
    (
        (x+1)*sum
        (
            /*
                mi * ai * p_legendre(i)
            */ 
            (mi) * ((2*i + 1) / sqrt((k + 1)*(k + 2))) * p_legendre(i),
            i,
            0,
            k
        )^2,
        x,
        -1,
        2*u - 1
    ))
);

/*
    print polynomial and TF for order n
*/

print_p(n) :=
(
    if mod(n, 2) = 0
    then
        block( Hs_denom: 1 + p_optimal_even(n, -s^2),
                Pw2: p_optimal_even(n, w^2))
    else
        block( Hs_denom: 1 + p_optimal_odd(n, -s^2),
                Pw2: p_optimal_odd(n, w^2)),
    
    print(""),
    print("=============================================================="),
    print("Order: ", n),
    print("P",n,"[x] = ", p_legendre(n)),
    print("L",n,"[",w^2,"] = ", Pw2),
    print(""),
    print("H(s)*H(-s) = 1 / (1 + L",n,"[-s^2]) = 1 / (", Hs_denom, ")")
    /* ,
    print(""),
    print("roots: ", allroots(Hs_denom)),
    print("H(s) = ", 1 / Hs_denom)
    */
)$

for p: 0 while p <= 6 do print_p(p)$

files:
http://tinyurl.com/4dl52bc

unless i have mistakes in the equations, a common factor has to be removed from Ln(w^2) before the forming of the polynomial for the transfer function denominator or the roots will be off.

EDIT 24.02.11: fixed the devisor issue, added calculations for the even order optimal.

---

[Music Engineer]

thanks for the maxima stuff, lubomir. i must digest this. so far, i used maxima myself only for much simpler stuff (solving systems of equations, obtaining some antiderivative and the like).

btw.: self-correction: i think the coeff should be divided rather than multiplied by its new index when forming the antiderivative. multiplying would occur when taking the derivative.

as for the binomial cofficient thing: i was thinking that the new (antiderivative) coeffs multiply terms of the general form (2u-1)^n but i may be missing something here in my superficial look at it.

[lubomir.ivanov]

thanks for the maxima stuff, lubomir. i must digest this. so far, i used maxima myself only for much simpler stuff (solving systems of equations, obtaining some antiderivative and the like).

just updated new files above with calculations for the even order (for completeness).

btw.: self-correction: i think the coeff should be divided rather than multiplied by its new index when forming the antiderivative. multiplying would occur when taking the derivative.

yes, thats right.

as for the binomial cofficient thing: i was thinking that the new (antiderivative) coeffs multiply terms of the general form (2u-1)^n but i may be missing something here in my superficial look at it.

ah i see, well it should work when substituting with the upper limit.

---

[Music Engineer]

ah i see, well it should work when substituting with the upper limit.

yep. and then we would *somehow* (probably messy) obtain another set of polynomial coefficients that represents P(2u-1) for any x, then we would add the constant number P(-1) to the constant coeff of that polynomial and would have the polynomial coefficients that represent that definite integral as a function of x. ...or something.

btw: does anyone know how this formula with the integral was derived in the first place? is, maybe, a link to the original papoulis derivation publically available? i couldn't find anything.