Confusion about how to use LPC, and the Numerical Error it gives

[Dimev]

Describe the bug
I don't understand what the Numerical error, input ill-conditioned?" error means in the LPC function, and I also find the documentation a bit confusing, as it doesn't really specify how I can use the LPC function.

To Reproduce

import numpy as np
import librosa

y = np.sin(np.arange(100) * 0.1) * 0.1

y_hat = librosa.lpc(y, 16) # throws said error here

Software versions*

>>> import platform; print(platform.platform())
Windows-10-10.0.19041-SP0
>>> import sys; print("Python", sys.version)
Python 3.7.4 (tags/v3.7.4:e09359112e, Jul  8 2019, 20:34:20) [MSC v.1916 64 bit (AMD64)]
>>> import numpy; print("NumPy", numpy.__version__)
NumPy 1.17.0
>>> import scipy; print("SciPy", scipy.__version__)
SciPy 1.4.1
>>> import librosa; print("librosa", librosa.__version__)
librosa 0.7.1
>>> 
>>> librosa.show_versions()
INSTALLED VERSIONS
------------------
python: 3.7.4 (tags/v3.7.4:e09359112e, Jul  8 2019, 20:34:20) [MSC v.1916 64 bit (AMD64)]

librosa: 0.7.1

audioread: 2.1.8
numpy: 1.17.0
scipy: 1.4.1
sklearn: 0.22
joblib: 0.14.1
decorator: 4.4.1
six: 1.12.0
soundfile: 0.10.3
resampy: 0.2.2
numba: 0.46.0

numpydoc: None
sphinx: None
sphinx_rtd_theme: None
sphinxcontrib.versioning: None
sphinx-gallery: None
pytest: None
pytest-mpl: None
pytest-cov: None
matplotlib: 3.1.2
presets: None

Additional context
I'm trying to use LPC to resynthesize speech from a given audio sample, however I struggle to understand how I should do this and how I have to use LPC in general.
From what I understand from here: https://www.fon.hum.uva.nl/rob/VocalTractExamples/
It's possible to use LPC to get the filter coefficients for source-filter synthesis, which is what I'm tying to do.

[bmcfee]

I don't have a ton of experience with LPC, but a few things come to mind:

1. Your test signal might be too short for the filter length you're using. I played with your example, and even going up to N=113 seems to fix the issue. Reducing the number of coefficients also seems to fix it, though I guess that's less relevant for your use case.
2. Your test signal might be too "simple". The check that triggers the error message in question is based on the combination of forward and backward prediction errors:

   librosa/librosa/core/audio.py

   Line 958 in a53fa56

   den = q * den - bwd_pred_error[-1] ** 2 - fwd_pred_error[0] ** 2

   if these go to zero, I could see there being trouble in the check. Just to test this out, I tried your code with random gaussian noise added to the signal, and it worked fine.

I suspect these issues shouldn't occur when analyzing real data, but I'd like to have a better understanding of what's going on here. Maybe @ajweiss can chime in?

[ajweiss]

I'll take a look at the numerical problem now, I'm guessing it may have to do with scaling. Yes, LPC is what you want, it will give you the filter coefficients for an IIR filter of the given order, given the input signal and an order.

[ajweiss]

Incidentally, it's long been on my todo list to put together an example notebook that demos source filter separation and synthesis by analysis using LPC, maybe it's time to pick that back up...

[Dimev]

Ah thanks for the explanation. Some examples on this would be great!

[Dimev]

I'm currently trying to port the algorithm to rust, with the source from here: http://www.emptyloop.com/technotes/A%20tutorial%20on%20Burg's%20method,%20algorithm%20and%20recursion.pdf
It seems to output different values however compared to the one in librosa. Is this because the method is different or because I'm doing something wrong?

[Dimev]

@ajweiss Do you have any sources about LPC I can read about? I understand it can be used to get the parameters of a linear model to predict the next value of a signal, but then why do the parameters also work as a filter, that also generates roughly the same signal if you put something through it?

[ajweiss]

Hey, sorry I haven't had much chance to look at this. To quickly answer questions:

1. LPC is used for both. It's origins (as I understand) come from compression, it was actually used in the popular American speak n' spell toy from the 1980s if I recall to compress all the word vocalizations down onto a cheap enough ROM to put in an affordable toy. The general idea there is to generate filter coefficients on for frames of input data (like 30ms or so) and then store those coefficients and the residual from the model (or an approximation of it) for compression and then decompress by using the residual and the sequence of filter coefficients to reconstruct an approximation of the original signal. This is what I'd like to create a notebook for.

2. Yes, computing LP coefficients is that same as finding the coefficients for an autoregressive model/filter, I think the difference is which side of the transfer fraction you put the coefficients. (numerator is predictive model, denominator is filter, I believe). There are a few methods for finding these coefficients, yule-walker and burg are common. There's also a fast method that was published by the folks who publish the open source Opus codec (I can find this for you if you like)

3. I believe there may be a nice chapter on LPC in "Speech Synthethesis and Recognition" by John and Wendy Holmes which is a wonderful little (if handwavey) introduction to all sorts of speech DSP methods.

4. I believe it is also documented in a popular DSP textbook by John Proakis, but it has been a little while.

5. If you'd like extra eyes on your implementation, I can try to help you, but this is probably not the right place for it. Feel free to e-mail directly at adam (at) signal11 (dot) com. I can't make huge time commitments, but I'm happy to help however I can.

6. Dan Ellis (who's research lab produced this library) has some slides on LPC compression on his website at Columbia for past courses he has taught. These resources are excellent and I highly recommend perusing them if you're interested. They are easily found via Google.

7. Regarding the numerical error, that came directly from the implementation in the Marple paper referenced in the code. Essentially the denominator being positive is an invariant. The derivation of the LPC algorithm involves an autocorrelation matrix that can go singular I believe which can cause the algorithm to fail, but I have not worked this derivation to see how exactly this happens. Empirically, reducing the model order or increasing the sample length can avoid the issue. I don't have a better answer for you now, but would be interested in digging into one.

Anyway, this is just a quick note and it's from memory on stuff I haven't looked at in a little while. Please feel free to e-mail me if you have additional questions that are not librosa related. Otherwise happy to continue to conversation here for librosa related stuff (specifically ensuring the LPC feature is usable, what a nice demo notebook may look like or any other problems suggestions for librosa).

Sound good?

[ajweiss]

I should also add (in case you're asking about general applications), LPC is also useful for "analysis by synthesis" which is often used to analyze speech. I think your original reference talks about this, where the vocal tract is modeled by a source-filter system, where the source is a model of "glottal pulses" that provides excitation to an IIR filter which models the resonance of the vocal tract. (In the compression example, you store the frequency of the glottal pulse excitations and then the parameters for the AR filter which shapes it's input to model the resonant modes of the vocal tract). An additional form of analysis is to look at the frequency response of this sequence of digital filters, where peaks in the spectrogram will show these resonant modes, or formants. This along with the vocoders seems like fun for the demo. :)

[Dimev]

Ah thanks again for the detailed response!
I'll look into the sources you cited, and probably also send an email about my implementation

As for a good demo, I think it would be nice to both give an example of a predictive model and a filter, and an explanation of how/why that works. (along with the other things you mentioned)

[bmcfee]

Chiming in here -- I'm extending the API to support LPC on multichannel inputs over in #1351 (prototype code buried in the comment thread), and ran into the stability issue when applying LPC to a framed signal. In the current implementation, we fail if den <= 0, and in the multichannel case, this extends to if np.any(den <= 0) (which causes the entire function to fail). This would be bad news if we wanted to implement something like LPCC, where failure in one frame means failure for the entire computation.

Since I'm still a relative newbie at LPC business, I thought I'd ask for opinions on modifying this line:

librosa/librosa/core/audio.py

Line 943 in d9b72af

reflect_coeff = dtype(-2) * np.dot(bwd_pred_error, fwd_pred_error) / dtype(den)

to something like

        reflect_coeff = dtype(-2) * np.dot(bwd_pred_error, fwd_pred_error) / (dtype(den) + epsilon)

where epsilon = librosa.util.tiny(dtype), and dropping the underflow condition in the loop. (Alternatively, we could do np.max(den, epsilon).) This way we'll avoid divide by zeros, and potentially produce some garbage features, but not cancel the entire computation.

Thoughts? (Esp. @ajweiss ?)

[ajweiss]

I think the concern here, at least theoretically, is that a failed computation could produce a filter with poles outside the unit circle. In this case, the filter becomes unstable and can cease responding to changes in input. This is theoretical though, and I'm not sure the actual impacts on a real digital filter. I'm just getting my day started, I'll take a look in depth later and get back to you!

…

On Thu, Aug 5, 2021, 7:39 AM Brian McFee ***@***.***> wrote: Chiming in here -- I'm extending the API to support LPC on multichannel inputs over in #1351 <#1351> (prototype code buried in the comment thread), and ran into the stability issue when applying LPC to a framed signal. In the current implementation, we fail if den <= 0, and in the multichannel case, this extends to if np.any(den <= 0) (which causes the entire function to fail). This would be bad news if we wanted to implement something like LPCC, where failure in one frame means failure for the entire computation. Since I'm still a relative newbie at LPC business, I thought I'd ask for opinions on modifying this line: https://github.com/librosa/librosa/blob/d9b72af964820fc796425c55d3dccd7f6b908bf9/librosa/core/audio.py#L943 to something like reflect_coeff = dtype(-2) * np.dot(bwd_pred_error, fwd_pred_error) / (dtype(den) + epsilon) where epsilon = librosa.util.tiny(dtype), and dropping the underflow condition in the loop. (Alternatively, we could do np.max(den, epsilon).) This way we'll avoid divide by zeros, and potentially produce some garbage features, but not cancel the entire computation. Thoughts? (Esp. @ajweiss <https://github.com/ajweiss> ?) — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1277 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AADJB232NU7WGVFNH6PCM4LT3KPDXANCNFSM4W5OWX7A> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email> .

[bmcfee]

I think the concern here, at least theoretically, is that a failed computation could produce a filter with poles outside the unit circle. In this case, the filter becomes unstable and can cease responding to changes in input. This is theoretical though, and I'm not sure the actual impacts on a real digital filter.

Sure. I think adding a floor to the divisor still achieves that goal, but we'd suffer a little in estimation error. I haven't thought this through deeply though, so perhaps I'm missing something.

[bmcfee]

Just a quick heads up that I tried out this idea, and it seems to work fine. An example that previously caused a numerical underflow seems relatively well behaved; the coefficients don't blow up on us (capping out around ~200 in this case). So, not a definitive proof that things will always come out okay, but it seems like a promising start.

[ajweiss]

Interesting. Does the spectrum of the input signal and the frequency response of the found filter look similar? (example of how to get the frequency response using freqz here on this old branch https://github.com/ajweiss/librosa/blob/lpcspec/librosa/core/spectrum.py#L1580

[ajweiss]

Looking at the code now, I think it's probably fine. If we wanted to be pedantic then maybe instead of throwing the FloatingPointError we could instead emit a warning (maybe a flag for the warning so it doesn't emit one for every loop iteration) and then add the eps in?

By doing this we do deviate from the paper definition some more, but it's not the first time someone just added an eps in to keep the machine running. :)

[bmcfee]

If we wanted to be pedantic then maybe instead of throwing the FloatingPointError we could instead emit a warning (maybe a flag for the warning so it doesn't emit one for every loop iteration) and then add the eps in?

Yes, flag for the warning is a good idea. This has the added benefit of removing a branch point from the inner loop, so it should be substantially faster as well.

I'll hold off on implementing that until after merging #1351 .

[miccio-dk]

Hey everyone!
I see that #1351 has been merged (hurray!) so I was wondering if any help was needed with getting the numerical error sorted out. I'd be happy to submit a PR. I'd also be very happy to share some code on how to use LPC for analysis and synthesis (@ajweiss).

PS: funnily enough, I get the FloatingPointError only if I low-pass my input signal 😅

[bmcfee]

so I was wondering if any help was needed with getting the numerical error sorted out. I'd be happy to submit a PR.

I'm actually not sure anything needs to be done - I got impatient because the exception was breaking some test cases that really ought to have worked.

I had commented out the exception in #1351, see here:

librosa/librosa/core/audio.py

Lines 944 to 947 in eb603e7

	for i in range(order):
	# can be removed if we keep the epsilon bias
	# if np.any(den <= 0):
	# raise FloatingPointError("numerical error, input ill-conditioned?")

and the epsilon floor is implemented below:

librosa/librosa/core/audio.py

Line 954 in eb603e7

reflect_coeff[0] /= den[0] + epsilon

I'd also be very happy to share some code on how to use LPC for analysis and synthesis

That would be great! Maybe as an example to add to the gallery?

[miccio-dk]

Oh, I see now! I didn't look into the source so I thought it hadn't been addressed yet.
Then I shall wait for v0.82 to come out (meanwhile a lower-order LP filter what all I needed to not trigger the error) :)

That would be great! Maybe as an example to add to the gallery?

Do you mean the examples here or here? I suppose it's the latter. I will find some time to add a PR.
Cheers!

[bmcfee]

Then I shall wait for v0.82 to come out

It'll be 0.9, but hopefully before end-of-year!

I suppose it's the latter.

Correct - my plan is to deprecate the examples folder #1328 and shift entirely to notebook-based examples built in the documentation.

I will find some time to add a PR.

🎉

[ajweiss]

In terms of analysis, I have a branch (https://github.com/ajweiss/librosa/tree/lpcspec) that implements both a short time LPC transform (windowed LPC coeffs) and also a function that then takes those coefficients and computes their spectrum. It hasn't been touched in quite a while, but is a pretty good start (I think) on the analysis side of things. What were you thinking in terms of analysis? I've been thinking about dusting this off for quite some time.

I'd also be interested in collaborating on examples. In general I was considering LPC spectral analysis and perhaps a basic compression demo... What did you have in mind?

[bmcfee]

In terms of analysis, I have a branch (https://github.com/ajweiss/librosa/tree/lpcspec) that implements both a short time LPC transform (windowed LPC coeffs) and also a function that then takes those coefficients and computes their spectrum. It hasn't been touched in quite a while, but is a pretty good start (I think) on the analysis side of things.

That sounds great! I imagine it's much easier to implement now that the core LPC function can directly operate along a target axis (eg after framing).

[miccio-dk]

THanks for the pointer, @ajweiss! I will take a look at your code.

[ajweiss]

Okay @miccio-dk keep me posted on the LPC spectrum and STLPC stuff specifically. If there isn't overlap there with the stuff you're hoping to do (or the stuff that is done is sufficiently complete that you don't feel like it would make for an interesting project), let me know and I'll pick it back up and create a PR from it. Otherwise, I'll get out of your way. :)

[bmcfee]

Is there anything left to do here? It seems to me that the original issue is now resolved (since 0.9.0). If we want to discuss LPCC, perhaps that should be spun into a new feature request issue.

[bmcfee]

I'll close this one out as resolved, but I'd still be keen to discuss extending our LPC functionality in other issues.