/* Copyright (c) 2005 Miller Puckette. BSD licensed. No warranties. */
/*
fix parameter settings
not to report pitch if evidence too scanty?
note-on detection triggered by falling envelope (a posteriori)
reentrancy bug setting loud flag (other parameters too?)
tweaked freqs still not stable enough
implement block ("-b") mode
*/
#include "m_pd.h"
#include <math.h>
#include <stdio.h>
#include <string.h>
#ifdef MSW
#include <malloc.h>
#else
#include <alloca.h>
#endif
#include <stdlib.h>
#ifdef NT
#pragma warning( disable : 4244 )
#pragma warning( disable : 4305 )
#endif
typedef struct peak
{
float p_freq;
float p_amp;
float p_ampreal;
float p_ampimag;
float p_pit;
float p_db;
float p_salience;
float p_tmp;
} t_peak;
/********************** service routines **************************/
static int sigmund_ilog2(int n)
{
int ret = -1;
while (n)
{
n >>= 1;
ret++;
}
return (ret);
}
/* parameters for von Hann window (change these to get Hamming if desired) */
#define W_ALPHA 0.5
#define W_BETA 0.5
#define NEGBINS 4 /* number of bins of negative frequency we'll need */
#define PI 3.14159265
#define LOG2 0.69314718
#define LOG10 2.30258509
static float sinx(float theta, float sintheta)
{
if (theta > -0.003 && theta < 0.003)
return (1);
else return (sintheta/theta);
}
static float window_hann_mag(float pidetune, float sinpidetune)
{
return (W_ALPHA * sinx(pidetune, sinpidetune)
- 0.5 * W_BETA *
(sinx(pidetune+PI, sinpidetune) + sinx(pidetune-PI, sinpidetune)));
}
static float window_mag(float pidetune, float cospidetune)
{
return (sinx(pidetune + (PI/2), cospidetune)
+ sinx(pidetune - (PI/2), -cospidetune));
}
/*********** Routines to analyze a window into sinusoidal peaks *************/
static int sigmund_cmp_freq(const void *p1, const void *p2)
{
if ((*(t_peak **)p1)->p_freq > (*(t_peak **)p2)->p_freq)
return (1);
else if ((*(t_peak **)p1)->p_freq < (*(t_peak **)p2)->p_freq)
return (-1);
else return (0);
}
static void sigmund_tweak(int npts, float *ftreal, float *ftimag,
int npeak, t_peak *peaks, float fperbin, int loud)
{
t_peak **peakptrs = (t_peak **)alloca(sizeof (*peakptrs) * (npeak+1));
t_peak negpeak;
int peaki, j, k;
float ampreal[3], ampimag[3];
float binperf = 1./fperbin;
float phaseperbin = (npts-0.5)/npts, oneovern = 1./npts;
if (npeak < 1)
return;
for (peaki = 0; peaki < npeak; peaki++)
peakptrs[peaki+1] = &peaks[peaki];
qsort(peakptrs+1, npeak, sizeof (*peakptrs), sigmund_cmp_freq);
peakptrs[0] = &negpeak;
negpeak.p_ampreal = peakptrs[1]->p_ampreal;
negpeak.p_ampimag = -peakptrs[1]->p_ampimag;
negpeak.p_freq = -peakptrs[1]->p_freq;
for (peaki = 1; peaki <= npeak; peaki++)
{
int cbin = peakptrs[peaki]->p_freq*binperf + 0.5;
int nsub = (peaki == npeak ? 1:2);
float windreal, windimag, windpower, detune, pidetune, sinpidetune,
cospidetune, ampcorrect, ampout, ampoutreal, ampoutimag, freqout;
/* post("3 nsub %d amp %f freq %f", nsub,
peakptrs[peaki]->p_amp, peakptrs[peaki]->p_freq); */
if (cbin < 0 || cbin > 2*npts - 3)
continue;
for (j = 0; j < 3; j++)
ampreal[j] = ftreal[cbin+2*j-2], ampimag[j] = ftimag[cbin+2*j-2];
/* post("a %f %f", ampreal[1], ampimag[1]); */
for (j = 0; j < nsub; j++)
{
t_peak *neighbor = peakptrs[(peaki-1) + 2*j];
float neighborreal = npts * neighbor->p_ampreal;
float neighborimag = npts * neighbor->p_ampimag;
for (k = 0; k < 3; k++)
{
float freqdiff = (0.5*PI) * ((cbin + 2*k-2)
-binperf * neighbor->p_freq);
float sx = sinx(freqdiff, sin(freqdiff));
float phasere = cos(freqdiff * phaseperbin);
float phaseim = sin(freqdiff * phaseperbin);
ampreal[k] -=
sx * (phasere * neighborreal - phaseim * neighborimag);
ampimag[k] -=
sx * (phaseim * neighborreal + phasere * neighborimag);
}
/* post("b %f %f", ampreal[1], ampimag[1]); */
}
windreal = W_ALPHA * ampreal[1] -
(0.5 * W_BETA) * (ampreal[0] + ampreal[2]);
windimag = W_ALPHA * ampimag[1] -
(0.5 * W_BETA) * (ampimag[0] + ampimag[2]);
windpower = windreal * windreal + windimag * windimag;
detune = (
W_BETA*(ampreal[0] - ampreal[2]) *
(2.0*W_ALPHA * ampreal[1] - W_BETA * (ampreal[0] + ampreal[2]))
+
W_BETA*(ampimag[0] - ampimag[2]) *
(2.0*W_ALPHA * ampimag[1] - W_BETA * (ampimag[0] + ampimag[2]))
) / (4.0 * windpower);
if (detune > 0.5)
detune = 0.5;
else if (detune < -0.5)
detune = -0.5;
/* if (loud > 0)
post("tweak: windpower %f, bin %d, detune %f",
windpower, cbin, detune); */
pidetune = PI * detune;
sinpidetune = sin(pidetune);
cospidetune = cos(pidetune);
ampcorrect = 1.0 / window_hann_mag(pidetune, sinpidetune);
ampout = oneovern * ampcorrect *sqrt(windpower);
ampoutreal = oneovern * ampcorrect *
(windreal * cospidetune - windimag * sinpidetune);
ampoutimag = oneovern * ampcorrect *
(windreal * sinpidetune + windimag * cospidetune);
freqout = (cbin + 2*detune) * fperbin;
if (loud > 1)
post("amp %f, freq %f", ampout, freqout);
peakptrs[peaki]->p_freq = freqout;
peakptrs[peaki]->p_amp = ampout;
peakptrs[peaki]->p_ampreal = ampoutreal;
peakptrs[peaki]->p_ampimag = ampoutimag;
}
}
static void sigmund_getrawpeaks(int npts, float *insamps,
int npeak, t_peak *peakv, int *nfound, float *power, float srate, int loud,
float param1, float param2, float param3, float hifreq)
{
float oneovern = 1.0/ (float)npts;
float fperbin = 0.5 * srate * oneovern;
int npts2 = 2*npts, i, bin;
int count, peakcount = 0;
float *fp1, *fp2;
float *rawpow, *rawreal, *rawimag, *maskbuf;
float *bigbuf = alloca(sizeof (float ) * (2*NEGBINS + 5*npts));
int maxbin = hifreq/fperbin;
int altwind = (param3 == 1);
int tweak = (param3 == 0);
if (maxbin > npts - NEGBINS)
maxbin = npts - NEGBINS;
if (loud) post("tweak %d", tweak);
maskbuf = bigbuf + npts2;
rawreal = maskbuf + npts+NEGBINS;
rawimag = rawreal+npts+NEGBINS;
for (i = 0; i < npts; i++)
maskbuf[i] = 0;
for (i = 0; i < npts; i++)
bigbuf[i] = insamps[i];
for (i = npts; i < 2*npts; i++)
bigbuf[i] = 0;
mayer_realfft(npts2, bigbuf);
for (i = 0; i < npts; i++)
rawreal[i] = bigbuf[i];
for (i = 1; i < npts-1; i++)
rawimag[i] = bigbuf[npts2-i];
if (loud && npts == 1024)
{
float bigbuf2[2048];
for (i = 0; i < 1024; i++)
bigbuf2[i] = insamps[i];
for (i = 1024; i < 2048; i++)
bigbuf2[i] = 0;
mayer_realfft(2048, bigbuf2);
for (i = 1; i < 10; i++)
post("(%10.2f, %10.2f) -> (%10.2f, %10.2f)",
bigbuf2[i], bigbuf2[2048-i], rawreal[i], rawimag[i]);
}
rawreal[-1] = rawreal[1];
rawreal[-2] = rawreal[2];
rawreal[-3] = rawreal[3];
rawreal[-4] = rawreal[4];
rawimag[0] = rawimag[npts-1] = 0;
rawimag[-1] = -rawimag[1];
rawimag[-2] = -rawimag[2];
rawimag[-3] = -rawimag[3];
rawimag[-4] = -rawimag[4];
for (peakcount = 0; peakcount < npeak; peakcount++)
{
float pow1, maxpower = 0, totalpower = 0, windreal, windimag, windpower,
detune, pidetune, sinpidetune, cospidetune, ampcorrect, ampout,
ampoutreal, ampoutimag, freqout, freqcount1, freqcount2, powmask;
int bestindex = -1;
for (bin = 2, fp1 = rawreal+2, fp2 = rawimag+2;
bin < maxbin; bin++, fp1++, fp2++)
{
float x1, x2, a1, a2, b1, b2, thresh;
if (altwind)
{
x1 = fp1[0] - 0.5*(fp1[-2] +fp1[2]);
x2 = fp2[0] - 0.5*(fp2[-2] +fp2[2]);
a1 = fp1[4] - 0.5*(fp1[2] +fp1[6]);
a2 = fp2[2] - 0.5*(fp2[2] +fp2[6]);
b1 = fp1[-4] - 0.5*(fp1[-2] +fp1[-6]);
b2 = fp2[-4] - 0.5*(fp2[-2] +fp2[-6]);
thresh = param2 * (a1*a1+a2*a2+b1*b1+b2*b2);
}
else
{
x1 = fp1[1] - fp1[-1];
x2 = fp2[1] - fp2[-1];
a1 = fp1[3] - fp1[1];
a2 = fp2[3] - fp2[1];
b1 = fp1[-1] - fp1[-3];
b2 = fp2[-1] - fp2[-3];
thresh = param2 * (a1*a1+a2*a2+b1*b1+b2*b2);
}
pow1 = x1*x1+x2*x2;
if (pow1 > maxpower && pow1 > maskbuf[bin])
{
if (pow1 > thresh)
maxpower = pow1, bestindex = bin;
}
totalpower += pow1;
}
if (totalpower <= 0 || maxpower < 1e-10*totalpower || bestindex < 0)
break;
fp1 = rawreal+bestindex;
fp2 = rawimag+bestindex;
*power = 0.5 * totalpower *oneovern * oneovern;
powmask = maxpower * exp(-param1 * log(10.) / 10.);
if (loud > 2)
post("maxpower %f, powmask %f, param1 %f",
maxpower, powmask, param1);
for (bin = 2; bin < maxbin; bin++)
{
float bindiff = bin - bestindex;
float mymask =
powmask/ (1. + bindiff * bindiff * bindiff * bindiff);
if (bindiff < 2 && bindiff > -2)
mymask = 2*maxpower;
if (mymask > maskbuf[bin])
maskbuf[bin] = mymask;
}
if (loud > 1)
post("best index %d, total power %f", bestindex, totalpower);
if (altwind)
{
windreal = W_ALPHA * fp1[0] - (0.5 * W_BETA) * (fp1[2] + fp1[-2]);
windimag = W_ALPHA * fp2[0] - (0.5 * W_BETA) * (fp2[2] + fp2[-2]);
windpower = windreal * windreal + windimag * windimag;
detune = (
(W_BETA*(rawreal[bestindex-2] - rawreal[bestindex+2])) *
(2.0 * W_ALPHA * rawreal[bestindex] -
W_BETA * (rawreal[bestindex-2] + rawreal[bestindex+2]))
+
(W_BETA*(rawimag[bestindex-2] - rawimag[bestindex+2])) *
(2.0 * W_ALPHA * rawimag[bestindex] -
W_BETA * (rawimag[bestindex-2] + rawimag[bestindex+2]))
) / (4.0 * windpower);
}
else
{
windreal = fp1[1] - fp1[-1];
windimag = fp2[1] - fp2[-1];
windpower = windreal * windreal + windimag * windimag;
detune = ((fp1[1] * fp1[1] - fp1[-1]*fp1[-1])
+ (fp2[1] * fp2[1] - fp2[-1]*fp2[-1])) / (2 * windpower);
if (loud > 2) post("(-1) %f %f; (1) %f %f",
fp1[-1], fp2[-1], fp1[1], fp2[1]);
if (loud > 2) post("peak %f %f",
fp1[0], fp2[0]);
}
if (detune > 0.5)
detune = 0.5;
else if (detune < -0.5)
detune = -0.5;
if (loud > 1)
post("windpower %f, index %d, detune %f",
windpower, bestindex, detune);
pidetune = PI * detune;
sinpidetune = sin(pidetune);
cospidetune = cos(pidetune);
if (altwind)
ampcorrect = 1.0 / window_hann_mag(pidetune, sinpidetune);
else ampcorrect = 1.0 / window_mag(pidetune, cospidetune);
ampout = ampcorrect *sqrt(windpower);
ampoutreal = ampcorrect *
(windreal * cospidetune - windimag * sinpidetune);
ampoutimag = ampcorrect *
(windreal * sinpidetune + windimag * cospidetune);
/* the frequency is the sum of the bin frequency and detuning */
peakv[peakcount].p_freq = (freqout = (bestindex + 2*detune)) * fperbin;
peakv[peakcount].p_amp = oneovern * ampout;
peakv[peakcount].p_ampreal = oneovern * ampoutreal;
peakv[peakcount].p_ampimag = oneovern * ampoutimag;
}
if (tweak)
{
sigmund_tweak(npts, rawreal, rawimag, peakcount, peakv, fperbin, loud);
sigmund_tweak(npts, rawreal, rawimag, peakcount, peakv, fperbin, loud);
}
for (i = 0; i < peakcount; i++)
{
peakv[i].p_pit = ftom(peakv[i].p_freq);
peakv[i].p_db = powtodb(peakv[i].p_amp);
}
*nfound = peakcount;
}
/*************** Routines for finding fundamental pitch *************/
#define PITCHNPEAK 12
#define PITCHUNCERTAINTY 0.3
#define HALFTONEINC 0.059
#define SUBHARMONICS 16
#define DBPERHALFTONE 0.5
static void sigmund_getpitch(int npeak, t_peak *peakv, float *freqp,
float npts, float srate, int loud)
{
float fperbin = 0.5 * srate / npts;
int npit = 48 * sigmund_ilog2(npts), i, j, k, nsalient;
float bestbin, bestweight, sumamp, sumweight, sumfreq, sumallamp,
freq;
float *weights = (float *)alloca(sizeof(float) * npit);
t_peak *bigpeaks[PITCHNPEAK];
int nbigpeaks;
if (npeak < 1)
{
freq = 0;
goto done;
}
for (i = 0; i < npit; i++)
weights[i] = 0;
for (i = 0; i < npeak; i++)
{
peakv[i].p_tmp = 0;
peakv[i].p_salience = peakv[i].p_db - DBPERHALFTONE * peakv[i].p_pit;
}
for (nsalient = 0; nsalient < PITCHNPEAK; nsalient++)
{
t_peak *bestpeak = 0;
float bestsalience = -1e20;
for (j = 0; j < npeak; j++)
if (peakv[j].p_tmp == 0 && peakv[j].p_salience > bestsalience)
{
bestsalience = peakv[j].p_salience;
bestpeak = &peakv[j];
}
if (!bestpeak)
break;
bigpeaks[nsalient] = bestpeak;
bestpeak->p_tmp = 1;
/* post("peak f=%f a=%f", bestpeak->p_freq, bestpeak->p_amp); */
}
sumweight = 0;
for (i = 0; i < nsalient; i++)
{
t_peak *thispeak = bigpeaks[i];
float pitchuncertainty =
4 * PITCHUNCERTAINTY * fperbin / (HALFTONEINC * thispeak->p_freq);
float weightindex = (48./LOG2) *
log(thispeak->p_freq/(2.*fperbin));
float loudness = sqrt(thispeak->p_amp);
/* post("index %f, uncertainty %f", weightindex, pitchuncertainty); */
for (j = 0; j < SUBHARMONICS; j++)
{
float subindex = weightindex -
(48./LOG2) * log(j + 1.);
int loindex = subindex - pitchuncertainty;
int hiindex = subindex + pitchuncertainty + 1;
if (hiindex < 0)
break;
if (hiindex >= npit)
continue;
if (loindex < 0)
loindex = 0;
for (k = loindex; k <= hiindex; k++)
weights[k] += loudness * 4. / (4. + j);
}
sumweight += loudness;
}
#if 0
for (i = 0; i < npit; i++)
{
postfloat(weights[i]);
if (!((i+1)%12)) post("");
}
#endif
bestbin = -1;
bestweight = -1e20;
for (i = 0; i < npit; i++)
if (weights[i] > bestweight)
bestweight = weights[i], bestbin = i;
if (bestweight < sumweight * 0.4)
bestbin = -1;
if (bestbin < 0)
{
freq = 0;
goto done;
}
for (i = bestbin+1; i < npit; i++)
{
if (weights[i] < bestweight)
break;
bestbin += 0.5;
}
freq = 2*fperbin * exp((LOG2/48.)*bestbin);
for (sumamp = sumweight = sumfreq = 0, i = 0; i < nsalient; i++)
{
t_peak *thispeak = bigpeaks[i];
float thisloudness = sqrt(thispeak->p_amp);
float thisfreq = thispeak->p_freq;
float harmonic = thisfreq/freq;
float intpart = (int)(0.5 + harmonic);
float inharm = freq * (harmonic - intpart);
if (harmonic < 1)
continue;
if (inharm < 0.25*fperbin && inharm > -0.25*fperbin)
{
float weight = thisloudness * intpart;
sumweight += weight;
sumfreq += weight*thisfreq/intpart;
}
}
if (sumweight > 0)
freq = sumfreq / sumweight;
done:
if (!(freq >= 0 || freq <= 0))
{
post("freq nan cancelled");
freq = 0;
}
*freqp = freq;
}
/*************** gather peak lists into sinusoidal tracks *************/
static void sigmund_peaktrack(int ninpeak, t_peak *inpeakv,
int noutpeak, t_peak *outpeakv, int loud)
{
int incnt, outcnt;
for (outcnt = 0; outcnt < noutpeak; outcnt++)
outpeakv[outcnt].p_tmp = -1;
/* first pass. Match each "in" peak with the closest previous
"out" peak, but no two to the same one. */
for (incnt = 0; incnt < ninpeak; incnt++)
{
float besterror = 1e20;
int bestcnt = -1;
inpeakv[incnt].p_tmp = -1;
for (outcnt = 0; outcnt < noutpeak; outcnt++)
{
float thiserror =
inpeakv[incnt].p_freq - outpeakv[outcnt].p_freq;
if (thiserror < 0)
thiserror = -thiserror;
if (thiserror < besterror)
{
besterror = thiserror;
bestcnt = outcnt;
}
}
if (outpeakv[bestcnt].p_tmp < 0)
{
outpeakv[bestcnt] = inpeakv[incnt];
inpeakv[incnt].p_tmp = 0;
outpeakv[bestcnt].p_tmp = 0;
}
}
/* second pass. Unmatched "in" peaks assigned to free "out"
peaks */
for (incnt = 0; incnt < ninpeak; incnt++)
if (inpeakv[incnt].p_tmp < 0)
{
for (outcnt = 0; outcnt < noutpeak; outcnt++)
if (outpeakv[outcnt].p_tmp < 0)
{
outpeakv[outcnt] = inpeakv[incnt];
inpeakv[incnt].p_tmp = 0;
outpeakv[outcnt].p_tmp = 1;
break;
}
}
for (outcnt = 0; outcnt < noutpeak; outcnt++)
if (outpeakv[outcnt].p_tmp == -1)
outpeakv[outcnt].p_amp = 0;
}
/**************** parse continuous pitch into note starts ***************/
#define NHISTPOINT 100
typedef struct _histpoint
{
float h_freq;
float h_power;
} t_histpoint;
typedef struct _notefinder
{
float n_age;
float n_hifreq;
float n_lofreq;
int n_peaked;
t_histpoint n_hist[NHISTPOINT];
int n_histphase;
} t_notefinder;
static void notefinder_init(t_notefinder *x)
{
int i;
x->n_peaked = x->n_age = 0;
x->n_hifreq = x->n_lofreq = 0;
x->n_histphase = 0;
for (i = 0; i < NHISTPOINT; i++)
x->n_hist[i].h_freq =x->n_hist[i].h_power = 0;
}
static void notefinder_doit(t_notefinder *x, float freq, float power,
float *note, float vibrato, int stableperiod, float powerthresh,
float growththresh, int loud)
{
/* calculate frequency ratio between allowable vibrato extremes
(equal to twice the vibrato deviation from center) */
float vibmultiple = exp((2*LOG2/12) * vibrato);
int oldhistphase, i, k;
if (stableperiod > NHISTPOINT - 1)
stableperiod = NHISTPOINT - 1;
if (++x->n_histphase == NHISTPOINT)
x->n_histphase = 0;
x->n_hist[x->n_histphase].h_freq = freq;
x->n_hist[x->n_histphase].h_power = power;
x->n_age++;
*note = 0;
if (loud)
{
post("stable %d, age %d, vibmultiple %f, powerthresh %f, hifreq %f",
stableperiod, (int)x->n_age ,vibmultiple, powerthresh, x->n_hifreq);
post("histfreq %f %f %f %f",
x->n_hist[x->n_histphase].h_freq,
x->n_hist[(x->n_histphase+NHISTPOINT-1)%NHISTPOINT].h_freq,
x->n_hist[(x->n_histphase+NHISTPOINT-2)%NHISTPOINT].h_freq,
x->n_hist[(x->n_histphase+NHISTPOINT-3)%NHISTPOINT].h_freq);
post("power %f %f %f %f",
x->n_hist[x->n_histphase].h_power,
x->n_hist[(x->n_histphase+NHISTPOINT-1)%NHISTPOINT].h_power,
x->n_hist[(x->n_histphase+NHISTPOINT-2)%NHISTPOINT].h_power,
x->n_hist[(x->n_histphase+NHISTPOINT-3)%NHISTPOINT].h_power);
for (i = 0, k = x->n_histphase; i < stableperiod; i++)
{
post("pit %5.1f pow %f", ftom(x->n_hist[k].h_freq),
x->n_hist[k].h_power);
if (--k < 0)
k = NHISTPOINT - 1;
}
}
/* look for shorter notes than "stableperiod" in length.
The amplitude must rise and then fall while the pitch holds
steady. */
if (x->n_hifreq <= 0 && x->n_age > stableperiod)
{
float maxpow = 0, freqatmaxpow = 0,
localhifreq = -1e20, locallofreq = 1e20;
int startphase = x->n_histphase - stableperiod + 1;
if (startphase < 0)
startphase += NHISTPOINT;
for (i = 0, k = startphase; i < stableperiod; i++)
{
if (x->n_hist[k].h_freq <= 0)
break;
if (x->n_hist[k].h_power > maxpow)
maxpow = x->n_hist[k].h_power,
freqatmaxpow = x->n_hist[k].h_freq;
if (x->n_hist[k].h_freq > localhifreq)
localhifreq = x->n_hist[k].h_freq;
if (x->n_hist[k].h_freq < locallofreq)
locallofreq = x->n_hist[k].h_freq;
if (localhifreq > locallofreq * vibmultiple)
break;
if (maxpow > power * growththresh &&
maxpow > x->n_hist[startphase].h_power * growththresh &&
localhifreq < vibmultiple * locallofreq
&& freqatmaxpow > 0 && maxpow > powerthresh)
{
x->n_hifreq = x->n_lofreq = *note = freqatmaxpow;
x->n_age = 0;
x->n_peaked = 0;
/* post("got short note"); */
return;
}
if (++k >= NHISTPOINT)
k = 0;
}
}
if (x->n_hifreq > 0)
{
/* test if we're within "vibrato" range, and if so update range */
if (freq * vibmultiple >= x->n_hifreq &&
x->n_lofreq * vibmultiple >= freq)
{
if (freq > x->n_hifreq)
x->n_hifreq = freq;
if (freq < x->n_lofreq)
x->n_lofreq = freq;
}
else if (x->n_hifreq > 0 && x->n_age > stableperiod)
{
/* if we've been out of range at least 1/2 the
last "stableperiod" analyses, clear the note */
int nbad = 0;
for (i = 0, k = x->n_histphase; i < stableperiod - 1; i++)
{
if (--k < 0)
k = NHISTPOINT - 1;
if (x->n_hist[k].h_freq * vibmultiple <= x->n_hifreq ||
x->n_lofreq * vibmultiple <= x->n_hist[k].h_freq)
nbad++;
}
if (2 * nbad >= stableperiod)
{
x->n_hifreq = x->n_lofreq = 0;
x->n_age = 0;
}
}
}
oldhistphase = x->n_histphase - stableperiod;
if (oldhistphase < 0)
oldhistphase += NHISTPOINT;
/* look for envelope attacks */
if (x->n_hifreq > 0 && x->n_peaked)
{
if (freq > 0 && power > powerthresh &&
power > x->n_hist[oldhistphase].h_power *
exp((LOG10*0.1)*growththresh))
{
/* clear it and fall through for new stable-note test */
x->n_peaked = 0;
x->n_hifreq = x->n_lofreq = 0;
x->n_age = 0;
}
}
else if (!x->n_peaked)
{
if (x->n_hist[oldhistphase].h_power > powerthresh &&
x->n_hist[oldhistphase].h_power > power)
x->n_peaked = 1;
}
/* test for a new note using a stability criterion. */
if (freq >= 0 &&
(x->n_hifreq <= 0 || freq > x->n_hifreq || freq < x->n_lofreq))
{
float testfhi, testflo, maxpow = 0;
for (i = 0, k = x->n_histphase, testfhi = testflo = freq;
i < stableperiod-1; i++)
{
if (--k < 0)
k = NHISTPOINT - 1;
if (x->n_hist[k].h_freq > testfhi)
testfhi = x->n_hist[k].h_freq;
if (x->n_hist[k].h_freq < testflo)
testflo = x->n_hist[k].h_freq;
if (x->n_hist[k].h_power > maxpow)
maxpow = x->n_hist[k].h_power;
}
if (testflo > 0 && testfhi <= vibmultiple * testflo
&& maxpow > powerthresh)
{
/* report new note */
float sumf = 0, sumw = 0, thisf, thisw;
for (i = 0, k = x->n_histphase; i < stableperiod; i++)
{
thisw = x->n_hist[k].h_power;
sumw += thisw;
sumf += thisw*x->n_hist[k].h_freq;
if (--k < 0)
k = NHISTPOINT - 1;
}
x->n_hifreq = x->n_lofreq = *note = (sumw > 0 ? sumf/sumw : 0);
#if 0
/* debugging printout */
for (i = 0; i < stableperiod; i++)
{
int k3 = x->n_histphase - i;
if (k3 < 0)
k3 += NHISTPOINT;
startpost("%5.1f ", ftom(x->n_hist[k3].h_freq));
}
post("");
#endif
x->n_age = 0;
x->n_peaked = 0;
return;
}
}
*note = 0;
return;
}
/*************************** Glue for Pd ************************/
static t_class *sigmund_class;
#define NHIST 100
#define MODE_STREAM 1
#define MODE_BLOCK 2 /* unimplemented */
#define MODE_TABLE 3
#define NPOINTS_DEF 1024
#define NPOINTS_MIN 128
#define HOP_DEF 512
#define NPEAK_DEF 20
#define VIBRATO_DEF 1
#define STABLETIME_DEF 50
#define MINPOWER_DEF 50
#define GROWTH_DEF 7
#define OUT_PITCH 0
#define OUT_ENV 1
#define OUT_NOTE 2
#define OUT_PEAKS 3
#define OUT_TRACKS 4
#define OUT_SMSPITCH 5
#define OUT_SMSNONPITCH 6
typedef struct _varout
{
t_outlet *v_outlet;
int v_what;
} t_varout;
typedef struct _sigmund
{
t_object x_obj;
t_varout *x_varoutv;
int x_nvarout;
t_clock *x_clock;
float x_f; /* for main signal inlet */
float x_sr; /* sample rate */
int x_mode; /* MODE_STREAM, etc. */
int x_npts; /* number of points in analysis window */
int x_npeak; /* number of peaks to find */
int x_loud; /* debug level */
t_sample *x_inbuf; /* input buffer */
int x_infill; /* number of points filled */
int x_countdown; /* countdown to start filling buffer */
int x_hop; /* samples between analyses */
float x_maxfreq; /* highest-frequency peak to report */
float x_vibrato; /* vibrato depth in half tones */
float x_stabletime; /* period of stability needed for note */
float x_growth; /* growth to set off a new note */
float x_minpower; /* minimum power, in DB, for a note */
float x_param1;
float x_param2;
float x_param3;
t_notefinder x_notefinder;
t_peak *x_trackv;
int x_ntrack;
unsigned int x_dopitch:1;
unsigned int x_donote:1;
unsigned int x_dotracks:1;
} t_sigmund;
static void sigmund_clock(t_sigmund *x);
static void sigmund_clear(t_sigmund *x);
static void sigmund_npts(t_sigmund *x, t_floatarg f);
static void sigmund_hop(t_sigmund *x, t_floatarg f);
static void sigmund_npeak(t_sigmund *x, t_floatarg f);
static void sigmund_maxfreq(t_sigmund *x, t_floatarg f);
static void sigmund_vibrato(t_sigmund *x, t_floatarg f);
static void sigmund_stabletime(t_sigmund *x, t_floatarg f);
static void sigmund_growth(t_sigmund *x, t_floatarg f);
static void sigmund_minpower(t_sigmund *x, t_floatarg f);
static void *sigmund_new(t_symbol *s, int argc, t_atom *argv)
{
t_sigmund *x = (t_sigmund *)pd_new(sigmund_class);
x->x_param1 = 0;
x->x_param2 = 0.6;
x->x_param3 = 0;
x->x_npts = NPOINTS_DEF;
x->x_hop = HOP_DEF;
x->x_mode = MODE_STREAM;
x->x_npeak = NPEAK_DEF;
x->x_vibrato = VIBRATO_DEF;
x->x_stabletime = STABLETIME_DEF;
x->x_growth = GROWTH_DEF;
x->x_minpower = MINPOWER_DEF;
x->x_maxfreq = 1000000;
x->x_loud = 0;
x->x_sr = 1;
x->x_nvarout = 0;
x->x_varoutv = (t_varout *)getbytes(0);
x->x_trackv = 0;
x->x_ntrack = 0;
x->x_dopitch = x->x_donote = x->x_dotracks = 0;
while (argc > 0)
{
t_symbol *firstarg = atom_getsymbolarg(0, argc, argv);
if (!strcmp(firstarg->s_name, "-t"))
{
x->x_mode = MODE_TABLE;
argc--, argv++;
}
else if (!strcmp(firstarg->s_name, "-s"))
{
x->x_mode = MODE_STREAM;
argc--, argv++;
}
#if 0
else if (!strcmp(firstarg->s_name, "-b"))
{
x->x_mode = MODE_BLOCK;
argc--, argv++;
}
#endif
else if (!strcmp(firstarg->s_name, "-npts") && argc > 1)
{
sigmund_npts(x, atom_getfloatarg(1, argc, argv));
argc -= 2; argv += 2;
}
else if (!strcmp(firstarg->s_name, "-hop") && argc > 1)
{
sigmund_hop(x, atom_getfloatarg(1, argc, argv));
argc -= 2; argv += 2;
}
else if (!strcmp(firstarg->s_name, "-npeak") && argc > 1)
{
sigmund_npeak(x, atom_getfloatarg(1, argc, argv));
argc -= 2; argv += 2;
}
else if (!strcmp(firstarg->s_name, "-maxfreq") && argc > 1)
{
sigmund_maxfreq(x, atom_getfloatarg(1, argc, argv));
argc -= 2; argv += 2;
}
else if (!strcmp(firstarg->s_name, "-vibrato") && argc > 1)
{
sigmund_vibrato(x, atom_getfloatarg(1, argc, argv));
argc -= 2; argv += 2;
}
else if (!strcmp(firstarg->s_name, "-stabletime") && argc > 1)
{
sigmund_stabletime(x, atom_getfloatarg(1, argc, argv));
argc -= 2; argv += 2;
}
else if (!strcmp(firstarg->s_name, "-growth") && argc > 1)
{
sigmund_growth(x, atom_getfloatarg(1, argc, argv));
argc -= 2; argv += 2;
}
else if (!strcmp(firstarg->s_name, "-minpower") && argc > 1)
{
sigmund_minpower(x, atom_getfloatarg(1, argc, argv));
argc -= 2; argv += 2;
}
else if (!strcmp(firstarg->s_name, "pitch"))
{
int n2 = x->x_nvarout+1;
x->x_varoutv = (t_varout *)t_resizebytes(x->x_varoutv,
x->x_nvarout*sizeof(t_varout), n2*sizeof(t_varout));
x->x_varoutv[x->x_nvarout].v_outlet =
outlet_new(&x->x_obj, &s_float);
x->x_varoutv[x->x_nvarout].v_what = OUT_PITCH;
x->x_nvarout = n2;
x->x_dopitch = 1;
argc--, argv++;
}
else if (!strcmp(firstarg->s_name, "env"))
{
int n2 = x->x_nvarout+1;
x->x_varoutv = (t_varout *)t_resizebytes(x->x_varoutv,
x->x_nvarout*sizeof(t_varout), n2*sizeof(t_varout));
x->x_varoutv[x->x_nvarout].v_outlet =
outlet_new(&x->x_obj, &s_float);
x->x_varoutv[x->x_nvarout].v_what = OUT_ENV;
x->x_nvarout = n2;
argc--, argv++;
}
else if (!strcmp(firstarg->s_name, "note")
|| !strcmp(firstarg->s_name, "notes"))
{
int n2 = x->x_nvarout+1;
x->x_varoutv = (t_varout *)t_resizebytes(x->x_varoutv,
x->x_nvarout*sizeof(t_varout), n2*sizeof(t_varout));
x->x_varoutv[x->x_nvarout].v_outlet =
outlet_new(&x->x_obj, &s_float);
x->x_varoutv[x->x_nvarout].v_what = OUT_NOTE;
x->x_nvarout = n2;
x->x_dopitch = x->x_donote = 1;
argc--, argv++;
}
else if (!strcmp(firstarg->s_name, "peaks"))
{
int n2 = x->x_nvarout+1;
x->x_varoutv = (t_varout *)t_resizebytes(x->x_varoutv,
x->x_nvarout*sizeof(t_varout), n2*sizeof(t_varout));
x->x_varoutv[x->x_nvarout].v_outlet =
outlet_new(&x->x_obj, &s_list);
x->x_varoutv[x->x_nvarout].v_what = OUT_PEAKS;
x->x_nvarout = n2;
argc--, argv++;
}
else if (!strcmp(firstarg->s_name, "tracks"))
{
int n2 = x->x_nvarout+1;
x->x_varoutv = (t_varout *)t_resizebytes(x->x_varoutv,
x->x_nvarout*sizeof(t_varout), n2*sizeof(t_varout));
x->x_varoutv[x->x_nvarout].v_outlet =
outlet_new(&x->x_obj, &s_list);
x->x_varoutv[x->x_nvarout].v_what = OUT_TRACKS;
x->x_nvarout = n2;
x->x_dotracks = 1;
argc--, argv++;
}
else
{
pd_error(x, "sigmund: %s: unknown flag or argument missing",
firstarg->s_name);
argc--, argv++;
}
}
if (!x->x_nvarout)
{
x->x_varoutv = (t_varout *)t_resizebytes(x->x_varoutv,
0, 2*sizeof(t_varout));
x->x_varoutv[0].v_outlet = outlet_new(&x->x_obj, &s_float);
x->x_varoutv[0].v_what = OUT_PITCH;
x->x_varoutv[1].v_outlet = outlet_new(&x->x_obj, &s_float);
x->x_varoutv[1].v_what = OUT_ENV;
x->x_nvarout = 2;
x->x_dopitch = 1;
}
if (x->x_dotracks)
{
x->x_ntrack = x->x_npeak;
x->x_trackv = (t_peak *)getbytes(x->x_ntrack * sizeof(*x->x_trackv));
}
x->x_clock = clock_new(&x->x_obj.ob_pd, (t_method)sigmund_clock);
/* check parameter ranges */
if (x->x_npts < NPOINTS_MIN)
post("sigmund~: minimum points %d", NPOINTS_MIN),
x->x_npts = NPOINTS_MIN;
if (x->x_npts != (1 << sigmund_ilog2(x->x_npts)))
post("sigmund~: adjusting analysis size to %d points",
(x->x_npts = (1 << sigmund_ilog2(x->x_npts))));
if (x->x_hop != (1 << sigmund_ilog2(x->x_hop)))
post("sigmund~: adjusting hop size to %d points",
(x->x_hop = (1 << sigmund_ilog2(x->x_hop))));
if (x->x_mode == MODE_STREAM)
x->x_inbuf = getbytes(sizeof(*x->x_inbuf) * x->x_npts);
else x->x_inbuf = 0;
x->x_infill = 0;
x->x_countdown = 0;
notefinder_init(&x->x_notefinder);
sigmund_clear(x);
return (x);
}
static void sigmund_doit(t_sigmund *x, int npts, float *arraypoints,
int loud, float srate)
{
t_peak *peakv = (t_peak *)alloca(sizeof(t_peak) * x->x_npeak);
int nfound, i, cnt;
float freq = 0, power, note = 0;
sigmund_getrawpeaks(npts, arraypoints, x->x_npeak, peakv,
&nfound, &power, srate, loud, x->x_param1, x->x_param2, x->x_param3,
x->x_maxfreq);
if (x->x_dopitch)
sigmund_getpitch(nfound, peakv, &freq, npts, srate, loud);
if (x->x_donote)
notefinder_doit(&x->x_notefinder, freq, power, ¬e, x->x_vibrato,
x->x_stabletime * 0.001f * x->x_sr / (float)x->x_hop,
exp(LOG10*0.1*(x->x_minpower - 100)), x->x_growth, loud);
if (x->x_dotracks)
sigmund_peaktrack(nfound, peakv, x->x_ntrack, x->x_trackv, loud);
for (cnt = x->x_nvarout; cnt--;)
{
t_varout *v = &x->x_varoutv[cnt];
switch (v->v_what)
{
case OUT_PITCH:
outlet_float(v->v_outlet, ftom(freq));
break;
case OUT_ENV:
outlet_float(v->v_outlet, powtodb(power));
break;
case OUT_NOTE:
if (note > 0)
outlet_float(v->v_outlet, ftom(note));
break;
case OUT_PEAKS:
for (i = 0; i < nfound; i++)
{
t_atom at[5];
SETFLOAT(at, (float)i);
SETFLOAT(at+1, peakv[i].p_freq);
SETFLOAT(at+2, 2*peakv[i].p_amp);
SETFLOAT(at+3, 2*peakv[i].p_ampreal);
SETFLOAT(at+4, 2*peakv[i].p_ampimag);
outlet_list(v->v_outlet, &s_list, 5, at);
}
break;
case OUT_TRACKS:
for (i = 0; i < x->x_ntrack; i++)
{
t_atom at[4];
SETFLOAT(at, (float)i);
SETFLOAT(at+1, x->x_trackv[i].p_freq);
SETFLOAT(at+2, 2*x->x_trackv[i].p_amp);
SETFLOAT(at+3, x->x_trackv[i].p_tmp);
outlet_list(v->v_outlet, &s_list, 4, at);
}
break;
}
}
}
static void sigmund_list(t_sigmund *x, t_symbol *s, int argc, t_atom *argv)
{
t_symbol *syminput = atom_getsymbolarg(0, argc, argv);
int npts = atom_getintarg(1, argc, argv);
int onset = atom_getintarg(2, argc, argv);
float srate = atom_getfloatarg(3, argc, argv);
int loud = atom_getfloatarg(4, argc, argv);
int arraysize, totstorage, nfound, i;
t_garray *a;
float *arraypoints, pit;
t_word *wordarray = 0;
if (argc < 5)
{
post(
"sigmund: array-name, npts, array-onset, samplerate, loud");
return;
}
if (npts < 64 || npts != (1 << ilog2(npts)))
{
error("sigmund: bad npoints");
return;
}
if (onset < 0)
{
error("sigmund: negative onset");
return;
}
arraypoints = alloca(sizeof(float)*npts);
if (!(a = (t_garray *)pd_findbyclass(syminput, garray_class)) ||
!garray_getfloatwords(a, &arraysize, &wordarray) ||
arraysize < onset + npts)
{
error("%s: array missing or too small", syminput->s_name);
return;
}
if (arraysize < npts)
{
error("sigmund~: too few points in array");
return;
}
for (i = 0; i < npts; i++)
arraypoints[i] = wordarray[i+onset].w_float;
sigmund_doit(x, npts, arraypoints, loud, srate);
}
static void sigmund_clear(t_sigmund *x)
{
if (x->x_trackv)
memset(x->x_trackv, 0, x->x_ntrack * sizeof(*x->x_trackv));
x->x_infill = x->x_countdown = 0;
}
/* these are for testing; their meanings vary... */
static void sigmund_param1(t_sigmund *x, t_floatarg f)
{
x->x_param1 = f;
}
static void sigmund_param2(t_sigmund *x, t_floatarg f)
{
x->x_param2 = f;
}
static void sigmund_param3(t_sigmund *x, t_floatarg f)
{
x->x_param3 = f;
}
static void sigmund_npts(t_sigmund *x, t_floatarg f)
{
x->x_npts = f;
/* check parameter ranges */
if (x->x_npts < NPOINTS_MIN)
post("sigmund~: minimum points %d", NPOINTS_MIN),
x->x_npts = NPOINTS_MIN;
if (x->x_npts != (1 << sigmund_ilog2(x->x_npts)))
post("sigmund~: adjusting analysis size to %d points",
(x->x_npts = (1 << sigmund_ilog2(x->x_npts))));
}
static void sigmund_hop(t_sigmund *x, t_floatarg f)
{
x->x_hop = f;
/* check parameter ranges */
if (x->x_hop != (1 << sigmund_ilog2(x->x_hop)))
post("sigmund~: adjusting analysis size to %d points",
(x->x_hop = (1 << sigmund_ilog2(x->x_hop))));
}
static void sigmund_npeak(t_sigmund *x, t_floatarg f)
{
if (f < 1)
f = 1;
x->x_npeak = f;
}
static void sigmund_maxfreq(t_sigmund *x, t_floatarg f)
{
x->x_maxfreq = f;
}
static void sigmund_vibrato(t_sigmund *x, t_floatarg f)
{
if (f < 0)
f = 0;
x->x_vibrato = f;
}
static void sigmund_stabletime(t_sigmund *x, t_floatarg f)
{
if (f < 0)
f = 0;
x->x_stabletime = f;
}
static void sigmund_growth(t_sigmund *x, t_floatarg f)
{
if (f < 0)
f = 0;
x->x_growth = f;
}
static void sigmund_minpower(t_sigmund *x, t_floatarg f)
{
if (f < 0)
f = 0;
x->x_minpower = f;
}
static void sigmund_print(t_sigmund *x)
{
post("sigmund~ settings:");
post("npts %d", (int)x->x_npts);
post("hop %d", (int)x->x_hop);
post("npeak %d", (int)x->x_npeak);
post("maxfreq %g", x->x_maxfreq);
post("vibrato %g", x->x_vibrato);
post("stabletime %g", x->x_stabletime);
post("growth %g", x->x_growth);
post("minpower %g", x->x_minpower);
}
static void sigmund_printnext(t_sigmund *x, t_float f)
{
x->x_loud = f;
}
static void sigmund_clock(t_sigmund *x)
{
if (x->x_infill == x->x_npts)
{
sigmund_doit(x, x->x_npts, x->x_inbuf, x->x_loud, x->x_sr);
if (x->x_hop >= x->x_npts)
{
x->x_infill = 0;
x->x_countdown = x->x_hop - x->x_npts;
}
else
{
memmove(x->x_inbuf, x->x_inbuf + x->x_hop,
(x->x_infill = x->x_npts - x->x_hop) * sizeof(*x->x_inbuf));
x->x_countdown = 0;
}
x->x_loud = 0;
}
}
static t_int *sigmund_perform(t_int *w)
{
t_sigmund *x = (t_sigmund *)(w[1]);
float *in = (float *)(w[2]);
int n = (int)(w[3]);
if (x->x_countdown > 0)
x->x_countdown -= n;
else if (x->x_infill != x->x_npts)
{
int i, j;
float *fp = x->x_inbuf + x->x_infill;
for (j = 0; j < n; j++)
*fp++ = *in++;
x->x_infill += n;
if (x->x_infill == x->x_npts)
clock_delay(x->x_clock, 0);
}
return (w+4);
}
static void sigmund_dsp(t_sigmund *x, t_signal **sp)
{
if (x->x_mode == MODE_STREAM)
{
if (x->x_hop % sp[0]->s_n)
post("sigmund: adjusting hop size to %d",
(x->x_hop = sp[0]->s_n * (x->x_hop / sp[0]->s_n)));
x->x_sr = sp[0]->s_sr;
dsp_add(sigmund_perform, 3, x, sp[0]->s_vec, sp[0]->s_n);
}
}
static void sigmund_free(t_sigmund *x)
{
if (x->x_inbuf)
freebytes(x->x_inbuf, x->x_npts * sizeof(*x->x_inbuf));
if (x->x_trackv)
freebytes(x->x_trackv, x->x_ntrack * sizeof(*x->x_trackv));
clock_free(x->x_clock);
}
void sigmund_tilde_setup(void)
{
sigmund_class = class_new(gensym("sigmund~"), (t_newmethod)sigmund_new,
(t_method)sigmund_free, sizeof(t_sigmund), 0, A_GIMME, 0);
class_addlist(sigmund_class, sigmund_list);
class_addmethod(sigmund_class, (t_method)sigmund_dsp, gensym("dsp"), 0);
CLASS_MAINSIGNALIN(sigmund_class, t_sigmund, x_f);
class_addmethod(sigmund_class, (t_method)sigmund_param1,
gensym("param1"), A_FLOAT, 0);
class_addmethod(sigmund_class, (t_method)sigmund_param2,
gensym("param2"), A_FLOAT, 0);
class_addmethod(sigmund_class, (t_method)sigmund_param3,
gensym("param3"), A_FLOAT, 0);
class_addmethod(sigmund_class, (t_method)sigmund_npts,
gensym("npts"), A_FLOAT, 0);
class_addmethod(sigmund_class, (t_method)sigmund_hop,
gensym("hop"), A_FLOAT, 0);
class_addmethod(sigmund_class, (t_method)sigmund_maxfreq,
gensym("maxfreq"), A_FLOAT, 0);
class_addmethod(sigmund_class, (t_method)sigmund_npeak,
gensym("npeak"), A_FLOAT, 0);
class_addmethod(sigmund_class, (t_method)sigmund_vibrato,
gensym("vibrato"), A_FLOAT, 0);
class_addmethod(sigmund_class, (t_method)sigmund_stabletime,
gensym("stabletime"), A_FLOAT, 0);
class_addmethod(sigmund_class, (t_method)sigmund_growth,
gensym("growth"), A_FLOAT, 0);
class_addmethod(sigmund_class, (t_method)sigmund_minpower,
gensym("minpower"), A_FLOAT, 0);
class_addmethod(sigmund_class, (t_method)sigmund_print,
gensym("print"), 0);
class_addmethod(sigmund_class, (t_method)sigmund_printnext,
gensym("printnext"), A_FLOAT, 0);
post("sigmund version 0.02");
}