Autocorrelation and levinson recursion (LPC)

[Michael]

Hi!

I found a code written in C (I think) that is supposed to do linear predictive coding using autocorrelation and Levinson-Durbin recursion. Is it possible to translate this code to Synthmaker Code?

Grateful for any opinions,
Michael

Here´s the code. (http://www.musicdsp.org/showone.php?id=137):



//find the order-P autocorrelation array, R, for the sequence x of length L and warping of lambda
//wAutocorrelate(&pfSrc[stIndex],siglen,R,P,0);
wAutocorrelate(float * x, unsigned int L, float * R, unsigned int P, float lambda)
{
double * dl = new double [L];
double * Rt = new double [L];
double r1,r2,r1t;
R[0]=0;
Rt[0]=0;
r1=0;
r2=0;
r1t=0;
for(unsigned int k=0; k<L;k++)
{
Rt[0]+=double(x[k])*double(x[k]);

dl[k]=r1-double(lambda)*double(x[k]-r2);
r1 = x[k];
r2 = dl[k];
}
for(unsigned int i=1; i<=P; i++)
{
Rt[i]=0;
r1=0;
r2=0;
for(unsigned int k=0; k<L;k++)
{
Rt[i]+=double(dl[k])*double(x[k]);

r1t = dl[k];
dl[k]=r1-double(lambda)*double(r1t-r2);
r1 = r1t;
r2 = dl[k];
}
}
for(i=0; i<=P; i++)
R[i]=float(Rt[i]);
delete[] dl;
delete[] Rt;
}

// Calculate the Levinson-Durbin recursion for the autocorrelation sequence R of length P+1 and return the autocorrelation coefficients a and reflection coefficients K
LevinsonRecursion(unsigned int P, float *R, float *A, float *K)
{
double Am1[62];

if(R[0]==0.0) {
for(unsigned int i=1; i<=P; i++)
{
K[i]=0.0;
A[i]=0.0;
}}
else {
double km,Em1,Em;
unsigned int k,s,m;
for (k=0;k<=P;k++){
A[0]=0;
Am1[0]=0; }
A[0]=1;
Am1[0]=1;
km=0;
Em1=R[0];
for (m=1;m<=P;m++) //m=2:N+1
{
double err=0.0f; //err = 0;
for (k=1;k<=m-1;k++) //for k=2:m-1
err += Am1[k]*R[m-k]; // err = err + am1(k)*R(m-k+1);
km = (R[m]-err)/Em1; //km=(R(m)-err)/Em1;
K[m-1] = -float(km);
A[m]=(float)km; //am(m)=km;
for (k=1;k<=m-1;k++) //for k=2:m-1
A[k]=float(Am1[k]-km*Am1[m-k]); // am(k)=am1(k)-km*am1(m-k+1);
Em=(1-km*km)*Em1; //Em=(1-km*km)*Em1;
for(s=0;s<=P;s++) //for s=1:N+1
Am1[s] = A[s]; // am1(s) = am(s)
Em1 = Em; //Em1 = Em;
}
}
return 0;
}

[trogluddite]

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]

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]

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]

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]

This code is really hard-coded.

[martinvicanek]

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]

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]

!

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]

Very nice and quite the mindblower!

[martinvicanek]

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

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]

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]

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?

[...] 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]

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

LPC_v3 with reg watcher.osm

(141.89 KiB) Downloaded 65 times

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]

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?