OXYGENHUB
◆ main openmsx / src / unittest / yin_test.cc
8.2 KiB · C++ · last changed chore(): initial commit 1 month ago source blame history raw
1 
#include "catch.hpp"
2 
#include "yin.hh"
3
4 
#include <cmath>
5 
#include <cstdlib>
6 
#include <numbers>
7 
#include <random>
8 
#include <vector>
9
10 
static constexpr float PI = std::numbers::pi_v<float>;
11 
static constexpr float TAU = 2.0f * PI;
12
13 
static constexpr float SAMPLE_RATE = 44100.0f;
14 
static constexpr size_t BUFFER_SIZE = 4096; // ~93ms at 44.1kHz
15 
static constexpr float MIN_FREQ = 50.0f;
16 
static constexpr float MAX_FREQ = 2000.0f;
17
18 
[[nodiscard]] static bool approx_rel(float value, float expected, float rel_err)
19 
{
20 
	return std::abs(value - expected) < (expected * rel_err);
21 
}
22 
[[nodiscard]] static bool approx_abs(float value, float expected, float abs_err)
23 
{
24 
	return std::abs(value - expected) < abs_err;
25 
}
26
27 
static void addSine(std::span<float> signal, float frequency, float amplitude = 1.0f, size_t start = 0, size_t end = size_t(-1))
28 
{
29 
	if (end == size_t(-1)) end = signal.size();
30 
	auto C = TAU * frequency / SAMPLE_RATE;
31 
	for (size_t i = start; i < end; ++i) {
32 
		signal[i] += amplitude * std::sin(float(i) * C);
33 
	}
34 
};
35 
static void addNoise(std::span<float> signal, float amplitude)
36 
{
37 
	std::mt19937 rng(42); // fixed seed
38 
	std::normal_distribution<float> dist(0.0f, 1.0f); // mean, stddev
39 
	for (size_t i = 0; i < signal.size(); ++i) {
40 
		signal[i] += amplitude * dist(rng);
41 
	}
42 
};
43
44 
[[nodiscard]] static auto run(std::span<const float> signal)
45 
{
46 
	return yin::detectPitch(signal, SAMPLE_RATE, MIN_FREQ, MAX_FREQ);
47 
}
48
49 
TEST_CASE("YIN: silence")
50 
{
51 
	auto signal = std::vector<float>(BUFFER_SIZE, 0.0f);
52
53 
	auto [freq, err] = run(signal);
54 
	CHECK(approx_abs(freq, 0.0f, 1.0f));
55 
	CHECK(approx_rel(err, 1.0f, 0.05f));
56 
}
57
58 
TEST_CASE("YIN: sine")
59 
{
60 
	auto generateSine = [](float frequency) {
61 
		std::vector<float> signal(BUFFER_SIZE);
62 
		addSine(signal, frequency);
63 
		return signal;
64 
	};
65 
	SECTION("50Hz") { // MIN_FREQ
66 
		auto [freq, err] = run(generateSine(50.0f));
67 
		CHECK(approx_rel(freq, 50.0f, 0.01f));
68 
		CHECK(err < 0.01f);
69 
	}
70 
	SECTION("100Hz") {
71 
		auto [freq, err] = run(generateSine(100.0f));
72 
		CHECK(approx_rel(freq, 100.0f, 0.01f));
73 
		CHECK(err < 0.01f);
74 
	}
75 
	SECTION("440Hz") {
76 
		auto [freq, err] = run(generateSine(440.0f));
77 
		CHECK(approx_rel(freq, 440.0f, 0.01f));
78 
		CHECK(err < 0.01f);
79 
	}
80 
	SECTION("1000Hz") {
81 
		auto [freq, err] = run(generateSine(1000.0f));
82 
		CHECK(approx_rel(freq, 1000.0f, 0.01f));
83 
		CHECK(err < 0.01f);
84 
	}
85 
	SECTION("2000Hz") { // MAX_FREQ
86 
		auto [freq, err] = run(generateSine(2000.0f));
87 
		CHECK(approx_rel(freq, 2000.0f, 0.01f));
88 
		CHECK(err < 0.01f);
89 
	}
90 
}
91
92 
TEST_CASE("YIN: square")
93 
{
94 
	auto generateSquare = [](float frequency) {
95 
		// Fourier series: square wave = sum of odd harmonics
96 
		// Square(t) = (4/π) * [sin(ωt) + (1/3)sin(3ωt) + (1/5)sin(5ωt) + ...]
97 
		std::vector<float> signal(BUFFER_SIZE, 0.0f);
98 
		for (int harmonic = 1; harmonic <= 11; harmonic += 2) {
99 
			float amplitude = 4.0f / (PI * harmonic);
100 
			addSine(signal, frequency * harmonic, amplitude);
101 
		}
102 
		return signal;
103 
	};
104 
	auto [freq, err] = run(generateSquare(321.0f));
105 
	CHECK(approx_rel(freq, 321.0f, 0.01f));
106 
	CHECK(err < 0.01f);
107 
}
108
109 
TEST_CASE("YIN: triangle")
110 
{
111 
	auto generateTriangle = [](float frequency) {
112 
		// Triangle wave = sum of odd harmonics with 1/n² amplitude
113 
		std::vector<float> signal(BUFFER_SIZE, 0.0f);
114 
		for (int harmonic = 1; harmonic <= 11; harmonic += 2) {
115 
			float amplitude = 8.0f / (PI * PI * harmonic * harmonic);
116 
			if ((harmonic - 1) / 2 % 2 == 1) amplitude = -amplitude; // alternating sign
117 
			addSine(signal, frequency * harmonic, amplitude);
118 
		}
119 
		return signal;
120 
	};
121 
	auto [freq, err] = run(generateTriangle(555.0f));
122 
	CHECK(approx_rel(freq, 555.0f, 0.01f));
123 
	CHECK(err < 0.01f);
124 
}
125
126 
TEST_CASE("YIN: sawtooth")
127 
{
128 
	auto generateSawtooth = [](float frequency) {
129 
		// Sawtooth = sum of all harmonics with 1/n amplitude
130 
		std::vector<float> signal(BUFFER_SIZE, 0.0f);
131 
		for (int harmonic = 1; harmonic <= 20; ++harmonic) {
132 
			float amplitude = 2.0f / (PI * harmonic);
133 
			addSine(signal, frequency * harmonic, amplitude);
134 
		}
135 
		return signal;
136 
	};
137 
	auto [freq, err] = run(generateSawtooth(222.0f));
138 
	CHECK(approx_rel(freq, 222.0f, 0.01f));
139 
	CHECK(err < 0.01f);
140 
}
141
142 
TEST_CASE("YIN: sine with harmonics")
143 
{
144 
	SECTION("Sine 440 Hz + decaying harmonics (2nd, 3rd)") {
145 
		std::vector<float> signal(BUFFER_SIZE, 0.0f);
146 
		addSine(signal, 1.0f * 440.0f, 1.0f);
147 
		addSine(signal, 2.0f * 440.0f, 0.5f);
148 
		addSine(signal, 3.0f * 440.0f, 0.25f);
149
150 
		auto [freq, err] = run(signal);
151 
		CHECK(approx_rel(freq, 440.0f, 0.01f));
152 
		CHECK(err < 0.01f);
153 
	}
154 
	SECTION("Sine 440 Hz + strong harmonics") {
155 
		std::vector<float> signal(BUFFER_SIZE, 0.0f);
156 
		addSine(signal, 1.0f * 440.0f, 1.0f);
157 
		addSine(signal, 2.0f * 440.0f, 2.0f);
158 
		addSine(signal, 3.0f * 440.0f, 1.5f);
159
160 
		auto [freq, err] = run(signal);
161 
		CHECK(approx_rel(freq, 440.0f, 0.01f));
162 
		CHECK(err < 0.01f);
163 
	}
164 
	SECTION("Sine 440 Hz + 880Hz equally strong") {
165 
		std::vector<float> signal(BUFFER_SIZE, 0.0f);
166 
		addSine(signal, 1.0f * 440.0f, 1.0f);
167 
		addSine(signal, 2.0f * 440.0f, 1.0f);
168
169 
		auto [freq, err] = run(signal);
170 
		CHECK(approx_rel(freq, 440.0f, 0.01f));
171 
		CHECK(err < 0.01f);
172 
	}
173 
}
174
175 
TEST_CASE("YIN: noise")
176 
{
177 
	std::mt19937 rng(12345); // fixed seed
178 
	SECTION("uniform") {
179 
		std::uniform_real_distribution<float> dist(-1.0f, 1.0f);
180 
		std::vector<float> signal(BUFFER_SIZE);
181 
		std::ranges::generate(signal, [&]{ return dist(rng); });
182
183 
		auto [freq, err] = run(signal);
184 
		// freq doesn't matter
185 
		CHECK(err > 0.8f);
186 
	}
187 
	SECTION("normal") {
188 
		std::normal_distribution<float> dist(0.0f, 1.0f); // mean, stddev
189 
		std::vector<float> signal(BUFFER_SIZE);
190 
		std::ranges::generate(signal, [&]{ return dist(rng); });
191
192 
		auto [freq, err] = run(signal);
193 
		// freq doesn't matter
194 
		CHECK(err > 0.8f);
195 
	}
196 
}
197
198 
TEST_CASE("YIN: chirp") // frequency sweep
199 
{
200 
	auto generateChirp = [](float startFreq, float endFreq) {
201 
		std::vector<float> signal(BUFFER_SIZE);
202 
		for (size_t i = 0; i < BUFFER_SIZE; ++i) {
203 
			float t = float(i) / SAMPLE_RATE;
204 
			float freq = startFreq + (endFreq - startFreq) * t / (float(BUFFER_SIZE) / SAMPLE_RATE);
205 
			signal[i] = std::sin(TAU * freq * t);
206 
		}
207 
		return signal;
208 
	};
209 
	auto [freq, err] = run(generateChirp(100.0f, 1000.0f));
210 
	// freq doesn't matter
211 
	CHECK(err > 0.8f);
212 
}
213
214 
TEST_CASE("YIN: sine decaying in amplitude")
215 
{
216 
	// Generate sine with decaying amplitude (exponential decay)
217 
	auto generateDecayingSine = [](float frequency, float decayTime) {
218 
		std::vector<float> signal(BUFFER_SIZE);
219 
		float decayConstant = -std::log(0.01f) / decayTime; // decay to 1% over decayTime
220 
		for (size_t i = 0; i < BUFFER_SIZE; ++i) {
221 
			float t = float(i) / SAMPLE_RATE;
222 
			float envelope = std::exp(-decayConstant * t);
223 
			signal[i] = envelope * std::sin(TAU * frequency * t);
224 
		}
225 
		return signal;
226 
	};
227 
	auto [freq, err] = run(generateDecayingSine(440.0f, 0.05f));
228 
	CHECK(approx_rel(freq, 440.0f, 0.01f));
229 
	CHECK(err < 0.05f); // still very good, but less confident than the 'pure' signals
230 
}
231
232 
TEST_CASE("YIN: split frequency")
233 
{
234 
	// split frequency: first half at 100Hz, second half at 123Hz (not a nice ratio)
235 
	std::vector<float> signal(BUFFER_SIZE);
236 
	addSine(signal, 100.0f, 1.0f, 0, BUFFER_SIZE / 2);
237 
	addSine(signal, 123.0f, 1.0f, BUFFER_SIZE / 2);
238
239 
	auto [freq, err] = run(signal);
240 
	CHECK(approx_rel(freq, 113.0f, 0.1f)); // something in between
241 
	CHECK(err < 0.25f); // but not a very confident detection
242 
}
243
244 
TEST_CASE("YIN: sine with noise")
245 
{
246 
	auto generateSineWithNoise = [](float freq, float noiseAmplitude) {
247 
		std::vector<float> signal(BUFFER_SIZE);
248 
		addSine(signal, freq);
249 
		addNoise(signal, noiseAmplitude);
250 
		return signal;
251 
	};
252 
	SECTION("0.1") {
253 
		auto [freq, err] = run(generateSineWithNoise(440.0f, 0.1f));
254 
		CHECK(approx_rel(freq, 440.0f, 0.01f));
255 
		CHECK(err < 0.02f); // still very good, but less confident than the 'pure' signals
256 
	}
257 
	SECTION("0.2") {
258 
		auto [freq, err] = run(generateSineWithNoise(440.0f, 0.2f));
259 
		CHECK(approx_rel(freq, 440.0f, 0.01f));
260 
		CHECK(err < 0.1f); // more noise, less confident
261 
	}
262 
	SECTION("0.3") {
263 
		auto [freq, err] = run(generateSineWithNoise(440.0f, 0.3f));
264 
		// This is actually wrong, detects freq 2x too low, but at the border of confident (so acceptable?)
265 
		CHECK(approx_rel(freq, 220.0f, 0.01f));
266 
		CHECK(err > 0.14f);
267 
	}
268 
	SECTION("0.4") {
269 
		auto [freq, err] = run(generateSineWithNoise(440.0f, 0.4f));
270 
		// freq doesn't matter
271 
		CHECK(err > 0.2f);
272 
	}
273 
	SECTION("0.5") {
274 
		auto [freq, err] = run(generateSineWithNoise(440.0f, 0.5f));
275 
		// freq doesn't matter
276 
		CHECK(err > 0.2f);
277 
	}
278 
}