Marsyas
0.6.0-alpha
|
00001 /* 00002 ** Copyright (C) 1998-2011 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 "MFCC.h" 00020 00021 00022 using std::ostringstream; 00023 00024 00025 using namespace Marsyas; 00026 00027 MFCC::MFCC(mrs_string name):MarSystem("MFCC",name) 00028 { 00029 addControls(); 00030 pfftSize_ = 0; 00031 psamplingRate_ = 0; 00032 mfcc_offsets_ = NULL; 00033 pcepstralCoefs_ = 0; 00034 cepstralCoefs_ = MFCC::cepstralCoefs_default; 00035 } 00036 00037 MFCC::MFCC(const MFCC& a) : MarSystem(a) 00038 { 00039 ctrl_cepstralCoefs_ = getctrl("mrs_natural/coefficients"); 00040 pfftSize_ = 0; 00041 psamplingRate_ = 0; 00042 mfcc_offsets_ = NULL; 00043 pcepstralCoefs_ = 0; 00044 cepstralCoefs_ = MFCC::cepstralCoefs_default; 00045 } 00046 00047 MFCC::~MFCC() 00048 { 00049 if (mfcc_offsets_!=NULL) 00050 delete [] mfcc_offsets_; 00051 } 00052 00053 MarSystem* 00054 MFCC::clone() const 00055 { 00056 return new MFCC(*this); 00057 } 00058 00059 00060 void 00061 MFCC::addControls() { 00063 addControl("mrs_natural/coefficients", MFCC::cepstralCoefs_default, ctrl_cepstralCoefs_); 00064 setControlState("mrs_natural/coefficients", true); 00065 } 00066 00067 void 00068 MFCC::myUpdate(MarControlPtr sender) 00069 { 00070 (void) sender; //suppress warning of unused parameter(s) 00071 00072 // Get the number of cepstral coefficients from the control 00073 cepstralCoefs_ = ctrl_cepstralCoefs_->to<mrs_natural>(); 00074 00075 ctrl_onSamples_->setValue((mrs_natural)1, NOUPDATE); 00076 ctrl_onObservations_->setValue((mrs_natural)cepstralCoefs_, NOUPDATE); 00077 ctrl_osrate_->setValue(ctrl_israte_, NOUPDATE); 00078 00079 // Initialize frequency boundaries for filters 00080 mrs_natural i,j; 00081 fftSize_ = 2 * (ctrl_inObservations_->to<mrs_natural>()-1); //PowerSpectrum outputs N/2+1 "magnitude" spectral points! 00082 if (fftSize_ == 0) return; // skip unnecessary updates 00083 00084 samplingRate_ = (mrs_natural) (ctrl_israte_->to<mrs_real>() * fftSize_); 00085 00086 // Set the observation names: use the first item of the 00087 // input observation names and prefix it with "MFCCxx" 00088 mrs_string inObsName = stringSplit(ctrl_inObsNames_->to<mrs_string>(), ",")[0]; 00089 ostringstream oss; 00090 for (i=0; i < cepstralCoefs_; ++i) 00091 { 00092 oss << "MFCC" << i << "_" << inObsName << ","; 00093 } 00094 ctrl_onObsNames_->setValue(oss.str(), NOUPDATE); 00095 00096 if ((pfftSize_ != fftSize_) || (psamplingRate_ != samplingRate_) || (pcepstralCoefs_ != cepstralCoefs_)) 00097 { 00098 00099 freqs_.create(42); 00100 lowestFrequency_ = 133.3333f; 00101 linearFilters_ = 13; 00102 linearSpacing_ = 66.66666f; 00103 logFilters_ = 27; 00104 logSpacing_ = 1.0711703f; 00105 00106 totalFilters_ = linearFilters_ + logFilters_; 00107 lower_.create(totalFilters_); 00108 center_.create(totalFilters_); 00109 upper_.create(totalFilters_); 00110 triangle_heights_.create(totalFilters_); 00111 00112 // Linear filter boundaries 00113 for (i=0; i< linearFilters_; ++i) 00114 { 00115 freqs_(i) = lowestFrequency_ + i * linearSpacing_; 00116 } 00117 00118 // Logarithmic filter boundaries 00119 mrs_real first_log = freqs_(linearFilters_-1); 00120 for (i=1; i<=logFilters_+2; ++i) 00121 { 00122 freqs_(linearFilters_-1+i) = first_log * pow(logSpacing_, (mrs_real)i); 00123 } 00124 00125 // Triangles information 00126 for (i=0; i<totalFilters_; ++i) 00127 { 00128 lower_(i) = freqs_(i); 00129 center_(i) = freqs_(i+1); 00130 upper_(i) = freqs_(i+2); 00131 triangle_heights_(i) = (mrs_real)(2.0 / (upper_(i) - lower_(i))); 00132 } 00133 00134 fftFreqs_.stretch(fftSize_); 00135 00136 for (i=0; i< fftSize_; ++i) 00137 { 00138 fftFreqs_(i) = (float)i / (float)fftSize_ * (float)samplingRate_; 00139 } 00140 00141 mfccFilterWeights_.create(totalFilters_, fftSize_); 00142 mfccDCT_.create(cepstralCoefs_, totalFilters_); 00143 00144 mrs_natural chan; 00145 00146 // NEIL's filter weight speedup 00147 if (pfftSize_!=fftSize_) 00148 { 00149 if (mfcc_offsets_!=NULL) 00150 { 00151 delete [] mfcc_offsets_; 00152 } 00153 mfcc_offsets_ = new int[totalFilters_*fftSize_*2]; 00154 } 00155 // Initialize mfccFilterWeights 00156 for (chan = 0; chan < totalFilters_; chan++) 00157 { 00158 // NEIL's filter weight speedup 00159 int len=-1; 00160 int pos=0; 00161 for (i=0; i< fftSize_; ++i) 00162 { 00163 if ((fftFreqs_(i) > lower_(chan))&& (fftFreqs_(i) <= center_(chan))) 00164 { 00165 mfccFilterWeights_(chan, i) = triangle_heights_(chan) * 00166 ((fftFreqs_(i) - lower_(chan))/(center_(chan) - lower_(chan))); 00167 // NEIL's filter weight speedup 00168 if (len==-1) 00169 { 00170 pos=i; 00171 } 00172 len=i; 00173 } 00174 if ((fftFreqs_(i) > center_(chan)) && (fftFreqs_(i) <= upper_(chan))) 00175 { 00176 mfccFilterWeights_(chan, i) = triangle_heights_(chan) * 00177 ((upper_(chan) - fftFreqs_(i))/(upper_(chan) - center_(chan))); 00178 // NEIL's filter weight speedup 00179 if (len==-1) 00180 { 00181 pos=i; 00182 } 00183 len=i; 00184 } 00185 } 00186 // NEIL's filter weight speedup 00187 mfcc_offsets_[chan] = pos; 00188 mfcc_offsets_[chan+totalFilters_] = len; 00189 } 00190 00191 // Initialize MFCC_DCT 00192 mrs_real scale_fac = (mrs_real)(1.0/ sqrt((mrs_real)(totalFilters_/2))); 00193 for (j = 0; j<cepstralCoefs_; j++) 00194 { 00195 for (i=0; i< totalFilters_; ++i) 00196 { 00197 mfccDCT_(j, i) = scale_fac * cos(j * (2*i +1) * PI/2/totalFilters_); 00198 if (j == 0) 00199 { 00200 mfccDCT_(j,i) *= (mrs_real)(sqrt(2.0)/2.0); 00201 } 00202 } 00203 } 00204 } 00205 00206 00207 pfftSize_ = fftSize_; 00208 psamplingRate_ = samplingRate_; 00209 00210 fmagnitude_.stretch(ctrl_inObservations_->to<mrs_natural>() * 2); 00211 earMagnitude_.stretch(totalFilters_); 00212 } 00213 00214 void 00215 MFCC::myProcess(realvec& in, realvec& out) 00216 { 00217 mrs_natural i,k,o; 00218 00219 // mirror the spectrum 00220 for (o=0; o < inObservations_; o++) 00221 { 00222 fmagnitude_(o) = in(o,0); 00223 } 00224 00225 for (o=0; o< inObservations_; o++) 00226 { 00227 fmagnitude_(o + inObservations_) = fmagnitude_(inObservations_ - o -1); 00228 } 00229 00230 mrs_real sum =0.0; 00231 // Calculate the filterbank responce 00232 for (i=0; i<totalFilters_; ++i) 00233 { 00234 sum = 0.0; 00235 // NEIL's filter weight speedup 00236 for (k=mfcc_offsets_[i]; k<=mfcc_offsets_[i+totalFilters_]; k++) 00237 { 00238 sum += (mfccFilterWeights_(i, k) * fmagnitude_(k)); 00239 } 00240 if (sum != 0.0) 00241 { 00242 earMagnitude_(i) = log10(sum); 00243 } 00244 else 00245 { 00246 earMagnitude_(i) = 0.0; 00247 } 00248 } 00249 /* The way it used to be : NEIL 00250 for (i=0; i<totalFilters_; ++i) { 00251 sum = 0.0; 00252 for (k=0; k<fftSize_; k++) { 00253 sum += (mfccFilterWeights_(i, k) * fmagnitude_(k)); 00254 } 00255 if (sum != 0.0) 00256 earMagnitude_(i) = log10(sum); 00257 else 00258 earMagnitude_(i) = 0.0; 00259 } 00260 */ 00261 00262 // Take the DCT 00263 for (o=0; o < cepstralCoefs_; o++) 00264 { 00265 sum =0.0; 00266 for (k=0; k < totalFilters_; k++) 00267 { 00268 sum += (mfccDCT_(o,k) * earMagnitude_(k)); 00269 } 00270 out(o,0) = sum; 00271 } 00272 }