Add more STFT tools (partial-sums), odds and ends

README.md

@@ -25,4 +25,4 @@ The goal (where possible) is for the actual audio characteristics of the tools (
 
 ### License
 
-This code is [MIT licensed](LICENSE.txt).
+This code is [MIT licensed](LICENSE.txt).  If you'd prefer something else, get in touch.

delay.h

@@ -416,53 +416,17 @@ namespace delay {
 		}
 	};
 
-//	template<typename Sample>
-//	using InterpolatorLagrange3 = InterpolatorLagrangeN<Sample, 3>;
+	template<typename Sample>
+	using InterpolatorLagrange3 = InterpolatorLagrangeN<Sample, 3>;
 	template<typename Sample>
 	using InterpolatorLagrange7 = InterpolatorLagrangeN<Sample, 7>;
 	template<typename Sample>
 	using InterpolatorLagrange19 = InterpolatorLagrangeN<Sample, 19>;
 
-	template<typename Sample>
-	struct InterpolatorLagrange3 {
-		static constexpr int inputLength = 4;
-		static constexpr int latency = 1;
-
-	protected:
-		template<class Data>
-		static Sample fractional(const Data &data, Sample x) {
-			Sample d0 = data[0];
-			Sample forwardFactor = (x + 1);
-			Sample d1 = data[1]*forwardFactor;
-			forwardFactor *= x;
-			Sample d2 = data[2]*forwardFactor;
-			Sample xm1 = (x - 1);
-			forwardFactor *= xm1;
-			Sample d3 = data[3]*forwardFactor;
-
-			Sample backwardsFactor = (x - 2);
-			Sample result = d3*Sample(1.0/6) + d2*backwardsFactor*Sample(1.0/-2);
-			backwardsFactor *= xm1;
-			result += d1*backwardsFactor*Sample(1.0/2);
-			backwardsFactor *= x;
-			result += d0*backwardsFactor*Sample(1.0/-6);
-			return result;
-
-//			Sample a = data[0], b = data[1], c = data[2], d = data[3];
-//
-//			constexpr Sample frac6 = 1/6.0;
-//			Sample frac6_da = frac6*(d - a);
-//			Sample k3 = frac6_da + 0.5*(b - c);
-//			Sample k2 = 0.5*(a + c) - b;
-//			Sample k1 = -0.5*(a + b) + c - frac6_da;
-//			return b + fractional*(k1 + fractional*(k2 + fractional*k3)); // 18 total
-		}
-	};
-
 	template<typename Sample, int n, bool minimumPhase=false>
 	struct InterpolatorKaiserSincN {
 		static constexpr int inputLength = n;
-		static constexpr Sample latency = minimumPhase ? 1 : (n*Sample(0.5) - 1);
+		static constexpr Sample latency = minimumPhase ? 0 : (n*Sample(0.5) - 1);
 
 		int subSampleSteps;
 		std::vector<Sample> coefficients;
@@ -479,12 +443,6 @@ namespace delay {
 			std::vector<Sample> windowedSinc(subSampleSteps*n + 1);
 
 			signalsmith::windows::fillKaiserBandwidth(windowedSinc, windowedSinc.size(), kaiserBandwidth, false);
-			if (minimumPhase) {
-				for (size_t i = 0; i < windowedSinc.size(); ++i) {
-//					double r = (i + 0.5)/windowedSinc.size();
-//					windowedSinc[i] = 0.5 - 0.5*std::cos(r*2*M_PI);
-				}
-			}
 
 			for (size_t i = 0; i < windowedSinc.size(); ++i) {
 				double x = (i - centreIndex)*scaleFactor;
@@ -526,9 +484,6 @@ namespace delay {
 				windowedSinc.resize(subSampleSteps*n + 1);
 				for (size_t i = 0; i < windowedSinc.size(); ++i) {
 					windowedSinc[i] = cepstrum[i].real()*scaling;
-////					 Taper off end
-//					double r = (i + 0.5)/windowedSinc.size();
-//					windowedSinc[i] *= 0.5 + 0.5*std::cos(r*M_PI);
 				}
 			}
 
@@ -583,8 +538,6 @@ namespace delay {
 	class Reader : public Interpolator<Sample> /* so we can get the empty-base-class optimisation */ {
 		using Super = Interpolator<Sample>;
 	public:
-		static constexpr Sample latency = Super::latency;
-
 		Reader () {}
 		/// Pass in a configured interpolator
 		Reader (const Interpolator<Sample> &interpolator) : Super(interpolator) {}

spectral.h

@@ -28,7 +28,7 @@ namespace spectral {
 		using Complex = std::complex<Sample>;
 		MRFFT mrfft{2};
 
-		std::vector<Sample> window;
+		std::vector<Sample> fftWindow;
 		std::vector<Sample> timeBuffer;
 		std::vector<Complex> freqBuffer;
 	public:
@@ -53,10 +53,10 @@ namespace spectral {
 		/// Sets the size, returning the window for modification (initially all 1s)
 		std::vector<Sample> & setSizeWindow(int size) {
 			mrfft.setSize(size);
-			window.resize(size, 1);
+			fftWindow.resize(size, 1);
 			timeBuffer.resize(size);
 			freqBuffer.resize(size);
-			return window;
+			return fftWindow;
 		}
 		/// Sets the FFT size, with a user-defined functor for the window
 		template<class WindowFn>
@@ -66,7 +66,7 @@ namespace spectral {
 			Sample invSize = 1/(Sample)size;
 			for (int i = 0; i < size; ++i) {
 				Sample r = (i + windowOffset)*invSize;
-				window[i] = fn(r);
+				fftWindow[i] = fn(r);
 			}
 		}
 		/// Sets the size (using the default Blackman-Harris window)
@@ -77,6 +77,13 @@ namespace spectral {
 				return 0.35875 + 0.48829*std::cos(phase) + 0.14128*std::cos(phase*2) + 0.1168*std::cos(phase*3);
 			});
 		}
+
+		const std::vector<Sample> & window() {
+			return this->fftWindow;
+		}
+		int size() {
+			return mrfft.size();
+		}
 
 		/// Performs an FFT (with windowing)
 		template<class Input, class Output>
@@ -89,7 +96,12 @@ namespace spectral {
 				}
 			};
 
-			mrfft.fft(WindowedInput{input, window}, output);
+			mrfft.fft(WindowedInput{input, fftWindow}, output);
+		}
+		/// Performs an FFT (no windowing)
+		template<class Input, class Output>
+		void fftRaw(Input &&input, Output &&output) {
+			mrfft.fft(input, output);
 		}
 
 		/// Inverse FFT, with windowing and 1/N scaling
@@ -99,7 +111,7 @@ namespace spectral {
 			int size = mrfft.size();
 			Sample norm = 1/(Sample)size;
 			for (int i = 0; i < size; ++i) {
-				output[i] *= norm*window[i];
+				output[i] *= norm*fftWindow[i];
 			}
 		}
 	};
@@ -218,6 +230,22 @@ namespace spectral {
 		int windowSize() {
 			return _windowSize;
 		}
+		/// Returns the (analysis and synthesis) window
+		decltype(fft.window()) window() {
+			return fft.window();
+		}
+		/// Calculates the effective window for the partially-summed future output (relative to the most recent block)
+		std::vector<Sample> partialSumWindow() {
+			const auto &w = window();
+			std::vector<Sample> result(_windowSize, 0);
+			for (int offset = 0; offset < _windowSize; offset += interval) {
+				for (int i = 0; i < _windowSize - offset; ++i) {
+					Sample value = w[i + offset];
+					result[i] += value*value;
+				}
+			}
+			return result;
+		}
 
 		/// Resets everything - since we clear the output sum, it will take `windowSize` samples to get proper output.
 		void reset() {
@@ -275,6 +303,17 @@ namespace spectral {
 		void analyse(int c, Data &&data) {
 			fft.fft(data, spectrum[c]);
 		}
+		/// Analyse without windowing
+		template<class Data>
+		void analyseRaw(Data &&data) {
+			for (int c = 0; c < channels; ++c) {
+				fft.fftRaw(data[c], spectrum[c]);
+			}
+		}
+		template<class Data>
+		void analyseRaw(int c, Data &&data) {
+			fft.fftRaw(data, spectrum[c]);
+		}
 
 		int bands() const {
 			return fftSize/2;