libs/Music/TimeAnalysis.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 "TimeAnalysis.h"
  21
  22 #include <cassert>
  23 #include <cmath>
  24 #include <deque>
  25 #include <iostream>
  26 #include <limits>
  27 using namespace std;
  28 #include <CppAddons/CAMath.h>
  29 using namespace Math;
  30
  31 #include "Music.h"
  32
  33 namespace Music
  34 {
  35         double InterpolatedPeriod(const std::deque<double>& queue, int left, int right)
  36         {
  37                 double l = left - queue[left]/(queue[left+1]-queue[left]);
  38
  39                 double r = right - queue[right]/(queue[right+1]-queue[right]);
  40
  41                 return r - l;
  42         }
  43
  44         /*
  45          * on peut imaginer des cas qui mettent en échec cette procédure:
  46          *      on selectionne un zéro qui n'en n'est pas un une periode plus
  47          *      tard et si un autre zéro se trouve dans la zone de tolérance la longeur
  48          *      ainsi calculée entre ces deux zéro (qui ne se correspondent donc pas) sera fausse.
  49          *      example: une fréquence très basse avec une seule harmonique très très
  50          *      haute.
  51          *      - il faut utiliser des zÃ©ros significatifs ... et ... et ... et voilà .
  52          *      - ou encore écarter les solutions trop élognées de la moyenne
  53          */
  54         double GetAveragePeriodFromApprox(const std::deque<double>& queue, int approx, int n)
  55         {
  56                 if(GetAFreq()<=0.0 || GetSamplingRate()<=0.0 || int(queue.size())<approx)
  57                         return 0.0;
  58
  59                 deque<int> ups;                                                                 // the upper peeks
  60
  61                 // parse the whole buffer, for n zeros
  62                 for(int i=0; int(ups.size())<n && i+1<int(queue.size()); i++)
  63                         if(queue[i]<0 && queue[i+1]>0)                          // if it cross the axis
  64                                 ups.push_back(i);
  65
  66 //              cerr << "approx=" << approx << " ups.size()=" << ups.size();
  67                 if(ups.empty())
  68                         return 0.0;
  69
  70                 double ht = f2hf(double(GetSamplingRate())/approx);
  71                 int period_low_bound = int(GetSamplingRate()/h2f(ht+1))-2;
  72                 int period_high_bound = int(GetSamplingRate()/h2f(ht-1))+2;
  73
  74 //              cerr << " ht=" << ht << " lb=" << period_low_bound << " rb=" << period_high_bound;
  75
  76 //              cerr << " periods=(";
  77
  78                 double period = 0.0;
  79                 int count = 0;
  80
  81                 for(int i=0; i<int(ups.size()) && count<n; i++)
  82                 {
  83                         int i_seek = ups[i] + approx;
  84
  85                         int lower_i_seek = i_seek;
  86                         int low_bound = ups[i] + period_low_bound;
  87                         int higher_i_seek = i_seek;
  88                         int high_bound = min(int(queue.size())-1, ups[i]+period_high_bound);
  89
  90 //                      cerr << "{" << low_bound << ":" << i_seek << ":" << high_bound << "}";
  91
  92                         if(low_bound+1>=int(queue.size()))
  93                                 i = ups.size();                                                                         // stop loop
  94                         else
  95                         {
  96                                 if(!(queue[i_seek]<=0.0 && queue[i_seek+1]>0.0))
  97                                 {
  98                                         while(lower_i_seek>low_bound &&
  99                                                 !(queue[lower_i_seek]<=0.0 && queue[lower_i_seek+1]>0.0))
 100                                                 lower_i_seek--;
 101
 102                                         while(higher_i_seek<high_bound &&
 103                                                 !(queue[higher_i_seek]<=0.0 && queue[higher_i_seek+1]>0.0))
 104                                                 higher_i_seek++;
 105
 106                                         if(i_seek-lower_i_seek < higher_i_seek-i_seek)  // take the nearest to i_seek
 107                                                 i_seek = lower_i_seek;
 108                                         else
 109                                                 i_seek = higher_i_seek;
 110                                 }
 111
 112 //                              cerr << i_seek << "=>";
 113
 114                                 if(low_bound<i_seek && i_seek<high_bound)
 115                                 {
 116                                         double per = InterpolatedPeriod(queue, ups[i], i_seek);
 117
 118 //                                      cerr << "f=" << GetSamplingRate()/per << " ";
 119
 120                                         period += per;
 121                                         count++;
 122                                 }
 123                         }
 124                 }
 125
 126                 if(count==0)
 127                         return 0.0;
 128
 129                 period /= count;
 130
 131 //              cerr << ")=" << GetSamplingRate()/period << "(" << count << ")" << endl;
 132
 133                 return period;
 134         }
 135
 136         void GetWaveSample(const std::deque<double>& queue, size_t wave_length, std::deque<double>& sample)
 137         {
 138                 assert(wave_length>0);
 139                 if(queue.size()<2*wave_length)  return;
 140
 141                 // find the highest peek in the second period
 142                 int left = 0;
 143                 double min_vol = 0;
 144                 for(int i=int(wave_length); i<int(queue.size()) && i<int(2*wave_length); i++)
 145                 {
 146                         if(queue[i]<min_vol)
 147                         {
 148                                 min_vol = queue[i];
 149                                 left = i;
 150                         }
 151                 }
 152
 153                 // adjust the right bound to the nearest minima
 154                 int left_right = int(left + 0.9*wave_length);
 155                 int right_right = int(left + wave_length/0.9);
 156
 157                 int right = int(left + wave_length);                            // init to a default value
 158
 159                 if(right_right>=int(queue.size()))
 160                         return;
 161
 162                 min_vol = 0.0;
 163                 for(int i=left_right; i<=right_right; i++)
 164                 {
 165                         if(queue[i]<min_vol)
 166                         {
 167                                 min_vol = queue[i];
 168                                 right = i;
 169                         }
 170                 }
 171
 172                 // fill in the sample
 173                 sample.clear();
 174                 for(int i=left; i<int(queue.size()) && i<right; i++)
 175                         sample.push_back(queue[i]);
 176         }
 177
 178         double GetAverageWaveLengthFromApproxEnergy(const std::deque<double>& queue, double approx, int n)
 179         {
 180                 assert(GetSamplingRate()>0);
 181
 182                 if(queue.size()<approx*1.5)
 183                         return 0.0;
 184
 185 //              cerr << queue.size() << "=>" << approx << " n=" << n << endl;
 186
 187                 double wave_length = 0.0;
 188                 int count = 0;
 189                 int seek = 0;
 190                 while(count<n && seek<int(queue.size()) && seek!=-1)
 191                 {
 192 //                      cerr << "new " << flush;
 193
 194                         // in one period, compute the energy over approx/4
 195                         int w = int(approx/4);          // TODO ptr un peu long
 196                         if(w<4) w=4;
 197                         vector<double> en;
 198                         for(int i=0; i<w/2; i++)
 199                                 en.push_back(0.0);
 200                         for(int i=w/2; int(en.size())<int(1.25*approx); i++)
 201                         {
 202                                 en.push_back(0.0);
 203                                 for(int j=-w/2; j<=w/2 && seek+i+j<int(queue.size()); j++)
 204                                         en.back() += queue[seek+i+j];   //*queue[seek+i+j]
 205                         }
 206
 207                         // find the highest energy peak
 208                         int i_max = 0;
 209                         double en_max = 0.0;
 210                         for(int i=0; i<int(en.size()); i++)
 211                         {
 212                                 if(en[i]>en_max)
 213                                 {
 214                                         en_max = en[i];
 215                                         i_max = i;
 216                                 }
 217                         }
 218                         seek += i_max;
 219
 220 //                      cerr << "max-seek=" << seek << " " << flush;
 221
 222                         int old_seek=seek;
 223                         // go back to the previous zero
 224                         while(seek>=0 && !(queue[seek]<=0 && queue[seek+1]>0) && seek>old_seek-approx/2)
 225                                 seek--;
 226
 227 //                      cerr << "zero-seek=" << seek << " " << flush;
 228
 229                         if(seek<0 || seek<=old_seek-approx/2)
 230                         {
 231                                 seek += int(approx);
 232 //                              cerr << endl;
 233                                 continue;
 234                         }
 235
 236                         int left = seek;
 237                         int right = int(left + approx);
 238                         int sright = right;
 239                         int downlimit = int(left + 0.75*approx);
 240                         int bright = right;
 241                         int uplimit = int(left + 1.25*approx);
 242
 243                         // look for the nearest zero of right
 244                         while(sright+1<int(queue.size()) && sright>downlimit &&
 245                                                    !(queue[sright]<=0 && queue[sright+1]>0))
 246                                 sright--;
 247                         while(bright+1<int(queue.size()) && bright<uplimit &&
 248                                                    !(queue[bright]<=0 && queue[bright+1]>0))
 249                                 bright++;
 250
 251                         if(sright>=int(queue.size()) || bright>=int(queue.size()))
 252                         {
 253                                 seek = -1;
 254 //                              cerr << endl;
 255                                 continue;
 256                         }
 257
 258                         double sw = InterpolatedPeriod(queue, left, sright);
 259                         double bw = InterpolatedPeriod(queue, left, bright);
 260                         // keep the nearest one after approx
 261                         double wl = 0.0;
 262                         if(abs(sw-approx)<abs(bw-approx))       wl = sw;
 263                         else                                                            wl = bw;
 264
 265 //                      cerr << "wl=" << wl << flush;
 266
 267                         wave_length += wl;
 268                         count++;
 269
 270                         seek += int(0.9*approx);
 271
 272 //                      cerr << endl;
 273                 }
 274
 275 //              cerr << "("<<count<<")"<< flush;
 276
 277                 if(count==0)    return 0.0;
 278
 279                 wave_length /= count;
 280
 281 //              cerr << GetSamplingRate()/wave_length << endl;
 282
 283                 return wave_length;
 284         }
 285 }