Add Kaiser performance-estimation methods

spectral.h

@@ -124,7 +124,7 @@ namespace spectral {
 
 	/** STFT synthesis/analysis/processing, built on a `MultiBuffer`.
  
-		Any window length and block interval is supported, but the FFT size may be rounded up to a faster size (by zero-padding).  It uses a Kaiser window modified for perfect-reconstruction, with shape chosen for almost-optimal aliasing (band-separation) performance.
+		Any window length and block interval is supported, but the FFT size may be rounded up to a faster size (by zero-padding).  It uses a heuristically-optimal Kaiser window modified for perfect-reconstruction.
 		
 		\diagram{stft-aliasing-simulated.svg,Simulated bad-case aliasing (random phase-shift for each band) for overlapping ratios}
 

windows.h

@@ -17,9 +17,8 @@ namespace windows {
 	*/
 
 	/** @brief The Kaiser window (almost) maximises the energy in the main-lobe compared to the side-lobes.
-		These can be specified by the shape-parameter (beta) or by the main lobe's bandwidth (with an optional heuristic: see: `.withBandwidth()`):
-		\diagram{kaiser-windows.svg,You can see that the main lobe matches the specified bandwidth.}
-	*/
+		
+		Kaiser windows can be constructing using the shape-parameter (beta) or using the static `with???()` methods.*/
 	class Kaiser {
 		inline static double bessel0(double x) {
 			constexpr double significanceLimit = 1e-6;
@@ -40,29 +39,96 @@ namespace windows {
 	public:
 		/// Set up a Kaiser window with a given shape.  `beta` is `pi*alpha` (since there is ambiguity about shape parameters)
 		Kaiser(double beta) : beta(beta), invB0(1/bessel0(beta)) {}
-
-		/** @brief Returns the Kaiser shape where the main lobe has the specified bandwidth (as a factor of 1/window-length).
-		If `approximateOptimal` is enabled, the main lobe width is _slightly_ wider, following an approximate compromise between total/peak energy ratios on either side of the boundary.
-		\diagram{kaiser-windows-heuristic.svg,Compare this to the diagram at the top\, where the main-lobe ends exactly at the boundary.}
-		*/
-		static double bandwidthToBeta(double bandwidth, bool approximateOptimal=false) {
-			if (approximateOptimal) {
-				// Heuristic based on numerical search
-				bandwidth = bandwidth + 2/(bandwidth*bandwidth); // mostly improves peak/avg, but gives slightly higher peaks for small bandwidth (in exchange for better average case)
+
+		/// @name Bandwidth methods
+		/// @{
+		static Kaiser withBandwidth(double bandwidth, bool heuristicOptimal=false) {
+			return Kaiser(bandwidthToBeta(bandwidth, heuristicOptimal));
+		}
+
+		/** Returns the Kaiser shape where the main lobe has the specified bandwidth (as a factor of 1/window-length).
+		\diagram{kaiser-windows.svg,You can see that the main lobe matches the specified bandwidth.}
+		If `heuristicOptimal` is enabled, the main lobe width is _slightly_ wider, improving both the peak and total energy - see `bandwidthToEnergyDb()` and `bandwidthToPeakDb()`.
+		\diagram{kaiser-windows-heuristic.svg, The main lobe extends to ±bandwidth/2.} */
+		static double bandwidthToBeta(double bandwidth, bool heuristicOptimal=false) {
+			if (heuristicOptimal) { // Heuristic based on numerical search
+				// Good peaks
+				//bandwidth += 8/((bandwidth + 3)*(bandwidth + 3));
+				// Good average
+				//bandwidth += 14/((bandwidth + 2.5)*(bandwidth + 2.5));
+				// Compromise
+				bandwidth += 8/((bandwidth + 3)*(bandwidth + 3)) + 0.25*std::max(3 - bandwidth, 0.0);
 			}
+			bandwidth = std::max(bandwidth, 2.0);
 			double alpha = std::sqrt(bandwidth*bandwidth*0.25 - 1);
 			return alpha*M_PI;
 		}
-		/** @brief Returns a Kaiser with the appropriate shape (according to `bandwidthToBeta()`). */
+
 		static double betaToBandwidth(double beta) {
 			double alpha = beta*(1.0/M_PI);
 			return 2*std::sqrt(alpha*alpha + 1);
 		}
-		/** @brief Returns a Kaiser with the appropriate shape (according to `bandwidthToBeta()`). */
-		static Kaiser withBandwidth(double bandwidth, bool approximateOptimal=false) {
-			return Kaiser(bandwidthToBeta(bandwidth, approximateOptimal));
-		}
+		/// @}
 
+		/// @name Performance methods
+		/// @{
+		/** @brief Total energy ratio (in dB) between side-lobes and the main lobe.
+			\diagram{kaiser-bandwidth-sidelobes-energy.svg,Measured main/side lobe energy ratio.  You can see that the heuristic improves performance for all bandwidth values.}
+			This function uses an approximation which is accurate to ±0.5dB for 2 ⩽ bandwidth ≤ 10, or 1 ⩽ bandwidth ≤ 10 when `heuristicOptimal`is enabled.
+		*/
+		static double bandwidthToEnergyDb(double bandwidth, bool heuristicOptimal=false) {
+			// Horrible heuristic fits
+			if (heuristicOptimal) {
+				if (bandwidth < 3) bandwidth += (3 - bandwidth)*0.5;
+				return 12.9 + -3/(bandwidth + 0.4) - 13.4*bandwidth + (bandwidth < 3)*-9.6*(bandwidth - 3);
+			}
+			return 10.5 + 15/(bandwidth + 0.4) - 13.25*bandwidth + (bandwidth < 2)*13*(bandwidth - 2);
+		}
+		static double energyDbToBandwidth(double energyDb, bool heuristicOptimal=false) {
+			double bw = 1;
+			while (bw < 20 && bandwidthToEnergyDb(bw, heuristicOptimal) > energyDb) {
+				bw *= 2;
+			}
+			double step = bw/2;
+			while (step > 0.0001) {
+				if (bandwidthToEnergyDb(bw, heuristicOptimal) > energyDb) {
+					bw += step;
+				} else {
+					bw -= step;
+				}
+				step *= 0.5;
+			}
+			return bw;
+		}
+		/** @brief Peak ratio (in dB) between side-lobes and the main lobe.
+			\diagram{kaiser-bandwidth-sidelobes-peak.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.}
+			This function uses an approximation which is accurate to ±0.5dB for 2 ⩽ bandwidth ≤ 9, or 0.5 ⩽ bandwidth ≤ 9 when `heuristicOptimal`is enabled.
+		*/
+		static double bandwidthToPeakDb(double bandwidth, bool heuristicOptimal=false) {
+			// Horrible heuristic fits
+			if (heuristicOptimal) {
+				return 14.2 - 20/(bandwidth + 1) - 13*bandwidth + (bandwidth < 3)*-6*(bandwidth - 3) + (bandwidth < 2.25)*5.8*(bandwidth - 2.25);
+			}
+			return 10 + 8/(bandwidth + 2) - 12.75*bandwidth + (bandwidth < 2)*4*(bandwidth - 2);
+		}
+		static double peakDbToBandwidth(double peakDb, bool heuristicOptimal=false) {
+			double bw = 1;
+			while (bw < 20 && bandwidthToPeakDb(bw, heuristicOptimal) > peakDb) {
+				bw *= 2;
+			}
+			double step = bw/2;
+			while (step > 0.0001) {
+				if (bandwidthToPeakDb(bw, heuristicOptimal) > peakDb) {
+					bw += step;
+				} else {
+					bw -= step;
+				}
+				step *= 0.5;
+			}
+			return bw;
+		}
+		/** @} */
+
 		/// Fills an arbitrary container with a Kaiser window
 		template<typename Data>
 		void fill(Data &data, int size) const {
@@ -73,22 +139,21 @@ namespace windows {
 				data[i] = bessel0(beta*arg)*invB0;
 			}
 		}
-
 	};
 
-	/**@brief Forces STFT perfect-reconstruction (WOLA) on an existing window by direct calculation.
+	/** Forces STFT perfect-reconstruction (WOLA) on an existing window, for a given STFT interval.
 	For example, here are perfect-reconstruction versions of the approximately-optimal @ref Kaiser windows:
 	\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.}
 	*/
 	template<typename Data>
-	void forcePerfectReconstruction(Data &data, int size, int stride) {
-		for (int i = 0; i < stride; ++i) {
+	void forcePerfectReconstruction(Data &data, int windowLength, int interval) {
+		for (int i = 0; i < interval; ++i) {
 			double sum2 = 0;
-			for (int index = i; index < size; index += stride) {
+			for (int index = i; index < windowLength; index += interval) {
 				sum2 += data[index]*data[index];
 			}
-			double factor = std::sqrt(1/sum2);
-			for (int index = i; index < size; index += stride) {
+			double factor = 1/std::sqrt(sum2);
+			for (int index = i; index < windowLength; index += interval) {
 				data[index] *= factor;
 			}
 		}