History

implemented basic DSP for Spectacle

Joel Matthys committed 12 years ago

Commit 4325a601 parent 0ecc72732 files · +3252 −1434Browse files

Showing the changes to Spectacle/SpectacleBase.cpp (1 of 32 files). Show all

Unified · Split

Spectacle/SpectacleBase.cpp+646−0
@@ -0,0 +1,646 @@
1/*
2 * Copyright (C) 2007 John Gibson
3 *
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License, version 2, as
6 * published by the Free Software Foundation.
7 *
8 * This program is distributed in the hope that it will be useful,
9 * but WITHOUT ANY WARRANTY; without even the implied warranty of
10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 * GNU General Public License for more details.
12 */
13
14// Base class for the various spectacle externals, derived from
15// my SPECTACLE2 RTcmix instruments.
16// John Gibson <johgibso at indiana dot edu>, 11/23/07.
17
18//#define DEBUG
19//#define DEBUG_BINGROUPS
20//#define CHECK_BINGROUPS
21
22#include "SpectacleBase.h"
23#include <float.h>
24#include <cstddef>
25
26#ifndef cosf
27#define cosf(x) cos((x))
28#endif
29#ifndef sinf
30#define sinf(x) sin((x))
31#endif
32#ifndef hypotf
33#define hypotf(x, y) hypot((x), (y))
34#endif
35#ifndef atan2f
36#define atan2f(y, x) atan2((y), (x))
37#endif
38
39const int kMaxWindowLen = kMaxFFTLen * 8;
40const int kMinOverlap = 1;
41const int kMaxOverlap = 64;
42
43
44// ------------------------------------------------------------ SpectacleBase --
45SpectacleBase::SpectacleBase() :
46 _control_table_size(0),
47 _bin_groups(NULL),
48 _binmaptable(NULL),
49 _print_stats(false),
50 _reading_input(true),
51 _posteq(false),
52 _prev_bg_ignorevals(-1),
53 _srate(-FLT_MAX),
54 _minfreq(-FLT_MAX), _maxfreq(-FLT_MAX),
55 _anal_window(NULL), _synth_window(NULL),
56 _input(NULL), _output(NULL),
57 _outbuf(NULL),
58 _cursamp(0L),
59 _fft(NULL),
60 _bucket(NULL)
61{
62}
63
64
65// ----------------------------------------------------------- ~SpectacleBase --
66SpectacleBase::~SpectacleBase()
67{
68 delete [] _outbuf;
69 delete [] _input;
70 delete [] _output;
71 delete [] _anal_window;
72 delete [] _synth_window;
73 delete [] _bin_groups;
74 delete _bucket;
75 delete _fft;
76}
77
78
79// ------------------------------------------------------------- make_windows --
80int SpectacleBase::make_windows()
81{
82 // Hamming window
83 for (int i = 0; i < _window_len; i++)
84 _anal_window[i] = _synth_window[i] = 0.54f - 0.46f
85 * cosf(2.0f * M_PI * i / (_window_len - 1));
86
87 // When _window_len > _fftlen, also apply interpolating (sinc) windows to
88 // ensure that window is 0 at increments of _fftlen away from the center
89 // of the analysis window and of decimation away from the center of the
90 // synthesis window.
91
92 if (_window_len > _fftlen) {
93 float x = -(_window_len - 1) / 2.0;
94 for (int i = 0; i < _window_len; i++, x += 1.0f)
95 if (x != 0.0f) {
96 _anal_window[i] *= _fftlen * sin(M_PI * x / _fftlen) / (M_PI * x);
97 if (_decimation)
98 _synth_window[i] *= _decimation * sin(M_PI * x / _decimation)
99 / (M_PI * x);
100 }
101 }
102
103 // Normalize windows for unity gain across unmodified
104 // analysis-synthesis procedure.
105
106 float sum = 0.0f;
107 for (int i = 0; i < _window_len; i++)
108 sum += _anal_window[i];
109
110 for (int i = 0; i < _window_len; i++) {
111 float afac = 2.0f / sum;
112 float sfac = _window_len > _fftlen ? 1.0f / afac : afac;
113 _anal_window[i] *= afac;
114 _synth_window[i] *= sfac;
115 }
116
117 if (_window_len <= _fftlen && _decimation) {
118 sum = 0.0f;
119 for (int i = 0; i < _window_len; i += _decimation)
120 sum += _synth_window[i] * _synth_window[i];
121 sum = 1.0f / sum;
122 for (int i = 0; i < _window_len; i++)
123 _synth_window[i] *= sum;
124 }
125
126 return 0;
127}
128
129
130// ---------------------------------------------------------------- set_srate --
131// Must call this before anything that updates bin groups in subclass.
132void SpectacleBase::set_srate(float srate)
133{
134 if (srate != _srate) {
135 _srate = srate;
136 if (_maxfreq == _nyquist)
137 _maxfreq = 0.0f;
138 _nyquist = srate / 2;
139 if (_maxfreq == 0.0f)
140 _maxfreq = _nyquist;
141 _fund_anal_freq = srate / float(_fftlen);
142 update_bin_groups(_bin_groups, NULL, _minfreq, _maxfreq,
143 _control_table_size);
144 }
145}
146
147
148// --------------------------------------------------------------------- init --
149int SpectacleBase::init(int fftlen, int windowlen, int overlap, float srate)
150{
151 _fftlen = fftlen;
152 _window_len = windowlen;
153 _overlap = overlap;
154
155 // Make sure FFT length is a power of 2 <= kMaxFFTLen.
156 bool valid = false;
157 for (int x = 1; x <= kMaxFFTLen; x *= 2) {
158 if (_fftlen == x) {
159 valid = true;
160 break;
161 }
162 }
163 if (!valid) {
164 //post("%s: FFT length must be a power of two <= %d. Setting to 1024...",
165 //instname(), kMaxFFTLen);
166 _fftlen = 1024;
167 }
168
169 _half_fftlen = _fftlen / 2;
170
171 // Make sure window length is a power of 2 >= FFT length.
172 valid = false;
173 for (int x = _fftlen; x <= kMaxWindowLen; x *= 2) {
174 if (_window_len == x) {
175 valid = true;
176 break;
177 }
178 }
179 if (!valid) {
180 //post("%s: Window length must be a power of two >= FFT length (%d)\n"
181 //"and <= %d. Setting to 2048...",
182 //instname(), _fftlen, kMaxWindowLen);
183 _window_len = _fftlen * 2;
184 }
185
186 // Make sure _overlap is a power of 2 in allowed range.
187 valid = false;
188 for (int x = kMinOverlap; x <= kMaxOverlap; x *= 2) {
189 if (_overlap == x) {
190 valid = true;
191 break;
192 }
193 }
194 if (!valid) {
195 //post("%s: Overlap must be a power of two between %d and %d. "
196 //"Setting to 2...",
197 //instname(), kMinOverlap, kMaxOverlap);
198 _overlap = 2;
199 }
200
201 _bin_groups = new int [_half_fftlen];
202
203 set_srate(srate);
204 set_freqrange(0.0f, 0.0f);
205
206 // derive decimation from overlap
207 _decimation = int(_fftlen / _overlap);
208
209 _window_len_minus_decimation = _window_len - _decimation;
210
211 DPRINT2("_fftlen=%d, _decimation=%d", _fftlen, _decimation);
212
213 _input = new float [_window_len]; // interior input buffer
214 _output = new float [_window_len]; // interior output buffer
215 if (_input == NULL || _output == NULL)
216 return -1;
217 for (int i = 0; i < _window_len; i++)
218 _input[i] = _output[i] = 0.0f;
219
220 // Read index chases write index by _decimation; add 2 extra locations to
221 // keep read point from stepping on write point. Verify with asserts in
222 // increment_out_*_index().
223 _outframes = _decimation + 2;
224 _out_read_index = _outframes - _decimation;
225 _out_write_index = 0;
226 _outbuf = new float [_outframes];
227 if (_outbuf == NULL)
228 return -1;
229 for (int i = 0; i < _outframes; i++)
230 _outbuf[i] = 0.0f;
231 DPRINT1("_outframes: %d", _outframes);
232
233 _anal_window = new float [_window_len];
234 _synth_window = new float [_window_len];
235 if (_anal_window == NULL || _synth_window == NULL)
236 return -1;
237 if (make_windows() != 0)
238 return -1;
239
240 _bucket = new Obucket(_decimation, process_wrapper, (void *) this);
241
242 _fft = new Offt(_fftlen);
243 _fft_buf = _fft->getbuf();
244
245 return 0;
246}
247
248
249// -------------------------------------------------------------------- clear --
250void SpectacleBase::clear()
251{
252 for (int i = 0; i < _fftlen; i++)
253 _fft_buf[i] = 0.0f;
254 for (int i = 0; i < _window_len; i++)
255 _input[i] = _output[i] = 0.0f;
256 for (int i = 0; i < _outframes; i++)
257 _outbuf[i] = 0.0f;
258 _bucket->clear();
259}
260
261
262// ------------------------------------------------------------ prepare_input --
263// <buf> contains <_decimation> samples (mono) from most recent input.
264
265void SpectacleBase::prepare_input(const float buf[])
266{
267 DPRINT("prepare_input");
268
269 // Shift samples in <_input> from right to left by <_decimation>,
270 // leaving a hole of <_decimation> slots at right end.
271
272 for (int i = 0; i < _window_len_minus_decimation; i++)
273 _input[i] = _input[i + _decimation];
274
275 // Copy <_decimation> samples from <buf> to right end of <_input>.
276
277 for (int i = _window_len_minus_decimation, j = 0; i < _window_len; i++, j++)
278 _input[i] = buf[j];
279
280 // Multiply input array by analysis window, both of length <_window_len>.
281 // Fold and rotate windowed real input into FFT buffer.
282
283 for (int i = 0; i < _fftlen; i++)
284 _fft_buf[i] = 0.0f;
285 int j = _cursamp % _fftlen;
286 for (int i = 0; i < _window_len; i++) {
287 _fft_buf[j] += _input[i] * _anal_window[i];
288 if (++j == _fftlen)
289 j = 0;
290 }
291}
292
293
294// ----------------------------------------------------------- prepare_output --
295void SpectacleBase::prepare_output()
296{
297 DPRINT("prepare_output");
298
299 // overlap-add <_fft_buf> real data into <_output>
300 int j = _cursamp % _fftlen;
301 for (int i = 0; i < _window_len; i++) {
302 _output[i] += _fft_buf[j] * _synth_window[i];
303 if (++j == _fftlen)
304 j = 0;
305 }
306
307 // transfer samples from <_output> to outer output buffer.
308 for (int i = 0; i < _decimation; i++) {
309 _outbuf[_out_write_index] = _output[i];
310 increment_out_write_index();
311 }
312
313 // shift samples in <_output> from right to left by <_decimation>
314 for (int i = 0; i < _window_len_minus_decimation; i++)
315 _output[i] = _output[i + _decimation];
316 for (int i = _window_len_minus_decimation; i < _window_len; i++)
317 _output[i] = 0.0f;
318}
319
320
321// ---------------------------------------------------------- process_wrapper --
322// Called by Obucket whenever the input bucket is full (i.e., has _decimation
323// samps). This static member wrapper function lets us call a non-static member
324// function, which we can't pass directly as a callback to Obucket.
325// See http://www.newty.de/fpt/callback.html for one explanation of this.
326
327void SpectacleBase::process_wrapper(const float buf[], const int len,
328 void *obj)
329{
330 SpectacleBase *myself = (SpectacleBase *) obj;
331 myself->process(buf, len);
332}
333
334
335// ------------------------------------------------------------------ process --
336void SpectacleBase::process(const float *buf, const int)
337{
338 DPRINT("SpectacleBase::process");
339
340 if (_reading_input) {
341 prepare_input(buf);
342 _fft->r2c();
343 }
344 modify_analysis(_reading_input);
345 _fft->c2r();
346 prepare_output();
347}
348
349
350// ------------------------------------------------------------ set_freqrange --
351// Update the frequency range within which control tables operate. Return true
352// if range has actually changed, resulting in an update_bin_groups call.
353// Otherwise, return false.
354
355bool SpectacleBase::set_freqrange(float min, float max)
356{
357 if (min != _minfreq || max != _maxfreq) {
358 if (max == 0.0f)
359 max = _nyquist;
360 if (max < min)
361 _fswap(&min, &max);
362 if (max > _nyquist)
363 max = _nyquist;
364 _maxfreq = max;
365 _minfreq = _fclamp(0.0f, min, max);
366 update_bin_groups(_bin_groups, NULL, _minfreq, max, _control_table_size);
367 return true;
368 }
369 return false;
370}
371
372
373// -------------------------------------------------------- update_bin_groups --
374// Update a bin groups array. (See .h for a description of bin groups.) Note
375// that this is called from set_freqrange and any subclass equivalents.
376// Three methods of updating bin groups array: (1) use a binmap table,
377// (2) use a one-to-one mapping of bin groups to bins, or (3) use a non-linear
378// mapping that starts as one-to-one, but combines bins in ever larger
379// groups as bin frequency increases. The <binmaptable> argument is an
380// optional binmap table to use in place of the base class <_binmaptable>.
381// If not overriding the base class table, pass NULL.
382
383void SpectacleBase::update_bin_groups(
384 int bin_groups[],
385 int binmaptable[],
386 const float minfreq,
387 const float maxfreq,
388 const int control_table_size)
389{
390 if (control_table_size == 0) // no control table
391 return;
392
393 const int nbins = _half_fftlen;
394 const int *bmtable = binmaptable ? binmaptable : _binmaptable;
395
396 // The optional <bmtable> has <control_table_size> elements, each giving
397 // the number of adjacent FFT bins controlled by that element in any of the
398 // control tables (e.g., EQ). Copy this information into the <bin_groups>
399 // array, which is sized to <nbins>. We ignore minfreq and maxfreq.
400
401 if (bmtable != NULL) {
402#ifdef DEBUG_BINGROUPS
403 //post("using binmap table -----------------------------");
404 for (int i = 0; i < control_table_size; i++)
405 //post("[%d] %d", i, bmtable[i]);
406#endif
407 int bidx = 0;
408 for (int i = 0; i < control_table_size; i++) {
409 int bincount = bmtable[i];
410 while (bincount-- > 0) {
411 if (bidx == nbins)
412 goto highcount;
413 bin_groups[bidx++] = i;
414 }
415 }
416 highcount:
417 if (bidx < nbins) {
418 const int last = control_table_size - 1;
419 while (bidx < nbins)
420 bin_groups[bidx++] = last;
421 }
422 }
423
424 // If there is no binmap table...
425 // If freq. range is full bandwidth, and control table size is >= to the
426 // number of bins, then use linear mapping of control table slots to bins,
427 // as in SPECTACLE v1. If table is larger than the number of bins, warn
428 // about ignoring the extra table slots.
429
430 else if (control_table_size >= nbins
431 && minfreq == 0 && maxfreq == _nyquist) {
432 for (int i = 0; i < nbins; i++)
433 bin_groups[i] = i;
434
435 if (control_table_size > nbins) {
436 const int ignorevals = control_table_size - nbins;
437 if (ignorevals != _prev_bg_ignorevals) {
438 //post("%s: Control table of size %d too large for "
439 //"frequency range...ignoring last %d values",
440 //instname(), control_table_size, ignorevals);
441 _prev_bg_ignorevals = ignorevals;
442 }
443 }
444 }
445
446 // Otherwise, there is either a freq. range that is not full bandwidth, or
447 // the control table size is less than the number of bins. If the former,
448 // we insure that the first array slot affects all bins whose frequencies
449 // are below (and possibly a little above) the min. freq., and that the last
450 // array slot affects all bins whose frequencies are above (and possibly a
451 // little below) the max. freq. If the control table is smaller than the
452 // number of bins within the freq. range, then we construct a one-to-one or
453 // one-to-many mapping of control table slots to bins, arranged so that we
454 // can control lower frequencies with greater resolution.
455 //
456 // Depending on the difference between the control table size and the number
457 // of FFT bins, we create a mapping that begins linearly and grows
458 // arithmetically after a certain point. The purpose is to get higher
459 // resolution in the lower frequencies of the range, where it matters most.
460 // So near the low end of the range, there is a one-to-one mapping of control
461 // array slots to FFT bins. Near the high end of the range, one control
462 // array slot affects many FFT bins. We attempt to transition smoothly
463 // between these two extremes. As an example, the number of FFT bins
464 // controlled by each slot of a control table array might look like this:
465 //
466 // number of bins spanned by each cntl slot
467 // (lowshelf) 1 1 1 1 1 1 1 1 1 1 1 2 3 4 5 6 7 8 9 10 14 (highshelf)
468 //
469 // The scheme isn't perfect -- e.g., the penultimate slot might control more
470 // than you would expect -- but it does guarantee that each control table
471 // slot affects at least one FFT bin.
472
473 else {
474 const int lowshelfbin = closest_bin(minfreq);
475 int highshelfbin = closest_bin(maxfreq);
476 bool maxfreq_is_nyquist = false;
477 if (highshelfbin == nbins) {
478 maxfreq_is_nyquist = true;
479 highshelfbin--; // we don't control the nyquist bin
480 }
481 int cntltablen = control_table_size;
482 const float endflatrange = minfreq + (cntltablen * _fund_anal_freq);
483 bool linear_map = true;
484 if (endflatrange > maxfreq - _fund_anal_freq) {
485 const int ignorevals = int((endflatrange - maxfreq)
486 / _fund_anal_freq);
487 if (ignorevals > 0) {
488 cntltablen = control_table_size - ignorevals;
489 if (ignorevals != _prev_bg_ignorevals) {
490 //post("%s: Control table of size %d too large for "
491 //"frequency range...ignoring last %d values",
492 //instname(), control_table_size, ignorevals);
493 _prev_bg_ignorevals = ignorevals;
494 }
495 }
496 }
497 else
498 linear_map = false;
499#ifdef DEBUG_BINGROUPS
500 //post("lowbin=%d, highbin=%d, endflatrange=%f, map=%s",
501 //lowshelfbin, highshelfbin, endflatrange,
502 //linear_map ? "linear" : "nonlinear");
503#endif
504
505 // assign low shelf group
506 int bidx = 0;
507 while (bidx <= lowshelfbin)
508 bin_groups[bidx++] = 0;
509
510 // assign interior groups
511 if (linear_map) {
512 const int lastbingroup = maxfreq_is_nyquist ? cntltablen - 1
513 : cntltablen - 2;
514 int bingroup = 1;
515 while (bidx < highshelfbin) {
516 bin_groups[bidx++] = bingroup;
517 if (bingroup < lastbingroup)
518 bingroup++;
519 }
520 }
521 else { // higher array slots control more and more bins
522 const int extrabins = (highshelfbin - lowshelfbin) - cntltablen;
523 //if (extrabins <= 0)
524 //error("%s: program error - contact johgibso@gmail.com", instname());
525 const int nsums = int(sqrt(2 * extrabins) + 2);
526
527 // Compute the sum of integers for n in [nsums-2, nsums);
528 // formula is the closed-form equation: x = n * (n + 1) / 2.
529 const int an = nsums - 2;
530 const int bn = nsums - 1;
531 const int a = (an * (an + 1)) / 2;
532 const int b = (bn * (bn + 1)) / 2;
533 int cntlspan = 0, binspan = 0;
534 if (a > extrabins) {
535 cntlspan = an;
536 binspan = a;
537 }
538 else if (b > extrabins) {
539 cntlspan = bn;
540 binspan = b;
541 }
542 //else // this should never happen
543 //error("%s: program error - contact johgibso@gmail.com", instname());
544
545#ifdef DEBUG_BINGROUPS
546 //post("nbins=%d, lowbin=%d, highbin=%d, tablen=%d, extrabins=%d, "
547 //"cntlspan=%d, binspan=%d, maxfreqisnyq=%d\n",
548 //nbins, lowshelfbin, highshelfbin, cntltablen, extrabins,
549 //cntlspan, binspan, maxfreq_is_nyquist);
550#endif
551
552 // first, linear mapping
553 const int end = (cntltablen - cntlspan) - 1;
554 int bingroup = 1;
555 while (bingroup < end)
556 bin_groups[bidx++] = bingroup++;
557
558 // then map using sum-of-integers series
559 const int lastbingroup = maxfreq_is_nyquist ? cntltablen - 1
560 : cntltablen - 2;
561 int incr = 1;
562 int count = 0;
563 while (bidx < highshelfbin) {
564 bin_groups[bidx++] = bingroup;
565 if (++count == incr) {
566 count = 0;
567 if (bingroup < lastbingroup)
568 bingroup++;
569 incr++;
570 }
571 }
572 }
573
574 // assign high shelf group
575 const int last = cntltablen - 1;
576 while (bidx < nbins)
577 bin_groups[bidx++] = last;
578}
579
580#ifdef DEBUG_BINGROUPS
581//post("bin_groups[] -----------------------------");
582for (int i = 0; i < nbins; i++)
583 //post("[%d] %d (%f)", i, bin_groups[i], i * _fund_anal_freq);
584#endif
585#ifdef CHECK_BINGROUPS
586 for (int i = 1; i < nbins; i++) {
587 int diff = bin_groups[i] - bin_groups[i - 1];
588 if (diff < 0 || diff > 1)
589 //error("%s: bin group (%p) index %d not 0 or 1 greater than prev entry",
590 instname(), bin_groups, i);
591}
592#endif
593}
594
595
596// --------------------------------------------------------- print_bin_groups --
597// print a user-readable listing of bin groups
598void SpectacleBase::print_bin_groups(
599 int bin_groups[],
600 const char *type) // to identify bin groups to user
601{
602 const int nbins = _half_fftlen;
603
604 //post("\n%s: %s table bin groups", instname(), type);
605 //post("------------------------------------------");
606 int cntltabslot = 0;
607 int startbin = 0;
608 for (int i = 0; i < nbins; i++) {
609 const int thisslot = bin_groups[i];
610 if (thisslot > cntltabslot) {
611 // print stats for previous control table slot <cntltabslot>
612 const int startfreq = int(_fund_anal_freq * startbin + 0.5f);
613 if (i - startbin > 1) {
614 const int endfreq = int(_fund_anal_freq * (i - 1) + 0.5f);
615 //post(" [%d]\t%d-%d Hz (%d bins)",
616 //cntltabslot, startfreq, endfreq, i - startbin);
617 }
618 else
619 //post(" [%d]\t%d Hz", cntltabslot, startfreq);
620 cntltabslot = thisslot;
621 startbin = i;
622 }
623}
624}
625
626
627// ---------------------------------------------------------------------- run --
628void SpectacleBase::run(float in[], float out[], int nframes)
629{
630 if (_reading_input) {
631 for (int i = 0; i < nframes; i++) {
632 _bucket->drop(in[i]); // may process <_decimation> input frames
633 out[i] = _outbuf[_out_read_index];
634 increment_out_read_index();
635 _cursamp++;
636 }
637 }
638 else {
639 for (int i = 0; i < nframes; i++) {
640 _bucket->drop(0.0f); // may process <_decimation> input frames
641 out[i] = _outbuf[_out_read_index];
642 increment_out_read_index();
643 _cursamp++;
644 }
645 }
646}

Keyboard shortcuts

General

Show keyboard shortcuts
?
Focus the current search
/

Issues and pull requests

Open filters
f

Issues, pull requests and Inbox

Select next item
j↓
Select previous item
k↑
Open selected item
Enter
Dismiss filters or clear selection
Esc

Repository code

Go to file
t