Diff: SpactrumAnalysisClasses/LPC_Analysis.cpp

Revision:

0:c5b026c2d07e

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SpactrumAnalysisClasses/LPC_Analysis.cpp	Sun Jul 26 02:48:23 2015 +0000
@@ -0,0 +1,62 @@
+//-------------------------------------------------------
+// Class for spectrum analysis using linear prediction
+// Copyright (c) 2014 MIKAMI, Naoki,  2014/12/30
+//-------------------------------------------------------
+
+#include "LPC_Analysis.hpp"
+
+namespace Mikami
+{
+    LpcAnalyzer::LpcAnalyzer(int nData, int order, int nFft)
+        : N_DATA_(nData), ORDER_(order), N_FFT_(nFft),
+          hm_(nData-1, nData-1), lp_(nData-1, order),
+          fft_(nFft)
+    {
+        pkHolder_ = new PeakHolder[nFft/2+1];
+        for (int n=0; n<=nFft/2; n++)
+            pkHolder_[n].SetCoefs(0.2f);    
+
+        xData_ = new float[nData];      // Data to be analyzed
+        an_ = new float[order];         // Linear-predictive coefficients
+        xFft_ = new float[nFft];        // Input for FFT
+        yFft_ = new Complex[nFft/2+1];  // Output of FFT
+        normY_ = new float[nFft/2+1];   // Powerspectrum
+    }
+
+    LpcAnalyzer::~LpcAnalyzer()
+    {
+        delete[] pkHolder_;
+        delete[] xData_;
+        delete[] an_;
+        delete[] xFft_;
+        delete[] yFft_;
+        delete[] normY_;
+    }
+
+    void LpcAnalyzer::Execute(float xn[], float db[])
+    {
+        // Differencing
+        for (int n=0; n<N_DATA_-1; n++)
+            xData_[n] = xn[n+1] - xn[n];
+
+        hm_.Execute(xData_, xData_); // Windowing
+        float em;
+        lp_.Execute(xData_, an_, em);
+
+        // To spectrum
+        xFft_[0] = 1.0f;
+        for (int n=0; n<ORDER_; n++) xFft_[n+1] = -an_[n];
+        for (int n=ORDER_+1; n<N_FFT_; n++) xFft_[n] = 0.0f;
+        fft_.Execute(xFft_, yFft_); // Execute FFT
+
+        // Smoothing
+        for (int n=0; n<=N_FFT_/2; n++)
+            normY_[n] = pkHolder_[n].Execute(Sqr(yFft_[n].real())
+                      + Sqr(yFft_[n].imag()));
+
+        // Translate to dB
+        float b0Db = 20.0f*log10f(em);  // b0 to dB
+        for (int n=0; n<=N_FFT_/2; n++)
+            db[n] = -10.0f*log10f(normY_[n]) + b0Db + 30;
+    }
+}