LPC_Analysis.cpp

00001 //-------------------------------------------------------
00002 // Class for spectrum analysis using linear prediction
00003 // Copyright (c) 2014 MIKAMI, Naoki,  2014/12/30
00004 //-------------------------------------------------------
00005 
00006 #include "LPC_Analysis.hpp"
00007 
00008 namespace Mikami
00009 {
00010     LpcAnalyzer::LpcAnalyzer(int nData, int order, int nFft)
00011         : N_DATA_(nData), ORDER_(order), N_FFT_(nFft),
00012           hm_(nData-1, nData-1), lp_(nData-1, order),
00013           fft_(nFft)
00014     {
00015         pkHolder_ = new PeakHolder[nFft/2+1];
00016         for (int n=0; n<=nFft/2; n++)
00017             pkHolder_[n].SetCoefs(0.2f);    
00018 
00019         xData_ = new float[nData];      // Data to be analyzed
00020         an_ = new float[order];         // Linear-predictive coefficients
00021         xFft_ = new float[nFft];        // Input for FFT
00022         yFft_ = new Complex[nFft/2+1];  // Output of FFT
00023         normY_ = new float[nFft/2+1];   // Powerspectrum
00024     }
00025 
00026     LpcAnalyzer::~LpcAnalyzer()
00027     {
00028         delete[] pkHolder_;
00029         delete[] xData_;
00030         delete[] an_;
00031         delete[] xFft_;
00032         delete[] yFft_;
00033         delete[] normY_;
00034     }
00035 
00036     void LpcAnalyzer::Execute(float xn[], float db[])
00037     {
00038         // Differencing
00039         for (int n=0; n<N_DATA_-1; n++)
00040             xData_[n] = xn[n+1] - xn[n];
00041 
00042         hm_.Execute(xData_, xData_); // Windowing
00043         float em;
00044         lp_.Execute(xData_, an_, em);
00045 
00046         // To spectrum
00047         xFft_[0] = 1.0f;
00048         for (int n=0; n<ORDER_; n++) xFft_[n+1] = -an_[n];
00049         for (int n=ORDER_+1; n<N_FFT_; n++) xFft_[n] = 0.0f;
00050         fft_.Execute(xFft_, yFft_); // Execute FFT
00051 
00052         // Smoothing
00053         for (int n=0; n<=N_FFT_/2; n++)
00054             normY_[n] = pkHolder_[n].Execute(Sqr(yFft_[n].real())
00055                       + Sqr(yFft_[n].imag()));
00056 
00057         // Translate to dB
00058         float b0Db = 20.0f*log10f(em);  // b0 to dB
00059         for (int n=0; n<=N_FFT_/2; n++)
00060             db[n] = -10.0f*log10f(normY_[n]) + b0Db + 30;
00061     }
00062 }