LinearPrediction.cpp

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