Split files, better peak detector

This commit is contained in:
Marcello Mamino 2015-09-29 22:28:12 +02:00
parent abfd160bb5
commit 5f836f6abf
5 changed files with 564 additions and 525 deletions

View file

@ -1,4 +1,4 @@
all: tg
tg: tg.c
gcc `pkg-config --libs --cflags gtk+-2.0 portaudio-2.0 fftw3f` -lm tg.c -O3 -fcx-limited-range -o tg
tg: interface.c algo.c tg.h
gcc -Wall `pkg-config --libs --cflags gtk+-2.0 portaudio-2.0 fftw3f` -lm -O3 -fcx-limited-range -o tg interface.c algo.c audio.c

380
algo.c Normal file
View file

@ -0,0 +1,380 @@
#include "tg.h"
int fl_cmp(const void *a, const void *b)
{
float x = *(float*)a;
float y = *(float*)b;
return x<y ? -1 : x>y ? 1 : 0;
}
void make_hp(struct filter *f, double freq)
{
double K = tan(M_PI * freq);
double norm = 1 / (1 + K * sqrt(2) + K * K);
f->a0 = 1 * norm;
f->a1 = -2 * f->a0;
f->a2 = f->a0;
f->b1 = 2 * (K * K - 1) * norm;
f->b2 = (1 - K * sqrt(2) + K * K) * norm;
}
void make_lp(struct filter *f, double freq)
{
double K = tan(M_PI * freq);
double norm = 1 / (1 + K * sqrt(2) + K * K);
f->a0 = K * K * norm;
f->a1 = 2 * f->a0;
f->a2 = f->a0;
f->b1 = 2 * (K * K - 1) * norm;
f->b2 = (1 - K * sqrt(2) + K * K) * norm;
}
void run_filter(struct filter *f, float *buff, int size)
{
int i;
double z1 = 0, z2 = 0;
for(i=0; i<size; i++) {
double in = buff[i];
double out = in * f->a0 + z1;
z1 = in * f->a1 + z2 - f->b1 * out;
z2 = in * f->a2 - f->b2 * out;
buff[i] = out;
}
}
void setup_buffers(struct processing_buffers *b)
{
b->samples = fftwf_malloc(2 * b->sample_count * sizeof(float));
b->samples_sc = malloc(2 * b->sample_count * sizeof(float));
b->waveform = malloc(2 * b->sample_rate * sizeof(float));
b->waveform_sc = malloc(2 * b->sample_rate * sizeof(float));
b->fft = malloc((b->sample_count + 1) * sizeof(fftwf_complex));
b->plan_a = fftwf_plan_dft_r2c_1d(2 * b->sample_count, b->samples, b->fft, FFTW_ESTIMATE);
b->plan_b = fftwf_plan_dft_c2r_1d(2 * b->sample_count, b->fft, b->samples_sc, FFTW_ESTIMATE);
b->plan_c = fftwf_plan_dft_r2c_1d(2 * b->sample_rate, b->waveform, b->fft, FFTW_ESTIMATE);
b->plan_d = fftwf_plan_dft_c2r_1d(2 * b->sample_rate, b->fft, b->waveform_sc, FFTW_ESTIMATE);
b->hpf = malloc(sizeof(struct filter));
make_hp(b->hpf,(double)FILTER_CUTOFF/b->sample_rate);
b->lpf = malloc(sizeof(struct filter));
make_lp(b->lpf,(double)FILTER_CUTOFF/b->sample_rate);
b->ready = 0;
}
struct processing_buffers *pb_clone(struct processing_buffers *p)
{
struct processing_buffers *new = malloc(sizeof(struct processing_buffers));
new->sample_rate = p->sample_rate;
new->waveform = malloc(new->sample_rate * sizeof(float));
memcpy(new->waveform, p->waveform, new->sample_rate * sizeof(float));
new->sample_count = p->sample_count;
new->samples_sc = malloc(new->sample_count * sizeof(float));
memcpy(new->samples_sc, p->samples_sc, new->sample_rate * sizeof(float));
new->period = p->period;
new->sigma = p->sigma;
new->be = p->be;
new->waveform_max = p->waveform_max;
new->tic = p->tic;
new->toc = p->toc;
new->ready = p->ready;
return new;
}
void pb_destroy(struct processing_buffers *p)
{
// Future BUG: we free only waveform
free(p->waveform);
free(p->samples_sc);
free(p);
}
void noise_suppressor(struct processing_buffers *p)
{
float *a = p->samples_sc;
float *b = p->samples_sc + p->sample_count;
int window = p->sample_rate / 50;
int i;
for(i = 0; i < p->sample_count; i++)
a[i] = p->samples[i] * p->samples[i];
double r_av = 0;
for(i = 0; i < window; i++)
r_av += a[i];
for(i = 0;; i++) {
b[i] = r_av;
if(i + window == p->sample_count) break;
r_av += a[i + window] - a[i];
}
int m = p->sample_count - window + 1;
float max = 0;
int j = 0;
for(i = 0; i < m; i++) {
if(b[i] > max) max = b[i];
if((i+1) % (p->sample_rate/2) == 0) {
a[j++] = max;
max = 0;
}
}
qsort(a, j, sizeof(float), fl_cmp);
float k = a[j/2];
for(i = 0; i < p->sample_count; i++) {
int j = i - window / 2;
j = j < 0 ? 0 : j > p->sample_count - window ? p->sample_count - window : j;
if(b[j] > 2*k) p->samples[i] = 0;
}
}
void prepare_data(struct processing_buffers *b)
{
int i;
int first_fft_size = b->sample_count/2 + 1;
run_filter(b->hpf, b->samples, b->sample_count);
noise_suppressor(b);
for(i=0; i < b->sample_count; i++)
b->samples[i] = fabs(b->samples[i]);
run_filter(b->lpf, b->samples, b->sample_count);
double average = 0;
for(i=0; i < b->sample_count; i++)
average += b->samples[i];
average /= b->sample_count;
for(i=0; i < b->sample_count; i++)
b->samples[i] -= average;
for(i=0; i < b->sample_rate/10; i++) {
double k = ( 1 - cos(i*M_PI/(b->sample_rate/10)) ) / 2;
b->samples[i] *= k;
b->samples[b->sample_count - i - 1] *= k;
}
fftwf_execute(b->plan_a);
for(i=0; i < b->sample_count+1; i++)
b->fft[i] = b->fft[i] * conj(b->fft[i]);
fftwf_execute(b->plan_b);
}
int peak_detector(float *buff, struct processing_buffers *p, int a, int b)
{
int i;
double max = buff[a];
int i_max = a;
for(i=a+1; i<=b; i++) {
if(buff[i] > max) {
max = buff[i];
i_max = i;
}
}
if(max <= 0) return -1;
float v[b-a+1];
for(i=a; i<=b; i++)
v[i-a] = buff[i];
qsort(v, b-a+1, sizeof(float), fl_cmp);
float x = v[(b-a+1)/10];
float y = v[(b-a+1)/2];
if(max < 5*y - 4*x)
return -1;
return i_max;
}
double estimate_period(struct processing_buffers *p)
{
int estimate = peak_detector(p->samples_sc, p, p->sample_rate / 12, p->sample_rate / 2);
if(estimate == -1) {
debug("no candidate period\n");
return -1;
}
int a = estimate*3/2 - p->sample_rate / 50;
int b = estimate*3/2 + p->sample_rate / 50;
double max = p->samples_sc[a];
int i;
for(i=a+1; i<=b; i++)
if(p->samples_sc[i] > max)
max = p->samples_sc[i];
if(max < 0.2 * p->samples_sc[estimate]) {
if(estimate * 2 < p->sample_rate / 2) {
debug("double triggered\n");
return peak_detector(p->samples_sc, p,
estimate*2 - p->sample_rate / 50,
estimate*2 + p->sample_rate / 50);
} else {
debug("period rejected (immense beat error?)");
return -1;
}
} else return estimate;
}
int compute_period(struct processing_buffers *b, int bph)
{
double estimate;
if(bph)
estimate = peak_detector(b->samples_sc, b,
7200 * b->sample_rate / bph - b->sample_rate / 100,
7200 * b->sample_rate / bph + b->sample_rate / 100);
else
estimate = estimate_period(b);
if(estimate == -1) {
debug("failed to estimate period\n");
return 1;
}
double delta = b->sample_rate * 0.01;
double new_estimate = estimate;
double sum = 0;
double sq_sum = 0;
int count = 0;
int cycle = 1;
for(;;) {
int inf = floor(new_estimate * cycle - delta);
int sup = ceil(new_estimate * cycle + delta);
if(sup > b->sample_count * 2 / 3)
break;
new_estimate = peak_detector(b->samples_sc,b,inf,sup);
if(new_estimate == -1) {
debug("cycle = %d peak not found\n",cycle);
return 1;
}
new_estimate /= cycle;
if(new_estimate < estimate - delta || new_estimate > estimate + delta) {
debug("cycle = %d new_estimate = %f invalid peak\n",cycle,new_estimate/b->sample_rate);
return 1;
} else
debug("cycle = %d new_estimate = %f\n",cycle,new_estimate/b->sample_rate);
if(inf > b->sample_count / 3) {
sum += new_estimate;
sq_sum += new_estimate * new_estimate;
count++;
}
cycle++;
}
if(count > 0) estimate = sum / count;
if(estimate >= b->sample_rate / 2) return 1;
b->period = estimate;
if(count > 1)
b->sigma = sqrt((sq_sum - count * estimate * estimate)/ (count-1));
else
b->sigma = b->period;
return 0;
}
float tmean(float *x, int n)
{
qsort(x,n,sizeof(float),fl_cmp);
int i;
double sum = 0;
for(i=0; i < n*4/5; i++)
sum += x[i];
return sum/(n*4/5);
}
int compute_parameters(struct processing_buffers *p)
{
int i;
double x = 0, y = 0;
for(i=0; i<p->sample_count; i++) {
double a = i * 4 * M_PI / p->period;
x += p->samples[i] * cos(a);
y += p->samples[i] * sin(a);
}
double s = p->period * (M_PI + atan2(y,x)) / (4 * M_PI);
for(i=0; i<2*p->sample_rate; i++)
p->waveform[i] = 0;
float bin[(int)ceil(1 + p->sample_count / p->period)];
for(i=0; i < p->period; i++) {
int j;
double k = fmod(i+s,p->period);
for(j=0;;j++) {
int n = round(k+j*p->period);
if(n >= p->sample_count) break;
bin[j] = p->samples[n];
}
p->waveform[i] = tmean(bin,j);
}
for(i=0; i<p->period; i++)
p->waveform_sc[i] = p->waveform[i];
qsort(p->waveform_sc,floor(p->period),sizeof(float),fl_cmp);
double nl = p->waveform_sc[(int)floor(p->period/2)];
for(i=0; i<p->period; i++) {
p->waveform[i] -= nl;
//p->waveform[i] = cbrt(p->waveform[i]);
//if(p->waveform[i] <= 0) p->waveform[i] = 0;
}
fftwf_execute(p->plan_c);
for(i=0; i < p->sample_rate+1; i++)
p->fft[i] = p->fft[i] * conj(p->fft[i]);
fftwf_execute(p->plan_d);
for(i=0;i<p->period;i++)
p->samples_sc[i] = p->waveform_sc[i];
int tic_to_toc = peak_detector(p->waveform_sc,p,
floor(p->period/2)-p->sample_rate/50,
floor(p->period/2)+p->sample_rate/50);
if(tic_to_toc < 0) {
p->tic = p->toc = -1;
p->be = -1;
debug("beat error = ---\n");
return 1;
} else {
p->be = p->period/2 - tic_to_toc;
debug("beat error = %.1f\n",fabs(p->be)*1000/p->sample_rate);
}
double max = 0;
int max_i = -1;
for(i=0;i<p->period;i++) {
if(p->waveform[i] > max) {
max = p->waveform[i];
max_i = i;
}
}
p->waveform_max = max;
if(max_i < p->period/2) {
p->tic = max_i;
p->toc = max_i + tic_to_toc;
if(p->toc > p->period)
p->toc = floor(p->period);
} else {
p->toc = max_i;
p->tic = max_i - tic_to_toc;
if(p->tic < 0)
p->tic = 0;
}
double phase = p->timestamp - p->last_tic;
double apparent_phase = p->sample_count - (s + p->tic);
double shift = fmod(apparent_phase - phase, p->period);
if(shift < 0) shift += p->period;
debug("shift = %.3f\n",shift / p->period);
if(fabs(shift - p->period/2) < p->period/4) {
p->last_toc = p->timestamp - (uint64_t)round(fmod(apparent_phase, p->period));
int t = p->tic;
p->tic = p->toc;
p->toc = t;
apparent_phase = p->sample_count - (s + p->tic);
} else
p->last_toc = p->timestamp - (uint64_t)round(fmod(p->sample_count - (s + p->toc), p->period));
p->last_tic = p->timestamp - (uint64_t)round(fmod(apparent_phase, p->period));
return 0;
}
void process(struct processing_buffers *p, int bph)
{
int i;
prepare_data(p);
p->ready = ! ( compute_period(p,bph) || compute_parameters(p) );
}

82
audio.c Normal file
View file

@ -0,0 +1,82 @@
#include "tg.h"
volatile float pa_buffers[2][PA_BUFF_SIZE];
volatile int write_pointer = 0;
volatile uint64_t timestamp = 0;
int paudio_callback(const void *input_buffer,
void *output_buffer,
unsigned long frame_count,
const PaStreamCallbackTimeInfo *time_info,
PaStreamCallbackFlags status_flags,
void *data)
{
unsigned long i;
for(i=0; i < frame_count; i++) {
pa_buffers[0][write_pointer] = ((float *)input_buffer)[2*i];
pa_buffers[1][write_pointer] = ((float *)input_buffer)[2*i + 1];
if(write_pointer < PA_BUFF_SIZE - 1) write_pointer++;
else write_pointer = 0;
timestamp++;
}
return 0;
}
int start_portaudio()
{
PaStream *stream;
PaStream **x = malloc(sizeof(PaStream*));
PaError err = Pa_Initialize();
if(err!=paNoError)
goto error;
err = Pa_OpenDefaultStream(&stream,2,0,paFloat32,PA_SAMPLE_RATE,paFramesPerBufferUnspecified,paudio_callback,x);
*x = stream;
if(err!=paNoError)
goto error;
err = Pa_StartStream(stream);
if(err!=paNoError)
goto error;
return 0;
error:
error("Error opening audio input: %s\n", Pa_GetErrorText(err));
return 1;
}
int analyze_pa_data(struct processing_buffers *p, int bph)
{
static uint64_t last_tic = 0;
int wp = write_pointer;
uint64_t ts = timestamp;
if(wp < 0 || wp >= PA_BUFF_SIZE) wp = 0;
int i;
for(i=0; i<NSTEPS; i++) {
int j,k;
memset(p[i].samples,0,2 * p[i].sample_count * sizeof(float));
k = wp - p[i].sample_count;
if(k < 0) k += PA_BUFF_SIZE;
for(j=0; j < p[i].sample_count; j++) {
// p[i].samples[j] = pa_buffers[0][k] + pa_buffers[1][k];
p[i].samples[j] = pa_buffers[1][k];
if(++k == PA_BUFF_SIZE) k = 0;
}
}
for(i=0; i<NSTEPS; i++) {
p[i].timestamp = ts;
p[i].last_tic = last_tic;
process(&p[i],bph);
if( !p[i].ready ) break;
debug("step %d : %f +- %f\n",i,p[i].period/p[i].sample_rate,p[i].sigma/p[i].sample_rate);
}
if(i) {
last_tic = p[i-1].last_tic;
debug("%f +- %f\n",p[i-1].period/p[i-1].sample_rate,p[i-1].sigma/p[i-1].sample_rate);
} else
debug("---\n");
return i;
}

View file

@ -1,46 +1,6 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <complex.h>
#include <fftw3.h>
#include <stdint.h>
#include <portaudio.h>
#include <stdarg.h>
#include <gtk/gtk.h>
#include "tg.h"
#include <assert.h>
#define FILTER_CUTOFF 3000
#define NSTEPS 4
#define FIRST_STEP 1
#define PA_SAMPLE_RATE 44100
#define PA_BUFF_SIZE (PA_SAMPLE_RATE << (NSTEPS + FIRST_STEP))
#define FPS 2
#define AMPLITUDE_WIDTH .07
#define OUTPUT_FONT 50
#define OUTPUT_WINDOW_HEIGHT 80
#define POSITIVE_SPAN 10
#define NEGATIVE_SPAN 25
#define EVENTS_COUNT 10000
#define MIN_BPH 12000
#define MAX_BPH 36000
#define DEFAULT_BPH 21600
#define MIN_LA 10
#define MAX_LA 90
#define DEFAULT_LA 52
int preset_bph[] = { 12000, 14400, 18000, 19800, 21600, 25200, 28800, 36000, 0 };
volatile float pa_buffers[2][PA_BUFF_SIZE];
volatile int write_pointer = 0;
volatile uint64_t timestamp = 0;
int preset_bph[] = PRESET_BPH;
void debug(char *format,...)
{
@ -58,478 +18,6 @@ void error(char *format,...)
va_end(args);
}
int paudio_callback(const void *input_buffer,
void *output_buffer,
unsigned long frame_count,
const PaStreamCallbackTimeInfo *time_info,
PaStreamCallbackFlags status_flags,
void *data)
{
unsigned long i;
for(i=0; i < frame_count; i++) {
pa_buffers[0][write_pointer] = ((float *)input_buffer)[2*i];
pa_buffers[1][write_pointer] = ((float *)input_buffer)[2*i + 1];
/*
if(timestamp % 4410 > 4000) {
pa_buffers[0][write_pointer] = 1; //timestamp % 2;
pa_buffers[1][write_pointer] = 1; //timestamp % 2;
} else {
pa_buffers[0][write_pointer] = 0;
pa_buffers[1][write_pointer] = 0;
}
*/
if(write_pointer < PA_BUFF_SIZE - 1) write_pointer++;
else write_pointer = 0;
timestamp++;
}
return 0;
}
int start_portaudio()
{
PaStream *stream;
PaStream **x = malloc(sizeof(PaStream*));
PaError err = Pa_Initialize();
if(err!=paNoError)
goto error;
err = Pa_OpenDefaultStream(&stream,2,0,paFloat32,PA_SAMPLE_RATE,paFramesPerBufferUnspecified,paudio_callback,x);
*x = stream;
if(err!=paNoError)
goto error;
err = Pa_StartStream(stream);
if(err!=paNoError)
goto error;
return 0;
error:
error("Error opening audio input: %s\n", Pa_GetErrorText(err));
return 1;
}
int fl_cmp(const void *a, const void *b)
{
float x = *(float*)a;
float y = *(float*)b;
return x<y ? -1 : x>y ? 1 : 0;
}
struct filter {
double a0,a1,a2,b1,b2;
};
void make_hp(struct filter *f, double freq)
{
double K = tan(M_PI * freq);
double norm = 1 / (1 + K * sqrt(2) + K * K);
f->a0 = 1 * norm;
f->a1 = -2 * f->a0;
f->a2 = f->a0;
f->b1 = 2 * (K * K - 1) * norm;
f->b2 = (1 - K * sqrt(2) + K * K) * norm;
}
void make_lp(struct filter *f, double freq)
{
double K = tan(M_PI * freq);
double norm = 1 / (1 + K * sqrt(2) + K * K);
f->a0 = K * K * norm;
f->a1 = 2 * f->a0;
f->a2 = f->a0;
f->b1 = 2 * (K * K - 1) * norm;
f->b2 = (1 - K * sqrt(2) + K * K) * norm;
}
void run_filter(struct filter *f, float *buff, int size)
{
int i;
double z1 = 0, z2 = 0;
for(i=0; i<size; i++) {
double in = buff[i];
double out = in * f->a0 + z1;
z1 = in * f->a1 + z2 - f->b1 * out;
z2 = in * f->a2 - f->b2 * out;
buff[i] = out;
}
}
struct processing_buffers {
int sample_rate;
int sample_count;
float *samples, *samples_sc, *waveform, *waveform_sc;
fftwf_complex *fft;
fftwf_plan plan_a, plan_b, plan_c, plan_d;
struct filter *hpf, *lpf;
double period,sigma,be,waveform_max;
int tic,toc;
int ready;
uint64_t timestamp, last_tic, last_toc;
};
void setup_buffers(struct processing_buffers *b)
{
b->samples = fftwf_malloc(2 * b->sample_count * sizeof(float));
b->samples_sc = malloc(2 * b->sample_count * sizeof(float));
b->waveform = malloc(2 * b->sample_rate * sizeof(float));
b->waveform_sc = malloc(2 * b->sample_rate * sizeof(float));
b->fft = malloc((b->sample_count + 1) * sizeof(fftwf_complex));
b->plan_a = fftwf_plan_dft_r2c_1d(2 * b->sample_count, b->samples, b->fft, FFTW_ESTIMATE);
b->plan_b = fftwf_plan_dft_c2r_1d(2 * b->sample_count, b->fft, b->samples_sc, FFTW_ESTIMATE);
b->plan_c = fftwf_plan_dft_r2c_1d(2 * b->sample_rate, b->waveform, b->fft, FFTW_ESTIMATE);
b->plan_d = fftwf_plan_dft_c2r_1d(2 * b->sample_rate, b->fft, b->waveform_sc, FFTW_ESTIMATE);
b->hpf = malloc(sizeof(struct filter));
make_hp(b->hpf,(double)FILTER_CUTOFF/b->sample_rate);
b->lpf = malloc(sizeof(struct filter));
make_lp(b->lpf,(double)FILTER_CUTOFF/b->sample_rate);
b->ready = 0;
}
struct processing_buffers *pb_clone(struct processing_buffers *p)
{
struct processing_buffers *new = malloc(sizeof(struct processing_buffers));
new->sample_rate = p->sample_rate;
new->waveform = malloc(new->sample_rate * sizeof(float));
memcpy(new->waveform, p->waveform, new->sample_rate * sizeof(float));
new->sample_count = p->sample_count;
new->samples_sc = malloc(new->sample_count * sizeof(float));
memcpy(new->samples_sc, p->samples_sc, new->sample_rate * sizeof(float));
new->period = p->period;
new->sigma = p->sigma;
new->be = p->be;
new->waveform_max = p->waveform_max;
new->tic = p->tic;
new->toc = p->toc;
new->ready = p->ready;
return new;
}
void pb_destroy(struct processing_buffers *p)
{
// Future BUG: we free only waveform
free(p->waveform);
free(p->samples_sc);
free(p);
}
void noise_suppressor(struct processing_buffers *p)
{
float *a = p->samples_sc;
float *b = p->samples_sc + p->sample_count;
int window = p->sample_rate / 50;
int i;
for(i = 0; i < p->sample_count; i++)
a[i] = p->samples[i] * p->samples[i];
double r_av = 0;
for(i = 0; i < window; i++)
r_av += a[i];
for(i = 0;; i++) {
b[i] = r_av;
if(i + window == p->sample_count) break;
r_av += a[i + window] - a[i];
}
int m = p->sample_count - window + 1;
float max = 0;
int j = 0;
for(i = 0; i < m; i++) {
if(b[i] > max) max = b[i];
if((i+1) % (p->sample_rate/2) == 0) {
a[j++] = max;
max = 0;
}
}
qsort(a, j, sizeof(float), fl_cmp);
float k = a[j/2];
for(i = 0; i < p->sample_count; i++) {
int j = i - window / 2;
j = j < 0 ? 0 : j > p->sample_count - window ? p->sample_count - window : j;
if(b[j] > 2*k) p->samples[i] = 0;
}
}
void prepare_data(struct processing_buffers *b)
{
int i;
int first_fft_size = b->sample_count/2 + 1;
run_filter(b->hpf, b->samples, b->sample_count);
noise_suppressor(b);
for(i=0; i < b->sample_count; i++)
b->samples[i] = fabs(b->samples[i]);
run_filter(b->lpf, b->samples, b->sample_count);
double average = 0;
for(i=0; i < b->sample_count; i++)
average += b->samples[i];
average /= b->sample_count;
for(i=0; i < b->sample_count; i++)
b->samples[i] -= average;
for(i=0; i < b->sample_rate/10; i++) {
double k = ( 1 - cos(i*M_PI/(b->sample_rate/10)) ) / 2;
b->samples[i] *= k;
b->samples[b->sample_count - i - 1] *= k;
}
fftwf_execute(b->plan_a);
for(i=0; i < b->sample_count+1; i++)
b->fft[i] = b->fft[i] * conj(b->fft[i]);
fftwf_execute(b->plan_b);
}
int peak_detector(float *buff, struct processing_buffers *p, int a, int b)
{
int i;
double max = buff[a];
int i_max = a;
for(i=a+1; i<=b; i++) {
if(buff[i] > max) {
max = buff[i];
i_max = i;
}
}
if(max <= 0) return -1;
int x,y;
for(x = i_max; x >= a && buff[x] > 0.7*max; x--);
for(y = i_max; y <= b && buff[y] > 0.7*max; y++);
if( x < a || y > b || y-x < p->sample_rate / FILTER_CUTOFF) return -1;
return i_max;
}
double estimate_period(struct processing_buffers *p)
{
int estimate = peak_detector(p->samples_sc, p, p->sample_rate / 12, p->sample_rate / 2);
if(estimate == -1) {
debug("no candidate period\n");
return -1;
}
int a = estimate*3/2 - p->sample_rate / 50;
int b = estimate*3/2 + p->sample_rate / 50;
double max = p->samples_sc[a];
int i;
for(i=a+1; i<=b; i++)
if(p->samples_sc[i] > max)
max = p->samples_sc[i];
if(max < 0.2 * p->samples_sc[estimate]) {
if(estimate * 2 < p->sample_rate / 2) {
debug("double triggered\n");
return peak_detector(p->samples_sc, p,
estimate*2 - p->sample_rate / 50,
estimate*2 + p->sample_rate / 50);
} else {
debug("period rejected (immense beat error?)");
return -1;
}
} else return estimate;
}
int compute_period(struct processing_buffers *b, int bph)
{
double estimate;
if(bph)
estimate = peak_detector(b->samples_sc, b,
7200 * b->sample_rate / bph - b->sample_rate / 100,
7200 * b->sample_rate / bph + b->sample_rate / 100);
else
estimate = estimate_period(b);
if(estimate == -1) {
debug("failed to estimate period\n");
return 1;
}
double delta = b->sample_rate * 0.01;
double new_estimate = estimate;
double sum = 0;
double sq_sum = 0;
int count = 0;
int cycle = 1;
for(;;) {
int inf = floor(new_estimate * cycle - delta);
int sup = ceil(new_estimate * cycle + delta);
if(sup > b->sample_count * 2 / 3)
break;
new_estimate = peak_detector(b->samples_sc,b,inf,sup);
if(new_estimate == -1) {
debug("cycle = %d peak not found\n",cycle);
return 1;
}
new_estimate /= cycle;
if(new_estimate < estimate - delta || new_estimate > estimate + delta) {
debug("cycle = %d new_estimate = %f invalid peak\n",cycle,new_estimate/b->sample_rate);
return 1;
} else
debug("cycle = %d new_estimate = %f\n",cycle,new_estimate/b->sample_rate);
if(inf > b->sample_count / 3) {
sum += new_estimate;
sq_sum += new_estimate * new_estimate;
count++;
}
cycle++;
}
if(count > 0) estimate = sum / count;
if(estimate >= b->sample_rate / 2) return 1;
b->period = estimate;
if(count > 1)
b->sigma = sqrt((sq_sum - count * estimate * estimate)/ (count-1));
else
b->sigma = b->period;
return 0;
}
float tmean(float *x, int n)
{
qsort(x,n,sizeof(float),fl_cmp);
int i;
double sum = 0;
for(i=0; i < n*4/5; i++)
sum += x[i];
return sum/(n*4/5);
}
int compute_parameters(struct processing_buffers *p)
{
int i;
double x = 0, y = 0;
for(i=0; i<p->sample_count; i++) {
double a = i * 4 * M_PI / p->period;
x += p->samples[i] * cos(a);
y += p->samples[i] * sin(a);
}
double s = p->period * (M_PI + atan2(y,x)) / (4 * M_PI);
for(i=0; i<2*p->sample_rate; i++)
p->waveform[i] = 0;
float bin[(int)ceil(1 + p->sample_count / p->period)];
for(i=0; i < p->period; i++) {
int j;
double k = fmod(i+s,p->period);
for(j=0;;j++) {
int n = round(k+j*p->period);
if(n >= p->sample_count) break;
bin[j] = p->samples[n];
}
p->waveform[i] = tmean(bin,j);
}
fftwf_execute(p->plan_c);
for(i=0; i < p->sample_rate+1; i++)
p->fft[i] = p->fft[i] * conj(p->fft[i]);
fftwf_execute(p->plan_d);
int tic_to_toc = peak_detector(p->waveform_sc,p,
floor(p->period/2)-p->sample_rate/50,
floor(p->period/2)+p->sample_rate/50);
if(tic_to_toc < 0) {
p->tic = p->toc = -1;
p->be = -1;
debug("beat error = ---\n");
return 1;
} else {
p->be = p->period/2 - tic_to_toc;
debug("beat error = %.1f\n",fabs(p->be)*1000/p->sample_rate);
}
for(i=0; i<p->period; i++)
p->waveform_sc[i] = p->waveform[i];
qsort(p->waveform_sc,floor(p->period),sizeof(float),fl_cmp);
double nl = p->waveform_sc[(int)floor(p->period/2)];
for(i=0; i<p->period; i++)
p->waveform[i] -= nl;
double max = 0;
int max_i = -1;
for(i=0;i<p->period;i++) {
if(p->waveform[i] > max) {
max = p->waveform[i];
max_i = i;
}
}
p->waveform_max = max;
if(max_i < p->period/2) {
p->tic = max_i;
p->toc = max_i + tic_to_toc;
if(p->toc > p->period)
p->toc = floor(p->period);
} else {
p->toc = max_i;
p->tic = max_i - tic_to_toc;
if(p->tic < 0)
p->tic = 0;
}
double phase = p->timestamp - p->last_tic;
double apparent_phase = p->sample_count - (s + p->tic);
double shift = fmod(apparent_phase - phase, p->period);
if(shift < 0) shift += p->period;
debug("shift = %.3f\n",shift / p->period);
if(fabs(shift - p->period/2) < p->period/4) {
p->last_toc = p->timestamp - (uint64_t)round(fmod(apparent_phase, p->period));
int t = p->tic;
p->tic = p->toc;
p->toc = t;
apparent_phase = p->sample_count - (s + p->tic);
} else
p->last_toc = p->timestamp - (uint64_t)round(fmod(p->sample_count - (s + p->toc), p->period));
p->last_tic = p->timestamp - (uint64_t)round(fmod(apparent_phase, p->period));
return 0;
}
void process(struct processing_buffers *p, int bph)
{
int i;
prepare_data(p);
p->ready = ! ( compute_period(p,bph) || compute_parameters(p) );
}
int analyze_pa_data(struct processing_buffers *p, int bph)
{
static uint64_t last_tic = 0;
int wp = write_pointer;
uint64_t ts = timestamp;
if(wp < 0 || wp >= PA_BUFF_SIZE) wp = 0;
int i;
for(i=0; i<NSTEPS; i++) {
int j,k;
memset(p[i].samples,0,2 * p[i].sample_count * sizeof(float));
k = wp - p[i].sample_count;
if(k < 0) k += PA_BUFF_SIZE;
for(j=0; j < p[i].sample_count; j++) {
// p[i].samples[j] = pa_buffers[0][k] + pa_buffers[1][k];
p[i].samples[j] = pa_buffers[1][k];
if(++k == PA_BUFF_SIZE) k = 0;
}
}
for(i=0; i<NSTEPS; i++) {
p[i].timestamp = ts;
p[i].last_tic = last_tic;
process(&p[i],bph);
if( !p[i].ready ) break;
debug("step %d : %f +- %f\n",i,p[i].period/p[i].sample_rate,p[i].sigma/p[i].sample_rate);
}
if(i) {
last_tic = p[i-1].last_tic;
debug("%f +- %f\n",p[i-1].period/p[i-1].sample_rate,p[i-1].sigma/p[i-1].sample_rate);
} else
debug("---\n");
return i;
}
struct main_window {
GtkWidget *window;
GtkWidget *bph_combo_box;
@ -622,8 +110,15 @@ void draw_amp_graph(double a, double b, cairo_t *c, struct processing_buffers *p
int bi = 1+round(b);
if(ai < 0) ai = 0;
if(bi > p->sample_count) bi = p->sample_count;
for(i=ai; i<bi; i++)
if(p->samples_sc[i] > max) max = p->samples_sc[i];
int qq = 0;
for(i=ai; i<bi; i++) {
if(p->samples_sc[i] > max) {
max = p->samples_sc[i];
qq = i;
}
}
debug("+++ max = %f (%d %d %d)\n",max,ai,qq,bi);
int first = 1;
for(i=0; i<width; i++) {
@ -666,9 +161,9 @@ int guess_bph(double period)
return preset_bph[ret];
}
void draw_watch_icon(cairo_t *c, int happy)
double draw_watch_icon(cairo_t *c, int signal)
{
happy = !!happy;
int happy = !!signal;
cairo_set_line_width(c,3);
cairo_set_source(c,happy?green:red);
cairo_move_to(c, OUTPUT_WINDOW_HEIGHT * 0.5, OUTPUT_WINDOW_HEIGHT * 0.5);
@ -678,6 +173,19 @@ void draw_watch_icon(cairo_t *c, int happy)
cairo_stroke(c);
cairo_arc(c, OUTPUT_WINDOW_HEIGHT * 0.5, OUTPUT_WINDOW_HEIGHT * 0.5, OUTPUT_WINDOW_HEIGHT * 0.4, 0, 2*M_PI);
cairo_stroke(c);
const int l = OUTPUT_WINDOW_HEIGHT * 0.8 / (2*NSTEPS - 1);
int i;
cairo_set_line_width(c,1);
for(i = 0; i < signal; i++) {
cairo_move_to(c, OUTPUT_WINDOW_HEIGHT + 0.5*l, OUTPUT_WINDOW_HEIGHT * 0.9 - 2*i*l);
cairo_line_to(c, OUTPUT_WINDOW_HEIGHT + 1.5*l, OUTPUT_WINDOW_HEIGHT * 0.9 - 2*i*l);
cairo_line_to(c, OUTPUT_WINDOW_HEIGHT + 1.5*l, OUTPUT_WINDOW_HEIGHT * 0.9 - (2*i+1)*l);
cairo_line_to(c, OUTPUT_WINDOW_HEIGHT + 0.5*l, OUTPUT_WINDOW_HEIGHT * 0.9 - (2*i+1)*l);
cairo_line_to(c, OUTPUT_WINDOW_HEIGHT + 0.5*l, OUTPUT_WINDOW_HEIGHT * 0.9 - 2*i*l);
cairo_stroke_preserve(c);
cairo_fill(c);
}
return OUTPUT_WINDOW_HEIGHT + 3*l;
}
gboolean output_expose_event(GtkWidget *widget, GdkEvent *event, struct main_window *w)
@ -690,11 +198,11 @@ gboolean output_expose_event(GtkWidget *widget, GdkEvent *event, struct main_win
cairo_set_source(c,black);
cairo_paint(c);
draw_watch_icon(c,w->signal);
int old;
struct processing_buffers *p = w->get_data(w,&old);
double x = draw_watch_icon(c,w->signal);
char rates[100];
char bphs[100];
@ -719,7 +227,6 @@ gboolean output_expose_event(GtkWidget *widget, GdkEvent *event, struct main_win
cairo_text_extents_t extents;
cairo_text_extents(c,"0",&extents);
double x = OUTPUT_WINDOW_HEIGHT + (double)OUTPUT_FONT/2;
double y = (double)OUTPUT_WINDOW_HEIGHT/2 - extents.y_bearing - extents.height/2;
cairo_move_to(c,x,y);
@ -769,8 +276,8 @@ gboolean amp_expose_event(GtkWidget *widget, GdkEvent *event, struct main_window
if(p) {
double span_time = p->period * span;
double a = p->period - span_time;
double b = p->period + span_time;
double a = p->period / 10; /// 2 - span_time;
double b = p->period - 1; // / 2 + span_time;
draw_amp_graph(a,b,c,p,w->amp_drawing_area);
@ -954,6 +461,8 @@ gboolean waveform_expose_event(GtkWidget *widget, GdkEvent *event, struct main_w
return FALSE;
}
extern volatile uint64_t timestamp;
gboolean paperstrip_expose_event(GtkWidget *widget, GdkEvent *event, struct main_window *w)
{
int i,old;
@ -1204,6 +713,7 @@ void int_recompute(struct main_window *w)
w->signal = analyze_pa_data(p, w->bph);
int old;
p = int_get_data(w,&old);
if(old) w->signal = -w->signal;
if(p)
w->guessed_bph = w->bph ? w->bph : guess_bph(p->period / p->sample_rate);
}

67
tg.h Normal file
View file

@ -0,0 +1,67 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <complex.h>
#include <fftw3.h>
#include <stdint.h>
#include <portaudio.h>
#include <stdarg.h>
#include <gtk/gtk.h>
#define FILTER_CUTOFF 3000
#define NSTEPS 4
#define FIRST_STEP 1
#define PA_SAMPLE_RATE 44100
#define PA_BUFF_SIZE (PA_SAMPLE_RATE << (NSTEPS + FIRST_STEP))
#define FPS 2
#define AMPLITUDE_WIDTH .07
#define OUTPUT_FONT 50
#define OUTPUT_WINDOW_HEIGHT 80
#define POSITIVE_SPAN 10
#define NEGATIVE_SPAN 25
#define EVENTS_COUNT 10000
#define MIN_BPH 12000
#define MAX_BPH 36000
#define DEFAULT_BPH 21600
#define MIN_LA 10
#define MAX_LA 90
#define DEFAULT_LA 52
#define PRESET_BPH { 12000, 14400, 18000, 19800, 21600, 25200, 28800, 36000, 0 };
/* algo.c */
struct filter {
double a0,a1,a2,b1,b2;
};
struct processing_buffers {
int sample_rate;
int sample_count;
float *samples, *samples_sc, *waveform, *waveform_sc;
fftwf_complex *fft;
fftwf_plan plan_a, plan_b, plan_c, plan_d;
struct filter *hpf, *lpf;
double period,sigma,be,waveform_max;
int tic,toc;
int ready;
uint64_t timestamp, last_tic, last_toc;
};
void setup_buffers(struct processing_buffers *b);
struct processing_buffers *pb_clone(struct processing_buffers *p);
void process(struct processing_buffers *p, int bph);
/* audio.c */
int start_portaudio();
int analyze_pa_data(struct processing_buffers *p, int bph);
/* interface.c */
void debug(char *format,...);
void error(char *format,...);