LinearPrediction.cpp

00001 //-----------------------------------------------------
00002 // Class for linear prediction
00003 // Copyright (c) 2014 MIKAMI, Naoki,  2014/12/30
00004 //-----------------------------------------------------
00005 
00006 #include "LinearPrediction.hpp"
00007 
00008 namespace Mikami
00009 {
00010     LinearPred::LinearPred(int nData, int order)
00011         : N_DATA_(nData), ORDER_(order)
00012     {
00013         r_ = new float[order+1];    // for auto-correlation
00014         k_ = new float[order];      // for PARCOR coefficients
00015         am_ = new float[order];     // working area
00016     }
00017 
00018     LinearPred::~LinearPred()
00019     {
00020         delete[] r_;
00021         delete[] k_;
00022         delete[] am_;
00023     }
00024 
00025     // Calculate linear-predictive coefficients
00026     bool LinearPred::Execute(const float x[], float a[],
00027                              float &em)
00028     {
00029         AutoCorr(x);
00030         return Durbin(a, em);
00031     }
00032 
00033     // Calculate auto-correlation
00034     void LinearPred::AutoCorr(const float x[])
00035     {
00036         for (int j=0; j<=ORDER_; j++)
00037         {
00038             r_[j] = 0.0;
00039             for (int n=0; n<N_DATA_-j; n++)
00040                 r_[j] = r_[j] + x[n]*x[n+j];
00041         }
00042     }
00043 
00044     // Levinson-Durbin algorithm
00045     bool LinearPred::Durbin(float a[], float &em)
00046     {
00047         // Initialization
00048         em = r_[0];
00049 
00050         // Repeat
00051         for (int m=0; m<ORDER_; m++)
00052         {
00053             float w = r_[m+1];
00054             for (int j=0; j<=m-1; j++)
00055                 w = w - r_[m-j]*a[j];
00056 
00057             k_[m] = w/em;
00058             em = em*(1 - k_[m]*k_[m]);
00059 
00060             if (em < 0) break;      // Error for negative squared sum of residual
00061 
00062             a[m] = k_[m];
00063             for (int j=0; j<=m-1; j++)
00064                 am_[j] = a[j];
00065             for (int j=0; j<=m-1; j++)
00066                 a[j] = am_[j] - k_[m]*am_[m-j-1];
00067         }
00068 
00069         if (em < 0) return false;
00070         else        return true;
00071     }
00072 }