Legendere / "Optimum-L" filter designs SOLVED!

[thevinn]

I see allusions to these filters but no transfer functions. Do they exist?

[Music Engineer]

this:

www.crbond.com/papers/optf2.pdf

has tabulated the poles for orders up to 10. how they were calculated - i don't know.

[thevinn]

I fired an email off to the author of that site.

[Music Engineer]

I fired an email off to the author of that site.

cool. i'd be interested in the outcome of that too.

[thevinn]

Here was the response

Code: Select all

Information on these filters is, as you found out, very difficult to come by. I've attached a C++ file which calculates the coefficients for "L" filters of arbitrary order. The algorithm differs depending on whether the order of the polynomial is even or odd and, if it's even, whether the order is divisible by 4 or not. You'll really need to find access to suitable discussions to make sense out of it, because the algorithm requires recursive calculation of internal coefficients and the integration of a polynomial.

In practice, the L(w^2) polynomials are generated by my program. The magnitude response squared is then:

                    1
      M^2(w) = ------------                    eq. 1
               1 + Ln(w^2)    
Note: Ln is the 'n'th order of the optimal polynomial --- not the natural log!

From this, form
                                 1
      h(s^2) = H(s)*H(-s) = ------------       eq.2
                             1 + Ln(-s^2)

Then find H(s) and scale it for DC gain of unity.

Example:

L3(w^2) = 3w^6 - 3w^4 +w^2                     eq.3

                  1
M^2(w) = --------------------                  eq.4
         1 + w^2 - 3w^4 +3w^6

                     1
H(s)H(-s) =  -------------------               eq.5
             1 - s^2 - s^4 - s^6

and

                  0.577
H(s) = ------------------------------          eq.6
       s^3 + 1.31s^2 + 1.359s + 0.577

I hope this is sufficient for you to put things
together the way you want. Let me know if you
need more...

[thevinn]

I'm a little bit befuddled about exactly how to go from eq.4 to eq.6

[thevinn]

I hard-coded the coefficients for the Legendre filter of order 3 into my root finder, and then calculated the poles. It worked. Here's an image. You can clearly see that the magnitude response is monotonically decreasing as expected.

[thevinn]

Without using any math, I tried shoving the coefficients of the Legendre polynomial into various places and after a lot of comparison of the output of the root finder, with the table of pole/zeros up to order 10 from the original paper, I totally guessed a working implementation. I can't justify it at all but here's a screenshot:



If you want to play with Legendre filters, its in the latest SVN source of DspFilters but I DID NOT update the /zip files yet.

Here's the algorithm:

Given D(s) = 1 + Ln(-s^2) where Ln() is the n'th order Legendre optimal polynomial, perform the following:

- Find the set of complex roots of D(s)
- Swap the real and imaginary part of each root
- Discard roots whose real part is positive

The remaining roots are the poles of the low pass analog prototype.

[thevinn]

You can view the source code for producing the Legendre polynomial coefficients, as well as my code to extract the poles, here:

http://code.google.com/p/dspfilterscpp/ ... gendre.cpp

[Music Engineer]

hey vincent, cool. having a working implementation of papoulis filter design around is ....NEW! thanks for figuring that out. so let's see if i got that right - the function:

void lopt(double *w,int n)

writes the desired polynomial coeffs (for order n) into the array "w" and then you apply your root finder to find the poles, yes?

[thevinn]

hey vincent, cool. having a working implementation of papoulis filter design around is ....NEW!

Thanks! Finally, something real after over a year and a half. I guess if I was pursuing an equivalent college level it would have taken something on the order of that time as well.

void lopt(double *w,int n)

writes the desired polynomial coeffs (for order n) into the array "w" and then you apply your root finder to find the poles, yes?

Yes, however, keep in mind that's based on H(s)*H(-s) = 1 / P(s^2) where P(s^2) is the Legendre polynomial.

Since the polynomial is in terms of s^2, any root of the form a+bi will also appear as a-bi (conjugate) but also -a+bi and -a-bi. I think...correct me if I'm wrong because I don't really have a good grasp on all this.

Anyway, if you use the root finder you will notice two things about your roots, in comparison to the numbers given for the poles in the original document. One, that the real and imaginary parts are transposed, and two that there are roots with positive reals (-a+bi and -a-bi where a+bi is a root with negative real).

My solution as I mentioned earlier was to first swap the real and imaginary parts of the roots, then discard the roots with positive real parts, finally use whatever is left as the poles.

I'm sure you can justify this treatment mathematically I just don't know how. Something having to do with solving for H(s) in terms of H(-s) and probably the complex identity conj(-i*(a+bi))==b+ai.

NOTE: I forgot to mention in the previous post, to swap the real and imaginary parts, I edited it to include that tidbit

[lubomir.ivanov]

hello vinnie,
nice to see these filters in your c++ library.

> In response to some of the questions raised by Vinnie, I prepared
> the attached file, "lopt.pdf" which does a better job of explaining
> how the transfer function is computed.
>
> The peculiar transforms used result from theorems in network theory,
> particularly focusing on realizability considerations. That is,
> given a function of s^2 it is necessary to extract those elements
> which produce a positive real function. What appears to be a
> 'swapping' of real and imaginary elements is just a result of
> certain symmetries in the pole locations for filters of this
> type. I hope the attached paper will clarify the method.
>
> I'll be posting the same file on my website over the weekend.
>
> Thanks for your interest in these remarkable and neglected
> filters.
>
> Best regards,
>
> C. Bond <email>
>
> P.S. Thanks for the link!

i also think you might be swapping some signs in there...for example:

Code: Select all

C3(w^2) = 0.27*w^2 + 0.46*w^4 + 0.27*w^6
P3(s^2) = 1 + C3(-s^2)

thats a mistake i did previously:
P3(s^2) = 1 + 0.27*s^2 + 0.46*s^4 + 0.27*s^6
instead of:
P3(s^2) = 1 - 0.27*s^2 + 0.46*s^4 - 0.27*s^6

which results in "real/imaginary parts swapped" in the conjugate pairs:
0 + 1.442*i
0 - 1.442*i
...
instead of:
1.442 + 0*i
-1.442 + 0*i
...

the brand new paper on legendre filters by c. bond:
http://www.crbond.com/papers/lopt.pdf

a small side note: the standard hurwitz polynomial in the denominator is formed by reverse factoring the set of binomials, containing the "stable" roots.

the second link in question that explains the cohen polynomial and filter derivatives:
http://www.ncbi.nlm.nih.gov/pmc/article ... 9-0524.pdf


cheers
lubomir


----

[thevinn]

i also think you might be swapping some signs in there...

Thanks for pointing this out. In fact, the bug is that I wasn't swapping signs, I just used the polynomial from H(w^2) directly. Some of that new paper from C.R.Bond sunk in and I finally understood w^2=-s^2. This essentially means reversing the sign of every other term starting with s^2 (s^2, s^6, etc).

I changed my code and then I was able to take out the hack of manually swapping the real and imaginary parts.

[thevinn]

http://www.ncbi.nlm.nih.gov/pmc/article ... 9-0524.pdf

Stupid question, is it possible to perform two dimensional convolution using a kernal exhibiting radial symmetry, by application of one dimensional IIR filters with suitable coefficients, first on one axis and then on the other axis (of the intermediate filtered data)?

I've never heard of IIR filters being used on image data...usually its done using FIRs (i.e. convolution kernels).

[lubomir.ivanov]

http://www.ncbi.nlm.nih.gov/pmc/article ... 9-0524.pdf

Stupid question, is it possible to perform two dimensional convolution using a kernal exhibiting radial symmetry, by application of one dimensional IIR filters with suitable coefficients, first on one axis and then on the other axis (of the intermediate filtered data)?

I've never heard of IIR filters being used on image data...usually its done using FIRs (i.e. convolution kernels).

i believe for these applications, firs are mostly used due to the potential of linear phase in the passband and greater simplicity. but i'm sure there are iir implementations in both software and hardware as well.

but yes...thats pretty much the theory behind a separable 2-d filtering - you can apply a 1-d filter over each axis.
of course, often times with mathematical simplicity come problems, as the quality of the 1-d filter would be of importance here. for example the passband ripple (from factor epsilon) of a chebyshev type-1 would be multiplied as (2*ripple).

on the topic of designing a 2-d filter, outer product expansion:
https://ccrma.stanford.edu/~jos/Variabl ... nsion.html
has to say, that a 2-d filter transfer function H2 can be expressed by two 1-d transfer functions H0, H1.

also consider that a 2-d filter might have problems which can only be solved if the filter is separated into two pairs of 1-d filters.
here is an example of outer product expansion with a 2-d notch iir:
http://www.ece.uvic.ca/~wslu/Publicatio ... al/J22.pdf


lubomir

---