yin.c

Audacity是一款用於錄音和編輯聲音的、免費的開放源碼軟體。它可以執行於Mac OS X、Microsoft Windows、GNU/Linux和其它作業系統

/* yin.c -- partial implementation of the YIN algorithm, with some 
 *   fixes by DM. This code should be replaced with the fall 2002
 *   intro to computer music implementation project.
 */

#include "stdio.h"
#ifdef UNIX
#include "sys/file.h"
#endif
#ifndef mips
#include "stdlib.h"
#endif
#include "snd.h"
#include "xlisp.h"
#include "sound.h"
#include "falloc.h"
#include "yin.h"

void yin_free();

/* for multiple channel results, one susp is shared by all sounds */
/* the susp in turn must point back to all sound list tails */

typedef struct yin_susp_struct {
    snd_susp_node susp;
    long terminate_cnt;
    boolean logically_stopped;
    sound_type s;
    long s_cnt;
    sample_block_values_type s_ptr;
    long blocksize;
    long stepsize;
    sample_type *block;
    float *temp;
    sample_type *fillptr;
    sample_type *endptr;
    snd_list_type chan[2];	/* array of back pointers */
    long cnt;	/* how many sample frames to read */
    long m;
    long middle;
} yin_susp_node, *yin_susp_type;


// Uses cubic interpolation to return the value of x such
// that the function defined by f(0), f(1), f(2), and f(3)
// is maximized.
//
float CubicMaximize(float y0, float y1, float y2, float y3)
{
  // Find coefficients of cubic

  float a, b, c, d;
  float da, db, dc;
  float discriminant;
  float x1, x2;
  float dda, ddb;
  
  a = (float) (y0/-6.0 + y1/2.0 - y2/2.0 + y3/6.0);
  b = (float) (y0 - 5.0*y1/2.0 + 2.0*y2 - y3/2.0);
  c = (float) (-11.0*y0/6.0 + 3.0*y1 - 3.0*y2/2.0 + y3/3.0);
  d = y0;

  // Take derivative

  da = 3*a;
  db = 2*b;
  dc = c;

  // Find zeroes of derivative using quadratic equation
  
  discriminant = db*db - 4*da*dc;
  if (discriminant < 0.0)
    return -1.0; // error
  
  x1 = (float) ((-db + sqrt(discriminant)) / (2 * da));
  x2 = (float) ((-db - sqrt(discriminant)) / (2 * da));
  
  // The one which corresponds to a local _maximum_ in the
  // cubic is the one we want - the one with a negative
  // second derivative
  
  dda = 2*da;
  ddb = db;
  
  if (dda*x1 + ddb < 0)
    return x1;
  else
    return x2;
}


void yin_compute(yin_susp_type susp, float *pitch, float *harmonicity)
{
    float *samples = susp->block;
    int middle = susp->middle;
    /* int n = middle * 2; */
    int m = susp->m;
    float threshold = 0.9F;
    float *results = susp->temp;

    /* samples is a buffer of samples */
    /* n is the number of samples, equals twice longest period, must be even */
    /* m is the shortest period in samples */
    /* results is an array of size n/2 - m + 1, the number of different lags */

    /* work from the middle of the buffer: */
    int i, j; /* loop counters */
    /* how many different lags do we compute? */
    /* int iterations = middle + 1 - m; */
    float left_energy = 0;
    float right_energy = 0;
    /* for each window, we keep the energy so we can compute the next one */
    /* incrementally. First, we need to compute the energies for lag m-1: */
    *pitch = 0;
    for (i = 0; i < m - 1; i++) {
        float left = samples[middle - 1 - i];
        float right = samples[middle + i];
        left_energy += left * left;
        right_energy += right * right;
    }
    for (i = m; i <= middle; i++) {
        /* i is the lag and the length of the window */
        /* compute the energy for left and right */
        float left, right, energy, a;
        float harmonic;
        left = samples[middle - i];
        left_energy += left * left;
        right = samples[middle - 1 + i];
        right_energy += right * right;
        /*  compute the autocorrelation */
        a = 0;
        for (j = 0; j < i; j++) {
            a += samples[middle - i + j] * samples[middle + j];
        }
        energy = left_energy + right_energy;
        harmonic = (2 * a) / energy;
        results[i - m] = harmonic;
    }
    for (i = m; i <= middle; i++) {
        if (results[i - m] > threshold) {
            float f_i = (i - 1) + 
                CubicMaximize(results[i - m - 1], results[i - m],
                              results[i - m + 1], results[i - m + 2]);
            if (f_i < i - m - 1 || f_i > i - m + 2) f_i = (float) i;
            *pitch = (float) hz_to_step((float) susp->susp.sr / f_i);
            *harmonicity = results[i - m];
            break;
        }
    }
}


/* yin_fetch - compute F0 and harmonicity using YIN approach.  */
/*
 * The pitch (F0) is determined by finding two periods whose
 * inner product accounts for almost all of the energy. Let X and Y
 * be adjacent vectors of length N in the sample stream. Then, 
 *    if 2X*Y > threshold * (X*X + Y*Y)
 *    then the period is given by N
 * In the algorithm, we compute different sizes until we find a
 * peak above threshold. Then, we use cubic interpolation to get
 * a precise value. If no peak above threshold is found, we return
 * the first peak. The second channel returns the value 2X*Y/(X*X+Y*Y)
 * which is refered to as the "harmonicity" -- the amount of energy
 * accounted for by periodicity.
 *
 * Low sample rates are advised because of the high cost of computing
 * inner products (fast autocorrelation is not used).
 *
 * The result is a 2-channel signal running at the requested rate.
 * The first channel is the estimated pitch, and the second channel
 * is the harmonicity.
 *
 * This code is adopted from multiread, currently the only other
 * multichannel suspension in Nyquist. Comments from multiread include:
 * The susp is shared by all channels.  The susp has backpointers
 * to the tail-most snd_list node of each channel, and it is by
 * extending the list at these nodes that sounds are read in.
 * To avoid a circularity, the reference counts on snd_list nodes
 * do not include the backpointers from this susp.  When a snd_list
 * node refcount goes to zero, the yin susp's free routine
 * is called.  This must scan the backpointers to find the node that
 * has a zero refcount (the free routine is called before the node
 * is deallocated, so this is safe).  The backpointer is then set
 * to NULL.  When all backpointers are NULL, the susp itself is
 * deallocated, because it can only be referenced through the
 * snd_list nodes to which there are backpointers.
 */
void yin_fetch(yin_susp_type susp, snd_list_type snd_list)
{
    int cnt = 0; /* how many samples computed */
    int togo = 0;
    int n;
    sample_block_type f0;
    sample_block_values_type f0_ptr = NULL;
    sample_block_type harmonicity;
    sample_block_values_type harmonicity_ptr = NULL;

    register sample_block_values_type s_ptr_reg;
    register sample_type *fillptr_reg;
    register sample_type *endptr_reg = susp->endptr;

    if (susp->chan[0]) {
        falloc_sample_block(f0, "yin_fetch");
        f0_ptr = f0->samples;
        /* Since susp->chan[i] exists, we want to append a block of samples.
         * The block, out, has been allocated.  Before we insert the block,
         * we must figure out whether to insert a new snd_list_type node for
         * the block.  Recall that before SND_get_next is called, the last
         * snd_list_type in the list will have a null block pointer, and the
         * snd_list_type's susp field points to the suspension (in this case,
         * susp).  When SND_get_next (in sound.c) is called, it appends a new
         * snd_list_type and points the previous one to internal_zero_block 
         * before calling this fetch routine.  On the other hand, since 
         * SND_get_next is only going to be called on one of the channels, the
         * other channels will not have had a snd_list_type appended.
         * SND_get_next does not tell us directly which channel it wants (it
         * doesn't know), but we can test by looking for a non-null block in the
         * snd_list_type pointed to by our back-pointers in susp->chan[].  If
         * the block is null, the channel was untouched by SND_get_next, and
         * we should append a snd_list_type.  If it is non-null, then it
         * points to internal_zero_block (the block inserted by SND_get_next)
         * and a new snd_list_type has already been appended.
         */
        /* Before proceeding, it may be that garbage collection ran when we
         * allocated out, so check again to see if susp->chan[j] is Null:
         */
        if (!susp->chan[0]) {
            ffree_sample_block(f0, "yin_fetch");
            f0 = NULL; /* make sure we don't free it again */
            f0_ptr = NULL; /* make sure we don't output f0 samples */
        } else if (!susp->chan[0]->block) {
            snd_list_type snd_list = snd_list_create((snd_susp_type) susp);
            /* Now we have a snd_list to append to the channel, but a very