Diff: MyClasses/LinearPrediction.cpp

Revision:

0:5c237fdcba23

Child:

4:99d4d5ea06a2

diff -r 000000000000 -r 5c237fdcba23 MyClasses/LinearPrediction.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/MyClasses/LinearPrediction.cpp	Wed Dec 09 05:05:00 2015 +0000
@@ -0,0 +1,71 @@
+//-----------------------------------------------------
+//  Class for linear prediction
+//
+//  2015/12/08, Copyright (c) 2015 MIKAMI, Naoki
+//-----------------------------------------------------
+
+#include "LinearPrediction.hpp"
+
+namespace Mikami
+{
+    LinearPred::LinearPred(int nData, int order)
+        : N_DATA_(nData), ORDER_(order),
+          r_(new float[order+1]), k_(new float[order]),
+          am_(new float[order]) {}
+
+    LinearPred::~LinearPred()
+    {
+        delete[] r_;
+        delete[] k_;
+        delete[] am_;
+    }
+
+    // Calculate linear-predictive coefficients
+    bool LinearPred::Execute(const float x[], float a[],
+                             float &em)
+    {
+        AutoCorr(x);
+        return Durbin(a, em);
+    }
+
+    // Calculate auto-correlation
+    void LinearPred::AutoCorr(const float x[])
+    {
+        for (int j=0; j<=ORDER_; j++)
+        {
+            r_[j] = 0.0;
+            for (int n=0; n<N_DATA_-j; n++)
+                r_[j] = r_[j] + x[n]*x[n+j];
+        }
+    }
+
+    // Levinson-Durbin algorithm
+    bool LinearPred::Durbin(float a[], float &em)
+    {
+        // Initialization
+        em = r_[0];
+
+        // Repeat
+        for (int m=0; m<ORDER_; m++)
+        {
+            float w = r_[m+1];
+            for (int j=0; j<=m-1; j++)
+                w = w - r_[m-j]*a[j];
+
+            k_[m] = w/em;
+            em = em*(1 - k_[m]*k_[m]);
+
+            if (em < 0) break;  // Error for negative squared sum of residual
+
+            a[m] = k_[m];
+            for (int j=0; j<=m-1; j++)
+                am_[j] = a[j];
+            for (int j=0; j<=m-1; j++)
+                a[j] = am_[j] - k_[m]*am_[m-j-1];
+        }
+
+        if (em < 0) return false;
+        else        return true;
+    }
+}
+