libs/Music/CombedFT.cpp

   1 // Copyright 2007 "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 #include "CombedFT.h"
  20
  21 #include <assert.h>
  22 #include <iostream>
  23 using namespace std;
  24 using namespace Math;
  25 #include "Music.h"
  26 #include "SPWindow.h"
  27 #include "FreqAnalysis.h"
  28
  29 namespace Music
  30 {
  31         CombedFT::CombedFT()
  32         {
  33                 m_use_audibility_treshold = false;
  34                 m_audib_ratio = 0.1;
  35
  36                 m_zp_factor = 1.0;
  37                 m_window_factor = 1.0;
  38
  39                 init();
  40         }
  41         void CombedFT::setZeroPaddingFactor(double zp)
  42         {
  43                 if(zp!=m_zp_factor)
  44                 {
  45                         m_zp_factor = zp;
  46                         init();
  47                 }
  48         }
  49         void CombedFT::setWindowFactor(double wf)
  50         {
  51                 if(wf!=m_window_factor)
  52                 {
  53                         m_window_factor = wf;
  54                         init();
  55                 }
  56         }
  57
  58         void CombedFT::init()
  59         {
  60                 if(GetSamplingRate()<=0)        return;
  61
  62                 m_f0 = 0.0;
  63
  64                 int win_size = int(m_window_factor*GetSamplingRate()/h2f(GetSemitoneMin()));// at least m_window_factor period of the lowest freq
  65
  66                 int best_size = 2;
  67                 while(best_size<int(win_size))
  68                         best_size *= 2;
  69
  70                 m_win = hann(best_size);
  71
  72                 while(best_size<int(m_zp_factor*win_size))
  73                         best_size *= 2;
  74
  75                 cerr << "CombedFT: INFO: window size=" << win_size << " FFT size=" << best_size << " window size factor=" << m_window_factor << " zero padding factor=" << m_zp_factor << endl;
  76
  77                 m_plan.resize(best_size);
  78                 m_components.resize(m_plan.size()/2);
  79                 m_comb.resize(m_plan.size()/2);
  80         }
  81         int CombedFT::getSampleAlgoLatency() const
  82         {
  83                 return 1000*m_win.size()/GetSamplingRate();
  84         }
  85         int CombedFT::getMinSize() const
  86         {
  87                 return m_win.size();
  88         }
  89         double CombedFT::getFondamentalFreq() const
  90         {
  91                 return m_f0;
  92         }
  93
  94         void CombedFT::apply(const deque<double>& buff)
  95         {
  96 //              cerr << getAmplitudeTreshold() << " " << getComponentTreshold() << " " << m_audib_ratio << endl;
  97
  98                 if(int(buff.size())<getMinSize() || m_win.empty())      return;
  99
 100                 m_max_amplitude = 0.0;
 101                 for(size_t i=0; i<m_win.size(); i++)
 102                 {
 103                         m_plan.in[i] = buff[i]*m_win[i];
 104
 105                         if(abs(buff[i])>m_max_amplitude)
 106                                 m_max_amplitude = abs(buff[i]);
 107                 }
 108                 for(int i=m_win.size(); i<m_plan.size(); i++)                   // padd with zeros
 109                         m_plan.in[i] = 0.0;
 110
 111                 m_plan.execute();
 112
 113                 for(size_t i=0; i<m_comb.in.size() && i<m_plan.out.size(); i++)
 114                         m_comb.in[i] = mod(m_plan.out[i]);
 115
 116                 // compute max with respect of the bounds
 117                 m_components_max = 0.0;
 118                 int max_index = -1;
 119                 double fmin = h2f(GetSemitoneMin());
 120                 double fmax = h2f(GetSemitoneMax());
 121                 for(size_t i=0; i<m_components.size(); i++)
 122                 {
 123                         double fi = i*double(GetSamplingRate())/m_plan.size();
 124                         if(fmin<=fi && fi<=fmax
 125                                                 && m_comb.in[i]>m_components_max)
 126                         {
 127                                 max_index = i;
 128                                 m_components_max = m_comb.in[i];
 129                         }
 130                 }
 131
 132                 m_f0 = 0.0;
 133
 134                 if(m_components_max>getComponentTreshold())
 135                 {
 136                         m_f0 = PeakRefinementLogParabola(m_plan.out, max_index)*double(GetSamplingRate())/m_plan.size();
 137
 138                         // TODO TEST *win[i]; // me semble qu'une trans de harm signal n'est pas trop discontinue aux extrémités
 139                         m_comb.execute();
 140
 141                         for(size_t i=0; i<m_components.size(); i++)
 142                                 m_components[i] = 0.0;
 143
 144                         for(size_t i=0; i<m_comb.out.size()/2; i++)
 145                                 m_components[i] = real(m_comb.out[i]);
 146
 147                         // avg is not interesting
 148                         m_components[0] = 0.0;
 149
 150                         // frequencies between 1 and g_res_factor are virtuals
 151 //                      for(int i=m_comb.out.size()-(g_res_factor-1); i<m_comb.out.size(); i++)
 152 //                              m_components[i] = 0.0;
 153
 154                         double step = (m_plan.size()/max_index)/2;
 155                         // hyp: the fund freq is not greater than the max amplitude harmonic
 156                         // keep only multiples of the max amplitude harmonic
 157                         for(int i=1; i<int(m_components.size())/step; i++)
 158                                 m_components[i] = m_components[int(step*i)];
 159                         for(int i=int(m_components.size()/step); i<int(m_components.size()); i++)
 160                                 m_components[i] = 0.0;
 161
 162                         if(m_use_audibility_treshold)
 163                         {
 164                                 for(int i=1; i<int(m_components.size()/step); i++)
 165                                         m_components[i] /= (1+m_audib_ratio*(i-1));
 166                         }
 167                         else
 168                         {
 169                                 vector<double> temp_comp(int(m_components.size()/step), 0.0);
 170                                 for(int i=1; i<int(temp_comp.size()); i++)
 171                                         temp_comp[i] =  (m_components[i]-m_components[i-1]) +
 172                                                                         (m_components[i]-m_components[i+1]);
 173                                 for(int i=1; i<int(temp_comp.size()); i++)
 174                                         m_components[i] = temp_comp[i];
 175                         }
 176
 177                         // find the max
 178                         double max_amp = 0.0;
 179                         max_index = -1;
 180                         for(size_t i=0; i<m_components.size(); i++)
 181                         {
 182 //                              double fi = i*double(GetSamplingRate())/m_win.size();
 183                                 if(m_components[i]>max_amp)
 184                                 {
 185                                         max_index = i;
 186                                         max_amp = m_components[i];
 187                                 }
 188                         }
 189
 190                         if(max_index>0)
 191                                 m_f0 /= max_index;
 192                 }
 193
 194 //              cerr << " final: " << GetSamplingRate() << ":" << m_f0 << endl;
 195         }
 196         CombedFT::~CombedFT()
 197         {
 198         }
 199 }