\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
\r
\*---------------------------------------------------------------------------*/\r
\r
-int lpc_to_lsp (float *a,int lpcrdr,float *freq,int nb,float delta)\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 the x domain */\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
\r
if((Q = (float *) malloc((m+1)*sizeof(float))) == NULL){\r
printf("not enough memory to allocate buffer\n");\r
- exit(1);\r
+ exit(1);\r
}\r
\r
if((P = (float *) malloc((m+1)*sizeof(float))) == NULL){\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
\r
-void lsp_to_lpc(float *freq,float *ak,int lpcrdr)\r
-/* float *freq array of LSP frequencies in the x domain */\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
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
#ifndef __LSP__
#define __LSP__
-int lpc_to_lsp (float *a,int lpcrdr,float *freq,int nb,float delta);
-void lsp_to_lpc(float *freq,float *ak,int lpcrdr);
+int lpc_to_lsp (float *a, int lpcrdr, float *freq, int nb, float delta);
+void lsp_to_lpc(float *freq, float *ak, int lpcrdr);
#endif
#include "sine.h"
#include "quantise.h"
#include "lpc.h"
+#include "lsp.h"
#include "dump.h"
#define MAX_ORDER 20
float lsp[MAX_ORDER];
float lsp_hz[MAX_ORDER];
float lsp_[MAX_ORDER];
- float lspd[MAX_ORDER];
int roots; /* number of LSP roots found */
int index;
float se;
- int l,k,m;
+ int k,m;
float *cb;
float wt[MAX_ORDER];
- float maxA, dB;
-
for(i=0; i<M; i++)
Wn[i] = Sn[i]*w[i];
autocorrelate(Wn,R,M,order);
for(i=0; i<=order; i++)
E += ak[i]*R[i];
+ for(i=0; i<order; i++)
+ wt[i] = 1.0;
+
if (lsp_quant) {
- roots = lpc_to_lsp(&ak[1], order, lsp, 5, LSP_DELTA1, NULL);
+ roots = lpc_to_lsp(ak, order, lsp, 5, LSP_DELTA1);
if (roots != order)
printf("LSP roots not found\n");
+ /* convert from radians to Hz to make quantisers more
+ human readable */
+
for(i=0; i<order; i++)
lsp_hz[i] = (4000.0/PI)*lsp[i];
+ /* simple uniform scalar quantisers */
+
for(i=0; i<10; i++) {
k = lsp_q[i].k;
m = lsp_q[i].m;
lsp_hz[i] = cb[index*k];
}
- /*
+ /* experiment: simulating uniform quantisation error
for(i=0; i<order; i++)
lsp[i] += PI*(12.5/4000.0)*(1.0 - 2.0*(float)rand()/RAND_MAX);
*/
for(i=0; i<order; i++)
lsp[i] = (PI/4000.0)*lsp_hz[i];
+ /* Bandwidth Expansion (BW). Prevents any two LSPs getting too
+ close together after quantisation. We know from experiment
+ that LSP quantisation errors < 12.5Hz (25Hz setp size) are
+ inaudible so we use that as the minimum LSP separation.
+ */
+
for(i=1; i<5; i++) {
if (lsp[i] - lsp[i-1] < PI*(12.5/4000.0))
lsp[i] = lsp[i-1] + PI*(12.5/4000.0);
/* as quantiser gaps increased, larger BW expansion was required
to prevent twinkly noises */
+
for(i=5; i<8; i++) {
if (lsp[i] - lsp[i-1] < PI*(25.0/4000.0))
lsp[i] = lsp[i-1] + PI*(25.0/4000.0);
lsp[i] = lsp[i-1] + PI*(75.0/4000.0);
}
- //#define OLD_VQ
-#ifdef OLD_VQ
- lspd[0] = lsp[0];
- for(i=1; i<order; i++)
- lspd[i] = lsp[i] - lsp[i-1];
- for(i=0; i<order; i++)
- wt[i] = 1.0;
-
- i = 0; /* i-th codebook */
- l = 0; /* which starts at l-th lsp */
- while(lsp_q[i].k) {
- k = lsp_q[i].k;
- m = lsp_q[i].m;
- cb = plsp_cb[i];
- index = quantise(cb, &lspd[l], wt, k, m, &se);
-
- for(j=0; j<k; j++)
- lspd[l+j] = cb[index*k+j];
-
- /* compute quantised lsp so we can adjust for quantisation error
- below */
-
- for(j=l; j<l+k; j++) {
- if (j==0)
- lsp_[0] = lspd[0];
- else
- lsp_[j] = lsp_[j-1] + lspd[j];
- }
-
- l += k;
- assert(l <= order);
-
- /* adjust next lspd to account for quantisation error */
-
- lspd[l] = lsp[l] - lsp_[l-1];
-
- i++;
- assert(i < MAX_CB);
- }
-#else
- l = 0;
-#endif
- /* used during development: copy remaining LSPs from orig if we haven't
- quantised all of them */
- for(j=l; j<order; j++)
+ for(j=0; j<order; j++)
lsp_[j] = lsp[j];
- lsp_to_lpc(lsp_, &ak[1], order, NULL);
+ lsp_to_lpc(lsp_, ak, order);
dump_lsp(lsp);
}
aks_to_M2(ak,order,model,E,&snr); /* {ak} -> {Am} LPC decode */
- #ifdef CLICKY
- /* Adding a random component to low energy harmonic phase seems to
- improve low pitch speakers. Adding a small random component to
- low energy harmonic amplitudes also helps low pitch speakers after
- LPC modelling (see LPC modelling/amplitude quantisation code).
- */
-
- maxA = 0.0;
- for(i=1; i<=model->L; i++) {
- if (model->A[i] > maxA) {
- maxA = model->A[i];
- }
- }
- for(i=1; i<=model->L; i++) {
- if (model->A[i] < 0.1*maxA) {
- dB = 3.0 - 6.0*(float)rand()/RAND_MAX;
- model->A[i] *= pow(10.0, dB/20.0);
- }
- }
- #endif
-
return snr;
}
projects / freetel-svn-tracking.git / commitdiff