Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/Filter.cpp
Go to the documentation of this file.
00001 /*
00002 ** Copyright (C) 1998-2010 George Tzanetakis <gtzan@cs.uvic.ca>
00003 **
00004 ** This program is free software; you can redistribute it and/or modify
00005 ** it under the terms of the GNU General Public License as published by
00006 ** the Free Software Foundation; either version 2 of the License, or
00007 ** (at your option) any later version.
00008 **
00009 ** This program is distributed in the hope that it will be useful,
00010 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012 ** GNU General Public License for more details.
00013 **
00014 ** You should have received a copy of the GNU General Public License
00015 ** along with this program; if not, write to the Free Software
00016 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00017 */
00018 
00019 #include "Filter.h"
00020 #include "../common_source.h"
00021 
00022 
00023 using std::ostringstream;
00024 using std::ofstream;
00025 using std::endl;
00026 
00027 
00028 using namespace Marsyas;
00029 
00030 Filter::Filter(mrs_string name)
00031   :
00032   MarSystem("Filter",name),
00033   norder_(2),
00034   dorder_(1),
00035   channels_(1),
00036   order_(2),
00037   fgain_(1.0)
00038 {
00039   ncoeffs_.create(norder_);
00040   dcoeffs_.create(dorder_);
00041   state_.create(channels_,order_-1);
00042 
00043   ncoeffs_(0) = 1.0;
00044   dcoeffs_(0) = 1.0;
00045 
00046   addControls();
00047 }
00048 
00049 Filter::~Filter()
00050 {
00051 
00052 }
00053 
00054 
00055 MarSystem*
00056 Filter::clone() const
00057 {
00058   return new Filter(*this);
00059 }
00060 
00061 void
00062 Filter::addControls()
00063 {
00064   addctrl("mrs_realvec/ncoeffs", ncoeffs_);
00065   addctrl("mrs_realvec/dcoeffs", dcoeffs_);
00066   addctrl("mrs_real/fgain", fgain_);
00067   addctrl("mrs_natural/stateUpdate", mrs_natural(0));
00068   addctrl("mrs_realvec/state", state_);
00069 
00070   setctrlState("mrs_realvec/ncoeffs", true);
00071   setctrlState("mrs_realvec/dcoeffs", true);
00072   setctrlState("mrs_realvec/state", true);
00073 }
00074 
00075 void Filter::myUpdate(MarControlPtr sender)
00076 {
00077   (void) sender;  //suppress warning of unused parameter(s)
00078   MRSDIAG("Filter.cpp - Filter:myUpdate");
00079 
00080   setctrl("mrs_natural/onSamples", getctrl("mrs_natural/inSamples"));
00081   setctrl("mrs_natural/onObservations", getctrl("mrs_natural/inObservations"));
00082   setctrl("mrs_real/osrate", getctrl("mrs_real/israte"));
00083 
00084   ctrl_onObsNames_->setValue("Filter_" + ctrl_inObsNames_->to<mrs_string>() + "," , NOUPDATE);
00085 
00086   if (getctrl("mrs_realvec/ncoeffs")->to<mrs_realvec>().getSize() != norder_)
00087   {
00088     ncoeffs_.create(getctrl("mrs_realvec/ncoeffs")->to<mrs_realvec>().getSize());
00089     norder_ = ncoeffs_.getSize();
00090     order_ = (norder_ > dorder_) ? norder_ : dorder_;
00091     channels_ = getctrl("mrs_natural/inObservations")->to<mrs_natural>();
00092     state_.create(channels_,order_-1);
00093     setctrl("mrs_realvec/state", state_);
00094   }
00095 
00096   if (getctrl("mrs_realvec/dcoeffs")->to<mrs_realvec>().getSize() != dorder_)
00097   {
00098 
00099     dcoeffs_.create(getctrl("mrs_realvec/dcoeffs")->to<mrs_realvec>().getSize());
00100     dorder_ = dcoeffs_.getSize();
00101     order_ = (norder_ > dorder_) ? norder_ : dorder_;
00102     channels_ = getctrl("mrs_natural/inObservations")->to<mrs_natural>();
00103     state_.create(channels_,order_-1);
00104     setctrl("mrs_realvec/state", state_);
00105   }
00106 
00107   if (getctrl("mrs_natural/inObservations")->to<mrs_natural>() != channels_)
00108   {
00109     channels_ = getctrl("mrs_natural/inObservations")->to<mrs_natural>();
00110     state_.create(channels_,order_-1);
00111   }
00112 
00113   ncoeffs_ = getctrl("mrs_realvec/ncoeffs")->to<mrs_realvec>();
00114   dcoeffs_ = getctrl("mrs_realvec/dcoeffs")->to<mrs_realvec>();
00115 
00116   if (getctrl("mrs_natural/stateUpdate")->to<mrs_natural>())
00117     state_ = getctrl("mrs_realvec/state")->to<mrs_realvec>();
00118 
00119   mrs_real d0 = dcoeffs_(0);
00120   if (d0 != 1.0) {
00121     for (mrs_natural i = 0; i < dorder_; ++i) {
00122       dcoeffs_(i) /= d0;
00123     }
00124 
00125     for (mrs_natural i = 0; i < norder_; ++i) {
00126       ncoeffs_(i) /= d0;
00127     }
00128   }
00129 
00130   fgain_ = 1.0f;
00131   setctrl("mrs_real/fgain", 1.0);
00132 }
00133 
00134 void
00135 Filter::write(mrs_string filename)
00136 {
00137   ofstream os(filename.c_str());
00138   os << (*this) << endl;
00139 }
00140 
00141 void
00142 Filter::myProcess(realvec& in, realvec& out)
00143 {
00144   //checkFlow(in,out);
00145 
00146   mrs_natural i,j,c;
00147   mrs_natural size = in.getCols();
00148   mrs_natural stateSize = state_.getCols();
00149   mrs_natural channels = in.getRows();
00150 
00151   mrs_real gain = getctrl("mrs_real/fgain")->to<mrs_real>();
00152 
00153   // State array holds the various delays for the difference equation
00154   // of the filter. Similar implementation as described in the manual
00155   // for MATLAB Signal Processing Toolbox. State corresponds to
00156   // the z, num_coefs to the b and denom_coefs to the a vector respectively
00157   // in_window is the input x(n) and out_window is the output y(n)
00158 
00159   //dcoeffs_/=10;
00160 
00161 
00162   // state_.setval(0);
00163   if (norder_ == dorder_) {
00164     for (c = 0; c < channels; ++c) {
00165       for (i = 0; i < size; ++i) {
00166         out(c,i) = ncoeffs_(0) * in(c,i) + state_(c,0);
00167         for (j = 0; j < stateSize - 1; j++)
00168         {
00169           state_(c,j) = ncoeffs_(j+1) * in(c,i) + state_(c,j+1) - dcoeffs_(j+1) * out(c,i);
00170         }
00171         state_(c,stateSize - 1) = ncoeffs_(order_-1) * in(c,i) - dcoeffs_(order_-1) * out(c,i);
00172       }
00173     }
00174   }
00175   else if (norder_ < dorder_) {
00176     for (c = 0; c < channels; ++c) {
00177       for (i = 0; i < size; ++i) {
00178         out(c,i) = ncoeffs_(0) * in(c,i) + state_(c,0);
00179         for (j = 0; j < norder_ - 1; j++)
00180         {
00181           state_(c,j) = ncoeffs_(j+1) * in(c,i) + state_(c,j+1) - dcoeffs_(j+1) * out(c,i);
00182         }
00183         for (j = norder_ - 1; j < stateSize - 1; j++)
00184         {
00185           state_(c,j) = state_(c,j+1) - dcoeffs_(j+1) * out(c,i);
00186         }
00187         state_(c,stateSize - 1) = -dcoeffs_(order_ - 1) * out(c,i);
00188       }
00189     }
00190   }
00191   else {
00192     for (c = 0; c < channels; ++c) {
00193       for (i = 0; i < size; ++i) {
00194         out(c,i) = ncoeffs_(0) * in(c,i) + state_(c,0);
00195         for (j = 0; j < dorder_ - 1; j++)
00196         {
00197           state_(c,j) = ncoeffs_(j+1) * in(c,i) + state_(c,j+1) - dcoeffs_(j+1) * out(c,i);
00198         }
00199         for (j = dorder_ - 1; j < stateSize - 1; j++)
00200         {
00201           state_(c,j) = ncoeffs_(j+1) * in(c,i) + state_(c,j+1);
00202         }
00203         state_(c,stateSize - 1) = ncoeffs_(order_-1) * in(c,i);
00204       }
00205     }
00206   }
00207   out *= gain;
00208 
00209 
00210   //        MATLAB_PUT(in, "Filter_in");
00211   //        MATLAB_PUT(out, "Filter_out");
00212   //        MATLAB_PUT(ncoeffs_, "ncoeffs_");
00213   //        MATLAB_PUT(dcoeffs_, "dcoeffs_");
00214   //        MATLAB_EVAL("MAT_out = filter(ncoeffs_, dcoeffs_, Filter_in)");
00215   //
00216   //        MATLAB_EVAL("spec_in = abs(fft(Filter_in));");
00217   //        MATLAB_EVAL("spec_out = abs(fft(Filter_out));");
00218   //        MATLAB_EVAL("spec_mat = abs(fft(MAT_out));");
00219   //
00220   //        MATLAB_EVAL("subplot(2,1,1);plot(Filter_in);hold on; plot(Filter_out, 'r'); plot(MAT_out, 'g');hold off");
00221   //        MATLAB_EVAL("subplot(2,1,2);plot(spec_in(1:end/2));hold on; plot(spec_out(1:end/2),'r');plot(spec_mat(1:end/2),'g');hold off;");
00222   //        MATLAB_EVAL("h = abs(fft([1 -.97], length(Filter_in)));");
00223   //        MATLAB_EVAL("hold on; plot(h(1:end/2), 'k'); hold off");
00224   //        //MATLAB_GET("MAT_out", out)
00225   //
00226 }