/* here is the warping filter, you give it a value 'd' between -1 and +1
/* and it does a non-linear warp on the formants, down or up.  the formula
/* is in our article in cmj, and I have a couple of routines on the IBM
/* that will translate freq to D and vice versa, I'll get them to you.
/* first the filter	*/

#include "sizes.h"

float cq,outold;

float warppol(sig,past,d,c)
float sig,*past,d,*c;
{
	float temp1,temp2;
	int n;

	temp1 = past[NPOLEM1];
	past[NPOLEM1] = cq * outold - d * past[NPOLEM1];
	for(n=NPOLE-2; n>=0; n--) {
		temp2 = past[n];
		past[n] = d * (past[n+1] - past[n]) + temp1;
		temp1 = temp2;
		}
	for(n=0;n<NPOLE;n++)  sig += c[n] * past[n];
	outold = sig;
	return(sig);
}

warpset(d,c)     
float d,*c;
{
	int m;
	float cl;

	for(m=1; m<NPOLE; m++)   c[m] += d * c[m-1];
	cl = 1./(1.-d * c[NPOLEM1]);
	cq = cl * (1. - d * d);
	return(cl);
}

warpinit()
{
	outold = 0;
}

and now a sample instrument that uses it.  notice that there is an 
extra call at the control rate time

#include "sizes.h"

int anal,pitches;
float sampspf;

play(p,n_args)
float *p;
{
	float c[FRAMSIZE],past[NPOLE*2],val,sig,dur,*cpoint;
	float pasto[NPOLE*2],val2;
	int kcount;
	float allpole(),allpoli();
	int jcount = 0;
	int i,j,framefirst,framelast;
	float d,cl,warpset(),warppol(),getfreq(),peak,temval,dtemp;
	float cps,curtime,cerror,error,thresh;
	int pframe;
	int time;

	anal = open("../lpc/analysis",0);
	pitches = open("pitches",0);
/* p0 is first frame (starting with 0) p1 last frame, p2 is inputskip
   p3 is samples per frame p4 is d p5 is outputskip p6 is thresh to go to d=0*/
	framefirst = p[0];
	framelast  = p[1];
	sampspf = p[3];
	d = p[4];
	thresh = p[6];
	dur = 1./(p[1] - p[0] + 1.);
	cpoint = c+4;
	setnot(p[2],dur,0); /* set input file */
	setnot(p[5],dur,1); /* set output file */
	for(j=0;j<NPOLE*2;j++)  past[j]=pasto[j]=0; 
	jcount = 0;
	kcount = 0;
	warpinit();
	for (i=framefirst; i<=framelast; i++)  {
			getfr(i,c);
			cerror = c[2];
			dtemp = (c[2] > thresh) ? 0. : d;
			cl = warpset(dtemp,c+4);
			curtime = (float)i * sampspf/SR;
			cps = getfreq(curtime,&error,&pframe);
printf("%d %f %d %f %f %f %f %d %f\n",i,peak,time,cl,cps,cerror,error,pframe,curtime);
peak = 0;
cl = 1;
		for(j=0;j<sampspf;j++) {
			getin(&sig,0);
			val = cl * warppol(sig,pasto,dtemp,c+4);
temval = (val > 0.) ? val : - val;
if(temval > peak) { time = j; peak = temval; }
			putout(&val,1);        
			}
	}
	endnot(1);
}