/* This code is placed in the public domain. */ #include <stdlib.h>#include <math.h>#include <string.h>#include "soundpipe.h" #ifndef M_PI#define M_PI 3.14159265358979323846#endif #include "voc.h" #ifndef MIN#define MIN(A,B) ((A) < (B) ? (A) : (B))#endif #ifndef MAX#define MAX(A,B) ((A) > (B) ? (A) : (B))#endif #define EPSILON 1.0e-38 #define MAX_TRANSIENTS 4   typedef struct {    SPFLOAT  freq;    SPFLOAT  tenseness;    SPFLOAT  Rd;    SPFLOAT  waveform_length;    SPFLOAT  time_in_waveform;     SPFLOAT  alpha;    SPFLOAT  E0;    SPFLOAT  epsilon;    SPFLOAT  shift;    SPFLOAT  delta;    SPFLOAT  Te;    SPFLOAT  omega;     SPFLOAT  T;} glottis;   typedef struct transient {    int  position;    SPFLOAT  time_alive;    SPFLOAT  lifetime;    SPFLOAT  strength;    SPFLOAT  exponent;    char is_free;    unsigned int id;    struct transient *next;} transient;  typedef struct {    transient pool[MAX_TRANSIENTS];    transient *root;    int size;    int next_free;} transient_pool;   typedef struct {    int n;     SPFLOAT  diameter[44];    SPFLOAT  rest_diameter[44];    SPFLOAT  target_diameter[44];    SPFLOAT  new_diameter[44];    SPFLOAT  R[44];    SPFLOAT  L[44];    SPFLOAT  reflection[45];    SPFLOAT  new_reflection[45];    SPFLOAT  junction_outL[45];    SPFLOAT  junction_outR[45];    SPFLOAT  A[44];     int nose_length;      int nose_start;     int tip_start;    SPFLOAT  noseL[28];    SPFLOAT  noseR[28];    SPFLOAT  nose_junc_outL[29];    SPFLOAT  nose_junc_outR[29];    SPFLOAT  nose_reflection[29];    SPFLOAT  nose_diameter[28];    SPFLOAT  noseA[28];     SPFLOAT  reflection_left;    SPFLOAT  reflection_right;    SPFLOAT  reflection_nose;     SPFLOAT  new_reflection_left;    SPFLOAT  new_reflection_right;    SPFLOAT  new_reflection_nose;     SPFLOAT  velum_target;     SPFLOAT  glottal_reflection;    SPFLOAT  lip_reflection;    int  last_obstruction;    SPFLOAT  fade;    SPFLOAT  movement_speed;    SPFLOAT  lip_output;    SPFLOAT  nose_output;    SPFLOAT  block_time;     transient_pool tpool;    SPFLOAT  T;} tract;  struct sp_voc {    glottis  glot; /*The Glottis*/    tract  tr; /*The Vocal Tract */    SPFLOAT  buf[512];    int counter;};     static void glottis_setup_waveform(glottis *glot, SPFLOAT lambda){    USED(lambda);     SPFLOAT Rd;    SPFLOAT Ra;    SPFLOAT Rk;    SPFLOAT Rg;     SPFLOAT Ta;    SPFLOAT Tp;    SPFLOAT Te;     SPFLOAT epsilon;    SPFLOAT shift;    SPFLOAT delta;    SPFLOAT rhs_integral;     SPFLOAT lower_integral;    SPFLOAT upper_integral;     SPFLOAT omega;    SPFLOAT s;    SPFLOAT y;    SPFLOAT z;     SPFLOAT alpha;    SPFLOAT E0;      glot->Rd = 3 * (1 - glot->tenseness);    glot->waveform_length = 1.0 / glot->freq;     Rd = glot->Rd;    if(Rd < 0.5) Rd = 0.5;    if(Rd > 2.7) Rd = 2.7;      Ra = -0.01 + 0.048*Rd;    Rk = 0.224 + 0.118*Rd;    Rg = (Rk/4)*(0.5 + 1.2*Rk)/(0.11*Rd-Ra*(0.5+1.2*Rk));      Ta = Ra;    Tp = (SPFLOAT)1.0 / (2*Rg);    Te = Tp + Tp*Rk;       epsilon = (SPFLOAT)1.0 / Ta;    shift = exp(-epsilon * (1 - Te));    delta = 1 - shift;       rhs_integral = (SPFLOAT)(1.0/epsilon) * (shift-1) + (1-Te)*shift;    rhs_integral = rhs_integral / delta;    lower_integral = - (Te - Tp) / 2 + rhs_integral;    upper_integral = -lower_integral;      omega = M_PI / Tp;    s = sin(omega * Te);     y = -M_PI * s * upper_integral / (Tp*2);    z = log(y);    alpha = z / (Tp/2 - Te);    E0 = -1 / (s * exp(alpha*Te));       glot->alpha = alpha;    glot->E0 = E0;    glot->epsilon = epsilon;    glot->shift = shift;    glot->delta = delta;    glot->Te = Te;    glot->omega = omega;}  static void glottis_init(glottis *glot, SPFLOAT sr){    glot->freq = 140; /* 140Hz frequency by default */    glot->tenseness = 0.6; /* value between 0 and 1 */    glot->T = 1.0/sr; /* big T */    glot->time_in_waveform = 0;    glottis_setup_waveform(glot, 0);}  static SPFLOAT glottis_compute(sp_data *sp, glottis *glot, SPFLOAT lambda){    SPFLOAT out;    SPFLOAT aspiration;    SPFLOAT noise;    SPFLOAT t;    SPFLOAT intensity;     intensity = 1.0;    glot->time_in_waveform += glot->T;     if(glot->time_in_waveform > glot->waveform_length) {        glot->time_in_waveform -= glot->waveform_length;        glottis_setup_waveform(glot, lambda);     }     t = (glot->time_in_waveform / glot->waveform_length);     if(t > glot->Te) {        out = (-exp(-glot->epsilon * (t-glot->Te)) + glot->shift) / glot->delta;    } else {        out = glot->E0 * exp(glot->alpha * t) * sin(glot->omega * t);    }        noise = 2.0 * ((SPFLOAT) sp_rand(sp) / SP_RANDMAX) - 1;      aspiration = intensity * (1 - sqrt(glot->tenseness)) * 0.3 * noise;     aspiration *= 0.2;     out += aspiration;     return out;}    static void tract_calculate_reflections(tract *tr){    int i;    SPFLOAT  sum;     for(i = 0; i < tr->n; i++) {        tr->A[i] = tr->diameter[i] * tr->diameter[i];        /* Calculate area from diameter squared*/    }     for(i = 1; i < tr->n; i++) {        tr->reflection[i] = tr->new_reflection[i];        if(tr->A[i] == 0) {            tr->new_reflection[i] = 0.999; /* to prevent bad behavior if 0 */        } else {            tr->new_reflection[i] =                (tr->A[i - 1] - tr->A[i]) / (tr->A[i - 1] + tr->A[i]);        }    }     tr->reflection_left = tr->new_reflection_left;    tr->reflection_right = tr->new_reflection_right;    tr->reflection_nose = tr->new_reflection_nose;     sum = tr->A[tr->nose_start] + tr->A[tr->nose_start + 1] + tr->noseA[0];    tr->new_reflection_left = (SPFLOAT)(2 * tr->A[tr->nose_start] - sum) / sum;    tr->new_reflection_right = (SPFLOAT)(2 * tr->A[tr->nose_start + 1] - sum) / sum;    tr->new_reflection_nose = (SPFLOAT)(2 * tr->noseA[0] - sum) / sum;}  static void tract_calculate_nose_reflections(tract *tr){    int i;     for(i = 0; i < tr->nose_length; i++) {        tr->noseA[i] = tr->nose_diameter[i] * tr->nose_diameter[i];    }     for(i = 1; i < tr->nose_length; i++) {        tr->nose_reflection[i] = (tr->noseA[i - 1] - tr->noseA[i]) /            (tr->noseA[i-1] + tr->noseA[i]);    }}   static int append_transient(transient_pool *pool, int position){    int i;    int free_id;    transient *t;     free_id = pool->next_free;    if(pool->size == MAX_TRANSIENTS) return 0;     if(free_id == -1) {        for(i = 0; i < MAX_TRANSIENTS; i++) {            if(pool->pool[i].is_free) {                free_id = i;                break;            }        }    }     if(free_id == -1) return 0;     t = &pool->pool[free_id];    t->next = pool->root;    pool->root = t;    pool->size++;    t->is_free = 0;    t->time_alive = 0;    t->lifetime = 0.2;    t->strength = 0.3;    t->exponent = 200;    t->position = position;    pool->next_free = -1;    return 0;}  static void remove_transient(transient_pool *pool, unsigned int id){    int i;    transient *n;     pool->next_free = id;    n = pool->root;    if(id == n->id) {        pool->root = n->next;        pool->size--;        return;    }     for(i = 0; i < pool->size; i++) {        if(n->next->id == id) {            pool->size--;            n->next->is_free = 1;            n->next = n->next->next;            break;        }        n = n->next;    }}   static SPFLOAT move_towards(SPFLOAT current, SPFLOAT target,        SPFLOAT amt_up, SPFLOAT amt_down){    SPFLOAT tmp;    if(current < target) {        tmp = current + amt_up;        return MIN(tmp, target);    }    tmp = current - amt_down;    return MAX(tmp, target);} static void tract_reshape(tract *tr){    SPFLOAT amount;    SPFLOAT slow_return;    SPFLOAT diameter;    SPFLOAT target_diameter;    int i;    int current_obstruction;     current_obstruction = -1;    amount = tr->block_time * tr->movement_speed;     for(i = 0; i < tr->n; i++) {        diameter = tr->diameter[i];        target_diameter = tr->target_diameter[i];         if(diameter < 0.001) current_obstruction = i;         if(i < tr->nose_start) slow_return = 0.6;        else if(i >= tr->tip_start) slow_return = 1.0;        else {            slow_return =                0.6+0.4*(i - tr->nose_start)/(tr->tip_start - tr->nose_start);        }         tr->diameter[i] = move_towards(diameter, target_diameter,                slow_return * amount, 2 * amount);     }     if(tr->last_obstruction > -1 && current_obstruction == -1 &&            tr->noseA[0] < 0.05) {        append_transient(&tr->tpool, tr->last_obstruction);    }    tr->last_obstruction = current_obstruction;     tr->nose_diameter[0] = move_towards(tr->nose_diameter[0], tr->velum_target,            amount * 0.25, amount * 0.1);    tr->noseA[0] = tr->nose_diameter[0] * tr->nose_diameter[0];}  static void tract_init(sp_data *sp, tract *tr){    int i;    SPFLOAT diameter, d; /* needed to set up diameter arrays */     tr->n = 44;    tr->nose_length = 28;    tr->nose_start = 17;     tr->reflection_left = 0.0;    tr->reflection_right = 0.0;    tr->reflection_nose = 0.0;    tr->new_reflection_left = 0.0;    tr->new_reflection_right= 0.0;    tr->new_reflection_nose = 0.0;    tr->velum_target = 0.01;    tr->glottal_reflection = 0.75;    tr->lip_reflection = -0.85;    tr->last_obstruction = -1;    tr->movement_speed = 15;    tr->lip_output = 0;    tr->nose_output = 0;    tr->tip_start = 32;      memset(tr->diameter, 0, tr->n * sizeof(SPFLOAT));    memset(tr->rest_diameter, 0, tr->n * sizeof(SPFLOAT));    memset(tr->target_diameter, 0, tr->n * sizeof(SPFLOAT));    memset(tr->new_diameter, 0, tr->n * sizeof(SPFLOAT));    memset(tr->L, 0, tr->n * sizeof(SPFLOAT));    memset(tr->R, 0, tr->n * sizeof(SPFLOAT));    memset(tr->reflection, 0, (tr->n + 1) * sizeof(SPFLOAT));    memset(tr->new_reflection, 0, (tr->n + 1) * sizeof(SPFLOAT));    memset(tr->junction_outL, 0, (tr->n + 1) * sizeof(SPFLOAT));    memset(tr->junction_outR, 0, (tr->n + 1) * sizeof(SPFLOAT));    memset(tr->A, 0, tr->n * sizeof(SPFLOAT));    memset(tr->noseL, 0, tr->nose_length * sizeof(SPFLOAT));    memset(tr->noseR, 0, tr->nose_length * sizeof(SPFLOAT));    memset(tr->nose_junc_outL, 0, (tr->nose_length + 1) * sizeof(SPFLOAT));    memset(tr->nose_junc_outR, 0, (tr->nose_length + 1) * sizeof(SPFLOAT));    memset(tr->nose_diameter, 0, tr->nose_length * sizeof(SPFLOAT));    memset(tr->noseA, 0, tr->nose_length * sizeof(SPFLOAT));      for(i = 0; i < tr->n; i++) {        if(i < 7 * (SPFLOAT)tr->n / 44 - 0.5) {            diameter = 0.6;        } else if( i < 12 * (SPFLOAT)tr->n / 44) {            diameter = 1.1;        } else {            diameter = 1.5;        }         tr->diameter[i] =            tr->rest_diameter[i] =            tr->target_diameter[i] =            tr->new_diameter[i] = diameter;     }          for(i = 0; i < tr->nose_length; i++) {            d = 2 * ((SPFLOAT)i / tr->nose_length);            if(d < 1) {                diameter = 0.4 + 1.6 * d;            } else {                diameter = 0.5 + 1.5*(2-d);            }            diameter = MIN(diameter, 1.9);            tr->nose_diameter[i] = diameter;        } 

Lines 1–500 of 788. Load more