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