Patch from Bruce:
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Fri, 1 Oct 2010 00:46:46 +0000 (00:46 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Fri, 1 Oct 2010 00:46:46 +0000 (00:46 +0000)
Generate const C tables for LSP codebook at compile time.
Replace mallocs with stack arrays.
ANSI C90 compatibility.
Conditional compile DUMP, give embedded systems a space break.
Conditional-define out unused VQ functions.
Initialize variables that are first passed by address with no
initialization, to keep the compiler from complaining.
Set GCC pedantic flag, we need to be that portable.
Get rid of any BSS or DATA segment data, there were a few unused static
variables. All compile-time initialized data is read-only and in the
text segment now.

git-svn-id: https://svn.code.sf.net/p/freetel/code@196 01035d8c-6547-0410-b346-abe4f91aad63

14 files changed:
codec2/src/Makefile
codec2/src/c2dec.c
codec2/src/c2enc.c
codec2/src/c2sim.c
codec2/src/defines.h
codec2/src/dump.c
codec2/src/lsp.c
codec2/src/nlp.c
codec2/src/pack.c
codec2/src/postfilter.c
codec2/src/quantise.c
codec2/src/quantise.h
codec2/unittest/Makefile
codec2/unittest/tnlp.c

index 66af79f456ae894a5c271f8837037e1cd22f275e..4254db8296ba97d5e67d6c3fa7e517fe7a754e64 100644 (file)
@@ -1,14 +1,27 @@
 CC=gcc
-CFLAGS=-g -Wall -I. -I../src -Wall -g -DFLOATING_POINT -DVAR_ARRAYS -O2
+CFLAGS=--pedantic -Wall -I. -I../src -Wall -g -DFLOATING_POINT -DVAR_ARRAYS -O2
 
 C2SIM_OBJ = sine.o nlp.o four1.o dump.o quantise.o lpc.o lsp.o phase.o \
-            pack.o postfilter.o interp.o codec2.o c2sim.o 
+            pack.o postfilter.o interp.o codec2.o c2sim.o codebook.o
 
 C2ENC_OBJ = sine.o nlp.o four1.o dump.o quantise.o lpc.o lsp.o phase.o \
-            pack.o postfilter.o interp.o codec2.o c2enc.o 
+            pack.o postfilter.o interp.o codec2.o c2enc.o codebook.o
 
 C2DEC_OBJ = sine.o nlp.o four1.o dump.o quantise.o lpc.o lsp.o phase.o \
-            pack.o postfilter.o interp.o codec2.o c2dec.o 
+            pack.o postfilter.o interp.o codec2.o c2dec.o codebook.o
+
+D=codebook
+CODEBOOKS= \
+  $D/lsp1.txt \
+  $D/lsp2.txt \
+  $D/lsp3.txt \
+  $D/lsp4.txt \
+  $D/lsp5.txt \
+  $D/lsp6.txt \
+  $D/lsp7.txt \
+  $D/lsp8.txt \
+  $D/lsp9.txt \
+  $D/lsp10.txt
 
 all: c2sim c2enc c2dec
 
@@ -21,8 +34,17 @@ c2enc: $(C2ENC_OBJ)
 c2dec: $(C2DEC_OBJ)
        $(CC) $(CFLAGS) $(C2DEC_OBJ) -o c2dec -lm
 
+codebook.c: generate_codebook $(CODEBOOKS)
+       ./generate_codebook $(CODEBOOKS) > codebook.c
+
+codebook.o: codebook.c
+       $(CC) -c $(CFLAGS) $< -o $@
+
+generate_codebook: generate_codebook.o
+       $(CC) $(CFLAGS) generate_codebook.o -o generate_codebook -lm
+
 %.o : %.c
        $(CC) -c $(CFLAGS) $< -o $@
 
 clean :
-       rm -f *.o *~ src/*~ c2enc c2dec c2sim
+       rm -f *.o *~ src/*~ c2enc c2dec c2sim generate_codebook src/codebook.c
index 3b876bcac035d6778eaa952d1a9cf1e0b6977c70..f64b38073a193cb8a02d116bbfea2a6506702ba5 100644 (file)
 #include <string.h>
 #include <errno.h>
 
+#define BITS_SIZE      ((CODEC2_BITS_PER_FRAME + 7) / 8)
+
 int main(int argc, char *argv[])
 {
-    static const int bitsSize = ((CODEC2_BITS_PER_FRAME + 7) / 8);
     void *codec2;
     FILE *fin;
     FILE *fout;
     short buf[CODEC2_SAMPLES_PER_FRAME];
-    unsigned char  bits[bitsSize];
+    unsigned char  bits[BITS_SIZE];
 
     if (argc != 3) {
        printf("usage: %s InputBitFile OutputRawSpeechFile\n", argv[0]);
@@ -66,7 +67,7 @@ int main(int argc, char *argv[])
 
     codec2 = codec2_create();
 
-    while(fread(bits, sizeof(char), bitsSize, fin) == bitsSize) {
+    while(fread(bits, sizeof(char), BITS_SIZE, fin) == BITS_SIZE) {
        codec2_decode(codec2, buf, bits);
        fwrite(buf, sizeof(short), CODEC2_SAMPLES_PER_FRAME, fout);
     }
index 8fd7c7778d20698bc7554a7d2b886139632fcb62..a2d02db65b645c447b59274ea5dd84111b326281 100644 (file)
 #include <string.h>
 #include <errno.h>
 
+#define BITS_SIZE      ((CODEC2_BITS_PER_FRAME + 7) / 8)
+
 int main(int argc, char *argv[])
 {
-    static const int bitsSize = ((CODEC2_BITS_PER_FRAME + 7) / 8);
     void *codec2;
     FILE *fin;
     FILE *fout;
     short buf[CODEC2_SAMPLES_PER_FRAME];
-    unsigned char  bits[bitsSize];
+    unsigned char  bits[BITS_SIZE];
 
     if (argc != 3) {
        printf("usage: %s InputRawspeechFile OutputBitFile\n", argv[0]);
@@ -70,7 +71,7 @@ int main(int argc, char *argv[])
     while(fread(buf, sizeof(short), CODEC2_SAMPLES_PER_FRAME, fin) ==
          CODEC2_SAMPLES_PER_FRAME) {
        codec2_encode(codec2, bits, buf);
-       fwrite(bits, sizeof(char), bitsSize, fout);
+       fwrite(bits, sizeof(char), BITS_SIZE, fout);
     }
 
     codec2_destroy(codec2);
index 72c5632fcd1c43510e8eb192aa2b70dda40114d3..1230425d350731bab8a56c2604e7eb9f9208aff8 100644 (file)
@@ -100,7 +100,7 @@ int main(int argc, char *argv[])
   float snr;
   float sum_snr;
 
-  int lpc_model, order;
+  int lpc_model, order = 0;
   int lsp, lsp_quantiser;
   float ak[LPC_MAX];
   COMP  Sw_[FFT_ENC];
@@ -114,7 +114,7 @@ int main(int argc, char *argv[])
   float bg_est;
 
   int   hand_voicing;
-  FILE *fvoicing;
+  FILE *fvoicing = 0;
 
   MODEL prev_model, interp_model;
   int decimate;
@@ -188,8 +188,10 @@ int main(int argc, char *argv[])
   }
 
   dump = switch_present("--dump",argc,argv);
+#ifdef DUMP
   if (dump) 
       dump_on(argv[dump+1]);
+#endif
 
   lsp = switch_present("--lsp",argc,argv);
   lsp_quantiser = 0;
@@ -241,7 +243,9 @@ int main(int argc, char *argv[])
     dft_speech(Sw, Sn, w); 
     two_stage_pitch_refinement(&model, Sw);
     estimate_amplitudes(&model, Sw, W);
+#ifdef DUMP 
     dump_Sn(Sn); dump_Sw(Sw); dump_model(&model);
+#endif
 
     /* optional zero-phase modelling */
 
@@ -249,7 +253,9 @@ int main(int argc, char *argv[])
        float Wn[M];                    /* windowed speech samples */
        float Rk[LPC_ORD+1];            /* autocorrelation coeffs  */
        
+#ifdef DUMP
        dump_phase(&model.phi[0], model.L);
+#endif
 
        /* find aks here, these are overwritten if LPC modelling is enabled */
 
@@ -261,13 +267,17 @@ int main(int argc, char *argv[])
        if (lpc_model)
            assert(order == LPC_ORD);
 
+#ifdef DUMP
        dump_ak(ak, LPC_ORD);
+#endif
        
        /* determine voicing */
 
        snr = est_voicing_mbe(&model, Sw, W, (FS/TWO_PI)*model.Wo, Sw_);
+#ifdef DUMP
        dump_Sw_(Sw_);
        dump_snr(snr);
+#endif
 
        /* just to make sure we are not cheating - kill all phases */
 
@@ -309,7 +319,9 @@ int main(int argc, char *argv[])
        aks_to_M2(ak, order, &model, e, &snr, 1); 
        apply_lpc_correction(&model, lpc_correction);
        sum_snr += snr;
+#ifdef DUMP
         dump_quantised_model(&model);
+#endif
     }
 
     /* option decimation to 20ms rate, which enables interpolation
@@ -380,8 +392,10 @@ int main(int argc, char *argv[])
   if (lpc_model)
       printf("SNR av = %5.2f dB\n", sum_snr/frames);
 
+#ifdef DUMP
   if (dump)
       dump_off();
+#endif
 
   if (hand_voicing)
     fclose(fvoicing);
index ef4899f70ac059bef53c310f177e823062ead3ea..76955e6526377d710723650c196396c9eb8a5e97 100644 (file)
@@ -81,4 +81,14 @@ typedef struct {
   int   voiced;                /* non-zero if this frame is voiced           */
 } MODEL;
 
+/* describes each codebook  */
+
+struct lsp_codebook {
+    int                        k;        /* dimension of vector        */
+    int                        log2m;    /* number of bits in m        */
+    int                        m;        /* elements in codebook       */
+    const float        *       cb;       /* The elements               */
+};
+extern const struct lsp_codebook lsp_cb[];
+
 #endif
index 2d187444838d219bb480fdbef9fa82b176d34b83..1c5d74ac9455e488608998a87d56bfd8144a660d 100644 (file)
@@ -32,6 +32,7 @@
 #include <string.h>
 #include <math.h>
 
+#ifdef DUMP
 static int dumpon = 0;
 
 static FILE *fsn = NULL;
@@ -400,3 +401,4 @@ void dump_E(float E) {
 
     fprintf(fE,"%f\n", 10.0*log10(E));
 }
+#endif
index feab4219abc7ecc38fe22b1d8170cc7a4aa620bc..158e9a192dd967b98c546fb472fe3b66a0808aad 100644 (file)
-/*---------------------------------------------------------------------------*\\r
-\r
-  FILE........: lsp.c\r
-  AUTHOR......: David Rowe\r
-  DATE CREATED: 24/2/93\r
-\r
-\r
-  This file contains functions for LPC to LSP conversion and LSP to\r
-  LPC conversion. Note that the LSP coefficients are not in radians\r
-  format but in the x domain of the unit circle.\r
-\r
-\*---------------------------------------------------------------------------*/\r
-\r
-#include "defines.h"\r
-#include "lsp.h"\r
-#include <math.h>\r
-#include <stdio.h>\r
-#include <stdlib.h>\r
-\r
-/*---------------------------------------------------------------------------*\\r
-\r
-  Introduction to Line Spectrum Pairs (LSPs)\r
-  ------------------------------------------\r
-\r
-  LSPs are used to encode the LPC filter coefficients {ak} for\r
-  transmission over the channel.  LSPs have several properties (like\r
-  less sensitivity to quantisation noise) that make them superior to\r
-  direct quantisation of {ak}.\r
-\r
-  A(z) is a polynomial of order lpcrdr with {ak} as the coefficients.\r
-\r
-  A(z) is transformed to P(z) and Q(z) (using a substitution and some\r
-  algebra), to obtain something like:\r
-\r
-    A(z) = 0.5[P(z)(z+z^-1) + Q(z)(z-z^-1)]  (1)\r
-\r
-  As you can imagine A(z) has complex zeros all over the z-plane. P(z)\r
-  and Q(z) have the very neat property of only having zeros _on_ the\r
-  unit circle.  So to find them we take a test point z=exp(jw) and\r
-  evaluate P (exp(jw)) and Q(exp(jw)) using a grid of points between 0\r
-  and pi.\r
-\r
-  The zeros (roots) of P(z) also happen to alternate, which is why we\r
-  swap coefficients as we find roots.  So the process of finding the\r
-  LSP frequencies is basically finding the roots of 5th order\r
-  polynomials.\r
-\r
-  The root so P(z) and Q(z) occur in symmetrical pairs at +/-w, hence\r
-  the name Line Spectrum Pairs (LSPs).\r
-\r
-  To convert back to ak we just evaluate (1), "clocking" an impulse\r
-  thru it lpcrdr times gives us the impulse response of A(z) which is\r
-  {ak}.\r
-\r
-\*---------------------------------------------------------------------------*/\r
-\r
-/*---------------------------------------------------------------------------*\\r
-\r
-  FUNCTION....: cheb_poly_eva()\r
-  AUTHOR......: David Rowe\r
-  DATE CREATED: 24/2/93\r
-\r
-  This function evalutes a series of chebyshev polynomials\r
-\r
-  FIXME: performing memory allocation at run time is very inefficient,\r
-  replace with stack variables of MAX_P size.\r
-\r
-\*---------------------------------------------------------------------------*/\r
-\r
-\r
-float cheb_poly_eva(float *coef,float x,int m)\r
-/*  float coef[]       coefficients of the polynomial to be evaluated  */\r
-/*  float x            the point where polynomial is to be evaluated   */\r
-/*  int m              order of the polynomial                         */\r
-{\r
-    int i;\r
-    float *T,*t,*u,*v,sum;\r
-\r
-    /* Allocate memory for chebyshev series formulation */\r
-\r
-    if((T = (float *)malloc((m/2+1)*sizeof(float))) == NULL){\r
-       fprintf(stderr, "not enough memory to allocate buffer\n");\r
-       exit(1);\r
-    }\r
-\r
-    /* Initialise pointers */\r
-\r
-    t = T;                             /* T[i-2]                       */\r
-    *t++ = 1.0;\r
-    u = t--;                           /* T[i-1]                       */\r
-    *u++ = x;\r
-    v = u--;                           /* T[i]                         */\r
-\r
-    /* Evaluate chebyshev series formulation using iterative approach  */\r
-\r
-    for(i=2;i<=m/2;i++)\r
-       *v++ = (2*x)*(*u++) - *t++;     /* T[i] = 2*x*T[i-1] - T[i-2]   */\r
-\r
-    sum=0.0;                           /* initialise sum to zero       */\r
-    t = T;                             /* reset pointer                */\r
-\r
-    /* Evaluate polynomial and return value also free memory space */\r
-\r
-    for(i=0;i<=m/2;i++)\r
-       sum+=coef[(m/2)-i]**t++;\r
-\r
-    free(T);\r
-    return sum;\r
-}\r
-\r
-\r
-/*---------------------------------------------------------------------------*\\r
-\r
-  FUNCTION....: lpc_to_lsp()\r
-  AUTHOR......: David Rowe\r
-  DATE CREATED: 24/2/93\r
-\r
-  This function converts LPC coefficients to LSP coefficients.\r
-\r
-\*---------------------------------------------------------------------------*/\r
-\r
-int lpc_to_lsp (float *a, int lpcrdr, float *freq, int nb, float delta)\r
-/*  float *a                   lpc coefficients                        */\r
-/*  int lpcrdr                 order of LPC coefficients (10)          */\r
-/*  float *freq                LSP frequencies in radians              */\r
-/*  int nb                     number of sub-intervals (4)             */\r
-/*  float delta                        grid spacing interval (0.02)            */\r
-{\r
-    float psuml,psumr,psumm,temp_xr,xl,xr,xm;\r
-    float temp_psumr;\r
-    int i,j,m,flag,k;\r
-    float *Q;                  /* ptrs for memory allocation           */\r
-    float *P;\r
-    float *px;                 /* ptrs of respective P'(z) & Q'(z)     */\r
-    float *qx;\r
-    float *p;\r
-    float *q;\r
-    float *pt;                 /* ptr used for cheb_poly_eval()\r
-                                  whether P' or Q'                     */\r
-    int roots=0;               /* number of roots found                */\r
-    flag = 1;                  \r
-    m = lpcrdr/2;              /* order of P'(z) & Q'(z) polynimials   */\r
-\r
-    /* Allocate memory space for polynomials */\r
-\r
-    Q = (float *) malloc((m+1)*sizeof(float));\r
-    P = (float *) malloc((m+1)*sizeof(float));\r
-    if( (P == NULL) || (Q == NULL) ) {\r
-       fprintf(stderr,"not enough memory to allocate buffer\n");\r
-       exit(1);\r
-    }\r
-\r
-    /* determine P'(z)'s and Q'(z)'s coefficients where\r
-      P'(z) = P(z)/(1 + z^(-1)) and Q'(z) = Q(z)/(1-z^(-1)) */\r
-\r
-    px = P;                      /* initilaise ptrs */\r
-    qx = Q;\r
-    p = px;\r
-    q = qx;\r
-    *px++ = 1.0;\r
-    *qx++ = 1.0;\r
-    for(i=1;i<=m;i++){\r
-       *px++ = a[i]+a[lpcrdr+1-i]-*p++;\r
-       *qx++ = a[i]-a[lpcrdr+1-i]+*q++;\r
-    }\r
-    px = P;\r
-    qx = Q;\r
-    for(i=0;i<m;i++){\r
-       *px = 2**px;\r
-       *qx = 2**qx;\r
-        px++;\r
-        qx++;\r
-    }\r
-    px = P;                    /* re-initialise ptrs                   */\r
-    qx = Q;\r
-\r
-    /* Search for a zero in P'(z) polynomial first and then alternate to Q'(z).\r
-    Keep alternating between the two polynomials as each zero is found         */\r
-\r
-    xr = 0;                    /* initialise xr to zero                */\r
-    xl = 1.0;                  /* start at point xl = 1                */\r
-\r
-\r
-    for(j=0;j<lpcrdr;j++){\r
-       if(j%2)                 /* determines whether P' or Q' is eval. */\r
-           pt = qx;\r
-       else\r
-           pt = px;\r
-\r
-       psuml = cheb_poly_eva(pt,xl,lpcrdr);    /* evals poly. at xl    */\r
-       flag = 1;\r
-       while(flag && (xr >= -1.0)){\r
-           xr = xl - delta ;                   /* interval spacing     */\r
-           psumr = cheb_poly_eva(pt,xr,lpcrdr);/* poly(xl-delta_x)     */\r
-           temp_psumr = psumr;\r
-           temp_xr = xr;\r
-\r
-        /* if no sign change increment xr and re-evaluate\r
-           poly(xr). Repeat til sign change.  if a sign change has\r
-           occurred the interval is bisected and then checked again\r
-           for a sign change which determines in which interval the\r
-           zero lies in.  If there is no sign change between poly(xm)\r
-           and poly(xl) set interval between xm and xr else set\r
-           interval between xl and xr and repeat till root is located\r
-           within the specified limits  */\r
-\r
-           if((psumr*psuml)<0.0){\r
-               roots++;\r
-\r
-               psumm=psuml;\r
-               for(k=0;k<=nb;k++){\r
-                   xm = (xl+xr)/2;             /* bisect the interval  */\r
-                   psumm=cheb_poly_eva(pt,xm,lpcrdr);\r
-                   if(psumm*psuml>0.){\r
-                       psuml=psumm;\r
-                       xl=xm;\r
-                   }\r
-                   else{\r
-                       psumr=psumm;\r
-                       xr=xm;\r
-                   }\r
-               }\r
-\r
-              /* once zero is found, reset initial interval to xr      */\r
-              freq[j] = (xm);\r
-              xl = xm;\r
-              flag = 0;                /* reset flag for next search   */\r
-           }\r
-           else{\r
-               psuml=temp_psumr;\r
-               xl=temp_xr;\r
-           }\r
-       }\r
-    }\r
-    free(P);                           /* free memory space            */\r
-    free(Q);\r
-\r
-    /* convert from x domain to radians */\r
-\r
-    for(i=0; i<lpcrdr; i++) {\r
-       freq[i] = acos(freq[i]);\r
-    }\r
-\r
-    return(roots);\r
-}\r
-\r
-/*---------------------------------------------------------------------------*\\r
-\r
-  FUNCTION....: lsp_to_lpc()\r
-  AUTHOR......: David Rowe\r
-  DATE CREATED: 24/2/93\r
-\r
-  This function converts LSP coefficients to LPC coefficients.  In the\r
-  Speex code we worked out a wayto simplify this significantly.\r
-\r
-\*---------------------------------------------------------------------------*/\r
-\r
-void lsp_to_lpc(float *freq, float *ak, int lpcrdr)\r
-/*  float *freq         array of LSP frequencies in radians            */\r
-/*  float *ak          array of LPC coefficients                       */\r
-/*  int lpcrdr         order of LPC coefficients                       */\r
-\r
-\r
-{\r
-    int i,j;\r
-    float xout1,xout2,xin1,xin2;\r
-    float *Wp;\r
-    float *pw,*n1,*n2,*n3,*n4;\r
-    int m = lpcrdr/2;\r
-\r
-    /* convert from radians to the x=cos(w) domain */\r
-\r
-    for(i=0; i<lpcrdr; i++)\r
-       freq[i] = cos(freq[i]);\r
-\r
-    if((Wp = (float *) malloc((4*m+2)*sizeof(float))) == NULL){\r
-       printf("not enough memory to allocate buffer\n");\r
-       exit(1);\r
-    }\r
-    pw = Wp;\r
-\r
-    /* initialise contents of array */\r
-\r
-    for(i=0;i<=4*m+1;i++){             /* set contents of buffer to 0 */\r
-       *pw++ = 0.0;\r
-    }\r
-\r
-    /* Set pointers up */\r
-\r
-    pw = Wp;\r
-    xin1 = 1.0;\r
-    xin2 = 1.0;\r
-\r
-    /* reconstruct P(z) and Q(z) by cascading second order polynomials\r
-      in form 1 - 2xz(-1) +z(-2), where x is the LSP coefficient */\r
-\r
-    for(j=0;j<=lpcrdr;j++){\r
-       for(i=0;i<m;i++){\r
-           n1 = pw+(i*4);\r
-           n2 = n1 + 1;\r
-           n3 = n2 + 1;\r
-           n4 = n3 + 1;\r
-           xout1 = xin1 - 2*(freq[2*i]) * *n1 + *n2;\r
-           xout2 = xin2 - 2*(freq[2*i+1]) * *n3 + *n4;\r
-           *n2 = *n1;\r
-           *n4 = *n3;\r
-           *n1 = xin1;\r
-           *n3 = xin2;\r
-           xin1 = xout1;\r
-           xin2 = xout2;\r
-       }\r
-       xout1 = xin1 + *(n4+1);\r
-       xout2 = xin2 - *(n4+2);\r
-       ak[j] = (xout1 + xout2)*0.5;\r
-       *(n4+1) = xin1;\r
-       *(n4+2) = xin2;\r
-\r
-       xin1 = 0.0;\r
-       xin2 = 0.0;\r
-    }\r
-    free(Wp);\r
-}\r
-\r
+/*---------------------------------------------------------------------------*\
+
+  FILE........: lsp.c
+  AUTHOR......: David Rowe
+  DATE CREATED: 24/2/93
+
+
+  This file contains functions for LPC to LSP conversion and LSP to
+  LPC conversion. Note that the LSP coefficients are not in radians
+  format but in the x domain of the unit circle.
+
+\*---------------------------------------------------------------------------*/
+
+#include "defines.h"
+#include "lsp.h"
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+/* Only 10 gets used, so far. */
+#define LSP_MAX_ORDER  20
+/*---------------------------------------------------------------------------*\
+
+  Introduction to Line Spectrum Pairs (LSPs)
+  ------------------------------------------
+
+  LSPs are used to encode the LPC filter coefficients {ak} for
+  transmission over the channel.  LSPs have several properties (like
+  less sensitivity to quantisation noise) that make them superior to
+  direct quantisation of {ak}.
+
+  A(z) is a polynomial of order lpcrdr with {ak} as the coefficients.
+
+  A(z) is transformed to P(z) and Q(z) (using a substitution and some
+  algebra), to obtain something like:
+
+    A(z) = 0.5[P(z)(z+z^-1) + Q(z)(z-z^-1)]  (1)
+
+  As you can imagine A(z) has complex zeros all over the z-plane. P(z)
+  and Q(z) have the very neat property of only having zeros _on_ the
+  unit circle.  So to find them we take a test point z=exp(jw) and
+  evaluate P (exp(jw)) and Q(exp(jw)) using a grid of points between 0
+  and pi.
+
+  The zeros (roots) of P(z) also happen to alternate, which is why we
+  swap coefficients as we find roots.  So the process of finding the
+  LSP frequencies is basically finding the roots of 5th order
+  polynomials.
+
+  The root so P(z) and Q(z) occur in symmetrical pairs at +/-w, hence
+  the name Line Spectrum Pairs (LSPs).
+
+  To convert back to ak we just evaluate (1), "clocking" an impulse
+  thru it lpcrdr times gives us the impulse response of A(z) which is
+  {ak}.
+
+\*---------------------------------------------------------------------------*/
+
+/*---------------------------------------------------------------------------*\
+
+  FUNCTION....: cheb_poly_eva()
+  AUTHOR......: David Rowe
+  DATE CREATED: 24/2/93
+
+  This function evalutes a series of chebyshev polynomials
+
+  FIXME: performing memory allocation at run time is very inefficient,
+  replace with stack variables of MAX_P size.
+
+\*---------------------------------------------------------------------------*/
+
+
+static float
+cheb_poly_eva(float *coef,float x,int m)
+/*  float coef[]       coefficients of the polynomial to be evaluated  */
+/*  float x            the point where polynomial is to be evaluated   */
+/*  int m              order of the polynomial                         */
+{
+    int i;
+    float *t,*u,*v,sum;
+    float T[(LSP_MAX_ORDER / 2) + 1];
+
+    /* Initialise pointers */
+
+    t = T;                             /* T[i-2]                       */
+    *t++ = 1.0;
+    u = t--;                           /* T[i-1]                       */
+    *u++ = x;
+    v = u--;                           /* T[i]                         */
+
+    /* Evaluate chebyshev series formulation using iterative approach  */
+
+    for(i=2;i<=m/2;i++)
+       *v++ = (2*x)*(*u++) - *t++;     /* T[i] = 2*x*T[i-1] - T[i-2]   */
+
+    sum=0.0;                           /* initialise sum to zero       */
+    t = T;                             /* reset pointer                */
+
+    /* Evaluate polynomial and return value also free memory space */
+
+    for(i=0;i<=m/2;i++)
+       sum+=coef[(m/2)-i]**t++;
+
+    return sum;
+}
+
+
+/*---------------------------------------------------------------------------*\
+
+  FUNCTION....: lpc_to_lsp()
+  AUTHOR......: David Rowe
+  DATE CREATED: 24/2/93
+
+  This function converts LPC coefficients to LSP coefficients.
+
+\*---------------------------------------------------------------------------*/
+
+int lpc_to_lsp (float *a, int lpcrdr, float *freq, int nb, float delta)
+/*  float *a                   lpc coefficients                        */
+/*  int lpcrdr                 order of LPC coefficients (10)          */
+/*  float *freq                LSP frequencies in radians              */
+/*  int nb                     number of sub-intervals (4)             */
+/*  float delta                        grid spacing interval (0.02)            */
+{
+    float psuml,psumr,psumm,temp_xr,xl,xr,xm = 0;
+    float temp_psumr;
+    int i,j,m,flag,k;
+    float *px;                 /* ptrs of respective P'(z) & Q'(z)     */
+    float *qx;
+    float *p;
+    float *q;
+    float *pt;                 /* ptr used for cheb_poly_eval()
+                                  whether P' or Q'                     */
+    int roots=0;               /* number of roots found                */
+    float Q[LSP_MAX_ORDER + 1];
+    float P[LSP_MAX_ORDER + 1];
+
+    flag = 1;                  
+    m = lpcrdr/2;              /* order of P'(z) & Q'(z) polynimials   */
+
+    /* Allocate memory space for polynomials */
+
+    /* determine P'(z)'s and Q'(z)'s coefficients where
+      P'(z) = P(z)/(1 + z^(-1)) and Q'(z) = Q(z)/(1-z^(-1)) */
+
+    px = P;                      /* initilaise ptrs */
+    qx = Q;
+    p = px;
+    q = qx;
+    *px++ = 1.0;
+    *qx++ = 1.0;
+    for(i=1;i<=m;i++){
+       *px++ = a[i]+a[lpcrdr+1-i]-*p++;
+       *qx++ = a[i]-a[lpcrdr+1-i]+*q++;
+    }
+    px = P;
+    qx = Q;
+    for(i=0;i<m;i++){
+       *px = 2**px;
+       *qx = 2**qx;
+        px++;
+        qx++;
+    }
+    px = P;                    /* re-initialise ptrs                   */
+    qx = Q;
+
+    /* Search for a zero in P'(z) polynomial first and then alternate to Q'(z).
+    Keep alternating between the two polynomials as each zero is found         */
+
+    xr = 0;                    /* initialise xr to zero                */
+    xl = 1.0;                  /* start at point xl = 1                */
+
+
+    for(j=0;j<lpcrdr;j++){
+       if(j%2)                 /* determines whether P' or Q' is eval. */
+           pt = qx;
+       else
+           pt = px;
+
+       psuml = cheb_poly_eva(pt,xl,lpcrdr);    /* evals poly. at xl    */
+       flag = 1;
+       while(flag && (xr >= -1.0)){
+           xr = xl - delta ;                   /* interval spacing     */
+           psumr = cheb_poly_eva(pt,xr,lpcrdr);/* poly(xl-delta_x)     */
+           temp_psumr = psumr;
+           temp_xr = xr;
+
+        /* if no sign change increment xr and re-evaluate
+           poly(xr). Repeat til sign change.  if a sign change has
+           occurred the interval is bisected and then checked again
+           for a sign change which determines in which interval the
+           zero lies in.  If there is no sign change between poly(xm)
+           and poly(xl) set interval between xm and xr else set
+           interval between xl and xr and repeat till root is located
+           within the specified limits  */
+
+           if((psumr*psuml)<0.0){
+               roots++;
+
+               psumm=psuml;
+               for(k=0;k<=nb;k++){
+                   xm = (xl+xr)/2;             /* bisect the interval  */
+                   psumm=cheb_poly_eva(pt,xm,lpcrdr);
+                   if(psumm*psuml>0.){
+                       psuml=psumm;
+                       xl=xm;
+                   }
+                   else{
+                       psumr=psumm;
+                       xr=xm;
+                   }
+               }
+
+              /* once zero is found, reset initial interval to xr      */
+              freq[j] = (xm);
+              xl = xm;
+              flag = 0;                /* reset flag for next search   */
+           }
+           else{
+               psuml=temp_psumr;
+               xl=temp_xr;
+           }
+       }
+    }
+
+    /* convert from x domain to radians */
+
+    for(i=0; i<lpcrdr; i++) {
+       freq[i] = acos(freq[i]);
+    }
+
+    return(roots);
+}
+
+/*---------------------------------------------------------------------------*\
+
+  FUNCTION....: lsp_to_lpc()
+  AUTHOR......: David Rowe
+  DATE CREATED: 24/2/93
+
+  This function converts LSP coefficients to LPC coefficients.  In the
+  Speex code we worked out a wayto simplify this significantly.
+
+\*---------------------------------------------------------------------------*/
+
+void lsp_to_lpc(float *freq, float *ak, int lpcrdr)
+/*  float *freq         array of LSP frequencies in radians            */
+/*  float *ak          array of LPC coefficients                       */
+/*  int lpcrdr         order of LPC coefficients                       */
+
+
+{
+    int i,j;
+    float xout1,xout2,xin1,xin2;
+    float *pw,*n1,*n2,*n3,*n4 = 0;
+    int m = lpcrdr/2;
+    float Wp[(LSP_MAX_ORDER * 4) + 2];
+
+    /* convert from radians to the x=cos(w) domain */
+
+    for(i=0; i<lpcrdr; i++)
+       freq[i] = cos(freq[i]);
+
+    pw = Wp;
+
+    /* initialise contents of array */
+
+    for(i=0;i<=4*m+1;i++){             /* set contents of buffer to 0 */
+       *pw++ = 0.0;
+    }
+
+    /* Set pointers up */
+
+    pw = Wp;
+    xin1 = 1.0;
+    xin2 = 1.0;
+
+    /* reconstruct P(z) and Q(z) by cascading second order polynomials
+      in form 1 - 2xz(-1) +z(-2), where x is the LSP coefficient */
+
+    for(j=0;j<=lpcrdr;j++){
+       for(i=0;i<m;i++){
+           n1 = pw+(i*4);
+           n2 = n1 + 1;
+           n3 = n2 + 1;
+           n4 = n3 + 1;
+           xout1 = xin1 - 2*(freq[2*i]) * *n1 + *n2;
+           xout2 = xin2 - 2*(freq[2*i+1]) * *n3 + *n4;
+           *n2 = *n1;
+           *n4 = *n3;
+           *n1 = xin1;
+           *n3 = xin2;
+           xin1 = xout1;
+           xin2 = xout2;
+       }
+       xout1 = xin1 + *(n4+1);
+       xout2 = xin2 - *(n4+2);
+       ak[j] = (xout1 + xout2)*0.5;
+       *(n4+1) = xin1;
+       *(n4+2) = xin2;
+
+       xin1 = 0.0;
+       xin2 = 0.0;
+    }
+}
+
index 193ca92109f0e25a8ff1900c464ba4be7b938051..8cf1fc7d70fc7a24338932fafdba960ab941c1ed 100644 (file)
@@ -60,7 +60,7 @@
 
 /* 48 tap 600Hz low pass FIR filter coefficients */
 
-float nlp_fir[] = {
+const float nlp_fir[] = {
   -1.0818124e-03,
   -1.1008344e-03,
   -9.2768838e-04,
@@ -249,13 +249,17 @@ float nlp(
     for(i=0; i<m/DEC; i++) {
        Fw[i].real = nlp->sq[i*DEC]*(0.5 - 0.5*cos(2*PI*i/(m/DEC-1)));
     }
+#ifdef DUMP
     dump_dec(Fw);
+#endif
     four1(&Fw[-1].imag,PE_FFT_SIZE,1);
     for(i=0; i<PE_FFT_SIZE; i++)
        Fw[i].real = Fw[i].real*Fw[i].real + Fw[i].imag*Fw[i].imag;
 
+#ifdef DUMP
     dump_sq(nlp->sq);
     dump_Fw(Fw);
+#endif
 
     /* find global peak */
 
index 2cbff4438a6dfdd2bae4a063da950b492be2ca2b..31551dfc43df82728ee6203c30252c1e419c0eca 100644 (file)
@@ -81,7 +81,8 @@ unpack(
  unsigned int          fieldWidth/* Width of the field in BITS, not bytes. */
  )
 {
-  unsigned int         field = 0;
+  unsigned int field = 0;
+  unsigned int t;
 
   do {
     unsigned int       bI = *bitIndex;
@@ -96,7 +97,7 @@ unpack(
   } while ( fieldWidth != 0 );
 
   /* Convert from Gray code to binary. Works for maximum 8-bit fields. */
-  unsigned int t = field ^ (field >> 8);
+  t = field ^ (field >> 8);
   t ^= (t >> 4);
   t ^= (t >> 2);
   t ^= (t >> 1);
index 6dad76b1e1ce6e6ea1a7a3704d73ff1cd4d98876..cd6c28b0a9193db1711ced6fd51735a55e1e00ae 100644 (file)
@@ -126,6 +126,8 @@ void postfilter(
              uv++;
          }
 
+#ifdef DUMP
   dump_bg(e, *bg_est, 100.0*uv/model->L);
+#endif
 
 }
index 90a251311ec38ac5b8b438ec9bc13f4b04b3099f..e386bf8905c8ac9d856db729c106d3cdc51ec817 100644 (file)
 #include "four1.h"
 
 #define LSP_DELTA1 0.01         /* grid spacing for LSP root searches */
-#define MAX_CB       20         /* max number of codebooks */
-
-/* describes each codebook  */
-
-typedef struct {
-    int   k;        /* dimension of vector                   */
-    int   log2m;    /* number of bits in m                   */
-    int   m;        /* elements in codebook                  */
-    char *fn;       /* file name of text file storing the VQ */
-} LSP_CB;
-
-/* lsp_q describes entire quantiser made up of several codebooks */
-
-#ifdef OLDER
-/* 10+10+6+6 = 32 bit LSP difference split VQ */
-
-LSP_CB lsp_q[] = {
-    {3,   1024, "../unittest/lspd123.txt"},
-    {3,   1024, "../unittest/lspd456.txt"},
-    {2,     64, "../unittest/lspd78.txt"},
-    {2,     64, "../unittest/lspd910.txt"},
-    {0,    0, ""}
-};
-#endif
-
-LSP_CB lsp_q[] = {
-    {1,4,16, "../unittest/lsp1.txt"},
-    {1,4,16, "../unittest/lsp2.txt"},
-    {1,4,16, "../unittest/lsp3.txt"},
-    {1,4,16, "../unittest/lsp4.txt"},
-    {1,4,16, "../unittest/lsp5.txt"},
-    {1,4,16, "../unittest/lsp6.txt"},
-    {1,4,16, "../unittest/lsp7.txt"},
-    {1,3,8, "../unittest/lsp8.txt"},
-    {1,3,8, "../unittest/lsp9.txt"},
-    {1,2,4, "../unittest/lsp10.txt"},
-    {0,0,0, ""}
-};
-
-/* ptr to each codebook */
-
-static float *plsp_cb[MAX_CB];
 
 /*---------------------------------------------------------------------------*\
                                                                              
@@ -98,9 +56,10 @@ float speech_to_uq_lsps(float lsp[], float ak[], float Sn[], float w[],
 \*---------------------------------------------------------------------------*/
 
 int lsp_bits(int i) {
-    return lsp_q[i].log2m;
+    return lsp_cb[i].log2m;
 }
 
+#if VECTOR_QUANTISATION
 /*---------------------------------------------------------------------------*\
                                                                              
   quantise_uniform
@@ -161,61 +120,7 @@ void lsp_quantise(
     for(i=1; i<order; i++)
        lsp_[i] = lsp_[i-1] + dlsp_[i];
 }
-
-/*---------------------------------------------------------------------------*\
-
-  scan_line()
-
-  This function reads a vector of floats from a line in a text file.
-
-\*---------------------------------------------------------------------------*/
-
-void scan_line(FILE *fp, float f[], int n)
-/*  FILE   *fp;                file ptr to text file           */
-/*  float  f[];        array of floats to return       */
-/*  int    n;          number of floats in line        */
-{
-    char   s[MAX_STR];
-    char   *ps,*pe;
-    int           i;
-
-    fgets(s,MAX_STR,fp);
-    ps = pe = s;
-    for(i=0; i<n; i++) {
-       while( isspace(*pe)) pe++;
-       while( !isspace(*pe)) pe++;
-       sscanf(ps,"%f",&f[i]);
-       ps = pe;
-    }
-}
-
-/*---------------------------------------------------------------------------*\
-
-  load_cb
-
-  Loads a single codebook (LSP vector quantiser) into memory.
-
-\*---------------------------------------------------------------------------*/
-
-void load_cb(char *filename, float *cb, int k, int m)
-{
-    FILE *ftext;
-    int   lines;
-    int   i;
-
-    ftext = fopen(filename,"rt");
-    if (ftext == NULL) {
-       printf("Error opening text file: %s\n",filename);
-       exit(1);
-    }
-
-    lines = 0;
-    for(i=0; i<m; i++) {
-       scan_line(ftext, &cb[k*lines++], k);
-    }
-
-    fclose(ftext);
-}
+#endif
 
 /*---------------------------------------------------------------------------*\
 
@@ -228,18 +133,6 @@ void load_cb(char *filename, float *cb, int k, int m)
 
 void quantise_init()
 {
-    int i,k,m;
-
-    i = 0;
-    while(lsp_q[i].k) {
-       k = lsp_q[i].k;
-               m = lsp_q[i].m;
-       plsp_cb[i] = (float*)malloc(sizeof(float)*k*m);
-       assert(plsp_cb[i] != NULL);
-       load_cb(lsp_q[i].fn, plsp_cb[i], k, m);
-       i++;
-       assert(i < MAX_CB);
-    }
 }
 
 /*---------------------------------------------------------------------------*\
@@ -252,7 +145,7 @@ void quantise_init()
 
 \*---------------------------------------------------------------------------*/
 
-long quantise(float cb[], float vec[], float w[], int k, int m, float *se)
+long quantise(const float * cb, float vec[], float w[], int k, int m, float *se)
 /* float   cb[][K];    current VQ codebook             */
 /* float   vec[];      vector to quantise              */
 /* float   w[];         weighting vector                */
@@ -283,19 +176,6 @@ long quantise(float cb[], float vec[], float w[], int k, int m, float *se)
    return(besti);
 }
 
-static float gmin=PI;
-
-float get_gmin(void) { return gmin; }
-
-void min_lsp_dist(float lsp[], int order)
-{
-    int   i;
-
-    for(i=1; i<order; i++)
-       if ((lsp[i]-lsp[i-1]) < gmin)
-           gmin = lsp[i]-lsp[i-1];
-}
-
 void check_lsp_order(float lsp[], int lpc_order)
 {
     int   i;
@@ -352,7 +232,7 @@ float lpc_model_amplitudes(
   int   index;
   float se;
   int   k,m;
-  float *cb;
+  const float * cb;
   float wt[LPC_MAX];
 
   for(i=0; i<M; i++)
@@ -381,9 +261,9 @@ float lpc_model_amplitudes(
     /* simple uniform scalar quantisers */
 
     for(i=0; i<10; i++) {
-       k = lsp_q[i].k;
-       m = lsp_q[i].m;
-       cb = plsp_cb[i];
+       k = lsp_cb[i].k;
+       m = lsp_cb[i].m;
+       cb = lsp_cb[i].cb;
        index = quantise(cb, &lsp_hz[i], wt, k, m, &se);
        lsp_hz[i] = cb[index*k];
     }
@@ -423,10 +303,14 @@ float lpc_model_amplitudes(
        lsp_[j] = lsp[j];
 
     lsp_to_lpc(lsp_, ak, order);
+#ifdef DUMP
     dump_lsp(lsp);
+#endif
   }
 
+#ifdef DUMP
   dump_E(E);
+#endif
   #ifdef SIM_QUANT
   /* simulated LPC energy quantisation */
   {
@@ -485,8 +369,10 @@ void aks_to_M2(
 
   for(i=0; i<FFT_DEC/2; i++)
     Pw[i].real = E/(Pw[i].real*Pw[i].real + Pw[i].imag*Pw[i].imag);
+#ifdef DUMP
   if (dump) 
       dump_Pw(Pw);
+#endif
 
   /* Determine magnitudes by linear interpolation of P(w) -------------------*/
 
@@ -610,7 +496,8 @@ void encode_lsps(int indexes[], float lsp[], int order)
     int    i,k,m;
     float  wt[1];
     float  lsp_hz[LPC_MAX];
-    float *cb, se;
+    const float * cb;
+    float se;
 
     /* convert from radians to Hz so we can use human readable
        frequencies */
@@ -622,9 +509,9 @@ void encode_lsps(int indexes[], float lsp[], int order)
 
     wt[0] = 1.0;
     for(i=0; i<order; i++) {
-       k = lsp_q[i].k;
-       m = lsp_q[i].m;
-       cb = plsp_cb[i];
+       k = lsp_cb[i].k;
+       m = lsp_cb[i].m;
+       cb = lsp_cb[i].cb;
        indexes[i] = quantise(cb, &lsp_hz[i], wt, k, m, &se);
     }
 }
@@ -644,11 +531,11 @@ void decode_lsps(float lsp[], int indexes[], int order)
 {
     int    i,k;
     float  lsp_hz[LPC_MAX];
-    float *cb;
+    const float * cb;
 
     for(i=0; i<order; i++) {
-       k = lsp_q[i].k;
-       cb = plsp_cb[i];
+       k = lsp_cb[i].k;
+       cb = lsp_cb[i].cb;
        lsp_hz[i] = cb[indexes[i]*k];
     }
 
index ded7645381fe41ba6e9d9ef51f65340b9faaf9e6..eeef8469955cf7f2a17cc626e3ba4116f35cdb8f 100644 (file)
@@ -39,7 +39,6 @@ float lpc_model_amplitudes(float Sn[], float w[], MODEL *model, int order,
                           int lsp,float ak[]);
 void aks_to_M2(float ak[], int order, MODEL *model, float E, float *snr, 
               int dump);
-float get_gmin(void);
 
 int   encode_Wo(float Wo);
 float decode_Wo(int index);
index d5f436a803c6a6e17a71e91190b854abefa72435..3de0ddcf59be32c3c893ef27d5137dcf2013762b 100644 (file)
@@ -14,10 +14,12 @@ TCONTPHASE_OBJ = tcontphase.o ../src/globals.o ../src/dump.o ../src/synth.o \
 
 TINTERP_OBJ    = tinterp.o ../src/sine.o ../src/four1.o ../src/interp.o
 
-TQUANT_OBJ     = tquant.o ../src/quantise.o ../src/lpc.o ../src/lsp.o \
+TQUANT_OBJ     = tquant.o ../src/quantise.o ../src/codebook.o \
+                ../src/lpc.o ../src/lsp.o \
                  ../src/dump.o ../src/four1.o
 
-TCODEC2_OBJ    = tcodec2.o ../src/quantise.o ../src/lpc.o ../src/lsp.o \
+TCODEC2_OBJ    = tcodec2.o ../src/quantise.o  ../src/codebook.o \
+               ../src/lpc.o ../src/lsp.o ../src/pack.o \
                  ../src/dump.o ../src/four1.o ../src/codec2.o ../src/sine.o \
                ../src/nlp.o ../src/postfilter.o ../src/phase.o ../src/interp.o
 
@@ -55,4 +57,5 @@ tcodec2: $(TCODEC2_OBJ)
        $(CC) -c $(CFLAGS) $< -o $@
 
 clean :
-       rm -f *.o *~ src/*~
+       rm -f *.o *~ src/*~ \
+       extract genlsp genres tcodec2 tinterp tnlp tquant vqtrain
index 4abf69c4efa3c3a288be41200e982c4a01037d38..19fa09408e5258c0d3414ef66685434ce0e3a504 100644 (file)
@@ -86,9 +86,11 @@ char *argv[];
     COMP  W[FFT_ENC];          /* DFT of w[] */
     float pitch;
     int   i; 
-    int   dump;
     float prev_Wo;
     void  *nlp_states;
+#ifdef DUMP
+    int   dump;
+#endif
 
     if (argc < 3) {
        printf("\nusage: tnlp InputRawSpeechFile OutputPitchTextFile "
@@ -110,9 +112,13 @@ char *argv[];
       exit(1);
     }
 
+#ifdef DUMP
     dump = switch_present("--dump",argc,argv);
     if (dump) 
       dump_on(argv[dump+1]);
+#else
+#warning "Compile with -DDUMP if you expect to dump anything."
+#endif
 
     nlp_states = nlp_create();
     make_analysis_window(w,W);
@@ -129,7 +135,9 @@ char *argv[];
       for(i=0; i<N; i++)
         Sn[i+M-N] = buf[i];
       dft_speech(Sw, Sn, w);
+#ifdef DUMP
       dump_Sn(Sn); dump_Sw(Sw); 
+#endif
 
       nlp(nlp_states,Sn,N,M,PITCH_MIN,PITCH_MAX,&pitch,Sw,&prev_Wo);
       prev_Wo = TWO_PI/pitch;
@@ -139,7 +147,9 @@ char *argv[];
 
     fclose(fin);
     fclose(fout);
+#ifdef DUMP
     if (dump) dump_off();
+#endif
     nlp_destroy(nlp_states);
 
     return 0;