1 ////////////////////////////////////////////////////////////////////////////////
2 ///
3 /// Peak detection routine.
4 ///
5 /// The routine detects highest value on an array of values and calculates the
6 /// precise peak location as a mass-center of the 'hump' around the peak value.
7 ///
8 /// Author : Copyright (c) Olli Parviainen
9 /// Author e-mail : oparviai 'at' iki.fi
10 /// SoundTouch WWW: http://www.surina.net/soundtouch
11 ///
12 ////////////////////////////////////////////////////////////////////////////////
13 //
14 // Last changed : $Date: 2015-05-18 18:22:02 +0300 (ma, 18 touko 2015) $
15 // File revision : $Revision: 4 $
16 //
17 // $Id: PeakFinder.cpp 213 2015-05-18 15:22:02Z oparviai $
18 //
19 ////////////////////////////////////////////////////////////////////////////////
20 //
21 // License :
22 //
23 // SoundTouch audio processing library
24 // Copyright (c) Olli Parviainen
25 //
26 // This library is free software; you can redistribute it and/or
27 // modify it under the terms of the GNU Lesser General Public
28 // License as published by the Free Software Foundation; either
29 // version 2.1 of the License, or (at your option) any later version.
30 //
31 // This library is distributed in the hope that it will be useful,
32 // but WITHOUT ANY WARRANTY; without even the implied warranty of
33 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
34 // Lesser General Public License for more details.
35 //
36 // You should have received a copy of the GNU Lesser General Public
37 // License along with this library; if not, write to the Free Software
38 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
39 //
40 ////////////////////////////////////////////////////////////////////////////////
41
42 #include <math.h>
43 #include <assert.h>
44
45 #include "PeakFinder.h"
46
47 using namespace soundtouch;
48
49 #define max(x, y) (((x) > (y)) ? (x) : (y))
50
51
PeakFinder()52 PeakFinder::PeakFinder()
53 {
54 minPos = maxPos = 0;
55 }
56
57
58 // Finds real 'top' of a peak hump from neighnourhood of the given 'peakpos'.
findTop(const float * data,int peakpos) const59 int PeakFinder::findTop(const float *data, int peakpos) const
60 {
61 int i;
62 int start, end;
63 float refvalue;
64
65 refvalue = data[peakpos];
66
67 // seek within �10 points
68 start = peakpos - 10;
69 if (start < minPos) start = minPos;
70 end = peakpos + 10;
71 if (end > maxPos) end = maxPos;
72
73 for (i = start; i <= end; i ++)
74 {
75 if (data[i] > refvalue)
76 {
77 peakpos = i;
78 refvalue = data[i];
79 }
80 }
81
82 // failure if max value is at edges of seek range => it's not peak, it's at slope.
83 if ((peakpos == start) || (peakpos == end)) return 0;
84
85 return peakpos;
86 }
87
88
89 // Finds 'ground level' of a peak hump by starting from 'peakpos' and proceeding
90 // to direction defined by 'direction' until next 'hump' after minimum value will
91 // begin
findGround(const float * data,int peakpos,int direction) const92 int PeakFinder::findGround(const float *data, int peakpos, int direction) const
93 {
94 int lowpos;
95 int pos;
96 int climb_count;
97 float refvalue;
98 float delta;
99
100 climb_count = 0;
101 refvalue = data[peakpos];
102 lowpos = peakpos;
103
104 pos = peakpos;
105
106 while ((pos > minPos+1) && (pos < maxPos-1))
107 {
108 int prevpos;
109
110 prevpos = pos;
111 pos += direction;
112
113 // calculate derivate
114 delta = data[pos] - data[prevpos];
115 if (delta <= 0)
116 {
117 // going downhill, ok
118 if (climb_count)
119 {
120 climb_count --; // decrease climb count
121 }
122
123 // check if new minimum found
124 if (data[pos] < refvalue)
125 {
126 // new minimum found
127 lowpos = pos;
128 refvalue = data[pos];
129 }
130 }
131 else
132 {
133 // going uphill, increase climbing counter
134 climb_count ++;
135 if (climb_count > 5) break; // we've been climbing too long => it's next uphill => quit
136 }
137 }
138 return lowpos;
139 }
140
141
142 // Find offset where the value crosses the given level, when starting from 'peakpos' and
143 // proceeds to direction defined in 'direction'
findCrossingLevel(const float * data,float level,int peakpos,int direction) const144 int PeakFinder::findCrossingLevel(const float *data, float level, int peakpos, int direction) const
145 {
146 float peaklevel;
147 int pos;
148
149 peaklevel = data[peakpos];
150 assert(peaklevel >= level);
151 pos = peakpos;
152 while ((pos >= minPos) && (pos < maxPos))
153 {
154 if (data[pos + direction] < level) return pos; // crossing found
155 pos += direction;
156 }
157 return -1; // not found
158 }
159
160
161 // Calculates the center of mass location of 'data' array items between 'firstPos' and 'lastPos'
calcMassCenter(const float * data,int firstPos,int lastPos) const162 double PeakFinder::calcMassCenter(const float *data, int firstPos, int lastPos) const
163 {
164 int i;
165 float sum;
166 float wsum;
167
168 sum = 0;
169 wsum = 0;
170 for (i = firstPos; i <= lastPos; i ++)
171 {
172 sum += (float)i * data[i];
173 wsum += data[i];
174 }
175
176 if (wsum < 1e-6) return 0;
177 return sum / wsum;
178 }
179
180
181
182 /// get exact center of peak near given position by calculating local mass of center
getPeakCenter(const float * data,int peakpos) const183 double PeakFinder::getPeakCenter(const float *data, int peakpos) const
184 {
185 float peakLevel; // peak level
186 int crosspos1, crosspos2; // position where the peak 'hump' crosses cutting level
187 float cutLevel; // cutting value
188 float groundLevel; // ground level of the peak
189 int gp1, gp2; // bottom positions of the peak 'hump'
190
191 // find ground positions.
192 gp1 = findGround(data, peakpos, -1);
193 gp2 = findGround(data, peakpos, 1);
194
195 peakLevel = data[peakpos];
196
197 if (gp1 == gp2)
198 {
199 // avoid rounding errors when all are equal
200 assert(gp1 == peakpos);
201 cutLevel = groundLevel = peakLevel;
202 } else {
203 // get average of the ground levels
204 groundLevel = 0.5f * (data[gp1] + data[gp2]);
205
206 // calculate 70%-level of the peak
207 cutLevel = 0.70f * peakLevel + 0.30f * groundLevel;
208 }
209
210 // find mid-level crossings
211 crosspos1 = findCrossingLevel(data, cutLevel, peakpos, -1);
212 crosspos2 = findCrossingLevel(data, cutLevel, peakpos, 1);
213
214 if ((crosspos1 < 0) || (crosspos2 < 0)) return 0; // no crossing, no peak..
215
216 // calculate mass center of the peak surroundings
217 return calcMassCenter(data, crosspos1, crosspos2);
218 }
219
220
221
detectPeak(const float * data,int aminPos,int amaxPos)222 double PeakFinder::detectPeak(const float *data, int aminPos, int amaxPos)
223 {
224
225 int i;
226 int peakpos; // position of peak level
227 double highPeak, peak;
228
229 this->minPos = aminPos;
230 this->maxPos = amaxPos;
231
232 // find absolute peak
233 peakpos = minPos;
234 peak = data[minPos];
235 for (i = minPos + 1; i < maxPos; i ++)
236 {
237 if (data[i] > peak)
238 {
239 peak = data[i];
240 peakpos = i;
241 }
242 }
243
244 // Calculate exact location of the highest peak mass center
245 highPeak = getPeakCenter(data, peakpos);
246 peak = highPeak;
247
248 // Now check if the highest peak were in fact harmonic of the true base beat peak
249 // - sometimes the highest peak can be Nth harmonic of the true base peak yet
250 // just a slightly higher than the true base
251
252 for (i = 3; i < 10; i ++)
253 {
254 double peaktmp, harmonic;
255 int i1,i2;
256
257 harmonic = (double)i * 0.5;
258 peakpos = (int)(highPeak / harmonic + 0.5f);
259 if (peakpos < minPos) break;
260 peakpos = findTop(data, peakpos); // seek true local maximum index
261 if (peakpos == 0) continue; // no local max here
262
263 // calculate mass-center of possible harmonic peak
264 peaktmp = getPeakCenter(data, peakpos);
265
266 // accept harmonic peak if
267 // (a) it is found
268 // (b) is within �4% of the expected harmonic interval
269 // (c) has at least half x-corr value of the max. peak
270
271 double diff = harmonic * peaktmp / highPeak;
272 if ((diff < 0.96) || (diff > 1.04)) continue; // peak too afar from expected
273
274 // now compare to highest detected peak
275 i1 = (int)(highPeak + 0.5);
276 i2 = (int)(peaktmp + 0.5);
277 if (data[i2] >= 0.4*data[i1])
278 {
279 // The harmonic is at least half as high primary peak,
280 // thus use the harmonic peak instead
281 peak = peaktmp;
282 }
283 }
284
285 return peak;
286 }
287