libs/Music/MultiCumulativeDiffAlgo.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 "MultiCumulativeDiffAlgo.h"
  21
  22 #include <cassert>
  23 #include <cmath>
  24 #include <iostream>
  25 #include <limits>
  26 using namespace std;
  27 #include <CppAddons/CAMath.h>
  28 using namespace Math;
  29
  30 #include "Music.h"
  31
  32 //#define MUSIC_DEBUG
  33 #ifdef MUSIC_DEBUG
  34 #define LOG(a)  a
  35 #else
  36 #define LOG(a)
  37 #endif
  38
  39 namespace Music
  40 {
  41         void MultiCumulativeDiffAlgo::init()
  42         {
  43                 for(size_t i=0; i<size(); i++)
  44                 {
  45                         if(m_diffs[i]!=NULL)    delete m_diffs[i];
  46
  47                         m_diffs[i] = new CumulativeDiff(1, i+GetSemitoneMin());
  48                 }
  49         }
  50
  51         MultiCumulativeDiffAlgo::MultiCumulativeDiffAlgo(int latency_factor, double test_complexity)
  52         {
  53                 assert(GetSamplingRate()>0);
  54
  55                 m_test_complexity = test_complexity;
  56
  57                 m_diffs.resize(size());
  58                 for(size_t i=0; i<size(); i++)
  59                         m_diffs[i] = NULL;
  60                 init();
  61         }
  62         bool MultiCumulativeDiffAlgo::is_minima(int ih)
  63         {
  64                 if(ih+1>=0 && ih+1<int(size()))
  65                         if(m_components[ih+1]<=m_components[ih])
  66                                 return false;
  67                 if(ih-1>=0 && ih-1<int(size()))
  68                         if(m_components[ih-1]<=m_components[ih])
  69                                 return false;
  70                 return true;
  71         }
  72         double MultiCumulativeDiffAlgo::getFondamentalWaveLength() const
  73         {
  74                 return int(GetSamplingRate()/h2f(m_first_fond+GetSemitoneMin(), GetAFreq()));
  75         }
  76         void MultiCumulativeDiffAlgo::apply(const deque<double>& buff)
  77         {
  78                 assert(GetSamplingRate()>0);
  79                 for(size_t i=0; i<size(); i++)
  80                 {
  81                         m_components[i] = 1.0;
  82                         m_is_fondamental[i] = false;
  83                 }
  84
  85                 m_first_fond = -1;
  86
  87                 if(buff.empty() || buff.size()<(m_test_complexity+1)*m_diffs[0]->m_s)
  88                         return;
  89
  90                 double v = 0.0;
  91                 for(size_t i=0; i<buff.size() && v<=getAmplitudeTreshold() && i<m_diffs[0]->m_s; i++)
  92                         v = max(v, abs(buff[i]));
  93
  94                 if(v>getAmplitudeTreshold())
  95                 {
  96                         // compute all components
  97                         m_components_max = 0.0;
  98                         double min_comp = 1000000;
  99                         double max_sum = 0.0;
 100                         for(int ih=int(size())-1; ih>=0; ih--)
 101                         {
 102                                 m_diffs[ih]->receive(buff, 0);
 103                                 m_components[ih] = m_diffs[ih]->m_error;
 104                                 m_components_max = max(m_components_max, m_components[ih]);
 105                         }
 106
 107                         // test components
 108                         for(int ih=int(size())-1; ih>=0; ih--)
 109                         {
 110                                 bool ok = true;
 111
 112                                 bool crit_min = true;
 113                                 // criteria: the fond and his first harmonics are minimas
 114                                 if(ok)  ok =
 115                                                         is_minima(ih) &&
 116                                                         is_minima(ih-12) &&
 117                                                         (is_minima(ih-19) || is_minima(ih-19-1) || is_minima(ih-19+1)) &&
 118                                                         is_minima(ih-24);
 119
 120                                 crit_min = ok;
 121
 122                                 bool crit_small_enough = true;
 123                                 if(ok)
 124                                 {
 125                                         if(m_components[ih]/m_components_max>getComponentTreshold())
 126                                                 crit_small_enough = false;
 127                                         ok = crit_small_enough;
 128                                 }
 129
 130 //                              bool crit_cross_zero = true;
 131                                 /*
 132                                 // criteria: wave should cross the zero value
 133                                 if(ok)
 134                                 {
 135                                         bool cross = true;
 136                                         for(int s=0; cross && s<int(m_diffs[ih]->m_s); s++)
 137                                         {
 138                                                 double sg = Math::sgn(buff[s]);
 139                                                 bool same_side = true;
 140                                                 for(int i=0; same_side && i<int(m_diffs[ih]->m_s); i++)
 141                                                         same_side = Math::sgn(buff[s+i])==sg;
 142
 143                                                 cross = !same_side;
 144                                         }
 145
 146                                         ok = crit_cross_zero = cross;
 147                                 }
 148                                 */
 149
 150 //                              bool crit_integ_cst = true;
 151                                 // criteria: integral should be nearly constant while shifting
 152                                 // TODO
 153                                 // (the previous criteria seems sufficient to remove high comp.)
 154
 155 //                              LOG(if(crit_min)
 156 //                                      cerr << "ih=" << ih <<
 157 //                                      " harm_min=(("<<is_minima(ih-12)<<","<<is_minima(ih-12-1)<<","<<is_minima(ih-12+1)<<"),("<<
 158 //                                              is_minima(ih-19)<<","<<is_minima(ih-19-1)<<","<<is_minima(ih-19+1)<<"),("<<
 159 //                                              is_minima(ih-24)<<","<<is_minima(ih-24-1)<<","<<is_minima(ih-24+1)<<"))"<< crit_min <<
 160 //                                      " increase=" << crit_incr_err <<
 161 //                                      " cross_zero=" << crit_cross_zero <<
 162 //                                      " integ_cst=" << crit_integ_cst <<
 163 //                                      " ok=" << ok << " c=" << m_components[ih] << endl;)
 164
 165                                 // if all criteria are ok
 166                                 if(ok)
 167                                 {
 168                                         double sum = 0.0;
 169                                         int n=0;
 170                                         int i=0;
 171                                         double wh = 1.0;
 172                                         sum += wh*(m_components_max-m_components[ih]); n++;
 173                                         i=12;   if(ih-i>=0)     {sum+=wh*(m_components_max-m_components[ih-i]); n++;}
 174                                         i=19;   if(ih-i>=0)     {sum+=wh*(m_components_max-m_components[ih-i]); n++;}
 175                                         i=24;   if(ih-i>=0)     {sum+=wh*(m_components_max-m_components[ih-i]); n++;}
 176
 177                                         LOG(cerr << "ih=" << ih << " sum=" << sum << endl;)
 178
 179                                         // get the "best"
 180                                         if(sum>max_sum)
 181                                         {
 182                                                 size_t step = size_t(m_diffs[ih]->m_s/m_test_complexity);
 183                                                 if(step<1)      step = 1;
 184                                                 for(size_t s=0; ok && s<m_diffs[ih]->m_s; s+=step)
 185                                                 {
 186                                                         if(ih-1>=0){
 187                                                                 m_diffs[ih-1]->receive(buff, s);
 188                                                                 m_components[ih-1] = m_diffs[ih-1]->m_error;
 189                                                         }
 190                                                         if(ih+1<int(size())){
 191                                                                 m_diffs[ih+1]->receive(buff, s);
 192                                                                 m_components[ih+1] = m_diffs[ih+1]->m_error;
 193                                                         }
 194                                                         m_diffs[ih]->receive(buff, s);
 195                                                         m_components[ih] = m_diffs[ih]->m_error;
 196                                                         ok = is_minima(ih);
 197                                                 }
 198
 199                                                 if(ok)
 200                                                 {
 201                                                         max_sum = sum;
 202                                                         min_comp = m_components[ih];
 203                                                         m_first_fond = ih;
 204                                                 }
 205                                         }
 206                                 }
 207                         }
 208
 209 //                      cerr << "ff: " << m_first_fond << endl;
 210
 211                         if(m_first_fond!=-1)
 212                                 m_is_fondamental[m_first_fond] = true;
 213
 214                         LOG(cerr << "m_first_fond=" << m_first_fond << endl;)
 215                 }
 216         }
 217         MultiCumulativeDiffAlgo::~MultiCumulativeDiffAlgo()
 218         {
 219                 for(size_t i=0; i<m_diffs.size(); i++)
 220                         delete m_diffs[i];
 221         }
 222 }