libs/Music/BubbleAlgo.cpp

   1 // Copyright 2004 "Gilles Degottex"
   2
   3 // This file is part of "Music"
   4
   5 // "Music" is free software; you can redistribute it and/or modify
   6 // it under the terms of the GNU Lesser General Public License as published by
   7 // the Free Software Foundation; either version 2.1 of the License, or
   8 // (at your option) any later version.
   9 //
  10 // "Music" is distributed in the hope that it will be useful,
  11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  13 // GNU Lesser General Public License for more details.
  14 //
  15 // You should have received a copy of the GNU Lesser General Public License
  16 // along with this program; if not, write to the Free Software
  17 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  18
  19
  20 #include "BubbleAlgo.h"
  21
  22 #include <cassert>
  23 #include <cmath>
  24 #include <deque>
  25 #include <iostream>
  26 #include <limits>
  27 using namespace std;
  28 #include <CppAddons/Math.h>
  29 using namespace Math;
  30
  31 #include "Music.h"
  32
  33 #define MUSIC_DEBUG
  34 #ifdef MUSIC_DEBUG
  35         #define LOG(a)  a
  36 #else
  37         #define LOG(a)
  38 #endif
  39
  40
  41 namespace Music
  42 {
  43         int gcd(int a, int b)
  44         {
  45                 if (b == 0) return a;
  46
  47                 return gcd(b, a%b);
  48         }
  49         int gcd(vector<int> v)
  50         {
  51                 if(v.empty())   return -1;
  52
  53                 int s = v[0];
  54                 for(size_t i=1; i<v.size(); i++)
  55                         s = gcd(s, v[i]);
  56
  57                 return s;
  58         }
  59
  60         void BubbleAlgo::init()
  61         {
  62                 m_min_length = int(GetSamplingRate()/h2f(GetSemitoneMax()));
  63                 m_max_length = int(GetSamplingRate()/h2f(GetSemitoneMin()));
  64                 setMinMaxLength(m_min_length, m_max_length);
  65                 m_bubbles.resize(m_max_length);
  66                 m_waves.resize(m_max_length);
  67                 m_waves_best.resize(m_max_length);
  68
  69                 cerr << "BubbleAlgo::init [" << m_min_length << ";" << m_max_length << "]" << endl;
  70
  71                 m_error_threshold = 0.1;
  72                 m_conv_threshold = 0.1;
  73                 double gauss_factor = 2.0;
  74                 size_t latency_factor = 2;
  75
  76                 double u = Usefull(Win_Sinc(gauss_factor));
  77
  78                 for(size_t s=m_min_length; s<m_max_length; s++)
  79                 {
  80                         m_bubbles[s].s = s;
  81                         m_bubbles[s].length = s;
  82                         while(m_bubbles[s].length<100)
  83                                 m_bubbles[s].length += s;
  84
  85                         m_waves[s].resize(latency_factor*s);
  86                         m_waves_best[s].resize(m_waves[s].size());
  87 //                      double c = - 2.0*Math::Pi * m_freq / sampling_rate;
  88                         double c = - 2.0*Math::Pi / s;
  89                         double d = (2.0/m_waves[s].size());
  90
  91                         for(size_t j=0; j<m_waves[s].size(); j++)
  92                                 m_waves[s][j] = exp(complex<double>(0.0, c*j)) * d * win_sinc(j/double(m_waves[s].size()), gauss_factor)/u;
  93
  94                         complex<double> b(0.0,0.0);
  95                         for(int i=s-1; i>=0; i--)
  96                         {
  97                                 for(size_t l=0; l<latency_factor; l++)
  98                                         b += m_waves[s][i+l*s]*sin(2.0*Math::Pi*(i+l*s)/s);
  99                                 m_waves_best[s][i] = normm(b);
 100 //                              if(s==109)
 101 //                                      cerr<<"nv="<<m_waves_best[s][i]<<endl;
 102                         }
 103                 }
 104
 105                 cerr << "BubbleAlgo::init " << m_waves.back().size() << endl;
 106         }
 107
 108         BubbleAlgo::BubbleAlgo()
 109         : Algorithm(0.1)
 110         {
 111                 init();
 112         }
 113
 114         //      double err = abs(buff[s] - buff[0])/((buff[s]+buff[0])/2);
 115         //      double err = abs(buff[s] - buff[0])/(max(abs(buff[s]),abs(buff[0])));
 116
 117         /*
 118          * * normalize errors ?
 119          * for too high results
 120          *   * longer correlation length
 121          *     6 is not statisticaly a valuable average size ...
 122          *   * random shifting
 123          *     + kill high errors
 124          *     - more latency (*1.5 or more)
 125          * for too low results
 126          *   * divide by num of present multiples in the sort
 127          *
 128          * max is not stable enough
 129          * difficult to use conv because there is sound with fondamental with zero energy
 130          */
 131         void BubbleAlgo::apply(const deque<double>& buff)
 132         {
 133                 m_wave_length = 0;
 134
 135                 if(buff.size()<m_waves.back().size())   return;
 136
 137                 LOG(cerr<<"BubbleAlgo::apply min_length="<<m_min_length<<" max_length="<<m_max_length<<endl;)
 138
 139                 Type max_vol = 0.0;
 140                 for(size_t i=0; i<m_max_length; i++)
 141                         max_vol = max(Type(max_vol), Type(buff[i]));
 142
 143                 m_error_threshold = 0.33;
 144                 m_conv_threshold = 0.0;
 145                 // use relative thresholds
 146                 double err_threshold = m_error_threshold*max_vol;
 147                 size_t err_drop = 0, err_drop_t = 0;
 148                 double conv_threshold = m_conv_threshold*max_vol;
 149                 size_t conv_drop = 0;
 150                 size_t score_drop = 0, score_drop_t = 0;
 151                 size_t sgn_drop = 0, sgn_drop_t = 0;
 152                 size_t mult_drop = 0, mult_drop_t = 0;
 153
 154                 // init
 155                 m_sort.clear();
 156                 m_best_set.clear();
 157                 size_t map_size = 0;
 158                 size_t max_n = 0;
 159                 for(size_t s=m_min_length; s<m_max_length; s++)
 160                 {
 161                         m_bubbles[s].score = 0.0;
 162                         m_bubbles[s].err = 0.0;
 163                         m_bubbles[s].conv = complex<double>(0.0,0.0);
 164                         m_bubbles[s].count = 0;
 165                         m_bubbles[s].todrop = false;
 166                         m_bubbles[s].gcd_count = 0;
 167
 168                         if(sgn(buff[s])==sgn(buff[0]) || sgn(buff[s+1])==sgn(buff[1]))
 169                         {
 170                                 m_bubbles[s].err = diff(buff, s, 0);
 171 //                              m_bubbles[s].conv = m_waves[s][0] * buff[0];
 172 //                              m_bubbles[s].conv += m_waves[s][s] * buff[s];
 173                                 m_bubbles[s].count++;
 174 //                              m_bubbles[s].score = max(m_bubbles[s].err, 1.0-normm(m_bubbles[s].conv));
 175                                 m_bubbles[s].score = m_bubbles[s].err;
 176
 177 //                              if(m_bubbles[s].err<err_threshold && normm(m_bubbles[s].conv)>conv_threshold/s)
 178                                 if(m_bubbles[s].err<err_threshold)
 179                                 {
 180                                         m_sort[m_bubbles[s].score].push_back(s);
 181                                         map_size++;
 182                                         max_n += s;
 183                                 }
 184                         }
 185                 }
 186                 size_t init_map_size = map_size;
 187                 LOG(cerr<<"map size="<<map_size<<" ("<<int(log(float(map_size))+1)<<")  droped="<<(m_max_length-m_min_length)-map_size<<endl;)
 188
 189                 if(map_size==0) return;
 190
 191                 // seek
 192                 bool search = true;
 193                 size_t n=0;
 194                 size_t f=0;
 195                 size_t current_best = 0;
 196                 while(!m_sort.empty() && search)
 197                 {
 198                         n++;
 199
 200                         // extract the next bubble to raise, lower or drop in the chart
 201                         size_t s = (*(m_sort.begin())).second.front();
 202                         (*(m_sort.begin())).second.pop_front();
 203                         map_size--;
 204                         if((*(m_sort.begin())).second.empty())
 205                                 m_sort.erase(m_sort.begin());
 206
 207                         // 2. should be dropped by multiplicity
 208                         if(m_bubbles[s].todrop)
 209                         {
 210                                 mult_drop++;
 211                                 mult_drop_t += n;
 212                         }
 213                         // 1. if sign is the same (or nearly)
 214                         else if(!(sgn(buff[s+m_bubbles[s].count])==sgn(buff[m_bubbles[s].count]) || sgn(buff[s+m_bubbles[s].count+1])==sgn(buff[m_bubbles[s].count+1])))
 215                         {
 216                                 sgn_drop++;
 217                                 sgn_drop_t += n;
 218                         }
 219                         else
 220                         {
 221                                 m_bubbles[s].err += diff(buff, s+m_bubbles[s].count, m_bubbles[s].count);
 222 //                              if(m_bubbles[s].count<s)
 223 //                              {
 224 //                                      m_bubbles[s].conv += m_waves[s][m_bubbles[s].count] * buff[m_bubbles[s].count];
 225 //                                      m_bubbles[s].conv += m_waves[s][m_bubbles[s].count+s] * buff[m_bubbles[s].count+s];
 226 //                              }
 227                                 m_bubbles[s].count++;
 228 //                              m_bubbles[s].score = max((m_bubbles[s].err/m_bubbles[s].count), 1.0-normm(m_bubbles[s].conv));
 229                                 m_bubbles[s].score = m_bubbles[s].err/m_bubbles[s].count;
 230
 231                                 // if finished
 232                                 if(m_bubbles[s].count>=m_bubbles[s].length)
 233                                 {
 234                                         if(m_bubbles[s].err/m_bubbles[s].count>err_threshold)
 235                                                 err_drop++;
 236 //                                      else if(normm(m_bubbles[s].conv)<conv_threshold)
 237 //                                              conv_drop++;
 238                                         else
 239                                         {
 240                                                 m_best_set[m_bubbles[s].score].push_back(s);
 241                                                 m_bubbles[s].n = n;
 242                                                 m_bubbles[s].f = f;
 243                                                 f++;
 244                                                 if(f==1 || m_bubbles[s].score<m_bubbles[current_best].score)
 245                                                         current_best = s;
 246
 247                                                 // drop nearby and all mult bubble
 248                                                 size_t ms=2*s;
 249                                                 if(s+1<m_max_length)    m_bubbles[s+1].todrop = true;
 250                                                 if(s-1>=m_min_length)   m_bubbles[s-1].todrop = true;
 251                                                 while(ms<m_max_length)
 252                                                 {
 253                                                         m_bubbles[ms].todrop = true;
 254                                                         if(ms+1<m_max_length)   m_bubbles[ms+1].todrop = true;
 255                                                         if(ms-1>=m_min_length)  m_bubbles[ms-1].todrop = true;
 256                                                         ms += s;
 257                                                 }
 258                                                 size_t d=2;
 259                                                 m_bubbles[s].gcd_count++;
 260 //                                              int old_ds = s;
 261                                                 while(s/d>=m_min_length)
 262                                                 {
 263                                                         size_t ds1 = s/d;
 264                                                         size_t ds2 = 1+ds1;
 265 //                                                      size_t ds2 = ds1;
 266 //                                                      for(int i=old_ds-1; i>ds2; i--)
 267 //                                                              m_bubbles[ds1].todrop = true;
 268 //                                                      old_ds = ds1;
 269                                                         m_bubbles[ds1].gcd_count++;
 270                                                         if(m_bubbles[ds1].gcd_count<f)
 271                                                                 m_bubbles[ds1].todrop = true;
 272                                                         m_bubbles[ds2].gcd_count++;
 273                                                         if(m_bubbles[ds2].gcd_count<f)
 274                                                                 m_bubbles[ds2].todrop = true;
 275
 276                                                         d++;
 277                                                 }
 278                                         }
 279                                 }
 280                                 // if there is no chance for s to be lower than error threshold
 281 //                              else if(m_bubbles[s].err/s>err_threshold)
 282 //                              {
 283 //                                      err_drop++;
 284 //                                      err_drop_t += n;
 285 //                              }
 286                                 else
 287                                 {
 288                                         // all test past, can continue
 289                                         m_sort[m_bubbles[s].score].push_back(s);
 290                                         map_size++;
 291                                 }
 292                         }
 293
 294 //                      if(f>0) search = false;
 295
 296                         // all possible multiples have been finished
 297                         if(f>m_max_length/m_min_length) search = false;
 298                 }
 299
 300                 int best_tot_c = 0;
 301
 302                 size_t gcds = 0;
 303
 304                 for(map<Type,list<size_t> >::iterator it=m_best_set.begin(); it!=m_best_set.end(); ++it)
 305                 {
 306                         for(list<size_t>::iterator itl=(*it).second.begin(); itl!=(*it).second.end(); ++itl)
 307                         {
 308                                 size_t s = *itl;
 309 //                              if(m_bubbles[s].gcd_count>1)
 310                                 {
 311                                         LOG(cerr<<m_bubbles[s].n<<"("<<int(100*float(m_bubbles[s].n)/n)<<"%) "<<m_bubbles[s].f<<": finished score="<<m_bubbles[s].score<<" ("<<m_bubbles[s].err/m_bubbles[s].count<<","<<1.0-normm(m_bubbles[s].conv)<<") s="<<s<<" "<<m_bubbles[s].gcd_count<<" "<<m_bubbles[s].todrop<<" :"<<h2n(f2h(48000.0f/s))<<endl;)
 312                                 best_tot_c += m_bubbles[s].count;
 313
 314                                 }
 315                                 if(m_bubbles[s].gcd_count>=f)
 316                                         if(gcds<s)      gcds = s;
 317                         }
 318                 }
 319
 320                 LOG(
 321                 cerr<<n<<": f=" << f << " bests count="<<int(10000*float(best_tot_c)/n)/100.0f<<"% max n="<<max_n<<endl;
 322                 cerr<<"finished map size="<<map_size<<" ("<<int(log(float(map_size))+1)<<")" << endl;
 323
 324                 cerr<<"drops=[err="<<err_drop<<"("<<100*err_drop/init_map_size<<"%)";
 325                 if(err_drop>0)  cerr<<"(pos:"<<100*err_drop_t/err_drop/n<<"%)";
 326                 cerr<<",conv="<<conv_drop;
 327                 cerr<<",score="<<score_drop<<"("<<100*score_drop/init_map_size<<"%)";
 328                 if(score_drop>0)        cerr<<"(pos:"<<100*score_drop_t/score_drop/n<<"%)";
 329                 cerr<<",sgn="<<sgn_drop<<"("<<100*sgn_drop/init_map_size<<"%)";
 330                 if(sgn_drop>0)  cerr<<"(pos:"<<100*sgn_drop_t/sgn_drop/n<<"%)";
 331                 cerr<<",mult="<<mult_drop<<"("<<100*mult_drop/init_map_size<<"%)";
 332                 if(mult_drop>0) cerr<<"(pos:"<<100*mult_drop_t/mult_drop/n<<"%)";
 333                 cerr<<"]"<<endl;
 334                 )
 335
 336                 if(f==0)
 337                         m_wave_length = 0;
 338                 else
 339                 {
 340 //                      m_wave_length = (*m_best_set.begin()).second.front();
 341 //
 342 //                      if(m_bubbles[m_wave_length].err/m_bubbles[m_wave_length].count>err_threshold
 343 //                                      || normm(m_bubbles[m_wave_length].conv)<conv_threshold)
 344 //                              m_wave_length = 0;
 345
 346                         m_wave_length = gcds;
 347                 }
 348
 349                 LOG(cerr << "Final " << n << ": wave length=" << m_wave_length << endl;)
 350         }
 351 }
 352