Working codec2 support

libs/libcodec2/unittest/lsptest.c

@@ -0,0 +1,177 @@
+/*---------------------------------------------------------------------------*\
+                                                                           
+  FILE........: lsptest.c   
+  AUTHOR......: David Rowe                                                      
+  DATE CREATED: 24/8/09                                                   
+                                                                          
+  Test Speech LPC to LSP conversion and quantisation.
+                                                                          
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2009 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2, as
+  published by the Free Software Foundation.  This program is
+  distributed in the hope that it will be useful, but WITHOUT ANY
+  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+  License for more details.
+
+  You should have received a copy of the GNU Lesser General Public License
+  along with this program; if not, write to the Free Software
+  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+*/
+
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <lpc.h>
+#include <lsp.h>
+#include <sd.h>
+
+#define N 160
+#define P 10
+
+#define LPC_FLOOR 0.0002        /* autocorrelation floor */
+#define LSP_DELTA1 0.2          /* grid spacing for LSP root searches */
+#define NDFT    256             /* DFT size for SD calculation          */
+
+/* Speex lag window */
+
+const float lag_window[11] = {
+   1.00000, 0.99716, 0.98869, 0.97474, 0.95554, 0.93140, 0.90273, 0.86998, 
+   0.83367, 0.79434, 0.75258
+};
+
+/*---------------------------------------------------------------------------*\
+                                                                            
+  find_aks_for_lsp()
+                                                          
+  This function takes a frame of samples, and determines the linear           
+  prediction coefficients for that frame of samples.  Modified version of
+  find_aks from lpc.c to include autocorrelation noise floor and lag window
+  to match Speex processing steps prior to LSP conversion.
+                                                                            
+\*---------------------------------------------------------------------------*/
+
+void find_aks_for_lsp(
+  float Sn[],	/* Nsam samples with order sample memory */
+  float a[],	/* order+1 LPCs with first coeff 1.0 */
+  int Nsam,	/* number of input speech samples */
+  int order,	/* order of the LPC analysis */
+  float *E	/* residual energy */
+)
+{
+  float Wn[N];	/* windowed frame of Nsam speech samples */
+  float R[P+1];	/* order+1 autocorrelation values of Sn[] */
+  int i;
+
+  hanning_window(Sn,Wn,Nsam);
+
+  autocorrelate(Wn,R,Nsam,order);
+  R[0] += LPC_FLOOR;
+  assert(order == 10); /* lag window only defined for order == 10 */
+  for(i=0; i<=order; i++)
+      R[i] *= lag_window[i];
+  levinson_durbin(R,a,order);
+
+  *E = 0.0;
+  for(i=0; i<=order; i++)
+    *E += a[i]*R[i];
+  if (*E < 0.0)
+    *E = 1E-12;
+}
+
+/*---------------------------------------------------------------------------*\
+                                                                            
+                                MAIN 
+                                   
+\*---------------------------------------------------------------------------*/
+
+int main(int argc, char *argv[])
+{
+  FILE *fin;            /* input speech files */
+  short buf[N];         /* buffer of 16 bit speech samples */
+  float Sn[P+N];        /* input speech samples */
+  float E;
+  float ak[P+1];        /* LP coeffs */
+  float ak_[P+1];       /* quantised LP coeffs */
+  float lsp[P];
+  float lsp_[P];        /* quantised LSPs */
+  int   roots;          /* number of LSP roots found */
+  int frames;           /* frames processed so far */
+  int i;                /* loop variables */
+
+  SpeexBits bits;
+
+  float  sd;            /* SD for this frame                */
+  float  totsd;         /* accumulated SD so far            */
+  int    gt2,gt4;       /* number of frames > 2 and 4 dB SD */
+  int    unstables;     /* number of unstable LSP frames    */
+
+  if (argc < 2) {
+    printf("usage: %s InputFile\n", argv[0]);
+    exit(0);
+  }
+
+  /* Open files */
+
+  if ((fin = fopen(argv[1],"rb")) == NULL) {
+    printf("Error opening input file: %s\n",argv[1]);
+    exit(0);
+  }
+
+  /* Initialise */
+
+  frames = 0;
+  for(i=0; i<P; i++) {
+    Sn[i] = 0.0;
+  }
+  ak_[0] = 1.0;
+
+  speex_bits_init(&bits);
+
+  totsd = 0.0;
+  unstables = 0;
+  gt2 = 0; gt4 = 0;
+
+  /* Main loop */
+
+  while( (fread(buf,sizeof(short),N,fin)) == N) {
+    frames++;
+    for(i=0; i<N; i++)
+      Sn[P+i] = (float)buf[i];
+
+    /* convert to LSP domain and back */
+
+    find_aks(&Sn[P], ak, N, P, &E);
+    roots = lpc_to_lsp(&ak[1], P , lsp, 10, LSP_DELTA1, NULL);
+    if (roots == P) {
+
+        speex_bits_reset(&bits);
+	lsp_quant_lbr(lsp, lsp_, P, &bits);	
+	lsp_to_lpc(lsp_, &ak_[1], P, NULL);
+	
+	/* measure spectral distortion */
+	sd = spectral_dist(ak, ak_, P, NDFT);
+	if (sd > 2.0) gt2++;
+	if (sd > 4.0) gt4++;
+	totsd += sd;
+    }
+    else
+	unstables++;
+  }
+
+  fclose(fin);
+
+  printf("frames = %d Av sd = %3.2f dB", frames, totsd/frames);
+  printf("  >2 dB %3.2f%%  >4 dB %3.2f%%  unstables: %d\n",gt2*100.0/frames,
+         gt4*100.0/frames, unstables);
+
+  return 0;
+}
+