CLAM-Development  1.4.0
CircularPeakPicking.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2001-2006 MUSIC TECHNOLOGY GROUP (MTG)
3  * UNIVERSITAT POMPEU FABRA
4  *
5  *
6  * This program is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 #ifndef CircularPeakPicking_hxx
23 #define CircularPeakPicking_hxx
24 #include <list>
25 #include <vector>
26 
27 namespace Simac
28 {
29 
46 {
47 public:
48  typedef std::vector<std::pair<double, double> > PeakList;
49 private:
50  PeakList _output;
51  unsigned int _maxSize;
52  double _binSize;
53  double _offset;
54 public:
55  CircularPeakPicking(unsigned chromagramSize, double binSize=1.0, double offset=0.0)
56  : _maxSize(chromagramSize), _binSize(binSize), _offset(offset)
57  {
58  _output.reserve(chromagramSize);
59  }
70  std::pair<double,double> interpolate(double y0, double y1, double y2)
71  {
72  // From the quadratic lagrange interpolation
73  // y= y0(x1-x)(x2-x)/(x1-x0)(x2-x0) +
74  // + y1(x-x0)(x2-x)/(x1-x0)(x2-x0) +
75  // + y1(x2-x)(x-x0)/(x2-x1)(x2-x0) +
76  // + y2(x-x1)(x-x0)/(x2-x1)(x2-x0) =
77  //
78  // considering x0=0, x1=1 and x2=2
79  //
80  // = y0(x-1)(x-2)/2 +
81  // - y1(x)(x-2)/2 +
82  // - y1(x-2)(x)/2 +
83  // + y2(x-1)(x)/2 =
84  //
85  // = y0(x-1)(x-2)/2 +
86  // - y1(x-2)(x) +
87  // + y2(x-1)(x)/2 =
88  //
89  // = y0/2 (x^2 - 3x + 2) - y1 (x^2-2x) + y2/2 (x^2-x)
90  // = (y0/2-y1+y2/2) x^2 + (-3*y0/2 + 2*y1 - y2/2) x + y0
91 
92  double a = y0/2 - y1 + y2/2;
93  double b = y1 -y0 -a; // = -3*y0/2 + 2*y1 -y2/2;
94  double c = y0;
95 
96  // From equating to zero the derivate of x*x*a + x*b + c
97  double xmax = -b/(a*2);
98  // ymax = xmax*xmax*a + b*xmax + c =
99  // = a*b*b/(4*a*a) -b*b/(2*a) + c =
100  // = b*b/(4*a) -b*b/(2*a) + c =
101  // = -b*b/(4*a) + c
102  double ymax = b*xmax/2 + y0;
103 
104  return std::make_pair(xmax, ymax);
105  }
106  void doIt(const std::vector<double> & chromagram)
107  {
108  _output.resize(0);
109  unsigned i0=_maxSize-2;
110  unsigned i1=_maxSize-1;
111  for (unsigned i=0; i<_maxSize; i0=i1, i1=i, i++)
112  {
113  // not equal to support plain two bins peaks
114  if (chromagram[i0] > chromagram[i1]) continue;
115  if (chromagram[i ] >= chromagram[i1]) continue;
116 
117  std::pair<double,double> interpoled=
118  interpolate(chromagram[i0],chromagram[i1],chromagram[i]);
119  // Adding the base bin
120  interpoled.first+=i0;
121  // Folding to [0,_maxSize) interval
122  while (interpoled.first<0) interpoled.first +=_maxSize;
123  while (interpoled.first>=_maxSize) interpoled.first -=_maxSize;
124  // Scaling
125  interpoled.first*=_binSize;
126  // Shifting to the first bin position
127  interpoled.first+=_offset;
128 
129  _output.push_back(interpoled);
130  }
131  }
132  const PeakList & output() const
133  {
134  return _output;
135  }
136 };
137 
138 } // namespace Simac
139 
140 #endif// CircularPeakPicking_hxx
141