Diff: SpactrumAnalysisClasses/LPC_Analysis.cpp

Revision:

0:c35b8a23a863

Child:

1:2ef356b500f2

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SpactrumAnalysisClasses/LPC_Analysis.cpp	Mon Oct 26 08:06:57 2015 +0000
@@ -0,0 +1,62 @@
+//-------------------------------------------------------
+// Class for spectrum analysis using linear prediction
+// Copyright (c) 2015 MIKAMI, Naoki,  2015/10/26
+//-------------------------------------------------------
+
+#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)
+    {
+        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[] 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
+
+        // Squared magnitude
+        for (int n=0; n<=N_FFT_/2; n++)
+            normY_[n] = 1.0f/(Sqr(yFft_[n].real()) + Sqr(yFft_[n].imag()));
+
+        // Searce maximum
+        float max = 0;
+        for (int n=0; n<=N_FFT_/2; n++)
+            max = (max > normY_[n]) ? max : normY_[n];
+        float invMax = 1.0f/max;
+
+        // Translate to dB
+        for (int n=0; n<=N_FFT_/2; n++)
+            db[n] = 10.0f*log10f(invMax*normY_[n]);
+    }
+}
+