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 }