Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/LSP.cpp
Go to the documentation of this file.
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 "../common_source.h"
00022 
00023 #include <marsyas/NumericLib.h>
00024 
00025 #include <algorithm>
00026 #include <vector>
00027 
00028 // #define _MATLAB_LSP_
00029 
00030 
00031 using std::ostringstream;
00032 using std::vector;
00033 using std::polar;
00034 
00035 using namespace Marsyas;
00036 
00037 LSP::LSP(mrs_string name):MarSystem("LSP",name)
00038 {
00039   addControls();
00040 }
00041 
00042 LSP::~LSP()
00043 {
00044 }
00045 
00046 MarSystem*
00047 LSP::clone() const
00048 {
00049   return new LSP(*this);
00050 }
00051 
00052 void
00053 LSP::addControls()
00054 {
00055   addctrl("mrs_natural/order", (mrs_natural)10);//read-only control
00056   addctrl("mrs_real/gamma", (mrs_real)1.0);
00057 }
00058 
00059 void
00060 LSP::myUpdate(MarControlPtr sender)
00061 {
00062   (void) sender;  //suppress warning of unused parameter(s)
00063   MRSDIAG("LSP.cpp - LSP:myUpdate");
00064 
00065   order_ = getctrl("mrs_natural/inObservations")->to<mrs_natural>() - 2;
00066   setctrl("mrs_natural/order", order_);//read-only control
00067 
00068   setctrl("mrs_natural/onObservations", order_);
00069   setctrl("mrs_natural/onSamples", getctrl("mrs_natural/inSamples"));
00070   setctrl("mrs_real/osrate", getctrl("mrs_real/israte"));
00071 
00072   //LSP features names
00073   ostringstream oss;
00074   for (mrs_natural i = 0; i < order_; ++i)
00075     oss << "LSP_" << i+1 << ",";
00076   setctrl("mrs_string/onObsNames", oss.str());
00077 }
00078 
00079 void
00080 LSP::myProcess(realvec& in, realvec& out)
00081 {
00082   NumericLib numLib;
00083 
00084   mrs_real gamma = getctrl("mrs_real/gamma")->to<mrs_real>();
00085   vector<mrs_real> ak(order_);
00086 
00087   if( gamma != 1.0)
00088     for(mrs_natural j = 0; j < order_ ; j++)
00089     {
00090       ak[j] = in(j)*pow((mrs_real)gamma, (mrs_real)j+1);//*(-1.0);//apply pole-shifting
00091     }
00092   else
00093     for(mrs_natural j = 0; j < order_ ; j++)
00094     {
00095       ak[j] = in(j);//*(-1.0); //no pole-shifting applied
00096     }
00097 
00098   vector<mrs_complex> P(order_+2);
00099   vector<mrs_complex> Q(order_+2);
00100   vector<mrs_complex> Proots(order_+1);
00101   vector<mrs_complex> Qroots(order_+1);
00102 
00103   P[order_+1] = polar(1.0, 0.0);
00104   Q[order_+1] = polar(1.0, 0.0);
00105   for(mrs_natural k = 0; k < order_; k++)
00106   {
00107     P[order_-k] = polar((double)(ak[k] + ak[order_-1-k]), 0.0);
00108     Q[order_-k] = polar((double)(ak[k] - ak[order_-1-k]), 0.0);
00109   }
00110   P[0] = polar(1.0, 0.0);
00111   Q[0] = polar(-1.0, 0.0);
00112 
00113   if (!numLib.polyRoots(P, false, order_+1, Proots)) {
00114     //P has only real coefs => complexCoefs = false
00115     MRSERR("LSP::myProcess() - numerical error in polynomial root calculation!");
00116   }
00117   if(!numLib.polyRoots(Q, false, order_+1, Qroots)) {
00118     //Q has only real coefs => complexCoefs = false
00119     MRSERR("LSP::myProcess() - numerical error in polynomial root calculation!");
00120   }
00121 
00122   mrs_real phase;
00123   vector<mrs_real> out_vec;
00124   for(mrs_natural k = 0; k <= order_; k++)
00125   {
00126     phase = arg(Proots[k]);
00127     if((phase > 0) && (phase < PI))
00128     {
00129       out_vec.push_back(phase);
00130     }
00131   }
00132   for(mrs_natural k = 0; k <= order_; k++)
00133   {
00134     phase = arg(Qroots[k]);
00135     if((phase > 0) && (phase < PI))
00136     {
00137       out_vec.push_back(phase);
00138     }
00139   }
00140   sort(out_vec.begin(), out_vec.end()); //sorts LSP freqs into ascending order
00141 
00142   //output sorted LSP frequencies
00143   for(mrs_natural i = 0; i < order_; ++i)
00144     out(i) = out_vec[i];
00145 
00146 #ifdef _MATLAB_LSP_
00147   MATLAB_PUT(order_, "LSP_order");
00148   MATLAB_PUT(in, "LSP_in");
00149   MATLAB_PUT(P, "LSP_P");
00150   MATLAB_PUT(Q, "LSP_Q");
00151   MATLAB_PUT(Proots, "LSP_Proots");
00152   MATLAB_PUT(Qroots, "LSP_Qroots");
00153   MATLAB_PUT(out_vec, "LSP_out1");
00154   MATLAB_PUT(out, "LSP_out2");
00155   MATLAB_EVAL("LSP_test(LSP_order, LSP_in, LSP_P, LSP_Q, LSP_Proots, LSP_Qroots, LSP_out1, LSP_out2);");
00156 #endif
00157 }