From 0182bbff2871114a4e93cc97942da621491f0e02 Mon Sep 17 00:00:00 2001 From: Thomas Grill Date: Tue, 7 Jan 2003 00:28:39 +0000 Subject: "" svn path=/trunk/; revision=325 --- externals/grill/fftease/src/burrow~.cpp | 255 +++++++++++++++++++++ externals/grill/fftease/src/convert.c | 55 +++++ externals/grill/fftease/src/convert_new.c | 64 ++++++ externals/grill/fftease/src/fft.c | 161 +++++++++++++ externals/grill/fftease/src/fft4.c | 340 ++++++++++++++++++++++++++++ externals/grill/fftease/src/fold.c | 24 ++ externals/grill/fftease/src/leanconvert.c | 19 ++ externals/grill/fftease/src/leanunconvert.c | 23 ++ externals/grill/fftease/src/main.cpp | 60 +++++ externals/grill/fftease/src/main.h | 40 ++++ externals/grill/fftease/src/makewindows.c | 202 +++++++++++++++++ externals/grill/fftease/src/overlapadd.c | 17 ++ externals/grill/fftease/src/pv.h | 60 +++++ externals/grill/fftease/src/unconvert.c | 36 +++ 14 files changed, 1356 insertions(+) create mode 100644 externals/grill/fftease/src/burrow~.cpp create mode 100644 externals/grill/fftease/src/convert.c create mode 100644 externals/grill/fftease/src/convert_new.c create mode 100644 externals/grill/fftease/src/fft.c create mode 100644 externals/grill/fftease/src/fft4.c create mode 100644 externals/grill/fftease/src/fold.c create mode 100644 externals/grill/fftease/src/leanconvert.c create mode 100644 externals/grill/fftease/src/leanunconvert.c create mode 100644 externals/grill/fftease/src/main.cpp create mode 100644 externals/grill/fftease/src/main.h create mode 100644 externals/grill/fftease/src/makewindows.c create mode 100644 externals/grill/fftease/src/overlapadd.c create mode 100644 externals/grill/fftease/src/pv.h create mode 100644 externals/grill/fftease/src/unconvert.c (limited to 'externals/grill/fftease/src') 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 +#include + + + +/* 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 +#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 + +#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 + +// ------------------------------------- +// modifications by Thomas Grill + +//#include +#include + +#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); + } + */ + } + +} -- cgit v1.2.1