=Norig) twidx-=Norig;
+ C_MUL(t,scratch[q] , twiddles[twidx] );
+ C_ADDTO( Fout[ k ] ,t);
+ }
+ k += m;
+ }
+ }
+ KISS_FFT_TMP_FREE(scratch);
+}
+
+static
+void kf_work(
+ kiss_fft_cpx * Fout,
+ const kiss_fft_cpx * f,
+ const size_t fstride,
+ int in_stride,
+ int * factors,
+ const kiss_fft_cfg st
+ )
+{
+ kiss_fft_cpx * Fout_beg=Fout;
+ const int p=*factors++; /* the radix */
+ const int m=*factors++; /* stage's fft length/p */
+ const kiss_fft_cpx * Fout_end = Fout + p*m;
+
+#ifdef _OPENMP
+ // use openmp extensions at the
+ // top-level (not recursive)
+ if (fstride==1 && p<=5)
+ {
+ int k;
+
+ // execute the p different work units in different threads
+# pragma omp parallel for
+ for (k=0;k floor_sqrt)
+ p = n; /* no more factors, skip to end */
+ }
+ n /= p;
+ *facbuf++ = p;
+ *facbuf++ = n;
+ } while (n > 1);
+}
+
+/*
+ *
+ * User-callable function to allocate all necessary storage space for the fft.
+ *
+ * The return value is a contiguous block of memory, allocated with malloc. As such,
+ * It can be freed with free(), rather than a kiss_fft-specific function.
+ * */
+kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem )
+{
+ kiss_fft_cfg st=NULL;
+ size_t memneeded = sizeof(struct kiss_fft_state)
+ + sizeof(kiss_fft_cpx)*(nfft-1); /* twiddle factors*/
+
+ if ( lenmem==NULL ) {
+ st = ( kiss_fft_cfg)KISS_FFT_MALLOC( memneeded );
+ }else{
+ if (mem != NULL && *lenmem >= memneeded)
+ st = (kiss_fft_cfg)mem;
+ *lenmem = memneeded;
+ }
+ if (st) {
+ int i;
+ st->nfft=nfft;
+ st->inverse = inverse_fft;
+
+ for (i=0;iinverse)
+ phase *= -1;
+ kf_cexp(st->twiddles+i, phase );
+ }
+
+ kf_factor(nfft,st->factors);
+ }
+ return st;
+}
+
+
+void kiss_fft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride)
+{
+ if (fin == fout) {
+ //NOTE: this is not really an in-place FFT algorithm.
+ //It just performs an out-of-place FFT into a temp buffer
+ kiss_fft_cpx * tmpbuf = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC( sizeof(kiss_fft_cpx)*st->nfft);
+ kf_work(tmpbuf,fin,1,in_stride, st->factors,st);
+ memcpy(fout,tmpbuf,sizeof(kiss_fft_cpx)*st->nfft);
+ KISS_FFT_TMP_FREE(tmpbuf);
+ }else{
+ kf_work( fout, fin, 1,in_stride, st->factors,st );
+ }
+}
+
+void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout)
+{
+ kiss_fft_stride(cfg,fin,fout,1);
+}
+
+
+void kiss_fft_cleanup(void)
+{
+ // nothing needed any more
+}
+
+int kiss_fft_next_fast_size(int n)
+{
+ while(1) {
+ int m=n;
+ while ( (m%2) == 0 ) m/=2;
+ while ( (m%3) == 0 ) m/=3;
+ while ( (m%5) == 0 ) m/=5;
+ if (m<=1)
+ break; /* n is completely factorable by twos, threes, and fives */
+ n++;
+ }
+ return n;
+}
diff --git a/libs/libcodec2/src/kiss_fft.h b/libs/libcodec2/src/kiss_fft.h
new file mode 100644
index 0000000000..64c50f4aae
--- /dev/null
+++ b/libs/libcodec2/src/kiss_fft.h
@@ -0,0 +1,124 @@
+#ifndef KISS_FFT_H
+#define KISS_FFT_H
+
+#include
+#include
+#include
+#include
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/*
+ ATTENTION!
+ If you would like a :
+ -- a utility that will handle the caching of fft objects
+ -- real-only (no imaginary time component ) FFT
+ -- a multi-dimensional FFT
+ -- a command-line utility to perform ffts
+ -- a command-line utility to perform fast-convolution filtering
+
+ Then see kfc.h kiss_fftr.h kiss_fftnd.h fftutil.c kiss_fastfir.c
+ in the tools/ directory.
+*/
+
+#ifdef USE_SIMD
+# include
+# define kiss_fft_scalar __m128
+#define KISS_FFT_MALLOC(nbytes) _mm_malloc(nbytes,16)
+#define KISS_FFT_FREE _mm_free
+#else
+#define KISS_FFT_MALLOC malloc
+#define KISS_FFT_FREE free
+#endif
+
+
+#ifdef FIXED_POINT
+#include
+# if (FIXED_POINT == 32)
+# define kiss_fft_scalar int32_t
+# else
+# define kiss_fft_scalar int16_t
+# endif
+#else
+# ifndef kiss_fft_scalar
+/* default is float */
+# define kiss_fft_scalar float
+# endif
+#endif
+
+typedef struct {
+ kiss_fft_scalar r;
+ kiss_fft_scalar i;
+}kiss_fft_cpx;
+
+typedef struct kiss_fft_state* kiss_fft_cfg;
+
+/*
+ * kiss_fft_alloc
+ *
+ * Initialize a FFT (or IFFT) algorithm's cfg/state buffer.
+ *
+ * typical usage: kiss_fft_cfg mycfg=kiss_fft_alloc(1024,0,NULL,NULL);
+ *
+ * The return value from fft_alloc is a cfg buffer used internally
+ * by the fft routine or NULL.
+ *
+ * If lenmem is NULL, then kiss_fft_alloc will allocate a cfg buffer using malloc.
+ * The returned value should be free()d when done to avoid memory leaks.
+ *
+ * The state can be placed in a user supplied buffer 'mem':
+ * If lenmem is not NULL and mem is not NULL and *lenmem is large enough,
+ * then the function places the cfg in mem and the size used in *lenmem
+ * and returns mem.
+ *
+ * If lenmem is not NULL and ( mem is NULL or *lenmem is not large enough),
+ * then the function returns NULL and places the minimum cfg
+ * buffer size in *lenmem.
+ * */
+
+kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem);
+
+/*
+ * kiss_fft(cfg,in_out_buf)
+ *
+ * Perform an FFT on a complex input buffer.
+ * for a forward FFT,
+ * fin should be f[0] , f[1] , ... ,f[nfft-1]
+ * fout will be F[0] , F[1] , ... ,F[nfft-1]
+ * Note that each element is complex and can be accessed like
+ f[k].r and f[k].i
+ * */
+void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout);
+
+/*
+ A more generic version of the above function. It reads its input from every Nth sample.
+ * */
+void kiss_fft_stride(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int fin_stride);
+
+/* If kiss_fft_alloc allocated a buffer, it is one contiguous
+ buffer and can be simply free()d when no longer needed*/
+#define kiss_fft_free free
+
+/*
+ Cleans up some memory that gets managed internally. Not necessary to call, but it might clean up
+ your compiler output to call this before you exit.
+*/
+void kiss_fft_cleanup(void);
+
+
+/*
+ * Returns the smallest integer k, such that k>=n and k has only "fast" factors (2,3,5)
+ */
+int kiss_fft_next_fast_size(int n);
+
+/* for real ffts, we need an even size */
+#define kiss_fftr_next_fast_size_real(n) \
+ (kiss_fft_next_fast_size( ((n)+1)>>1)<<1)
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
diff --git a/libs/libcodec2/src/listensim.sh b/libs/libcodec2/src/listensim.sh
index 64f7455ab3..b296cac588 100755
--- a/libs/libcodec2/src/listensim.sh
+++ b/libs/libcodec2/src/listensim.sh
@@ -4,6 +4,6 @@
#
# Listen to files processed with sim.sh
-../script/menu.sh ../raw/$1.raw $1_uq.raw $1_phase0.raw $1_lpc10.raw $1_lsp.raw $1_phase0_lpc10.raw $1_phase0_lsp.raw $1_phase0_lsp_dec.raw $2 $3
+../script/menu.sh $1_uq.raw $1_lpc10.raw $1_lpcpf.raw $1_phase0.raw $1_phase0_lpcpf.raw $2 $3 $4 $5
diff --git a/libs/libcodec2/src/lpc.c b/libs/libcodec2/src/lpc.c
index 1f9ff2bf10..a253289a46 100644
--- a/libs/libcodec2/src/lpc.c
+++ b/libs/libcodec2/src/lpc.c
@@ -2,14 +2,14 @@
FILE........: lpc.c
AUTHOR......: David Rowe
- DATE CREATED: 30/9/90
+ DATE CREATED: 30 Sep 1990 (!)
Linear Prediction functions written in C.
\*---------------------------------------------------------------------------*/
/*
- Copyright (C) 2009 David Rowe
+ Copyright (C) 2009-2012 David Rowe
All rights reserved.
@@ -22,18 +22,74 @@
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.
+ along with this program; if not, see .
*/
#define LPC_MAX_N 512 /* maximum no. of samples in frame */
#define PI 3.141592654 /* mathematical constant */
+#define ALPHA 1.0
+#define BETA 0.94
+
#include
#include
#include "defines.h"
#include "lpc.h"
+/*---------------------------------------------------------------------------*\
+
+ pre_emp()
+
+ Pre-emphasise (high pass filter with zero close to 0 Hz) a frame of
+ speech samples. Helps reduce dynamic range of LPC spectrum, giving
+ greater weight and hensea better match to low energy formants.
+
+ Should be balanced by de-emphasis of the output speech.
+
+\*---------------------------------------------------------------------------*/
+
+void pre_emp(
+ float Sn_pre[], /* output frame of speech samples */
+ float Sn[], /* input frame of speech samples */
+ float *mem, /* Sn[-1]single sample memory */
+ int Nsam /* number of speech samples to use */
+)
+{
+ int i;
+
+ for(i=0; i.
*/
#ifndef __LPC__
@@ -31,6 +30,8 @@
#define LPC_MAX_ORDER 20
+void pre_emp(float Sn_pre[], float Sn[], float *mem, int Nsam);
+void de_emp(float Sn_se[], float Sn[], float *mem, int Nsam);
void hanning_window(float Sn[], float Wn[], int Nsam);
void autocorrelate(float Sn[], float Rn[], int Nsam, int order);
void levinson_durbin(float R[], float lpcs[], int order);
diff --git a/libs/libcodec2/src/lsp.c b/libs/libcodec2/src/lsp.c
index feab4219ab..47001c1efd 100644
--- a/libs/libcodec2/src/lsp.c
+++ b/libs/libcodec2/src/lsp.c
@@ -1,323 +1,325 @@
-/*---------------------------------------------------------------------------*\
-
- 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
-#include
-#include
-
-/*---------------------------------------------------------------------------*\
-
- 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.
-
-\*---------------------------------------------------------------------------*/
-
-
-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,*t,*u,*v,sum;
-
- /* Allocate memory for chebyshev series formulation */
-
- if((T = (float *)malloc((m/2+1)*sizeof(float))) == NULL){
- fprintf(stderr, "not enough memory to allocate buffer\n");
- exit(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++;
-
- free(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;
- float temp_psumr;
- int i,j,m,flag,k;
- float *Q; /* ptrs for memory allocation */
- float *P;
- 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 */
- flag = 1;
- m = lpcrdr/2; /* order of P'(z) & Q'(z) polynimials */
-
- /* Allocate memory space for polynomials */
-
- Q = (float *) malloc((m+1)*sizeof(float));
- P = (float *) malloc((m+1)*sizeof(float));
- if( (P == NULL) || (Q == NULL) ) {
- fprintf(stderr,"not enough memory to allocate buffer\n");
- exit(1);
- }
-
- /* 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= -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;
- }
- }
- }
- free(P); /* free memory space */
- free(Q);
-
- /* convert from x domain to radians */
-
- for(i=0; i.
+*/
+
+#include "defines.h"
+#include "lsp.h"
+#include
+#include
+#include
+
+/* 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= -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.
+*/
+
#ifndef __LSP__
#define __LSP__
diff --git a/libs/libcodec2/src/nlp.c b/libs/libcodec2/src/nlp.c
index 193ca92109..3214578e84 100644
--- a/libs/libcodec2/src/nlp.c
+++ b/libs/libcodec2/src/nlp.c
@@ -4,8 +4,8 @@
AUTHOR......: David Rowe
DATE CREATED: 23/3/93
- Non Linear Pitch (NLP) estimation functions.
-
+ Non Linear Pitch (NLP) estimation functions.
+
\*---------------------------------------------------------------------------*/
/*
@@ -22,14 +22,13 @@
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.
+ along with this program; if not, see .
*/
#include "defines.h"
#include "nlp.h"
#include "dump.h"
-#include "four1.h"
+#include "kiss_fft.h"
#include
#include
@@ -60,7 +59,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,
@@ -112,12 +111,14 @@ float nlp_fir[] = {
};
typedef struct {
- float sq[PMAX_M]; /* squared speech samples */
- float mem_x,mem_y; /* memory for notch filter */
- float mem_fir[NLP_NTAP]; /* decimation FIR filter memory */
+ float sq[PMAX_M]; /* squared speech samples */
+ float mem_x,mem_y; /* memory for notch filter */
+ float mem_fir[NLP_NTAP]; /* decimation FIR filter memory */
+ kiss_fft_cfg fft_cfg; /* kiss FFT config */
} NLP;
-float post_process_mbe(COMP Fw[], int pmin, int pmax, float gmax);
+float test_candidate_mbe(COMP Sw[], COMP W[], float f0);
+float post_process_mbe(COMP Fw[], int pmin, int pmax, float gmax, COMP Sw[], COMP W[], float *prev_Wo);
float post_process_sub_multiples(COMP Fw[],
int pmin, int pmax, float gmax, int gmax_bin,
float *prev_Wo);
@@ -146,20 +147,27 @@ void *nlp_create()
for(i=0; imem_fir[i] = 0.0;
+ nlp->fft_cfg = kiss_fft_alloc (PE_FFT_SIZE, 0, NULL, NULL);
+ assert(nlp->fft_cfg != NULL);
+
return (void*)nlp;
}
/*---------------------------------------------------------------------------*\
- nlp_destory()
+ nlp_destroy()
- Initialisation function for NLP pitch estimator.
+ Shut down function for NLP pitch estimator.
\*---------------------------------------------------------------------------*/
void nlp_destroy(void *nlp_state)
{
+ NLP *nlp;
assert(nlp_state != NULL);
+ nlp = (NLP*)nlp_state;
+
+ KISS_FFT_FREE(nlp->fft_cfg);
free(nlp_state);
}
@@ -198,27 +206,30 @@ float nlp(
float Sn[], /* input speech vector */
int n, /* frames shift (no. new samples in Sn[]) */
int m, /* analysis window size */
- int pmin, /* minimum pitch value */
+ int pmin, /* minimum pitch value */
int pmax, /* maximum pitch value */
float *pitch, /* estimated pitch period in samples */
COMP Sw[], /* Freq domain version of Sn[] */
+ COMP W[], /* Freq domain window */
float *prev_Wo
)
{
NLP *nlp;
- float notch; /* current notch filter output */
- COMP Fw[PE_FFT_SIZE]; /* DFT of squared signal */
+ float notch; /* current notch filter output */
+ COMP fw[PE_FFT_SIZE]; /* DFT of squared signal (input) */
+ COMP Fw[PE_FFT_SIZE]; /* DFT of squared signal (output) */
float gmax;
int gmax_bin;
int i,j;
float best_f0;
assert(nlp_state != NULL);
+ assert(m <= PMAX_M);
nlp = (NLP*)nlp_state;
/* Square, notch filter at DC, and LP filter vector */
- for(i=m-n; isq[i] = Sn[i]*Sn[i];
for(i=m-n; imem_y;
nlp->mem_x = nlp->sq[i];
nlp->mem_y = notch;
- nlp->sq[i] = notch;
+ nlp->sq[i] = notch + 1.0; /* With 0 input vectors to codec,
+ kiss_fft() would take a long
+ time to execute when running in
+ real time. Problem was traced
+ to kiss_fft function call in
+ this function. Adding this small
+ constant fixed problem. Not
+ exactly sure why. */
}
for(i=m-n; isq[i*DEC]*(0.5 - 0.5*cos(2*PI*i/(m/DEC-1)));
+ fw[i].real = nlp->sq[i*DEC]*(0.5 - 0.5*cos(2*PI*i/(m/DEC-1)));
}
+ #ifdef DUMP
dump_dec(Fw);
- four1(&Fw[-1].imag,PE_FFT_SIZE,1);
+ #endif
+
+ kiss_fft(nlp->fft_cfg, (kiss_fft_cpx *)fw, (kiss_fft_cpx *)Fw);
for(i=0; isq);
dump_Fw(Fw);
+ #endif
/* find global peak */
@@ -267,9 +290,13 @@ float nlp(
gmax_bin = i;
}
}
-
- best_f0 = post_process_sub_multiples(Fw, pmin, pmax, gmax, gmax_bin,
- prev_Wo);
+
+ //#define POST_PROCESS_MBE
+ #ifdef POST_PROCESS_MBE
+ best_f0 = post_process_mbe(Fw, pmin, pmax, gmax, Sw, W, prev_Wo);
+ #else
+ best_f0 = post_process_sub_multiples(Fw, pmin, pmax, gmax, gmax_bin, prev_Wo);
+ #endif
/* Shift samples in buffer to make room for new samples */
@@ -286,7 +313,7 @@ float nlp(
post_process_sub_multiples()
- Given the global maximma of Fw[] we search interger submultiples for
+ Given the global maximma of Fw[] we search integer submultiples for
local maxima. If local maxima exist and they are above an
experimentally derived threshold (OK a magic number I pulled out of
the air) we choose the submultiple as the F0 estimate.
@@ -317,10 +344,10 @@ float post_process_sub_multiples(COMP Fw[],
/* post process estimate by searching submultiples */
mult = 2;
- min_bin = PE_FFT_SIZE*DEC/pmax;
+ min_bin = PE_FFT_SIZE*DEC/pmax;
cmax_bin = gmax_bin;
prev_f0_bin = *prev_Wo*(4000.0/PI)*(PE_FFT_SIZE*DEC)/SAMPLE_RATE;
-
+
while(gmax_bin/mult >= min_bin) {
b = gmax_bin/mult; /* determine search interval */
@@ -339,7 +366,7 @@ float post_process_sub_multiples(COMP Fw[],
lmax = 0;
lmax_bin = bmin;
- for (b=bmin; b<=bmax; b++) /* look for maximum in interval */
+ for (b=bmin; b<=bmax; b++) /* look for maximum in interval */
if (Fw[b].real > lmax) {
lmax = Fw[b].real;
lmax_bin = b;
@@ -359,3 +386,158 @@ float post_process_sub_multiples(COMP Fw[],
return best_f0;
}
+/*---------------------------------------------------------------------------*\
+
+ post_process_mbe()
+
+ Use the MBE pitch estimation algorithm to evaluate pitch candidates. This
+ works OK but the accuracy at low F0 is affected by NW, the analysis window
+ size used for the DFT of the input speech Sw[]. Also favours high F0 in
+ the presence of background noise which causes periodic artifacts in the
+ synthesised speech.
+
+\*---------------------------------------------------------------------------*/
+
+float post_process_mbe(COMP Fw[], int pmin, int pmax, float gmax, COMP Sw[], COMP W[], float *prev_Wo)
+{
+ float candidate_f0;
+ float f0,best_f0; /* fundamental frequency */
+ float e,e_min; /* MBE cost function */
+ int i;
+ float e_hz[F0_MAX];
+ int bin;
+ float f0_min, f0_max;
+ float f0_start, f0_end;
+
+ f0_min = (float)SAMPLE_RATE/pmax;
+ f0_max = (float)SAMPLE_RATE/pmin;
+
+ /* Now look for local maxima. Each local maxima is a candidate
+ that we test using the MBE pitch estimation algotithm */
+
+ for(i=0; i Fw[i-1].real) && (Fw[i].real > Fw[i+1].real)) {
+
+ /* local maxima found, lets test if it's big enough */
+
+ if (Fw[i].real > T*gmax) {
+
+ /* OK, sample MBE cost function over +/- 10Hz range in 2.5Hz steps */
+
+ candidate_f0 = (float)i*SAMPLE_RATE/(PE_FFT_SIZE*DEC);
+ f0_start = candidate_f0-20;
+ f0_end = candidate_f0+20;
+ if (f0_start < f0_min) f0_start = f0_min;
+ if (f0_end > f0_max) f0_end = f0_max;
+
+ for(f0=f0_start; f0<=f0_end; f0+= 2.5) {
+ e = test_candidate_mbe(Sw, W, f0);
+ bin = floor(f0); assert((bin > 0) && (bin < F0_MAX));
+ e_hz[bin] = e;
+ if (e < e_min) {
+ e_min = e;
+ best_f0 = f0;
+ }
+ }
+
+ }
+ }
+ }
+
+ /* finally sample MBE cost function around previous pitch estimate
+ (form of pitch tracking) */
+
+ candidate_f0 = *prev_Wo * SAMPLE_RATE/TWO_PI;
+ f0_start = candidate_f0-20;
+ f0_end = candidate_f0+20;
+ if (f0_start < f0_min) f0_start = f0_min;
+ if (f0_end > f0_max) f0_end = f0_max;
+
+ for(f0=f0_start; f0<=f0_end; f0+= 2.5) {
+ e = test_candidate_mbe(Sw, W, f0);
+ bin = floor(f0); assert((bin > 0) && (bin < F0_MAX));
+ e_hz[bin] = e;
+ if (e < e_min) {
+ e_min = e;
+ best_f0 = f0;
+ }
+ }
+
+ #ifdef DUMP
+ dump_e(e_hz);
+ #endif
+
+ return best_f0;
+}
+
+/*---------------------------------------------------------------------------*\
+
+ test_candidate_mbe()
+
+ Returns the error of the MBE cost function for the input f0.
+
+ Note: I think a lot of the operations below can be simplified as
+ W[].imag = 0 and has been normalised such that den always equals 1.
+
+\*---------------------------------------------------------------------------*/
+
+float test_candidate_mbe(
+ COMP Sw[],
+ COMP W[],
+ float f0
+)
+{
+ COMP Sw_[FFT_ENC]; /* DFT of all voiced synthesised signal */
+ int l,al,bl,m; /* loop variables */
+ COMP Am; /* amplitude sample for this band */
+ int offset; /* centers Hw[] about current harmonic */
+ float den; /* denominator of Am expression */
+ float error; /* accumulated error between originl and synthesised */
+ float Wo; /* current "test" fundamental freq. */
+ int L;
+
+ L = floor((SAMPLE_RATE/2.0)/f0);
+ Wo = f0*(2*PI/SAMPLE_RATE);
+
+ error = 0.0;
+
+ /* Just test across the harmonics in the first 1000 Hz (L/4) */
+
+ for(l=1; l.
*/
#ifndef __NLP__
#define __NLP__
+#include "comp.h"
+
void *nlp_create();
void nlp_destroy(void *nlp_state);
float nlp(void *nlp_state, float Sn[], int n, int m, int pmin, int pmax,
- float *pitch, COMP Sw[], float *prev_Wo);
-float test_candidate_mbe(COMP Sw[], float f0, COMP Sw_[]);
+ float *pitch, COMP Sw[], COMP W[], float *prev_Wo);
#endif
diff --git a/libs/libcodec2/src/octave.c b/libs/libcodec2/src/octave.c
new file mode 100644
index 0000000000..2ff5ad1413
--- /dev/null
+++ b/libs/libcodec2/src/octave.c
@@ -0,0 +1,85 @@
+/*---------------------------------------------------------------------------*\
+
+ FILE........: octave.c
+ AUTHOR......: David Rowe
+ DATE CREATED: April 28 2012
+
+ Functions to save C arrays in GNU Octave matrix format. The output text
+ file can be directly read into Octave using "load filename".
+
+\*---------------------------------------------------------------------------*/
+
+
+/*
+ Copyright (C) 2012 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, see .
+*/
+
+#include
+#include "octave.h"
+
+void octave_save_int(FILE *f, char name[], int data[], int rows, int cols)
+{
+ int r,c;
+
+ fprintf(f, "# name: %s\n", name);
+ fprintf(f, "# type: matrix\n");
+ fprintf(f, "# rows: %d\n", rows);
+ fprintf(f, "# columns: %d\n", cols);
+
+ for(r=0; r.
+*/
+
+#ifndef __OCTAVE__
+#define __OCTAVE__
+
+#include "comp.h"
+
+void octave_save_int(FILE *f, char name[], int data[], int rows, int cols);
+void octave_save_float(FILE *f, char name[], float data[], int rows, int cols, int col_len);
+void octave_save_complex(FILE *f, char name[], COMP data[], int rows, int cols, int col_len);
+
+#endif
diff --git a/libs/libcodec2/src/os.h b/libs/libcodec2/src/os.h
new file mode 100644
index 0000000000..0dae9bfd24
--- /dev/null
+++ b/libs/libcodec2/src/os.h
@@ -0,0 +1,53 @@
+/* Generate using fir1(47,1/6) in Octave */
+
+const float fdmdv_os_filter[]= {
+ -3.55606818e-04,
+ -8.98615286e-04,
+ -1.40119781e-03,
+ -1.71713852e-03,
+ -1.56471179e-03,
+ -6.28128960e-04,
+ 1.24522223e-03,
+ 3.83138676e-03,
+ 6.41309478e-03,
+ 7.85893186e-03,
+ 6.93514929e-03,
+ 2.79361991e-03,
+ -4.51051400e-03,
+ -1.36671853e-02,
+ -2.21034939e-02,
+ -2.64084653e-02,
+ -2.31425052e-02,
+ -9.84218694e-03,
+ 1.40648474e-02,
+ 4.67316298e-02,
+ 8.39615986e-02,
+ 1.19925275e-01,
+ 1.48381174e-01,
+ 1.64097819e-01,
+ 1.64097819e-01,
+ 1.48381174e-01,
+ 1.19925275e-01,
+ 8.39615986e-02,
+ 4.67316298e-02,
+ 1.40648474e-02,
+ -9.84218694e-03,
+ -2.31425052e-02,
+ -2.64084653e-02,
+ -2.21034939e-02,
+ -1.36671853e-02,
+ -4.51051400e-03,
+ 2.79361991e-03,
+ 6.93514929e-03,
+ 7.85893186e-03,
+ 6.41309478e-03,
+ 3.83138676e-03,
+ 1.24522223e-03,
+ -6.28128960e-04,
+ -1.56471179e-03,
+ -1.71713852e-03,
+ -1.40119781e-03,
+ -8.98615286e-04,
+ -3.55606818e-04
+};
+
diff --git a/libs/libcodec2/src/pack.c b/libs/libcodec2/src/pack.c
index 2cbff4438a..5d67c3296e 100644
--- a/libs/libcodec2/src/pack.c
+++ b/libs/libcodec2/src/pack.c
@@ -1,20 +1,20 @@
/*
Copyright (C) 2010 Perens LLC
- This program is free software: you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
+ All rights reserved.
- 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.
+ 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.1, 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 General Public License
- along with this program. If not, see .
+ You should have received a copy of the GNU Lesser General Public License
+ along with this program; if not, see .
+*/
- */
#include "defines.h"
#include "quantise.h"
#include
@@ -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);
diff --git a/libs/libcodec2/src/phase.c b/libs/libcodec2/src/phase.c
index 83fd680e79..d41303b73c 100644
--- a/libs/libcodec2/src/phase.c
+++ b/libs/libcodec2/src/phase.c
@@ -22,20 +22,22 @@
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.
+ along with this program; if not,see .
*/
#include "defines.h"
#include "phase.h"
-#include "four1.h"
+#include "kiss_fft.h"
+#include "comp.h"
+#include "glottal.c"
#include
+#include
#include
#include
#include
-#define VTHRESH 4.0
+#define GLOTTAL_FFT_SIZE 512
/*---------------------------------------------------------------------------*\
@@ -47,6 +49,7 @@
\*---------------------------------------------------------------------------*/
void aks_to_H(
+ kiss_fft_cfg fft_fwd_cfg,
MODEL *model, /* model parameters */
float aks[], /* LPC's */
float G, /* energy term */
@@ -54,7 +57,8 @@ void aks_to_H(
int order
)
{
- COMP Pw[FFT_DEC]; /* power spectrum */
+ COMP pw[FFT_ENC]; /* power spectrum (input) */
+ COMP Pw[FFT_ENC]; /* power spectrum (output) */
int i,m; /* loop variables */
int am,bm; /* limits of current band */
float r; /* no. rads/bin */
@@ -63,19 +67,19 @@ void aks_to_H(
int b; /* centre bin of harmonic */
float phi_; /* phase of LPC spectra */
- r = TWO_PI/(FFT_DEC);
+ r = TWO_PI/(FFT_ENC);
/* Determine DFT of A(exp(jw)) ------------------------------------------*/
- for(i=0; iWo)*N/2;
+ ex_phase[0] += (*prev_Wo+model->Wo)*N/2;
*/
ex_phase[0] += (model->Wo)*N;
ex_phase[0] -= TWO_PI*floor(ex_phase[0]/TWO_PI + 0.5);
+ r = TWO_PI/GLOTTAL_FFT_SIZE;
for(m=1; m<=model->L; m++) {
-
+
/* generate excitation */
-
+
if (model->voiced) {
- /* This method of adding jitter really helped remove the clicky
- sound in low pitched makes like hts1a. This moves the onset
- of each harmonic over at +/- 0.25 of a sample.
+ //float rnd;
+
+ b = floor(m*model->Wo/r + 0.5);
+ if (b > ((GLOTTAL_FFT_SIZE/2)-1)) {
+ b = (GLOTTAL_FFT_SIZE/2)-1;
+ }
+
+ /* I think adding a little jitter helps improve low pitch
+ males like hts1a. This moves the onset of each harmonic
+ over +/- 0.25 of a sample.
*/
- jitter = 0.25*(1.0 - 2.0*rand()/RAND_MAX);
- Ex[m].real = cos(ex_phase[0]*m - jitter*model->Wo*m);
- Ex[m].imag = sin(ex_phase[0]*m - jitter*model->Wo*m);
+ //jitter = 0.25*(1.0 - 2.0*rand()/RAND_MAX);
+ jitter = 0;
+
+ //rnd = (PI/8)*(1.0 - 2.0*rand()/RAND_MAX);
+ Ex[m].real = cos(ex_phase[0]*m/* - jitter*model->Wo*m + glottal[b]*/);
+ Ex[m].imag = sin(ex_phase[0]*m/* - jitter*model->Wo*m + glottal[b]*/);
}
else {
@@ -252,3 +270,4 @@ void phase_synth_zero_order(
}
}
+
diff --git a/libs/libcodec2/src/phase.h b/libs/libcodec2/src/phase.h
index 6dbf3fa2d6..367948dffb 100644
--- a/libs/libcodec2/src/phase.h
+++ b/libs/libcodec2/src/phase.h
@@ -22,13 +22,18 @@
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.
+ along with this program; if not, see .
*/
#ifndef __PHASE__
#define __PHASE__
-void phase_synth_zero_order(MODEL *model, float aks[], float *ex_phase);
+#include "kiss_fft.h"
+
+void phase_synth_zero_order(kiss_fft_cfg fft_dec_cfg,
+ MODEL *model,
+ float aks[],
+ float *ex_phase,
+ int order);
#endif
diff --git a/libs/libcodec2/src/phaseexp.c b/libs/libcodec2/src/phaseexp.c
new file mode 100644
index 0000000000..57db0f0919
--- /dev/null
+++ b/libs/libcodec2/src/phaseexp.c
@@ -0,0 +1,1574 @@
+/*---------------------------------------------------------------------------*\
+
+ FILE........: phaseexp.c
+ AUTHOR......: David Rowe
+ DATE CREATED: June 2012
+
+ Experimental functions for quantising, modelling and synthesising phase.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+ Copyright (C) 2012 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.1, 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,see .
+*/
+
+#include "defines.h"
+#include "phase.h"
+#include "kiss_fft.h"
+#include "comp.h"
+
+#include
+#include
+#include
+#include
+#include
+
+/* Bruce Perens' funcs to load codebook files */
+
+struct codebook {
+ unsigned int k;
+ unsigned int log2m;
+ unsigned int m;
+ COMP *cb;
+ unsigned int offset;
+};
+
+static const char format[] =
+"The table format must be:\n"
+"\tTwo integers describing the dimensions of the codebook.\n"
+"\tThen, enough numbers to fill the specified dimensions.\n";
+
+float get_float(FILE * in, const char * name, char * * cursor, char * buffer, int size)
+{
+ for ( ; ; ) {
+ char * s = *cursor;
+ char c;
+
+ while ( (c = *s) != '\0' && !isdigit(c) && c != '-' && c != '.' )
+ s++;
+
+ /* Comments start with "#" and continue to the end of the line. */
+ if ( c != '\0' && c != '#' ) {
+ char * end = 0;
+ float f = 0;
+
+ f = strtod(s, &end);
+
+ if ( end != s )
+ *cursor = end;
+ return f;
+ }
+
+ if ( fgets(buffer, size, in) == NULL ) {
+ fprintf(stderr, "%s: Format error. %s\n", name, format);
+ exit(1);
+ }
+ *cursor = buffer;
+ }
+}
+
+static struct codebook *load(const char * name)
+{
+ FILE *file;
+ char line[2048];
+ char *cursor = line;
+ struct codebook *b = malloc(sizeof(struct codebook));
+ int i;
+ int size;
+ float angle;
+
+ file = fopen(name, "rt");
+ assert(file != NULL);
+
+ *cursor = '\0';
+
+ b->k = (int)get_float(file, name, &cursor, line, sizeof(line));
+ b->m = (int)get_float(file, name ,&cursor, line, sizeof(line));
+ size = b->k * b->m;
+
+ b->cb = (COMP *)malloc(size * sizeof(COMP));
+
+ for ( i = 0; i < size; i++ ) {
+ angle = get_float(file, name, &cursor, line, sizeof(line));
+ b->cb[i].real = cos(angle);
+ b->cb[i].imag = sin(angle);
+ }
+
+ fclose(file);
+
+ return b;
+}
+
+
+/* states for phase experiments */
+
+struct PEXP {
+ float phi1;
+ float phi_prev[MAX_AMP];
+ float Wo_prev;
+ int frames;
+ float snr;
+ float var;
+ int var_n;
+ struct codebook *vq1,*vq2,*vq3,*vq4,*vq5;
+ float vq_var;
+ int vq_var_n;
+ MODEL prev_model;
+ int state;
+};
+
+
+/*---------------------------------------------------------------------------* \
+
+ phase_experiment_create()
+
+ Inits states for phase quantisation experiments.
+
+\*---------------------------------------------------------------------------*/
+
+struct PEXP * phase_experiment_create() {
+ struct PEXP *pexp;
+ int i;
+
+ pexp = (struct PEXP *)malloc(sizeof(struct PEXP));
+ assert (pexp != NULL);
+
+ pexp->phi1 = 0;
+ for(i=0; iphi_prev[i] = 0.0;
+ pexp->Wo_prev = 0.0;
+ pexp->frames = 0;
+ pexp->snr = 0.0;
+ pexp->var = 0.0;
+ pexp->var_n = 0;
+
+ /* smoothed 10th order for 1st 1 khz */
+ //pexp->vq1 = load("../unittest/ph1_10_1024.txt");
+ //pexp->vq1->offset = 0;
+
+ /* load experimental phase VQ */
+
+ //pexp->vq1 = load("../unittest/testn1_20_1024.txt");
+ pexp->vq1 = load("../unittest/test.txt");
+ //pexp->vq2 = load("../unittest/testn21_40_1024.txt");
+ pexp->vq2 = load("../unittest/test11_20_1024.txt");
+ pexp->vq3 = load("../unittest/test21_30_1024.txt");
+ pexp->vq4 = load("../unittest/test31_40_1024.txt");
+ pexp->vq5 = load("../unittest/test41_60_1024.txt");
+ pexp->vq1->offset = 0;
+ pexp->vq2->offset = 10;
+ pexp->vq3->offset = 20;
+ pexp->vq4->offset = 30;
+ pexp->vq5->offset = 40;
+
+ pexp->vq_var = 0.0;
+ pexp->vq_var_n = 0;
+
+ pexp->state = 0;
+
+ return pexp;
+}
+
+
+/*---------------------------------------------------------------------------* \
+
+ phase_experiment_destroy()
+
+\*---------------------------------------------------------------------------*/
+
+void phase_experiment_destroy(struct PEXP *pexp) {
+ assert(pexp != NULL);
+ if (pexp->snr != 0.0)
+ printf("snr: %4.2f dB\n", pexp->snr/pexp->frames);
+ if (pexp->var != 0.0)
+ printf("var...: %4.3f std dev...: %4.3f (%d non zero phases)\n",
+ pexp->var/pexp->var_n, sqrt(pexp->var/pexp->var_n), pexp->var_n);
+ if (pexp->vq_var != 0.0)
+ printf("vq var: %4.3f vq std dev: %4.3f (%d non zero phases)\n",
+ pexp->vq_var/pexp->vq_var_n, sqrt(pexp->vq_var/pexp->vq_var_n), pexp->vq_var_n);
+ free(pexp);
+}
+
+
+/*---------------------------------------------------------------------------* \
+
+ Various test and experimental functions ................
+
+\*---------------------------------------------------------------------------*/
+
+/* Bubblesort to find highest amplitude harmonics */
+
+struct AMPINDEX {
+ float amp;
+ int index;
+};
+
+static void bubbleSort(struct AMPINDEX numbers[], int array_size)
+{
+ int i, j;
+ struct AMPINDEX temp;
+
+ for (i = (array_size - 1); i > 0; i--)
+ {
+ for (j = 1; j <= i; j++)
+ {
+ //printf("i %d j %d %f %f \n", i, j, numbers[j-1].amp, numbers[j].amp);
+ if (numbers[j-1].amp < numbers[j].amp)
+ {
+ temp = numbers[j-1];
+ numbers[j-1] = numbers[j];
+ numbers[j] = temp;
+ }
+ }
+ }
+}
+
+
+static void print_pred_error(struct PEXP *pexp, MODEL *model, int start, int end, float mag_thresh) {
+ int i;
+ float mag;
+
+ mag = 0.0;
+ for(i=start; i<=end; i++)
+ mag += model->A[i]*model->A[i];
+ mag = 10*log10(mag/(end-start));
+
+ if (mag > mag_thresh) {
+ for(i=start; i<=end; i++) {
+ float pred = pexp->phi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0;
+ float err = pred - model->phi[i];
+ err = atan2(sin(err),cos(err));
+ printf("%f\n",err);
+ }
+ //printf("\n");
+ }
+
+}
+
+
+static void predict_phases(struct PEXP *pexp, MODEL *model, int start, int end) {
+ int i;
+
+ for(i=start; i<=end; i++) {
+ model->phi[i] = pexp->phi_prev[i] + N*i*model->Wo;
+ }
+
+}
+static float refine_Wo(struct PEXP *pexp,
+ MODEL *model,
+ int start,
+ int end);
+
+/* Fancy state based phase prediction. Actually works OK on most utterances,
+ but could use some tuning. Breaks down a bit on mmt1. */
+
+static void predict_phases_state(struct PEXP *pexp, MODEL *model, int start, int end) {
+ int i, next_state;
+ float best_Wo, dWo;
+
+ //best_Wo = refine_Wo(pexp, model, start, end);
+ //best_Wo = (model->Wo + pexp->Wo_prev)/2.0;
+ best_Wo = model->Wo;
+
+ dWo = fabs(model->Wo - pexp->Wo_prev)/model->Wo;
+ next_state = pexp->state;
+ switch(pexp->state) {
+ case 0:
+ if (dWo < 0.1) {
+ /* UV -> V transition, so start with phases in lock. They will
+ drift a bit over voiced track which is kinda what we want, so
+ we don't get clicky speech.
+ */
+ next_state = 1;
+ for(i=start; i<=end; i++)
+ pexp->phi_prev[i] = i*pexp->phi1;
+ }
+
+ break;
+ case 1:
+ if (dWo > 0.1)
+ next_state = 0;
+ break;
+ }
+ pexp->state = next_state;
+
+ if (pexp->state == 0)
+ for(i=start; i<=end; i++) {
+ model->phi[i] = PI*(1.0 - 2.0*rand()/RAND_MAX);
+ }
+ else
+ for(i=start; i<=end; i++) {
+ model->phi[i] = pexp->phi_prev[i] + N*i*best_Wo;
+ }
+ printf("state %d\n", pexp->state);
+}
+
+static void struct_phases(struct PEXP *pexp, MODEL *model, int start, int end) {
+ int i;
+
+ for(i=start; i<=end; i++)
+ model->phi[i] = pexp->phi1*i;
+
+}
+
+
+static void predict_phases2(struct PEXP *pexp, MODEL *model, int start, int end) {
+ int i;
+ float pred, str, diff;
+
+ for(i=start; i<=end; i++) {
+ pred = pexp->phi_prev[i] + N*i*model->Wo;
+ str = pexp->phi1*i;
+ diff = str - pred;
+ diff = atan2(sin(diff), cos(diff));
+ if (diff > 0)
+ pred += PI/16;
+ else
+ pred -= PI/16;
+ model->phi[i] = pred;
+ }
+
+}
+
+static void skip_phases(struct PEXP *pexp, MODEL *model, int start, int end) {
+ int i;
+
+ for(i=start; i<=end; i+=2)
+ model->phi[i] = model->phi[i-1] - model->phi[i-2];
+
+}
+
+static void rand_phases(MODEL *model, int start, int end) {
+ int i;
+
+ for(i=start; i<=end; i++)
+ model->phi[i] = PI*(1.0 - 2.0*(float)rand()/RAND_MAX);
+
+}
+
+static void quant_phase(float *phase, float min, float max, int bits) {
+ int levels = 1 << bits;
+ int index;
+ float norm, step;
+
+ norm = (*phase - min)/(max - min);
+ index = floor(levels*norm);
+
+ //printf("phase %f norm %f index %d ", *phase, norm, index);
+ if (index < 0 ) index = 0;
+ if (index > (levels-1)) index = levels-1;
+ //printf("index %d ", index);
+ step = (max - min)/levels;
+ *phase = min + step*index + 0.5*step;
+ //printf("step %f phase %f\n", step, *phase);
+}
+
+static void quant_phases(MODEL *model, int start, int end, int bits) {
+ int i;
+
+ for(i=start; i<=end; i++) {
+ quant_phase(&model->phi[i], -PI, PI, bits);
+ }
+}
+
+static void fixed_bits_per_frame(struct PEXP *pexp, MODEL *model, int m, int budget) {
+ int res, finished;
+
+ res = 3;
+ finished = 0;
+
+ while(!finished) {
+ if (m > model->L/2)
+ res = 2;
+ if (((budget - res) < 0) || (m > model->L))
+ finished = 1;
+ else {
+ quant_phase(&model->phi[m], -PI, PI, res);
+ budget -= res;
+ m++;
+ }
+ }
+ printf("m: %d L: %d budget: %d\n", m, model->L, budget);
+ predict_phases(pexp, model, m, model->L);
+ //rand_phases(model, m, model->L);
+}
+
+/* used to plot histogram of quantisation error, for 3 bits, 8 levels,
+ should be uniform between +/- PI/8 */
+
+static void check_phase_quant(MODEL *model, float tol)
+{
+ int m;
+ float phi_before[MAX_AMP];
+
+ for(m=1; m<=model->L; m++)
+ phi_before[m] = model->phi[m];
+
+ quant_phases(model, 1, model->L, 3);
+
+ for(m=1; m<=model->L; m++) {
+ float err = phi_before[m] - model->phi[m];
+ printf("%f\n", err);
+ if (fabs(err) > tol)
+ exit(0);
+ }
+}
+
+
+static float est_phi1(MODEL *model, int start, int end)
+{
+ int m;
+ float delta, s, c, phi1_est;
+
+ if (end > model->L)
+ end = model->L;
+
+ s = c = 0.0;
+ for(m=start; mphi[m+1] - model->phi[m];
+ s += sin(delta);
+ c += cos(delta);
+ }
+
+ phi1_est = atan2(s,c);
+
+ return phi1_est;
+}
+
+static void print_phi1_pred_error(MODEL *model, int start, int end)
+{
+ int m;
+ float phi1_est;
+
+ phi1_est = est_phi1(model, start, end);
+
+ for(m=start; mphi[m+1] - model->phi[m] - phi1_est;
+ err = atan2(sin(err),cos(err));
+ printf("%f\n", err);
+ }
+}
+
+
+static void first_order_band(MODEL *model, int start, int end, float phi1_est)
+{
+ int m;
+ float pred_err, av_pred_err;
+ float c,s;
+
+ s = c = 0.0;
+ for(m=start; mphi[m] - phi1_est*m;
+ s += sin(pred_err);
+ c += cos(pred_err);
+ }
+
+ av_pred_err = atan2(s,c);
+ for(m=start; mphi[m] = av_pred_err + phi1_est*m;
+ model->phi[m] = atan2(sin(model->phi[m]), cos(model->phi[m]));
+ }
+
+}
+
+
+static void sub_linear(MODEL *model, int start, int end, float phi1_est)
+{
+ int m;
+
+ for(m=start; mphi[m] = m*phi1_est;
+ }
+}
+
+
+static void top_amp(struct PEXP *pexp, MODEL *model, int start, int end, int n_harm, int pred)
+{
+ int removed = 0, not_removed = 0;
+ int top, i, j;
+ struct AMPINDEX sorted[MAX_AMP];
+
+ /* sort into ascending order of amplitude */
+
+ printf("\n");
+ for(i=start,j=0; iA[i];
+ sorted[j].index = i;
+ printf("%f ", model->A[i]);
+ }
+ bubbleSort(sorted, end-start);
+
+ printf("\n");
+ for(j=0; jA[i] == sorted[j].amp) {
+ top = 1;
+ assert(i == sorted[j].index);
+ }
+ }
+
+ #define ALTTOP
+ #ifdef ALTTOP
+ model->phi[i] = 0.0; /* make sure */
+ if (top) {
+ model->phi[i] = i*pexp->phi1;
+ removed++;
+ }
+ else {
+ model->phi[i] = PI*(1.0 - 2.0*(float)rand()/RAND_MAX); // note: try rand for higher harms
+ removed++;
+ }
+ #else
+ if (!top) {
+ model->phi[i] = 0.0; /* make sure */
+ if (pred) {
+ //model->phi[i] = pexp->phi_prev[i] + i*N*(model->Wo + pexp->Wo_prev)/2.0;
+ model->phi[i] = i*model->phi[1];
+ }
+ else
+ model->phi[i] = PI*(1.0 - 2.0*(float)rand()/RAND_MAX); // note: try rand for higher harms
+ removed++;
+ }
+ else {
+ /* need to make this work thru budget of bits */
+ quant_phase(&model->phi[i], -PI, PI, 3);
+ not_removed++;
+ }
+ #endif
+ }
+ printf("dim: %d rem %d not_rem %d\n", end-start, removed, not_removed);
+
+}
+
+
+static void limit_prediction_error(struct PEXP *pexp, MODEL *model, int start, int end, float limit)
+{
+ int i;
+ float pred, pred_error, error;
+
+ for(i=start; i<=end; i++) {
+ pred = pexp->phi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0;
+ pred_error = pred - model->phi[i];
+ pred_error -= TWO_PI*floor((pred_error+PI)/TWO_PI);
+ quant_phase(&pred_error, -limit, limit, 2);
+
+ error = pred - pred_error - model->phi[i];
+ error -= TWO_PI*floor((error+PI)/TWO_PI);
+ printf("%f\n", pred_error);
+ model->phi[i] = pred - pred_error;
+ }
+}
+
+
+static void quant_prediction_error(struct PEXP *pexp, MODEL *model, int start, int end, float limit)
+{
+ int i;
+ float pred, pred_error;
+
+ for(i=start; i<=end; i++) {
+ pred = pexp->phi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0;
+ pred_error = pred - model->phi[i];
+ pred_error -= TWO_PI*floor((pred_error+PI)/TWO_PI);
+
+ printf("%f\n", pred_error);
+ model->phi[i] = pred - pred_error;
+ }
+}
+
+
+static void print_sparse_pred_error(struct PEXP *pexp, MODEL *model, int start, int end, float mag_thresh)
+{
+ int i, index;
+ float mag, pred, error;
+ float sparse_pe[MAX_AMP];
+
+ mag = 0.0;
+ for(i=start; i<=end; i++)
+ mag += model->A[i]*model->A[i];
+ mag = 10*log10(mag/(end-start));
+
+ if (mag > mag_thresh) {
+ for(i=0; iphi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0;
+ error = pred - model->phi[i];
+ error = atan2(sin(error),cos(error));
+
+ index = MAX_AMP*i*model->Wo/PI;
+ assert(index < MAX_AMP);
+ sparse_pe[index] = error;
+ }
+
+ /* dump spare phase vector in polar format */
+
+ for(i=0; iL; m++) {
+ signal += model->A[m]*model->A[m];
+ diff = cos(model->phi[m]) - cos(before[m]);
+ noise += pow(model->A[m]*diff, 2.0);
+ diff = sin(model->phi[m]) - sin(before[m]);
+ noise += pow(model->A[m]*diff, 2.0);
+ //printf("%f %f\n", before[m], model->phi[m]);
+ }
+ //printf("%f %f snr = %f\n", signal, noise, 10.0*log10(signal/noise));
+ pexp->snr += 10.0*log10(signal/noise);
+}
+
+
+static void update_variance_calc(struct PEXP *pexp, MODEL *model, float before[])
+{
+ int m;
+ float diff;
+
+ for(m=1; mL; m++) {
+ diff = model->phi[m] - before[m];
+ diff = atan2(sin(diff), cos(diff));
+ pexp->var += diff*diff;
+ }
+ pexp->var_n += model->L;
+}
+
+void print_vec(COMP cb[], int d, int e)
+{
+ int i,j;
+
+ for(j=0; jWo + pexp->Wo_prev)/2.0;
+ best_var = 1E32;
+ for(Wo=0.97*Wo_est; Wo<=1.03*Wo_est; Wo+=0.001*Wo_est) {
+
+ /* predict phase and sum differences between harmonics */
+
+ var = 0.0;
+ for(i=start; i<=end; i++) {
+ pred = pexp->phi_prev[i] + N*i*Wo;
+ error = pred - model->phi[i];
+ error = atan2(sin(error),cos(error));
+ var += error*error;
+ }
+
+ if (var < best_var) {
+ best_var = var;
+ best_Wo = Wo;
+ }
+ }
+
+ return best_Wo;
+}
+
+
+static void split_vq(COMP sparse_pe_out[], struct PEXP *pexp, struct codebook *vq, float weights[], COMP sparse_pe_in[])
+{
+ int i, j, non_zero, vq_ind;
+
+ //printf("\n offset %d k %d m %d j: ", vq->offset, vq->k, vq->m);
+ vq_ind = vq_phase(vq->cb, &sparse_pe_in[vq->offset], &weights[vq->offset], vq->k, vq->m, &pexp->vq_var);
+
+ non_zero = 0;
+ for(i=0, j=vq->offset; ik; i++,j++) {
+ //printf("%f ", atan2(sparse_pe[i].imag, sparse_pe[i].real));
+ if ((sparse_pe_in[j].real != 0.0) && (sparse_pe_in[j].imag != 0.0)) {
+ //printf("%d ", j);
+ sparse_pe_out[j] = vq->cb[vq->k * vq_ind + i];
+ non_zero++;
+ }
+ }
+ pexp->vq_var_n += non_zero;
+}
+
+
+static void sparse_vq_pred_error(struct PEXP *pexp,
+ MODEL *model
+)
+{
+ int i, index;
+ float pred, error, error_q_angle, best_Wo;
+ COMP sparse_pe_in[MAX_AMP], sparse_pe_out[MAX_AMP];
+ float weights[MAX_AMP];
+ COMP error_q_rect;
+
+ best_Wo = refine_Wo(pexp, model, 1, model->L);
+ //best_Wo = (model->Wo + pexp->Wo_prev)/2.0;
+
+ /* transform to sparse pred error vector */
+
+ for(i=0; iL; i++) {
+ pred = pexp->phi_prev[i] + N*i*best_Wo;
+ error = pred - model->phi[i];
+
+ index = MAX_AMP*i*model->Wo/PI;
+ assert(index < MAX_AMP);
+ sparse_pe_in[index].real = cos(error);
+ sparse_pe_in[index].imag = sin(error);
+ sparse_pe_out[index] = sparse_pe_in[index];
+ weights[index] = model->A[i];
+ //printf("%d ", index);
+ }
+
+ /* vector quantise */
+
+ split_vq(sparse_pe_out, pexp, pexp->vq1, weights, sparse_pe_in);
+ split_vq(sparse_pe_out, pexp, pexp->vq2, weights, sparse_pe_in);
+ split_vq(sparse_pe_out, pexp, pexp->vq3, weights, sparse_pe_in);
+ split_vq(sparse_pe_out, pexp, pexp->vq4, weights, sparse_pe_in);
+ split_vq(sparse_pe_out, pexp, pexp->vq5, weights, sparse_pe_in);
+
+ /* transform quantised phases back */
+
+ for(i=1; i<=model->L; i++) {
+ pred = pexp->phi_prev[i] + N*i*best_Wo;
+
+ index = MAX_AMP*i*model->Wo/PI;
+ assert(index < MAX_AMP);
+ error_q_rect = sparse_pe_out[index];
+ error_q_angle = atan2(error_q_rect.imag, error_q_rect.real);
+ model->phi[i] = pred - error_q_angle;
+ model->phi[i] = atan2(sin(model->phi[i]), cos(model->phi[i]));
+ }
+}
+
+
+/*
+ est delta (in Wo)
+ see if this gives us a better (smaller variance) prediction error
+*/
+
+static void print_pred_error_sparse_wo_correction(struct PEXP *pexp,
+ MODEL *model,
+ int start, int end,
+ float mag_thresh)
+{
+ int i, index;
+ float mag, pred, error[MAX_AMP], diff, c, s, delta, err;
+ float sparse_pe[MAX_AMP];
+
+ mag = 0.0;
+ for(i=start; i<=end; i++)
+ mag += model->A[i]*model->A[i];
+ mag = 10*log10(mag/(end-start));
+
+ if (mag > mag_thresh) {
+ for(i=0; iphi[i] = pexp->phi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0 + 0.01*N*i;
+ pred = pexp->phi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0;
+ error[i] = pred - model->phi[i];
+ }
+
+ /* estimate delta Wo */
+
+ c = s = 0.0;
+ for(i=start+1; i<=end; i++) {
+ diff = error[i] - error[i-1];
+ c += log(model->A[i])*cos(diff);
+ s += log(model->A[i])*sin(diff);
+ }
+ delta = atan2(s,c)/N;
+ //printf("delta %f\n",delta);
+ delta = 0;
+ /* now predict phases using corrected Wo */
+
+ for(i=start; i<=end; i++) {
+ pred = pexp->phi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0 - N*i*delta;
+ err = pred - model->phi[i];
+ err = atan2(sin(err),cos(err));
+
+ index = MAX_AMP*i*model->Wo/PI;
+ assert(index < MAX_AMP);
+ sparse_pe[index] = err;
+ }
+
+ /* dump spare phase vector in polar format */
+
+ for(i=0; iA[i]*model->A[i];
+ mag = 10*log10(mag/(end-start));
+
+ if (mag > mag_thresh) {
+
+ best_Wo = refine_Wo(pexp, model, start, end);
+
+ /* now predict phases using corrected Wo */
+
+ for(i=0; iphi_prev[i] + N*i*best_Wo;
+ err = pred - model->phi[i];
+ err = atan2(sin(err),cos(err));
+
+ index = MAX_AMP*i*model->Wo/PI;
+ assert(index < MAX_AMP);
+ sparse_pe[index] = err;
+ }
+
+ /* dump spare phase vector in polar format */
+
+ for(i=0; iL);
+
+ for(i=start; i<=end; i++) {
+ model->phi[i] = pexp->phi_prev[i] + N*i*best_Wo;
+ }
+}
+
+
+/*
+ This functions tests theory that some bands can be combined together
+ due to less frequency resolution at higher frequencies. This will
+ reduce the amount of information we need to encode.
+*/
+
+void smooth_phase(struct PEXP *pexp, MODEL *model, int mode)
+{
+ int m, i, j, index, step, v, en, nav, st;
+ COMP sparse_pe_in[MAX_AMP], av;
+ COMP sparse_pe_out[MAX_AMP];
+ COMP smoothed[MAX_AMP];
+ float best_Wo, pred, err;
+ float weights[MAX_AMP];
+ float avw, smoothed_weights[MAX_AMP];
+ COMP smoothed_in[MAX_AMP], smoothed_out[MAX_AMP];
+
+ best_Wo = refine_Wo(pexp, model, 1, model->L);
+
+ for(m=0; mL; m++) {
+ pred = pexp->phi_prev[m] + N*m*best_Wo;
+ err = model->phi[m] - pred;
+ err = atan2(sin(err),cos(err));
+
+ index = MAX_AMP*m*model->Wo/PI;
+ assert(index < MAX_AMP);
+ sparse_pe_in[index].real = model->A[m]*cos(err);
+ sparse_pe_in[index].imag = model->A[m]*sin(err);
+ sparse_pe_out[index] = sparse_pe_in[index];
+ weights[index] = model->A[m];
+ }
+
+ /* now combine samples at high frequencies to reduce dimension */
+
+ step = 2;
+ st = 0;
+ for(i=st,v=0; i (MAX_AMP-1))
+ en = MAX_AMP-1;
+ for(j=i; jvq1, smoothed_weights, smoothed_in);
+ for(i=0; i (MAX_AMP-1))
+ en = MAX_AMP-1;
+ for(j=i; jL; m++) {
+ index = MAX_AMP*m*model->Wo/PI;
+ assert(index < MAX_AMP);
+ pred = pexp->phi_prev[m] + N*m*best_Wo;
+ err = atan2(sparse_pe_out[index].imag, sparse_pe_out[index].real);
+ model->phi[m] = pred + err;
+ }
+
+}
+
+/*
+ Another version of a functions that tests the theory that some bands
+ can be combined together due to less frequency resolution at higher
+ frequencies. This will reduce the amount of information we need to
+ encode.
+*/
+
+void smooth_phase2(struct PEXP *pexp, MODEL *model) {
+ float m;
+ float step;
+ int a,b,h,i;
+ float best_Wo, pred, err, s,c, phi1_;
+
+ best_Wo = refine_Wo(pexp, model, 1, model->L);
+
+ step = (float)model->L/30;
+ printf("\nL: %d step: %3.2f am,bm: ", model->L, step);
+ for(m=(float)model->L/4; m<=model->L; m+=step) {
+ a = floor(m);
+ b = floor(m+step);
+ if (b > model->L) b = model->L;
+ h = b-a;
+
+ printf("%d,%d,(%d) ", a, b, h);
+ c = s = 0.0;
+ if (h>1) {
+ for(i=a; iphi_prev[i] + N*i*best_Wo;
+ err = model->phi[i] - pred;
+ c += cos(err); s += sin(err);
+ }
+ phi1_ = atan2(s,c);
+ for(i=a; iphi_prev[i] + N*i*best_Wo;
+ printf("%d: %4.3f -> ", i, model->phi[i]);
+ model->phi[i] = pred + phi1_;
+ model->phi[i] = atan2(sin(model->phi[i]),cos(model->phi[i]));
+ printf("%4.3f ", model->phi[i]);
+ }
+ }
+ }
+}
+
+
+#define MAX_BINS 40
+//static float bins[] = {2600.0, 2800.0, 3000.0, 3200.0, 3400.0, 3600.0, 3800.0, 4000.0};
+static float bins[] = {/*
+
+ 1000.0, 1100.0, 1200.0, 1300.0, 1400.0,
+ 1500.0, 1600.0, 1700.0, 1800.0, 1900.0,*/
+
+ 2000.0, 2400.0, 2800.0,
+ 3000.0, 3400.0, 3600.0, 4000.0};
+
+void smooth_phase3(struct PEXP *pexp, MODEL *model) {
+ int m, i;
+ int nbins;
+ int b;
+ float f, best_Wo, pred, err;
+ COMP av[MAX_BINS];
+
+ nbins = sizeof(bins)/sizeof(float);
+ best_Wo = refine_Wo(pexp, model, 1, model->L);
+
+ /* clear all bins */
+
+ for(i=0; iL; m++) {
+ f = m*model->Wo*FS/TWO_PI;
+ if (f > bins[0]) {
+
+ /* find bin */
+
+ for(i=0; i bins[i]) && (f <= bins[i+1]))
+ b = i;
+ assert(b < MAX_BINS);
+
+ /* est predicted phase from average */
+
+ pred = pexp->phi_prev[m] + N*m*best_Wo;
+ err = model->phi[m] - pred;
+ av[b].real += cos(err); av[b].imag += sin(err);
+ }
+
+ }
+
+ /* use averages to est phases */
+
+ for(m=1; m<=model->L; m++) {
+ f = m*model->Wo*FS/TWO_PI;
+ if (f > bins[0]) {
+
+ /* find bin */
+
+ for(i=0; i bins[i]) && (f <= bins[i+1]))
+ b = i;
+ assert(b < MAX_BINS);
+
+ /* add predicted phase error to this bin */
+
+ printf("L %d m %d f %4.f b %d\n", model->L, m, f, b);
+
+ pred = pexp->phi_prev[m] + N*m*best_Wo;
+ err = atan2(av[b].imag, av[b].real);
+ printf(" %d: %4.3f -> ", m, model->phi[m]);
+ model->phi[m] = pred + err;
+ model->phi[m] = atan2(sin(model->phi[m]),cos(model->phi[m]));
+ printf("%4.3f\n", model->phi[m]);
+ }
+ }
+ printf("\n");
+}
+
+
+/*
+ Try to code the phase of the largest amplitude in each band. Randomise the
+ phase of the other harmonics. The theory is that only the largest harmonic
+ will be audible.
+*/
+
+void cb_phase1(struct PEXP *pexp, MODEL *model) {
+ int m, i;
+ int nbins;
+ int b;
+ float f, best_Wo;
+ float max_val[MAX_BINS];
+ int max_ind[MAX_BINS];
+
+ nbins = sizeof(bins)/sizeof(float);
+ best_Wo = refine_Wo(pexp, model, 1, model->L);
+
+ for(i=0; iL; m++) {
+ f = m*model->Wo*FS/TWO_PI;
+ if (f > bins[0]) {
+
+ /* find bin */
+
+ for(i=0; i bins[i]) && (f <= bins[i+1]))
+ b = i;
+ assert(b < MAX_BINS);
+
+ if (model->A[m] > max_val[b]) {
+ max_val[b] = model->A[m];
+ max_ind[b] = m;
+ }
+ }
+
+ }
+
+ /* randomise phase of other harmonics */
+
+ for(m=1; m<=model->L; m++) {
+ f = m*model->Wo*FS/TWO_PI;
+ if (f > bins[0]) {
+
+ /* find bin */
+
+ for(i=0; i bins[i]) && (f <= bins[i+1]))
+ b = i;
+ assert(b < MAX_BINS);
+
+ if (m != max_ind[b])
+ model->phi[m] = pexp->phi_prev[m] + N*m*best_Wo;
+ }
+ }
+}
+
+
+/*
+ Theory is only the phase of the envelope of signal matters within a
+ Critical Band. So we estimate the position of an impulse that
+ approximates the envelope of the signal.
+*/
+
+void cb_phase2(struct PEXP *pexp, MODEL *model) {
+ int st, m, i, a, b, step;
+ float diff,w,c,s,phi1_;
+ float A[MAX_AMP];
+
+ for(m=1; m<=model->L; m++) {
+ A[m] = model->A[m];
+ model->A[m] = 0;
+ }
+
+ st = 2*model->L/4;
+ step = 3;
+ model->phi[1] = pexp->phi_prev[1] + (pexp->Wo_prev+model->Wo)*N/2.0;
+
+ printf("L=%d ", model->L);
+ for(m=st; m model->L)
+ b = model->L;
+
+ c = s = 0;
+ for(i=a; iphi[i+1] - model->phi[i];
+ //w = (model->A[i+1] + model->A[i])/2;
+ w = 1.0;
+ c += w*cos(diff); s += w*sin(diff);
+ }
+ phi1_ = atan2(s,c);
+ printf("replacing: ");
+ for(i=a; iphi[i] = i*phi1_;
+ //model->phi[i] = i*model->phi[1];
+ //model->phi[i] = m*(pexp->Wo_prev+model->Wo)*N/2.0;
+ model->A[m] = A[m];
+ printf("%d ", i);
+ }
+ printf(" . ");
+ }
+ printf("\n");
+}
+
+
+static void smooth_phase4(MODEL *model) {
+ int m;
+ float phi_m, phi_m_1;
+
+ if (model->L > 25) {
+ printf("\nL %d\n", model->L);
+ for(m=model->L/2; m<=model->L; m+=2) {
+ if ((m+1) <= model->L) {
+ phi_m = (model->phi[m] - model->phi[m+1])/2.0;
+ phi_m_1 = (model->phi[m+1] - model->phi[m])/2.0;
+ model->phi[m] = phi_m;
+ model->phi[m+1] = phi_m_1;
+ printf("%d %4.3f %4.3f ", m, phi_m, phi_m_1);
+ }
+ }
+ }
+
+}
+
+/* try repeating last frame, just advance phases to account for time shift */
+
+static void repeat_phases(struct PEXP *pexp, MODEL *model) {
+ int m;
+
+ *model = pexp->prev_model;
+ for(m=1; m<=model->L; m++)
+ model->phi[m] += N*m*model->Wo;
+
+}
+
+/*---------------------------------------------------------------------------*\
+
+ phase_experiment()
+
+ Phase quantisation experiments.
+
+\*---------------------------------------------------------------------------*/
+
+void phase_experiment(struct PEXP *pexp, MODEL *model, char *arg) {
+ int m;
+ float before[MAX_AMP], best_Wo;
+
+ assert(pexp != NULL);
+ memcpy(before, &model->phi[0], sizeof(float)*MAX_AMP);
+
+ if (strcmp(arg,"q3") == 0) {
+ quant_phases(model, 1, model->L, 3);
+ update_snr_calc(pexp, model, before);
+ update_variance_calc(pexp, model, before);
+ }
+
+ if (strcmp(arg,"dec2") == 0) {
+ if ((pexp->frames % 2) != 0) {
+ predict_phases(pexp, model, 1, model->L);
+ update_snr_calc(pexp, model, before);
+ update_variance_calc(pexp, model, before);
+ }
+ }
+
+ if (strcmp(arg,"repeat") == 0) {
+ if ((pexp->frames % 2) != 0) {
+ repeat_phases(pexp, model);
+ update_snr_calc(pexp, model, before);
+ update_variance_calc(pexp, model, before);
+ }
+ }
+
+ if (strcmp(arg,"vq") == 0) {
+ sparse_vq_pred_error(pexp, model);
+ update_snr_calc(pexp, model, before);
+ update_variance_calc(pexp, model, before);
+ }
+
+ if (strcmp(arg,"pred") == 0)
+ predict_phases_state(pexp, model, 1, model->L);
+
+ if (strcmp(arg,"pred1k") == 0)
+ predict_phases(pexp, model, 1, model->L/4);
+
+ if (strcmp(arg,"smooth") == 0) {
+ smooth_phase(pexp, model,0);
+ update_snr_calc(pexp, model, before);
+ }
+ if (strcmp(arg,"smoothtrain") == 0)
+ smooth_phase(pexp, model,1);
+ if (strcmp(arg,"smoothvq") == 0) {
+ smooth_phase(pexp, model,2);
+ update_snr_calc(pexp, model, before);
+ }
+
+ if (strcmp(arg,"smooth2") == 0)
+ smooth_phase2(pexp, model);
+ if (strcmp(arg,"smooth3") == 0)
+ smooth_phase3(pexp, model);
+ if (strcmp(arg,"smooth4") == 0)
+ smooth_phase4(model);
+ if (strcmp(arg,"vqsmooth3") == 0) {
+ sparse_vq_pred_error(pexp, model);
+ smooth_phase3(pexp, model);
+ }
+
+ if (strcmp(arg,"cb1") == 0) {
+ cb_phase1(pexp, model);
+ update_snr_calc(pexp, model, before);
+ }
+
+ if (strcmp(arg,"top") == 0) {
+ //top_amp(pexp, model, 1, model->L/4, 4, 1);
+ //top_amp(pexp, model, model->L/4, model->L/3, 4, 1);
+ //top_amp(pexp, model, model->L/3+1, model->L/2, 4, 1);
+ //top_amp(pexp, model, model->L/2, model->L, 6, 1);
+ //rand_phases(model, model->L/2, 3*model->L/4);
+ //struct_phases(pexp, model, model->L/2, 3*model->L/4);
+ //update_snr_calc(pexp, model, before);
+ }
+
+ if (strcmp(arg,"pred23") == 0) {
+ predict_phases2(pexp, model, model->L/2, model->L);
+ update_snr_calc(pexp, model, before);
+ }
+
+ if (strcmp(arg,"struct23") == 0) {
+ struct_phases(pexp, model, model->L/2, 3*model->L/4 );
+ update_snr_calc(pexp, model, before);
+ }
+
+ if (strcmp(arg,"addnoise") == 0) {
+ int m;
+ float max;
+
+ max = 0;
+ for(m=1; m<=model->L; m++)
+ if (model->A[m] > max)
+ max = model->A[m];
+ max = 20.0*log10(max);
+ for(m=1; m<=model->L; m++)
+ if (20.0*log10(model->A[m]) < (max-20)) {
+ model->phi[m] += (PI/4)*(1.0 -2.0*rand()/RAND_MAX);
+ //printf("m %d\n", m);
+ }
+ }
+
+ /* normalise phases */
+
+ for(m=1; m<=model->L; m++)
+ model->phi[m] = atan2(sin(model->phi[m]), cos(model->phi[m]));
+
+ /* update states */
+
+ //best_Wo = refine_Wo(pexp, model, model->L/2, model->L);
+ pexp->phi1 += N*model->Wo;
+
+ for(m=1; m<=model->L; m++)
+ pexp->phi_prev[m] = model->phi[m];
+ pexp->Wo_prev = model->Wo;
+ pexp->frames++;
+ pexp->prev_model = *model;
+}
+
+#ifdef OLD_STUFF
+ //quant_phases(model, 1, model->L, 3);
+ //update_variance_calc(pexp, model, before);
+ //print_sparse_pred_error(pexp, model, 1, model->L, 40.0);
+
+ //sparse_vq_pred_error(pexp, model);
+
+ //quant_phases(model, model->L/4+1, model->L, 3);
+
+ //predict_phases1(pexp, model, 1, model->L/4);
+ //quant_phases(model, model->L/4+1, model->L, 3);
+
+ //quant_phases(model, 1, model->L/8, 3);
+
+ //update_snr_calc(pexp, model, before);
+ //update_variance_calc(pexp, model, before);
+
+ //fixed_bits_per_frame(pexp, model, 40);
+ //struct_phases(pexp, model, 1, model->L/4);
+ //rand_phases(model, 10, model->L);
+ //for(m=1; m<=model->L; m++)
+ // model->A[m] = 0.0;
+ //model->A[model->L/2] = 1000;
+ //repeat_phases(model, 20);
+ //predict_phases(pexp, model, 1, model->L/4);
+ //quant_phases(model, 1, 10, 3);
+ //quant_phases(model, 10, 20, 2);
+ //repeat_phases(model, 20);
+ //rand_phases(model, 3*model->L/4, model->L);
+ // print_phi1_pred_error(model, 1, model->L);
+ //predict_phases(pexp, model, 1, model->L/4);
+ //first_order_band(model, model->L/4, model->L/2);
+ //first_order_band(model, model->L/2, 3*model->L/4);
+ //if (fabs(model->Wo - pexp->Wo_prev)< 0.1*model->Wo)
+
+ //print_pred_error(pexp, model, 1, model->L, 40.0);
+ //print_sparse_pred_error(pexp, model, 1, model->L, 40.0);
+
+ //phi1_est = est_phi1(model, 1, model->L/4);
+ //print_phi1_pred_error(model, 1, model->L/4);
+
+ //first_order_band(model, 1, model->L/4, phi1_est);
+ //sub_linear(model, 1, model->L/4, phi1_est);
+
+ //top_amp(pexp, model, 1, model->L/4, 4);
+ //top_amp(pexp, model, model->L/4, model->L/2, 4);
+
+ //first_order_band(model, 1, model->L/4, phi1_est);
+ //first_order_band(model, model->L/4, model->L/2, phi1_est);
+
+ //if (fabs(model->Wo - pexp->Wo_prev) > 0.2*model->Wo)
+ // rand_phases(model, model->L/2, model->L);
+
+ //top_amp(pexp, model, 1, model->L/4, 4);
+ //top_amp(pexp, model, model->L/4, model->L/2, 8);
+ //top_amp(pexp, model, model->L/4+1, model->L/2, 10, 1);
+ //top_amp(pexp, model, 1, model->L/4, 10, 1);
+ //top_amp(pexp, model, model->L/4+1, 3*model->L/4, 10, 1);
+ //top_amp(pexp, model, 1, 3*model->L/4, 20, 1);
+
+ #ifdef REAS_CAND1
+ predict_phases(pexp, model, 1, model->L/4);
+ top_amp(pexp, model, model->L/4+1, 3*model->L/4, 10, 1);
+ rand_phases(model, 3*model->L/4+1, model->L);
+ #endif
+
+ #ifdef REAS_CAND2
+ if ((pexp->frames % 2) == 0) {
+ //printf("quant\n");
+ predict_phases(pexp, model, 1, model->L/4);
+ //top_amp(pexp, model, model->L/4+1, 3*model->L/4, 20, 1);
+ top_amp(pexp, model, model->L/4+1, 7*model->L/8, 20, 1);
+ rand_phases(model, 7*model->L/8+1, model->L);
+ }
+ else {
+ //printf("predict\n");
+ predict_phases(pexp, model, 1, model->L);
+ }
+ #endif
+
+ //#define REAS_CAND3
+ #ifdef REAS_CAND3
+ if ((pexp->frames % 3) != 0) {
+ printf("pred\n");
+ predict_phases(pexp, model, 1, model->L);
+ }
+ else {
+ predict_phases(pexp, model, 1, model->L/4);
+ fixed_bits_per_frame(pexp, model, model->L/4+1, 60);
+ }
+ #endif
+ //predict_phases(pexp, model, model->L/4, model->L);
+
+
+ //print_pred_error(pexp, model, 1, model->L);
+ //limit_prediction_error(pexp, model, model->L/2, model->L, PI/2);
+#endif
diff --git a/libs/libcodec2/src/phaseexp.h b/libs/libcodec2/src/phaseexp.h
new file mode 100644
index 0000000000..b43db75e83
--- /dev/null
+++ b/libs/libcodec2/src/phaseexp.h
@@ -0,0 +1,39 @@
+/*---------------------------------------------------------------------------*\
+
+ FILE........: phaseexp.h
+ AUTHOR......: David Rowe
+ DATE CREATED: June 2012
+
+ Experimental functions for quantising, modelling and synthesising phase.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+ Copyright (C) 2012 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.1, 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, see .
+*/
+
+#ifndef __PHASEEXP__
+#define __PHASEEXP__
+
+#include "kiss_fft.h"
+
+struct PEXP;
+
+struct PEXP * phase_experiment_create();
+void phase_experiment_destroy(struct PEXP *pexp);
+void phase_experiment(struct PEXP *pexp, MODEL *model, char *arg);
+
+#endif
diff --git a/libs/libcodec2/src/pilot_coeff.h b/libs/libcodec2/src/pilot_coeff.h
new file mode 100644
index 0000000000..66e7501d8f
--- /dev/null
+++ b/libs/libcodec2/src/pilot_coeff.h
@@ -0,0 +1,34 @@
+/* Generated by pilot_coeff_file() Octave function */
+
+const float pilot_coeff[]={
+ 0.00204705,
+ 0.00276339,
+ 0.00432595,
+ 0.00697042,
+ 0.0108452,
+ 0.0159865,
+ 0.0223035,
+ 0.029577,
+ 0.0374709,
+ 0.045557,
+ 0.0533491,
+ 0.0603458,
+ 0.0660751,
+ 0.070138,
+ 0.0722452,
+ 0.0722452,
+ 0.070138,
+ 0.0660751,
+ 0.0603458,
+ 0.0533491,
+ 0.045557,
+ 0.0374709,
+ 0.029577,
+ 0.0223035,
+ 0.0159865,
+ 0.0108452,
+ 0.00697042,
+ 0.00432595,
+ 0.00276339,
+ 0.00204705
+};
diff --git a/libs/libcodec2/src/postfilter.c b/libs/libcodec2/src/postfilter.c
index 6dad76b1e1..c78f495bed 100644
--- a/libs/libcodec2/src/postfilter.c
+++ b/libs/libcodec2/src/postfilter.c
@@ -24,15 +24,16 @@
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.
+ along with this program; if not, see .
*/
+#include
#include
#include
#include
#include "defines.h"
+#include "comp.h"
#include "dump.h"
#include "postfilter.h"
@@ -44,6 +45,11 @@
#define BG_THRESH 40.0 /* only consider low levels signals for bg_est */
#define BG_BETA 0.1 /* averaging filter constant */
+#define BG_MARGIN 6.0 /* harmonics this far above BG noise are
+ randomised. Helped make bg noise less
+ spikey (impulsive) for mmt1, but speech was
+ perhaps a little rougher.
+ */
/*---------------------------------------------------------------------------*\
@@ -61,7 +67,7 @@
(5-12) are required to transmit the frequency selective voicing
information. Mixed excitation also requires accurate voicing
estimation (parameter estimators always break occasionally under
- exceptional condition).
+ exceptional conditions).
In our case we use a post filter approach which requires no
additional bits to be transmitted. The decoder measures the average
@@ -105,6 +111,7 @@ void postfilter(
for(m=1; m<=model->L; m++)
e += model->A[m]*model->A[m];
+ assert(e > 0.0);
e = 10.0*log10(e/model->L);
/* If beneath threhold, update bg estimate. The idea
@@ -121,11 +128,13 @@ void postfilter(
uv = 0;
if (model->voiced)
for(m=1; m<=model->L; m++)
- if (20.0*log10(model->A[m]) < *bg_est) {
+ if (20.0*log10(model->A[m]) < (*bg_est + BG_MARGIN)) {
model->phi[m] = TWO_PI*(float)rand()/RAND_MAX;
uv++;
}
+#ifdef DUMP
dump_bg(e, *bg_est, 100.0*uv/model->L);
+#endif
}
diff --git a/libs/libcodec2/src/postfilter.h b/libs/libcodec2/src/postfilter.h
index d4dd4ae053..bf080b1b65 100644
--- a/libs/libcodec2/src/postfilter.h
+++ b/libs/libcodec2/src/postfilter.h
@@ -22,8 +22,7 @@
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.
+ along with this program; if not, see .
*/
#ifndef __POSTFILTER__
diff --git a/libs/libcodec2/src/quantise.c b/libs/libcodec2/src/quantise.c
index a1cd728112..1153943b9f 100644
--- a/libs/libcodec2/src/quantise.c
+++ b/libs/libcodec2/src/quantise.c
@@ -20,8 +20,8 @@
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.
+ along with this program; if not, see .
+
*/
#include
@@ -36,52 +36,9 @@
#include "quantise.h"
#include "lpc.h"
#include "lsp.h"
-#include "four1.h"
-#include "codebook.h"
+#include "kiss_fft.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 */
- float *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, "/usr/src/freeswitch/libs/libcodec2-1.0/unittest/lspd123.txt"},
- {3, 1024, "/usr/src/freeswitch/libs/libcodec2-1.0/unittest/lspd456.txt"},
- {2, 64, "/usr/src/freeswitch/libs/libcodec2-1.0/unittest/lspd78.txt"},
- {2, 64, "/usr/src/freeswitch/libs/libcodec2-1.0/unittest/lspd910.txt"},
- {0, 0, ""}
-};
-#endif
-
-LSP_CB lsp_q[] = {
- {1,4,16, codebook_lsp1 },
- {1,4,16, codebook_lsp2 },
- {1,4,16, codebook_lsp3 },
- {1,4,16, codebook_lsp4 },
- {1,4,16, codebook_lsp5 },
- {1,4,16, codebook_lsp6 },
- {1,4,16, codebook_lsp7 },
- {1,3,8, codebook_lsp8 },
- {1,3,8, codebook_lsp9 },
- {1,2,4, codebook_lsp10 },
- {0,0,0, NULL },
-};
-
-/* ptr to each codebook */
-
-static float *plsp_cb[MAX_CB];
/*---------------------------------------------------------------------------*\
@@ -99,114 +56,19 @@ 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;
}
-/*---------------------------------------------------------------------------*\
-
- quantise_uniform
-
- Simulates uniform quantising of a float.
-
-\*---------------------------------------------------------------------------*/
-
-void quantise_uniform(float *val, float min, float max, int bits)
-{
- int levels = 1 << (bits-1);
- float norm;
- int index;
-
- /* hard limit to quantiser range */
-
- printf("min: %f max: %f val: %f ", min, max, val[0]);
- if (val[0] < min) val[0] = min;
- if (val[0] > max) val[0] = max;
-
- norm = (*val - min)/(max-min);
- printf("%f norm: %f ", val[0], norm);
- index = fabs(levels*norm + 0.5);
-
- *val = min + index*(max-min)/levels;
-
- printf("index %d val_: %f\n", index, val[0]);
+int lspd_bits(int i) {
+ return lsp_cbd[i].log2m;
}
-/*---------------------------------------------------------------------------*\
-
- lspd_quantise
-
- Simulates differential lsp quantiser
-
-\*---------------------------------------------------------------------------*/
-
-void lsp_quantise(
- float lsp[],
- float lsp_[],
- int order
-)
-{
- int i;
- float dlsp[LPC_MAX];
- float dlsp_[LPC_MAX];
-
- dlsp[0] = lsp[0];
- for(i=1; i 0);
+ mbest = (struct MBEST *)malloc(sizeof(struct MBEST));
+ assert(mbest != NULL);
+
+ mbest->entries = entries;
+ mbest->list = (struct MBEST_LIST *)malloc(entries*sizeof(struct MBEST_LIST));
+ assert(mbest->list != NULL);
+
+ for(i=0; ientries; i++) {
+ for(j=0; jlist[i].index[j] = 0;
+ mbest->list[i].error = 1E32;
+ }
+
+ return mbest;
+}
+
+
+static void mbest_destroy(struct MBEST *mbest) {
+ assert(mbest != NULL);
+ free(mbest->list);
+ free(mbest);
+}
+
+
+/*---------------------------------------------------------------------------*\
+
+ mbest_insert
+
+ Insert the results of a vector to codebook entry comparison. The
+ list is ordered in order or error, so those entries with the
+ smallest error will be first on the list.
+
+\*---------------------------------------------------------------------------*/
+
+static void mbest_insert(struct MBEST *mbest, int index[], float error) {
+ int i, j, found;
+ struct MBEST_LIST *list = mbest->list;
+ int entries = mbest->entries;
+
+ found = 0;
+ for(i=0; ii; j--)
+ list[j] = list[j-1];
+ for(j=0; jentries; i++) {
+ for(j=0; jlist[i].index[j]);
+ printf(" %f\n", mbest->list[i].error);
+ }
+}
+
+
+/*---------------------------------------------------------------------------*\
+
+ mbest_search
+
+ Searches vec[] to a codebbook of vectors, and maintains a list of the mbest
+ closest matches.
+
+\*---------------------------------------------------------------------------*/
+
+static void mbest_search(
+ const float *cb, /* VQ codebook to search */
+ float vec[], /* target vector */
+ float w[], /* weighting vector */
+ int k, /* dimension of vector */
+ int m, /* number on entries in codebook */
+ struct MBEST *mbest, /* list of closest matches */
+ int index[] /* indexes that lead us here */
+)
+{
+ float e;
+ int i,j;
+ float diff;
+
+ for(j=0; jlist[j].index[0];
+ for(i=0; ilist[j].index[1];
+ index[1] = n2 = mbest_stage2->list[j].index[0];
+ for(i=0; i