LPC_Analysis.cpp

00001 //-------------------------------------------------------
00002 // Class for spectrum analysis using linear prediction
00003 // Copyright (c) 2015 MIKAMI, Naoki,  2015/10/26
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         xData_ = new float[nData];      // Data to be analyzed
00016         an_ = new float[order];         // Linear-predictive coefficients
00017         xFft_ = new float[nFft];        // Input for FFT
00018         yFft_ = new Complex[nFft/2+1];  // Output of FFT
00019         normY_ = new float[nFft/2+1];   // Powerspectrum
00020     }
00021 
00022     LpcAnalyzer::~LpcAnalyzer()
00023     {
00024         delete[] xData_;
00025         delete[] an_;
00026         delete[] xFft_;
00027         delete[] yFft_;
00028         delete[] normY_;
00029     }
00030 
00031     void LpcAnalyzer::Execute(float xn[], float db[])
00032     {
00033         // Differencing
00034         for (int n=0; n<N_DATA_-1; n++)
00035             xData_[n] = xn[n+1] - xn[n];
00036 
00037         hm_.Execute(xData_, xData_); // Windowing
00038         float em;
00039         lp_.Execute(xData_, an_, em);
00040 
00041         // To spectrum
00042         xFft_[0] = 1.0f;
00043         for (int n=0; n<ORDER_; n++) xFft_[n+1] = -an_[n];
00044         for (int n=ORDER_+1; n<N_FFT_; n++) xFft_[n] = 0.0f;
00045         fft_.Execute(xFft_, yFft_); // Execute FFT
00046 
00047         // Squared magnitude
00048         for (int n=0; n<=N_FFT_/2; n++)
00049             normY_[n] = 1.0f/(Sqr(yFft_[n].real()) + Sqr(yFft_[n].imag()));
00050 
00051         // Search maximum
00052         float max = 0;
00053         for (int n=0; n<=N_FFT_/2; n++)
00054             max = (max > normY_[n]) ? max : normY_[n];
00055         float invMax = 1.0f/max;
00056 
00057         // Translate to dB
00058         for (int n=0; n<=N_FFT_/2; n++)
00059             db[n] = 10.0f*log10f(invMax*normY_[n]);
00060     }
00061 }
00062