/*
* lingot, a musical instrument tuner.
*
* Copyright (C) 2004-2011 Ibán Cereijo Graña, Jairo Chapela Martínez.
*
* This file is part of lingot.
*
* lingot is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* lingot is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with lingot; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
#include <stdio.h>
#include <math.h>
#include <sys/soundcard.h>
#include <string.h>
#include <errno.h>
#include <sys/time.h>
#include <time.h>
#include <stdlib.h>
#include "lingot-complex.h"
#include "lingot-fft.h"
#include "lingot-signal.h"
#include "lingot-core.h"
#include "lingot-config.h"
#include "lingot-i18n.h"
#include "lingot-msg.h"
#ifndef timeradd
#define timeradd(a, b, result) \
do { \
(result)->tv_sec = (a)->tv_sec + (b)->tv_sec; \
(result)->tv_usec = (a)->tv_usec + (b)->tv_usec; \
if ((result)->tv_usec >= 1000000) \
{ \
++(result)->tv_sec; \
(result)->tv_usec -= 1000000; \
} \
} while (0)
#endif
int
lingot_core_read_callback(FLT* read_buffer, int read_buffer_size, void *arg);
void lingot_core_run_reading_thread(LingotCore* core);
void lingot_core_run_computation_thread(LingotCore* core);
int decimation_input_index = 0;
LingotCore* lingot_core_new(LingotConfig* conf) {
char buff[1000];
LingotCore* core = malloc(sizeof(LingotCore));
core->conf = conf;
core->running = 0;
core->audio = NULL;
core->spd_fft = NULL;
core->X = NULL;
core->spd_dft = NULL;
core->diff2_spd_fft = NULL;
core->flt_read_buffer = NULL;
core->temporal_buffer = NULL;
core->windowed_temporal_buffer = NULL;
core->windowed_fft_buffer = NULL;
core->hamming_window_temporal = NULL;
core->hamming_window_fft = NULL;
core->antialiasing_filter = NULL;
int requested_sample_rate = conf->sample_rate;
if (conf->sample_rate <= 0) {
conf->sample_rate = 0;
}
core->audio = lingot_audio_new(conf->audio_system,
conf->audio_dev[conf->audio_system], conf->sample_rate,
(LingotAudioProcessCallback) lingot_core_read_callback, core);
if (core->audio != NULL) {
if (requested_sample_rate != core->audio->real_sample_rate) {
conf->sample_rate = core->audio->real_sample_rate;
lingot_config_update_internal_params(conf);
sprintf(
buff,
_("The requested sample rate is not available, the real sample rate has been set to %d Hz"),
core->audio->real_sample_rate);
lingot_msg_add_warning(buff);
}
if (conf->temporal_buffer_size < conf->fft_size) {
conf->temporal_window = ((double) conf->fft_size
* conf->oversampling) / conf->sample_rate;
conf->temporal_buffer_size = conf->fft_size;
lingot_config_update_internal_params(conf);
sprintf(
buff,
_(
"The temporal buffer is smaller than FFT size. It has been increased to %0.3f seconds"),
conf->temporal_window);
lingot_msg_add_warning(buff);
}
// since the SPD is simmetrical, we only store the 1st half.
if (core->conf->fft_size > 256) {
core->spd_fft = malloc((core->conf->fft_size >> 1) * sizeof(FLT));
core->X = malloc((core->conf->fft_size >> 1) * sizeof(FLT));
memset(core->spd_fft, 0, (core->conf->fft_size >> 1) * sizeof(FLT));
memset(core->X, 0, (core->conf->fft_size >> 1) * sizeof(FLT));
} else { // if the fft size is 256, we store the whole signal for representation.
core->spd_fft = malloc((core->conf->fft_size) * sizeof(FLT));
core->X = malloc((core->conf->fft_size) * sizeof(FLT));
memset(core->spd_fft, 0, core->conf->fft_size * sizeof(FLT));
memset(core->X, 0, (core->conf->fft_size) * sizeof(FLT));
}
core->spd_dft = malloc((core->conf->dft_size) * sizeof(FLT));
memset(core->spd_dft, 0, core->conf->dft_size * sizeof(FLT));
core->diff2_spd_fft = malloc(core->conf->fft_size * sizeof(FLT)); // 2nd derivative from SPD.
memset(core->diff2_spd_fft, 0, core->conf->fft_size * sizeof(FLT));
memset(core->spd_dft, 0, core->conf->dft_size * sizeof(FLT));
lingot_fft_create_phase_factors(conf); // creates the phase factors for FFT.
core->fft_out = malloc((core->conf->fft_size) * sizeof(LingotComplex)); // complex signal in freq domain.
memset(core->fft_out, 0, core->conf->fft_size * sizeof(LingotComplex));
// audio source read in floating point format.
core->flt_read_buffer = malloc(core->audio->read_buffer_size
* sizeof(FLT));
memset(core->flt_read_buffer, 0, core->audio->read_buffer_size
* sizeof(FLT));
// stored samples.
core->temporal_buffer = malloc((core->conf->temporal_buffer_size)
* sizeof(FLT));
memset(core->temporal_buffer, 0, core->conf->temporal_buffer_size
* sizeof(FLT));
core->hamming_window_temporal = NULL;
core->hamming_window_fft = NULL;
if (conf->window_type != NONE) {
core->hamming_window_temporal = malloc(
(core->conf->temporal_buffer_size) * sizeof(FLT));
core->hamming_window_fft = malloc((core->conf->fft_size)
* sizeof(FLT));
lingot_signal_window(core->conf->temporal_buffer_size,
core->hamming_window_temporal, conf->window_type);
lingot_signal_window(core->conf->fft_size,
core->hamming_window_fft, conf->window_type);
}
core->windowed_temporal_buffer = malloc(
(core->conf->temporal_buffer_size) * sizeof(FLT));
memset(core->windowed_temporal_buffer, 0,
core->conf->temporal_buffer_size * sizeof(FLT));
core->windowed_fft_buffer
= malloc((core->conf->fft_size) * sizeof(FLT));
memset(core->windowed_fft_buffer, 0, core->conf->fft_size * sizeof(FLT));
/*
* 8 order Chebyshev filters, with wc=0.9/i (normalized respect to
* Pi). We take 0.9 instead of 1 to leave a 10% of safety margin,
* in order to avoid aliased frequencies near to w=Pi, due to non
* ideality of the filter.
*
* The cut frequencies wc=Pi/i, with i=1..20, correspond with the
* oversampling factor, avoiding aliasing at decimation.
*
* Why Chebyshev filters?, for a given order, those filters yield
* abrupt falls than other ones as Butterworth, making the most of
* the order. Although Chebyshev filters affects more to the phase,
* it doesn't matter due to the analysis is made on the signal
* power distribution (only module).
*/
core->antialiasing_filter = lingot_filter_cheby_design(8, 0.5, 0.9
/ core->conf->oversampling);
pthread_mutex_init(&core->temporal_buffer_mutex, NULL);
// ------------------------------------------------------------
core->running = 1;
}
core->freq = 0.0;
return core;
}
// -----------------------------------------------------------------------
/* Deallocate resources */
void lingot_core_destroy(LingotCore* core) {
if (core->audio != NULL) {
lingot_fft_destroy_phase_factors(); // destroy phase factors.
free(core->fft_out);
lingot_audio_destroy(core->audio);
free(core->spd_fft);
free(core->X);
free(core->spd_dft);
free(core->diff2_spd_fft);
free(core->flt_read_buffer);
free(core->temporal_buffer);
if (core->hamming_window_fft != NULL) {
free(core->hamming_window_temporal);
}
if (core->windowed_temporal_buffer != NULL) {
free(core->windowed_temporal_buffer);
}
if (core->hamming_window_fft != NULL) {
free(core->hamming_window_fft);
}
if (core->windowed_fft_buffer != NULL) {
free(core->windowed_fft_buffer);
}
if (core->antialiasing_filter != NULL)
lingot_filter_destroy(core->antialiasing_filter);
pthread_mutex_destroy(&core->temporal_buffer_mutex);
}
free(core);
}
// -----------------------------------------------------------------------
// reads a new piece of signal from audio source, apply filtering and
// decimation and appends it to the buffer
int lingot_core_read_callback(FLT* read_buffer, int read_buffer_size, void *arg) {
unsigned int i, decimation_output_index; // loop variables.
int decimation_output_len;
FLT* decimation_in;
FLT* decimation_out;
LingotCore* core = (LingotCore*) arg;
LingotConfig* conf = core->conf;
memcpy(core->flt_read_buffer, read_buffer, read_buffer_size * sizeof(FLT));
// double freq = 100.0;
// for (i = 0; i < read_buffer_size; i++) {
// core->flt_read_buffer[i] = 5e2 * cos(2.0 * M_PI * freq * i
// / conf->sample_rate);
// }
if (conf->gain_nu != 1.0) {
for (i = 0; i < read_buffer_size; i++)
core->flt_read_buffer[i] *= conf->gain_nu;
}
//
// just readed:
//
// ----------------------------
// |bxxxbxxxbxxxbxxxbxxxbxxxbxxx|
// ----------------------------
//
// <----------------------------> read_buffer_size*oversampling
//
decimation_output_len = 1 + (read_buffer_size
- (decimation_input_index + 1)) / conf->oversampling;
pthread_mutex_lock(&core->temporal_buffer_mutex);
/* we shift the temporal window to leave a hollow where place the new piece
of data read. The buffer is actually a queue. */
if (conf->temporal_buffer_size > decimation_output_len) {
memmove(core->temporal_buffer,
&core->temporal_buffer[decimation_output_len],
(conf->temporal_buffer_size - decimation_output_len)
* sizeof(FLT));
}
//
// previous buffer situation:
//
// ------------------------------------------
// | xxxxxxxxxxxxxxxxxxxxxx | yyyyy | aaaaaaa |
// ------------------------------------------
// <------> read_buffer_size
// <---------------> fft_size
// <----------------------------------------> temporal_buffer_size
//
// new situation:
//
// ------------------------------------------
// | xxxxxxxxxxxxxxxxyyyyaa | aaaaa | |
// ------------------------------------------
//
// decimation with lowpass filtering
/* we decimate the read signal and put it at the end of the buffer. */
if (conf->oversampling > 1) {
decimation_in = core->flt_read_buffer;
decimation_out = &core->temporal_buffer[conf->temporal_buffer_size
- decimation_output_len];
// low pass filter to avoid aliasing.
lingot_filter_filter(core->antialiasing_filter, read_buffer_size,
decimation_in, decimation_in);
// compression.
for (decimation_output_index = 0; decimation_input_index
< read_buffer_size; decimation_output_index++, decimation_input_index
+= conf->oversampling)
decimation_out[decimation_output_index]
= decimation_in[decimation_input_index];
decimation_input_index -= read_buffer_size;
} else
memcpy(&core->temporal_buffer[conf->temporal_buffer_size
- decimation_output_len], core->flt_read_buffer,
decimation_output_len * sizeof(FLT));
//
// ------------------------------------------
// | xxxxxxxxxxxxxxxxyyyyaa | aaaaa | bbbbbbb |
// ------------------------------------------
//
pthread_mutex_unlock(&core->temporal_buffer_mutex);
return 0;
}
void lingot_core_compute_fundamental_fequency(LingotCore* core) {
register unsigned int i, k; // loop variables.
LingotConfig* conf = core->conf;
FLT delta_w_FFT = 2.0 * M_PI / conf->fft_size; // FFT resolution in rads.
// ----------------- TRANSFORMATION TO FREQUENCY DOMAIN ----------------
FLT _1_N2 = 1.0 / (conf->fft_size * conf->fft_size);
// SPD normalization constant
//printf("%d samples transformed of a total of %d\n", conf->fft_size,
// conf->temporal_buffer_size);
pthread_mutex_lock(&core->temporal_buffer_mutex);
// windowing
if (conf->window_type != NONE) {
for (i = 0; i < conf->fft_size; i++) {
core->windowed_fft_buffer[i]
= core->temporal_buffer[conf->temporal_buffer_size
- conf->fft_size + i] * core->hamming_window_fft[i];
}
} else {
memmove(core->windowed_fft_buffer,
&core->temporal_buffer[conf->temporal_buffer_size
- conf->fft_size], conf->fft_size * sizeof(FLT));
}
// transformation.
lingot_fft_fft(core->windowed_fft_buffer, core->fft_out, conf->fft_size);
// esteem of SPD from FFT. (normalized squared module)
for (i = 0; i < ((conf->fft_size > 256) ? (conf->fft_size >> 1) : 256); i++)
core->spd_fft[i] = (core->fft_out[i].r * core->fft_out[i].r
+ core->fft_out[i].i * core->fft_out[i].i) * _1_N2;
// representable piece
memcpy(core->X, core->spd_fft, ((conf->fft_size > 256) ? (conf->fft_size
>> 1) : 256) * sizeof(FLT));
// truncated 2nd derivative esteem, to enhance peaks
core->diff2_spd_fft[0] = 0.0;
for (i = 1; i < (conf->fft_size >> 1) - 1; i++) {
core->diff2_spd_fft[i] = 2.0 * core->spd_fft[i] - core->spd_fft[i - 1]
- core->spd_fft[i + 1]; // centred 2nd order derivative, to avoid group delay.
if (core->diff2_spd_fft[i] < 0.0)
core->diff2_spd_fft[i] = 0.0; // truncation
}
// peaks searching in that signal.
int Mi = lingot_signal_get_fundamental_peak(conf, core->spd_fft,
core->diff2_spd_fft, (conf->fft_size >> 1)); // take the fundamental peak.
if (Mi == (signed) (conf->fft_size >> 1)) {
core->freq = 0.0;
pthread_mutex_unlock(&core->temporal_buffer_mutex);
return;
}
FLT w = (Mi - 1) * delta_w_FFT;
// frequency of sample previous to the peak
// Approximation to fundamental frequency by selective DFTs
// ---------------------------------------------------------
FLT d_w = delta_w_FFT;
for (k = 0; k < conf->dft_number; k++) {
d_w = 2.0 * d_w / (conf->dft_size - 1); // resolution in rads.
if (k == 0) {
lingot_fft_spd(core->windowed_fft_buffer, conf->fft_size, w + d_w,
d_w, &core->spd_dft[1], conf->dft_size - 2);
core->spd_dft[0] = core->spd_fft[Mi - 1];
core->spd_dft[conf->dft_size - 1] = core->spd_fft[Mi + 1]; // 2 samples known.
} else
lingot_fft_spd(core->windowed_fft_buffer, conf->fft_size, w, d_w,
core->spd_dft, conf->dft_size);
lingot_signal_get_max(core->spd_dft, conf->dft_size, &Mi); // search the maximum.
w += (Mi - 1) * d_w; // previous sample to the peak.
}
w += d_w; // DFT approximation.
// windowing
if (conf->window_type != NONE) {
for (i = 0; i < conf->temporal_buffer_size; i++) {
core->windowed_temporal_buffer[i] = core->temporal_buffer[i]
* core->hamming_window_temporal[i];
}
} else {
memmove(core->windowed_temporal_buffer, core->temporal_buffer,
conf->temporal_buffer_size * sizeof(FLT));
}
pthread_mutex_unlock(&core->temporal_buffer_mutex);
// Maximum finding by Newton-Raphson
// -----------------------------------
FLT wk = -1.0e5;
FLT wkm1 = w;
// first iterator set to the current approximation.
FLT d1_SPD, d2_SPD;
for (k = 0; (k < conf->max_nr_iter) && (fabs(wk - wkm1) > 1.0e-8); k++) {
wk = wkm1;
// ! we use the WHOLE temporal window for greater precission.
lingot_fft_spd_diffs(core->windowed_temporal_buffer,
conf->temporal_buffer_size, wk, &d1_SPD, &d2_SPD);
//printf("%f %f %f\n", wk, d1_SPD, d2_SPD);
wkm1 = wk - d1_SPD / d2_SPD;
}
w = wkm1; // frequency in rads.
core->freq = (w * conf->sample_rate) / (2.0 * M_PI * conf->oversampling); // analog frequency.
}
/* start running the core in another thread */
void lingot_core_start(LingotCore* core) {
int audio_status = 0;
decimation_input_index = 0;
if (core->audio != NULL) {
audio_status = lingot_audio_start(core->audio);
if (audio_status == 0) {
pthread_mutex_init(&core->thread_computation_mutex, NULL);
pthread_cond_init(&core->thread_computation_cond, NULL);
pthread_attr_init(&core->thread_computation_attr);
pthread_create(&core->thread_computation,
&core->thread_computation_attr,
(void* (*)(void*)) lingot_core_run_computation_thread, core);
}
core->running = 1;
}
}
/* stop running the core */
void lingot_core_stop(LingotCore* core) {
void* thread_result;
if (core->running == 1) {
core->running = 0;
// threads cancelation
pthread_cancel(core->thread_computation);
pthread_join(core->thread_computation, &thread_result);
// printf("%p %p %i\n", thread_result, PTHREAD_CANCELED, thread_result
// == PTHREAD_CANCELED);
pthread_attr_destroy(&core->thread_computation_attr);
memset(core->X, 0,
((core->conf->fft_size > 256) ? (core->conf->fft_size >> 1)
: core->conf->fft_size) * sizeof(FLT));
core->freq = 0.0;
}
if (core->audio != NULL)
lingot_audio_stop(core->audio);
}
/* run the core */
void lingot_core_run_computation_thread(LingotCore* core) {
struct timeval t, tout;
struct timespec tspec;
gettimeofday(&tout, NULL);
t.tv_sec = 0;
t.tv_usec = 1e6 / core->conf->calculation_rate;
while (core->running) {
lingot_core_compute_fundamental_fequency(core);
timeradd(&t, &tout, &tout);
tspec.tv_sec = tout.tv_sec;
tspec.tv_nsec = 1000 * tout.tv_usec;
pthread_cond_timedwait(&core->thread_computation_cond,
&core->thread_computation_mutex, &tspec);
if (core->audio != NULL) {
if (core->audio->interrupted) {
memset(core->X, 0,
((core->conf->fft_size > 256) ? (core->conf->fft_size
>> 1) : core->conf->fft_size) * sizeof(FLT));
core->freq = 0.0;
core->running = 0;
}
}
}
}