Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/TimeFreqPeakConnectivity.cpp
Go to the documentation of this file.
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 }