Test code and output plots comparing two very-close BLAMPs vs a single BLEP, with single-precision float

main.cpp

#include "elliptic-blep.h"
#include "plot.h"
#include "simple-fft.h"

#include <vector>

void generateWaveform(float freq, float shape, std::vector<float> &result, bool sawtooth) {
	float phase = 0;
	
	float gradUp = freq*2/shape;
	float gradDown = freq*-2/(1 - shape);
	float corner1 = 0.5*shape, corner2 = 1 - corner1;
	signalsmith::blep::EllipticBlep<float> blepper;
	
	// Initial gradient change to kick us off
	blepper.add(gradUp, 2);
	
	for (auto &v : result) {
		// Move the oscillator forwards
		float prevPhase = phase;
		phase += freq;
		// And the blepper
		blepper.step();

		if (sawtooth) {
			if (phase >= corner1 && prevPhase < corner1) {
				float samplesInPast = (phase - corner1)/freq;
				// -2 step jump
				blepper.add(-2, 1, samplesInPast);
			}
			if (phase >= 1) {
				phase -= 1;
			}
		} else {
			if (phase >= corner2 && prevPhase < corner2) {
				float samplesInPast = (phase - corner2)/freq;
				blepper.add(gradUp - gradDown, 2, samplesInPast);
			}
			if (phase >= 1) {
				phase -= 1;
				prevPhase -= 1;
			}
			// Discontinuities
			if (phase >= corner1 && prevPhase < corner1) {
				float samplesInPast = (phase - corner1)/freq;
				blepper.add(gradDown - gradUp, 2, samplesInPast);
			}
		}

		// Naive waveform
		if (phase < corner1) {
			v = 2*phase/shape;
		} else if (phase < corner2) {
			v = (2*phase - 1)/(shape - 1);
		} else {
			v = 2*(phase - 1)/shape;
		}
		// Add BLEP
		v += blepper.get();
	}
}

int main() {
	size_t length = 16384;
	float freq = 0.02373020996; // C6 at 44.1kHz

	signalsmith::plot::Figure figure;
	auto &timePlot = figure(0, 0).plot(600, 150);
	timePlot.x.linear(9.9/freq, 14.6/freq);
	timePlot.y.major(0).linear(-1.1, 1.1).minors(-1, 1);
	auto &freqPlot = figure(0, 1).plot(600, 200);
	freqPlot.x.linear(0, 0.5).major(0).minor(0.5, "Nyquist").label("frequency");
	freqPlot.y.linear(-130, 0).major(0).minors(-30, -60, -90, -120).label("dB");
	
	auto &timeLine = timePlot.line();
	auto &freqLine = freqPlot.line();
	
	std::vector<float> signal(length);
	signalsmith::linear::SimpleFFT<float> fft(length);

	using Complex = std::complex<float>;
	std::vector<Complex> fftTime(length), fftFreq(length);
	
	auto plotWaveform = [&](float shape, bool sawtooth){
		generateWaveform(freq, shape, signal, sawtooth);

		// Spectrum analysis
		for (size_t i = 0; i < length; ++i) {
			float r = (i + 0.5f)/length;
			fftTime[i] = signal[i]*(1 - std::cos(2*M_PI*r));
		}
		fft.fft(length, fftTime.data(), fftFreq.data());

		// Plot waveform and spectrum
		for (size_t i = 0; i < length; ++i) {
			timeLine.add(i, signal[i]);
			
			float db = 10*std::log10(std::norm(fftFreq[i]/float(length)) + 1e-30f);
			freqLine.add(i/float(length), db);
		}
	};
	
	for (double r = 0; r <= 1; r += 0.01) {
		float shape = 0.5 + 0.499*std::sin(r*2*M_PI);
		plotWaveform(shape, false);
		figure.toFrame(r*10);
	}
	figure.loopFrame(10);
	
	figure.write("animation.svg");

	figure.toFrame(0);
	figure.clearFrames();
	timePlot.title("shape = 0.5");
	plotWaveform(0.5, false);
	figure.write("shape-5.svg");

	figure.toFrame(0);
	figure.clearFrames();
	timePlot.title("shape = 0.9");
	plotWaveform(0.9, false);
	figure.write("shape-9.svg");

	figure.toFrame(0);
	figure.clearFrames();
	timePlot.title("shape = 0.999");
	plotWaveform(0.999, false);
	figure.write("shape-999.svg");
	
	figure.toFrame(0);
	figure.clearFrames();
	timePlot.title("shape = 0.9999");
	plotWaveform(0.9999, false);
	figure.write("shape-9999.svg");
	
	figure.toFrame(0);
	figure.clearFrames();
	timePlot.title("shape = 0.99999");
	plotWaveform(0.99999, false);
	figure.write("shape-99999.svg");

	figure.toFrame(0);
	figure.clearFrames();
	timePlot.title("shape = sawtooth");
	plotWaveform(1, true);
	figure.write("shape-sawtooth.svg");
}

shape-5.svg

shape-9.svg

shape-999.svg

shape-9999.svg

shape-99999.svg

shape-sawtooth.svg

simple-fft.h

#ifndef SIGNALSMITH_DSP_LINEAR_SIMPLEFFT_H
#define SIGNALSMITH_DSP_LINEAR_SIMPLEFFT_H

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

#include <complex>
#include <vector>

namespace signalsmith { namespace linear {

/// A self-contained, reasonably fast power-of-2 FFT template
template<typename Sample>
struct SimpleFFT {
	using Complex = std::complex<Sample>;
	
	SimpleFFT(int maxSize=0) {
		resize(maxSize);
	}
	
	void resize(int maxSize) {
		twiddles.resize(maxSize/2);
		for (int i = 0; i < maxSize/2; ++i) {
			double twiddlePhase = -2*M_PI*i/maxSize;
			twiddles[i] = {
				Sample(std::cos(twiddlePhase)),
				Sample(std::sin(twiddlePhase))
			};
		}
		working.resize(maxSize);
	}
	
	void fft(int size, const Complex *time, Complex *freq) {
		if (size <= 1) {
			*freq = *time;
			return;
		}
		fftPass<false>(size, 1, time, freq, working.data());
	}

	void ifft(int size, const Complex *freq, Complex *time) {
		if (size <= 1) {
			*time = *freq;
			return;
		}
		fftPass<true>(size, 1, freq, time, working.data());
	}
private:
	std::vector<Complex> twiddles;
	std::vector<Complex> working;

	// Calculate a [size]-point FFT, where each element is a block of [stride] values
	template<bool inverse>
	void fftPass(int size, int stride, const Complex *input, Complex *output, Complex *working) const {
		if (size > 2) {
			// Calculate the two half-size FFTs (odd and even) by doubling the stride
			fftPass<inverse>(size/2, stride*2, input, working, output);
			combine2<inverse>(size, stride, working, output);
		} else {
			// The input can already be considered a 1-point FFT
			combine2<inverse>(size, stride, input, output);
		}
	}
	
	// Combine interleaved even/odd results into a single spectrum
	template<bool inverse>
	void combine2(int size, int stride, const Complex *input, Complex *output) const {
		auto twiddleStep = twiddles.size()*2/size;
		for (int i = 0; i < size/2; ++i) {
			Complex twiddle = twiddles[i*twiddleStep];
			
			const Complex *inputA = input + 2*i*stride;
			const Complex *inputB = input + (2*i + 1)*stride;
			Complex *outputA = output + i*stride;
			Complex *outputB = output + (i + size/2)*stride;
			for (int s = 0; s < stride; ++s) {
				Complex a = inputA[s];
				Complex b = inputB[s]*(inverse ? std::conj(twiddle) : twiddle);
				outputA[s] = a + b;
				outputB[s] = a - b;
			}
		}
	}
};

}} // namespace
#endif // include guard