aboutsummaryrefslogtreecommitdiff
path: root/pvoc~.c
blob: c6c76ce0cb2087c0f7f7e2079397327e158f9d65 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
/* Copyright 2003 Benjamin R. Saylor <bensaylor@fastmail.fm>
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

#include <math.h>
#include <fftw3.h>
#include "m_pd.h"

// FIXME:
// set array when dsp is turned on
// get rid of shiftbuf, just save values that will be needed next before overwriting them
// cubic interp
// use float fftw?
// performance testing
// what if there are 2 transients less than fftsize apart?  second one might get smeared.
// compare sound with phaselockedvoc.pd
// detect transients
// peaks + noise
// other phase locking methods
// use floats?
// use in-place?
// if don't have an array, call setarray(x->arrayname)
// window size and fft size independent (what is gained by zero-padding?)
// error if parent blocksize is larger than hopsize
// slowly return window to true position after centering around a transient?
// use fewer fft arrays?

// DONE:
// use FFTW_MEASURE

static t_class *pvoc_class;

typedef struct _pvoc {
	t_object x_obj;
	t_symbol *arrayname;
	t_garray *arrayobj;
	t_float *array;
	int arraysize;
	double *window;
	int fftsize;
	int overlap;
	int hopsize;			// = fftsize / overlap
	int trans[256];			// sample indices of transients
	int ntrans;			// number of transients
	int wastrans;			// there was a transient in the left half of the window during the previous frame
	double phaselocking;
	fftw_plan fftplan;
	fftw_plan fft2plan;
	fftw_plan ifftplan;
	double *fftin;
	double *fft2in;
	double *ifftout;
	fftw_complex *fftout;
	fftw_complex *fft2out;
	fftw_complex *ifftin;
	fftw_complex *shiftbuf;
	double *outbuf;
	int outbufpos;
} t_pvoc;

// if there is a transient between samples a and b, return its position, else return -1
static inline int transient_between(t_pvoc *x, int a, int b)
{
	// linear search for now: FIXME
	int i;
	for (i = 0; i < x->ntrans; i++)
		if (a <= x->trans[i] && b >= x->trans[i])
			return x->trans[i];
	return -1;
}

#if 1
static inline double interpolate(t_pvoc *x, double t)
{
	// linear interpolation for now: FIXME
	if (t < 0 || t > (x->arraysize - 1))
		return 0.0;
	else {
		int x_1 = t;
		double y_1 = x->array[x_1];
		double y_2 = x->array[x_1 + 1];

		return (y_2 - y_1) * (t - x_1) + y_1;
	}
}
#else
static inline double interpolate(t_pvoc *x, double t)
{
	// FIXME check bounds (can't think now)
	int truncphase = (int) x->phase;
	double fr = x->phase - ((double) truncphase);
	double inm1 = x->ifftout[truncphase - 1];
	double in   = x->ifftout[truncphase + 0];
	double inp1 = x->ifftout[truncphase + 1];
	double inp2 = x->ifftout[truncphase + 2];

	// taken from swh-plugins-0.4.0/ladspa-util.h cube_interp, made to use doubles instead since doubles are what i'm using for some reason
	return in + 0.5 * fr * (inp1 - inm1 +
	 fr * (4.0 * inp1 + 2.0 * inm1 - 5.0 * in - inp2 +
	 fr * (3.0 * (in - inp1) - inm1 + inp2)));
}
#endif

static t_int *pvoc_perform(t_int *w)
{
	t_pvoc *x = (t_pvoc *)(w[1]);
	t_float *in1 = (t_float *)(w[2]);
	t_float *in2 = (t_float *)(w[3]);
	t_float *out = (t_float *)(w[4]);
	int n = (int)(w[5]);
	double t;
	double pitchshift;
	int transientpos;
	int desmear;
	double framestart;
	double frameend;
	int i;	// 0 to n -type iterator
	int j;	// start to end -type iterator
	int k;	// bin iterator
	double xlook;		// iterator for interpolated table lookup

	// if we are at the start of a new frame...
	if (x->outbufpos % x->hopsize == 0) {

		// don't desmear this frame by default
		desmear = 0;

		// sample the input signals (FIXME just sample these in the beginning..)
		t = in1[0];		// time position
		pitchshift = in2[0];	// pitch shift

		// set the frame boundaries with the desired time pos in the middle
		framestart = t - (pitchshift * x->fftsize / 2);
		frameend = framestart + pitchshift * x->fftsize;

		// prepare to de-smear transients
		transientpos = transient_between(x, (int) framestart, (int) frameend);
		if (transientpos != -1) {
			// there is a transient in this frame
#if 0
			if (transientpos > t) {
				// there is a transient in the right half of the window:
				// --> move the window left until the transient is outside it
				frameend = transientpos;
				framestart = frameend - x->fftsize;
				x->wastrans = 0;
			} else if ( ! x->wastrans) {
				// there is a transient in the left half of the window,
				// and there was no transient there during the previous frame:
				// --> center the window around the transient and remember to desmear this frame
				framestart = transientpos - (x->fftsize / 2);
				frameend = framestart + x->fftsize;
				desmear = 1;
				x->wastrans = 1;
			} else
				x->wastrans = 1;
#else
			// this simpler method turns out to sound better (timing sounds more accurate, no "frozen" sound preceding transients)
			if ( ! x->wastrans) {
				// there is a transient in the window,
				// and there wasn't during the previous frame:
				// --> center the window around the transient and remember to desmear this frame
				framestart = transientpos - (pitchshift * x->fftsize / 2);
				desmear = 1;
			}
			x->wastrans = 1;
#endif
		} else
			x->wastrans = 0;

		// interpolate-read the array from framestart to frameend into fftin, windowing it
		for (i = 0, xlook = framestart; i < x->fftsize; xlook += pitchshift, i++) {
			x->fftin[i] = interpolate(x, xlook) * x->window[i];
		}

		// hop forward and read the second frame into fft2in
		// FIXME merge the two loops?
		framestart += pitchshift * x->hopsize;
		for (i = 0, xlook = framestart; i < x->fftsize; xlook += pitchshift, i++) {
			x->fft2in[i] = interpolate(x, xlook) * x->window[i];
		}

		// do the ffts
		fftw_execute(x->fftplan);
		fftw_execute(x->fft2plan);

		if ( ! desmear) {
			// Miller Puckette's phase modification math (translation from 09.pvoc.pd and 10.phaselockedvoc.pd)

			double a, b, r, c, d;

			// propagate phase
			for (k = 0; k < (x->fftsize / 2 + 1); k++) {
				a = x->ifftin[k][0] * x->fftout[k][0] + x->ifftin[k][1] * x->fftout[k][1] + 0.00000000000000000001;
				b = x->ifftin[k][1] * x->fftout[k][0] - x->ifftin[k][0] * x->fftout[k][1];
				r = 1 / sqrt(a * a + b * b);
				c = a * r;
				d = b * r;
				x->shiftbuf[k][0] = c * x->fft2out[k][0] - d * x->fft2out[k][1];
				x->shiftbuf[k][1] = c * x->fft2out[k][1] + d * x->fft2out[k][0];
			}

			// don't phase-lock the first bin
			x->ifftin[0][0] = x->shiftbuf[0][0];
			x->ifftin[0][1] = x->shiftbuf[0][1];

			// phase-lock
			for (k = 1; k < (x->fftsize / 2); k++) {
				x->ifftin[k][0] = x->shiftbuf[k][0] - x->phaselocking * (x->shiftbuf[k - 1][0] + x->shiftbuf[k + 1][0]);
				x->ifftin[k][1] = x->shiftbuf[k][1] - x->phaselocking * (x->shiftbuf[k - 1][1] + x->shiftbuf[k + 1][1]);
			}

			// don't phase-lock the last bin
			x->ifftin[x->fftsize / 2][0] = x->shiftbuf[x->fftsize / 2][0];
			x->ifftin[x->fftsize / 2][1] = x->shiftbuf[x->fftsize / 2][1];

		} else {
			// this frame is to be de-smeared, which means don't modify the phases, just preserve the original phases
			for (k = 0; k < (x->fftsize / 2 + 1); k++) {
				x->ifftin[k][0] = x->fftout[k][0];
				x->ifftin[k][1] = x->fftout[k][1];
			}
		}

		// do the ifft
		fftw_execute(x->ifftplan);

		// add into output buffer, windowing and normalizing first (divide by blocksize)
		for (i = 0, j = x->outbufpos; i < x->fftsize; i++, j++) {
			x->outbuf[j % x->fftsize] += x->ifftout[i] / x->fftsize * x->window[i];
		}
	}

	// output one block of the output buffer
	for (i = 0, j = x->outbufpos; i < n; i++, j++) {
		out[i] = x->outbuf[j % x->fftsize];
		x->outbuf[j % x->fftsize] = 0;		// zero the part of the buffer that was just output
	}

	// move the output buffer pointer forward by one block
	x->outbufpos = (x->outbufpos + n) % x->fftsize;

	return (w+6);
}

static void pvoc_dsp(t_pvoc *x, t_signal **sp)
{
	    dsp_add(pvoc_perform, 5, x, sp[0]->s_vec, sp[1]->s_vec, sp[2]->s_vec, sp[0]->s_n);
}

// adapted from jsarlo's windowing library
// Hanning
static void makewindow(double *w, int n)
{
	int i;
	double xshift =  n / 2.0;
	double x;
	for (i = 0; i < n; i++) {
		x = (i - xshift) / xshift;
		w[i] = 0.5 * (1 + cos(M_PI * x));
	}
}

static void setarray(t_pvoc *x, t_symbol *s)
{
	x->arrayname = s;
	if ( ! (x->arrayobj = (t_garray *)pd_findbyclass(x->arrayname, garray_class))) {
 		if (*x->arrayname->s_name) pd_error(x, "pvoc~: %s: no such array", x->arrayname->s_name);
		x->array = NULL;
		x->arraysize = 0;
	} else if ( ! garray_getfloatarray(x->arrayobj, &x->arraysize, &x->array)) {
 		error("%s: bad template", x->arrayname->s_name);
		x->array = NULL;
		x->arraysize = 0;
	} else {
		garray_usedindsp(x->arrayobj);
	}
}

static void locking(t_pvoc *x, t_floatarg f)
{
	x->phaselocking = f;
}

// takes a list of sample positions of transients to be de-smeared
static void transients(t_pvoc *x, t_symbol *s, int argc, t_atom *argv)
{
	int i;

	x->ntrans = argc;
	for (i = 0; i < x->ntrans; i++)
		x->trans[i] = atom_getfloatarg(i, argc, argv);
}

// for clarity (same as "transients" with no args)
static void notransients(t_pvoc *x)
{
	x->ntrans = 0;
}

static void *pvoc_new(t_symbol *s, int argc, t_atom *argv)
{
	t_pvoc *x = (t_pvoc *)pd_new(pvoc_class);
	int i;

	inlet_new(&x->x_obj, &x->x_obj.ob_pd, &s_signal, &s_signal);	// pitch-shift inlet
	outlet_new(&x->x_obj, gensym("signal"));

	if (argc != 3) {
		post("argc = %d", argc);
		error("pvoc~: usage: [pvoc~ <arrayname> <fftsize> <overlap>]");
		return NULL;
	}

	x->fftsize = atom_getfloatarg(1, argc, argv);
	x->overlap = atom_getfloatarg(2, argc, argv);
	x->hopsize = x->fftsize / x->overlap;
	x->ntrans = 0;
	x->wastrans = 0;
	x->phaselocking = 0;

	// get the source array
	setarray(x, atom_getsymbol(argv));

	// set up output ring buffer
	x->outbuf = getbytes(sizeof(double) * x->fftsize);
	x->outbufpos = 0;
	for (i = 0; i < x->fftsize; i++)
		x->outbuf[i] = 0;

	// make table for window function
	x->window = getbytes(sizeof(double) * x->fftsize);
	makewindow(x->window, x->fftsize);

	// set up fftw stuff
	x->fftin = fftw_malloc(sizeof(double) * x->fftsize);
	x->fft2in = fftw_malloc(sizeof(double) * x->fftsize);
	x->ifftout = fftw_malloc(sizeof(double) * x->fftsize);
	x->fftout = fftw_malloc(sizeof(fftw_complex) * (x->fftsize / 2 + 1));
	x->fft2out = fftw_malloc(sizeof(fftw_complex) * (x->fftsize / 2 + 1));
	x->ifftin = fftw_malloc(sizeof(fftw_complex) * (x->fftsize / 2 + 1));
	x->shiftbuf = fftw_malloc(sizeof(fftw_complex) * (x->fftsize / 2 + 1));
	for (i = 0; i < (x->fftsize / 2 + 1); i++) {
		x->ifftin[i][0] = 0;	// need to start the phases from zero
		x->ifftin[i][1] = 0;
	}
	x->fftplan = fftw_plan_dft_r2c_1d(x->fftsize, x->fftin, x->fftout, FFTW_MEASURE);
	x->fft2plan = fftw_plan_dft_r2c_1d(x->fftsize, x->fft2in, x->fft2out, FFTW_MEASURE);
	x->ifftplan = fftw_plan_dft_c2r_1d(x->fftsize, x->ifftin, x->ifftout, FFTW_MEASURE | FFTW_PRESERVE_INPUT);
	
	return (x);
}

static void pvoc_free(t_pvoc *x)
{
	freebytes(x->outbuf, sizeof(double) * x->fftsize);
	freebytes(x->window, sizeof(double) * x->fftsize);
	fftw_free(x->fftin);
	fftw_free(x->fft2in);
	fftw_free(x->ifftout);
	fftw_free(x->fftout);
	fftw_free(x->fft2out);
	fftw_free(x->ifftin);
	fftw_free(x->shiftbuf);
	fftw_destroy_plan(x->fftplan);
	fftw_destroy_plan(x->fft2plan);
	fftw_destroy_plan(x->ifftplan);
}

void pvoc_tilde_setup(void)
{
	pvoc_class = class_new(gensym("pvoc~"), (t_newmethod)pvoc_new, (t_method)pvoc_free, sizeof(t_pvoc), 0, A_GIMME, 0);
	class_addmethod(pvoc_class, nullfn, gensym("signal"), 0);
	class_addmethod(pvoc_class, (t_method) pvoc_dsp, gensym("dsp"), 0);
	class_addmethod(pvoc_class, (t_method) setarray, gensym("setarray"), A_DEFSYMBOL, 0);
	class_addmethod(pvoc_class, (t_method) locking, gensym("locking"), A_DEFFLOAT, 0);
	class_addmethod(pvoc_class, (t_method) transients, gensym("transients"), A_GIMME, 0);
	class_addmethod(pvoc_class, (t_method) notransients, gensym("notransients"), 0);
}