""

svn path=/trunk/; revision=325

14 files changed, 1356 insertions, 0 deletions

diff --git a/externals/grill/fftease/src/burrow~.cpp b/externals/grill/fftease/src/burrow~.cpp
new file mode 100644
index 00000000..7fbc45b8
--- /dev/null
+++ b/externals/grill/fftease/src/burrow~.cpp
@@ -0,0 +1,255 @@
+/*
+
+FFTease - A set of Live Spectral Processors
+Originally written by Eric Lyon and Christopher Penrose for the Max/MSP platform
+
+This flext port is based on the jMax port of Christian Klippel
+
+Copyright (c)Thomas Grill (xovo@gmx.net)
+For information on usage and redistribution, and for a DISCLAIMER OF ALL
+WARRANTIES, see the file, "license.txt," in this distribution.  
+
+*/
+
+#include "main.h"
+
+
+class burrow:
+	public flext_dsp
+{
+	FLEXT_HEADER_S(burrow,flext_dsp,setup)
+	
+public:
+	burrow(I argc,const t_atom *argv);
+	~burrow();
+
+protected:
+
+	virtual V m_dsp(I n,S *const *in,S *const *out);
+	virtual V m_signal(I n,S *const *in,S *const *out);
+
+    BL _invert;
+    I _inCount;
+    I *_bitshuffle;
+
+    F _threshold,_multiplier;
+    F _thresh_dB,_mult_dB;
+
+    F *_Wanal;
+    F *_Wsyn;
+    F *_inputOne,*_inputTwo;
+    F *_Hwin;
+    F *_bufferOne,*_bufferTwo;
+    F *_channelOne,*_channelTwo;
+    F *_output;
+    F *_trigland;
+
+private:
+	V Reset();
+	V Clear();
+	
+	V ms_thresh(F v) { _threshold = (float) (pow( 10., ((_thresh_dB = v) * .05))); }
+	V ms_mult(F v) { _multiplier = (float) (pow( 10., ((_mult_dB = v) * .05))); }
+
+
+	static V setup(t_classid c);
+
+	FLEXT_ATTRGET_F(_thresh_dB)
+	FLEXT_CALLSET_F(ms_thresh)
+	FLEXT_ATTRGET_F(_mult_dB)
+	FLEXT_CALLSET_F(ms_mult)
+	FLEXT_ATTRVAR_B(_invert)
+};
+
+FLEXT_LIB_DSP_V("fftease, burrow~",burrow)
+
+
+V burrow::setup(t_classid c)
+{
+	FLEXT_CADDATTR_VAR(c,"thresh",_thresh_dB,ms_thresh);
+	FLEXT_CADDATTR_VAR(c,"mult",_mult_dB,ms_mult);
+	FLEXT_CADDATTR_VAR1(c,"invert",_invert);
+}
+
+
+burrow::burrow(I argc,const t_atom *argv):
+	_thresh_dB(-30),_mult_dB(-18),
+	_invert(false)
+{
+	/* parse and set object's options given */
+	if(argc >= 1) {
+		if(CanbeFloat(argv[0]))
+			_thresh_dB = GetAFloat(argv[0]);
+		else
+			post("%s - Threshold must be a float value - set to %0f",thisName(),_thresh_dB);
+	}
+	if(argc >= 2) {
+		if(CanbeFloat(argv[1]))
+			_mult_dB = GetAFloat(argv[1]);
+		else
+			post("%s - Multiplier must be a float value - set to %0f",thisName(),_mult_dB);
+	}
+	if(argc >= 3) {
+		if(CanbeBool(argv[2]))
+			_invert = GetABool(argv[2]);
+		else
+			post("%s - Invert flags must be a boolean value - set to %i",thisName(),_invert?1:0);
+	}
+
+	ms_thresh(_thresh_dB);
+	ms_mult(_mult_dB);
+
+	Reset();
+
+	AddInSignal(2);
+	AddOutSignal();
+}
+
+burrow::~burrow()
+{
+	Clear();
+}
+
+V burrow::Reset()
+{
+	_bitshuffle = NULL;
+	_trigland = NULL;
+	_inputOne = _inputTwo = NULL;
+	_Hwin = NULL;
+	_Wanal = _Wsyn = NULL;
+	_bufferOne = _bufferTwo = NULL;
+	_channelOne = _channelTwo = NULL;
+	_output = NULL;
+}
+
+V burrow::Clear() 
+{
+	if(_bitshuffle) delete[] _bitshuffle;
+	if(_trigland) delete[] _trigland;
+	if(_inputOne) delete[] _inputOne;
+	if(_inputTwo) delete[] _inputTwo;
+	if(_Hwin) delete[] _Hwin;
+	if(_Wanal) delete[] _Wanal;
+	if(_Wsyn) delete[] _Wsyn;
+	if(_bufferOne) delete[] _bufferOne;
+	if(_bufferTwo) delete[] _bufferTwo;
+	if(_channelOne) delete[] _channelOne;
+	if(_channelTwo) delete[] _channelTwo;
+	if(_output) delete[] _output;
+}
+
+
+
+V burrow::m_dsp(I n,S *const *in,S *const *out)
+{
+	Clear();
+
+	/* preset the objects data */
+	const I _D = Blocksize();
+	const I _N = _D* 4,_Nw = _N,_N2 = _N / 2,_Nw2 = _Nw / 2; 
+
+	_inCount = -_Nw;
+
+	/* assign memory to the buffers */
+	_bitshuffle = new I[_N*2];
+	_trigland = new F[_N*2];
+	_inputOne = new F[_Nw];
+	_inputTwo = new F[_Nw];
+	_Hwin = new F[_Nw];
+	_Wanal = new F[_Nw];
+	_Wsyn = new F[_Nw];
+	_bufferOne = new F[_N];
+	_bufferTwo = new F[_N];
+	_channelOne = new F[_N+2];
+	_channelTwo = new F[_N+2];
+	_output = new F[_Nw];
+
+	/* initialize pv-lib functions */
+	init_rdft( _N, _bitshuffle, _trigland);
+	makewindows( _Hwin, _Wanal, _Wsyn, _Nw, _N, _D, 0);
+}
+
+V burrow::m_signal(I n,S *const *in,S *const *out)
+{
+	const S *inOne = in[0],*inTwo = in[1];
+	S *outOne = out[0];
+
+	/* declare working variables */
+	I i, j, even, odd; 
+	const I _D = Blocksize();
+	const I _N = _D* 4,_Nw = _N,_N2 = _N / 2,_Nw2 = _Nw / 2; 
+
+	/* fill our retaining buffers */
+	_inCount += _D;
+
+	for(i = 0; i < _N-_D ; i++ ) {
+		_inputOne[i] = _inputOne[i+_D];
+		_inputTwo[i] = _inputTwo[i+_D];
+	}
+	for(j = 0; i < _N; i++,j++) {
+		_inputOne[i] = inOne[j];
+		_inputTwo[i] = inTwo[j];
+	}
+
+	/* apply hamming window and fold our window buffer into the fft buffer */
+	fold( _inputOne, _Wanal, _Nw, _bufferOne, _N, _inCount );
+	fold( _inputTwo, _Wanal, _Nw, _bufferTwo, _N, _inCount );
+
+	/* do an fft */
+	rdft( _N, 1, _bufferOne, _bitshuffle, _trigland );
+	rdft( _N, 1, _bufferTwo, _bitshuffle, _trigland );
+
+
+	/* convert to polar coordinates from complex values */
+	for ( i = 0; i <= _N2; i++ ) {
+		register F a,b;
+
+		odd = ( even = i<<1 ) + 1;
+
+		a = ( i == _N2 ? _bufferOne[1] : _bufferOne[even] );
+		b = ( i == 0 || i == _N2 ? 0. : _bufferOne[odd] );
+
+		_channelOne[even] = hypot( a, b );
+		_channelOne[odd] = -atan2( b, a );
+
+		a = ( i == _N2 ? _bufferTwo[1] : _bufferTwo[even] );
+		b = ( i == 0 || i == _N2 ? 0. : _bufferTwo[odd] );
+
+		_channelTwo[even] = hypot( a, b );
+
+		/* use simple threshold from second signal to trigger filtering */
+		if (_invert?(_channelTwo[even] < _threshold):(_channelTwo[even] > _threshold) )
+			_channelOne[even] *= _multiplier;
+	}
+
+	/* convert back to complex form, read for the inverse fft */
+	for ( i = 0; i <= _N2; i++ ) {
+		odd = ( even = i<<1 ) + 1;
+
+		*(_bufferOne+even) = *(_channelOne+even) * cos( *(_channelOne+odd) );
+
+		if ( i != _N2 )
+			*(_bufferOne+odd) = -(*(_channelOne+even)) * sin( *(_channelOne+odd) );
+	}
+
+	/* do an inverse fft */
+	rdft( _N, -1, _bufferOne, _bitshuffle, _trigland );
+
+	/* dewindow our result */
+	overlapadd( _bufferOne, _N, _Wsyn, _output, _Nw, _inCount);
+
+	/* set our output and adjust our retaining output buffer */
+	F mult = 1./_N;
+	for ( j = 0; j < _D; j++ )
+		outOne[j] = _output[j] * mult;
+
+	for ( j = 0; j < _N-_D; j++ )
+		_output[j] = _output[j+_D];
+	for (; j < _N; j++ )
+		_output[j] = 0.;
+
+	/* restore state variables */
+}
+
+
+
diff --git a/externals/grill/fftease/src/convert.c b/externals/grill/fftease/src/convert.c
new file mode 100644
index 00000000..8eb238f0
--- /dev/null
+++ b/externals/grill/fftease/src/convert.c
@@ -0,0 +1,55 @@
+#include "pv.h"
+
+/* S is a spectrum in rfft format, i.e., it contains N real values
+   arranged as real followed by imaginary values, except for first
+   two values, which are real parts of 0 and Nyquist frequencies;
+   convert first changes these into N/2+1 PAIRS of magnitude and
+   phase values to be stored in output array C; the phases are then
+   unwrapped and successive phase differences are used to compute
+   estimates of the instantaneous frequencies for each phase vocoder
+   analysis channel; decimation rate D and sampling rate R are used
+   to render these frequency values directly in Hz. */
+
+void convert(float *S, float *C, int N2, float *lastphase,  float fundamental, float factor )
+{
+  float 	phase,
+		phasediff;
+  int 		real,
+		imag,
+		amp,
+		freq;
+  float 	a,
+		b;
+  int 		i;
+
+  float myTWOPI, myPI;
+
+  myTWOPI = 8.*atan(1.);
+  myPI = 4.*atan(1.);
+
+
+    for ( i = 0; i <= N2; i++ ) {
+      imag = freq = ( real = amp = i<<1 ) + 1;
+      a = ( i == N2 ? S[1] : S[real] );
+      b = ( i == 0 || i == N2 ? 0. : S[imag] );
+
+      C[amp] = hypot( a, b );
+      if ( C[amp] == 0. )
+	phasediff = 0.;
+      else {
+	phasediff = ( phase = -atan2( b, a ) ) - lastphase[i];
+	lastphase[i] = phase;
+	
+	while ( phasediff > myPI )
+	  phasediff -= myTWOPI;
+	while ( phasediff < -myPI )
+	  phasediff += myTWOPI;
+      }
+      C[freq] = phasediff*factor + i*fundamental;
+      /*
+      if( i > 8 && i < 12 ) {
+	fprintf(stderr,"convert freq %d: %f\n",i, C[freq]);
+      }
+      */
+    }
+}
diff --git a/externals/grill/fftease/src/convert_new.c b/externals/grill/fftease/src/convert_new.c
new file mode 100644
index 00000000..5df717f8
--- /dev/null
+++ b/externals/grill/fftease/src/convert_new.c
@@ -0,0 +1,64 @@
+#include "pv.h"  // T.Grill
+
+// #include <stdio.h>
+#include <math.h>
+
+
+
+/* S is a spectrum in rfft format, i.e., it contains N real values
+   arranged as real followed by imaginary values, except for first
+   two values, which are real parts of 0 and Nyquist frequencies;
+   convert first changes these into N/2+1 PAIRS of magnitude and
+   phase values to be stored in output array C; the phases are then
+   unwrapped and successive phase differences are used to compute
+   estimates of the instantaneous frequencies for each phase vocoder
+   analysis channel; decimation rate D and sampling rate R are used
+   to render these frequency values directly in Hz. */
+
+
+
+
+
+void convert_new(float *S, float *C, int N2, float *lastphase,  float fundamental, float factor )
+{
+  float 	phase,
+		phasediff;
+  int 		real,
+		imag,
+		amp,
+		freq;
+  float 	a,
+		b;
+  int 		i;
+
+  float myTWOPI, myPI;
+  
+  myTWOPI = 8.*atan(1.);
+  myPI = 4.*atan(1.);
+
+
+    for ( i = 0; i <= N2; i++ ) {
+      imag = freq = ( real = amp = i<<1 ) + 1;
+      a = ( i == N2 ? S[1] : S[real] );
+      b = ( i == 0 || i == N2 ? 0. : S[imag] );
+
+      C[amp] = hypot( a, b );
+      if ( C[amp] == 0. )
+	phasediff = 0.;
+      else {
+	phasediff = ( phase = -atan2( b, a ) ) - lastphase[i];
+	lastphase[i] = phase;
+	
+	while ( phasediff > myPI )
+	  phasediff -= myTWOPI;
+	while ( phasediff < -myPI )
+	  phasediff += myTWOPI;
+      }
+      C[freq] = phasediff*factor + i*fundamental;
+      /*
+      if( i > 8 && i < 12 ) {
+	fprintf(stderr,"convert freq %d: %f\n",i, C[freq]);
+      }
+      */
+    }
+}
diff --git a/externals/grill/fftease/src/fft.c b/externals/grill/fftease/src/fft.c
new file mode 100644
index 00000000..f920e231
--- /dev/null
+++ b/externals/grill/fftease/src/fft.c
@@ -0,0 +1,161 @@
+#include "pv.h"
+
+
+//float TWOPI;
+
+/* If forward is true, rfft replaces 2*N real data points in x with
+   N complex values representing the positive frequency half of their
+   Fourier spectrum, with x[1] replaced with the real part of the Nyquist
+   frequency value.  If forward is false, rfft expects x to contain a
+   positive frequency spectrum arranged as before, and replaces it with
+   2*N real values.  N MUST be a power of 2. */
+static void 
+bitreverse( float *x, int N );
+
+void rfft( float *x, int N, int forward )
+{
+  float 	c1,c2,
+  		h1r,h1i,
+		h2r,h2i,
+		wr,wi,
+		wpr,wpi,
+  		temp,
+		theta;
+  float 	xr,xi;
+  int 		i,
+		i1,i2,i3,i4,
+		N2p1;
+  static int 	first = 1;
+float PI, TWOPI;
+
+PI = 3.141592653589793115997963468544185161590576171875;
+TWOPI = 6.28318530717958623199592693708837032318115234375;
+    if ( first ) {
+
+	first = 0;
+    }
+    theta = PI/N;
+    wr = 1.;
+    wi = 0.;
+    c1 = 0.5;
+    if ( forward ) {
+	c2 = -0.5;
+	cfft( x, N, forward );
+	xr = x[0];
+	xi = x[1];
+    } else {
+	c2 = 0.5;
+	theta = -theta;
+	xr = x[1];
+	xi = 0.;
+	x[1] = 0.;
+    }
+    wpr = -2.*pow( sin( 0.5*theta ), 2. );
+    wpi = sin( theta );
+    N2p1 = (N<<1) + 1;
+    for ( i = 0; i <= N>>1; i++ ) {
+	i1 = i<<1;
+	i2 = i1 + 1;
+	i3 = N2p1 - i2;
+	i4 = i3 + 1;
+	if ( i == 0 ) {
+	    h1r =  c1*(x[i1] + xr );
+	    h1i =  c1*(x[i2] - xi );
+	    h2r = -c2*(x[i2] + xi );
+	    h2i =  c2*(x[i1] - xr );
+	    x[i1] =  h1r + wr*h2r - wi*h2i;
+	    x[i2] =  h1i + wr*h2i + wi*h2r;
+	    xr =  h1r - wr*h2r + wi*h2i;
+	    xi = -h1i + wr*h2i + wi*h2r;
+	} else {
+	    h1r =  c1*(x[i1] + x[i3] );
+	    h1i =  c1*(x[i2] - x[i4] );
+	    h2r = -c2*(x[i2] + x[i4] );
+	    h2i =  c2*(x[i1] - x[i3] );
+	    x[i1] =  h1r + wr*h2r - wi*h2i;
+	    x[i2] =  h1i + wr*h2i + wi*h2r;
+	    x[i3] =  h1r - wr*h2r + wi*h2i;
+	    x[i4] = -h1i + wr*h2i + wi*h2r;
+	}
+	wr = (temp = wr)*wpr - wi*wpi + wr;
+	wi = wi*wpr + temp*wpi + wi;
+    }
+    if ( forward )
+	x[1] = xr;
+    else
+	cfft( x, N, forward );
+}
+
+/* cfft replaces float array x containing NC complex values
+   (2*NC float values alternating real, imagininary, etc.)
+   by its Fourier transform if forward is true, or by its
+   inverse Fourier transform if forward is false, using a
+   recursive Fast Fourier transform method due to Danielson
+   and Lanczos.  NC MUST be a power of 2. */
+
+void cfft( float *x, int NC, int forward )
+{
+  float 	wr,wi,
+		wpr,wpi,
+		theta,
+		scale;
+  int 		mmax,
+		ND,
+		m,
+		i,j,
+		delta;
+float TWOPI;
+TWOPI = 6.28318530717958623199592693708837032318115234375;
+    ND = NC<<1;
+    bitreverse( x, ND );
+    for ( mmax = 2; mmax < ND; mmax = delta ) {
+	delta = mmax<<1;
+	theta = TWOPI/( forward? mmax : -mmax );
+	wpr = -2.*pow( sin( 0.5*theta ), 2. );
+	wpi = sin( theta );
+	wr = 1.;
+	wi = 0.;
+	for ( m = 0; m < mmax; m += 2 ) {
+	 register float rtemp, itemp;
+	    for ( i = m; i < ND; i += delta ) {
+		j = i + mmax;
+		rtemp = wr*x[j] - wi*x[j+1];
+		itemp = wr*x[j+1] + wi*x[j];
+		x[j] = x[i] - rtemp;
+		x[j+1] = x[i+1] - itemp;
+		x[i] += rtemp;
+		x[i+1] += itemp;
+	    }
+	    wr = (rtemp = wr)*wpr - wi*wpi + wr;
+	    wi = wi*wpr + rtemp*wpi + wi;
+	}
+    }
+
+/* scale output */
+
+    scale = forward ? 1./ND : 2.;
+    { register float *xi=x, *xe=x+ND;
+	while ( xi < xe )
+	    *xi++ *= scale;
+    }
+}
+
+/* bitreverse places float array x containing N/2 complex values
+   into bit-reversed order */
+
+void bitreverse( float *x, int N )
+{
+  float 	rtemp,itemp;
+  int 		i,j,
+		m;
+
+    for ( i = j = 0; i < N; i += 2, j += m ) {
+	if ( j > i ) {
+	    rtemp = x[j]; itemp = x[j+1]; /* complex exchange */
+	    x[j] = x[i]; x[j+1] = x[i+1];
+	    x[i] = rtemp; x[i+1] = itemp;
+	}
+	for ( m = N>>1; m >= 2 && j >= m; m >>= 1 )
+	    j -= m;
+    }
+}
diff --git a/externals/grill/fftease/src/fft4.c b/externals/grill/fftease/src/fft4.c
new file mode 100644
index 00000000..ffb3574d
--- /dev/null
+++ b/externals/grill/fftease/src/fft4.c
@@ -0,0 +1,340 @@
+/*
+ *  This file is part of the pv-lib
+ *
+ *  The pv-lib can be used by everyone as desired
+ *
+ *  (c) Eric Lyon and Christopher Penrose
+ *
+ */
+  
+#include <math.h>
+#include "pv.h"
+
+static void bitrv2(int n, int *ip, float *a);
+static void cftsub(int n, float *a, float *w);
+static void rftsub(int n, float *a, int nc, float *c);
+static void makewt(int nw, int *ip, float *w);
+static void makect(int nc, int *ip, float *c);
+
+void init_rdft(int n, int *ip, float *w)
+{
+
+  int	nw,
+	nc;
+
+  nw = n >> 2;
+  makewt(nw, ip, w);
+
+  nc = n >> 2;
+  makect(nc, ip, w + nw);
+
+  return;
+}
+
+
+void rdft(int n, int isgn, float *a, int *ip, float *w)
+{
+
+  int		j,
+		nw,
+		nc;
+
+  float		xi;
+
+  nw = ip[0];
+  nc = ip[1];
+
+  if (isgn < 0) {
+    a[1] = 0.5 * (a[1] - a[0]);
+    a[0] += a[1];
+
+    for (j = 3; j <= n - 1; j += 2) {
+      a[j] = -a[j];
+    }
+
+    if (n > 4) {
+      rftsub(n, a, nc, w + nw);
+      bitrv2(n, ip + 2, a);
+    }
+
+    cftsub(n, a, w);
+
+    for (j = 1; j < n; j += 2) {
+      a[j] = -a[j];
+    }
+  }
+
+  else {
+
+    if (n > 4) {
+      bitrv2(n, ip + 2, a);
+    }
+
+    cftsub(n, a, w);
+
+    if (n > 4) {
+      rftsub(n, a, nc, w + nw);
+    }
+
+    xi = a[0] - a[1];
+    a[0] += a[1];
+    a[1] = xi;
+  }
+}
+
+
+static void bitrv2(int n, int *ip, float *a)
+{
+  int j, j1, k, k1, l, m, m2;
+  float xr, xi;
+
+  ip[0] = 0;
+  l = n;
+  m = 1;
+
+  while ((m << 2) < l) {
+    l >>= 1;
+    for (j = 0; j < m; j++) {
+      ip[m + j] = ip[j] + l;
+    }
+    m <<= 1;
+  }
+
+  if ((m << 2) > l) {
+
+    for (k = 1; k < m; k++) {
+
+      for (j = 0; j <= k - 1; j++) {
+	j1 = (j << 1) + ip[k];
+	k1 = (k << 1) + ip[j];
+	xr = a[j1];
+	xi = a[j1 + 1];
+	a[j1] = a[k1];
+	a[j1 + 1] = a[k1 + 1];
+	a[k1] = xr;
+	a[k1 + 1] = xi;
+      }
+    }
+  }
+
+  else {
+    m2 = m << 1;
+
+    for (k = 1; k < m; k++) {
+
+      for (j = 0; j < k; j++) {
+	j1 = (j << 1) + ip[k];
+	k1 = (k << 1) + ip[j];
+	xr = a[j1];
+	xi = a[j1 + 1];
+	a[j1] = a[k1];
+	a[j1 + 1] = a[k1 + 1];
+	a[k1] = xr;
+	a[k1 + 1] = xi;
+	j1 += m2;
+	k1 += m2;
+	xr = a[j1];
+	xi = a[j1 + 1];
+	a[j1] = a[k1];
+	a[j1 + 1] = a[k1 + 1];
+	a[k1] = xr;
+	a[k1 + 1] = xi;
+      }
+    }
+  }
+}
+
+
+static void cftsub(int n, float *a, float *w)
+{
+  int j, j1, j2, j3, k, k1, ks, l, m;
+  float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
+  float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
+    
+  l = 2;
+
+  while ((l << 1) < n) {
+    m = l << 2;
+
+    for (j = 0; j <= l - 2; j += 2) {
+      j1 = j + l;
+      j2 = j1 + l;
+      j3 = j2 + l;
+      x0r = a[j] + a[j1];
+      x0i = a[j + 1] + a[j1 + 1];
+      x1r = a[j] - a[j1];
+      x1i = a[j + 1] - a[j1 + 1];
+      x2r = a[j2] + a[j3];
+      x2i = a[j2 + 1] + a[j3 + 1];
+      x3r = a[j2] - a[j3];
+      x3i = a[j2 + 1] - a[j3 + 1];
+      a[j] = x0r + x2r;
+      a[j + 1] = x0i + x2i;
+      a[j2] = x0r - x2r;
+      a[j2 + 1] = x0i - x2i;
+      a[j1] = x1r - x3i;
+      a[j1 + 1] = x1i + x3r;
+      a[j3] = x1r + x3i;
+      a[j3 + 1] = x1i - x3r;
+    }
+
+    if (m < n) {
+      wk1r = w[2];
+
+      for (j = m; j <= l + m - 2; j += 2) {
+	j1 = j + l;
+	j2 = j1 + l;
+	j3 = j2 + l;
+	x0r = a[j] + a[j1];
+	x0i = a[j + 1] + a[j1 + 1];
+	x1r = a[j] - a[j1];
+	x1i = a[j + 1] - a[j1 + 1];
+	x2r = a[j2] + a[j3];
+	x2i = a[j2 + 1] + a[j3 + 1];
+	x3r = a[j2] - a[j3];
+	x3i = a[j2 + 1] - a[j3 + 1];
+	a[j] = x0r + x2r;
+	a[j + 1] = x0i + x2i;
+	a[j2] = x2i - x0i;
+	a[j2 + 1] = x0r - x2r;
+	x0r = x1r - x3i;
+	x0i = x1i + x3r;
+	a[j1] = wk1r * (x0r - x0i);
+	a[j1 + 1] = wk1r * (x0r + x0i);
+	x0r = x3i + x1r;
+	x0i = x3r - x1i;
+	a[j3] = wk1r * (x0i - x0r);
+	a[j3 + 1] = wk1r * (x0i + x0r);
+      }
+
+      k1 = 1;
+      ks = -1;
+
+      for (k = (m << 1); k <= n - m; k += m) {
+	k1++;
+	ks = -ks;
+	wk1r = w[k1 << 1];
+	wk1i = w[(k1 << 1) + 1];
+	wk2r = ks * w[k1];
+	wk2i = w[k1 + ks];
+	wk3r = wk1r - 2 * wk2i * wk1i;
+	wk3i = 2 * wk2i * wk1r - wk1i;
+
+	for (j = k; j <= l + k - 2; j += 2) {
+	  j1 = j + l;
+	  j2 = j1 + l;
+	  j3 = j2 + l;
+	  x0r = a[j] + a[j1];
+	  x0i = a[j + 1] + a[j1 + 1];
+	  x1r = a[j] - a[j1];
+	  x1i = a[j + 1] - a[j1 + 1];
+	  x2r = a[j2] + a[j3];
+	  x2i = a[j2 + 1] + a[j3 + 1];
+	  x3r = a[j2] - a[j3];
+	  x3i = a[j2 + 1] - a[j3 + 1];
+	  a[j] = x0r + x2r;
+	  a[j + 1] = x0i + x2i;
+	  x0r -= x2r;
+	  x0i -= x2i;
+	  a[j2] = wk2r * x0r - wk2i * x0i;
+	  a[j2 + 1] = wk2r * x0i + wk2i * x0r;
+	  x0r = x1r - x3i;
+	  x0i = x1i + x3r;
+	  a[j1] = wk1r * x0r - wk1i * x0i;
+	  a[j1 + 1] = wk1r * x0i + wk1i * x0r;
+	  x0r = x1r + x3i;
+	  x0i = x1i - x3r;
+	  a[j3] = wk3r * x0r - wk3i * x0i;
+	  a[j3 + 1] = wk3r * x0i + wk3i * x0r;
+	}
+      }
+    }
+
+    l = m;
+  }
+
+  if (l < n) {
+
+    for (j = 0; j <= l - 2; j += 2) {
+      j1 = j + l;
+      x0r = a[j] - a[j1];
+      x0i = a[j + 1] - a[j1 + 1];
+      a[j] += a[j1];
+      a[j + 1] += a[j1 + 1];
+      a[j1] = x0r;
+      a[j1 + 1] = x0i;
+    }
+  }
+}
+
+
+static void rftsub(int n, float *a, int nc, float *c)
+{
+  int j, k, kk, ks;
+  float wkr, wki, xr, xi, yr, yi;
+
+  ks = (nc << 2) / n;
+  kk = 0;
+
+  for (k = (n >> 1) - 2; k >= 2; k -= 2) {
+    j = n - k;
+    kk += ks;
+    wkr = 0.5 - c[kk];
+    wki = c[nc - kk];
+    xr = a[k] - a[j];
+    xi = a[k + 1] + a[j + 1];
+    yr = wkr * xr - wki * xi;
+    yi = wkr * xi + wki * xr;
+    a[k] -= yr;
+    a[k + 1] -= yi;
+    a[j] += yr;
+    a[j + 1] -= yi;
+  }
+}
+
+
+static void makewt(int nw, int *ip, float *w)
+{
+    int nwh, j;
+    float delta, x, y;
+
+    ip[0] = nw;
+    ip[1] = 1;
+    if (nw > 2) {
+        nwh = nw >> 1;
+        delta = atan(1.0) / nwh;
+        w[0] = 1;
+        w[1] = 0;
+        w[nwh] = cos(delta * nwh);
+        w[nwh + 1] = w[nwh];
+        for (j = 2; j <= nwh - 2; j += 2) {
+            x = cos(delta * j);
+            y = sin(delta * j);
+            w[j] = x;
+            w[j + 1] = y;
+            w[nw - j] = y;
+            w[nw - j + 1] = x;
+        }
+        bitrv2(nw, ip + 2, w);
+    }
+}
+
+
+static void makect(int nc, int *ip, float *c)
+{
+    int nch, j;
+    float delta;
+
+    ip[1] = nc;
+    if (nc > 1) {
+        nch = nc >> 1;
+        delta = atan(1.0) / nch;
+        c[0] = 0.5;
+        c[nch] = 0.5 * cos(delta * nch);
+        for (j = 1; j <= nch - 1; j++) {
+            c[j] = 0.5 * cos(delta * j);
+            c[nc - j] = 0.5 * sin(delta * j);
+        }
+    }
+}
+
diff --git a/externals/grill/fftease/src/fold.c b/externals/grill/fftease/src/fold.c
new file mode 100644
index 00000000..0ecee5d6
--- /dev/null
+++ b/externals/grill/fftease/src/fold.c
@@ -0,0 +1,24 @@
+/*
+ * multiply current input I by window W (both of length Nw);
+ * using modulus arithmetic, fold and rotate windowed input
+ * into output array O of (FFT) length N according to current
+ * input time n
+ */
+void fold( float *I, float  *W, int Nw, float *O, int N, int n )
+{
+ 
+    int i;
+
+    for ( i = 0; i < N; i++ )
+	O[i] = 0.;
+
+    while ( n < 0 )
+	n += N;
+    n %= N;
+    for ( i = 0; i < Nw; i++ ) {
+	O[n] += I[i]*W[i];
+	if ( ++n == N )
+	    n = 0;
+    }
+}
+
diff --git a/externals/grill/fftease/src/leanconvert.c b/externals/grill/fftease/src/leanconvert.c
new file mode 100644
index 00000000..36b1452a
--- /dev/null
+++ b/externals/grill/fftease/src/leanconvert.c
@@ -0,0 +1,19 @@
+#include "pv.h"
+
+void leanconvert( float *S, float *C, int N2 )
+{
+
+ int		real, imag,
+		amp, phase;
+ float		a, b;
+  int		i;
+	
+ for ( i = 0; i <= N2; i++ ) {
+   imag = phase = ( real = amp = i<<1 ) + 1;
+   a = ( i == N2 ? S[1] : S[real] );
+   b = ( i == 0 || i == N2 ? 0. : S[imag] );
+   C[amp] = hypot( a, b );
+   C[phase] = -atan2( b, a );
+ }
+}
+
diff --git a/externals/grill/fftease/src/leanunconvert.c b/externals/grill/fftease/src/leanunconvert.c
new file mode 100644
index 00000000..f8575168
--- /dev/null
+++ b/externals/grill/fftease/src/leanunconvert.c
@@ -0,0 +1,23 @@
+#include "pv.h"
+
+/* unconvert essentially undoes what convert does, i.e., it
+  turns N2+1 PAIRS of amplitude and frequency values in
+  C into N2 PAIR of complex spectrum data (in rfft format)
+  in output array S; sampling rate R and interpolation factor
+  I are used to recompute phase values from frequencies */
+
+void leanunconvert( float *C, float *S, int N2 )
+{
+  int		real, imag,
+		amp, phase;
+  float		a, b;
+  register int		i;
+  
+  for ( i = 0; i <= N2; i++ ) {
+    imag = phase = ( real = amp = i<<1 ) + 1;
+    S[real] = *(C+amp) * cos( *(C+phase) );
+    if ( i != N2 )
+      S[imag] = -*(C+amp) * sin( *(C+phase) );
+  }
+}
+
diff --git a/externals/grill/fftease/src/main.cpp b/externals/grill/fftease/src/main.cpp
new file mode 100644
index 00000000..1d627cf2
--- /dev/null
+++ b/externals/grill/fftease/src/main.cpp
@@ -0,0 +1,60 @@
+/*
+
+FFTease - A set of Live Spectral Processors
+Originally written by Eric Lyon and Christopher Penrose for the Max/MSP platform
+
+This flext port is based on the jMax port of Christian Klippel
+
+Copyright (c)Thomas Grill (xovo@gmx.net)
+For information on usage and redistribution, and for a DISCLAIMER OF ALL
+WARRANTIES, see the file, "license.txt," in this distribution.  
+
+*/
+
+#include "main.h"
+
+
+// Initialization function for xsample library
+V lib_setup()
+{
+	post("");
+	post("-------------------------------------------------------");
+	post("FFTease - A set of Live Spectral Processors");
+	post("Originally written by Eric Lyon and Christopher Penrose");
+	post("for the MAX/MSP platform.");
+	post("");
+	post("flext port done by Thomas Grill, (C)2003");
+	post("-------------------------------------------------------");
+	post("");
+
+	// call the objects' setup routines
+	FLEXT_DSP_SETUP(burrow);
+
+/*
+	FLEXT_DSP_SETUP(cross);
+	FLEXT_DSP_SETUP(dentist);
+	FLEXT_DSP_SETUP(disarray);
+	FLEXT_DSP_SETUP(drown);
+	FLEXT_DSP_SETUP(ether);
+
+	FLEXT_DSP_SETUP(morphine);
+	FLEXT_DSP_SETUP(pvcompand);
+	FLEXT_DSP_SETUP(pvoc);
+	FLEXT_DSP_SETUP(scrape);
+	FLEXT_DSP_SETUP(shapee);
+
+	FLEXT_DSP_SETUP(swinger);
+	FLEXT_DSP_SETUP(taint);
+	FLEXT_DSP_SETUP(thresher);
+	FLEXT_DSP_SETUP(vacancy);
+	FLEXT_DSP_SETUP(xsyn);
+*/
+
+#if FLEXT_SYS == FLEXT_SYS_MAX
+//	finder_addclass((C *)"FFTease",(C *)"xxxxx");
+#endif
+}
+
+// setup the library
+FLEXT_LIB_SETUP(fftease,lib_setup)
+
diff --git a/externals/grill/fftease/src/main.h b/externals/grill/fftease/src/main.h
new file mode 100644
index 00000000..c0fb56f8
--- /dev/null
+++ b/externals/grill/fftease/src/main.h
@@ -0,0 +1,40 @@
+/*
+
+FFTease - A set of Live Spectral Processors
+Originally written by Eric Lyon and Christopher Penrose for the Max/MSP platform
+
+This flext port is based on the jMax port of Christian Klippel
+
+Copyright (c)Thomas Grill (xovo@gmx.net)
+For information on usage and redistribution, and for a DISCLAIMER OF ALL
+WARRANTIES, see the file, "license.txt," in this distribution.  
+
+*/
+
+#ifndef __FFTEASE_H
+#define __FFTEASE_H
+
+#define FFTEASE_VERSION "0.0.1"
+
+#define FLEXT_ATTRIBUTES 1 
+
+#include <flext.h>
+
+#if !defined(FLEXT_VERSION) || (FLEXT_VERSION < 401)
+#error You need at least flext version 0.4.1
+#endif
+
+#include "pv.h"
+
+// lazy me
+#define F float
+#define D double
+#define I int
+#define L long
+#define C char
+#define V void
+#define BL bool
+#define S t_sample
+
+
+#endif
diff --git a/externals/grill/fftease/src/makewindows.c b/externals/grill/fftease/src/makewindows.c
new file mode 100644
index 00000000..1287390e
--- /dev/null
+++ b/externals/grill/fftease/src/makewindows.c
@@ -0,0 +1,202 @@
+#include "pv.h"
+#include "math.h"
+
+/*
+ * make balanced pair of analysis (A) and synthesis (S) windows;
+ * window lengths are Nw, FFT length is N, synthesis interpolation
+ * factor is I, and osc is true (1) if oscillator bank resynthesis 
+ * is specified
+ */
+void makewindows( float *H, float *A, float *S, int Nw, int N, int I, int osc )
+{
+ int i ;
+ float sum ;
+float PI, TWOPI;
+
+PI = 3.141592653589793115997963468544185161590576171875;
+TWOPI = 6.28318530717958623199592693708837032318115234375;
+
+/*
+ * basic Hamming windows
+ */
+    for ( i = 0 ; i < Nw ; i++ )
+	H[i] = A[i] = S[i] = 0.54 - 0.46*cos( TWOPI*i/(Nw - 1) ) ;
+/*
+ * when Nw > N, also apply interpolating (sinc) windows to
+ * ensure that window are 0 at increments of N (the FFT length)
+ * away from the center of the analysis window and of I away
+ * from the center of the synthesis window
+ */
+    if ( Nw > N ) {
+     float x ;
+
+/*
+ * take care to create symmetrical windows
+ */
+	x = -(Nw - 1)/2. ;
+	for ( i = 0 ; i < Nw ; i++, x += 1. )
+	    if ( x != 0. ) {
+		A[i] *= N*sin( PI*x/N )/(PI*x) ;
+		if ( I )
+		    S[i] *= I*sin( PI*x/I )/(PI*x) ;
+	    }
+    }
+/*
+ * normalize windows for unity gain across unmodified
+ * analysis-synthesis procedure
+ */
+    for ( sum = i = 0 ; i < Nw ; i++ )
+	sum += A[i] ;
+
+    for ( i = 0 ; i < Nw ; i++ ) {
+     float afac = 2./sum ;
+     float sfac = Nw > N ? 1./afac : afac ;
+	A[i] *= afac ;
+	S[i] *= sfac ;
+    }
+
+    if ( Nw <= N && I ) {
+	for ( sum = i = 0 ; i < Nw ; i += I )
+	    sum += S[i]*S[i] ;
+	for ( sum = 1./sum, i = 0 ; i < Nw ; i++ )
+	    S[i] *= sum ;
+    }
+}
+
+void makehamming( float *H, float *A, float *S, int Nw, int N, int I, int osc,int  odd )
+{
+ int i;
+ float sum ;
+float PI, TWOPI;
+
+PI = 3.141592653589793115997963468544185161590576171875;
+TWOPI = 6.28318530717958623199592693708837032318115234375;
+
+/*
+ * basic Hamming windows
+ */
+ 
+ 
+ if (odd) {
+    for ( i = 0 ; i < Nw ; i++ )
+	  H[i] = A[i] = S[i] = sqrt(0.54 - 0.46*cos( TWOPI*i/(Nw - 1) ));
+ }
+	
+ else {
+
+   for ( i = 0 ; i < Nw ; i++ )
+	  H[i] = A[i] = S[i] = 0.54 - 0.46*cos( TWOPI*i/(Nw - 1) );
+
+ }
+ 	
+/*
+ * when Nw > N, also apply interpolating (sinc) windows to
+ * ensure that window are 0 at increments of N (the FFT length)
+ * away from the center of the analysis window and of I away
+ * from the center of the synthesis window
+ */
+    if ( Nw > N ) {
+     float x ;
+
+/*
+ * take care to create symmetrical windows
+ */
+	x = -(Nw - 1)/2. ;
+	for ( i = 0 ; i < Nw ; i++, x += 1. )
+	    if ( x != 0. ) {
+		A[i] *= N*sin( PI*x/N )/(PI*x) ;
+		if ( I )
+		    S[i] *= I*sin( PI*x/I )/(PI*x) ;
+	    }
+    }
+/*
+ * normalize windows for unity gain across unmodified
+ * analysis-synthesis procedure
+ */
+    for ( sum = i = 0 ; i < Nw ; i++ )
+	sum += A[i] ;
+
+    for ( i = 0 ; i < Nw ; i++ ) {
+     float afac = 2./sum ;
+     float sfac = Nw > N ? 1./afac : afac ;
+	A[i] *= afac ;
+	S[i] *= sfac ;
+    }
+
+    if ( Nw <= N && I ) {
+	for ( sum = i = 0 ; i < Nw ; i += I )
+	    sum += S[i]*S[i] ;
+	for ( sum = 1./sum, i = 0 ; i < Nw ; i++ )
+	    S[i] *= sum ;
+    }
+}
+
+
+void makehanning( float *H, float *A, float *S, int Nw, int N, int I, int osc, int odd )
+{
+ int i;
+ float sum ;
+float PI, TWOPI;
+
+PI = 3.141592653589793115997963468544185161590576171875;
+TWOPI = 6.28318530717958623199592693708837032318115234375;
+
+/*
+ * basic Hanning windows
+ */
+ 
+ 
+ if (odd) {
+    for ( i = 0 ; i < Nw ; i++ )
+	  H[i] = A[i] = S[i] = sqrt(0.5 * (1. + cos(PI + TWOPI * i / (Nw - 1))));
+ }
+	
+ else {
+
+   for ( i = 0 ; i < Nw ; i++ )
+	  H[i] = A[i] = S[i] = 0.5 * (1. + cos(PI + TWOPI * i / (Nw - 1)));
+
+ }
+ 	
+/*
+ * when Nw > N, also apply interpolating (sinc) windows to
+ * ensure that window are 0 at increments of N (the FFT length)
+ * away from the center of the analysis window and of I away
+ * from the center of the synthesis window
+ */
+    if ( Nw > N ) {
+     float x ;
+
+/*
+ * take care to create symmetrical windows
+ */
+	x = -(Nw - 1)/2. ;
+	for ( i = 0 ; i < Nw ; i++, x += 1. )
+	    if ( x != 0. ) {
+		A[i] *= N*sin( PI*x/N )/(PI*x) ;
+		if ( I )
+		    S[i] *= I*sin( PI*x/I )/(PI*x) ;
+	    }
+    }
+/*
+ * normalize windows for unity gain across unmodified
+ * analysis-synthesis procedure
+ */
+    for ( sum = i = 0 ; i < Nw ; i++ )
+	sum += A[i] ;
+
+    for ( i = 0 ; i < Nw ; i++ ) {
+     float afac = 2./sum ;
+     float sfac = Nw > N ? 1./afac : afac ;
+	A[i] *= afac ;
+	S[i] *= sfac ;
+    }
+
+    if ( Nw <= N && I ) {
+	for ( sum = i = 0 ; i < Nw ; i += I )
+	    sum += S[i]*S[i] ;
+	for ( sum = 1./sum, i = 0 ; i < Nw ; i++ )
+	    S[i] *= sum ;
+    }
+}
+
diff --git a/externals/grill/fftease/src/overlapadd.c b/externals/grill/fftease/src/overlapadd.c
new file mode 100644
index 00000000..7e832a00
--- /dev/null
+++ b/externals/grill/fftease/src/overlapadd.c
@@ -0,0 +1,17 @@
+/*
+ * input I is a folded spectrum of length N; output O and
+ * synthesis window W are of length Nw--overlap-add windowed,
+ * unrotated, unfolded input data into output O
+ */
+void overlapadd( float *I, int N, float *W, float *O, int Nw, int n )
+{
+ int i ;
+    while ( n < 0 )
+	n += N ;
+    n %= N ;
+    for ( i = 0 ; i < Nw ; i++ ) {
+	O[i] += I[n]*W[i] ;
+	if ( ++n == N )
+	    n = 0 ;
+    }
+}
diff --git a/externals/grill/fftease/src/pv.h b/externals/grill/fftease/src/pv.h
new file mode 100644
index 00000000..b945e72e
--- /dev/null
+++ b/externals/grill/fftease/src/pv.h
@@ -0,0 +1,60 @@
+/*
+ *  This file is part of the pv-lib
+ *
+ *  The pv-lib can be used by everyone as desired
+ *
+ *  (c) Eric Lyon and Christopher Penrose
+ */
+
+#include <stdio.h>
+
+// -------------------------------------
+// modifications by Thomas Grill
+
+//#include <fts/fts.h>
+#include <math.h>
+
+#ifdef _MSC_VER
+#pragma warning(disable: 4305)
+#pragma warning(disable: 4244)
+#pragma warning(disable: 4101)
+#endif
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+// -------------------------------------
+
+
+#define FORWARD 1
+#define INVERSE 0
+
+typedef struct {
+  float    min;
+  float    max;
+} Bound;
+
+void init_rdft(int n, int *ip, float *w);
+void rdft(int n, int isgn, float *a, int *ip, float *w);
+void fold( float *I, float *W, int Nw, float *O, int N, int n );
+void overlapadd(float *I, int N, float *W, float *O, int Nw, int n );
+void makehanning( float *H, float *A, float *S, int Nw, int N, int I, int osc, int odd );
+void makewindows( float *H, float *A, float *S, int Nw, int N, int I, int osc );
+void leanconvert( float *S, float *C, int N2 );
+void leanunconvert( float *C, float *S, int N2 );
+void rfft( float *x, int N, int forward );
+void cfft( float *x, int NC, int forward );
+void convert_new(float *S, float *C, int N2, float *lastphase,  float fundamental, float factor );
+void convert(float *S, float *C, int N2, float *lastphase,  float fundamental, float factor );
+void unconvert(float  *C, float *S, int N2, float *lastphase, float fundamental,  float factor );
+
+
+// -------------------------------------
+
+#ifdef __cplusplus
+}
+#endif
+
+// -------------------------------------
+
diff --git a/externals/grill/fftease/src/unconvert.c b/externals/grill/fftease/src/unconvert.c
new file mode 100644
index 00000000..54bb6ca1
--- /dev/null
+++ b/externals/grill/fftease/src/unconvert.c
@@ -0,0 +1,36 @@
+#include "pv.h"
+
+
+
+void unconvert(float  *C, float *S, int N2, float *lastphase, float fundamental,  float factor )
+{
+  int 		i,
+		real,
+		imag,
+		amp,
+		freq;
+  float 	mag,
+		phase;
+
+    for ( i = 0; i <= N2; i++ ) {
+
+	imag = freq = ( real = amp = i<<1 ) + 1;
+
+	if ( i == N2 )
+	  real = 1;
+
+	mag = C[amp];
+	lastphase[i] += C[freq] - i*fundamental;
+	phase = lastphase[i]*factor;
+	S[real] = mag*cos( phase );
+
+	if ( i != N2 )
+	  S[imag] = -mag*sin( phase );
+	/*
+      if( i == 10 ) {
+	fprintf(stderr,"unconvert: amp: %f freq: %f funda %f fac %f\n", C[amp],C[freq],fundamental,factor);
+      }
+	*/
+    }
+
+}