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