Audacity  3.0.3
Spectrum.cpp
Go to the documentation of this file.
1 /**********************************************************************
2 
3  Audacity: A Digital Audio Editor
4 
5  Spectrum.cpp
6 
7  Dominic Mazzoni
8 
9 *******************************************************************//*******************************************************************/
15 
16 #include "Spectrum.h"
17 
18 #include <math.h>
19 
20 #include "SampleFormat.h"
21 
22 bool ComputeSpectrum(const float * data, size_t width,
23  size_t windowSize,
24  double WXUNUSED(rate), float *output,
25  bool autocorrelation, int windowFunc)
26 {
27  if (width < windowSize)
28  return false;
29 
30  if (!data || !output)
31  return true;
32 
33  Floats processed{ windowSize };
34 
35  for (size_t i = 0; i < windowSize; i++)
36  processed[i] = float(0.0);
37  auto half = windowSize / 2;
38 
39  Floats in{ windowSize };
40  Floats out{ windowSize };
41  Floats out2{ windowSize };
42 
43  size_t start = 0;
44  unsigned windows = 0;
45  while (start + windowSize <= width) {
46  for (size_t i = 0; i < windowSize; i++)
47  in[i] = data[start + i];
48 
49  WindowFunc(windowFunc, windowSize, in.get());
50 
51  if (autocorrelation) {
52  // Take FFT
53  RealFFT(windowSize, in.get(), out.get(), out2.get());
54  // Compute power
55  for (size_t i = 0; i < windowSize; i++)
56  in[i] = (out[i] * out[i]) + (out2[i] * out2[i]);
57 
58  // Tolonen and Karjalainen recommend taking the cube root
59  // of the power, instead of the square root
60 
61  for (size_t i = 0; i < windowSize; i++)
62  in[i] = powf(in[i], 1.0f / 3.0f);
63 
64  // Take FFT
65  RealFFT(windowSize, in.get(), out.get(), out2.get());
66  }
67  else
68  PowerSpectrum(windowSize, in.get(), out.get());
69 
70  // Take real part of result
71  for (size_t i = 0; i < half; i++)
72  processed[i] += out[i];
73 
74  start += half;
75  windows++;
76  }
77 
78  if (autocorrelation) {
79 
80  // Peak Pruning as described by Tolonen and Karjalainen, 2000
81  /*
82  Combine most of the calculations in a single for loop.
83  It should be safe, as indexes refer only to current and previous elements,
84  that have already been clipped, etc...
85  */
86  for (size_t i = 0; i < half; i++) {
87  // Clip at zero, copy to temp array
88  if (processed[i] < 0.0)
89  processed[i] = float(0.0);
90  out[i] = processed[i];
91  // Subtract a time-doubled signal (linearly interp.) from the original
92  // (clipped) signal
93  if ((i % 2) == 0)
94  processed[i] -= out[i / 2];
95  else
96  processed[i] -= ((out[i / 2] + out[i / 2 + 1]) / 2);
97 
98  // Clip at zero again
99  if (processed[i] < 0.0)
100  processed[i] = float(0.0);
101  }
102 
103  // Reverse and scale
104  for (size_t i = 0; i < half; i++)
105  in[i] = processed[i] / (windowSize / 4);
106  for (size_t i = 0; i < half; i++)
107  processed[half - 1 - i] = in[i];
108  } else {
109  // Convert to decibels
110  // But do it safely; -Inf is nobody's friend
111  for (size_t i = 0; i < half; i++){
112  float temp=(processed[i] / windowSize / windows);
113  if (temp > 0.0)
114  processed[i] = 10 * log10(temp);
115  else
116  processed[i] = 0;
117  }
118  }
119 
120  for(size_t i = 0; i < half; i++)
121  output[i] = processed[i];
122 
123  return true;
124 }
125 
WindowFunc
void WindowFunc(int whichFunction, size_t NumSamples, float *in)
Definition: FFT.cpp:515
PowerSpectrum
void PowerSpectrum(size_t NumSamples, const float *In, float *Out)
Definition: FFT.cpp:304
ComputeSpectrum
bool ComputeSpectrum(const float *data, size_t width, size_t windowSize, double WXUNUSED(rate), float *output, bool autocorrelation, int windowFunc)
Definition: Spectrum.cpp:22
RealFFT
void RealFFT(size_t NumSamples, const float *RealIn, float *RealOut, float *ImagOut)
Definition: FFT.cpp:230
Spectrum.h
SampleFormat.h
ArrayOf< float >