Merge pull request #83 from drowe67/ss-wreal

Change W[] vector to float as imag always 0.0

Author: drowe67

misc/tnlp.c

@@ -91,7 +91,7 @@ int main(int argc, char *argv[])
     kiss_fft_cfg  fft_fwd_cfg;
     COMP  Sw[FFT_ENC];	        /* DFT of Sn[] */
     float w[m];	                /* time domain hamming window */
-    COMP  W[FFT_ENC];	        /* DFT of w[] */
+    float W[FFT_ENC];	        /* DFT of w[] */
     float pitch_samples;
     int   i;
     float f0, prev_f0;

src/c2sim.c

@@ -404,7 +404,7 @@ int main(int argc, char *argv[])
     codec2_fftr_cfg  fftr_fwd_cfg;
     codec2_fftr_cfg  fftr_inv_cfg;
     float w[m_pitch];	        /* time domain hamming window            */
-    COMP  W[FFT_ENC];	/* DFT of w[]                            */
+    float W[FFT_ENC];	/* DFT of w[]                            */
     MODEL model;
     float Pn[2*N_SAMP];	/* trapezoidal synthesis window          */
     float Sn_[2*N_SAMP];	/* synthesised speech */

src/codec2_internal.h

@@ -42,7 +42,7 @@ struct CODEC2 {
     codec2_fft_cfg  fft_fwd_cfg;           /* forward FFT config                        */
     codec2_fftr_cfg fftr_fwd_cfg;          /* forward real FFT config                   */
     float        *w;	                   /* [m_pitch] time domain hamming window      */
-    COMP          W[FFT_ENC];	           /* DFT of w[]                                */
+    float         W[FFT_ENC];	           /* DFT of w[]                                */
     float        *Pn;	                   /* [2*n_samp] trapezoidal synthesis window   */
     float        *bpf_buf;                 /* buffer for band pass filter               */
     float        *Sn;                      /* [m_pitch] input speech                    */

src/nlp.c

@@ -53,7 +53,6 @@
 #define F0_MAX      500
 #define CNLP        0.3	        /* post processor constant              */
 #define NLP_NTAP 48	        /* Decimation LPF order */
-#undef  POST_PROCESS_MBE        /* choose post processor                */
 
 /* 8 to 16 kHz sample rate conversion */
 
@@ -132,10 +131,6 @@ typedef struct {
     FILE         *f;
 } NLP;
 
-#ifdef POST_PROCESS_MBE
-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);
-#endif
 float post_process_sub_multiples(COMP Fw[],
 				 int pmin, int pmax, float gmax, int gmax_bin,
 				 float *prev_f0);
@@ -258,7 +253,7 @@ float nlp(
   int    n,			/* frames shift (no. new samples in Sn[])             */
   float *pitch,			/* estimated pitch period in samples at current Fs    */
   COMP   Sw[],                  /* Freq domain version of Sn[]                        */
-  COMP   W[],                   /* Freq domain window                                 */
+  float  W[],                   /* Freq domain window                                 */
   float *prev_f0                /* previous pitch f0 in Hz, memory for pitch tracking */
 )
 {
@@ -389,11 +384,7 @@ float nlp(
 
     PROFILE_SAMPLE_AND_LOG(peakpick, magsq, "      peak pick");
 
-    #ifdef POST_PROCESS_MBE
-    best_f0 = post_process_mbe(Fw, pmin, pmax, gmax, Sw, W, prev_f0);
-    #else
     best_f0 = post_process_sub_multiples(Fw, pmin, pmax, gmax, gmax_bin, prev_f0);
-    #endif
 
     PROFILE_SAMPLE_AND_LOG(shiftmem, peakpick,  "      post process");
 
@@ -491,178 +482,6 @@ float post_process_sub_multiples(COMP Fw[],
     return best_f0;
 }
 
-#ifdef POST_PROCESS_MBE
-
-/*---------------------------------------------------------------------------*\
-
-  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;
-  #ifdef DUMP
-  float e_hz[F0_MAX];
-  #endif
-  #if !defined(NDEBUG) || defined(DUMP)
-  int   bin;
-  #endif
-  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 */
-
-  #ifdef DUMP
-  for(i=0; i<F0_MAX; i++)
-      e_hz[i] = -1;
-  #endif
-  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);
-	    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);
-		#if !defined(NDEBUG) || defined(DUMP)
-		bin = floorf(f0); assert((bin > 0) && (bin < F0_MAX));
-		#endif
-		#ifdef DUMP
-                e_hz[bin] = e;
-                #endif
-		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);
-      #if !defined(NDEBUG) || defined(DUMP)
-      bin = floorf(f0); assert((bin > 0) && (bin < F0_MAX));
-      #endif
-      #ifdef DUMP
-      e_hz[bin] = e;
-      #endif
-      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 = floorf((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<L/4; l++) {
-	Am.real = 0.0;
-	Am.imag = 0.0;
-	den = 0.0;
-	al = ceilf((l - 0.5)*Wo*FFT_ENC/TWO_PI);
-	bl = ceilf((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;
-}
-
-#endif
-
 /*---------------------------------------------------------------------------*\
 
   FUNCTION....: fdmdv_16_to_8()

src/nlp.h

@@ -33,6 +33,6 @@
 void *nlp_create(C2CONST *c2const);
 void nlp_destroy(void *nlp_state);
 float nlp(void *nlp_state, float Sn[], int n, 
-	  float *pitch_samples, COMP Sw[], COMP W[], float *prev_f0);
+	  float *pitch_samples, COMP Sw[], float W[], float *prev_f0);
 
 #endif

src/sine.c

@@ -97,11 +97,10 @@ C2CONST c2const_create(int Fs, float framelength_s) {
 
 \*---------------------------------------------------------------------------*/
 
-void make_analysis_window(C2CONST *c2const, codec2_fft_cfg fft_fwd_cfg, float w[], COMP W[])
+void make_analysis_window(C2CONST *c2const, codec2_fft_cfg fft_fwd_cfg, float w[], float W[])
 {
   float m;
   COMP  wshift[FFT_ENC];
-  COMP  temp;
   int   i,j;
   int   m_pitch = c2const->m_pitch;
   int   nw      = c2const->nw;
@@ -156,6 +155,8 @@ void make_analysis_window(C2CONST *c2const, codec2_fft_cfg fft_fwd_cfg, float w[
        nw/2              nw/2
   */
 
+  COMP temp[FFT_ENC];
+
   for(i=0; i<FFT_ENC; i++) {
     wshift[i].real = 0.0;
     wshift[i].imag = 0.0;
@@ -165,7 +166,7 @@ void make_analysis_window(C2CONST *c2const, codec2_fft_cfg fft_fwd_cfg, float w[
   for(i=FFT_ENC-nw/2,j=m_pitch/2-nw/2; i<FFT_ENC; i++,j++)
    wshift[i].real = w[j];
 
-  codec2_fft(fft_fwd_cfg, wshift, W);
+  codec2_fft(fft_fwd_cfg, wshift, temp);
 
   /*
       Re-arrange W[] to be symmetrical about FFT_ENC/2.  Makes later
@@ -192,12 +193,8 @@ void make_analysis_window(C2CONST *c2const, codec2_fft_cfg fft_fwd_cfg, float w[
 
 
   for(i=0; i<FFT_ENC/2; i++) {
-    temp.real = W[i].real;
-    temp.imag = W[i].imag;
-    W[i].real = W[i+FFT_ENC/2].real;
-    W[i].imag = W[i+FFT_ENC/2].imag;
-    W[i+FFT_ENC/2].real = temp.real;
-    W[i+FFT_ENC/2].imag = temp.imag;
+      W[i] = temp[i + FFT_ENC / 2].real;
+      W[i + FFT_ENC / 2] = temp[i].real;
   }
 
 }
@@ -402,39 +399,30 @@ void hs_pitch_refinement(MODEL *model, COMP Sw[], float pmin, float pmax, float
 
 \*---------------------------------------------------------------------------*/
 
-void estimate_amplitudes(MODEL *model, COMP Sw[], COMP W[], int est_phase)
+void estimate_amplitudes(MODEL *model, COMP Sw[], float W[], int est_phase)
 {
   int   i,m;		/* loop variables */
   int   am,bm;		/* bounds of current harmonic */
-  int   b;		/* DFT bin of centre of current harmonic */
   float den;		/* denominator of amplitude expression */
-  float r, one_on_r;	/* number of rads/bin */
-  int   offset;
-  COMP  Am;
 
-  r = TWO_PI/FFT_ENC;
-  one_on_r = 1.0/r;
+  float r = TWO_PI/FFT_ENC;
+  float one_on_r = 1.0/r;
 
   for(m=1; m<=model->L; m++) {
+    /* Estimate ampltude of harmonic */
+
     den = 0.0;
     am = (int)((m - 0.5)*model->Wo*one_on_r + 0.5);
     bm = (int)((m + 0.5)*model->Wo*one_on_r + 0.5);
-    b = (int)(m*model->Wo/r + 0.5);
-
-    /* Estimate ampltude of harmonic */
 
-    den = 0.0;
-    Am.real = Am.imag = 0.0;
-    offset = FFT_ENC/2 - (int)(m*model->Wo*one_on_r + 0.5);
     for(i=am; i<bm; i++) {
       den += Sw[i].real*Sw[i].real + Sw[i].imag*Sw[i].imag;
-      Am.real += Sw[i].real*W[i + offset].real;
-      Am.imag += Sw[i].imag*W[i + offset].real;
     }
 
     model->A[m] = sqrtf(den);
 
     if (est_phase) {
+        int b = (int)(m*model->Wo/r + 0.5); /* DFT bin of centre of current harmonic */
 
         /* Estimate phase of harmonic, this is expensive in CPU for
            embedded devicesso we make it an option */
@@ -459,7 +447,7 @@ float est_voicing_mbe(
                       C2CONST *c2const,
                       MODEL *model,
                       COMP   Sw[],
-                      COMP   W[]
+                      float  W[]
                       )
 {
     int   l,al,bl,m;    /* loop variables */
@@ -497,20 +485,19 @@ float est_voicing_mbe(
 
         offset = FFT_ENC/2 - l*Wo*FFT_ENC/TWO_PI + 0.5;
 	for(m=al; m<bl; m++) {
-	    Am.real += Sw[m].real*W[offset+m].real;
-	    Am.imag += Sw[m].imag*W[offset+m].real;
-	    den += W[offset+m].real*W[offset+m].real;
+	    Am.real += Sw[m].real*W[offset+m];
+	    Am.imag += Sw[m].imag*W[offset+m];
+	    den += W[offset+m]*W[offset+m];
         }
 
         Am.real = Am.real/den;
         Am.imag = Am.imag/den;
 
         /* Determine error between estimated harmonic and original */
 
-// Redundant!        offset = FFT_ENC/2 - l*Wo*FFT_ENC/TWO_PI + 0.5;
         for(m=al; m<bl; m++) {
-	    Ew.real = Sw[m].real - Am.real*W[offset+m].real;
-	    Ew.imag = Sw[m].imag - Am.imag*W[offset+m].real;
+	    Ew.real = Sw[m].real - Am.real*W[offset+m];
+	    Ew.imag = Sw[m].imag - Am.imag*W[offset+m];
 	    error += Ew.real*Ew.real;
 	    error += Ew.imag*Ew.imag;
 	}

src/sine.h

@@ -34,12 +34,12 @@
 
 C2CONST c2const_create(int Fs, float framelength_ms);
 
-void make_analysis_window(C2CONST *c2const, codec2_fft_cfg fft_fwd_cfg, float w[], COMP W[]);
+void make_analysis_window(C2CONST *c2const, codec2_fft_cfg fft_fwd_cfg, float w[], float W[]);
 float hpf(float x, float states[]);
 void dft_speech(C2CONST *c2const, codec2_fft_cfg fft_fwd_cfg, COMP Sw[], float Sn[], float w[]);
 void two_stage_pitch_refinement(C2CONST *c2const, MODEL *model, COMP Sw[]);
-void estimate_amplitudes(MODEL *model, COMP Sw[], COMP W[], int est_phase);
-float est_voicing_mbe(C2CONST *c2const, MODEL *model, COMP Sw[], COMP W[]);
+void estimate_amplitudes(MODEL *model, COMP Sw[], float W[], int est_phase);
+float est_voicing_mbe(C2CONST *c2const, MODEL *model, COMP Sw[], float W[]);
 void make_synthesis_window(C2CONST *c2const, float Pn[]);
 void synthesise(int n_samp, codec2_fftr_cfg fftr_inv_cfg, float Sn_[], MODEL *model, float Pn[], int shift);