Audacity 3.2.0
Public Types | Public Member Functions | Private Member Functions | Private Attributes | List of all members
SpectrumAnalyst Class Reference

Used for finding the peaks, for snapping to peaks. More...

#include <SpectrumAnalyst.h>

Collaboration diagram for SpectrumAnalyst:
[legend]

Public Types

enum  Algorithm {
  Spectrum , Autocorrelation , CubeRootAutocorrelation , EnhancedAutocorrelation ,
  Cepstrum , NumAlgorithms
}
 

Public Member Functions

 SpectrumAnalyst ()
 
 ~SpectrumAnalyst ()
 
bool Calculate (Algorithm alg, int windowFunc, size_t windowSize, double rate, const float *data, size_t dataLen, float *pYMin=NULL, float *pYMax=NULL, FreqGauge *progress=NULL)
 
const float * GetProcessed () const
 
int GetProcessedSize () const
 
float GetProcessedValue (float freq0, float freq1) const
 
float FindPeak (float xPos, float *pY) const
 

Private Member Functions

float CubicInterpolate (float y0, float y1, float y2, float y3, float x) const
 
float CubicMaximize (float y0, float y1, float y2, float y3, float *max) const
 

Private Attributes

Algorithm mAlg
 
double mRate
 
size_t mWindowSize
 
std::vector< float > mProcessed
 

Detailed Description

Used for finding the peaks, for snapping to peaks.

This class is used to do the 'find peaks' snapping both in FreqPlot and in the spectrogram spectral selection.

Definition at line 20 of file SpectrumAnalyst.h.

Member Enumeration Documentation

◆ Algorithm

Enumerator
Spectrum 
Autocorrelation 
CubeRootAutocorrelation 
EnhancedAutocorrelation 
Cepstrum 
NumAlgorithms 

Definition at line 24 of file SpectrumAnalyst.h.

Constructor & Destructor Documentation

◆ SpectrumAnalyst()

SpectrumAnalyst::SpectrumAnalyst ( )

Definition at line 80 of file SpectrumAnalyst.cpp.

82, mRate(0.0)
83, mWindowSize(0)
84{
85}

◆ ~SpectrumAnalyst()

SpectrumAnalyst::~SpectrumAnalyst ( )

Definition at line 87 of file SpectrumAnalyst.cpp.

88{
89}

Member Function Documentation

◆ Calculate()

bool SpectrumAnalyst::Calculate ( Algorithm  alg,
int  windowFunc,
size_t  windowSize,
double  rate,
const float *  data,
size_t  dataLen,
float *  pYMin = NULL,
float *  pYMax = NULL,
FreqGauge progress = NULL 
)

Definition at line 91 of file SpectrumAnalyst.cpp.

96{
97 // Wipe old data
98 mProcessed.resize(0);
99 mRate = 0.0;
100 mWindowSize = 0;
101
102 // Validate inputs
103 int f = NumWindowFuncs();
104
105 if (!(windowSize >= 32 && windowSize <= 131072 &&
108 windowFunc >= 0 && windowFunc < f)) {
109 return false;
110 }
111
112 if (dataLen < windowSize) {
113 return false;
114 }
115
116 // Now repopulate
117 mRate = rate;
118 mWindowSize = windowSize;
119 mAlg = alg;
120
121 auto half = mWindowSize / 2;
122 mProcessed.resize(mWindowSize);
123
124 Floats in{ mWindowSize };
125 Floats out{ mWindowSize };
126 Floats out2{ mWindowSize };
127 Floats win{ mWindowSize };
128
129 for (size_t i = 0; i < mWindowSize; i++) {
130 mProcessed[i] = 0.0f;
131 win[i] = 1.0f;
132 }
133
134 WindowFunc(windowFunc, mWindowSize, win.get());
135
136 // Scale window such that an amplitude of 1.0 in the time domain
137 // shows an amplitude of 0dB in the frequency domain
138 double wss = 0;
139 for (size_t i = 0; i<mWindowSize; i++)
140 wss += win[i];
141 if(wss > 0)
142 wss = 4.0 / (wss*wss);
143 else
144 wss = 1.0;
145
146 if (progress) {
147 progress->SetRange(dataLen);
148 }
149
150 size_t start = 0;
151 int windows = 0;
152 while (start + mWindowSize <= dataLen) {
153 for (size_t i = 0; i < mWindowSize; i++)
154 in[i] = win[i] * data[start + i];
155
156 switch (alg) {
157 case Spectrum:
158 PowerSpectrum(mWindowSize, in.get(), out.get());
159
160 for (size_t i = 0; i < half; i++)
161 mProcessed[i] += out[i];
162 break;
163
164 case Autocorrelation:
167
168 // Take FFT
169 RealFFT(mWindowSize, in.get(), out.get(), out2.get());
170 // Compute power
171 for (size_t i = 0; i < mWindowSize; i++)
172 in[i] = (out[i] * out[i]) + (out2[i] * out2[i]);
173
174 if (alg == Autocorrelation) {
175 for (size_t i = 0; i < mWindowSize; i++)
176 in[i] = sqrt(in[i]);
177 }
178 if (alg == CubeRootAutocorrelation ||
180 // Tolonen and Karjalainen recommend taking the cube root
181 // of the power, instead of the square root
182
183 for (size_t i = 0; i < mWindowSize; i++)
184 in[i] = pow(in[i], 1.0f / 3.0f);
185 }
186 // Take FFT
187 RealFFT(mWindowSize, in.get(), out.get(), out2.get());
188
189 // Take real part of result
190 for (size_t i = 0; i < half; i++)
191 mProcessed[i] += out[i];
192 break;
193
194 case Cepstrum:
195 RealFFT(mWindowSize, in.get(), out.get(), out2.get());
196
197 // Compute log power
198 // Set a sane lower limit assuming maximum time amplitude of 1.0
199 {
200 float power;
201 float minpower = 1e-20*mWindowSize*mWindowSize;
202 for (size_t i = 0; i < mWindowSize; i++)
203 {
204 power = (out[i] * out[i]) + (out2[i] * out2[i]);
205 if(power < minpower)
206 in[i] = log(minpower);
207 else
208 in[i] = log(power);
209 }
210 // Take IFFT
211 InverseRealFFT(mWindowSize, in.get(), NULL, out.get());
212
213 // Take real part of result
214 for (size_t i = 0; i < half; i++)
215 mProcessed[i] += out[i];
216 }
217
218 break;
219
220 default:
221 wxASSERT(false);
222 break;
223 } //switch
224
225 // Update the progress bar
226 if (progress) {
227 progress->SetValue(start);
228 }
229
230 start += half;
231 windows++;
232 }
233
234 if (progress) {
235 // Reset for next time
236 progress->Reset();
237 }
238
239 float mYMin = 1000000, mYMax = -1000000;
240 double scale;
241 switch (alg) {
242 case Spectrum:
243 // Convert to decibels
244 mYMin = 1000000.;
245 mYMax = -1000000.;
246 scale = wss / (double)windows;
247 for (size_t i = 0; i < half; i++)
248 {
249 mProcessed[i] = 10 * log10(mProcessed[i] * scale);
250 if(mProcessed[i] > mYMax)
251 mYMax = mProcessed[i];
252 else if(mProcessed[i] < mYMin)
253 mYMin = mProcessed[i];
254 }
255 break;
256
257 case Autocorrelation:
259 for (size_t i = 0; i < half; i++)
260 mProcessed[i] = mProcessed[i] / windows;
261
262 // Find min/max
263 mYMin = mProcessed[0];
264 mYMax = mProcessed[0];
265 for (size_t i = 1; i < half; i++)
266 if (mProcessed[i] > mYMax)
267 mYMax = mProcessed[i];
268 else if (mProcessed[i] < mYMin)
269 mYMin = mProcessed[i];
270 break;
271
273 for (size_t i = 0; i < half; i++)
274 mProcessed[i] = mProcessed[i] / windows;
275
276 // Peak Pruning as described by Tolonen and Karjalainen, 2000
277
278 // Clip at zero, copy to temp array
279 for (size_t i = 0; i < half; i++) {
280 if (mProcessed[i] < 0.0)
281 mProcessed[i] = float(0.0);
282 out[i] = mProcessed[i];
283 }
284
285 // Subtract a time-doubled signal (linearly interp.) from the original
286 // (clipped) signal
287 for (size_t i = 0; i < half; i++)
288 if ((i % 2) == 0)
289 mProcessed[i] -= out[i / 2];
290 else
291 mProcessed[i] -= ((out[i / 2] + out[i / 2 + 1]) / 2);
292
293 // Clip at zero again
294 for (size_t i = 0; i < half; i++)
295 if (mProcessed[i] < 0.0)
296 mProcessed[i] = float(0.0);
297
298 // Find NEW min/max
299 mYMin = mProcessed[0];
300 mYMax = mProcessed[0];
301 for (size_t i = 1; i < half; i++)
302 if (mProcessed[i] > mYMax)
303 mYMax = mProcessed[i];
304 else if (mProcessed[i] < mYMin)
305 mYMin = mProcessed[i];
306 break;
307
308 case Cepstrum:
309 for (size_t i = 0; i < half; i++)
310 mProcessed[i] = mProcessed[i] / windows;
311
312 // Find min/max, ignoring first and last few values
313 {
314 size_t ignore = 4;
315 mYMin = mProcessed[ignore];
316 mYMax = mProcessed[ignore];
317 for (size_t i = ignore + 1; i + ignore < half; i++)
318 if (mProcessed[i] > mYMax)
319 mYMax = mProcessed[i];
320 else if (mProcessed[i] < mYMin)
321 mYMin = mProcessed[i];
322 }
323 break;
324
325 default:
326 wxASSERT(false);
327 break;
328 }
329
330 if (pYMin)
331 *pYMin = mYMin;
332 if (pYMax)
333 *pYMax = mYMax;
334
335 return true;
336}
void WindowFunc(int whichFunction, size_t NumSamples, float *in)
Definition: FFT.cpp:513
int NumWindowFuncs()
Definition: FFT.cpp:327
void PowerSpectrum(size_t NumSamples, const float *In, float *Out)
Definition: FFT.cpp:302
void InverseRealFFT(size_t NumSamples, const float *RealIn, const float *ImagIn, float *RealOut)
Definition: FFT.cpp:266
void RealFFT(size_t NumSamples, const float *RealIn, float *RealOut, float *ImagOut)
Definition: FFT.cpp:228
void SetRange(int range, int bar=12, int gap=3)
void SetValue(int value)
std::vector< float > mProcessed
constexpr fastfloat_really_inline int32_t power(int32_t q) noexcept
Definition: fast_float.h:1454
__finl float_x4 __vecc sqrt(const float_x4 &a)

References Autocorrelation, Cepstrum, CubeRootAutocorrelation, EnhancedAutocorrelation, InverseRealFFT(), mAlg, mProcessed, mRate, mWindowSize, NumAlgorithms, NumWindowFuncs(), fast_float::detail::power(), PowerSpectrum(), RealFFT(), FreqGauge::Reset(), FreqGauge::SetRange(), FreqGauge::SetValue(), Spectrum, staffpad::audio::simd::sqrt(), and WindowFunc().

Referenced by SelectHandle::StartSnappingFreqSelection().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ CubicInterpolate()

float SpectrumAnalyst::CubicInterpolate ( float  y0,
float  y1,
float  y2,
float  y3,
float  x 
) const
private

Definition at line 449 of file SpectrumAnalyst.cpp.

450{
451 float a, b, c, d;
452
453 a = y0 / -6.0 + y1 / 2.0 - y2 / 2.0 + y3 / 6.0;
454 b = y0 - 5.0 * y1 / 2.0 + 2.0 * y2 - y3 / 2.0;
455 c = -11.0 * y0 / 6.0 + 3.0 * y1 - 3.0 * y2 / 2.0 + y3 / 3.0;
456 d = y0;
457
458 float xx = x * x;
459 float xxx = xx * x;
460
461 return (a * xxx + b * xx + c * x + d);
462}

Referenced by GetProcessedValue().

Here is the caller graph for this function:

◆ CubicMaximize()

float SpectrumAnalyst::CubicMaximize ( float  y0,
float  y1,
float  y2,
float  y3,
float *  max 
) const
private

Definition at line 464 of file SpectrumAnalyst.cpp.

465{
466 // Find coefficients of cubic
467
468 float a, b, c, d;
469
470 a = y0 / -6.0 + y1 / 2.0 - y2 / 2.0 + y3 / 6.0;
471 b = y0 - 5.0 * y1 / 2.0 + 2.0 * y2 - y3 / 2.0;
472 c = -11.0 * y0 / 6.0 + 3.0 * y1 - 3.0 * y2 / 2.0 + y3 / 3.0;
473 d = y0;
474
475 // Take derivative
476
477 float da, db, dc;
478
479 da = 3 * a;
480 db = 2 * b;
481 dc = c;
482
483 // Find zeroes of derivative using quadratic equation
484
485 float discriminant = db * db - 4 * da * dc;
486 if (discriminant < 0.0)
487 return float(-1.0); // error
488
489 float x1 = (-db + sqrt(discriminant)) / (2 * da);
490 float x2 = (-db - sqrt(discriminant)) / (2 * da);
491
492 // The one which corresponds to a local _maximum_ in the
493 // cubic is the one we want - the one with a negative
494 // second derivative
495
496 float dda = 2 * da;
497 float ddb = db;
498
499 if (dda * x1 + ddb < 0)
500 {
501 *max = a*x1*x1*x1+b*x1*x1+c*x1+d;
502 return x1;
503 }
504 else
505 {
506 *max = a*x2*x2*x2+b*x2*x2+c*x2+d;
507 return x2;
508 }
509}

References staffpad::audio::simd::sqrt().

Referenced by FindPeak().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ FindPeak()

float SpectrumAnalyst::FindPeak ( float  xPos,
float *  pY 
) const

Definition at line 397 of file SpectrumAnalyst.cpp.

398{
399 float bestpeak = 0.0f;
400 float bestValue = 0.0;
401 if (GetProcessedSize() > 1) {
402 bool up = (mProcessed[1] > mProcessed[0]);
403 float bestdist = 1000000;
404 for (int bin = 3; bin < GetProcessedSize() - 1; bin++) {
405 bool nowUp = mProcessed[bin] > mProcessed[bin - 1];
406 if (!nowUp && up) {
407 // Local maximum. Find actual value by cubic interpolation
408 int leftbin = bin - 2;
409 /*
410 if (leftbin < 1)
411 leftbin = 1;
412 */
413 float valueAtMax = 0.0;
414 float max = leftbin + CubicMaximize(mProcessed[leftbin],
415 mProcessed[leftbin + 1],
416 mProcessed[leftbin + 2],
417 mProcessed[leftbin + 3],
418 &valueAtMax);
419
420 float thispeak;
421 if (mAlg == Spectrum)
422 thispeak = max * mRate / mWindowSize;
423 else
424 thispeak = max / mRate;
425
426 if (fabs(thispeak - xPos) < bestdist) {
427 bestpeak = thispeak;
428 bestdist = fabs(thispeak - xPos);
429 bestValue = valueAtMax;
430 // Should this test come after the enclosing if?
431 if (thispeak > xPos)
432 break;
433 }
434 }
435 up = nowUp;
436 }
437 }
438
439 if (pY)
440 *pY = bestValue;
441 return bestpeak;
442}
float CubicMaximize(float y0, float y1, float y2, float y3, float *max) const
int GetProcessedSize() const

References CubicMaximize(), GetProcessedSize(), mAlg, mProcessed, mRate, mWindowSize, and Spectrum.

Referenced by SelectHandle::SnapCenterOnce().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetProcessed()

const float * SpectrumAnalyst::GetProcessed ( ) const

Definition at line 338 of file SpectrumAnalyst.cpp.

339{
340 return &mProcessed[0];
341}

References mProcessed.

◆ GetProcessedSize()

int SpectrumAnalyst::GetProcessedSize ( ) const

Definition at line 343 of file SpectrumAnalyst.cpp.

344{
345 return mProcessed.size() / 2;
346}

References mProcessed.

Referenced by FindPeak(), and GetProcessedValue().

Here is the caller graph for this function:

◆ GetProcessedValue()

float SpectrumAnalyst::GetProcessedValue ( float  freq0,
float  freq1 
) const

Definition at line 348 of file SpectrumAnalyst.cpp.

349{
350 float bin0, bin1, binwidth;
351
352 if (mAlg == Spectrum) {
353 bin0 = freq0 * mWindowSize / mRate;
354 bin1 = freq1 * mWindowSize / mRate;
355 } else {
356 bin0 = freq0 * mRate;
357 bin1 = freq1 * mRate;
358 }
359 binwidth = bin1 - bin0;
360
361 float value = float(0.0);
362
363 if (binwidth < 1.0) {
364 float binmid = (bin0 + bin1) / 2.0;
365 int ibin = (int)(binmid) - 1;
366 if (ibin < 1)
367 ibin = 1;
368 if (ibin >= GetProcessedSize() - 3)
369 ibin = std::max(0, GetProcessedSize() - 4);
370
371 value = CubicInterpolate(mProcessed[ibin],
372 mProcessed[ibin + 1],
373 mProcessed[ibin + 2],
374 mProcessed[ibin + 3], binmid - ibin);
375
376 } else {
377 if (bin0 < 0)
378 bin0 = 0;
379 if (bin1 >= GetProcessedSize())
380 bin1 = GetProcessedSize() - 1;
381
382 if ((int)(bin1) > (int)(bin0))
383 value += mProcessed[(int)(bin0)] * ((int)(bin0) + 1 - bin0);
384 bin0 = 1 + (int)(bin0);
385 while (bin0 < (int)(bin1)) {
386 value += mProcessed[(int)(bin0)];
387 bin0 += 1.0;
388 }
389 value += mProcessed[(int)(bin1)] * (bin1 - (int)(bin1));
390
391 value /= binwidth;
392 }
393
394 return value;
395}
float CubicInterpolate(float y0, float y1, float y2, float y3, float x) const

References CubicInterpolate(), GetProcessedSize(), mAlg, mProcessed, mRate, mWindowSize, and Spectrum.

Here is the call graph for this function:

Member Data Documentation

◆ mAlg

Algorithm SpectrumAnalyst::mAlg
private

Definition at line 56 of file SpectrumAnalyst.h.

Referenced by Calculate(), FindPeak(), and GetProcessedValue().

◆ mProcessed

std::vector<float> SpectrumAnalyst::mProcessed
private

◆ mRate

double SpectrumAnalyst::mRate
private

Definition at line 57 of file SpectrumAnalyst.h.

Referenced by Calculate(), FindPeak(), and GetProcessedValue().

◆ mWindowSize

size_t SpectrumAnalyst::mWindowSize
private

Definition at line 58 of file SpectrumAnalyst.h.

Referenced by Calculate(), FindPeak(), and GetProcessedValue().


The documentation for this class was generated from the following files: