Audacity  3.0.3
SpectrumAnalyst.cpp
Go to the documentation of this file.
1 /**********************************************************************
2 
3  Audacity: A Digital Audio Editor
4 
5  SpectrumAnalyst.cpp
6 
7  Dominic Mazzoni
8  Paul Licameli split from FreqWindow.cpp
9 
10 *******************************************************************//*******************************************************************/
19 
20 /*
21  Salvo Ventura - November 2006
22  Extended range check for additional FFT windows
23 */
24 
25 #include "SpectrumAnalyst.h"
26 #include "FFT.h"
27 
28 #include "SampleFormat.h"
29 #include <wx/dcclient.h>
30 
31 FreqGauge::FreqGauge(wxWindow * parent, wxWindowID winid)
32 : wxStatusBar(parent, winid, wxST_SIZEGRIP)
33 {
34  mRange = 0;
35 }
36 
37 void FreqGauge::SetRange(int range, int bar, int gap)
38 {
39  mRange = range;
40  mBar = bar;
41  mGap = gap;
42 
43  GetFieldRect(0, mRect);
44  mRect.Inflate(-1);
45 
46  mInterval = mRange / (mRect.width / (mBar + mGap));
47  mRect.width = mBar;
48  mMargin = mRect.x;
49  mLast = -1;
50 
51  Update();
52 }
53 
54 void FreqGauge::SetValue(int value)
55 {
56  mCur = value / mInterval;
57 
58  if (mCur != mLast)
59  {
60  wxClientDC dc(this);
61  dc.SetPen(*wxTRANSPARENT_PEN);
62  dc.SetBrush(wxColour(100, 100, 220));
63 
64  while (mLast < mCur)
65  {
66  mLast++;
67  mRect.x = mMargin + mLast * (mBar + mGap);
68  dc.DrawRectangle(mRect);
69  }
70  Update();
71  }
72 }
73 
75 {
76  mRange = 0;
77  Refresh(true);
78 }
79 
81 : mAlg(Spectrum)
82 , mRate(0.0)
83 , mWindowSize(0)
84 {
85 }
86 
88 {
89 }
90 
91 bool SpectrumAnalyst::Calculate(Algorithm alg, int windowFunc,
92  size_t windowSize, double rate,
93  const float *data, size_t dataLen,
94  float *pYMin, float *pYMax,
95  FreqGauge *progress)
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 }
337 
338 const float *SpectrumAnalyst::GetProcessed() const
339 {
340  return &mProcessed[0];
341 }
342 
344 {
345  return mProcessed.size() / 2;
346 }
347 
348 float SpectrumAnalyst::GetProcessedValue(float freq0, float freq1) const
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 }
396 
397 float SpectrumAnalyst::FindPeak(float xPos, float *pY) const
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 }
443 
444 // If f(0)=y0, f(1)=y1, f(2)=y2, and f(3)=y3, this function finds
445 // the degree-three polynomial which best fits these points and
446 // returns the value of this polynomial at a value x. Usually
447 // 0 < x < 3
448 
449 float SpectrumAnalyst::CubicInterpolate(float y0, float y1, float y2, float y3, float x) const
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 }
463 
464 float SpectrumAnalyst::CubicMaximize(float y0, float y1, float y2, float y3, float * max) const
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 }
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
FreqGauge::mRect
wxRect mRect
Definition: SpectrumAnalyst.h:72
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::mMargin
int mMargin
Definition: SpectrumAnalyst.h:79
FreqGauge::mInterval
int mInterval
Definition: SpectrumAnalyst.h:76
SpectrumAnalyst.h
gap
static const int gap
Definition: MeterPanel.cpp:257
FreqGauge::mRange
int mRange
Definition: SpectrumAnalyst.h:73
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
FFT.h
SpectrumAnalyst::GetProcessedSize
int GetProcessedSize() const
Definition: SpectrumAnalyst.cpp:343
FreqGauge::FreqGauge
FreqGauge(wxWindow *parent, wxWindowID winid)
Definition: SpectrumAnalyst.cpp:31
SpectrumAnalyst::mWindowSize
size_t mWindowSize
Definition: SpectrumAnalyst.h:58
FreqGauge::mGap
int mGap
Definition: SpectrumAnalyst.h:78
InverseRealFFT
void InverseRealFFT(size_t NumSamples, const float *RealIn, const float *ImagIn, float *RealOut)
Definition: FFT.cpp:268
SpectrumAnalyst::Calculate
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)
Definition: SpectrumAnalyst.cpp:91
FreqGauge::SetRange
void SetRange(int range, int bar=12, int gap=3)
Definition: SpectrumAnalyst.cpp:37
SpectrumAnalyst::~SpectrumAnalyst
~SpectrumAnalyst()
Definition: SpectrumAnalyst.cpp:87
FreqGauge
Definition: SpectrumAnalyst.h:63
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::FindPeak
float FindPeak(float xPos, float *pY) const
Definition: SpectrumAnalyst.cpp:397
SpectrumAnalyst::SpectrumAnalyst
SpectrumAnalyst()
Definition: SpectrumAnalyst.cpp:80
FreqGauge::mCur
int mCur
Definition: SpectrumAnalyst.h:74
SpectrumAnalyst::GetProcessedValue
float GetProcessedValue(float freq0, float freq1) const
Definition: SpectrumAnalyst.cpp:348
FreqGauge::mBar
int mBar
Definition: SpectrumAnalyst.h:77
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
WaveTrackViewConstants::Spectrum
@ Spectrum
Definition: WaveTrackViewConstants.h:29
SpectrumAnalyst::GetProcessed
const float * GetProcessed() const
Definition: SpectrumAnalyst.cpp:338
SpectrumAnalyst::NumAlgorithms
@ NumAlgorithms
Definition: SpectrumAnalyst.h:31
SpectrumAnalyst::Algorithm
Algorithm
Definition: SpectrumAnalyst.h:24
SpectrumAnalyst::mRate
double mRate
Definition: SpectrumAnalyst.h:57
SampleFormat.h
ArrayOf< float >
FreqGauge::mLast
int mLast
Definition: SpectrumAnalyst.h:75