[Clam-devel] [patch] patch for finding polynomial roots: the route to formants

[David García Garzón]

On Monday 09 July 2007 22:43:53 abe wrote:
> Hi David,
> Thanks for the comments.  My replies are below...
>
> David García Garzón wrote:
> > On Thursday 05 July 2007 19:18:33 abe kazemzadeh wrote:
> >> Hi All,
> >>
> >> Here's a patch of what I've been working on.  I added a method to
> >> LPModel to get the roots (poles) of the LPC coefficients.  Since it is a
> >> general method,
> >
> > Not sure if this is not already on Sandra's code, did you checked it?
>
> The idea of making a polynomial class I got from Sandra's code.  There
> is a root finding algorithm in Sandra's code as well, but I didn't use
> that one.  Now that I have a working prototype I can try testing both of
> ours. Also, like you say, since the algorithm I implemented uses
> matrices, I should try testing some matrix libraries.
>
> >> not specific to LPC, I made a new class, Polynomial, in
> >> CLAM/src/Standard. I'm not sure if this is the best location to put it.
> >
> > Standard? It is.
> >
> >> Also, in this class
> >> there are
> >> methods that might be better in separate classes (eg, finding the
> >> eigenvalues
> >> of a matrix:  there seem to be some stubs for this purpose in
> >> CLAM/src/Standard/MatrixTmplDec.hxx:217), so let me know if anyone has
> >> ideas about this.
> >
> > I think that the method is ok where you put it. I don't like at all the
> > Matrix object. I would like to reduce any remaining use of Matrix to drop
> > it. We could eventually adopt any good matrix library as dependency. Will
> > you mind using a plain std::vector instead a matrix?
>
> Sure I can give this a try.  Also, if I find a library that is good I
> can use that.  By the way, what is it that you don't like about the
> Matrix object?   I also don't like it, mainly b/c it was hard to figure
> out how to use (at first, I thought that I could just declare a Matrix
> object, but it turned out that I needed a MatrixTmplDef).  Is that what
> you dislike about it or are there other reasons?

I agree with you with the non obvious usage. Also that it has too much 
functionality without testing that is has not been used. And, as most of the 
legacy classes is over designed and hard to maintain.


> >> Overview:
> >> I added a method to LPModel, LPModel::ToRoots(), which just calls the
> >> root solving function of the new class, Polynomial::PolyRoots() on the
> >> lpc coefficients.
> >> PolyRoots takes the lpc coefficients, creates a companion matrix
> >> (Polynomial::BuildCompanion()), balances it,and then gets the
> >> eigenvalues of
> >> this matrix (Polynomial::EigenHessenberg()).  This output is the roots,
> >> which
> >> are the formants of speech.  However, I haven't converted them from
> >> their complex
> >> representation to the frequency yet.  Also, for actual formant tracking,
> >> I still need to add
> >> ways to smooth the output.  It's still a bit rough, but I wanted to get
> >> feedback before
> >> further work.
> >
> > I have no knowledge on the subject, but seems ok to me.
> >
> >> Some issues:
> >>
> >> -right now the default lpc order is 11.  This seems like the textbook
> >> value that
> >> people quote, but I think that it applies to speech coding in 8kHz, so
> >> for higher
> >> sampling rates it might be better to have a higher
> >> order.  I dug out the notes from
> >> the speech processing class that I took and there is a nice derivation
> >> of how to pick
> >> the order (based on the sampling rate, the length of the vocal tract,
> >> and the speed of
> >> sound).  I'm not sure how feasible it would be to get the LPModel class
> >> to configure
> >> itself based on the sampling rate, but this would make it convenient for
> >> the user.  Either
> >> that or downsampling to 8kHz before calculating the LPC/formants (it
> >> seems that that
> >> might be what Sandra Gilabert did).
> >
> > If you know how to estimate the order taking into account the sampling
> > rate i think it is the nicer way. But also the resampling can be
> > convenient. Sandra's code included a resampler but we were thinking about
> > adding resampling processing module to clam by using the libsamplerate
> > [1] library.
> >
> > [1] http://www.mega-nerd.com/SRC/index.html
>
> Cool. I'll look into this some more and see what's best.
> Down/resampling might be reasonable for speech purposes since the
> formants needed for intelligible speech (f1-f3) are fine for the range
> of 100-4000 Hz.   Having a higher LPC/higher sampling rates/more
> formants may have some benefits if higher formants are needed, but could
> have a drawback of spurious formants.  As far as performance is
> concerned, I don't have much of an idea which method would be  best.

I have received very good references about what libsamplerate is able to do in 
terms of quality and efficiency. If you decide to resample, please, use it. 
If not, at some point we will implement such processing for other purposes.

> >> -I implemented one of the algorithms that I found online that seemed
> >> good. I
> >> considered LAPACK++, but didn't use it b/c it seemed like a lot to learn
> >> when I
> >> was already learning clam.  In retrospect, it took me a fair amount of
> >> time to
> >> translate, debug, and test the algorithm I used, making LAPACK seem good
> >> in retrospect,
> >> so I wanted to see if anyone on the list is familiar with it LAPACK.
> >
> > Not me. But isn't there any FOSS library which does that?
>
> I think LAPACK++ is FOSS.  It seems like there are 2 versions of
> LAPACK++ (not to mention the fortran version LAPACK w/o the ++... sorry
> if that was unclear) so I'm not too sure about the licence.

Sorry, i though LAPACK was an algorithm to solve linear systems not a library. 
As i said i am not familiar to it.

> >> -one of the constants in the algorithm is epsilon, a tiny value such
> >> that anything
> >> less is negligible.  I was trying to see if CLAM has such a constant and
> >> I found
> >> something in CLAM/test/UnitTests/cppUnitHelper.hxx:174, but I couldn't
> >> figure out
> >> what it was or how to use it.
> >
> > Epsilon in tests is used to accept two double values as equal. Does it
> > give you any conflict?
>
> The problem I had with this is that I didn't want to use it in tests but
> rather in the code itself... I couldn't figure out a way to do this.
> For now, I'm just using an arbitrary very small number and it seems to
> work okay.

There is such a constant in the <limits> header.
http://www.tacc.utexas.edu/services/userguides/pgi/pgC++_lib/stdlibcr/num_5679.htm#Member 
Fields and Functionsepsilon()

"Returns the machine epsilon (the difference between 1 and the least value 
greater than 1 that is representable). This function is meaningful for 
floating point types only."

Not exactly what you look for. Take a look to the other functions there.


> >> -I tested out the algorithm manually, but I was wondering about
> >> automatic tests: would that be
> >> good or necessary, and if so, could someone point me in the right
> >> direction for doing this.
> >
> > An algorithm can be covered by Back to Back tests or Unit tests:
> > - If some (simple) input has some known output you can manually generate
> > them and test as in a back to back.
> > - If you have a reference implementation generate I/O data for the B2B
> > test. - If you hand checked for the results on real data, do a back to
> > back test to be warned when the results change. You should provide some
> > criteria to validate further changes.
> > - That can be automated when your algorithm has some fitness criteria. I
> > mean, if you are testing a segment extractor and you have a hand
> > annotated wave and a fitness function.
> >
> > I have to rewrite this as a wiki page.
>
> Thanks, this and the new wiki page are helpful... Also, if you can put
> pointers to some good examples in clam that would be good (where the
> tests are built into clam library/applications).

Of course, that text is mostly a cut and paste some abstract guidelines we had 
to write on what kind of tests we did in CLAM to review the descriptors. Some 
concrete instructions on what are the tests should be given in order to 
implement them. Yesterday i gave a Roman some guidelines on that. I could 
part from that point.