Marsyas
0.6.0-alpha
|
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 }