Sample rate normalised most noise powers

Updated the gain of white, pink, red, violet, and blue noise sources so that they depend on sample rate, and the perceived loudness of each source is independent as much as possible between sample rates. Left grey noise alone since I don't understand the algorithm, but it also needs fixing for sample rate independence.

src/Noise.cpp

@@ -4,20 +4,21 @@
 /** Based on "The Voss algorithm"
 http://www.firstpr.com.au/dsp/pink-noise/
 */
-template <int QUALITY = 8>
+template <int MAX_OCTAVES = 16>
 struct PinkNoiseGenerator {
 	int frame = -1;
-	float values[QUALITY] = {};
+	float values[MAX_OCTAVES] = {};
+	int octaves = 8;
 
 	float process() {
 		int lastFrame = frame;
 		frame++;
-		if (frame >= (1 << QUALITY))
+		if (frame >= (1 << octaves))
 			frame = 0;
 		int diff = lastFrame ^ frame;
 
 		float sum = 0.f;
-		for (int i = 0; i < QUALITY; i++) {
+		for (int i = 0; i < octaves; i++) {
 			if (diff & (1 << i)) {
 				values[i] = random::uniform() - 0.5f;
 			}
@@ -38,7 +39,7 @@ struct InverseAWeightingFFTFilter {
 
 	InverseAWeightingFFTFilter() : fft(BUFFER_LEN) {}
 
-	float process(float deltaTime, float x) {
+	float process(float sampleRate, float x) {
 		inputBuffer[frame] = x;
 		if (++frame >= BUFFER_LEN) {
 			frame = 0;
@@ -46,7 +47,7 @@ struct InverseAWeightingFFTFilter {
 			fft.rfft(inputBuffer, freqBuffer);
 
 			for (int i = 0; i < BUFFER_LEN; i++) {
-				float f = 1 / deltaTime / 2 / BUFFER_LEN * i;
+				float f = sampleRate / 2 / BUFFER_LEN * i;
 				float amp = 0.f;
 				if (80.f <= f && f <= 20000.f) {
 					float f2 = f * f;
@@ -64,6 +65,12 @@ struct InverseAWeightingFFTFilter {
 };
 
 
+template <typename T>
+static T softclip(T x) {
+	// rational approximation of 12*tanh(x/12) with derivative 0 at x = 24
+	x = simd::clamp(x, -24.f, 24.f);
+	return x * (3888 + x * x) / (3888 + 9 * x * x);
+}
 
 
 struct Noise : Module {
@@ -88,10 +95,18 @@ struct Noise : Module {
 	};
 
 	dsp::ClockDivider blackDivider;
-	PinkNoiseGenerator<8> pinkNoiseGenerator;
+	PinkNoiseGenerator<16> pinkNoiseGenerator;
+	dsp::IIRFilter<2, 2> whiteFilter;
+	dsp::IIRFilter<2, 2> pinkFilter;
 	dsp::IIRFilter<2, 2> redFilter;
+	dsp::IIRFilter<2, 2> violetFilter;
+	dsp::IIRFilter<2, 2> blueFilter;
 	float lastWhite = 0.f;
 	float lastPink = 0.f;
+	float gainWhite = 0.f;
+	float gainWhiteDifferenced = 0.f;
+	float gainPink = 0.f;
+	float gainPinkDifferenced = 0.f;
 	InverseAWeightingFFTFilter grayFilter;
 
 	// For calibrating levels
@@ -114,52 +129,82 @@ struct Noise : Module {
 		configOutput(BLACK_OUTPUT, "Black noise");
 		outputInfos[BLACK_OUTPUT]->description = "Uniform random numbers";
 
-		// Hard-code coefficients for Butterworth lowpass with cutoff 20 Hz @ 44.1kHz.
-		const float b[] = {0.00425611, 0.00425611};
-		const float a[] = {-0.99148778};
-		redFilter.setCoefficients(b, a);
+		SampleRateChangeEvent e;
+		e.sampleRate = 44100.f;
+		e.sampleTime = 1.f/44100.f;
+		onSampleRateChange(e);
 	}
 
-	void process(const ProcessArgs& args) override {
+	void onSampleRateChange(const SampleRateChangeEvent& e) override {
 		// All noise is calibrated to 1 RMS.
 		// Then they should be scaled to match the RMS of a sine wave with 5V amplitude.
-		const float gain = 5.f / std::sqrt(2.f);
+		// Also normalised by bandwidth of sample rate
+		gainWhite = 5.f * std::sqrt(e.sampleRate)/(std::sqrt(44100.f)*std::sqrt(2.f));
+		gainWhiteDifferenced = gainWhite*e.sampleRate/44100.f;
+		gainPink = 5.f/std::sqrt(2.f);
+		gainPinkDifferenced = gainPink*e.sampleRate/44100.f;
+
+		// DF1 Butterworth lowpass with cutoff 18000 Hz to keep amplitude down
+		{
+			const float g = std::tan(M_PI*fmin(18000.f*e.sampleTime, 0.49f));
+			const float gnrm = 1.f/(1.f + g);
+			const float b[] = {g*gnrm, g*gnrm};
+			const float a[] = {(g - 1.f)*gnrm};
+			whiteFilter.setCoefficients(b, a);
+			pinkFilter.setCoefficients(b, a);
+			violetFilter.setCoefficients(b, a);
+			blueFilter.setCoefficients(b, a);
+		}
+
+		// DF1 Butterworth lowpass with cutoff 60 Hz
+		{
+			const float g = std::tan(M_PI*fmin(60.f*e.sampleTime, 0.49f));
+			const float gnrm = 1.f/(1.f + g);
+			const float b[] = {g*gnrm, g*gnrm};
+			const float a[] = {(g - 1.f)*gnrm};
+			redFilter.setCoefficients(b, a);
+		}
+		const int octaves = clamp(int(std::log2(e.sampleRate) + 1) - 8, 4, 16);
+		pinkNoiseGenerator.octaves = octaves;
+	}
+
+	void process(const ProcessArgs& args) override {
 
 		if (outputs[WHITE_OUTPUT].isConnected() || outputs[RED_OUTPUT].isConnected() || outputs[VIOLET_OUTPUT].isConnected() || outputs[GRAY_OUTPUT].isConnected()) {
 			// White noise: equal power density
-			float white = random::normal();
-			outputs[WHITE_OUTPUT].setVoltage(white * gain);
+			float white = whiteFilter.process(random::normal());
+			outputs[WHITE_OUTPUT].setVoltage(softclip(white * gainWhite));
 
 			// Red/Brownian noise: -6dB/oct
 			if (outputs[RED_OUTPUT].isConnected()) {
 				float red = redFilter.process(white) / 0.0645f;
-				outputs[RED_OUTPUT].setVoltage(red * gain);
+				outputs[RED_OUTPUT].setVoltage(softclip(red * gainWhite));
 			}
 
 			// Violet/purple noise: 6dB/oct
 			if (outputs[VIOLET_OUTPUT].isConnected()) {
-				float violet = (white - lastWhite) / 1.41f;
+				float violet = violetFilter.process((white - lastWhite) / 1.41f);
 				lastWhite = white;
-				outputs[VIOLET_OUTPUT].setVoltage(violet * gain);
+				outputs[VIOLET_OUTPUT].setVoltage(softclip(violet * gainWhiteDifferenced));
 			}
 
 			// Gray noise: psychoacoustic equal loudness curve, specifically inverted A-weighted
 			if (outputs[GRAY_OUTPUT].isConnected()) {
-				float gray = grayFilter.process(args.sampleTime, white) / 1.67f;
-				outputs[GRAY_OUTPUT].setVoltage(gray * gain);
+				float gray = grayFilter.process(args.sampleRate, white) / 1.67f;
+				outputs[GRAY_OUTPUT].setVoltage(softclip(gray * gainWhite));
 			}
 		}
 
 		if (outputs[PINK_OUTPUT].isConnected() || outputs[BLUE_OUTPUT].isConnected()) {
 			// Pink noise: -3dB/oct
-			float pink = pinkNoiseGenerator.process() / 0.816f;
-			outputs[PINK_OUTPUT].setVoltage(pink * gain);
+			float pink = pinkFilter.process(pinkNoiseGenerator.process() / 0.816f);
+			outputs[PINK_OUTPUT].setVoltage(softclip(pink * gainPink));
 
 			// Blue noise: 3dB/oct
 			if (outputs[BLUE_OUTPUT].isConnected()) {
-				float blue = (pink - lastPink) / 0.705f;
+				float blue = blueFilter.process((pink - lastPink) / 0.705f);
 				lastPink = pink;
-				outputs[BLUE_OUTPUT].setVoltage(blue * gain);
+				outputs[BLUE_OUTPUT].setVoltage(softclip(blue * gainPinkDifferenced));
 			}
 		}