Last commit for src/lingot-core.c: 9bc63f42d03ef071132142c2462b03b992b95f3b

- make OSS backend optional

Piotr Pawlow [2014-03-13 20:52:10]
- make OSS backend optional
/*
 * 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>
#ifdef OSS
#include <sys/soundcard.h>
#endif
#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;
			}
		}
	}
}
ViewGit