From 5eba9ab7a48d25c40b85e43b6a373ae76f6faf49 Mon Sep 17 00:00:00 2001
From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Tue, 8 Sep 2009 04:42:10 +0000
Subject: [PATCH] added MBE pitch post processor to nlp.c, still debugging

git-svn-id: https://svn.code.sf.net/p/freetel/code@43 01035d8c-6547-0410-b346-abe4f91aad63
---
 codec2/src/defines.h  |   3 +-
 codec2/src/dump.c     |  47 ++++++++-
 codec2/src/dump.h     |   2 +
 codec2/src/globals.c  |   2 +-
 codec2/src/lpc.c      |   2 +-
 codec2/src/lpc.h      |   2 +
 codec2/src/nlp.c      | 225 ++++++++++++++++++++++++++++--------------
 codec2/src/nlp.h      |  39 ++++++++
 codec2/src/phase.c    |   1 +
 codec2/src/quantise.c |   6 +-
 codec2/src/quantise.h |   2 +-
 codec2/src/refine.c   |  16 ++-
 codec2/src/sinedec.c  |  45 +++++----
 codec2/src/sinenc.c   |  23 ++---
 14 files changed, 292 insertions(+), 123 deletions(-)
 create mode 100644 codec2/src/nlp.h

diff --git a/codec2/src/defines.h b/codec2/src/defines.h
index 13898e2e..c07e7a8f 100644
--- a/codec2/src/defines.h
+++ b/codec2/src/defines.h
@@ -51,7 +51,8 @@
 
 /* Encoder defines */
 
-#define AW_ENC 512		/* maximum encoder analysis window size */
+#define NW      220             /* analysis window size */
+#define AW_ENC  512		/* maximum encoder analysis window size */
 #define FFT_ENC 512		/* size of FFT used for encoder analysis */
 
 /* Decoder defines */
diff --git a/codec2/src/dump.c b/codec2/src/dump.c
index 1b8195ed..1295074b 100644
--- a/codec2/src/dump.c
+++ b/codec2/src/dump.c
@@ -41,6 +41,8 @@ static FILE *fpw = NULL;
 static FILE *flsp = NULL;
 static FILE *fphase = NULL;
 static FILE *fphase_ = NULL;
+static FILE *ffw = NULL;
+static FILE *fe = NULL;
 
 static char  prefix[MAX_STR];
 
@@ -68,6 +70,10 @@ void dump_off(){
 	fclose(fphase);
     if (fphase_ != NULL)
 	fclose(fphase_);
+    if (ffw != NULL)
+	fclose(ffw);
+    if (fe != NULL)
+	fclose(fe);
 }
 
 void dump_Sn(float Sn[]) {
@@ -85,10 +91,10 @@ void dump_Sn(float Sn[]) {
     /* split across two lines to avoid max line length problems */
     /* reconstruct in Octave */
 
-    for(i=0; i<AW_ENC/2; i++)
+    for(i=0; i<M/2; i++)
 	fprintf(fsn,"%f\t",Sn[i]);
     fprintf(fsn,"\n");    
-    for(i=AW_ENC/2; i<AW_ENC; i++)
+    for(i=M/2; i<M; i++)
 	fprintf(fsn,"%f\t",Sn[i]);
     fprintf(fsn,"\n");    
 }
@@ -245,4 +251,41 @@ void dump_lsp(float lsp[]) {
     fprintf(flsp,"\n");    
 }
 
+void dump_Fw(COMP Fw[]) {
+    int i;
+    char s[MAX_STR];
+
+    if (!dumpon) return;
+
+    if (ffw == NULL) {
+	sprintf(s,"%s_fw.txt", prefix);
+	ffw = fopen(s, "wt");
+	assert(ffw != NULL);
+    }
+
+    for(i=0; i<256; i++)
+	fprintf(ffw,"%f\t",Fw[i].real);
+    fprintf(ffw,"\n");    
+}
+
+void dump_e(float e_hz[]) {
+    int i;
+    char s[MAX_STR];
+
+    if (!dumpon) return;
+
+    if (fe == NULL) {
+	sprintf(s,"%s_e.txt", prefix);
+	fe = fopen(s, "wt");
+	assert(fe != NULL);
+    }
+
+    for(i=0; i<500/2; i++)
+	fprintf(fe,"%f\t",e_hz[i]);
+    fprintf(fe,"\n");    
+    for(i=500/2; i<500; i++)
+	fprintf(fe,"%f\t",e_hz[i]);
+    fprintf(fe,"\n");    
+}
+
 
diff --git a/codec2/src/dump.h b/codec2/src/dump.h
index db55a0a4..b309430c 100644
--- a/codec2/src/dump.h
+++ b/codec2/src/dump.h
@@ -40,5 +40,7 @@ void dump_Pw(COMP Pw[]);
 void dump_lsp(float lsp[]);
 void dump_phase(float phase[]);
 void dump_phase_(float phase[]);
+void dump_Fw(COMP Fw[]);
+void dump_e(float e_hz[]);
 
 #endif
diff --git a/codec2/src/globals.c b/codec2/src/globals.c
index da4ecb97..ad92851f 100644
--- a/codec2/src/globals.c
+++ b/codec2/src/globals.c
@@ -31,7 +31,7 @@
 /* Globals used in encoder and decoder */
 
 int frames;		/* number of frames processed so far */
-float Sn[M+AW_ENC/2];	/* float input speech samples */
+float Sn[M];	        /* float input speech samples */
 MODEL model;		/* model parameters for the current frame */
 int Nw;			/* number of samples in analysis window */
 float sig;		/* energy of current frame */
diff --git a/codec2/src/lpc.c b/codec2/src/lpc.c
index 62bbed17..b95a131d 100644
--- a/codec2/src/lpc.c
+++ b/codec2/src/lpc.c
@@ -244,7 +244,7 @@ void weight(
 {
   int i;
   
-  for(i=0; i<=order; i++)
+  for(i=1; i<=order; i++)
     akw[i] = ak[i]*pow(gamma,(float)i);
 }
     
diff --git a/codec2/src/lpc.h b/codec2/src/lpc.h
index dd133340..114add39 100644
--- a/codec2/src/lpc.h
+++ b/codec2/src/lpc.h
@@ -29,6 +29,8 @@
 #ifndef __LPC__
 #define __LPC__
 
+#define LPC_MAX_ORDER 20
+
 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/codec2/src/nlp.c b/codec2/src/nlp.c
index 2793dce5..c63f9d4f 100644
--- a/codec2/src/nlp.c
+++ b/codec2/src/nlp.c
@@ -26,20 +26,25 @@
   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
 
+#include "nlp.h"
+#include "sine.h"
+#include "dump.h"
+#include <assert.h>
+
 /*---------------------------------------------------------------------------*\
                                                                              
  				DEFINES                                       
                                                                              
 \*---------------------------------------------------------------------------*/
 
-#define PMAX_M 600		/* maximum NLP analysis window size */
-#define COEFF 0.95		/* noth filter parameter */
-#define NTAP 48			/* Decimation LPF order */
-#define PE_FFT_SIZE 512		/* DFT size for pitch estimation */
-#define DEC 5			/* decimation factor */
+#define PMAX_M      600		/* maximum NLP analysis window size     */
+#define COEFF       0.95	/* notch filter parameter               */
+#define PE_FFT_SIZE 512		/* DFT size for pitch estimation        */
+#define DEC         5		/* decimation factor                    */
 #define SAMPLE_RATE 8000
-#define PI 3.141592654		/* mathematical constant */
-#define CNLP 0.5		/* post processor constant */
+#define PI          3.141592654	/* mathematical constant                */
+#define T           0.1         /* threshold for local minima candidate */
+#define F0_MAX      500
 
 /*---------------------------------------------------------------------------*\
                                                                             
@@ -100,44 +105,46 @@ float nlp_fir[] = {
   -1.0818124e-03
 };
 
+float test_candidate_mbe(COMP Sw[], float f0);
+extern int frames;
+
 /*---------------------------------------------------------------------------*\
                                                                              
   void nlp()                                                                  
                                                                              
-  Determines the pitch in samples using the NLP algorithm. Returns the 	      
-  fundamental in Hz.						   	      
+  Determines the pitch in samples using the NLP algorithm. Returns the
+  fundamental in Hz.  Note that the actual pitch estimate is for the
+  centre of the M sample Sn[] vector, not the current N sample input
+  vector.  This is (I think) a delay of 2.5 frames with N=80 samples.
+  You should align further analysis using this pitch estimate to be
+  centred on the middle of Sn[].
                                                                              
 \*---------------------------------------------------------------------------*/
 
-float nlp(Sn,n,m,d,pmin,pmax,pitch)
-float Sn[];			/* input speech vector */
-int n;				/* frames shift (no. new samples in Sn[]) */
-int m;				/* analysis window size */
-int d;				/* additional delay (used for testing) */
-int pmin;			/* minimum pitch value */
-int pmax;			/* maximum pitch value */
-float *pitch;			/* estimated pitch */
+float nlp(
+  float  Sn[],			/* input speech vector */
+  int    n,			/* frames shift (no. new samples in Sn[]) */
+  int    m,			/* analysis window size */
+  int    d,			/* additional delay (used for testing) */
+  int    pmin,			/* minimum pitch value */
+  int    pmax,			/* maximum pitch value */
+  float *pitch,			/* estimated pitch period in samples */
+  COMP   Sw[]                   /* Freq domain version of Sn[] */
+)
 {
   static float sq[PMAX_M];	/* squared speech samples */
-  float notch;			/* current notch filter output */
+  float  notch;			/* current notch filter output */
   static float mem_x,mem_y;     /* memory for notch filter */
-  static float mem_fir[NTAP];	/* decimation FIR filter memory */
-  COMP Fw[PE_FFT_SIZE];		/* DFT of squared signal */
-
-  int gmax_bin;			/* DFT bin where global maxima occurs */
-  float gmax;			/* global maxima value */
-  float lmax;			/* current local maxima value */
-  int lmax_bin;			/* bin of current local maxima */
-  float cmax;			/* chosen local maxima value */
-  int cmax_bin;			/* bin of chosen local maxima */
-
-  int mult;			/* current submultiple */
-  int min_bin;			/* lowest possible bin */
-  int bmin,bmax;		/* range of local maxima search */
-  float thresh;			/* threshold for submultiple selection */
-
-  float F0;			/* fundamental frequency */
-  int i,j,b;
+  static float mem_fir[NLP_NTAP];/* decimation FIR filter memory */
+  COMP   Fw[PE_FFT_SIZE];	/* DFT of squared signal */
+  float  gmax;
+ 
+  float candidate_f0;
+  float f0,best_f0;		/* fundamental frequency */
+  float e,e_min;                /* MBE cost function */
+  int i,j;
+  float e_hz[F0_MAX];
+  int bin;
 
   /* Square, notch filter at DC, and LP filter vector */
 
@@ -154,12 +161,12 @@ float *pitch;			/* estimated pitch */
 
   for(i=m-n+d; i<m+d; i++) {	/* FIR filter vector */
 
-    for(j=0; j<NTAP-1; j++)
+    for(j=0; j<NLP_NTAP-1; j++)
       mem_fir[j] = mem_fir[j+1];
-    mem_fir[NTAP-1] = sq[i];
+    mem_fir[NLP_NTAP-1] = sq[i];
 
     sq[i] = 0.0;
-    for(j=0; j<NTAP; j++)
+    for(j=0; j<NLP_NTAP; j++)
       sq[i] += mem_fir[j]*nlp_fir[j];
   }
 
@@ -175,56 +182,128 @@ float *pitch;			/* estimated pitch */
   for(i=0; i<PE_FFT_SIZE; i++)
     Fw[i].real = Fw[i].real*Fw[i].real + Fw[i].imag*Fw[i].imag;
 
-  /* find global peak within limits, this corresponds to F0 estimate */
+  dump_Fw(Fw);
+
+  /* find global peak */
 
   gmax = 0.0;
   for(i=PE_FFT_SIZE*DEC/pmax; i<=PE_FFT_SIZE*DEC/pmin; i++) {
     if (Fw[i].real > gmax) {
       gmax = Fw[i].real;
-      gmax_bin = i;
     }
   }
 
-  /* Now post process estimate by searching submultiples */
-
-  mult = 2;
-  min_bin = PE_FFT_SIZE*DEC/pmax;
-  thresh = CNLP*gmax;
-  cmax_bin = gmax_bin;
-
-  while(gmax_bin/mult >= min_bin) {
-
-    b = gmax_bin/mult;			/* determine search interval */
-    bmin = 0.8*b;
-    bmax = 1.2*b;
-    if (bmin < min_bin)
-      bmin = min_bin;
-      
-    lmax = 0;
-    for (b=bmin; b<=bmax; b++) 		/* look for maximum in interval */
-      if (Fw[b].real > lmax) {
-        lmax = Fw[b].real;
-	lmax_bin = b;
-      }
-
-    if (lmax > thresh)
-      if (lmax > Fw[lmax_bin-1].real && lmax > Fw[lmax_bin+1].real) {
-	cmax = lmax;
-	cmax_bin = lmax_bin;
-      }
-
-    mult++;
-  }
+  /* Now look for local maxima.  Each local maxima is a candidate
+     that we test using the MBE pitch estimation algotithm */
 
-  F0 = (float)cmax_bin*SAMPLE_RATE/(PE_FFT_SIZE*DEC);
-  *pitch = SAMPLE_RATE/F0;
+  for(i=0; i<F0_MAX; i++)
+      e_hz[i] = -1;
+  e_min = 1E32;
+  best_f0 = 50;
+  for(i=PE_FFT_SIZE*DEC/pmax; i<=PE_FFT_SIZE*DEC/pmin; i++) {
+    if ((Fw[i].real > 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);
+	    if (frames == 29)
+	      printf("candidate F0: %f\n", candidate_f0);
+	    for(f0=candidate_f0-20; f0<=candidate_f0+20; f0+= 2.5) {
+		e = test_candidate_mbe(Sw, f0);
+		bin = floor(f0); assert((bin > 0) && (bin < F0_MAX));
+		e_hz[bin] = e;
+		if (frames == 29)
+		    printf("f0: %f e: %f e_min: %f best_f0: %f\n", 
+			   f0, e, e_min, best_f0);
+		if (e < e_min) {
+		    e_min = e;
+		    best_f0 = f0;
+		}
+	    }
+
+	}
+    }
+  }
+  dump_e(e_hz);
 
   /* Shift samples in buffer to make room for new samples */
 
   for(i=0; i<m-n+d; i++)
     sq[i] = sq[i+n];
 
-  return(F0);  
+  /* return pitch and F0 estimate */
+
+  *pitch = (float)SAMPLE_RATE/best_f0;
+  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[],
+    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;              /* number of bands */
+
+    Wo = f0*(2*PI/SAMPLE_RATE);
+    L = floor(PI/Wo);;   
+
+    error = 0.0;
+
+    /* Just test across the harmonics in the first 1000 Hz (L/4) */
+
+    for(l=1; l<L/4; l++) {
+	Am.real = 0.0;
+	Am.imag = 0.0;
+	den = 0.0;
+	al = ceil((l - 0.5)*Wo*FFT_ENC/TWO_PI);
+	bl = ceil((l + 0.5)*Wo*FFT_ENC/TWO_PI);
+
+	/* Estimate amplitude of harmonic assuming harmonic is totally voiced */
+
+	for(m=al; m<bl; m++) {
+	    offset = FFT_ENC/2 + m - l*Wo*FFT_ENC/TWO_PI + 0.5;
+	    Am.real += Sw[m].real*W[offset].real + Sw[m].imag*W[offset].imag;
+	    Am.imag += Sw[m].imag*W[offset].real - Sw[m].real*W[offset].imag;
+	    den += W[offset].real*W[offset].real + W[offset].imag*W[offset].imag;
+        }
+
+        Am.real = Am.real/den;
+        Am.imag = Am.imag/den;
+
+        /* Determine error between estimated harmonic and original */
+
+        for(m=al; m<bl; m++) {
+	    offset = FFT_ENC/2 + m - l*Wo*FFT_ENC/TWO_PI + 0.5;
+	    Sw_[m].real = Am.real*W[offset].real - Am.imag*W[offset].imag;
+	    Sw_[m].imag = Am.real*W[offset].imag + Am.imag*W[offset].real;
+	    error += (Sw[m].real - Sw_[m].real)*(Sw[m].real - Sw_[m].real);
+	    error += (Sw[m].imag - Sw_[m].imag)*(Sw[m].imag - Sw_[m].imag);
+	}
+    }
+
+    return error;
+}
+
 
diff --git a/codec2/src/nlp.h b/codec2/src/nlp.h
new file mode 100644
index 00000000..5479af37
--- /dev/null
+++ b/codec2/src/nlp.h
@@ -0,0 +1,39 @@
+/*---------------------------------------------------------------------------*\
+                                                 
+  FILE........: nlp.c                                                   
+  AUTHOR......: David Rowe                                      
+  DATE CREATED: 23/3/93                                    
+                                                         
+  Non Linear Pitch (NLP) estimation functions.			  
+                                                               
+\*---------------------------------------------------------------------------*/
+
+/*
+  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 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 General Public License
+  along with this program; if not, write to the Free Software
+  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+*/
+
+#ifndef __NLP__
+#define __NLP__
+
+#include "sine.h"
+
+#define NLP_NTAP 48	 /* Decimation LPF order */
+
+float nlp(float Sn[], int n, int m, int d, int pmin, int pmax, float *pitch,
+	  COMP  Sw[]);
+
+#endif
diff --git a/codec2/src/phase.c b/codec2/src/phase.c
index a88dd2ff..42b1c658 100644
--- a/codec2/src/phase.c
+++ b/codec2/src/phase.c
@@ -71,6 +71,7 @@ int    order;
 
   for(i=0; i<=order; i++)
     Pw[i].real = aks[i];
+
   four1(&Pw[-1].imag,FFT_DEC,-1);
 
   /* Sample magnitude and phase at harmonics */
diff --git a/codec2/src/quantise.c b/codec2/src/quantise.c
index a8ea3d5a..633812a9 100644
--- a/codec2/src/quantise.c
+++ b/codec2/src/quantise.c
@@ -176,7 +176,6 @@ void load_cb(char *filename, float *cb, int k, int m)
     for(i=0; i<m; i++) {
 	scan_line(ftext, &cb[k*lines++], k);
     }
-    printf("%d lines\n",lines);
 
     fclose(ftext);
 }
@@ -244,17 +243,16 @@ float lpc_model_amplitudes(
   float  Sn[],			/* Input frame of speech samples */
   MODEL *model,			/* sinusoidal model parameters */
   int    order,                 /* LPC model order */
-  int    lsp_quantisation       /* optional LSP quantisation if non-zero */
+  int    lsp_quantisation,      /* optional LSP quantisation if non-zero */
+  float  ak[]                   /* output aks */
 )
 {
   float Wn[AW_ENC];
   float R[MAX_ORDER+1];
-  float ak[MAX_ORDER+1];
   float E;
   int   i;
   float snr;	
   float lsp[MAX_ORDER];
-  float lsp_[MAX_ORDER];  /* quantised LSPs */
   int   roots;            /* number of LSP roots found */
   int   index;
   float se;
diff --git a/codec2/src/quantise.h b/codec2/src/quantise.h
index 762c0dfc..46ce79db 100644
--- a/codec2/src/quantise.h
+++ b/codec2/src/quantise.h
@@ -30,7 +30,7 @@
 #include "sine.h"
 
 void quantise_init();
-float lpc_model_amplitudes(float Sn[], MODEL *model, int order, int lsp);
+float lpc_model_amplitudes(float Sn[], MODEL *model, int order,int lsp,float ak[]);
 void aks_to_M2(float ak[], int order, MODEL *model, float E, float *snr);
 
 #endif
diff --git a/codec2/src/refine.c b/codec2/src/refine.c
index b050d8a3..3eb62984 100644
--- a/codec2/src/refine.c
+++ b/codec2/src/refine.c
@@ -50,12 +50,18 @@ void dft_speech()
     Sw[i].imag = 0.0;
   }
 
-  /* centre analysis window on time axis */
+  /* Centre analysis window on time axis, we need to arrange input
+     to FFT this way to make FFT phases correct */
   
-  for(i=0; i<AW_ENC/2; i++)
-    Sw[i].real = Sn[i+AW_ENC/2]*w[i+AW_ENC/2];
-  for(i=FFT_ENC-AW_ENC/2; i<FFT_ENC; i++)
-    Sw[i].real = Sn[i-FFT_ENC+AW_ENC/2]*w[i-FFT_ENC+AW_ENC/2];
+  /* move 2nd half to start of FFT input vector */
+
+  for(i=0; i<NW/2; i++)
+    Sw[i].real = Sn[i+M/2]*w[i+AW_ENC/2];
+
+  /* move 1st half to end of FFT input vector */
+
+  for(i=0; i<NW/2; i++)
+    Sw[FFT_ENC-NW/2+i].real = Sn[i+M/2-NW/2]*w[i+AW_ENC/2-NW/2];
 
   four1(&Sw[-1].imag,FFT_ENC,-1);
 }
diff --git a/codec2/src/sinedec.c b/codec2/src/sinedec.c
index a01a8e8d..54083330 100644
--- a/codec2/src/sinedec.c
+++ b/codec2/src/sinedec.c
@@ -80,6 +80,8 @@ int main(int argc, char *argv[])
 
   int lpc_model, order;
   int lsp, lsp_quantiser;
+  float ak[LPC_MAX_ORDER+1];
+
   int dump;
   
   int phase, phase_model;
@@ -152,9 +154,10 @@ int main(int argc, char *argv[])
   /* phase_model 0: zero phase
      phase_model 1: 1st order polynomial */
   phase = switch_present("--phase",argc,argv);
-  if (phase)
+  if (phase) {
       phase_model = atoi(argv[phase+1]);
-  assert((phase_model == 0) || (phase_model == 1));
+      assert((phase_model == 0) || (phase_model == 1));
+  }
 
   /* Initialise ------------------------------------------------------------*/
 
@@ -183,28 +186,39 @@ int main(int argc, char *argv[])
 
     dump_model(&model);
 
+    /* optional LPC model amplitudes */
+
+    if (lpc_model) {
+	snr = lpc_model_amplitudes(Sn, &model, order, lsp_quantiser, ak);
+	sum_snr += snr;
+        dump_quantised_model(&model);
+    }
+
     /* optional phase modelling */
 
     if (phase) {
 	float Wn[AW_ENC];		/* windowed speech samples */
 	float Rk[PHASE_LPC_ORD+1];	/* autocorrelation coeffs  */
-	float aks[PHASE_LPC_ORD+1];
         COMP  H[MAX_AMP];               /* LPC freq domain samples */
 	int   i_min;
 	COMP  min_Am;
 	
 	dump_phase(&model.phi[0]);
 
-	/* Determine LPC model using time domain LPC.  A little
-	   further down the development track optionally LPCs from lpc
-	   modelling/LSP quant for phase modelling */
+	if (!lpc_model) {
+	    /* Determine LPC model using time domain LPC if we don't have
+	       any LPCs yet */
 
-	for(i=0; i<AW_ENC; i++)
-	    Wn[i] = Sn[i]*w[i];
-	autocorrelate(Wn,Rk,AW_ENC,PHASE_LPC_ORD);
-	levinson_durbin(Rk,aks,PHASE_LPC_ORD);
+	    for(i=0; i<AW_ENC; i++)
+		Wn[i] = Sn[i]*w[i];
+	    autocorrelate(Wn,Rk,AW_ENC,PHASE_LPC_ORD);
+	    levinson_durbin(Rk,ak,PHASE_LPC_ORD);
+	}
+	else
+	    assert(order == PHASE_LPC_ORD);
 
-	snr = phase_model_first_order(aks, H, &i_min, &min_Am);
+	snr = phase_model_first_order(ak, H, &i_min, &min_Am);
+	snr = 5;
 	if (phase_model == 0)
 	    phase_synth_zero_order(snr, H, &prev_Wo, &ex_phase);
 	if (phase_model == 1)
@@ -213,15 +227,6 @@ int main(int argc, char *argv[])
         dump_phase_(&model.phi[0]);
     }
 
-    /* optional LPC model amplitudes */
-
-    if (lpc_model) {
-	snr = lpc_model_amplitudes(Sn, &model, order, lsp_quantiser);
-	sum_snr += snr;
-        dump_quantised_model(&model);
-    }
-
-
     /* Synthesise speech */
 
     if (fout != NULL) {
diff --git a/codec2/src/sinenc.c b/codec2/src/sinenc.c
index be266bd5..70f098da 100644
--- a/codec2/src/sinenc.c
+++ b/codec2/src/sinenc.c
@@ -41,11 +41,11 @@
 \*---------------------------------------------------------------------------*/
 
 int switch_present(sw,argc,argv)
-register char sw[];     /* switch in string form */
-register int argc;      /* number of command line arguments */
-register char *argv[];  /* array of command line arguments in string form */
+  char sw[];     /* switch in string form */
+  int argc;      /* number of command line arguments */
+  char *argv[];  /* array of command line arguments in string form */
 {
-  register int i;       /* loop variable */
+  int i;       /* loop variable */
 
   for(i=1; i<argc; i++)
     if (!strcmp(sw,argv[i]))
@@ -112,8 +112,7 @@ int main(int argc, char *argv[])
     fref = NULL;
 
   init_encoder();
-  Nw = 220;
-  make_window(Nw);
+  make_window(NW);
 
   /* Main loop ------------------------------------------------------------*/
 
@@ -123,20 +122,14 @@ int main(int argc, char *argv[])
 
     /* Update input speech buffers */
 
-    for(i=0; i<N+AW_ENC/2; i++)
+    for(i=0; i<M-N; i++)
       Sn[i] = Sn[i+N];
     for(i=0; i<N; i++)
-      Sn[i+N+AW_ENC/2] = buf[i];
+      Sn[i+M-N] = buf[i];
 
     /* Estimate pitch */
 
-    if (frames > 2) {
-      fscanf(fp,"%f\n",&pitch);
-      if (pitch > P_MAX) pitch = P_MAX;
-      if (pitch < P_MIN) pitch = P_MIN;
-    }
-    else
-      pitch = P_MIN;
+    fscanf(fp,"%f\n",&pitch);
 
     /* construct analysis window */
 
-- 
2.25.1