Marsyas  0.2
/home/gperciva/src/marsyas/src/marsyas/LSP.cpp
00001 
00002 /*
00003 ** Copyright (C) 1998-2010 George Tzanetakis <gtzan@cs.uvic.ca>
00004 **  
00005 ** This program is free software; you can redistribute it and/or modify
00006 ** it under the terms of the GNU General Public License as published by
00007 ** the Free Software Foundation; either version 2 of the License, or
00008 ** (at your option) any later version.
00009 ** 
00010 ** This program is distributed in the hope that it will be useful,
00011 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
00012 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013 ** GNU General Public License for more details.
00014 ** 
00015 ** You should have received a copy of the GNU General Public License
00016 ** along with this program; if not, write to the Free Software 
00017 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00018 */
00019 
00020 #include "LSP.h"
00021 #include "NumericLib.h"
00022 
00023 #include <algorithm>
00024 #include <vector>
00025 
00026 // #define _MATLAB_LSP_
00027 
00028  
00029 using std::ostringstream;
00030 using std::vector;
00031 using std::polar;
00032 
00033 using namespace Marsyas;
00034 
00035 LSP::LSP(mrs_string name):MarSystem("LSP",name)
00036 {
00037     addControls();
00038 }
00039 
00040 LSP::~LSP()
00041 {
00042 }
00043 
00044 MarSystem* 
00045 LSP::clone() const 
00046 {
00047     return new LSP(*this);
00048 }
00049 
00050 void 
00051 LSP::addControls()
00052 {
00053     addctrl("mrs_natural/order", (mrs_natural)10);//read-only control
00054     addctrl("mrs_real/gamma", (mrs_real)1.0);
00055 }
00056 
00057 void
00058 LSP::myUpdate(MarControlPtr sender)
00059 { 
00060     (void) sender;
00061     MRSDIAG("LSP.cpp - LSP:myUpdate");
00062 
00063     order_ = getctrl("mrs_natural/inObservations")->to<mrs_natural>() - 2;
00064     setctrl("mrs_natural/order", order_);//read-only control
00065 
00066     setctrl("mrs_natural/onObservations", order_);
00067     setctrl("mrs_natural/onSamples", getctrl("mrs_natural/inSamples"));
00068     setctrl("mrs_real/osrate", getctrl("mrs_real/israte"));
00069 
00070     //LSP features names
00071     ostringstream oss;
00072     for (mrs_natural i = 0; i < order_; ++i)
00073         oss << "LSP_" << i+1 << ",";
00074     setctrl("mrs_string/onObsNames", oss.str());
00075 }
00076 
00077 void 
00078 LSP::myProcess(realvec& in, realvec& out)
00079 {
00080     NumericLib numLib;
00081     
00082     mrs_real gamma = getctrl("mrs_real/gamma")->to<mrs_real>();
00083     vector<mrs_real> ak(order_);
00084 
00085     if( gamma != 1.0)
00086         for(mrs_natural j = 0; j < order_ ; j++)
00087         {
00088             ak[j] = in(j)*pow((mrs_real)gamma, (mrs_real)j+1);//*(-1.0);//apply pole-shifting
00089         }
00090     else
00091         for(mrs_natural j = 0; j < order_ ; j++)
00092         {
00093             ak[j] = in(j);//*(-1.0); //no pole-shifting applied
00094         }
00095 
00096     vector<mrs_complex> P(order_+2);
00097     vector<mrs_complex> Q(order_+2);
00098     vector<mrs_complex> Proots(order_+1); 
00099     vector<mrs_complex> Qroots(order_+1); 
00100 
00101     P[order_+1] = polar(1.0, 0.0);
00102     Q[order_+1] = polar(1.0, 0.0);
00103     for(mrs_natural k = 0; k < order_; k++)
00104     {
00105         P[order_-k] = polar((double)(ak[k] + ak[order_-1-k]), 0.0);
00106         Q[order_-k] = polar((double)(ak[k] - ak[order_-1-k]), 0.0);
00107     }
00108     P[0] = polar(1.0, 0.0);
00109     Q[0] = polar(-1.0, 0.0);
00110 
00111     if (!numLib.polyRoots(P, false, order_+1, Proots))//P has only real coefs => complexCoefs = false
00112         MRSERR("LSP::myProcess() - numerical error in polynomial root calculation!");
00113     if(!numLib.polyRoots(Q, false, order_+1, Qroots))//Q has only real coefs => complexCoefs = false
00114         MRSERR("LSP::myProcess() - numerical error in polynomial root calculation!");
00115         
00116 
00117     mrs_real phase;
00118     vector<mrs_real> out_vec;
00119     for(mrs_natural k = 0; k <= order_; k++)
00120     {
00121         phase = arg(Proots[k]);
00122         if((phase > 0) && (phase < PI))
00123         {
00124             out_vec.push_back(phase);
00125         }
00126     }
00127     for(mrs_natural k = 0; k <= order_; k++)
00128     {
00129         phase = arg(Qroots[k]);
00130         if((phase > 0) && (phase < PI))
00131         {
00132             out_vec.push_back(phase);
00133         }
00134     }
00135     sort(out_vec.begin(), out_vec.end()); //sorts LSP freqs into ascending order
00136 
00137     //output sorted LSP frequencies
00138     for(mrs_natural i = 0; i < order_; ++i)
00139         out(i) = out_vec[i];
00140 
00141 #ifdef _MATLAB_LSP_
00142     MATLAB_PUT(order_, "LSP_order");
00143     MATLAB_PUT(in, "LSP_in");
00144     MATLAB_PUT(P, "LSP_P");
00145     MATLAB_PUT(Q, "LSP_Q");
00146     MATLAB_PUT(Proots, "LSP_Proots");
00147     MATLAB_PUT(Qroots, "LSP_Qroots");
00148     MATLAB_PUT(out_vec, "LSP_out1");
00149     MATLAB_PUT(out, "LSP_out2");
00150     MATLAB_EVAL("LSP_test(LSP_order, LSP_in, LSP_P, LSP_Q, LSP_Proots, LSP_Qroots, LSP_out1, LSP_out2);");
00151 #endif
00152 }