On Monday, September 9, 2013 2:27:37 PM UTC-4, glen herrmannsfeldt wrote:> clay wrote: > > > > (snip) > > > Finding the LPCs directly in fixed point is somewhat tricky > > > since they often have a large dynamic range; however, an end > > > run around the problem can be had by finding the PARCORS via > > > the Leroux-Gueguen algorithm. This works very well in fixed point. > > > > How large is large? Floating point is designed for problems > > having a dynamic range of hundred of orders of magnitude, but > > it a reasonable relative error. > > > > The mass of a proton and the mass of the earth differ by over 50 > > orders of magnitude, the mass of the earth is known to about six > > digits and of the proton to about nine digits. Note that the relative > > uncertainties are about the same, where the absolute uncertainties > > are very different. > > > > Seems to me that in DSP 24 bits (about 7 decimal digits) is a large > > dynamic range. > > > > -- glenThe OP specifically wanted a fixed point implimentation, so I pointed him towards a way of getting his result. How large is large? That depends on the data? I didn't say it can't be done but rather it would be somewhat tricky. On the other hand the LG algo I mentioned is not tricky. One could read Le Roux and Gueguen, "A Fixed Point Computation of Partial Correlation Coefficients." and get more details. Clay
Fixed-point implementation of levinson durbin algorithm
Started by ●January 3, 2011
Reply by ●September 9, 20132013-09-09
Reply by ●September 9, 20132013-09-09
On 9/9/13 11:27 AM, glen herrmannsfeldt wrote:> Seems to me that in DSP 24 bits (about 7 decimal digits) is a large > dynamic range. >i think that there is some 56K Levinson-Durbin code laying around somewhere (it came from Dr. Bub from the previous millennium and i cleaned it up once). if prodded, i might be able to dig something up. if no one cares, i won't bother. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●September 9, 20132013-09-09
On 9/9/2013 2:19 PM, robert bristow-johnson wrote:> On 9/9/13 11:27 AM, glen herrmannsfeldt wrote: > >> Seems to me that in DSP 24 bits (about 7 decimal digits) is a large >> dynamic range. >> > > i think that there is some 56K Levinson-Durbin code laying around > somewhereThe original GSM 6.10 codec used Shur recursion. IIRC all math in there was designed for 12 bit machine. It is also possible to converge to LPCs directly from data using LMS SG algorithm. This is very straightforward, however it could take more computation compared to solution methods. Vladimir Vassilevsky DSP and Mixed Signal Designs www.abvolt.com
Reply by ●September 9, 20132013-09-09
On Monday, September 9, 2013 3:19:33 PM UTC-4, robert bristow-johnson wrote:> On 9/9/13 11:27 AM, glen herrmannsfeldt wrote: > > > > > Seems to me that in DSP 24 bits (about 7 decimal digits) is a large > > > dynamic range. > > > > > > i think that there is some 56K Levinson-Durbin code laying around > > somewhere (it came from Dr. Bub from the previous millennium and i > > cleaned it up once). if prodded, i might be able to dig something up. > > if no one cares, i won't bother. > > > > > > -- > > > > r b-j rbj> > > > "Imagination is more important than knowledge."I have both Durbin and the LG code (560000). It would take a little while to find and dust off - it is from 20 or so years ago. I likely have c code to go with it. Clay
Reply by ●September 10, 20132013-09-10
clay@claysturner.com writes:> On Monday, September 9, 2013 3:19:33 PM UTC-4, robert bristow-johnson wrote: >> On 9/9/13 11:27 AM, glen herrmannsfeldt wrote: >> >> >> >> > Seems to me that in DSP 24 bits (about 7 decimal digits) is a large >> >> > dynamic range. >> >> > >> >> >> i think that there is some 56K Levinson-Durbin code laying around >> >> somewhere (it came from Dr. Bub from the previous millennium and i >> >> cleaned it up once). if prodded, i might be able to dig something up. >> >> if no one cares, i won't bother. >> >> >> >> >> >> -- >> >> >> >> r b-j rbj> >> >> >> "Imagination is more important than knowledge." > > I have both Durbin and the LG code (560000). It would take a little > while to find and dust off - it is from 20 or so years ago. I likely > have c code to go with it.I could just write it by the time you and Robert go through your old floppies. -- Randy Yates Digital Signal Labs http://www.digitalsignallabs.com
Reply by ●September 10, 20132013-09-10
On 9/10/13 3:36 PM, Randy Yates wrote:> ... > I could just write it by the time you and Robert go through your old > floppies.sorry Randy. you need to be more explicit (that you want something) for me to get off my bum. ya gotta unwrap the lines. have phun. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge." PAGE 255 IF 0 /********************************************************************** * * * Durbin Algorithm for LPC Coefficients * * * * input: * * n = LPC order * * r -> r[0] to r[n] autocorrelation values * * * * output: * * a -> a[0] to a[n] LPC coefficients, a[0] = 1.0 * * k -> k[0] to k[n] reflection coefficients, k[0] = unused * * * * * **********************************************************************/ #define N 17 durbin(const int n, const double* r, double* a, double* k) { double a_temp[N], alpha, epsilon; /* n <= N = constant */ int i, j; k[0] = 0.0; /* unused */ a[0] = 1.0; a_new[0] = 1.0; /* unnecessary but consistent */ alpha = r[0]; for(i=1; i<=n; i++) { epsilon = r[i]; /* epsilon = a[0]*r[i]; */ for(j=1; j<i; j++) { epsilon += a[j]*r[i-j]; } a[i] = k[i] = -epsilon/alpha; alpha = alpha*(1.0 - k[i]*k[i]); for(j=1; j<i; j++) { a_temp[j] = a[j] + k[i]*a[i-j]; /* update a[] array into temporary array */ } for(j=1; j<i; j++) { a[j] = a_temp[j]; /* update a[] array */ } } } /* Less concise version with inside for() loops executed at least once in the same manner as Fortran version. */ durbin(const int n, const double* r, double* a, double* k) { double a_temp[N], alpha, epsilon; /* n <= N = constant */ int i, j; k[0] = 0.0; /* unused */ a[0] = 1.0; a_new[0] = 1.0; /* unnecessary but consistent */ a[1] = k[1] = -r[1]/r[0]; alpha = r[0]*(1.0 - k[1]*k[1]); for(i=2; i<=n; i++) { epsilon = 0.0; for(j=0; j<i; j++) { epsilon += a[j]*r[i-j]; } a[i] = k[i] = -epsilon/alpha; alpha = alpha*(1.0 - k[i]*k[i]); for(j=1; j<i; j++) { a_temp[j] = a[j] + k[i]*a[i-j]; /* update a[..] array first into temporary array */ } for(j=i-1; j>=1; j--) { a[j] = a_temp[j]; /* update a[..] array */ } } } (Nearly) original FORTRAN source program: [1] SUBROUTINE DURBIN(N,R,K,A) C C THIS ROUTINE USES THE DURBIN ALGORITHM TO C TRANSFORM THE AUTOCORRELATION COEFFICIENTS R(0) TO R(N) C TO THE REFLECTION COEFFICIENTS K(1) TO K(N) AND TO THE FILTER C COEFFICIENTS A(0) TO A(N) (A(0)=1). C REAL K(1..N),R(0..N),A(0..N) C C REAL A_TEMP(0..N),EPSILON C A(0) = 1.0 A_TEMP(0) = 1.0 K(1) = -R(1)/R(0) A(1) = K(1) ALPHA = R(0)*(1.0 - K(1)*K(1)) DO 400 I = 2,N EPSILON = 0.0 DO 100 J = 0,I-1 EPSILON = EPSILON + A(J)*R(I-J) 100 CONTINUE K(I) = - EPSILON / ALPHA A(I) = K(I) ALPHA = ALPHA*(1.0 - K(I)*K(I)) DO 200 J = 1,I-1 A_TEMP(J) = A(J) + K(I)*A(I-J) 200 CONTINUE DO 300 J = 1,I-1 A(J) = A_TEMP(J) 300 CONTINUE 400 CONTINUE RETURN END References [1] Papamichalis, Panos E., Practical Approaches to Speech Coding, Englewood Cliffs, NJ: Prentice - Hall, 1987. [2] Rabiner, L.R. and R.W. Schafer, Digital Processing of Speech Signals, Englewood Cliffs, NJ: Prentice - Hall, 1978. [3] Makhoul, John, "Linear Prediction: A Tutorial Review," Proceedings of the IEEE, Volume 63 (April,1975): pp. 561 - 580. ENDIF ; ; i/o file definitions ; file_data_in EQU y:$FFFF file_data_out EQU y:$FFFE LPC_order EQU 10 ; 10th order LPC analysis ORG x: r DS LPC_order+1 ; autocorrelation coefficients input ORG x: k DC 0.0 ; define k[0] = 0.0 (unused) DS LPC_order ; reflection coefficients output ORG y: a_coef DC 0.99999988 ; define a[0] ~= 1.0 DS LPC_order ; LPC filter coef array ORG x: a_temp DC 0.99999988 ; define a_temp[0] ~= 1.0 DS LPC_order ; temporary LPC filter coef array ; ; The following registers correspond to the C variables: ; ; r0 - r[i] r4 - a_coef[i] ; r1 - k[i-1] r5 - a_coef[j-1] ; r2 - loop counter r6 - a_coef[i-j+1] ; r3 - r[i-j+1], a_new[j-1] ; ; move #r,r0 ; r0 -> r[0] input array to algorithm move #k+1,r1 ; r1 -> k[1] output array of algorithm move #a_coef+1,n4 ; n4 -> a[1] move #a_temp,n3 ; n3 -> a_temp[0] move #LPC_order-1,n2 jsr durbin durbin move n4,r4 ; r4 -> a[1] move n4,r5 ; r5 -> a[1] move x:(r0)+,y1 ; r[0] move x:(r0)+,a ; r[1] move #2,n6 ; n6 = 2 ; signed division abs a a,b ; a = abs(r[1]), b = r[1] and #$FE,ccr do #24,_div1 div y1,a ; compute abs(r[1])/r[0] _div1 neg a a0,x0 eor y1,b a0,a ; sign bit, N = 1 if y1 and b are different signs tmi x0,a ; negate quotient if r[1] and r[0] have same sign ; a = -r[1]/r[0] move #<-1.0,b move a,x0 a,y:(r4)+ ; x0 = k[1], a[1] = k[1] = -r[1]/r[0] macr x0,x0,b ; b = -1.0 + k[1]^2 move #0,r2 ; r2 = 0 move b,x1 ; x1 = -(1.0 - k[1]^2) mpyr -x1,y1,b x0,x:(r1)+ ; b = r[0]*(1.0 - k[1]^2), write k[1] move (r5)- ; r5 -> a[0] move b,y1 ; alpha = r[0]*(1.0 - k[1]^2) do n2,_outer_loop ; for (i=2; i<=LPC_order; i++) { move r0,r3 ; r3 -> r[i] move (r0)+ ; r0 -> r[i+1] move (r2)+ ; r2 = i-1 nop clr a x:(r3)-,x0 y:(r5)+,y0 ; epsilon = 0; do r2,_loop1 ; for (j=0; j<i-1; j++) { // i-1 iterations mac x0,y0,a x:(r3)-,x0 y:(r5)+,y0 ; epsilon += a[j]*r[i-j]; _loop1 ; } mac x0,y0,a ; epsilon += a[i-1]*r[1]; ; signed division abs a a,b ; a = abs(epsilon), b = epsilon and #$FE,ccr do #24,_div2 div y1,a ; compute abs(epsilon)/alpha _div2 neg a a0,x0 eor y1,b a0,a ; sign bit, N = 1 if y1 and b are different signs tmi x0,a ; negate quotient if epsilon and alpha have same sign ; a = k[i] = -epsilon/alpha move #<-1.0,b move a,x0 a,y:(r4)+ ; x0 = k[i], a[i] = k[i] macr x0,x0,b r4,r6 ; b = -1.0 + k[i]^2, r6 -> a[i+1] move n4,r5 ; r5 -> a[1] move b,x1 ; x1 = -(1.0 - k[i]^2) mpyr -x1,y1,b x0,x:(r1)+ ; b = alpha*(1.0 - k[i]^2), store k[i] move (r6)-n6 ; r6 -> a[i-1] move b,y1 ; alpha <- alpha*(1.0 - k[i]^2) move y:(r6)-,y0 ; y0 = a[i-1] move y:(r5)+,a ; a = a[1] move n3,r3 ; r3 -> a_temp[0] do r2,_loop2 ; for (j=1; j<=i-1; j++) { // i-1 iterations macr x0,y0,a y:(r6)-,y0 ; // a[j] + k[i]*a[i-j], y0 = a[i-(j+1)] move (r3)+ ; // r3 -> a_temp[j] move a,x:(r3) y:(r5)+,a ; a_temp[j] = a[j] + k[i]*a[i-j] // a = a[j+1] _loop2 ; } move (r5)- ; r5 -> a[i-1] do r2,_loop3 ; for (j=i-1; j>=1; j--) { // i-1 iterations move x:(r3)-,y0 move y0,y:(r5)- ; a[j] = a_temp[j] _loop3 ; } nop ; r5 -> a[0] outer_loop ; } rts
Reply by ●September 10, 20132013-09-10
On 9/10/13 6:12 PM, robert bristow-johnson wrote:> > k[0] = 0.0; /* unused */ > > a[0] = 1.0; > a_new[0] = 1.0; /* unnecessary but consistent */should be a_temp[0] = 1.0; the C code was put into a comment of the asm program. i know the asm program worked, but not certain the C code in the comments worked, but i thought i tried it out once. obviously not with this "a_new" which should have been renamed lest the compiler complains. again, good luck with it, Randy. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●September 10, 20132013-09-10
robert bristow-johnson <rbj@audioimagination.com> writes:> On 9/10/13 3:36 PM, Randy Yates wrote: >> ... >> I could just write it by the time you and Robert go through your old >> floppies. > > sorry Randy. you need to be more explicit (that you want something) > for me to get off my bum. ya gotta unwrap the lines. > > have phun.Robert, I am not asking you for anything. Perhaps I was a bit too terse. What I intended to communicate is that I am so skilled in the art of DSP and such an awesome coder that, by the time you or Clay could FIND your old code, I could WRITE it! -- Randy Yates Digital Signal Labs http://www.digitalsignallabs.com
