From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Fri, 1 Apr 2011 00:20:56 +0000 (+0000)
Subject: support for new amp VQ scheme, aplogies if this breaks anything
X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=7b1ee1e7ef20bd9f358fb716cb3497d436aa27a6;p=freetel-svn-tracking.git

support for new amp VQ scheme, aplogies if this breaks anything

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

diff --git a/codec2/src/c2sim.c b/codec2/src/c2sim.c
index 8c2657d3..bb49c789 100644
--- a/codec2/src/c2sim.c
+++ b/codec2/src/c2sim.c
@@ -125,6 +125,11 @@ int main(int argc, char *argv[])
 
   void *nlp_states;
   float hpf_states[2];
+  int   resample;
+  float AresdB_prev[MAX_AMP];
+
+  for(i=0; i<MAX_AMP; i++)
+      AresdB_prev[i] = 0.0; 
 
   for(i=0; i<M; i++)
       Sn[i] = 1.0;
@@ -234,6 +239,9 @@ int main(int argc, char *argv[])
 
   decimate = switch_present("--dec",argc,argv);
 
+  arg = switch_present("--resample",argc,argv);
+  resample = atoi(argv[arg+1]);
+
   /* Initialise ------------------------------------------------------------*/
 
   make_analysis_window(w,W);
@@ -252,8 +260,10 @@ int main(int argc, char *argv[])
 
     for(i=0; i<M-N; i++)
       Sn[i] = Sn[i+N];
-    for(i=0; i<N; i++)
-	Sn[i+M-N] = buf[i];
+    for(i=0; i<N; i++) {
+	//Sn[i+M-N] = hpf((float)buf[i], hpf_states);
+	Sn[i+M-N] = (float)buf[i];
+    }
  
     /* Estimate pitch */
 
@@ -348,6 +358,17 @@ int main(int argc, char *argv[])
 #endif
     }
 
+    /* optional resampling of model amplitudes */
+
+    printf("frames=%d\n", frames);
+    if (resample) {
+	snr = resample_amp_nl(&model, resample, AresdB_prev);
+	sum_snr += snr;
+#ifdef DUMP
+        dump_quantised_model(&model);
+#endif
+    }
+
     /* option decimation to 20ms rate, which enables interpolation
        routine to synthesise in between frame */
   
@@ -414,7 +435,7 @@ int main(int argc, char *argv[])
   if (fout != NULL)
     fclose(fout);
 
-  if (lpc_model)
+  if (lpc_model || resample)
       printf("SNR av = %5.2f dB\n", sum_snr/frames);
 
 #ifdef DUMP
diff --git a/codec2/src/codec2.c b/codec2/src/codec2.c
index 4b669153..92708ee3 100644
--- a/codec2/src/codec2.c
+++ b/codec2/src/codec2.c
@@ -42,31 +42,7 @@
 #include "interp.h"
 #include "postfilter.h"
 #include "codec2.h"
-
-typedef struct {
-    float  w[M];	       /* time domain hamming window                */
-    COMP   W[FFT_ENC];	       /* DFT of w[]                                */
-    float  Pn[2*N];	       /* trapezoidal synthesis window              */
-    float  Sn[M];              /* input speech                              */
-    float  hpf_states[2];      /* high pass filter states                   */
-    void  *nlp;                /* pitch predictor states                    */
-    float  Sn_[2*N];	       /* synthesised output speech                 */
-    float  ex_phase;           /* excitation model phase track              */
-    float  bg_est;             /* background noise estimate for post filter */
-    float  prev_Wo;            /* previous frame's pitch estimate           */
-    MODEL  prev_model;         /* previous frame's model parameters         */
-    float  prev_lsps[LPC_ORD]; /* previous frame's LSPs                     */
-    float  prev_energy;        /* previous frame's LPC energy               */
-} CODEC2;
-
-/*---------------------------------------------------------------------------*\
-                                                       
-                             FUNCTION HEADERS
-
-\*---------------------------------------------------------------------------*/
-
-void analyse_one_frame(CODEC2 *c2, MODEL *model, short speech[]);
-void synthesise_one_frame(CODEC2 *c2, short speech[], MODEL *model,float ak[]);
+#include "codec2_internal.h"
 
 /*---------------------------------------------------------------------------*\
                                                        
diff --git a/codec2/src/codec2_internal.h b/codec2/src/codec2_internal.h
new file mode 100644
index 00000000..3943ac29
--- /dev/null
+++ b/codec2/src/codec2_internal.h
@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+
+  FILE........: codec2_internal.h
+  AUTHOR......: David Rowe
+  DATE CREATED: 22 March 2011
+
+  Some internal structures and states broken out here as they are useful for
+  testing and development.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2011 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 <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef __CODEC2_INTERNAL__
+#define  __CODEC2_INTERNAL__
+
+/*---------------------------------------------------------------------------*\
+                                                       
+                             STATES
+
+\*---------------------------------------------------------------------------*/
+
+typedef struct {
+    float  w[M];	       /* time domain hamming window                */
+    COMP   W[FFT_ENC];	       /* DFT of w[]                                */
+    float  Pn[2*N];	       /* trapezoidal synthesis window              */
+    float  Sn[M];              /* input speech                              */
+    float  hpf_states[2];      /* high pass filter states                   */
+    void  *nlp;                /* pitch predictor states                    */
+    float  Sn_[2*N];	       /* synthesised output speech                 */
+    float  ex_phase;           /* excitation model phase track              */
+    float  bg_est;             /* background noise estimate for post filter */
+    float  prev_Wo;            /* previous frame's pitch estimate           */
+    MODEL  prev_model;         /* previous frame's model parameters         */
+    float  prev_lsps[LPC_ORD]; /* previous frame's LSPs                     */
+    float  prev_energy;        /* previous frame's LPC energy               */
+} CODEC2;
+
+/*---------------------------------------------------------------------------*\
+                                                       
+                             FUNCTION HEADERS
+
+\*---------------------------------------------------------------------------*/
+
+void analyse_one_frame(CODEC2 *c2, MODEL *model, short speech[]);
+void synthesise_one_frame(CODEC2 *c2, short speech[], MODEL *model,float ak[]);
+
+#endif
diff --git a/codec2/src/dump.c b/codec2/src/dump.c
index 2215bf27..73a378e2 100644
--- a/codec2/src/dump.c
+++ b/codec2/src/dump.c
@@ -54,6 +54,7 @@ static FILE *fak = NULL;
 static FILE *fbg = NULL;
 static FILE *fE = NULL;
 static FILE *frk = NULL;
+static FILE *fres = NULL;
 
 static char  prefix[MAX_STR];
 
@@ -101,6 +102,8 @@ void dump_off(){
 	fclose(fE);
     if (frk != NULL)
 	fclose(frk);
+    if (fres != NULL)
+	fclose(fres);
 }
 
 void dump_Sn(float Sn[]) {
@@ -221,6 +224,26 @@ void dump_quantised_model(MODEL *model) {
     fprintf(fqmodel,"\n");    
 }
 
+void dump_resample(float w[], float A[], int n) {
+    int l;
+    char s[MAX_STR];
+
+    if (!dumpon) return;
+
+    if (fres == NULL) {
+	sprintf(s,"%s_res.txt", prefix);
+	fres = fopen(s, "wt");
+	assert(fres != NULL);
+    }
+
+    fprintf(fres,"%d\t",n);
+    for(l=0; l<n; l++)
+	fprintf(fres,"%f\t",w[l]);
+    for(l=0; l<n; l++)
+	fprintf(fres,"%f\t",A[l]);
+    fprintf(fres,"\n");    
+}
+
 void dump_phase(float phase[], int L) {
     int l;
     char s[MAX_STR];
diff --git a/codec2/src/dump.h b/codec2/src/dump.h
index 2c7de189..eeddd340 100644
--- a/codec2/src/dump.h
+++ b/codec2/src/dump.h
@@ -44,6 +44,7 @@ void dump_Pw(COMP Pw[]);
 void dump_lsp(float lsp[]);
 void dump_ak(float ak[], int order);
 void dump_E(float E);
+void dump_resample(float w[], float A[], int n);
 
 /* phase modelling */
 
diff --git a/codec2/src/interp.c b/codec2/src/interp.c
index dbdb5bf5..10332055 100644
--- a/codec2/src/interp.c
+++ b/codec2/src/interp.c
@@ -28,11 +28,14 @@
 #include <assert.h>
 #include <math.h>
 #include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
 
 #include "defines.h"
 #include "interp.h"
 #include "lsp.h"
 #include "quantise.h"
+#include "dump.h"
 
 float sample_log_amp(MODEL *model, float w);
 
@@ -110,8 +113,9 @@ float sample_log_amp(MODEL *model, float w)
 
     assert(w > 0.0); assert (w <= PI);
 
-    m = floor(w/model->Wo + 0.5);
-    f = (w - m*model->Wo)/w;
+    m = 0;
+    while ((m+1)*model->Wo < w) m++;
+    f = (w - m*model->Wo)/model->Wo;
     assert(f <= 1.0);
 
     if (m < 1) {
@@ -123,11 +127,287 @@ float sample_log_amp(MODEL *model, float w)
     else {
 	log_amp = (1.0-f)*log10(model->A[m] + 1E-6) + 
                   f*log10(model->A[m+1] + 1E-6);
+	//printf("m=%d A[m] %f A[m+1] %f x %f %f %f\n", m, model->A[m], 
+	//       model->A[m+1], pow(10.0, log_amp),
+	//       (1-f), f);
     }
 
     return log_amp;
 }
 
+/*---------------------------------------------------------------------------*\
+
+  FUNCTION....: sample_log_amp_quad()
+  AUTHOR......: David Rowe			      
+  DATE CREATED: 9 March 2011 
+        
+  Samples the amplitude envelope at an arbitrary frequency w.  Uses
+  quadratic interpolation in the log domain to sample between harmonic
+  amplitudes.
+  
+  y(x) = ax*x + bx + c
+
+  We assume three points are x=-1, x=0, x=1, which we map to m-1,m,m+1
+
+  c = y(0)
+  b = (y(1) - y(-1))/2
+  a = y(-1) + b - y(0)
+
+\*---------------------------------------------------------------------------*/
+
+float sample_log_amp_quad(MODEL *model, float w)
+{
+    int   m;
+    float a,b,c,x, log_amp;
+
+    assert(w > 0.0); assert (w <= PI);
+
+    m = floor(w/model->Wo + 0.5);
+    if (m < 2) m = 2;
+    if (m > (model->L-1)) m = model->L-1;
+    c = log10(model->A[m]+1E-6);
+    b = (log10(model->A[m+1]+1E-6) - log10(model->A[m-1]+1E-6))/2.0;
+    a = log10(model->A[m-1]+1E-6) + b - c;
+    x = (w - m*model->Wo)/model->Wo;
+
+    log_amp = a*x*x + b*x + c;
+    //printf("m=%d A[m-1] %f A[m] %f A[m+1] %f w %f x %f log_amp %f\n", m,
+    //	   model->A[m-1], 
+    //	   model->A[m], model->A[m+1], w, x, pow(10.0, log_amp));
+    return log_amp;
+}
+
+/*---------------------------------------------------------------------------*\
+
+  FUNCTION....: sample_log_amp_quad_nl()
+  AUTHOR......: David Rowe			      
+  DATE CREATED: 10 March 2011
+        
+  Samples the amplitude envelope at an arbitrary frequency w.  Uses
+  quadratic interpolation in the log domain to sample between harmonic
+  amplitudes.  This version can handle non-linear steps along a freq
+  axis defined by arbitrary steps.
+  
+  y(x) = ax*x + bx + c
+
+  We assume three points are (x_1,y_1), (0,y0) and (x1,y1).
+
+\*---------------------------------------------------------------------------*/
+
+float sample_log_amp_quad_nl(
+			     float w[],     /* frequency points            */
+			     float A[],     /* for these amplitude samples */
+			     int np,        /* number of frequency points  */
+			     float w_sample /* frequency of new samples    */
+)
+{
+    int   m,i;
+    float a,b,c,x, log_amp, best_dist;
+    float x_1, x1;
+    float y_1, y0, y1;
+
+    //printf("w_sample  %f\n", w_sample);
+    assert(w_sample >= 0.0); assert (w_sample <= 1.1*PI);
+
+    /* find closest point to centre quadratic interpolator */
+
+    best_dist = 1E32;
+    for (i=0; i<np; i++)
+	if (fabs(w[i] - w_sample) < best_dist) {
+	    best_dist = fabs(w[i] - w_sample);
+	    m = i; 
+	}
+    
+    /* stay one point away from edge of array */
+
+    if (m < 1) m = 1;
+    if (m > (np-2)) m = np - 2;
+
+    /* find polynomial coeffs */
+
+    x_1 = w[m-1]- w[m]; x1 = w[m+1] - w[m];
+    y_1 = log10(A[m-1]+1E-6);
+    y0  = log10(A[m]+1E-6);
+    y1  = log10(A[m+1]+1E-6);
+
+    c = y0;
+    a = (y_1*x1 - y1*x_1 + c*x_1 - c*x1)/(x_1*x_1*x1 - x1*x1*x_1);
+    b = (y1 -a*x1*x1 - c)/x1;
+    x = w_sample - w[m];
+    
+    //printf("%f   %f  %f\n", w[0], w[1], w[2]);
+    //printf("%f %f  %f %f  %f %f\n", x_1, y_1, 0.0, y0, x1, y1);
+    log_amp = a*x*x + b*x + c;
+    //printf("a %f  b %f  c %f\n", a, b, c);
+    //printf("m=%d A[m-1] %f A[m] %f A[m+1] %f w_sample %f w[m] %f x %f log_amp %f\n", m,
+    //	   A[m-1], 
+    //	   A[m], A[m+1], w_sample, w[m], x, log_amp);
+    //exit(0);
+    return log_amp;
+}
+
+#define M_MAX 40
+
+float fres[] = {100,   200,  300,  400,  500,  600,  700,  800,  900, 1000,
+		1200, 1400, 1600, 1850, 2100, 2350, 2600, 2900, 3400, 3800};
+
+/*---------------------------------------------------------------------------*\
+
+  FUNCTION....: resample_amp_nl()
+  AUTHOR......: David Rowe			      
+  DATE CREATED: 7 March 2011
+        
+  Converts the current model with L {Am} samples spaced Wo apart to 
+  RES_POINTS samples spaced Wo/RES_POINTS apart.  Then subtracts
+  from the previous frames samples to get the delta.
+
+\*---------------------------------------------------------------------------*/
+
+void resample_amp_fixed(MODEL *model, 
+			float w[], float A[],
+			float wres[], float Ares[],
+			float AresdB_prev[], 
+			float AresdB[], 
+			float deltat[])
+{
+    int   i;
+
+    for(i=1; i<=model->L; i++) {
+	w[i-1] = i*model->Wo;
+	A[i-1] = model->A[i];
+    }
+
+    for(i=0; i<RES_POINTS; i++) {
+	wres[i] = fres[i]*PI/4000.0;
+    }
+    
+    for(i=0; i<RES_POINTS; i++) {
+	Ares[i] = pow(10.0,sample_log_amp_quad_nl(w, A, model->L, wres[i]));
+    }
+
+    /* work out delta T vector for this frame */
+
+    for(i=0; i<RES_POINTS; i++) {
+	AresdB[i] = 20.0*log10(Ares[i]);
+	deltat[i] = AresdB[i] - AresdB_prev[i];
+    }
+
+}
+
+/*---------------------------------------------------------------------------*\
+
+  FUNCTION....: resample_amp_nl()
+  AUTHOR......: David Rowe			      
+  DATE CREATED: 7 March 2011
+        
+  Converts the current model with L {Am} samples spaced Wo apart to M
+  samples spaced Wo/M apart.  Then converts back to L {Am} samples.
+  used to prototype constant rate Amplitude encoding ideas.
+  
+  Returns the SNR in dB.
+
+\*---------------------------------------------------------------------------*/
+
+float resample_amp_nl(MODEL *model, int m, float AresdB_prev[])
+{
+    int   i;
+    float w[MAX_AMP], A[MAX_AMP];
+    float wres[MAX_AMP], Ares[MAX_AMP], AresdB[MAX_AMP];
+    float signal, noise, snr;
+    float new_A;
+    float deltat[MAX_AMP], deltat_q[MAX_AMP], AresdB_q[MAX_AMP];
+
+    resample_amp_fixed(model, w, A, wres, Ares, AresdB_prev, AresdB, deltat);
+
+    /* quantise delta T vector */
+
+    for(i=0; i<RES_POINTS; i++) {
+	noise = 3.0*(1.0 - 2.0*rand()/RAND_MAX);
+	//noise = 0.0;
+	deltat_q[i] = deltat[i] + noise;
+    }
+
+    /* recover Ares vector */
+
+    for(i=0; i<RES_POINTS; i++) {
+	AresdB_q[i] = AresdB_prev[i] + deltat_q[i];
+	Ares[i] = pow(10.0, AresdB_q[i]/20.0);
+	//printf("%d %f %f\n", i, AresdB[i], AresdB_q[i]);
+    }
+
+    /* update memory based on version at decoder */
+
+    for(i=0; i<RES_POINTS; i++) {
+	AresdB_prev[i] = AresdB_q[i];
+    }
+
+
+    dump_resample(wres,Ares,M_MAX);
+
+    signal = noise = 0.0;
+    
+    for(i=1; i<model->L; i++) {
+	new_A = pow(10.0,sample_log_amp_quad_nl(wres, Ares, RES_POINTS, model->Wo*i));
+	signal += pow(model->A[i], 2.0);
+	noise  += pow(model->A[i] - new_A, 2.0);
+	//printf("%f %f\n", model->A[i], new_A);
+	model->A[i] = new_A;
+    }
+
+    snr = 10.0*log10(signal/noise);
+    printf("snr = %3.2f\n", snr);
+    //exit(0);
+    return snr;
+}
+
+/*---------------------------------------------------------------------------*\
+
+  FUNCTION....: resample_amp()
+  AUTHOR......: David Rowe			      
+  DATE CREATED: 10 March 2011
+        
+  Converts the current model with L {Am} samples spaced Wo apart to M
+  samples with a non-linear spacing.  Then converts back to L {Am}
+  samples.  used to prototype constant rate Amplitude encoding ideas.
+  
+  Returns the SNR in dB.
+
+\*---------------------------------------------------------------------------*/
+
+float resample_amp(MODEL *model, int m)
+{
+    int   i;
+    MODEL model_m;
+    float new_A, signal, noise, snr, log_amp_dB;
+    float n_db = 0.0;
+
+    model_m.Wo = PI/(float)m;
+    model_m.L = PI/model_m.Wo;
+
+    for(i=1; i<=model_m.L; i++) {
+	log_amp_dB    = 20.0*sample_log_amp_quad(model, i*model_m.Wo);
+	log_amp_dB   += n_db*(1.0 - 2.0*rand()/RAND_MAX);
+	model_m.A[i]  = pow(10,log_amp_dB/20.0);
+    }
+
+    //dump_resample(&model_m);
+
+    signal = noise = 0.0;
+    
+    for(i=1; i<model->L/4; i++) {
+	new_A = pow(10,sample_log_amp_quad(&model_m, i*model->Wo));
+	signal += pow(model->A[i], 2.0);
+	noise  += pow(model->A[i] - new_A, 2.0);
+	//printf("%f %f\n", model->A[i], new_A);
+	model->A[i] = new_A;
+    }
+
+    snr = 10.0*log10(signal/noise);
+    //printf("snr = %3.2f\n", snr);
+    //exit(0);
+    return snr;
+}
+
 /*---------------------------------------------------------------------------*\
 
   FUNCTION....: interp_lsp()	     
diff --git a/codec2/src/interp.h b/codec2/src/interp.h
index 3dd5f837..d41eac3f 100644
--- a/codec2/src/interp.h
+++ b/codec2/src/interp.h
@@ -28,10 +28,14 @@
 #ifndef __INTERP__
 #define __INTERP__
 
+#define RES_POINTS 20
+
 void interpolate(MODEL *interp, MODEL *prev, MODEL *next);
 void interpolate_lsp(MODEL *interp, MODEL *prev, MODEL *next, 
 		     float *prev_lsps, float  prev_e,
 		     float *next_lsps, float  next_e,
 		     float *ak_interp);
+float resample_amp(MODEL *model, int m);
+float resample_amp_nl(MODEL *model, int m, float Ares_prev[]);
 
 #endif
diff --git a/codec2/src/lpc.c b/codec2/src/lpc.c
index e38bbd54..ba801137 100644
--- a/codec2/src/lpc.c
+++ b/codec2/src/lpc.c
@@ -78,6 +78,33 @@ void autocorrelate(
   }
 }
 
+/*---------------------------------------------------------------------------*\
+                                                                           
+  autocorrelate_freq()                                                          
+                                                                          
+  Finds the first P autocorrelation values from an array of frequency domain
+  power samples.                                                            
+                                                                          
+\*---------------------------------------------------------------------------*/
+
+void autocorrelate_freq(
+  float Pw[],	/* Nsam frequency domain power spectrum samples    */
+  float  w[],	/* frequency of each sample in Pw[]                */
+  float  R[],	/* array of order+1 autocorrelation coefficients   */
+  int Nsam,	/* number of windowed samples to use               */
+  int order	/* order of LPC analysis                           */
+)
+{
+  int i,j;	/* loop variables */
+
+  for(j=0; j<order+1; j++) {
+    R[j] = 0.0;
+    for(i=0; i<Nsam; i++)
+	R[j] += Pw[i]*cos(j*w[i]);
+  }
+  R[j] /= Nsam;
+}
+
 /*---------------------------------------------------------------------------*\
                                                                             
   levinson_durbin()                                                        
diff --git a/codec2/src/lpc.h b/codec2/src/lpc.h
index dc0cc40d..ead05e1b 100644
--- a/codec2/src/lpc.h
+++ b/codec2/src/lpc.h
@@ -32,6 +32,7 @@
 
 void hanning_window(float Sn[],	float Wn[], int Nsam);
 void autocorrelate(float Sn[], float Rn[], int Nsam, int order);
+void autocorrelate_freq(float Pw[], float w[], float  R[], int Nsam, int order);
 void levinson_durbin(float R[],	float lpcs[], int order);
 void inverse_filter(float Sn[], float a[], int Nsam, float res[], int order);
 void synthesis_filter(float res[], float a[], int Nsam,	int order, float Sn_[]);
diff --git a/codec2/src/phase.c b/codec2/src/phase.c
index d1a88fa3..0e1a14a6 100644
--- a/codec2/src/phase.c
+++ b/codec2/src/phase.c
@@ -221,42 +221,42 @@ void phase_synth_zero_order(
   r = TWO_PI/GLOTTAL_FFT_SIZE;
 
   for(m=1; m<=model->L; m++) {
-
-    /* generate excitation */
-
-    if (model->voiced) {
-	/* I think adding a little jitter helps improve low pitch
-	   males like hts1a. This moves the onset of each harmonic
-	   over at +/- 0.25 of a sample.
-	*/
-        jitter = 0.25*(1.0 - 2.0*rand()/RAND_MAX);
-        b = floor(m*model->Wo/r + 0.5);
-	if (b > ((GLOTTAL_FFT_SIZE/2)-1)) {
-		b = (GLOTTAL_FFT_SIZE/2)-1;
-	}
-	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 {
-
-	/* When a few samples were tested I found that LPC filter
-	   phase is not needed in the unvoiced case, but no harm in
-	   keeping it.
-        */
-	float phi = TWO_PI*(float)rand()/RAND_MAX;
-        Ex[m].real = cos(phi);
-	Ex[m].imag = sin(phi);
-    }
-
-    /* filter using LPC filter */
-
-    A_[m].real = H[m].real*Ex[m].real - H[m].imag*Ex[m].imag;
-    A_[m].imag = H[m].imag*Ex[m].real + H[m].real*Ex[m].imag;
-
-    /* modify sinusoidal phase */
+	  
+      /* generate excitation */
+
+      if (model->voiced) {
+	  /* I think adding a little jitter helps improve low pitch
+	     males like hts1a. This moves the onset of each harmonic
+	     over at +/- 0.25 of a sample.
+	  */
+	  jitter = 0.25*(1.0 - 2.0*rand()/RAND_MAX);
+	  b = floor(m*model->Wo/r + 0.5);
+	  if (b > ((GLOTTAL_FFT_SIZE/2)-1)) {
+	      b = (GLOTTAL_FFT_SIZE/2)-1;
+	  }
+	  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 {
+
+	  /* When a few samples were tested I found that LPC filter
+	     phase is not needed in the unvoiced case, but no harm in
+	     keeping it.
+	  */
+	  float phi = TWO_PI*(float)rand()/RAND_MAX;
+	  Ex[m].real = cos(phi);
+	  Ex[m].imag = sin(phi);
+      }
+
+      /* filter using LPC filter */
+
+      A_[m].real = H[m].real*Ex[m].real - H[m].imag*Ex[m].imag;
+      A_[m].imag = H[m].imag*Ex[m].real + H[m].real*Ex[m].imag;
+
+      /* modify sinusoidal phase */
    
-    new_phi = atan2(A_[m].imag, A_[m].real+1E-12);
-    model->phi[m] = new_phi;
+      new_phi = atan2(A_[m].imag, A_[m].real+1E-12);
+      model->phi[m] = new_phi;
   }
 
 }