1#include "./common.h"
2
3#ifndef SIGNALSMITH_DSP_WINDOWS_H
4#define SIGNALSMITH_DSP_WINDOWS_H
5
6#include <cmath>
7#include <algorithm>
8
9namespace signalsmith {
10namespace windows {
11 /** @defgroup Windows Window functions
12 @brief Windows for spectral analysis
13
14 These are generally double-precision, because they are mostly calculated during setup/reconfiguring, not real-time code.
15
16 @{
17 @file
18 */
19
20 /** @brief The Kaiser window (almost) maximises the energy in the main-lobe compared to the side-lobes.
21
22 Kaiser windows can be constructing using the shape-parameter (beta) or using the static `with???()` methods.*/
23 class Kaiser {
24 // I_0(x)=\sum_{k=0}^{N}\frac{x^{2k}}{(k!)^2\cdot4^k}
25 inline static double bessel0(double x) {
26 const double significanceLimit = 1e-4;
27 double result = 0;
28 double term = 1;
29 double m = 0;
30 while (term > significanceLimit) {
31 result += term;
32 ++m;
33 term *= (x*x)/(4*m*m);
34 }
35
36 return result;
37 }
38 double beta;
39 double invB0;
40
41 static double heuristicBandwidth(double bandwidth) {
42 // Good peaks
43 //return bandwidth + 8/((bandwidth + 3)*(bandwidth + 3));
44 // Good average
45 //return bandwidth + 14/((bandwidth + 2.5)*(bandwidth + 2.5));
46 // Compromise
47 return bandwidth + 8/((bandwidth + 3)*(bandwidth + 3)) + 0.25*std::max(a: 3 - bandwidth, b: 0.0);
48 }
49 public:
50 /// Set up a Kaiser window with a given shape. `beta` is `pi*alpha` (since there is ambiguity about shape parameters)
51 Kaiser(double beta) : beta(beta), invB0(1/bessel0(x: beta)) {}
52
53 /// @name Bandwidth methods
54 /// @{
55 static Kaiser withBandwidth(double bandwidth, bool heuristicOptimal=false) {
56 return Kaiser(bandwidthToBeta(bandwidth, heuristicOptimal));
57 }
58
59 /** Returns the Kaiser shape where the main lobe has the specified bandwidth (as a factor of 1/window-length).
60 \diagram{kaiser-windows.svg,You can see that the main lobe matches the specified bandwidth.}
61 If `heuristicOptimal` is enabled, the main lobe width is _slightly_ wider, improving both the peak and total energy - see `bandwidthToEnergyDb()` and `bandwidthToPeakDb()`.
62 \diagram{kaiser-windows-heuristic.svg, The main lobe extends to ±bandwidth/2.} */
63 static double bandwidthToBeta(double bandwidth, bool heuristicOptimal=false) {
64 if (heuristicOptimal) { // Heuristic based on numerical search
65 bandwidth = heuristicBandwidth(bandwidth);
66 }
67 bandwidth = std::max(a: bandwidth, b: 2.0);
68 double alpha = std::sqrt(x: bandwidth*bandwidth*0.25 - 1);
69 return alpha*M_PI;
70 }
71
72 static double betaToBandwidth(double beta) {
73 double alpha = beta*(1.0/M_PI);
74 return 2*std::sqrt(x: alpha*alpha + 1);
75 }
76 /// @}
77
78 /// @name Performance methods
79 /// @{
80 /** @brief Total energy ratio (in dB) between side-lobes and the main lobe.
81 \diagram{windows-kaiser-sidelobe-energy.svg,Measured main/side lobe energy ratio. You can see that the heuristic improves performance for all bandwidth values.}
82 This function uses an approximation which is accurate to ±0.5dB for 2 ⩽ bandwidth ≤ 10, or 1 ⩽ bandwidth ≤ 10 when `heuristicOptimal`is enabled.
83 */
84 static double bandwidthToEnergyDb(double bandwidth, bool heuristicOptimal=false) {
85 // Horrible heuristic fits
86 if (heuristicOptimal) {
87 if (bandwidth < 3) bandwidth += (3 - bandwidth)*0.5;
88 return 12.9 + -3/(bandwidth + 0.4) - 13.4*bandwidth + (bandwidth < 3)*-9.6*(bandwidth - 3);
89 }
90 return 10.5 + 15/(bandwidth + 0.4) - 13.25*bandwidth + (bandwidth < 2)*13*(bandwidth - 2);
91 }
92 static double energyDbToBandwidth(double energyDb, bool heuristicOptimal=false) {
93 double bw = 1;
94 while (bw < 20 && bandwidthToEnergyDb(bandwidth: bw, heuristicOptimal) > energyDb) {
95 bw *= 2;
96 }
97 double step = bw/2;
98 while (step > 0.0001) {
99 if (bandwidthToEnergyDb(bandwidth: bw, heuristicOptimal) > energyDb) {
100 bw += step;
101 } else {
102 bw -= step;
103 }
104 step *= 0.5;
105 }
106 return bw;
107 }
108 /** @brief Peak ratio (in dB) between side-lobes and the main lobe.
109 \diagram{windows-kaiser-sidelobe-peaks.svg,Measured main/side lobe peak ratio. You can see that the heuristic improves performance, except in the bandwidth range 1-2 where peak ratio was sacrificed to improve total energy ratio.}
110 This function uses an approximation which is accurate to ±0.5dB for 2 ⩽ bandwidth ≤ 9, or 0.5 ⩽ bandwidth ≤ 9 when `heuristicOptimal`is enabled.
111 */
112 static double bandwidthToPeakDb(double bandwidth, bool heuristicOptimal=false) {
113 // Horrible heuristic fits
114 if (heuristicOptimal) {
115 return 14.2 - 20/(bandwidth + 1) - 13*bandwidth + (bandwidth < 3)*-6*(bandwidth - 3) + (bandwidth < 2.25)*5.8*(bandwidth - 2.25);
116 }
117 return 10 + 8/(bandwidth + 2) - 12.75*bandwidth + (bandwidth < 2)*4*(bandwidth - 2);
118 }
119 static double peakDbToBandwidth(double peakDb, bool heuristicOptimal=false) {
120 double bw = 1;
121 while (bw < 20 && bandwidthToPeakDb(bandwidth: bw, heuristicOptimal) > peakDb) {
122 bw *= 2;
123 }
124 double step = bw/2;
125 while (step > 0.0001) {
126 if (bandwidthToPeakDb(bandwidth: bw, heuristicOptimal) > peakDb) {
127 bw += step;
128 } else {
129 bw -= step;
130 }
131 step *= 0.5;
132 }
133 return bw;
134 }
135 /** @} */
136
137 /** Equivalent noise bandwidth (ENBW), a measure of frequency resolution.
138 \diagram{windows-kaiser-enbw.svg,Measured ENBW\, with and without the heuristic bandwidth adjustment.}
139 This approximation is accurate to ±0.05 up to a bandwidth of 22.
140 */
141 static double bandwidthToEnbw(double bandwidth, bool heuristicOptimal=false) {
142 if (heuristicOptimal) bandwidth = heuristicBandwidth(bandwidth);
143 double b2 = std::max<double>(a: bandwidth - 2, b: 0);
144 return 1 + b2*(0.2 + b2*(-0.005 + b2*(-0.000005 + b2*0.0000022)));
145 }
146
147 /// Return the window's value for position in the range [0, 1]
148 double operator ()(double unit) {
149 double r = 2*unit - 1;
150 double arg = std::sqrt(x: 1 - r*r);
151 return bessel0(x: beta*arg)*invB0;
152 }
153
154 /// Fills an arbitrary container with a Kaiser window
155 template<typename Data>
156 void fill(Data &&data, int size) const {
157 double invSize = 1.0/size;
158 for (int i = 0; i < size; ++i) {
159 double r = (2*i + 1)*invSize - 1;
160 double arg = std::sqrt(x: 1 - r*r);
161 data[i] = bessel0(x: beta*arg)*invB0;
162 }
163 }
164 };
165
166 /** @brief The Approximate Confined Gaussian window is (almost) optimal
167
168 ACG windows can be constructing using the shape-parameter (sigma) or using the static `with???()` methods.*/
169 class ApproximateConfinedGaussian {
170 double gaussianFactor;
171
172 double gaussian(double x) const {
173 return std::exp(x: -x*x*gaussianFactor);
174 }
175 public:
176 /// Heuristic map from bandwidth to the appropriately-optimal sigma
177 static double bandwidthToSigma(double bandwidth) {
178 return 0.3/std::sqrt(x: bandwidth);
179 }
180 static ApproximateConfinedGaussian withBandwidth(double bandwidth) {
181 return ApproximateConfinedGaussian(bandwidthToSigma(bandwidth));
182 }
183
184 ApproximateConfinedGaussian(double sigma) : gaussianFactor(0.0625/(sigma*sigma)) {}
185
186 /// Fills an arbitrary container
187 template<typename Data>
188 void fill(Data &&data, int size) const {
189 double invSize = 1.0/size;
190 double offsetScale = gaussian(x: 1)/(gaussian(x: 3) + gaussian(x: -1));
191 double norm = 1/(gaussian(x: 0) - 2*offsetScale*(gaussian(x: 2)));
192 for (int i = 0; i < size; ++i) {
193 double r = (2*i + 1)*invSize - 1;
194 data[i] = norm*(gaussian(x: r) - offsetScale*(gaussian(x: r - 2) + gaussian(x: r + 2)));
195 }
196 }
197 };
198
199 /** Forces STFT perfect-reconstruction (WOLA) on an existing window, for a given STFT interval.
200 For example, here are perfect-reconstruction versions of the approximately-optimal @ref Kaiser windows:
201 \diagram{kaiser-windows-heuristic-pr.svg,Note the lower overall energy\, and the pointy top for 2x bandwidth. Spectral performance is about the same\, though.}
202 */
203 template<typename Data>
204 void forcePerfectReconstruction(Data &&data, int windowLength, int interval) {
205 for (int i = 0; i < interval; ++i) {
206 double sum2 = 0;
207 for (int index = i; index < windowLength; index += interval) {
208 sum2 += data[index]*data[index];
209 }
210 double factor = 1/std::sqrt(x: sum2);
211 for (int index = i; index < windowLength; index += interval) {
212 data[index] *= factor;
213 }
214 }
215 }
216
217/** @} */
218}} // signalsmith::windows
219#endif // include guard
220

source code of qtmultimedia/src/3rdparty/signalsmith-stretch/dsp/windows.h