Marsyas
0.6.0-alpha
|
00001 /* 00002 ** Copyright (C) 1998-2006 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 "../common_source.h" 00020 #include "TimeFreqPeakConnectivity.h" 00021 #include <algorithm> 00022 00023 using std::max; 00024 using namespace Marsyas; 00025 00026 //#define MATLAB_DBG_OUT 00027 //#define SANITY_CHECK 00028 00029 static const mrs_real costInit = 1e30; 00030 static const mrs_natural listLenInit = 16; 00031 class DoubleListEntries 00032 { 00033 public: 00034 DoubleListEntries () : maxListLength(0),currListLength(0) 00035 { 00036 list = 0; 00037 AllocMem (listLenInit); 00038 }; 00039 ~DoubleListEntries () 00040 { 00041 FreeMem (); 00042 }; 00043 mrs_natural GetNumIndices (mrs_natural row, mrs_natural col) 00044 { 00045 mrs_natural count = 0; 00046 for (mrs_natural i = 0; i < currListLength; i++) 00047 { 00048 if (row == list[i][kRow]) 00049 if (col == list[i][kCol]) 00050 count++; 00051 } 00052 return count; 00053 }; 00054 mrs_natural GetIndex (mrs_natural row, mrs_natural col, mrs_natural index) 00055 { 00056 for (mrs_natural i = 0; i < currListLength; i++) 00057 { 00058 if (row == list[i][kRow]) { 00059 if (col == list[i][kCol]) { 00060 if (index == 0) { 00061 return list[i][kIdx]; 00062 } else { 00063 index--; 00064 } 00065 } 00066 } 00067 } 00068 MRSASSERT(true); 00069 return -1; 00070 }; 00071 mrs_bool IsInList (mrs_natural row, mrs_natural col, mrs_natural index) 00072 { 00073 for (mrs_natural i = 0; i < currListLength; i++) 00074 { 00075 if (row == list[i][kRow]) 00076 if (col == list[i][kCol]) 00077 if (index == list[i][kIdx]) 00078 return true; 00079 } 00080 return false; 00081 }; 00082 void AddIndex (mrs_natural row, mrs_natural col, mrs_natural index) 00083 { 00084 if (IsInList (row,col,index)) 00085 return; 00086 if (currListLength == maxListLength -1) 00087 AllocMem (maxListLength<<1); 00088 list[currListLength][kRow] = row; 00089 list[currListLength][kCol] = col; 00090 list[currListLength][kIdx] = index; 00091 currListLength++; 00092 }; 00093 void Reset () 00094 { 00095 currListLength = 0; 00096 }; 00097 private: 00098 void AllocMem(mrs_natural numNewListItems) 00099 { 00100 mrs_natural i; 00101 MRSASSERT(numNewListItems >= maxListLength); 00102 mrs_natural **tmpList = new mrs_natural* [numNewListItems]; 00103 for (i = 0; i < maxListLength; i++) 00104 tmpList[i] = list[i]; 00105 for (i = maxListLength; i < numNewListItems; i++) 00106 tmpList[i] = new mrs_natural [kNumListEntries]; 00107 00108 delete [] list; 00109 list = tmpList; 00110 maxListLength = numNewListItems; 00111 return; 00112 }; 00113 void FreeMem () 00114 { 00115 for (mrs_natural i = 0; i < maxListLength; i++) 00116 delete [] list[i]; 00117 delete [] list; 00118 }; 00119 00120 enum ListEntries_t 00121 { 00122 kRow, 00123 kCol, 00124 kIdx, 00125 00126 kNumListEntries 00127 }; 00128 mrs_natural **list; 00129 mrs_natural maxListLength, 00130 currListLength; 00131 }; 00132 00133 00134 static inline mrs_natural MyMax (mrs_natural a , mrs_natural b) 00135 { 00136 return (a >= b) ? a : b; 00137 } 00138 static inline mrs_natural MyMin (mrs_natural a , mrs_natural b) 00139 { 00140 return (a <= b) ? a : b; 00141 00142 } 00143 00144 TimeFreqPeakConnectivity::TimeFreqPeakConnectivity(mrs_string name) : MarSystem("TimeFreqPeakConnectivity", name) 00145 { 00146 addControls(); 00147 00148 tracebackIdx_ = 0; 00149 peakIndices_ = 0; 00150 path_ = 0; 00151 numBands_ = 0; 00152 multipleIndices = 0; 00153 } 00154 00155 TimeFreqPeakConnectivity::TimeFreqPeakConnectivity(const TimeFreqPeakConnectivity& a) : MarSystem(a) 00156 { 00157 // IMPORTANT! 00160 // Otherwise this would result in trying to deallocate them twice! 00161 ctrl_reso_ = getctrl("mrs_real/freqResolution"); 00162 00163 tracebackIdx_ = 0; 00164 peakIndices_ = 0; 00165 path_ = 0; 00166 numBands_ = 0; 00167 multipleIndices = 0; 00168 } 00169 00170 00171 TimeFreqPeakConnectivity::~TimeFreqPeakConnectivity() 00172 { 00173 FreeMemory (); 00174 if (multipleIndices) 00175 delete multipleIndices; 00176 } 00177 00178 MarSystem* 00179 TimeFreqPeakConnectivity::clone() const 00180 { 00181 // Every MarSystem should do this. 00182 return new TimeFreqPeakConnectivity(*this); 00183 } 00184 00185 void 00186 TimeFreqPeakConnectivity::addControls() 00187 { 00188 addctrl("mrs_string/frequencyIntervalInHz", "MARSYAS_EMPTY"); 00189 setctrlState("mrs_string/frequencyIntervalInHz", true); 00190 00191 addctrl("mrs_bool/inBark", false); 00192 addctrl("mrs_real/freqResolution", 25., ctrl_reso_); 00193 //addctrl("mrs_bool/inBark", true); 00194 //addctrl("mrs_real/freqResolution", .15, ctrl_reso_); 00195 00196 addctrl("mrs_natural/textureWindowSize", 0); 00197 } 00198 00199 void 00200 TimeFreqPeakConnectivity::AllocMemory (mrs_natural numSamples) 00201 { 00202 tracebackIdx_ = new unsigned char* [numBands_]; 00203 peakIndices_ = new mrs_natural* [numBands_]; 00204 for (mrs_natural i = 0; i < numBands_; i++) 00205 { 00206 tracebackIdx_[i] = new unsigned char [numSamples]; 00207 peakIndices_[i] = new mrs_natural[numSamples]; 00208 } 00209 path_ = new mrs_natural [numSamples]; 00210 00211 if (multipleIndices == 0) 00212 multipleIndices = new DoubleListEntries (); 00213 else 00214 multipleIndices->Reset (); 00215 } 00216 00217 void 00218 TimeFreqPeakConnectivity::FreeMemory () 00219 { 00220 if (tracebackIdx_) 00221 { 00222 for (mrs_natural i = 0; i < numBands_; i++) 00223 delete [] tracebackIdx_[i]; 00224 delete [] tracebackIdx_; 00225 } 00226 if (peakIndices_) 00227 { 00228 for (mrs_natural i = 0; i < numBands_; i++) 00229 delete [] peakIndices_[i]; 00230 delete [] peakIndices_; 00231 } 00232 00233 delete [] path_; 00234 } 00235 00236 void 00237 TimeFreqPeakConnectivity::myUpdate(MarControlPtr sender) 00238 { 00239 MRSDIAG("TimeFreqPeakConnectivity.cpp - TimeFreqPeakConnectivity:myUpdate"); 00240 00242 MarSystem::myUpdate(sender); 00243 00244 const mrs_bool isInBark = getctrl("mrs_bool/inBark")->to<mrs_bool>(); 00245 00246 FreeMemory (); 00247 00248 // compute matrix dimensions 00249 if(getctrl("mrs_string/frequencyIntervalInHz")->to<mrs_string>() != "MARSYAS_EMPTY") 00250 { 00251 realvec conv(2); 00252 string2parameters(getctrl("mrs_string/frequencyIntervalInHz")->to<mrs_string>(), conv, '_'); //[!] 00253 downFreq_ = (isInBark)? hertz2bark (conv(0)) : conv(0); 00254 upFreq_ = (isInBark)? hertz2bark (conv(1)) : conv(1); 00255 numBands_ = (mrs_natural)((upFreq_-downFreq_)/ctrl_reso_->to<mrs_real>() + .5); 00256 } 00257 else 00258 { 00259 numBands_ = 0; 00260 downFreq_ = 0; 00261 upFreq_ = 0; 00262 } 00263 textWinSize_ = getControl ("mrs_natural/textureWindowSize")->to<mrs_natural>(); 00264 00265 // create peak matrix 00266 peakMatrix_.create (numBands_, textWinSize_); 00267 costMatrix_.create (numBands_, textWinSize_); 00268 00269 updControl ("mrs_natural/onObservations", inSamples_); 00270 updControl ("mrs_natural/onSamples", inSamples_); 00271 00272 AllocMemory (textWinSize_); 00273 } 00274 00275 void 00276 TimeFreqPeakConnectivity::myProcess(realvec& in, realvec& out) 00277 { 00278 // get pitch resolution 00279 const mrs_real reso = ctrl_reso_->to<mrs_real>(); 00280 const mrs_bool isInBark = getctrl("mrs_bool/inBark")->to<mrs_bool>(); 00281 00282 // a matrix indicating where peaks are in the time frequency plane 00283 MRSASSERT(textWinSize_ >= in(1,inSamples_-1)-in(1,0)); 00284 peakMatrix_.stretch(numBands_, textWinSize_); 00285 00286 // init 00287 peakMatrix_.setval (1); 00288 for (mrs_natural t = 0; t < textWinSize_; t++) 00289 for (mrs_natural o=0; o < numBands_; o++) 00290 peakIndices_[o][t] = -1; 00291 00292 // initialized pseudo spectrogram representation 00293 for (mrs_natural t = 0; t < inSamples_; t++) 00294 { 00295 mrs_natural row = (isInBark)? Freq2RowIdx (in(0,t),reso) : Freq2RowIdx (bark2hertz(in(0,t)),reso), // input is in bark frequency, so we might have to convert it back 00296 col = (mrs_natural)(in(1,t)-in(1,0)+.1); 00297 MRSASSERT(col >= 0 && col < textWinSize_); 00298 MRSASSERT(row >= 0 && row < numBands_); 00299 // check whether more than one peak are at that matrix pos, i.e. we already have set the entry 00300 if (peakIndices_[row][col] != -1) 00301 { 00302 multipleIndices->AddIndex (row,col,peakIndices_[row][col]); 00303 multipleIndices->AddIndex (row,col,t); 00304 peakIndices_[row][col] = -2; 00305 } 00306 else 00307 peakIndices_[row][col] = t; // a matrix indicating which matrix bin corresponds to which input index 00308 peakMatrix_(row, col) = 0; 00309 } 00310 00311 // initialize output 00312 out.setval (costInit); 00313 00314 #ifdef MATLAB_DBG_OUT 00315 #ifdef MARSYAS_MATLAB 00316 MATLAB_PUT(peakMatrix_, "peakMatrix"); 00317 MATLAB_EVAL ("figure(1),imagesc(peakMatrix),colorbar"); 00318 #endif 00319 #endif 00320 00321 // iteration over all pairs 00322 for (mrs_natural t = 0; t < inSamples_; t++) 00323 { 00324 for (mrs_natural o=inSamples_-1; o >= t; o--) 00325 { 00326 // don't compute distance if we already have it 00327 if (out(t,o) != costInit) 00328 continue; 00329 00330 // get peak matrix indices 00331 mrs_natural rowt = (isInBark)? Freq2RowIdx (in(0,t),reso) : Freq2RowIdx (bark2hertz(in(0,t)),reso), // input is in bark frequency, so we might have to convert it back 00332 rowo = (isInBark)? Freq2RowIdx (in(0,o),reso) : Freq2RowIdx (bark2hertz(in(0,o)),reso), // input is in bark frequency, so we might have to convert it back 00333 colt = (mrs_natural)(in(1,t)-in(1,0)+.1), 00334 colo = (mrs_natural)(in(1,o)-in(1,0)+.1), 00335 pathLength; 00336 00337 MRSASSERT(colt >= 0 && colt < textWinSize_); 00338 MRSASSERT(colo >= 0 && colo < textWinSize_); 00339 MRSASSERT(rowt >= 0 && rowt < numBands_); 00340 MRSASSERT(rowo >= 0 && rowo < numBands_); 00341 00342 00343 // self similarity and similarity with overlapping peaks 00344 if ((t == o) || (rowt == rowo && colt == colo)) 00345 { 00346 SetOutput(out, 0, rowt, colt, rowo, colo); 00347 continue; 00348 } 00349 00350 // check if path calculation makes sense with the current dp step size 00351 if (abs(rowt - rowo) > abs(colt-colo)) 00352 { 00353 SetOutput(out, 1, rowt, colt, rowo, colo); 00354 continue; 00355 } 00356 00357 // let's calculate only from left to right 00358 if (colo < colt) 00359 continue; 00360 00361 // dynamic programming 00362 CalcDp (peakMatrix_, rowt, colt, rowo, colo); 00363 pathLength = colo-colt+1; 00364 00365 #ifdef MATLAB_DBG_OUT 00366 #ifdef MARSYAS_MATLAB 00367 MATLAB_PUT(costMatrix_, "cost"); 00368 MATLAB_EVAL ("figure(2),imagesc(cost,[0 10]),colorbar"); 00369 #endif 00370 #endif 00371 00372 // set cost for this path and all subpaths 00373 for (mrs_natural i = 0; i < pathLength; i++) 00374 { 00375 if (peakMatrix_(path_[i],colt + i) > 0) 00376 { 00377 MRSASSERT(i>0); 00378 continue; 00379 } 00380 for (mrs_natural j = i; j < pathLength; j++) 00381 { 00382 if (peakMatrix_(path_[j],colt + j) > 0) 00383 continue; // this path entry is not a peak 00384 00385 mrs_real cost = (costMatrix_(path_[j],colt + j) - costMatrix_(path_[i],colt + i)) / (j-i+1); 00386 MRSASSERT (cost >= 0 && cost <=1); 00387 00388 // assign this (and related) output 00389 SetOutput(out, cost, path_[i], colt + i, path_[j], colt + j); 00390 } 00391 } 00392 } 00393 } 00394 multipleIndices->Reset (); 00395 00396 // scale output because of RBF without std?? 00397 //out *= 10.; 00398 00399 // convert distance to similarity 00400 //out *= -1; 00401 //out += 1; 00402 #ifdef SANITY_CHECK 00403 for (mrs_natural t=0; t < inSamples_; t++) 00404 for (mrs_natural o=0; o < inSamples_; o++) 00405 MRSASSERT (out(t,o) >= 0 && out(t,o) <=1); 00406 #endif 00407 #ifdef MATLAB_DBG_OUT 00408 #ifdef MARSYAS_MATLAB 00409 MATLAB_PUT(out, "out"); 00410 MATLAB_EVAL ("figure(3),imagesc(out),colorbar"); 00411 #endif 00412 #endif 00413 00414 } 00415 00416 00417 void TimeFreqPeakConnectivity::SetOutput(mrs_realvec &out, const mrs_real cost, mrs_natural idxRowA, mrs_natural idxColA, mrs_natural idxRowB, mrs_natural idxColB) 00418 { 00419 mrs_natural ml = 0, 00420 nl = 0, 00421 rowIdx = peakIndices_[idxRowA][idxColA], 00422 colIdx = peakIndices_[idxRowB][idxColB]; 00423 00424 if (rowIdx == -2) 00425 { 00426 ml = multipleIndices->GetNumIndices (idxRowA, idxColA); 00427 MRSASSERT(ml>0); 00428 rowIdx = multipleIndices->GetIndex (idxRowA, idxColA, 0); 00429 } 00430 if (colIdx == -2) 00431 { 00432 nl = multipleIndices->GetNumIndices (idxRowB, idxColB); 00433 MRSASSERT(nl>0); 00434 colIdx = multipleIndices->GetIndex (idxRowB, idxColB, 0); 00435 } 00436 00437 if (out(rowIdx, colIdx) != costInit) 00438 { 00439 MRSASSERT(out(rowIdx, colIdx) == cost); 00440 return; // cost has already been computed - nothing to do 00441 } 00442 //if ((rowIdx==53 && colIdx==55) || (rowIdx==55 && colIdx==55) || (rowIdx==55 && colIdx==53)) 00443 // cout << out(t,o) << endl; 00444 00445 out(rowIdx,colIdx) = cost; 00446 out(colIdx,rowIdx) = cost; 00447 00448 // only for the case of overlapping peaks 00449 if (ml > 0 || nl > 0) 00450 { 00451 mrs_natural m,n; 00452 00453 // both > 0 00454 if (ml * nl) 00455 { 00456 for (m = 0; m < ml; m++) 00457 { 00458 rowIdx = multipleIndices->GetIndex (idxRowA, idxColA, m); 00459 for (n = 0; n < nl; n++) 00460 { 00461 colIdx = multipleIndices->GetIndex (idxRowB, idxColB, n); 00462 //if ((rowIdx==53 && colIdx==55) || (rowIdx==55 && colIdx==55) || (rowIdx==55 && colIdx==53)) 00463 // cout << out(t,o) << endl; 00464 out(rowIdx,colIdx) = cost; 00465 out(colIdx,rowIdx) = cost; 00466 } 00467 } 00468 } 00469 else if (ml > 0) 00470 { 00471 for (m = 0; m < ml; m++) 00472 { 00473 rowIdx = multipleIndices->GetIndex (idxRowA, idxColA, m); 00474 //if ((rowIdx==53 && colIdx==55) || (rowIdx==55 && colIdx==55) || (rowIdx==55 && colIdx==53)) 00475 // cout << out(t,o) << endl; 00476 out(rowIdx,colIdx) = cost; 00477 out(colIdx,rowIdx) = cost; 00478 } 00479 00480 } 00481 else if (nl > 0) 00482 { 00483 for (n = 0; n < nl; n++) 00484 { 00485 colIdx = multipleIndices->GetIndex (idxRowB, idxColB, n); 00486 //if ((rowIdx==53 && colIdx==55) || (rowIdx==55 && colIdx==55) || (rowIdx==55 && colIdx==53)) 00487 // cout << out(t,o) << endl; 00488 out(rowIdx,colIdx) = cost; 00489 out(colIdx,rowIdx) = cost; 00490 } 00491 } 00492 } 00493 } 00494 00495 mrs_natural TimeFreqPeakConnectivity::Freq2RowIdx (mrs_real barkFreq, mrs_real bres) 00496 { 00497 mrs_natural result = (mrs_natural)((barkFreq-downFreq_)/bres +.5); 00498 00499 if (result < 0) 00500 result = 0; 00501 if (result >= numBands_) 00502 result =numBands_-1; 00503 00504 return result; 00505 } 00506 00507 void TimeFreqPeakConnectivity::InitMatrix (mrs_realvec &Matrix, unsigned char **traceback, mrs_natural rowIdx, mrs_natural colIdx) 00508 { 00509 mrs_natural i,j, 00510 numRows = Matrix.getRows (), 00511 numCols = Matrix.getCols (); 00512 mrs_natural iCol; 00513 00514 Matrix.setval(0); 00515 00516 traceback[rowIdx][colIdx] = kFromLeft; 00517 00518 // left of point of interest 00519 for (i = 0; i < numRows; i++) 00520 { 00521 for (j = 0; j < colIdx; j++) 00522 { 00523 Matrix(i,j) = costInit; 00524 traceback[i][j] = kFromLeft; 00525 } 00526 } 00527 //cout << Matrix << endl; 00528 //upper left corner 00529 for (i = 0; i < rowIdx; i++) 00530 { 00531 iCol = MyMin (rowIdx - i + colIdx, numCols); 00532 for (j = colIdx; j < iCol; j++) 00533 { 00534 Matrix(i,j) = costInit; 00535 traceback[i][j] = kFromLeft; 00536 } 00537 } 00538 //cout << Matrix << endl; 00539 // lower left corner 00540 for (i = rowIdx + 1; i < numRows; i++) 00541 { 00542 iCol = MyMin (i - rowIdx + colIdx, numCols); 00543 for (j = colIdx; j < iCol; j++) 00544 { 00545 Matrix(i,j) = costInit; 00546 traceback[i][j] = kFromLeft; 00547 } 00548 } 00549 //cout << Matrix << endl; 00550 } 00551 00552 void TimeFreqPeakConnectivity::CalcDp (mrs_realvec &Matrix, mrs_natural startr, mrs_natural startc, mrs_natural stopr, mrs_natural stopc) 00553 { 00554 mrs_natural i,j, 00555 numRows = Matrix.getRows (), 00556 numCols = Matrix.getCols (); 00557 mrs_real prevCost[kNumDirections] = {0,0,0}; 00558 00559 costMatrix_.stretch ( numRows, numCols); 00560 00561 // initialize cost and traceback matrix 00562 // upper and lower left corner 00563 InitMatrix (costMatrix_, tracebackIdx_, startr, startc); 00564 costMatrix_(startr, startc) = Matrix(startr, startc); 00565 00566 // compute cost matrix 00567 for (j = startc+1; j <= stopc; j++) 00568 { 00569 mrs_natural rowLoopStart = MyMax (0, startr - (j - startc)), 00570 rowLoopStop = MyMin (numRows, startr + (j-startc) + 1); 00571 for (i = rowLoopStart; i < rowLoopStop; i++) 00572 { 00573 prevCost[kFromLeft] = costMatrix_(i,j-1); 00574 prevCost[kFromUp] = (i >= numRows-1) ? costInit : costMatrix_(i+1, j-1); // the if isn't very nice here.... 00575 prevCost[kFromDown] = (i <= 0) ? costInit : costMatrix_(i-1, j-1); 00576 // assign cost 00577 costMatrix_(i,j) = Matrix(i,j) + dtwFindMin (prevCost, tracebackIdx_[i][j]); 00578 } 00579 } 00580 00581 // compute path 00582 i = stopr; 00583 for (j = stopc; j >= startc; j--) 00584 { 00585 path_[j-startc] = i; 00586 i -= (kFromLeft - tracebackIdx_[i][j]); // note: does work only for kFromUp, kFromLeft, kFromDown 00587 } 00588 }