Exhale LPC implementation, in case it actually proves useful

exhale_lpc.patch

diff --git a/src/app/CMakeLists.txt b/src/app/CMakeLists.txt
index 1ebe9b8..aaf821e 100644
--- a/src/app/CMakeLists.txt
+++ b/src/app/CMakeLists.txt
@@ -20,6 +20,8 @@ add_executable(exhaleApp
     basicWavReader.h
     basicWavReader.cpp
     exhaleAppPch.cpp
+    lpc.h
+    lpc.cpp
     ${PROJECT_SOURCE_DIR}/include/exhaleDecl.h
     ${PROJECT_SOURCE_DIR}/include/version.h)
 
diff --git a/src/app/exhaleAppPch.cpp b/src/app/exhaleAppPch.cpp
index 21d9ea0..6a90176 100644
--- a/src/app/exhaleAppPch.cpp
+++ b/src/app/exhaleAppPch.cpp
@@ -9,6 +9,7 @@
  */
 
 #include "exhaleAppPch.h"
+#include "lpc.h"
 
 // ISO/IEC 23003-3 USAC Table 67
 static const unsigned supportedSamplingRates[16] = {
@@ -22,16 +23,38 @@ void eaExtrapolate (int32_t* const pcmBuffer, const uint16_t pcmOffset, // start
 {
   const int32_t delta = (fadeIn ? -1 : 1) * numChannels;
   const uint16_t size = (fadeIn ? pcmOffset : frameSize - pcmOffset);
+  const int order=16;
 
   if ((pcmOffset == 0 && fadeIn) || (pcmOffset >= frameSize) || !pcmBuffer) return;
 
+  float lpc[order];
+  float work[frameSize];
+
   for (uint16_t ch = 0; ch < numChannels; ch++)
   {
-    int32_t* chPcmBuf = pcmBuffer + ch + (pcmOffset - (fadeIn ? 0 : 1)) * numChannels;
-    int32_t  result32 = (pcmOffset == 0 ? 0 : *chPcmBuf << 8); // input is known to be 24-bit PCM
-    const int32_t s32 = result32 / size;
+    if (frameSize - size > order * 2)
+    {
+      int32_t* chPcmBuf = pcmBuffer + ch + (fadeIn ? frameSize : -1) * numChannels;
+      for (uint16_t i = 0; i < frameSize; i++) work[i] = *(chPcmBuf += delta);
+
+      vorbis_lpc_from_data(work, lpc, frameSize - size, order);
+
+      vorbis_lpc_predict(lpc, work + frameSize - size - order,
+                         order,
+                         work + frameSize - size,
+                         size);
+
+      chPcmBuf = pcmBuffer + ch + (fadeIn ? frameSize : -1) * numChannels;
+      for (uint16_t i = 0; i < frameSize; i++) *(chPcmBuf += delta) = (int32_t)work[i];
+    }
+    else
+    {
+      int32_t* chPcmBuf = pcmBuffer + ch + (pcmOffset - (fadeIn ? 0 : 1)) * numChannels;
+      int32_t  result32 = (pcmOffset == 0 ? 0 : *chPcmBuf << 8); // input is known to be 24-bit PCM
+      const int32_t s32 = result32 / size;
 
-    for (uint16_t i = size; i > 0; i--) *(chPcmBuf += delta) = ((result32 -= s32) + 128) >> 8;
+      for (uint16_t i = size; i > 0; i--) *(chPcmBuf += delta) = ((result32 -= s32) + 128) >> 8;
+    }
   }
 }
 
diff --git a/src/app/exhaleApp_vs2012.vcxproj b/src/app/exhaleApp_vs2012.vcxproj
index 78ec24c..c9a6d69 100644
--- a/src/app/exhaleApp_vs2012.vcxproj
+++ b/src/app/exhaleApp_vs2012.vcxproj
@@ -160,6 +160,7 @@
     <ClInclude Include="basicWavReader.h" />
     <ClInclude Include="exhaleAppPch.h" />
     <ClInclude Include="loudnessEstim.h" />
+    <ClInclude Include="lpc.h" />
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="basicMP4Writer.cpp" />
@@ -172,6 +173,7 @@
       <PrecompiledHeader Condition="'$(Configuration)|$(Platform)'=='Release|x64'">Create</PrecompiledHeader>
     </ClCompile>
     <ClCompile Include="loudnessEstim.cpp" />
+    <ClCompile Include="lpc.cpp" />
   </ItemGroup>
   <ItemGroup>
     <ResourceCompile Include="exhaleApp.rc" />
diff --git a/src/app/exhaleApp_vs2012.vcxproj.filters b/src/app/exhaleApp_vs2012.vcxproj.filters
index bb3601c..9ad5aa0 100644
--- a/src/app/exhaleApp_vs2012.vcxproj.filters
+++ b/src/app/exhaleApp_vs2012.vcxproj.filters
@@ -33,6 +33,9 @@
     <ClInclude Include="loudnessEstim.h">
       <Filter>Header Files</Filter>
     </ClInclude>
+    <ClInclude Include="lpc.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="basicMP4Writer.cpp">
@@ -50,6 +53,9 @@
     <ClCompile Include="loudnessEstim.cpp">
       <Filter>Source Files</Filter>
     </ClCompile>
+    <ClCompile Include="lpc.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ResourceCompile Include="exhaleApp.rc">
diff --git a/src/app/lpc.cpp b/src/app/lpc.cpp
new file mode 100644
index 0000000..571614a
--- /dev/null
+++ b/src/app/lpc.cpp
@@ -0,0 +1,155 @@
+/********************************************************************
+ *                                                                  *
+ * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
+ *                                                                  *
+ * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009             *
+ * by the Xiph.Org Foundation https://xiph.org/                     *
+ *                                                                  *
+ ********************************************************************
+
+  function: LPC low level routines
+
+ ********************************************************************/
+
+/* Some of these routines (autocorrelator, LPC coefficient estimator)
+   are derived from code written by Jutta Degener and Carsten Bormann;
+   thus we include their copyright below.  The entirety of this file
+   is freely redistributable on the condition that both of these
+   copyright notices are preserved without modification.  */
+
+/* Preserved Copyright: *********************************************/
+
+/* Copyright 1992, 1993, 1994 by Jutta Degener and Carsten Bormann,
+Technische Universita"t Berlin
+
+Any use of this software is permitted provided that this notice is not
+removed and that neither the authors nor the Technische Universita"t
+Berlin are deemed to have made any representations as to the
+suitability of this software for any purpose nor are held responsible
+for any defects of this software. THERE IS ABSOLUTELY NO WARRANTY FOR
+THIS SOFTWARE.
+
+As a matter of courtesy, the authors request to be informed about uses
+this software has found, about bugs in this software, and about any
+improvements that may be of general interest.
+
+Berlin, 28.11.1994
+Jutta Degener
+Carsten Bormann
+
+*********************************************************************/
+
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "lpc.h"
+
+/* Autocorrelation LPC coeff generation algorithm invented by
+   N. Levinson in 1947, modified by J. Durbin in 1959. */
+
+/* Input : n elements of time doamin data
+   Output: m lpc coefficients, excitation energy */
+
+float vorbis_lpc_from_data(float *data,float *lpci,int n,int m){
+  double aut[m+1];
+  double lpc[m];
+  double error;
+  double epsilon;
+  int i,j;
+
+  /* autocorrelation, p+1 lag coefficients */
+  j=m+1;
+  while(j--){
+    double d=0; /* double needed for accumulator depth */
+    for(i=j;i<n;i++)d+=(double)data[i]*data[i-j];
+    aut[j]=d;
+  }
+
+  /* Generate lpc coefficients from autocorr values */
+
+  /* set our noise floor to about -100dB */
+  error=aut[0] * (1. + 1e-10);
+  epsilon=1e-9*aut[0]+1e-10;
+
+  for(i=0;i<m;i++){
+    double r= -aut[i+1];
+
+    if(error<epsilon){
+      memset(lpc+i,0,(m-i)*sizeof(*lpc));
+      goto done;
+    }
+
+    /* Sum up this iteration's reflection coefficient; note that in
+       Vorbis we don't save it.  If anyone wants to recycle this code
+       and needs reflection coefficients, save the results of 'r' from
+       each iteration. */
+
+    for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
+    r/=error;
+
+    /* Update LPC coefficients and total error */
+
+    lpc[i]=r;
+    for(j=0;j<i/2;j++){
+      double tmp=lpc[j];
+
+      lpc[j]+=r*lpc[i-1-j];
+      lpc[i-1-j]+=r*tmp;
+    }
+    if(i&1)lpc[j]+=lpc[j]*r;
+
+    error*=1.-r*r;
+
+  }
+
+ done:
+
+  /* slightly damp the filter */
+  {
+    double g = .99;
+    double damp = g;
+    for(j=0;j<m;j++){
+      lpc[j]*=damp;
+      damp*=g;
+    }
+  }
+
+  for(j=0;j<m;j++)lpci[j]=(float)lpc[j];
+
+  /* we need the error value to know how big an impulse to hit the
+     filter with later */
+
+  return error;
+}
+
+void vorbis_lpc_predict(float *coeff,float *prime,int m,
+                     float *data,long n){
+
+  /* in: coeff[0...m-1] LPC coefficients
+         prime[0...m-1] initial values (allocated size of n+m-1)
+    out: data[0...n-1] data samples */
+
+  long i,j,o,p;
+  float y;
+  float work[m+n];
+
+  if(!prime)
+    for(i=0;i<m;i++)
+      work[i]=0.f;
+  else
+    for(i=0;i<m;i++)
+      work[i]=prime[i];
+
+  for(i=0;i<n;i++){
+    y=0;
+    o=i;
+    p=m;
+    for(j=0;j<m;j++)
+      y-=work[o++]*coeff[--p];
+
+    data[i]=work[o]=y;
+  }
+}
diff --git a/src/app/lpc.h b/src/app/lpc.h
new file mode 100644
index 0000000..ddfdf85
--- /dev/null
+++ b/src/app/lpc.h
@@ -0,0 +1,26 @@
+/********************************************************************
+ *                                                                  *
+ * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
+ *                                                                  *
+ * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2007             *
+ * by the Xiph.Org Foundation https://xiph.org/                     *
+ *                                                                  *
+ ********************************************************************
+
+  function: LPC low level routines
+
+ ********************************************************************/
+
+#ifndef _V_LPC_H_
+#define _V_LPC_H_
+
+/* simple linear scale LPC code */
+extern float vorbis_lpc_from_data(float *data,float *lpc,int n,int m);
+
+extern void vorbis_lpc_predict(float *coeff,float *prime,int m,
+                               float *data,long n);
+
+#endif
diff --git a/src/app/makefile b/src/app/makefile
index 8625a0c..142baad 100644
--- a/src/app/makefile
+++ b/src/app/makefile
@@ -31,6 +31,7 @@ OBJS      = \
 	$(DIR_OBJ)/exhaleApp.o \
 	$(DIR_OBJ)/exhaleAppPch.o \
 	$(DIR_OBJ)/loudnessEstim.o \
+	$(DIR_OBJ)/lpc.o \
 
 # define libraries to link with
 LIBS      = -ldl