--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/autocorr.c	Wed Apr 03 05:31:37 2024 +0000
@@ -0,0 +1,99 @@
+/*************************************************************************
+ *
+ *  FUNCTION:  autocorr
+ *
+ *  PURPOSE:   Compute autocorrelations of signal with windowing
+ *
+ *  DESCRIPTION:
+ *       - Windowing of input speech:   s'[n] = s[n] * w[n]
+ *       - Autocorrelations of input speech:
+ *             r[k] = sum_{i=k}^{239} s'[i]*s'[i-k]    k=0,...,10
+ *         The autocorrelations are expressed in normalized double precision
+ *         format.
+ *
+ *************************************************************************/
+
+#include "typedef.h"
+#include "basic_op.h"
+#include "oper_32b.h"
+#include "count.h"
+#include "cnst.h"
+
+Word16 Autocorr (
+    Word16 x[],         /* (i)    : Input signal                    */
+    Word16 m,           /* (i)    : LPC order                       */
+    Word16 r_h[],       /* (o)    : Autocorrelations  (msb)         */
+    Word16 r_l[],       /* (o)    : Autocorrelations  (lsb)         */
+    Word16 wind[]       /* (i)    : window for LPC analysis         */
+)
+{
+    Word16 i, j, norm;
+    Word16 y[L_WINDOW];
+    Word32 sum;
+    Word16 overfl, overfl_shft;
+
+    /* Windowing of signal */
+
+    for (i = 0; i < L_WINDOW; i++)
+    {
+        y[i] = mult_r (x[i], wind[i]); move16 (); 
+    }
+
+    /* Compute r[0] and test for overflow */
+
+    overfl_shft = 0;                   move16 (); 
+
+    do
+    {
+        overfl = 0;                    move16 (); 
+        sum = 0L;                      move32 ();
+
+        for (i = 0; i < L_WINDOW; i++)
+        {
+            sum = L_mac (sum, y[i], y[i]);
+        }
+
+        /* If overflow divide y[] by 4 */
+
+        test (); 
+        if (L_sub (sum, MAX_32) == 0L)
+        {
+            overfl_shft = add (overfl_shft, 4);
+            overfl = 1;                move16 (); /* Set the overflow flag */
+
+            for (i = 0; i < L_WINDOW; i++)
+            {
+                y[i] = shr (y[i], 2);  move16 (); 
+            }
+        }
+        test (); 
+    }
+    while (overfl != 0);
+
+    sum = L_add (sum, 1L);             /* Avoid the case of all zeros */
+
+    /* Normalization of r[0] */
+
+    norm = norm_l (sum);
+    sum = L_shl (sum, norm);
+    L_Extract (sum, &r_h[0], &r_l[0]); /* Put in DPF format (see oper_32b) */
+
+    /* r[1] to r[m] */
+
+    for (i = 1; i <= m; i++)
+    {
+        sum = 0;                       move32 (); 
+
+        for (j = 0; j < L_WINDOW - i; j++)
+        {
+            sum = L_mac (sum, y[j], y[j + i]);
+        }
+
+        sum = L_shl (sum, norm);
+        L_Extract (sum, &r_h[i], &r_l[i]);
+    }
+
+    norm = sub (norm, overfl_shft);
+
+    return norm;
+}