Marsyas  0.5.0-beta1
/Users/jleben/code/marsyas/src/marsyas/marsystems/F0Analysis.cpp
Go to the documentation of this file.
00001 #include <map>
00002 #include <algorithm>
00003 #include "F0Analysis.h"
00004 
00005 #define ROUND(x)    (mrs_natural)floor(x+0.5)
00006 
00007 
00008 using std::ostringstream;
00009 using std::vector;
00010 using std::abs;
00011 
00012 using namespace Marsyas;
00013 
00014 F0Analysis::F0Analysis(mrs_string inName)
00015   :MarSystem("F0Analysis",inName)
00016 {
00017   addControls();
00018 }
00019 
00020 F0Analysis::F0Analysis(const F0Analysis& inToCopy)
00021   :MarSystem(inToCopy)
00022 {
00023   ctrl_SampleRate_ = getctrl("mrs_real/SampleRate");
00024   ctrl_NrOfHarmonics_ = getctrl("mrs_natural/NrOfHarmonics");
00025   ctrl_F0Weight_ = getctrl("mrs_real/F0Weight");
00026   ctrl_Attenuation_ = getctrl("mrs_real/Attenuation");
00027   ctrl_Tolerance_ = getctrl("mrs_real/Tolerance");
00028   ctrl_LowestF0_ = getctrl("mrs_real/LowestF0");
00029   ctrl_Compression_ = getctrl("mrs_real/Compression");
00030 
00031   SampleRate_ = inToCopy.SampleRate_;
00032   NrOfHarmonics_ = inToCopy.NrOfHarmonics_;
00033   F0Weight_ = inToCopy.F0Weight_;
00034   Attenuation_ = inToCopy.Attenuation_;
00035   Tolerance_ = inToCopy.Tolerance_;
00036   LowestF0_ = inToCopy.LowestF0_;
00037   Compression_ = inToCopy.Compression_;
00038 }
00039 
00040 F0Analysis::~F0Analysis() {}
00041 
00042 MarSystem* F0Analysis::clone() const
00043 {
00044   return new F0Analysis(*this);
00045 }
00046 
00047 void F0Analysis::addControls()
00048 {
00049   addctrl("mrs_real/SampleRate",8000.f,ctrl_SampleRate_);
00050   addctrl("mrs_natural/NrOfHarmonics",5,ctrl_NrOfHarmonics_);
00051   addctrl("mrs_real/F0Weight",0.5,ctrl_F0Weight_);
00052   addctrl("mrs_real/Attenuation",0.75,ctrl_Attenuation_);
00053   addctrl("mrs_real/Tolerance",0.03,ctrl_Tolerance_);
00054   addctrl("mrs_real/LowestF0",100.,ctrl_LowestF0_);
00055   addctrl("mrs_real/Compression",0.5,ctrl_Compression_);
00056   addctrl("mrs_real/ChordEvidence",0.);
00057 
00058   ctrl_SampleRate_->setState(true);
00059   ctrl_NrOfHarmonics_->setState(true);
00060   ctrl_F0Weight_->setState(true);
00061   ctrl_Attenuation_->setState(true);
00062   ctrl_Tolerance_->setState(true);
00063   ctrl_LowestF0_->setState(true);
00064   ctrl_Compression_->setState(true);
00065 
00066   SampleRate_ = 8000.f;
00067   NrOfHarmonics_ = 5;
00068   F0Weight_ = 0.5;
00069   Attenuation_ = 0.75;
00070   Tolerance_ = 0.03;
00071   LowestF0_ = 100.;
00072   Compression_ = 0.5;
00073 }
00074 
00075 void F0Analysis::myUpdate(MarControlPtr inSender)
00076 {
00077   SampleRate_ = ctrl_SampleRate_->to<mrs_real>();
00078   NrOfHarmonics_ = ctrl_NrOfHarmonics_->to<mrs_natural>();
00079   F0Weight_ = ctrl_F0Weight_->to<mrs_real>();
00080   Attenuation_ = ctrl_Attenuation_->to<mrs_real>();
00081   Tolerance_ = ctrl_Tolerance_->to<mrs_real>();
00082   LowestF0_ = ctrl_LowestF0_->to<mrs_real>();
00083   Compression_ = ctrl_Compression_->to<mrs_real>();
00084 
00085   MarSystem::myUpdate(inSender);
00086 }
00087 
00088 void F0Analysis::myProcess(realvec& inVec, realvec& outVec)
00089 {
00090   F2Fs theF0ToFks;      // Map between F0 and higher harmonics
00091   HarmMap theHarmSums;
00092   FindCandidateF0s(inVec, theHarmSums, theF0ToFks);
00093   SelectUnrelatedF0s(inVec, theHarmSums, theF0ToFks, outVec);
00094   updControl("mrs_real/ChordEvidence",ChordEvidence_);
00095 }
00096 
00097 bool F0Analysis::FindCandidateF0s(const realvec& inPeaks,
00098                                   HarmMap& outHarmSums, F2Fs& outF0ToFks) const
00099 {
00100   /* For each F0 > FLower, search for harmonically related freqs Fks
00101      1. Compute harmonic sum (S) and store (S,F0) in the map outHarmSums
00102      2. Store all Fks assigned to F0 in m_F0ToFks */
00103   outHarmSums.clear();
00104   outF0ToFks.clear();
00105 
00106   mrs_real theFreqStep = SampleRate_/(2.*inPeaks.getSize());
00107   for (mrs_natural i=0; i<inPeaks.getSize(); ++i) {
00108     mrs_real theF0 = (mrs_real)i*theFreqStep;
00109 
00110     // F0 > FLower
00111     if (inPeaks(i)>0 && theF0 >= LowestF0_) {
00112 
00113       // Compute harmonic sum & energy
00114       vector<double> theAssignedFks;
00115       mrs_real theSum = 0.0f, theNormFactor = 0.0f;
00116       for (mrs_natural j=i+1; j<inPeaks.getSize(); j++) {
00117         if (inPeaks(j)>0) {
00118           mrs_real theFk = (mrs_real)j*theFreqStep;
00119 
00120           /* Check whether Fk is one of the considered harmonics of F0:
00121              1. Compute k, the closest integer to Fk/F0
00122              2. Check whether |Fk/k-F0| <= tolerance x F0 */
00123           int k = ROUND((mrs_real)j/(mrs_real)i);
00124           if (k > 1 && k <= NrOfHarmonics_ &&
00125               abs(theFk/(mrs_real)k-theF0) <= Tolerance_*theF0) {
00126 
00127             theAssignedFks.push_back(theFk);
00128             double tmp = pow(Attenuation_,(double)k);
00129             theSum += pow(inPeaks(j),Compression_)*tmp;
00130             theNormFactor += tmp;
00131           }
00132         }
00133       }
00134 
00135       // Add F0
00136       if (theAssignedFks.size()>0) {
00137         outHarmSums[pow(inPeaks(i),Compression_*F0Weight_)
00138                     *pow(theSum/theNormFactor,1.-F0Weight_)] = theF0;
00139         outF0ToFks[theF0] = theAssignedFks;
00140       }
00141     }
00142   }
00143   return true;
00144 }
00145 
00146 bool F0Analysis::SelectUnrelatedF0s(const realvec& inPeaks,
00147                                     const HarmMap inHarmSums, const F2Fs& inF0ToFks, realvec& outNoteEvidence)
00148 {
00149   outNoteEvidence.setval(0);
00150 
00151   FreqMap thePeaks;
00152   mrs_real theFreqStep = SampleRate_/(2.*inPeaks.getSize());
00153   for (mrs_natural i=0; i<inPeaks.getSize(); ++i)
00154     if (inPeaks(i)>0)
00155       thePeaks[(mrs_real)i*theFreqStep] = inPeaks(i);
00156 
00157   ChordEvidence_ = 0.0f;
00158   mrs_natural theNrOfPitches = 0;
00159   if (!inHarmSums.empty()) {
00160     /* Select candidates F0 (called Fc) that are not harmonically related.
00161        At first, select Fc with the highest harmonic sum (S) value.
00162        At second, select the other Fcs in decreasing order of S while
00163        Fc is not related to a selected Fc (HARMONIC or SUBHARMONIC)*/
00164 
00165     HarmMap::const_iterator Cand;   // Candidate F0
00166     FreqMap::iterator Sel;          // Selected F0
00167 
00168     // Add first candidate
00169     Cand = inHarmSums.begin();
00170     double F0c = Cand->second;
00171     outNoteEvidence(ROUND(F0c/theFreqStep)) =
00172       ComputePowerOfF0(thePeaks, inF0ToFks, F0c);
00173     theNrOfPitches++;
00174     Cand++;
00175 
00176     // Proceed with other candidates
00177     mrs_real theHypPower = ComputePowerOfF0(thePeaks, inF0ToFks, F0c);
00178     mrs_real theAllPower = ComputePowerOfInput(thePeaks);
00179     while(Cand!=inHarmSums.end()) {
00180 
00181       // Check relationship with selected F0
00182       F0c = Cand->second;
00183       bool theRelFlag = false;
00184 
00185       for (mrs_natural i=0; i<outNoteEvidence.getSize(); ++i)
00186       {
00187         if (outNoteEvidence(i) > 0)
00188         {
00189           double F0 = (mrs_real)i*theFreqStep;
00190           int k = ROUND(F0c/F0);
00191           int l = ROUND(F0/F0c);
00192 
00193           // theRelFlag = true if Sel and Cand are harmonically related
00194           theRelFlag = (k > 0 && k <= NrOfHarmonics_ &&
00195                         abs(F0c/(double)k-F0) <= Tolerance_*F0) ||
00196                        (l > 0 && l <= NrOfHarmonics_ &&
00197                         abs((double)l*F0c-F0) <= Tolerance_*F0);
00198 
00199           if (theRelFlag)
00200             break;
00201         }
00202       }
00203 
00204       if (!theRelFlag)
00205       {
00206         outNoteEvidence(ROUND(F0c/theFreqStep)) = ComputePowerOfF0(thePeaks, inF0ToFks, F0c);
00207         theHypPower = ComputePowerOfHyp(thePeaks, inF0ToFks, outNoteEvidence);
00208         theNrOfPitches++;
00209       }
00210       Cand++;
00211     }
00212 
00213     // Normalize note relevances
00214     mrs_real theFactor = 0.0f;
00215     for (mrs_natural i=0; i<outNoteEvidence.getSize(); ++i)
00216       theFactor += outNoteEvidence(i);
00217     for (mrs_natural i=0; i<outNoteEvidence.getSize(); ++i)
00218       outNoteEvidence(i) /= theFactor;
00219 
00220     // Compute chord evidence if nr of notes >= 2
00221     if (theNrOfPitches>=2)
00222       ChordEvidence_= theHypPower/theAllPower;
00223   }
00224   return true;
00225 }
00226 
00227 mrs_real F0Analysis::ComputePowerOfF0(const FreqMap inPeaks, const F2Fs& inF0ToFks, double inF0) const {
00228   FreqMap::const_iterator iter1 = inPeaks.find(inF0);
00229   mrs_real thePower = pow(iter1->second,Compression_);
00230 
00231   F2Fs::const_iterator iter2 = inF0ToFks.find(inF0);
00232   vector<double> theFks = iter2->second;
00233   for (unsigned long i=0; i<theFks.size(); ++i) {
00234     iter1 = inPeaks.find(theFks[i]);
00235     thePower += pow(iter1->second,Compression_);
00236   }
00237   return thePower;
00238 }
00239 
00240 mrs_real F0Analysis::ComputePowerOfInput(const FreqMap inPeaks) const {
00241   mrs_real thePower = 0.0f;
00242   FreqMap::const_iterator iter;
00243   for (iter=inPeaks.begin(); iter!=inPeaks.end(); iter++)
00244     thePower += iter->second * iter->second;
00245   return thePower;
00246 }
00247 
00248 mrs_real F0Analysis::ComputePowerOfHyp(const FreqMap inPeaks,
00249                                        const F2Fs& inF0ToFks, realvec& inNoteEvidence) const
00250 {
00251   /* For each selected candidate F0, search the assigned higher
00252      harmonics and store them in the vector Tmp*/
00253   vector<double> Tmp;
00254   mrs_real theFreqStep = SampleRate_/(2*inNoteEvidence.getSize());
00255   for (mrs_natural i=0; i<inNoteEvidence.getSize(); ++i)
00256   {
00257     if (inNoteEvidence(i) > 0)
00258     {
00259       F2Fs::const_iterator iter2 =
00260         inF0ToFks.find((mrs_real)i*theFreqStep);
00261       vector<double> theHarm = iter2->second;
00262       for (unsigned long i=0; i<theHarm.size(); ++i)
00263         Tmp.push_back(theHarm[i]);
00264     }
00265   }
00266 
00267   // Remove duplicate frequencies
00268   sort(Tmp.begin(), Tmp.end());
00269   unique(Tmp.begin(), Tmp.end());
00270 
00271   // Compute power of unique frequencies
00272   mrs_real thePower = 0.0f;
00273   FreqMap::const_iterator iter;
00274   for (size_t i=0; i<Tmp.size(); ++i) {
00275     iter = inPeaks.find(Tmp[i]);
00276     thePower += iter->second * iter->second;
00277   }
00278   return thePower;
00279 }