Technical discussions related to Speech Coding (all itu and other vocoders, ACELP, CELP, AMR, etc)
|
Hi I am working on G729 implementation in c54x. Can any one give me algorithm to find log2 for fixed point arithmetic. it will be very helpful for me thanx gautam |
|
|
|
> Message: 2 > Date: Tue, 6 Aug 2002 05:14:27 -0700 (PDT) > From: gautam moharir <> > Subject: logarithm algorithm > > Hi > I am working on G729 implementation in c54x. > Can any one give me algorithm to find log2 for fixed > point arithmetic. You normalize the number to be in the form m*2^e m in [0.5,1) with log2(m*2^e)=log2(m) + e You can calculate log2(m) with a table interpolation. I annex as reference a matlab file to generate a table, the c code, and the DSP56000 code for a similar calculation (is loge instead of log2), I provide it as it is, without support (is really old stuff that I have not used in a while). It may help you as reference to adapt it for your implementation. The C codes includes quantification simulation, but I hope its readable, as I cannot provide the quantization and normalization routines used in it. For the table interpolation, a reference is in appendix B of: Robert C. Boman,"Integer Implementation of a Perceptual-Based Acoustic Analysis for Speech Recognition Robust to Additive and Convolutional Noise", http://www.icspat.com/db_area/1999/348.PDF I hope this help, Sara matlab code (.m file) ---- clear; NN=5;%number of elements of the table 2^NN+1 %generate the values xx=0.5:1/2^NN:1.0; table=log(xx); table=round(2^15*table)/2^15; %write the table fid = fopen('table_log.h','w'); fprintf(fid,['#define WP ' num2str(2^NN) '\n']); fprintf(fid,['#define WP2 ' num2str(2^(NN-1)) '\n']); fprintf(fid,['double tablog[' num2str(2^(NN-1)+1) '] = {\n']); fprintf(fid,'%1.20e,\n',table(1:length(table)-1)); fprintf(fid,'%1.20e\n',table(length(table))); fprintf(fid,'};\n'); fclose(fid); --- generated table (.h file) --- #define WP 32 #define WP2 16 double tablog[17] = { -6.93145751953125000000e-001, -6.32507324218750000000e-001, -5.75378417968750000000e-001, -5.21301269531250000000e-001, -4.70001220703125000000e-001, -4.21203613281250000000e-001, -3.74694824218750000000e-001, -3.30230712890625000000e-001, -2.87689208984375000000e-001, -2.46856689453125000000e-001, -2.07641601562500000000e-001, -1.69891357421875000000e-001, -1.33544921875000000000e-001, -9.84497070312500000000e-002, -6.45446777343750000000e-002, -3.17382812500000000000e-002, 0.00000000000000000000e+000 }; --- piece of c code that simulates algorithm --- int j,i,morder; double LOG2, c0, mz1, x1, x2, y1, y2, D; int e1, nn, TABADDR=0; if (gain==0) gain=pow(2.0,-15.0); /* normalization */ normN(gain, &mz1, &e1, 16); mz1=t_qnt(mz1,0); /* c0=log(mz1); */ nn=floor(mz1*WP); /* in the DSP the lower bits of x are masked and transferred to rx */ x1=(double) nn / (double) WP; y1=tablog[(TABADDR-WP2)+nn]; y2=tablog[(TABADDR-WP2)+nn+1]; D=y2-y1; D=t_qnt(D,0); c0=(mz1-x1); c0=t_qnt(c0,0); c0=c0* WP; c0=t_qnt(c0,1); c0=c0*D; c0=t_qnt(c0,1); c0=c0+y1; c0=t_qnt(c0,1); --- DSP56000 code (in assembler, draft not tested and not well commented, sorry) --- ;calculate log(mz1) ;Note: ;NN=5, WP=2^NN=32, WP2=2^NN-1=16 ;#LOGTAB = (address of the LOG table -WP2) First element of the 17 element table! move a,x1 ;put argument mz1 in x1 ;table interpolation move #>32,y0 ; y0=2^-10 mpy x1,y0,a #LOGTAB,r4 ; table address - WP2 in r4 ; a=mz1*2^-10, 5 high bits and "x1" in a1 ; x-"x1" in a0 (?) move a1,n4 ; extract addressing bits in n4 and n5 move n4,n5 move a0,a1 ; (x-"x1")*2^-10*2^16 in a (?) ; (it works fine because positive) asr a #LOGTAB+1,r5 ;inc table address in r5 ; Why a=a/2 ? a=(x-"x1")*2^-10*2^16*2^-1 ; =(x-"x1")*2^5=(x-"x1")*WP move a0,y0 ; (x-"x1")*WP in y0 move y:(r4+n4),x1 ; get "y1" in x1 move y:(r5+n5),a ; get "y2" in a sub x1,a ; "y2-y1" in a move a,x1 y:(r4+n4),a ; "y2-y1" in x1 ; get "y1" in a macr x1,y0,a y:(LOG2),x1 ; a= "y1"+(x-"x1")*WP*("y2-y1") ; get LOG2 in x1 |