Autocorrelation and levinson recursion (LPC)

trogluddite · 21 Oct 2011

Hi Michael, welcome to the forum.

Now, please bear in mind that I'm not really a C++ programmer, and I have no idea what a Levinson-Durbin recursion is, but I would say that the code would next to impossible in SM code

SM code is designed for very fast linear number crunching, but it is very limited when it comes to conditional statements and loops. This is intentional, as it allows SM's clever way of creating polyphonic voices for synthesisers by using SSE. In fact there are no conventional conditional statements! Only very simple ones can be realised using combinations of maths and bitwise logic. Loops are only possible if the number of iterations is fixed - and sadly the code has several instances of loops which are controlled by other variables. The code also has nested loops that are within conditional statements, further compounding the problem.

There might be an alternative method to achieve the same thing in SM, but this looks to me like the kind of very optimised matrix algorithm that would be uselessly inefficient if done any other way (SM lacks 'real time' FFT for similar reasons). There may be more chance if you used assembly rather than SM code - but that is far from a trivial task, and without a lot more flowcharting etc., I can't say for sure that you could overcome all of the limitations of the SM code that you would need to.

Sorry to sound so negative - but as I say, I am far from an expert in this kind of maths - there may be others here who have some insights into possible ways around the problems that I've highlighted.

martinvicanek · 21 Oct 2011

Hi Michael,

interesting. Not long ago I was considering LPC as a means of manipulating formats in voice processing. I finally dropped it because I figured there would be too much latency in a real-time application. I agree with Trog that it would be tough to implement the Levinson-Durbin algorithm in SM, although I think it would not be entirely impossible. The complexity is O(N^2) which means you have (only) two nested loops. You could unroll the loops for a fixed N. Although you might find it easier to design it from scratch rather than try and adapt the code from musicdsp. It is by no means a straight-forward job but I think it could be done. What do you want to use it for anyway?

Michael · 22 Oct 2011

Thanks for your swift replies!

I was afraid that the ifs and loops make things difficult....

Anyway, I try to use Synthmaker to extract vowel data from speech.

Maybe a little explanation of my thinking would be good: Linear prediction, LP, is widely used in speech sciences to extract formant data (vowel data), and it does it well .

The voice can be thought as a voice source (vocal folds) and an all-pole filter (vocall tract). Theory behind this is called "linear source-filter theory".

Although voice source and vocal tract also couple non-linearly, LP works well in many cases.

LP is also used to compress speech data, for example it is used in mobile phone data packing algorithms.

Levinson method is used to solve a matrix that is of "Toeplitz" form (matrix in which each descending diagonal from left to right is constant). When using autocorrelation method for formant analysis, a Toeplitz matrix is created. Levinson figured out to simplify the calculation for Toeplitz matrix using less time and computational power.

I am no mathematician or coder, just a musician doing some research on singing, and I would like to create a software to aid me.

The FFT, Float Array Section, and Float Max componensts work nicely when getting the intensities of vocal fold vibration harmonics. FFT method is seemingly not accurate enough for formant tracking (=finding out the vocal tract resonances, vocal tract being the space between vocal folds in the larynx, and lips), for example, the voice source harmonics interfere. The peaks detected by Float Array Max are not necessarily precisely the vocat tract resonances (formants).

Michael

infuzion · 24 Oct 2011

Hi Michael, & welcome to SM, There is a C++ to SM mini-tutorial in the WiKi, but as Trog said you have to rework your code. How many "reflection coefficients K" per sample are there please?

CoreStyler · 24 Oct 2011

This code is really hard-coded. 😠

martinvicanek · 3 Mar 2012

I finally got round to implementing this algo, and it turns out to work quite well. I have been struggling some parts of it and would like to share what I have learned along the road. The method analyses the audio input in blocks of 1024 samples (23 ms @ 44.1 kHz sample rate) and provides filter coefficients at the end of each frame. These can be used to manipulate a bearer signal, which actually gives a nice vocoder. The inherent 23 ms latency is not as bad as one would expect. There are design rules as to how many coefficients are required for good intelligibility, and of course there is a trade-off regarding CPU demand. The rule of thumb is that you need 1 formant per kHz bandwidth, which leaves the following alternatives: 1. Brute force: allow for 40 or so filter coefficients (20 formants) to cover the entire audible range. This would not seem to be an entirely hopeless route, however bear in mind the O(n^2) complexity and also the potential instability of such a high-order filter. 2. Downsample to reduce the requirement (as in a telephone transmission). I have done this, it works, however it sounds, uhm well, like a telephone. 3. Warp the frequency scale to give more importance to the relevant range (500 Hz to 4 kHz). This is a bit tricky and adds some computational overhead, but the results are quite convincing already with as few as 9 coefficients (4 formants). The schematic below implements the third alternative. You can change the warp parameter to see its effect (0 means no warping). If you use different warp for analysis and synthesis, respectively, you get Klingon or chipmunk sounds. It is interesting to note that the Lewinson-Durbin recursion is, despite its O(n^2) complexity, not the most CPU demanding part in the processing chain because it is executed only each 1024 samples. Two obvious improvements include (a) overlapping frames and (b) more formants. Have fun!

tektoog · 3 Mar 2012

whao... this is really fun ! and sounds really good... Surely way above my coding skills.... but I will study this... for sure... Good job there !

stw · 4 Mar 2012

martinvicanek wrote: !

Just took a short look at it and i guess its another masterpiece from you. Looking forward for deeper diving into it! Again some awesome work Martin! And yet again...Thanks for sharing!! 😍

aliasant · 4 Mar 2012

Very nice and quite the mindblower!

martinvicanek · 4 Mar 2012

Thanks for your kudos guys! 😊 I am just rebuilding what others have invented long ago. I am having fun anyway.

martinvicanek wrote: Two obvious improvements include (a) overlapping frames and (b) more formants. Have fun!

Regarding (b) I have managed to increase the number of filter coefficients from 9 to 15, resulting in up to 7 formants. That turned out to be a substantial improvement. The result is actually close to the "ultimate vocoder", one that not only distinguishes between "oooh" and "aaah", but so faithfully reproduces the spectral envelope that it is possible to distinguish the singer's identity. Indeed, to recognize the synthesized output as my own voice, even though it is merely a filtered saw, is a fascinating, even a bit scary experience. More fun below 😉

aliasant · 4 Mar 2012

martin. That does sound incredible. I have my own voice on a track in iTunes and as you said. It sounds like my voice even though its "just" a saw.

but. Theres some -1#IND's screwing it up and sometimes even killing the ASIO stream. I had a look at the synthesis module and found them on xmmreg 0,1 and 2

Another thing. I tried to, just to see the cpu rise, change the hops from 1024 to 256 and cpu climbed from 1.8 to 2.0 barely noticable but sound was a lot worse. Then I tried ( anyone getting the pattern? 😉 ) to find anything that i could alter to make it work at 256 but found nothing that made any sense in changing. I guess a filter freq somewere that needs to be 4 times "faster"? The Q is. Could it work using smaller hops and if so, would we really get a smoother sound? What about the other way around?

martinvicanek · 4 Mar 2012

aliasant wrote: Theres some -1#IND's screwing it up and sometimes even killing the ASIO stream. I had a look at the synthesis module and found them on xmmreg 0,1 and 2

That's interesting. How do you analyze that? I also noticed an odd accumulation of glitches over time. After a few minutes the sound gets so bad that I have to restart the ASIO stream then everything is OK again. I thought that problem was particular to my installation (ASIO4ALL and my new TASCAM US-144 MKII do not seem to like each other.). It does help somewhat to increase the ASIO buffer size and switch off the FFT displays. Would be interesting to have other's feedback on that, and maybe find a remedy someone?

aliasant wrote: [...] change the hops from 1024 to 256

Shorter frames have less spectral detail, you need a certain time for the formants to evolve (there's Heisenberg's uncertainty principle again). Hop(1024) seems to be a good compromise between spectral detail and latency, and the corresponding value of 23 ms matches the frame size commonly used in the literature (20 to 25 ms). Is 256 a magic number in SM's architecture? Who knows. 🙂 Unfortunately I have hardcoded the 1024 hop size, so if you change it in my OSM, you have to do it in all modules consistently, including the constants for the Hanning window. But I doubt that that has anything to do with the problems above.

aliasant · 4 Mar 2012

Heres your osm but with the added xmm reg watcher. DWB made this some time ago and its a nice tool for sure!

I played my voice again and after a little time I start getting these high amped annomalies and after another 15 sec or so ASIO died and xmm 0,1 and 2 had #IND's

I dont think 256 is a magic number other then that it is a very common buffer value but as usual Im probably wrong 🙂

martinvicanek · 4 Mar 2012

aliasant wrote: Heres your osm but with the added xmm reg watcher. DWB made this some time ago and its a nice tool for sure!

Thanks for that! Although I don't quite understand what it does exactly. 🙂 Anyway, I had it running for some 10 minutes without #INDs. Oh, and I solved my problems (see my last post) by reinstalling ASIO4ALL, so now all is good here. That does not help you, I am aware. Is there a difference between the first (9 pole) and the second (15 pole) OSM?

aliasant · 4 Mar 2012

martinvicanek wrote: aliasant wrote: Heres your osm but with the added xmm reg watcher. DWB made this some time ago and its a nice tool for sure!

Thanks for that! Although I don't quite understand what it does exactly. 🙂 Anyway, I had it running for some 10 minutes without #INDs. Oh, and I solved my problems (see my last post) by reinstalling ASIO4ALL, so now all is good here. That does not help you, I am aware. Is there a difference between the first (9 pole) and the second (15 pole) OSM?

I didnt have these problems with the first 9pole version. Not at all. Would be nice if others could test it and see if they get IND's or glitches.

I dont use ASIO4all since I have a real soundcard and I havnt had any similar problems before so theres something going on in your code and probably something happend when you added those extra poles.

martinvicanek · 4 Mar 2012

aliasant wrote: I didnt have these problems with the first 9pole version. Not at all. [...] so theres something going on in your code and probably something happend when you added those extra poles.

OK, that's a strong indication that the filter gets unstable with 15 poles. That would not be surprising, just weird that it's not happening here.

aliasant wrote: Would be nice if others could test it and see if they get IND's or glitches.

Yes, definitely! Meanwhile here is a quick hack for 12 poles, if you care to try.

aliasant · 4 Mar 2012

Crap.

I got the same problems so i downloaded your first one again and .....

😊 SORRY.

I get the same in that one too. Not sure why I didnt notice it the first time. The only difference that time was that I played music ( Depeche Mode , critical detail!....) and now I play my recorded voice and that has many pauses and music doesnt really have that.

The loud amped artifacts sounds like res peaks that sort of appears, explodes and disaperas really fast.

Could it be something with the feedback in the filters going berserk?

I just tried to connect my "tapeplayer" and loaded a drumloop then played that into the vocoder instead of an external ASIO signal. If i pause the drumloop while still holding a midinote I get a loud crack and gone is the audio. I have to deactivate ASIO and reactivate again to get it to work again. If i connect the tapeplayer directly to the amp after the voc alls well.

Uploaded this osm for testing.

martinvicanek · 4 Mar 2012

Thanks aliasant, that helps. Now I can reproduce it here, too. So far I have been feeding the analyser with a mic, so there was always some noise there. Will investigate why the schematic chokes on silence. The good news is that it is NOT caused by the increase from 9 to 15 poles. 😉

aliasant · 4 Mar 2012

martinvicanek wrote: Thanks aliasant, that helps. Now I can reproduce it here, too. So far I have been feeding the analyser with a mic, so there was always some noise there. Will investigate why the schematic chokes on silence. The good news is that it is NOT caused by the increase from 9 to 15 poles. 😉

Great ! 🙂 Carry on then...

btw. I made a Vocoder a couple of years ago and the build was huge. A huge load of filters and monopacks.... This sounds so much more promising at a 50th of the size.

martinvicanek · 4 Mar 2012

Here is a quick and dirty fix: place a denormal remover in the voice path before the autocorrelator and set noise level to -120 dB. (you have to go that high because the autocorrelator computes squares which would vanish otherwise) That way the subsequent stages have always something to chew on which will keep them from going berserk.

aliasant · 4 Mar 2012

martinvicanek wrote: Here is a quick and dirty fix: place a denormal remover in the voice path before the autocorrelator and set noise level to -120 dB. (you have to go that high because the autocorrelator computes squares which would vanish otherwise) That way the subsequent stages have always something to chew on which will keep them from going berserk.

I tried that and the tapeplayer behaved better but I still hat severe dropouts and ASIO shutdowns depending on what I fed it.

Then I looked in the Lewinson-Durbin code and the first line of code tries to add a very very small number to avoid dividing by zero. I dont think that small a number can be used like that so I changed it from 1e-19 to 1e-8

Now it seems much much more stable.

tektoog · 5 Mar 2012

Great improvment with the 15 poles version; It sounds rounder and more realisitic... I have the same problem as Aliasant... All goes really well with external sound input thru Asio (I also use a real soundcard routing, no virtual ASIO), but I also experience ASIO crashes when using the wavplayer. After feeding it with different wavfiles; it seems that the crash occurs when the sounds has lot of consonnance (SSS sounds, I don't know how to put that) Maybe a de-esser somehere in the chain could help...I'm throwing that without really understanding what's going on... but that might be a track to look into...

Try this wav with the player, in all three versions of your schematic, I have an ASIO crash with it...

I tried the quick fix you suggested, but that doesn't solve the ASIO crash... This is definitely one of the most exciting project I've seen... Once again... Amazing job...It's opening a wide new field to look into for cool new projects ...

martinvicanek · 5 Mar 2012

Thanks, tektoog! Looks like I have some more work. 😉 Ah, these recursive filters are so unforgiving! You say just one wrong word (#IND) to them and they never speak to you again! They are a curse, probably that's why they are called recursive, haha.

aliasant · 5 Mar 2012

Thanx Tektog. I tested your wav and even with my little fix changing the supersmall number to a small one I get ASI= crash almost emediatly.

What about the "Pre-Emphasis" filter then? That one is dangerous for sure. I placed a Tahn Softclipper just after it to tame it and no more ASIO crash using your testwav. That together with the divby small number fix seems to do miracles 🙂 Now that the highend of the signal is prevented from going berzerk I can hear the Sssssss distort when its hits the sat but it wont crash. btw. I did remove both the "Denormal Removers"

Heres that version for testing.

stw · 5 Mar 2012

I could reproduce these #IND errors only with the 15 pole so far. Chnaging the warp factor from 0.75 to 0.7 seemed to work at first but as i found out it only maxes the time until the error occurs. Maybe this help tracing the error. However i think there's another error in the filter module. You clip the (stream) "out" variable which feeds back into the filter. That leads to some nasty results depending on the input. I guess the clip (if needed there anyway) should be seperated from the filter code or you have to introduce a temp output for further use in the filter.

aliasant · 5 Mar 2012

I did more testing replacing the Emphasis Filter with a bandpass at about 900Hz and that works much better too even without a sat after but its still a "quick and dirty fix"

martinvicanek · 5 Mar 2012

stw, the clippers in the filter are there as a last resort to prevent a runaway in case the filter coefficients happen to produce poles outside the unit cirlce for some frame. For that frame, the sound won't be nice, I know, but it should normally not happen anyway (provided you drive the whole schematic at a reasonable dB level). If you place the clippers outside the filter then it won't prevent the runaway and the filter will get locked up. I think that the problems reported are not due to filter coefficients that make the filter unstable, but rather INF filter coefficients for some frames. The clippers won't save the situation in that case, that's why we are having these lock-ups. I need to understand EXACTLY what's happening and devise a strategy from there, as opposed to inventing ad hoc cures for symptoms. I am really grateful for the feedback and trouble reports so far. Is there a way to check for INFs and NANs etc. in stream, like bool expressions that could be used for masks?

stw · 5 Mar 2012

martinvicanek wrote: Is there a way to check for INFs and NANs etc. in stream, like bool expressions that could be used for masks?

check this useful toy from MyCo: A geek toy: Stream inspector

stw · 5 Mar 2012

martinvicanek wrote: stw, the clippers in the filter are there as a last resort to prevent a runaway in case the filter coefficients happen to produce poles outside the unit cirlce for some frame.

Hmmm... my filter design knowledge is very limited so i just throw in my 2 cents to let you shed some light on me... 😉 How could you prevent poles from running out of the unit circle by simply clipping the filter out? I guessed if the coefficiants are wrong and these wrong poles are generated the filter becomes unstable anyway? And on the other hand in guessed that levels outside the normal range in IIR filters (especially if you deal with resonance) are nothing uncommon?

martinvicanek wrote: If you place the clippers outside the filter then it won't prevent the runaway and the filter will get locked up.

I never could create that scenario. But i can safely reproduce glitches resulting from that clipping if the input contains rich higher harmonics (e.g. the Vocal file from tektoog). I'd presume the glitches result from the created pulse from the clipper.

martinvicanek wrote: I think that the problems reported are not due to filter coefficients that make the filter unstable, but rather INF filter coefficients for some frames. The clippers won't save the situation in that case, that's why we are having these lock-ups.

I found some NaN's in the G(ain) values but didn't take a deeper look to that.

EDIT: Anyway... i sorted it out. The #IND errors came from some negative sqrt calcs in the LD algo. I limited that and that solved the prob on my side.

EDIT EDIT: What the heck... i tested it for several minutes with all known conditions which caused the errors. And right after upping the osm and changing input it appears again. So i take my "solution" (which is one for at least that one value) works properly. Sorry folks... 😠

aliasant · 5 Mar 2012

stw.

That seemed to work better but when i move the Pitch on the tapeplayer between 1 and 3 ASIO breaks down. Maybe not emediatly but when you hit the right pitch at the right time it will and it didnt take me more then a few turns back and forth to break it.

I think that pitch is at 2.16 ( using samplerate 48kHz if that matters)

If I leave it playing the wave at that pitch I can see the "Warped Linear" Graph peak oddly at about 12kHz and even dive deep and then come back.

Edit: Theres another one at pitch 0.55

Could it be just one of the windows/channels that has a faulty codeline ?

stw · 5 Mar 2012

aliasant wrote: Could it be just one of the windows/channels that has a faulty codeline ?

After some more testing i guess i produced exactly what martin predicted when he put in the clipper. I commented the clipper out in my last version and the errors occured in different cases. I now put the clipper back and the #IND errors are gone. At least i found a setting in which the #IND error comes up immediatly when the clipper is out. And other wise you clearly hear the glitches if the clipper is in. Since i'm running out of time for today maybe that osm can be used for further testing...

trogluddite · 5 Mar 2012

martinvicanek wrote: Is there a way to check for INFs and NANs etc. in stream, like bool expressions that could be used for masks?

There is, but only in assembly - SM code has no equivalent. See HERE the definition for the cmpss opcode - by setting the third operand to 3 or 7 you can test for ORD(ered) or UNORD(ered) relationships. Two values are said to be unordered if either (or both) of them are NaN (Not A Number - usually #INF[inity] or #IND[eterminate])- i.e. you can't arrange two numbers in order if one of them isn't a number at all. Conversely, if both inputs are valid, the comparison is said to be ordered.

So this little chunk of code turns NaNs into zero (so far as I have tested)

//assume input value is in xmm0 movaps xmm1,xmm0; cmpps xmm1,xmm1,7; //True if xmm1 is not NaN andps xmm0,xmm1; //Set NaNs to zero, else allow the value through //value still in xmm0

This only works because of the bitwise AND - most opcodes will propagate a NaN straight through even if it sometimes seems counter-intuitive. For example multiplying a NaN by zero will not output zero - it will output a NaN. This ensures that even very lengthy calculations will reliably show whenever there has been an invalid step in the calculation, no matter how far back. This is the reason why it only takes one NaN to lock up the audio - the NaN just propagates ad infinitum around the filter feedback loops.

trogluddite · 5 Mar 2012

trogluddite wrote: SM code has no equivalent

Changed my mind! 😊 I just remembered that any comparison with a NaN (except "unord") is, by definition, false - including comparing a NaN to another NaN, or even itself. So, all we need to do is this...

NanTest = (Value != Value); //True if "Value" is a NaN ValidTest = (Value == Value); //True if "Value" is a real number

Looks nonsensical - but it works (unlike my memory and reasoning after a long day at work!) So to make all NaNs into zeros...

Value = Value & (Value == Value);

martinvicanek · 5 Mar 2012

Thank you all for the overwhelming amount of responses, all valid and helpful! Hey but let's not overdo it, there must be other things in life apart from LPC. 😀 I'll pick it up later this week, I am confident that we will finally tame this Levinson guy.

Value = Value & (Value == Value);

I like that! I think from now on this will be my favorite line of code.

DigitalWhiteByte · 5 Mar 2012

😀

trogluddite · 5 Mar 2012

martinvicanek wrote: Value = Value & (Value == Value); I like that! I think from now on this will be my favorite line of code.

He he, straight from the Lewis Carroll* Academy of Programming! 😉

All that thinking has left me with one of those annoying pop tunes going round and round in my head though (recursively, you might say); sure you know the one - goes a little like this...

...NaN, NaN, NaNaNaN, NaN, NaN, NananaNaN, NaN, NaN, NaNaNaN, NaN.... ...think it must be by one of those annoying #INDy bands! 😕

* Trivia fact of the day: As well as the Alice in Wonderland stories, and surreal verse, Carroll (Charles Dodgson) really did publish some serious academic tomes such as Symbolic Logic Parts I & II, Curiosa Mathematica...

nix · 5 Mar 2012

hehe

martinvicanek · 10 Mar 2012

Here is an update. I have had a closer look to trace the various glitches back to their origin. The problem occurs in the Levinson-Durbin routine when one of the so-called reflection coefficients k_i falls outside the (-1,1) interval. Theoretically this should never happen, but with single precision it does occur sometimes, and as a result all subsequent k_i are completely off, as are the filter coefficients. In particular, the variable EF can become negative, which results in a NaN for the square root in G (good catch, stw!), which in turn will break the stream when fed into the synthesis filter. So the trick is to keep the k_i within the (-1,1) interval. This can be achieved simply by increasing VERY slightly the diagonal matrix elements r0, which is more or less equivalent to adding white noise at -80 dB relative level to the modulator signal. So I am attaching the OSM with this fix. I tested thoroughly with a wide range of input signals and could not produce any instability. (The most challenging input is a sine - this would break all my previous schematics but not this one.) Appreciate your feedback!

tektoog · 11 Mar 2012

Hey Martin, It seems you nailed it.... I quickly tried different wavs with the player (havn't tried a sine) and it went flawless... I'll test it extensively later... Good job man !

RunBeerRun2 · 11 Mar 2012

I'm piecing together a couple of Martin's things, it's an autotune with autotuning vocoder as well. So far I've just added stuff to my toolbox and got it basically set up and working, it't at 24% cpu on my quad core so I'll try to widdle it down some.