1 /*
2 TiMidity++ -- MIDI to WAVE converter and player
3 Copyright (C) 1999-2002 Masanao Izumo <mo@goice.co.jp>
4 Copyright (C) 1995 Tuukka Toivonen <tt@cgs.fi>
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 #ifdef HAVE_CONFIG_H
22 #include "config.h"
23 #endif /* HAVE_CONFIG_H */
24
25 #include <stdio.h>
26 #include <stdlib.h>
27 #include <math.h>
28 #include "timidity.h"
29 #include "common.h"
30 #include "instrum.h"
31 #include "freq.h"
32 #include "fft4g.h"
33
34 float *floatdata, *magdata, *prunemagdata;
35 int *ip;
36 float *w;
37 uint32 oldfftsize = 0;
38 float pitchmags[129];
39 double pitchbins[129];
40 double new_pitchbins[129];
41 int *fft1_bin_to_pitch;
42
43 /* middle C = pitch 60 = 261.6 Hz
44 freq = 13.75 * exp((pitch - 9) / 12 * log(2))
45 pitch = 9 - log(13.75 / freq) * 12 / log(2)
46 = -36.37631656 + 17.31234049 * log(freq)
47 */
48 float pitch_freq_table[129] = {
49 8.17579892, 8.66195722, 9.17702400, 9.72271824, 10.3008612, 10.9133822,
50 11.5623257, 12.2498574, 12.9782718, 13.7500000, 14.5676175, 15.4338532,
51
52 16.3515978, 17.3239144, 18.3540480, 19.4454365, 20.6017223, 21.8267645,
53 23.1246514, 24.4997147, 25.9565436, 27.5000000, 29.1352351, 30.8677063,
54
55 32.7031957, 34.6478289, 36.7080960, 38.8908730, 41.2034446, 43.6535289,
56 46.2493028, 48.9994295, 51.9130872, 55.0000000, 58.2704702, 61.7354127,
57
58 65.4063913, 69.2956577, 73.4161920, 77.7817459, 82.4068892, 87.3070579,
59 92.4986057, 97.9988590, 103.826174, 110.000000, 116.540940, 123.470825,
60
61 130.812783, 138.591315, 146.832384, 155.563492, 164.813778, 174.614116,
62 184.997211, 195.997718, 207.652349, 220.000000, 233.081881, 246.941651,
63
64 261.625565, 277.182631, 293.664768, 311.126984, 329.627557, 349.228231,
65 369.994423, 391.995436, 415.304698, 440.000000, 466.163762, 493.883301,
66
67 523.251131, 554.365262, 587.329536, 622.253967, 659.255114, 698.456463,
68 739.988845, 783.990872, 830.609395, 880.000000, 932.327523, 987.766603,
69
70 1046.50226, 1108.73052, 1174.65907, 1244.50793, 1318.51023, 1396.91293,
71 1479.97769, 1567.98174, 1661.21879, 1760.00000, 1864.65505, 1975.53321,
72
73 2093.00452, 2217.46105, 2349.31814, 2489.01587, 2637.02046, 2793.82585,
74 2959.95538, 3135.96349, 3322.43758, 3520.00000, 3729.31009, 3951.06641,
75
76 4186.00904, 4434.92210, 4698.63629, 4978.03174, 5274.04091, 5587.65170,
77 5919.91076, 6271.92698, 6644.87516, 7040.00000, 7458.62018, 7902.13282,
78
79 8372.01809, 8869.84419, 9397.27257, 9956.06348, 10548.0818, 11175.3034,
80 11839.8215, 12543.8540, 13289.7503
81 };
82
83
84
85 /* center_pitch + 0.49999 */
86 float pitch_freq_ub_table[129] = {
87 8.41536325, 8.91576679, 9.44592587, 10.0076099, 10.6026933, 11.2331623,
88 11.9011208, 12.6087983, 13.3585565, 14.1528976, 14.9944727, 15.8860904,
89 16.8307265, 17.8315336, 18.8918517, 20.0152197, 21.2053866, 22.4663245,
90 23.8022417, 25.2175966, 26.7171129, 28.3057952, 29.9889453, 31.7721808,
91 33.6614530, 35.6630672, 37.7837035, 40.0304394, 42.4107732, 44.9326490,
92 47.6044834, 50.4351932, 53.4342259, 56.6115903, 59.9778907, 63.5443616,
93 67.3229060, 71.3261343, 75.5674070, 80.0608788, 84.8215464, 89.8652980,
94 95.2089667, 100.870386, 106.868452, 113.223181, 119.955781, 127.088723,
95 134.645812, 142.652269, 151.134814, 160.121758, 169.643093, 179.730596,
96 190.417933, 201.740773, 213.736904, 226.446361, 239.911563, 254.177446,
97 269.291624, 285.304537, 302.269628, 320.243515, 339.286186, 359.461192,
98 380.835867, 403.481546, 427.473807, 452.892723, 479.823125, 508.354893,
99 538.583248, 570.609074, 604.539256, 640.487030, 678.572371, 718.922384,
100 761.671734, 806.963092, 854.947614, 905.785445, 959.646250, 1016.70979,
101 1077.16650, 1141.21815, 1209.07851, 1280.97406, 1357.14474, 1437.84477,
102 1523.34347, 1613.92618, 1709.89523, 1811.57089, 1919.29250, 2033.41957,
103 2154.33299, 2282.43630, 2418.15702, 2561.94812, 2714.28948, 2875.68954,
104 3046.68693, 3227.85237, 3419.79046, 3623.14178, 3838.58500, 4066.83914,
105 4308.66598, 4564.87260, 4836.31405, 5123.89624, 5428.57897, 5751.37907,
106 6093.37387, 6455.70474, 6839.58092, 7246.28356, 7677.17000, 8133.67829,
107 8617.33197, 9129.74519, 9672.62809, 10247.7925, 10857.1579, 11502.7581,
108 12186.7477, 12911.4095, 13679.1618
109 };
110
111
112
113 /* center_pitch - 0.49999 */
114 float pitch_freq_lb_table[129] = {
115 7.94305438, 8.41537297, 8.91577709, 9.44593678, 10.0076214, 10.6027055,
116 11.2331752, 11.9011346, 12.6088129, 13.3585719, 14.1529139, 14.9944900,
117 15.8861088, 16.8307459, 17.8315542, 18.8918736, 20.0152428, 21.2054111,
118 22.4663505, 23.8022692, 25.2176258, 26.7171438, 28.3058279, 29.9889800,
119 31.7722175, 33.6614919, 35.6631084, 37.7837471, 40.0304857, 42.4108222,
120 44.9327009, 47.6045384, 50.4352515, 53.4342876, 56.6116557, 59.9779599,
121 63.5444350, 67.3229838, 71.3262167, 75.5674943, 80.0609713, 84.8216444,
122 89.8654018, 95.2090767, 100.870503, 106.868575, 113.223311, 119.955920,
123 127.088870, 134.645968, 142.652433, 151.134989, 160.121943, 169.643289,
124 179.730804, 190.418153, 201.741006, 213.737151, 226.446623, 239.911840,
125 254.177740, 269.291935, 285.304867, 302.269977, 320.243885, 339.286578,
126 359.461607, 380.836307, 403.482012, 427.474301, 452.893246, 479.823680,
127 508.355480, 538.583870, 570.609734, 604.539954, 640.487770, 678.573155,
128 718.923215, 761.672614, 806.964024, 854.948602, 905.786491, 959.647359,
129 1016.71096, 1077.16774, 1141.21947, 1209.07991, 1280.97554, 1357.14631,
130 1437.84643, 1523.34523, 1613.92805, 1709.89720, 1811.57298, 1919.29472,
131 2033.42192, 2154.33548, 2282.43893, 2418.15982, 2561.95108, 2714.29262,
132 2875.69286, 3046.69045, 3227.85610, 3419.79441, 3623.14597, 3838.58944,
133 4066.84384, 4308.67096, 4564.87787, 4836.31963, 5123.90216, 5428.58524,
134 5751.38572, 6093.38091, 6455.71219, 6839.58882, 7246.29193, 7677.17887,
135 8133.68768, 8617.34192, 9129.75574, 9672.63927, 10247.8043, 10857.1705,
136 11502.7714, 12186.7618, 12911.4244
137 };
138
139
140
141 /* (M)ajor, rotate back 1, rotate back 2
142 (m)inor, rotate back 1, rotate back 2
143 (d)iminished minor, rotate back 1, rotate back 2
144 (f)ifth, rotate back 1, rotate back 2
145 */
146 int chord_table[4][3][3] = {
147 { { 0, 4, 7, }, { -5, 0, 4, }, { -8, -5, 0, }, },
148 { { 0, 3, 7, }, { -5, 0, 3, }, { -9, -5, 0, }, },
149 { { 0, 3, 6, }, { -6, 0, 3, }, { -9, -6, 0, }, },
150 { { 0, 5, 7, }, { -5, 0, 5, }, { -7, -5, 0, }, },
151 };
152
153 /* write the chord type to *chord, returns the root note of the chord */
assign_chord(double * pitchbins,int * chord,int min_guesspitch,int max_guesspitch,int root_pitch)154 int assign_chord(double *pitchbins, int *chord,
155 int min_guesspitch, int max_guesspitch, int root_pitch)
156 {
157
158 int type, subtype;
159 int pitches[19] = {0};
160 int prune_pitches[10] = {0};
161 int i, j, k, n, n2;
162 double val, cutoff, max;
163 int root_flag;
164
165 *chord = -1;
166
167 if (root_pitch - 9 > min_guesspitch)
168 min_guesspitch = root_pitch - 9;
169
170 if (min_guesspitch <= LOWEST_PITCH)
171 min_guesspitch = LOWEST_PITCH + 1;
172
173 if (root_pitch + 9 < max_guesspitch)
174 max_guesspitch = root_pitch + 9;
175
176 if (max_guesspitch >= HIGHEST_PITCH)
177 max_guesspitch = HIGHEST_PITCH - 1;
178
179 /* keep only local maxima */
180 for (i = min_guesspitch, n = 0; i <= max_guesspitch; i++)
181 {
182 val = pitchbins[i];
183 if (val)
184 {
185 if (pitchbins[i-1] < val && pitchbins[i+1] < val)
186 pitches[n++] = i;
187 }
188 }
189
190 if (n < 3)
191 return -1;
192
193 /* find largest peak */
194 max = -1;
195 for (i = 0; i < n; i++)
196 {
197 val = pitchbins[pitches[i]];
198 if (val > max)
199 max = val;
200 }
201
202 /* discard any peaks below cutoff */
203 cutoff = 0.2 * max;
204 for (i = 0, n2 = 0, root_flag = 0; i < n; i++)
205 {
206 val = pitchbins[pitches[i]];
207 if (val >= cutoff)
208 {
209 prune_pitches[n2++] = pitches[i];
210 if (pitches[i] == root_pitch)
211 root_flag = 1;
212 }
213 }
214
215 if (!root_flag || n2 < 3)
216 return -1;
217
218 /* search for a chord, must contain root pitch */
219 for (i = 0; i < n2; i++)
220 {
221 for (subtype = 0; subtype < 3; subtype++)
222 {
223 if (i + subtype >= n2)
224 continue;
225
226 for (type = 0; type < 4; type++)
227 {
228 for (j = 0, n = 0, root_flag = 0; j < 3; j++)
229 {
230 k = i + j;
231
232 if (k >= n2)
233 continue;
234
235 if (prune_pitches[k] == root_pitch)
236 root_flag = 1;
237
238 if (prune_pitches[k] - prune_pitches[i+subtype] ==
239 chord_table[type][subtype][j])
240 n++;
241 }
242 if (root_flag && n == 3)
243 {
244 *chord = 3 * type + subtype;
245 return prune_pitches[i+subtype];
246 }
247 }
248 }
249 }
250
251 return -1;
252 }
253
254
255
256 /* initialize FFT arrays for the frequency analysis */
freq_initialize_fft_arrays(Sample * sp)257 int freq_initialize_fft_arrays(Sample *sp)
258 {
259
260 uint32 i;
261 uint32 length, newlength;
262 unsigned int rate;
263 sample_t *origdata;
264
265 rate = sp->sample_rate;
266 length = sp->data_length >> FRACTION_BITS;
267 origdata = sp->data;
268
269 /* copy the sample to a new float array */
270 floatdata = (float *) safe_malloc(length * sizeof(float));
271 for (i = 0; i < length; i++)
272 floatdata[i] = origdata[i];
273
274 /* length must be a power of 2 */
275 /* set it to smallest power of 2 >= 1.4*rate */
276 /* at least 1.4*rate is required for decent resolution of low notes */
277 newlength = pow(2, ceil(log(1.4*rate) / log(2)));
278 if (length < newlength)
279 {
280 floatdata = safe_realloc(floatdata, newlength * sizeof(float));
281 memset(floatdata + length, 0, (newlength - length) * sizeof(float));
282 }
283 length = newlength;
284
285 /* allocate FFT arrays */
286 /* calculate sin/cos and fft1_bin_to_pitch tables */
287 if (length != oldfftsize)
288 {
289 float f0;
290
291 if (oldfftsize > 0)
292 {
293 free(magdata);
294 free(prunemagdata);
295 free(ip);
296 free(w);
297 free(fft1_bin_to_pitch);
298 }
299 magdata = (float *) safe_malloc(length * sizeof(float));
300 prunemagdata = (float *) safe_malloc(length * sizeof(float));
301 ip = (int *) safe_malloc(2 + sqrt(length) * sizeof(int));
302 *ip = 0;
303 w = (float *) safe_malloc((length >> 1) * sizeof(float));
304 fft1_bin_to_pitch = safe_malloc((length >> 1) * sizeof(float));
305
306 for (i = 1, f0 = (float) rate / length; i < (length >> 1); i++) {
307 fft1_bin_to_pitch[i] = assign_pitch_to_freq(i * f0);
308 }
309 }
310 oldfftsize = length;
311
312 /* zero out arrays that need it */
313 memset(pitchmags, 0, 129 * sizeof(float));
314 memset(pitchbins, 0, 129 * sizeof(double));
315 memset(new_pitchbins, 0, 129 * sizeof(double));
316 memset(prunemagdata, 0, length * sizeof(float));
317
318 return(length);
319 }
320
321
322
323 /* return the frequency of the sample */
324 /* max of 1.4 - 2.0 seconds of audio is analyzed, depending on sample rate */
325 /* samples < 1.4 seconds are padded to max length for higher fft accuracy */
freq_fourier(Sample * sp,int * chord)326 float freq_fourier(Sample *sp, int *chord)
327 {
328
329 uint32 length, length0;
330 int32 maxoffset, minoffset, minoffset1, minoffset2;
331 int32 minbin, maxbin;
332 int32 bin;
333 int32 i, j, n, total;
334 unsigned int rate;
335 int pitch, bestpitch, minpitch, maxpitch, maxpitch2;
336 sample_t *origdata;
337 float f0, mag, maxmag;
338 int16 amp, oldamp, maxamp;
339 int32 maxpos = 0;
340 double sum, weightsum, maxsum;
341 double f0_inv;
342 float freq, newfreq, bestfreq, freq_inc;
343 float minfreq, maxfreq, minfreq2, maxfreq2;
344 float min_guessfreq, max_guessfreq;
345
346 rate = sp->sample_rate;
347 length = length0 = sp->data_length >> FRACTION_BITS;
348 origdata = sp->data;
349
350 length = freq_initialize_fft_arrays(sp);
351
352 /* base frequency of the FFT */
353 f0 = (float) rate / length;
354 f0_inv = 1.0 / f0;
355
356 /* get maximum amplitude */
357 maxamp = -1;
358 for (i = 0; i < length0; i++)
359 {
360 amp = abs(origdata[i]);
361 if (amp >= maxamp)
362 {
363 maxamp = amp;
364 maxpos = i;
365 }
366 }
367
368 /* go out 2 zero crossings in both directions, starting at maxpos */
369 /* find the peaks after the 2nd crossing */
370 minoffset1 = 0;
371 for (n = 0, oldamp = origdata[maxpos], i = maxpos - 1; i >= 0 && n < 2; i--)
372 {
373 amp = origdata[i];
374 if ((oldamp && amp == 0) || (oldamp > 0 && amp < 0) ||
375 (oldamp < 0 && amp > 0))
376 n++;
377 oldamp = amp;
378 }
379 minoffset1 = i;
380 maxamp = labs(origdata[i]);
381 while (i >= 0)
382 {
383 amp = origdata[i];
384 if ((oldamp && amp == 0) || (oldamp > 0 && amp < 0) ||
385 (oldamp < 0 && amp > 0))
386 {
387 break;
388 }
389 oldamp = amp;
390
391 amp = labs(amp);
392 if (amp > maxamp)
393 {
394 maxamp = amp;
395 minoffset1 = i;
396 }
397 i--;
398 }
399
400 minoffset2 = 0;
401 for (n = 0, oldamp = origdata[maxpos], i = maxpos + 1; i < length0 && n < 2; i++)
402 {
403 amp = origdata[i];
404 if ((oldamp && amp == 0) || (oldamp > 0 && amp < 0) ||
405 (oldamp < 0 && amp > 0))
406 n++;
407 oldamp = amp;
408 }
409 minoffset2 = i;
410 maxamp = labs(origdata[i]);
411 while (i < length0)
412 {
413 amp = origdata[i];
414 if ((oldamp && amp == 0) || (oldamp > 0 && amp < 0) ||
415 (oldamp < 0 && amp > 0))
416 {
417 break;
418 }
419 oldamp = amp;
420
421 amp = labs(amp);
422 if (amp > maxamp)
423 {
424 maxamp = amp;
425 minoffset2 = i;
426 }
427 i++;
428 }
429
430 /* upper bound on the detected frequency */
431 /* distance between the two peaks is at most 2 periods */
432 minoffset = (minoffset2 - minoffset1);
433 if (minoffset < 4)
434 minoffset = 4;
435 max_guessfreq = (float) rate / (minoffset * 0.5);
436 if (max_guessfreq >= (rate >> 1)) max_guessfreq = (rate >> 1) - 1;
437
438 /* lower bound on the detected frequency */
439 maxoffset = rate / pitch_freq_lb_table[LOWEST_PITCH] + 0.5;
440 if (maxoffset > (length >> 1))
441 maxoffset = (length >> 1);
442 min_guessfreq = (float) rate / maxoffset;
443
444 /* perform the in place FFT */
445 rdft(length, 1, floatdata, ip, w);
446
447 /* calc the magnitudes */
448 for (i = 2; i < length; i++)
449 {
450 mag = floatdata[i++];
451 mag *= mag;
452 mag += floatdata[i] * floatdata[i];
453 magdata[i >> 1] = sqrt(mag);
454 }
455
456 /* find max mag */
457 maxmag = 0;
458 for (i = 1; i < (length >> 1); i++)
459 {
460 mag = magdata[i];
461
462 pitch = fft1_bin_to_pitch[i];
463 if (pitch && mag > maxmag)
464 maxmag = mag;
465 }
466
467 /* Apply non-linear scaling to the magnitudes
468 * I don't know why this improves the pitch detection, but it does
469 * The best choice of power seems to be between 1.64 - 1.68
470 */
471 for (i = 1; i < (length >> 1); i++)
472 magdata[i] = maxmag * pow(magdata[i] / maxmag, 1.66);
473
474 /* bin the pitches */
475 for (i = 1; i < (length >> 1); i++)
476 {
477 mag = magdata[i];
478
479 pitch = fft1_bin_to_pitch[i];
480 pitchbins[pitch] += mag;
481
482 if (mag > pitchmags[pitch])
483 pitchmags[pitch] = mag;
484 }
485
486 /* zero out lowest pitch, since it contains all lower frequencies too */
487 pitchbins[LOWEST_PITCH] = 0;
488
489 /* find the largest peak */
490 for (i = LOWEST_PITCH + 1, maxsum = -42; i <= HIGHEST_PITCH; i++)
491 {
492 sum = pitchbins[i];
493 if (sum > maxsum)
494 {
495 maxsum = sum;
496 }
497 }
498
499 minpitch = assign_pitch_to_freq(min_guessfreq);
500 if (minpitch > HIGHEST_PITCH) minpitch = HIGHEST_PITCH;
501
502 /* zero out any peak below minpitch */
503 for (i = LOWEST_PITCH + 1; i < minpitch; i++)
504 pitchbins[i] = 0;
505
506 /* remove all pitches below threshold */
507 for (i = minpitch; i <= HIGHEST_PITCH; i++)
508 {
509 if (pitchbins[i] / maxsum < 0.01 && pitchmags[i] / maxmag < 0.01)
510 pitchbins[i] = 0;
511 }
512
513 /* keep local maxima */
514 for (i = LOWEST_PITCH + 1; i < HIGHEST_PITCH; i++)
515 {
516 double temp;
517
518 temp = pitchbins[i];
519
520 /* also keep significant bands to either side */
521 if (temp && pitchbins[i-1] < temp && pitchbins[i+1] < temp)
522 {
523 new_pitchbins[i] = temp;
524
525 temp *= 0.5;
526 if (pitchbins[i-1] >= temp)
527 new_pitchbins[i-1] = pitchbins[i-1];
528 if (pitchbins[i+1] >= temp)
529 new_pitchbins[i+1] = pitchbins[i-1];
530 }
531 }
532 memcpy(pitchbins, new_pitchbins, 129 * sizeof(double));
533
534 /* find lowest and highest pitches */
535 minpitch = LOWEST_PITCH;
536 while (minpitch < HIGHEST_PITCH && !pitchbins[minpitch])
537 minpitch++;
538 maxpitch = HIGHEST_PITCH;
539 while (maxpitch > LOWEST_PITCH && !pitchbins[maxpitch])
540 maxpitch--;
541
542 /* uh oh, no pitches left...
543 * best guess is middle C
544 * return 260 Hz, since exactly 260 Hz is never returned except on error
545 * this should only occur on blank/silent samples
546 */
547 if (maxpitch < minpitch)
548 {
549 free(floatdata);
550 return 260.0;
551 }
552
553 /* pitch assignment bounds based on zero crossings and pitches kept */
554 if (pitch_freq_lb_table[minpitch] > min_guessfreq)
555 min_guessfreq = pitch_freq_lb_table[minpitch];
556 if (pitch_freq_ub_table[maxpitch] < max_guessfreq)
557 max_guessfreq = pitch_freq_ub_table[maxpitch];
558
559 minfreq = pitch_freq_lb_table[minpitch];
560 if (minfreq >= (rate >> 1)) minfreq = (rate >> 1) - 1;
561
562 maxfreq = pitch_freq_ub_table[maxpitch];
563 if (maxfreq >= (rate >> 1)) maxfreq = (rate >> 1) - 1;
564
565 minbin = minfreq / f0;
566 if (!minbin)
567 minbin = 1;
568 maxbin = ceil(maxfreq / f0);
569 if (maxbin >= (length >> 1))
570 maxbin = (length >> 1) - 1;
571
572 /* filter out all "noise" from magnitude array */
573 for (i = minbin, n = 0; i <= maxbin; i++)
574 {
575 pitch = fft1_bin_to_pitch[i];
576 if (pitchbins[pitch])
577 {
578 prunemagdata[i] = magdata[i];
579 n++;
580 }
581 }
582
583 /* whoa!, there aren't any strong peaks at all !!! bomb early
584 * best guess is middle C
585 * return 260 Hz, since exactly 260 Hz is never returned except on error
586 * this should only occur on blank/silent samples
587 */
588 if (!n)
589 {
590 free(floatdata);
591 return 260.0;
592 }
593
594 memset(new_pitchbins, 0, 129 * sizeof(double));
595
596 maxsum = -1;
597 minpitch = assign_pitch_to_freq(min_guessfreq);
598 maxpitch = assign_pitch_to_freq(max_guessfreq);
599 maxpitch2 = assign_pitch_to_freq(max_guessfreq) + 9;
600 if (maxpitch2 > HIGHEST_PITCH) maxpitch2 = HIGHEST_PITCH;
601
602 /* initial guess is first local maximum */
603 bestfreq = pitch_freq_table[minpitch];
604 if (minpitch < HIGHEST_PITCH &&
605 pitchbins[minpitch+1] > pitchbins[minpitch])
606 bestfreq = pitch_freq_table[minpitch+1];
607
608 /* find best fundamental */
609 for (i = minpitch; i <= maxpitch2; i++)
610 {
611 if (!pitchbins[i])
612 continue;
613
614 minfreq2 = pitch_freq_lb_table[i];
615 maxfreq2 = pitch_freq_ub_table[i];
616 freq_inc = (maxfreq2 - minfreq2) * 0.1;
617 if (minfreq2 >= (rate >> 1)) minfreq2 = (rate >> 1) - 1;
618 if (maxfreq2 >= (rate >> 1)) maxfreq2 = (rate >> 1) - 1;
619
620 /* look for harmonics */
621 for (freq = minfreq2; freq <= maxfreq2; freq += freq_inc)
622 {
623 n = total = 0;
624 sum = weightsum = 0;
625
626 for (j = 1; j <= 32 && (newfreq = j*freq) <= maxfreq; j++)
627 {
628 pitch = assign_pitch_to_freq(newfreq);
629
630 if (pitchbins[pitch])
631 {
632 sum += pitchbins[pitch];
633 n++;
634 total = j;
635 }
636 }
637
638 /* only pitches with good harmonics are assignment candidates */
639 if (n > 1)
640 {
641 double ratio;
642
643 ratio = (double) n / total;
644 if (ratio >= 0.333333)
645 {
646 weightsum = ratio * sum;
647 pitch = assign_pitch_to_freq(freq);
648
649 /* use only these pitches for chord detection */
650 if (pitch <= HIGHEST_PITCH && pitchbins[pitch])
651 new_pitchbins[pitch] = weightsum;
652
653 if (pitch > maxpitch)
654 continue;
655
656 if (n < 2 || weightsum > maxsum)
657 {
658 maxsum = weightsum;
659 bestfreq = freq;
660 }
661 }
662 }
663 }
664 }
665
666 bestpitch = assign_pitch_to_freq(bestfreq);
667
668 /* assign chords */
669 if ((pitch = assign_chord(new_pitchbins, chord,
670 bestpitch - 9, maxpitch2, bestpitch)) >= 0)
671 bestpitch = pitch;
672
673 bestfreq = pitch_freq_table[bestpitch];
674
675 /* tune based on the fundamental and harmonics up to +5 octaves */
676 sum = weightsum = 0;
677 for (i = 1; i <= 32 && (freq = i*bestfreq) <= maxfreq; i++)
678 {
679 double tune;
680
681 minfreq2 = pitch_freq_lb_table[bestpitch];
682 maxfreq2 = pitch_freq_ub_table[bestpitch];
683
684 minbin = minfreq2 * f0_inv;
685 if (!minbin) minbin = 1;
686 maxbin = ceil(maxfreq2 * f0_inv);
687 if (maxbin >= (length>>1))
688 maxbin = (length>>1) - 1;
689
690 for (bin = minbin; bin <= maxbin; bin++)
691 {
692 tune = -36.37631656 + 17.31234049 * log(bin*f0) - bestpitch;
693 sum += magdata[bin];
694 weightsum += magdata[bin] * tune;
695 }
696 }
697
698 bestfreq = 13.75 * exp(((bestpitch + weightsum / sum) - 9) /
699 12 * log(2));
700
701 /* Since we are using exactly 260 Hz as an error code, fudge the freq
702 * on the extremely unlikely chance that the detected pitch is exactly
703 * 260 Hz.
704 */
705 if (bestfreq == 260.0)
706 bestfreq += 1E-5;
707
708 free(floatdata);
709
710 return bestfreq;
711 }
712
713
714
assign_pitch_to_freq(float freq)715 int assign_pitch_to_freq(float freq)
716 {
717
718 int pitch;
719
720 /* round to nearest integer using: ceil(fraction - 0.5) */
721 /* -0.5 is already added into the first constant below */
722 pitch = ceil(-36.87631656f + 17.31234049f * log(freq));
723
724 /* min and max pitches */
725 if (pitch < LOWEST_PITCH) pitch = LOWEST_PITCH;
726 else if (pitch > HIGHEST_PITCH) pitch = HIGHEST_PITCH;
727
728 return pitch;
729 }
730