1 /*
2 * lingot, a musical instrument tuner.
3 *
4 * Copyright (C) 2004-2018 Iban Cereijo.
5 * Copyright (C) 2004-2008 Jairo Chapela.
6
7 *
8 * This file is part of lingot.
9 *
10 * lingot is free software; you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation; either version 2 of the License, or
13 * (at your option) any later version.
14 *
15 * lingot is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with lingot; if not, write to the Free Software
22 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 */
24
25 #include <stdio.h>
26 #include <math.h>
27 #include <sys/soundcard.h>
28 #include <string.h>
29 #include <errno.h>
30 #include <sys/time.h>
31 #include <time.h>
32 #include <stdlib.h>
33 #include "lingot-fft.h"
34 #include "lingot-signal.h"
35 #include "lingot-core.h"
36 #include "lingot-config.h"
37 #include "lingot-i18n.h"
38 #include "lingot-msg.h"
39
40 int
41 lingot_core_read_callback(FLT* read_buffer, unsigned int samples_read, void *arg);
42
43 void lingot_core_run_computation_thread(LingotCore* core);
44
45 unsigned int decimation_input_index = 0;
46
lingot_core_new(LingotCore * core,LingotConfig * conf)47 void lingot_core_new(LingotCore* core, LingotConfig* conf) {
48
49 char buff[1000];
50
51 lingot_config_copy(&core->conf, conf);
52 core->running = 0;
53 core->spd_fft = NULL;
54 core->noise_level = NULL;
55 core->SPL = NULL;
56 core->flt_read_buffer = NULL;
57 core->temporal_buffer = NULL;
58 core->windowed_temporal_buffer = NULL;
59 core->windowed_fft_buffer = NULL;
60 core->hamming_window_temporal = NULL;
61 core->hamming_window_fft = NULL;
62
63 #ifdef DRAW_MARKERS
64 core->markers_size = 0;
65 core->markers_size2 = 0;
66 #endif
67
68 unsigned int requested_sample_rate = core->conf.sample_rate;
69
70 if (core->conf.sample_rate <= 0) {
71 core->conf.sample_rate = 0;
72 }
73
74 lingot_audio_new(&core->audio, core->conf.audio_system,
75 core->conf.audio_dev[core->conf.audio_system], core->conf.sample_rate,
76 (LingotAudioProcessCallback) lingot_core_read_callback, core);
77
78 if (core->audio.audio_system != -1) {
79
80 if (requested_sample_rate != core->audio.real_sample_rate) {
81 core->conf.sample_rate = core->audio.real_sample_rate;
82 lingot_config_update_internal_params(conf);
83 // sprintf(buff,
84 // _("The requested sample rate is not available, the real sample rate has been set to %d Hz"),
85 // core->audio.real_sample_rate);
86 // lingot_msg_add_warning(buff);
87 }
88
89 if (core->conf.temporal_buffer_size < core->conf.fft_size) {
90 core->conf.temporal_window = ((double) core->conf.fft_size
91 * core->conf.oversampling) / core->conf.sample_rate;
92 core->conf.temporal_buffer_size = core->conf.fft_size;
93 lingot_config_update_internal_params(conf);
94 snprintf(buff, sizeof(buff),
95 _(
96 "The temporal buffer is smaller than FFT size. It has been increased to %0.3f seconds"),
97 core->conf.temporal_window);
98 lingot_msg_add_warning(buff);
99 }
100
101 // Since the SPD is symmetrical, we only store the 1st half.
102 const unsigned int spd_size = (core->conf.fft_size / 2);
103
104 core->spd_fft = malloc(spd_size * sizeof(FLT));
105 core->noise_level = malloc(spd_size * sizeof(FLT));
106 core->SPL = malloc(spd_size * sizeof(FLT));
107
108 memset(core->spd_fft, 0, spd_size * sizeof(FLT));
109 memset(core->noise_level, 0, spd_size * sizeof(FLT));
110 memset(core->SPL, 0, spd_size * sizeof(FLT));
111
112 // audio source read in floating point format.
113 core->flt_read_buffer = malloc(
114 core->audio.read_buffer_size_samples * sizeof(FLT));
115 memset(core->flt_read_buffer, 0,
116 core->audio.read_buffer_size_samples * sizeof(FLT));
117
118 // stored samples.
119 core->temporal_buffer = malloc(
120 (core->conf.temporal_buffer_size) * sizeof(FLT));
121 memset(core->temporal_buffer, 0,
122 core->conf.temporal_buffer_size * sizeof(FLT));
123
124 core->hamming_window_temporal = NULL;
125 core->hamming_window_fft = NULL;
126
127 if (core->conf.window_type != NONE) {
128 core->hamming_window_temporal = malloc(
129 (core->conf.temporal_buffer_size) * sizeof(FLT));
130 core->hamming_window_fft = malloc(
131 (core->conf.fft_size) * sizeof(FLT));
132
133 lingot_signal_window(core->conf.temporal_buffer_size,
134 core->hamming_window_temporal, core->conf.window_type);
135 lingot_signal_window(core->conf.fft_size, core->hamming_window_fft,
136 core->conf.window_type);
137 }
138
139 core->windowed_temporal_buffer = malloc(
140 (core->conf.temporal_buffer_size) * sizeof(FLT));
141 memset(core->windowed_temporal_buffer, 0,
142 core->conf.temporal_buffer_size * sizeof(FLT));
143 core->windowed_fft_buffer = malloc(
144 (core->conf.fft_size) * sizeof(FLT));
145 memset(core->windowed_fft_buffer, 0,
146 core->conf.fft_size * sizeof(FLT));
147
148 lingot_fft_plan_create(&core->fftplan, core->windowed_fft_buffer,
149 core->conf.fft_size);
150
151 /*
152 * 8 order Chebyshev filters, with wc=0.9/i (normalised respect to
153 * Pi). We take 0.9 instead of 1 to leave a 10% of safety margin,
154 * in order to avoid aliased frequencies near to w=Pi, due to non
155 * ideality of the filter.
156 *
157 * The cut frequencies wc=Pi/i, with i=1..20, correspond with the
158 * oversampling factor, avoiding aliasing at decimation.
159 *
160 * Why Chebyshev filters?, for a given order, those filters yield
161 * abrupt falls than other ones as Butterworth, making the most of
162 * the order. Although Chebyshev filters affect more to the phase,
163 * it doesn't matter due to the analysis is made on the signal
164 * power distribution (only magnitude).
165 */
166 lingot_filter_cheby_design(&core->antialiasing_filter, 8, 0.5,
167 0.9 / core->conf.oversampling);
168
169 pthread_mutex_init(&core->temporal_buffer_mutex, NULL);
170
171 // ------------------------------------------------------------
172
173 core->running = 1;
174 }
175
176 core->freq = 0.0;
177 }
178
179 // -----------------------------------------------------------------------
180
181 /* Deallocate resources */
lingot_core_destroy(LingotCore * core)182 void lingot_core_destroy(LingotCore* core) {
183
184 if (core->audio.audio_system != -1) {
185 lingot_fft_plan_destroy(&core->fftplan);
186 lingot_audio_destroy(&core->audio);
187
188 free(core->spd_fft);
189 free(core->noise_level);
190 free(core->SPL);
191 free(core->flt_read_buffer);
192 free(core->temporal_buffer);
193
194 free(core->hamming_window_temporal);
195 free(core->windowed_temporal_buffer);
196 free(core->hamming_window_fft);
197 free(core->windowed_fft_buffer);
198
199 lingot_filter_destroy(&core->antialiasing_filter);
200
201 pthread_mutex_destroy(&core->temporal_buffer_mutex);
202 }
203 }
204
205 // -----------------------------------------------------------------------
206
207 // reads a new piece of signal from audio source, applies filtering and
208 // decimation and appends it to the buffer
lingot_core_read_callback(FLT * read_buffer,unsigned int samples_read,void * arg)209 int lingot_core_read_callback(FLT* read_buffer, unsigned int samples_read, void *arg) {
210
211 unsigned int decimation_output_index; // loop variables.
212 unsigned int decimation_output_len;
213 FLT* decimation_in;
214 FLT* decimation_out;
215 LingotCore* core = (LingotCore*) arg;
216 const LingotConfig* const conf = &core->conf;
217
218 memcpy(core->flt_read_buffer, read_buffer, samples_read * sizeof(FLT));
219 // double omega = 2.0 * M_PI * 100.0;
220 // double T = 1.0 / conf->sample_rate;
221 // static double t = 0.0;
222 // for (i = 0; i < read_buffer_size; i++) {
223 // core->flt_read_buffer[i] = 1e4
224 // * (cos(omega * t) + (0.5 * rand()) / RAND_MAX);
225 // t += T;
226 // }
227
228 // if (conf->gain_nu != 1.0) {
229 // for (i = 0; i < read_buffer_size; i++)
230 // core->flt_read_buffer[i] *= conf->gain_nu;
231 // }
232
233 //
234 // just read:
235 //
236 // ----------------------------
237 // |bxxxbxxxbxxxbxxxbxxxbxxxbxxx|
238 // ----------------------------
239 //
240 // <----------------------------> samples_read
241 //
242
243 decimation_output_len = 1
244 + (samples_read - (decimation_input_index + 1))
245 / conf->oversampling;
246
247 //#define DUMP
248
249 #ifdef DUMP
250 static FILE* fid0 = 0x0;
251 if (fid0 == 0x0) {
252 fid0 = fopen("/tmp/dump_pre_filter.txt", "w");
253 }
254
255 for (i = 0; i < samples_read; i++) {
256 fprintf(fid0, "%f ", core->flt_read_buffer[i]);
257 }
258 #endif
259
260 pthread_mutex_lock(&core->temporal_buffer_mutex);
261
262 /* we shift the temporal window to leave a hollow where place the new piece
263 of data read. The buffer is actually a queue. */
264 if (conf->temporal_buffer_size > decimation_output_len) {
265 memmove(core->temporal_buffer,
266 &core->temporal_buffer[decimation_output_len],
267 (conf->temporal_buffer_size - decimation_output_len)
268 * sizeof(FLT));
269 }
270
271 //
272 // previous buffer situation:
273 //
274 // ------------------------------------------
275 // | xxxxxxxxxxxxxxxxxxxxxx | yyyyy | aaaaaaa |
276 // ------------------------------------------
277 // <------> samples_read/oversampling
278 // <---------------> fft_size
279 // <----------------------------------------> temporal_buffer_size
280 //
281 // new situation:
282 //
283 // ------------------------------------------
284 // | xxxxxxxxxxxxxxxxyyyyaa | aaaaa | |
285 // ------------------------------------------
286 //
287
288 // decimation with low-pass filtering
289
290 /* we decimate the signal and append it at the end of the buffer. */
291 if (conf->oversampling > 1) {
292
293 decimation_in = core->flt_read_buffer;
294 decimation_out = &core->temporal_buffer[conf->temporal_buffer_size
295 - decimation_output_len];
296
297 // low pass filter to avoid aliasing.
298 lingot_filter_filter(&core->antialiasing_filter, samples_read,
299 decimation_in, decimation_in);
300
301 // downsampling.
302 for (decimation_output_index = 0; decimation_input_index < samples_read;
303 decimation_output_index++, decimation_input_index +=
304 conf->oversampling) {
305 decimation_out[decimation_output_index] =
306 decimation_in[decimation_input_index];
307 }
308 decimation_input_index -= samples_read;
309 } else {
310 memcpy(
311 &core->temporal_buffer[conf->temporal_buffer_size
312 - decimation_output_len], core->flt_read_buffer,
313 decimation_output_len * sizeof(FLT));
314 }
315 //
316 // ------------------------------------------
317 // | xxxxxxxxxxxxxxxxyyyyaa | aaaaa | bbbbbbb |
318 // ------------------------------------------
319 //
320
321 pthread_mutex_unlock(&core->temporal_buffer_mutex);
322
323 #ifdef DUMP
324 static FILE* fid1 = 0x0;
325 static FILE* fid2 = 0x0;
326
327 if (fid1 == 0x0) {
328 fid1 = fopen("/tmp/dump_post_filter.txt", "w");
329 fid2 = fopen("/tmp/dump_post.txt", "w");
330 }
331
332 for (i = 0; i < samples_read; i++) {
333 fprintf(fid1, "%f ", core->flt_read_buffer[i]);
334 }
335 decimation_out = &core->temporal_buffer[conf->temporal_buffer_size
336 - decimation_output_len];
337 for (i = 0; i < decimation_output_len; i++) {
338 fprintf(fid2, "%f ", decimation_out[i]);
339 }
340 #endif
341
342 return 0;
343 }
344
lingot_core_frequencies_related(FLT freq1,FLT freq2,FLT minFrequency,FLT * mulFreq1ToFreq,FLT * mulFreq2ToFreq)345 int lingot_core_frequencies_related(FLT freq1, FLT freq2, FLT minFrequency,
346 FLT* mulFreq1ToFreq, FLT* mulFreq2ToFreq) {
347
348 int result = 0;
349 static const FLT tol = 5e-2; // TODO: tune
350 static const int maxDivisor = 4;
351
352 if ((freq1 != 0.0) && (freq2 != 0.0)) {
353
354 FLT smallFreq = freq1;
355 FLT bigFreq = freq2;
356
357 if (bigFreq < smallFreq) {
358 smallFreq = freq2;
359 bigFreq = freq1;
360 }
361
362 int divisor;
363 FLT frac;
364 FLT error = -1.0;
365 for (divisor = 1; divisor <= maxDivisor; divisor++) {
366 if (minFrequency * divisor > smallFreq) {
367 break;
368 }
369
370 frac = bigFreq * divisor / smallFreq;
371 error = fabs(frac - round(frac));
372 if (error < tol) {
373 if (smallFreq == freq1) {
374 *mulFreq1ToFreq = 1.0 / divisor;
375 *mulFreq2ToFreq = 1.0 / round(frac);
376 } else {
377 *mulFreq1ToFreq = 1.0 / round(frac);
378 *mulFreq2ToFreq = 1.0 / divisor;
379 }
380 result = 1;
381 break;
382 }
383 }
384 } else {
385 *mulFreq1ToFreq = 0.0;
386 *mulFreq2ToFreq = 0.0;
387 }
388
389 // printf("relation %f, %f = %i\n", freq1, freq2, result);
390
391 return result;
392 }
393
lingot_core_frequency_locker(FLT freq,FLT minFrequency)394 static FLT lingot_core_frequency_locker(FLT freq, FLT minFrequency) {
395
396 static int locked = 0;
397 static FLT current_frequency = -1.0;
398 static int hits_counter = 0;
399 static int rehits_counter = 0;
400 static int rehits_up_counter = 0;
401 static const int nhits_to_lock = 4;
402 static const int nhits_to_unlock = 5;
403 static const int nhits_to_relock = 6;
404 static const int nhits_to_relock_up = 8;
405 FLT multiplier = 0.0;
406 FLT multiplier2 = 0.0;
407 static FLT old_multiplier = 0.0;
408 static FLT old_multiplier2 = 0.0;
409 int fail = 0;
410 FLT result = 0.0;
411
412 #ifdef DRAW_MARKERS
413 printf("f = %f\n", freq);
414 #endif
415 int consistent_with_current_frequency = 0;
416 consistent_with_current_frequency = lingot_core_frequencies_related(freq,
417 current_frequency, minFrequency, &multiplier, &multiplier2);
418
419 if (!locked) {
420
421 if ((freq > 0.0) && (current_frequency == 0.0)) {
422 consistent_with_current_frequency = 1;
423 multiplier = 1.0;
424 multiplier2 = 1.0;
425 }
426
427 // printf("filtering frequency %f, current %f\n", freq, current_frequency);
428
429 if (consistent_with_current_frequency && (multiplier == 1.0)
430 && (multiplier2 == 1.0)) {
431 current_frequency = freq * multiplier;
432
433 if (++hits_counter >= nhits_to_lock) {
434 locked = 1;
435 #ifdef DRAW_MARKERS
436 printf("locked to frequency %f\n", current_frequency);
437 #endif
438 hits_counter = 0;
439 }
440 } else {
441 hits_counter = 0;
442 current_frequency = 0.0;
443 }
444
445 // result = freq;
446 } else {
447 // printf("c = %i, f = %f, cf = %f, multiplier = %f, multiplier2 = %f\n",
448 // consistent_with_current_frequency, freq, current_frequency,
449 // multiplier, multiplier2);
450
451 if (consistent_with_current_frequency) {
452 if (fabs(multiplier2 - 1.0) < 1e-5) {
453 result = freq * multiplier;
454 current_frequency = result;
455 rehits_counter = 0;
456
457 if (fabs(multiplier - 1.0) > 1e-5) {
458 if (fabs(multiplier - old_multiplier) < 1e-5) {
459 #ifdef DRAW_MARKERS
460 printf("SEIN!!!! %f!\n", multiplier);
461 #endif
462 if (++rehits_up_counter >= nhits_to_relock_up) {
463 result = freq;
464 current_frequency = result;
465 #ifdef DRAW_MARKERS
466 printf("relock UP!! to %f\n\n\n", freq);
467 #endif
468 rehits_up_counter = 0;
469 fail = 0;
470 }
471 } else {
472 rehits_up_counter = 0;
473 }
474 } else {
475 rehits_up_counter = 0;
476 }
477 } else {
478 rehits_up_counter = 0;
479 #ifdef DRAW_MARKERS
480 printf("%f!\n", multiplier2);
481 #endif
482 if (fabs(multiplier2 - 0.5) < 1e-5) {
483 hits_counter--;
484 }
485 fail = 1;
486 if (freq * multiplier < minFrequency) {
487 #ifdef DRAW_MARKERS
488 printf("warning freq * mul = %f < min\n",
489 freq * multiplier);
490 #endif
491 } else {
492 // result = freq * multiplier;
493 // printf("hop detected!\n");
494 // current_frequency = result;
495
496 #ifdef DRAW_MARKERS
497 printf("(%f == %f)?\n", multiplier2, old_multiplier2);
498 #endif
499 if (fabs(multiplier2 - old_multiplier2) < 1e-5) {
500 #ifdef DRAW_MARKERS
501 printf("match for relock, %f == %f\n", multiplier2,
502 old_multiplier2);
503 #endif
504 if (++rehits_counter >= nhits_to_relock) {
505 result = freq * multiplier;
506 current_frequency = result;
507 #ifdef DRAW_MARKERS
508 printf("relock!! to %f\n", freq);
509 #endif
510 rehits_counter = 0;
511 fail = 0;
512 }
513 }
514 }
515 }
516
517 } else {
518 fail = 1;
519 }
520
521 if (fail) {
522 result = current_frequency;
523 hits_counter++;
524 if (hits_counter >= nhits_to_unlock) {
525 current_frequency = 0.0;
526 locked = 0;
527 hits_counter = 0;
528 #ifdef DRAW_MARKERS
529 printf("unlocked\n");
530 #endif
531 result = 0.0;
532 }
533 } else {
534 hits_counter = 0;
535 }
536 }
537
538 old_multiplier = multiplier;
539 old_multiplier2 = multiplier2;
540
541 // if (result != 0.0)
542 // printf("result = %f\n", result);
543 return result;
544 }
545
lingot_core_compute_fundamental_fequency(LingotCore * core)546 void lingot_core_compute_fundamental_fequency(LingotCore* core) {
547
548 register unsigned int i, k; // loop variables.
549 const LingotConfig* const conf = &core->conf;
550 const FLT index2f = ((FLT) conf->sample_rate)
551 / (conf->oversampling * conf->fft_size); // FFT resolution in Hz.
552 // const FLT index2w = 2.0 * M_PI / conf->fft_size; // FFT resolution in rads.
553 // const FLT f2w = 2 * M_PI * conf->oversampling / conf->sample_rate;
554 // const FLT w2f = 1.0 / f2w;
555
556 // ----------------- TRANSFORMATION TO FREQUENCY DOMAIN ----------------
557
558 pthread_mutex_lock(&core->temporal_buffer_mutex);
559
560 // windowing
561 if (conf->window_type != NONE) {
562 for (i = 0; i < conf->fft_size; i++) {
563 core->windowed_fft_buffer[i] =
564 core->temporal_buffer[conf->temporal_buffer_size
565 - conf->fft_size + i] * core->hamming_window_fft[i];
566 }
567 } else {
568 memmove(core->windowed_fft_buffer,
569 &core->temporal_buffer[conf->temporal_buffer_size
570 - conf->fft_size], conf->fft_size * sizeof(FLT));
571 }
572
573 const unsigned int spd_size = (conf->fft_size / 2);
574
575 // FFT
576 lingot_fft_compute_dft_and_spd(&core->fftplan, core->spd_fft, spd_size);
577
578 static const FLT minSPL = -200;
579 for (i = 0; i < spd_size; i++) {
580 core->SPL[i] = 10.0 * log10(core->spd_fft[i]);
581 if (core->SPL[i] < minSPL) {
582 core->SPL[i] = minSPL;
583 }
584 }
585
586 FLT noise_filter_width = 150.0; // hz
587 unsigned int noise_filter_width_samples = ceil(
588 noise_filter_width * conf->fft_size * conf->oversampling
589 / conf->sample_rate);
590
591 lingot_signal_compute_noise_level(core->SPL, spd_size,
592 noise_filter_width_samples, core->noise_level);
593 for (i = 0; i < spd_size; i++) {
594 core->SPL[i] -= core->noise_level[i];
595 }
596
597 unsigned int lowest_index = (unsigned int) ceil(
598 conf->internal_min_frequency
599 * (1.0 * conf->oversampling / conf->sample_rate)
600 * conf->fft_size);
601 unsigned int highest_index = (unsigned int) ceil(0.95 * spd_size);
602
603 short divisor = 1;
604 FLT f0 = lingot_signal_estimate_fundamental_frequency(core->SPL,
605 0.5 * core->freq, core->fftplan.fft_out, spd_size,
606 conf->peak_number, lowest_index, highest_index,
607 conf->peak_half_width, index2f, conf->min_SNR,
608 conf->min_overall_SNR, conf->internal_min_frequency, core,
609 &divisor);
610
611 FLT w;
612 FLT w0 =
613 (f0 == 0.0) ?
614 0.0 :
615 2 * M_PI * f0 * conf->oversampling / conf->sample_rate;
616 w = w0;
617 // int Mi = floor(w / index2w);
618
619 if (w != 0.0) {
620 // windowing
621 if (conf->window_type != NONE) {
622 for (i = 0; i < conf->temporal_buffer_size; i++) {
623 core->windowed_temporal_buffer[i] = core->temporal_buffer[i]
624 * core->hamming_window_temporal[i];
625 }
626 } else {
627 memmove(core->windowed_temporal_buffer, core->temporal_buffer,
628 conf->temporal_buffer_size * sizeof(FLT));
629 }
630 }
631
632 pthread_mutex_unlock(&core->temporal_buffer_mutex); // we don't need the read buffer anymore
633
634 if (w != 0.0) {
635
636 // Maximum finding by Newton-Raphson
637 // -----------------------------------
638
639 FLT wk = -1.0e5;
640 FLT wkm1 = w;
641 // first iterator set to the current approximation.
642 FLT d0_SPD = 0.0;
643 FLT d1_SPD = 0.0;
644 FLT d2_SPD = 0.0;
645 FLT d0_SPD_old = 0.0;
646
647 // printf("NR iter: %f ", w * w2f);
648
649 for (k = 0; (k < conf->max_nr_iter) && (fabs(wk - wkm1) > 1.0e-4);
650 k++) {
651 wk = wkm1;
652
653 d0_SPD_old = d0_SPD;
654 lingot_fft_spd_diffs_eval(core->windowed_fft_buffer, // TODO: iterate over this buffer?
655 conf->fft_size, wk, &d0_SPD, &d1_SPD, &d2_SPD);
656
657 wkm1 = wk - d1_SPD / d2_SPD;
658 // printf(" -> (%f,%g,%g,%g)", wkm1 * w2f, d0_SPD, d1_SPD, d2_SPD);
659
660 if (d0_SPD < d0_SPD_old) {
661 // printf("!!!", d0_SPD, d0_SPD_old);
662 wkm1 = 0.0;
663 break;
664 }
665
666 }
667 // printf("\n");
668
669 if (wkm1 > 0.0) {
670 w = wkm1; // frequency in rads.
671 wk = -1.0e5;
672 d0_SPD = 0.0;
673 // printf("NR2 iter: %f ", w * w2f);
674
675 for (k = 0;
676 (k <= 1)
677 || ((k < conf->max_nr_iter)
678 && (fabs(wk - wkm1) > 1.0e-4)); k++) {
679 wk = wkm1;
680
681 // ! we use the WHOLE temporal window for bigger precision.
682 d0_SPD_old = d0_SPD;
683 lingot_fft_spd_diffs_eval(core->windowed_temporal_buffer,
684 conf->temporal_buffer_size, wk, &d0_SPD, &d1_SPD,
685 &d2_SPD);
686
687 wkm1 = wk - d1_SPD / d2_SPD;
688 // printf(" -> (%f,%g,%g,%g)", wkm1 * w2f, d0_SPD, d1_SPD, d2_SPD);
689
690 if (d0_SPD < d0_SPD_old) {
691 // printf("!!!");
692 wkm1 = 0.0;
693 break;
694 }
695
696 }
697 // printf("\n");
698
699 if (wkm1 > 0.0) {
700 w = wkm1; // frequency in rads.
701 }
702 }
703 }
704
705 FLT freq =
706 (w == 0.0) ?
707 0.0 :
708 w * conf->sample_rate
709 / (divisor * 2.0 * M_PI * conf->oversampling); // analog frequency in Hz.
710 // core->freq = freq;
711 core->freq = lingot_core_frequency_locker(freq,
712 core->conf.internal_min_frequency);
713 // printf("-> %f\n", core->freq);
714 }
715
716 /* start running the core in another thread */
lingot_core_start(LingotCore * core)717 void lingot_core_start(LingotCore* core) {
718
719 int audio_status = 0;
720 decimation_input_index = 0;
721
722 if (core->audio.audio_system != -1) {
723 audio_status = lingot_audio_start(&core->audio);
724
725 if (audio_status == 0) {
726 pthread_mutex_init(&core->thread_computation_mutex, NULL);
727 pthread_cond_init(&core->thread_computation_cond, NULL);
728
729 pthread_attr_init(&core->thread_computation_attr);
730 pthread_create(&core->thread_computation,
731 &core->thread_computation_attr,
732 (void* (*)(void*)) lingot_core_run_computation_thread,
733 core);
734 core->running = 1;
735 } else {
736 core->running = 0;
737 lingot_audio_destroy(&core->audio);
738 }
739
740 }
741 }
742
743 /* stop running the core */
lingot_core_stop(LingotCore * core)744 void lingot_core_stop(LingotCore* core) {
745 void* thread_result;
746
747 int result;
748 struct timeval tout, tout_abs;
749 struct timespec tout_tspec;
750
751 gettimeofday(&tout_abs, NULL);
752 tout.tv_sec = 0;
753 tout.tv_usec = 300000;
754
755 if (core->running == 1) {
756 core->running = 0;
757
758 timeradd(&tout, &tout_abs, &tout_abs);
759 tout_tspec.tv_sec = tout_abs.tv_sec;
760 tout_tspec.tv_nsec = 1000 * tout_abs.tv_usec;
761
762 // watchdog timer
763 pthread_mutex_lock(&core->thread_computation_mutex);
764 result = pthread_cond_timedwait(&core->thread_computation_cond,
765 &core->thread_computation_mutex, &tout_tspec);
766 pthread_mutex_unlock(&core->thread_computation_mutex);
767
768 if (result == ETIMEDOUT) {
769 fprintf(stderr, "warning: cancelling computation thread\n");
770 pthread_cancel(core->thread_computation);
771 } else {
772 pthread_join(core->thread_computation, &thread_result);
773 }
774 pthread_attr_destroy(&core->thread_computation_attr);
775 pthread_mutex_destroy(&core->thread_computation_mutex);
776 pthread_cond_destroy(&core->thread_computation_cond);
777
778 int spd_size = core->conf.fft_size / 2;
779 memset(core->SPL, 0, spd_size * sizeof(FLT));
780 core->freq = 0.0;
781 }
782
783 if (core->audio.audio_system != -1) {
784 lingot_audio_stop(&core->audio);
785 }
786 }
787
788 /* run the core */
lingot_core_run_computation_thread(LingotCore * core)789 void lingot_core_run_computation_thread(LingotCore* core) {
790 struct timeval tout, tout_abs;
791 struct timespec tout_tspec;
792
793 gettimeofday(&tout_abs, NULL);
794 tout.tv_sec = 0;
795 tout.tv_usec = 1e6 / core->conf.calculation_rate;
796
797 while (core->running) {
798 lingot_core_compute_fundamental_fequency(core);
799 timeradd(&tout, &tout_abs, &tout_abs);
800 tout_tspec.tv_sec = tout_abs.tv_sec;
801 tout_tspec.tv_nsec = 1000 * tout_abs.tv_usec;
802 pthread_mutex_lock(&core->thread_computation_mutex);
803 pthread_cond_timedwait(&core->thread_computation_cond,
804 &core->thread_computation_mutex, &tout_tspec);
805 pthread_mutex_unlock(&core->thread_computation_mutex);
806
807 if (core->audio.audio_system != -1) {
808 const unsigned int spd_size = core->conf.fft_size / 2;
809 if (core->audio.interrupted) {
810 memset(core->SPL, 0, spd_size * sizeof(FLT));
811 core->freq = 0.0;
812 core->running = 0;
813 }
814 }
815 }
816
817 pthread_mutex_lock(&core->thread_computation_mutex);
818 pthread_cond_broadcast(&core->thread_computation_cond);
819 pthread_mutex_unlock(&core->thread_computation_mutex);
820 }
821