Audacity  3.0.3
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.

24  {
25  Spectrum,
29  Cepstrum,
30 
32  };

Constructor & Destructor Documentation

◆ SpectrumAnalyst()

SpectrumAnalyst::SpectrumAnalyst ( )

Definition at line 80 of file SpectrumAnalyst.cpp.

81 : mAlg(Spectrum)
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 &&
106  alg >= SpectrumAnalyst::Spectrum &&
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 ||
179  alg == EnhancedAutocorrelation) {
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 }

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, 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 }

Referenced by FindPeak().

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 }

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 }

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:
WindowFunc
void WindowFunc(int whichFunction, size_t NumSamples, float *in)
Definition: FFT.cpp:515
SpectrumAnalyst::Spectrum
@ Spectrum
Definition: SpectrumAnalyst.h:25
PowerSpectrum
void PowerSpectrum(size_t NumSamples, const float *In, float *Out)
Definition: FFT.cpp:304
fast_float::detail::power
constexpr fastfloat_really_inline int32_t power(int32_t q) noexcept
Definition: fast_float.h:1454
SpectrumAnalyst::CubicInterpolate
float CubicInterpolate(float y0, float y1, float y2, float y3, float x) const
Definition: SpectrumAnalyst.cpp:449
NumWindowFuncs
int NumWindowFuncs()
Definition: FFT.cpp:329
FreqGauge::SetValue
void SetValue(int value)
Definition: SpectrumAnalyst.cpp:54
RealFFT
void RealFFT(size_t NumSamples, const float *RealIn, float *RealOut, float *ImagOut)
Definition: FFT.cpp:230
FreqGauge::Reset
void Reset()
Definition: SpectrumAnalyst.cpp:74
SpectrumAnalyst::EnhancedAutocorrelation
@ EnhancedAutocorrelation
Definition: SpectrumAnalyst.h:28
SpectrumAnalyst::GetProcessedSize
int GetProcessedSize() const
Definition: SpectrumAnalyst.cpp:343
SpectrumAnalyst::mWindowSize
size_t mWindowSize
Definition: SpectrumAnalyst.h:58
InverseRealFFT
void InverseRealFFT(size_t NumSamples, const float *RealIn, const float *ImagIn, float *RealOut)
Definition: FFT.cpp:268
FreqGauge::SetRange
void SetRange(int range, int bar=12, int gap=3)
Definition: SpectrumAnalyst.cpp:37
SpectrumAnalyst::Cepstrum
@ Cepstrum
Definition: SpectrumAnalyst.h:29
SpectrumAnalyst::CubicMaximize
float CubicMaximize(float y0, float y1, float y2, float y3, float *max) const
Definition: SpectrumAnalyst.cpp:464
SpectrumAnalyst::mProcessed
std::vector< float > mProcessed
Definition: SpectrumAnalyst.h:59
SpectrumAnalyst::Autocorrelation
@ Autocorrelation
Definition: SpectrumAnalyst.h:26
SpectrumAnalyst::CubeRootAutocorrelation
@ CubeRootAutocorrelation
Definition: SpectrumAnalyst.h:27
SpectrumAnalyst::mAlg
Algorithm mAlg
Definition: SpectrumAnalyst.h:56
SpectrumAnalyst::NumAlgorithms
@ NumAlgorithms
Definition: SpectrumAnalyst.h:31
SpectrumAnalyst::mRate
double mRate
Definition: SpectrumAnalyst.h:57
ArrayOf< float >