BogaudioModules

BogaudioModules for VCV Rack
Log | Files | Refs | README | LICENSE

analyzer.cpp (5610B)


      1 
      2 #include "ffft/FFTRealFixLen.h"
      3 
      4 #include "buffer.hpp"
      5 #include "analyzer.hpp"
      6 
      7 using namespace bogaudio::dsp;
      8 
      9 void Window::apply(float* in, float* out) {
     10 	for (int i = 0; i < _size; ++i) {
     11 		out[i] = in[i] * _window[i];
     12 	}
     13 }
     14 
     15 
     16 HanningWindow::HanningWindow(int size, float alpha) : Window(size) {
     17 	const float invAlpha = 1.0 - alpha;
     18 	const float twoPIEtc = 2.0 * M_PI / (float)size;
     19 	for (int i = 0; i < _size; ++i) {
     20 		_sum += _window[i] = invAlpha*cos(twoPIEtc*(float)i + M_PI) + alpha;
     21 	}
     22 }
     23 
     24 
     25 KaiserWindow::KaiserWindow(int size, float alpha) : Window(size) {
     26 	float ii0a = 1.0f / i0(alpha);
     27 	float ism1 = 1.0f / (float)(size - 1);
     28 	for (int i = 0, n = _size - 1; i < n; ++i) {
     29 		float x = i * 2.0f;
     30 		x *= ism1;
     31 		x -= 1.0f;
     32 		x *= x;
     33 		x = 1.0f - x;
     34 		x = sqrtf(x);
     35 		x *= alpha;
     36 		_sum += _window[i] = i0(x) * ii0a;
     37 	}
     38 	_window[_size - 1] = 0;
     39 }
     40 
     41 // Rabiner, Gold: "The Theory and Application of Digital Signal Processing", 1975, page 103.
     42 float KaiserWindow::i0(float x) {
     43 	assert(x >= 0.0f);
     44 	assert(x < 20.0f);
     45 	float y = 0.5f * x;
     46 	float t = .1e-8f;
     47 	float e = 1.0f;
     48 	float de = 1.0f;
     49 	for (int i = 1; i <= 25; ++i) {
     50 		de = de * y / (float)i;
     51 		float sde = de * de;
     52 		e += sde;
     53 		if (e * t - sde > 0.0f) {
     54 			break;
     55 		}
     56 	}
     57 	return e;
     58 }
     59 
     60 
     61 PlanckTaperWindow::PlanckTaperWindow(int size, int taperSamples) : Window(size) {
     62 	_window[0] = 0.0f;
     63 	_sum += _window[size - 1] = 1.0f;
     64 
     65 	for (int i = 1; i < taperSamples; ++i) {
     66 		float x = ((float)taperSamples / (float)i) - ((float)taperSamples / (float)(taperSamples - i));
     67 		x = 1.0f + exp(x);
     68 		x = 1.0f / x;
     69 		_sum += _window[i] = x;
     70 	}
     71 	int nOnes = size - 2 * taperSamples;
     72 	std::fill_n(_window + taperSamples, nOnes, 1.0f);
     73 	_sum += nOnes;
     74 	for (int i = 0; i < taperSamples; ++i) {
     75 		_sum += _window[size - 1 - i] = _window[i];
     76 	}
     77 }
     78 
     79 
     80 typedef ffft::FFTRealFixLen<10> FIXED_FFT1024;
     81 
     82 FFT1024::FFT1024() {
     83 	_fft = new FIXED_FFT1024();
     84 }
     85 
     86 FFT1024::~FFT1024() {
     87 	delete (FIXED_FFT1024*)_fft;
     88 }
     89 
     90 void FFT1024::do_fft(float* out, float* in) {
     91 	((FIXED_FFT1024*)_fft)->do_fft(out, in);
     92 }
     93 
     94 
     95 typedef ffft::FFTRealFixLen<12> FIXED_FFT4096;
     96 
     97 FFT4096::FFT4096() {
     98 	_fft = new FIXED_FFT4096();
     99 }
    100 
    101 FFT4096::~FFT4096() {
    102 	delete (FIXED_FFT4096*)_fft;
    103 }
    104 
    105 void FFT4096::do_fft(float* out, float* in) {
    106 	((FIXED_FFT4096*)_fft)->do_fft(out, in);
    107 }
    108 
    109 
    110 typedef ffft::FFTRealFixLen<13> FIXED_FFT8192;
    111 
    112 FFT8192::FFT8192() {
    113 	_fft = new FIXED_FFT8192();
    114 }
    115 
    116 FFT8192::~FFT8192() {
    117 	delete (FIXED_FFT8192*)_fft;
    118 }
    119 
    120 void FFT8192::do_fft(float* out, float* in) {
    121 	((FIXED_FFT8192*)_fft)->do_fft(out, in);
    122 }
    123 
    124 
    125 typedef ffft::FFTRealFixLen<14> FIXED_FFT16384;
    126 
    127 FFT16384::FFT16384() {
    128 	_fft = new FIXED_FFT16384();
    129 }
    130 
    131 FFT16384::~FFT16384() {
    132 	delete (FIXED_FFT16384*)_fft;
    133 }
    134 
    135 void FFT16384::do_fft(float* out, float* in) {
    136 	((FIXED_FFT16384*)_fft)->do_fft(out, in);
    137 }
    138 
    139 
    140 typedef ffft::FFTRealFixLen<15> FIXED_FFT32768;
    141 
    142 FFT32768::FFT32768() {
    143 	_fft = new FIXED_FFT32768();
    144 }
    145 
    146 FFT32768::~FFT32768() {
    147 	delete (FIXED_FFT32768*)_fft;
    148 }
    149 
    150 void FFT32768::do_fft(float* out, float* in) {
    151 	((FIXED_FFT32768*)_fft)->do_fft(out, in);
    152 }
    153 
    154 
    155 SpectrumAnalyzer::SpectrumAnalyzer(
    156 	Size size,
    157 	Overlap overlap,
    158 	WindowType windowType,
    159 	float sampleRate,
    160 	bool autoProcess
    161 )
    162 : OverlappingBuffer(size, overlap, autoProcess)
    163 , _sampleRate(sampleRate)
    164 {
    165 	assert(size <= maxSize);
    166 	assert(_sampleRate > size);
    167 
    168 	switch (size) {
    169 		case SIZE_1024: {
    170 			_fft1024 = new FFT1024();
    171 			break;
    172 		}
    173 		case SIZE_4096: {
    174 			_fft4096 = new FFT4096();
    175 			break;
    176 		}
    177 		case SIZE_8192: {
    178 			_fft8192 = new FFT8192();
    179 			break;
    180 		}
    181 		case SIZE_16384: {
    182 			_fft16384 = new FFT16384();
    183 			break;
    184 		}
    185 		case SIZE_32768: {
    186 			_fft32768 = new FFT32768();
    187 			break;
    188 		}
    189 		default: {
    190 			_fft = new ffft::FFTReal<float>(size);
    191 		}
    192 	}
    193 
    194 	switch (windowType) {
    195 		case WINDOW_NONE: {
    196 			break;
    197 		}
    198 		case WINDOW_HANNING: {
    199 			_window = new HanningWindow(size);
    200 			_windowOut = new float[size];
    201 			break;
    202 		}
    203 		case WINDOW_HAMMING: {
    204 			_window = new HammingWindow(size);
    205 			_windowOut = new float[size];
    206 			break;
    207 		}
    208 		case WINDOW_KAISER: {
    209 			_window = new KaiserWindow(size);
    210 			_windowOut = new float[size];
    211 			break;
    212 		}
    213 	}
    214 
    215 	_fftOut = new float[_size];
    216 }
    217 
    218 SpectrumAnalyzer::~SpectrumAnalyzer() {
    219 	if (_fft) {
    220 		delete _fft;
    221 	}
    222 	if (_fft1024) {
    223 		delete _fft1024;
    224 	}
    225 	if (_fft4096) {
    226 		delete _fft4096;
    227 	}
    228 	if (_fft8192) {
    229 		delete _fft8192;
    230 	}
    231 	if (_fft16384) {
    232 		delete _fft16384;
    233 	}
    234 	if (_fft32768) {
    235 		delete _fft32768;
    236 	}
    237 
    238 	if (_window) {
    239 		delete _window;
    240 		delete[] _windowOut;
    241 	}
    242 
    243 	delete[] _fftOut;
    244 }
    245 
    246 void SpectrumAnalyzer::processBuffer(float* samples) {
    247 	float* input = samples;
    248 	if (_window) {
    249 		_window->apply(samples, _windowOut);
    250 		input = _windowOut;
    251 	}
    252 	if (_fft1024) {
    253 		_fft1024->do_fft(_fftOut, input);
    254 	}
    255 	else if (_fft4096) {
    256 		_fft4096->do_fft(_fftOut, input);
    257 	}
    258 	else if (_fft8192) {
    259 		_fft8192->do_fft(_fftOut, input);
    260 	}
    261 	else if (_fft16384) {
    262 		_fft16384->do_fft(_fftOut, input);
    263 	}
    264 	else if (_fft32768) {
    265 		_fft32768->do_fft(_fftOut, input);
    266 	}
    267 	else {
    268 		_fft->do_fft(_fftOut, input);
    269 	}
    270 }
    271 
    272 void SpectrumAnalyzer::getMagnitudes(float* bins, int nBins) {
    273 	assert(nBins <= _size / 2);
    274 	assert(_size % nBins == 0);
    275 
    276 	const int bands = _size / 2;
    277 	const int binWidth = bands / nBins;
    278 	const float invBinWidth = 1.0 / (float)binWidth;
    279 	const float normalization = 2.0 / powf(_window ? _window->sum() : _size, 2.0);
    280 	for (int bin = 0; bin < nBins; ++bin) {
    281 		float sum = 0.0;
    282 		int binEnd = bin * binWidth + binWidth;
    283 		for (int i = binEnd - binWidth; i < binEnd; ++i) {
    284 			sum += (_fftOut[i]*_fftOut[i] + _fftOut[i + bands]*_fftOut[i + bands]) * normalization;
    285 		}
    286 		bins[bin] = sum * invBinWidth;
    287 	}
    288 }