OXYGENHUB
◆ main openmsx / src / unittest / FFTReal_test.cc
7.8 KiB · C++ · last changed chore(): initial commit 1 month ago source blame history raw
e8555c89 Andy McCall 1 
#include "catch.hpp"
2 
#include "FFTReal.hh"
3
4 
#include <array>
5 
#include <cmath>
6
7 
// Compare two arrays with tolerance, return max error
8 
static float compareArrays(const float* a, const float* b, size_t size) {
9 
	float maxError = 0.0f;
10 
	for (size_t i = 0; i < size; ++i) {
11 
		float error = std::abs(a[i] - b[i]);
12 
		maxError = std::max(maxError, error);
13 
	}
14 
	return maxError;
15 
}
16
17 
TEST_CASE("FFTReal: Full mode - known patterns")
18 
{
19 
	constexpr float tolerance = 1e-5f;
20
21 
	SECTION("All zeros") {
22 
		constexpr unsigned FFT_LEN_L2 = 5;
23 
		constexpr unsigned FFT_LEN = 1 << FFT_LEN_L2; // 32
24
25 
		std::array<float, FFT_LEN> input = {};
26 
		std::array<float, FFT_LEN> output, tmpBuf;
27
28 
		auto result = FFTReal<FFT_LEN_L2>::execute(input, output, tmpBuf);
29
30 
		// FFT of all zeros should be all zeros
31 
		for (unsigned i = 0; i < FFT_LEN; ++i) {
32 
			CHECK(std::abs(result[i]) < tolerance);
33 
		}
34 
	}
35
36 
	SECTION("Impulse at t=0") {
37 
		constexpr unsigned FFT_LEN_L2 = 5;
38 
		constexpr unsigned FFT_LEN = 1 << FFT_LEN_L2; // 32
39
40 
		std::array<float, FFT_LEN> input = {};
41 
		input[0] = 1.0f;
42 
		std::array<float, FFT_LEN> output, tmpBuf;
43
44 
		auto result = FFTReal<FFT_LEN_L2>::execute(input, output, tmpBuf);
45
46 
		// FFT of impulse should be constant (DC = 1, all others = 0 imaginary)
47 
		// Real parts should all be 1.0
48 
		CHECK(std::abs(result[0] - 1.0f) < tolerance); // DC
49 
		for (unsigned i = 1; i <= FFT_LEN / 2; ++i) {
50 
			CHECK(std::abs(result[i] - 1.0f) < tolerance);
51 
		}
52 
		// Imaginary parts should all be ~0
53 
		for (unsigned i = FFT_LEN / 2 + 1; i < FFT_LEN; ++i) {
54 
			CHECK(std::abs(result[i]) < tolerance);
55 
		}
56 
	}
57
58 
	SECTION("Constant signal") {
59 
		constexpr unsigned FFT_LEN_L2 = 5;
60 
		constexpr unsigned FFT_LEN = 1 << FFT_LEN_L2; // 32
61
62 
		std::array<float, FFT_LEN> input;
63 
		for (unsigned i = 0; i < FFT_LEN; ++i) {
64 
			input[i] = 1.0f;
65 
		}
66 
		std::array<float, FFT_LEN> output, tmpBuf;
67
68 
		auto result = FFTReal<FFT_LEN_L2>::execute(input, output, tmpBuf);
69
70 
		// FFT of constant should have all energy in DC bin
71 
		CHECK(std::abs(result[0] - float(FFT_LEN)) < tolerance);
72 
		for (unsigned i = 1; i < FFT_LEN; ++i) {
73 
			CHECK(std::abs(result[i]) < tolerance);
74 
		}
75 
	}
76 
}
77
78 
TEST_CASE("FFTReal: Full mode - fixed random pattern")
79 
{
80 
	// This test uses a fixed random pattern (seed=42, size=32) and expected
81 
	// output calculated with numpy.fft.rfft([...])
82 
	constexpr float tolerance = 1e-4f;
83 
	constexpr unsigned FFT_LEN_L2 = 5;
84 
	constexpr unsigned FFT_LEN = 1 << FFT_LEN_L2; // 32
85
86 
	// Input: 32 random values (seed=42, range [-1, 1])
87 
	// Generated with: np.random.seed(42); np.random.uniform(-1, 1, 32)
88 
	std::array<float, FFT_LEN> input = {
89 
		-0.2509198f,  0.9014286f,  0.4639879f,  0.1973170f,
90 
		-0.6879627f, -0.6880109f, -0.8838328f,  0.7323523f,
91 
		 0.2022300f,  0.4161451f, -0.9588310f,  0.9398197f,
92 
		 0.6648853f, -0.5753218f, -0.6363501f, -0.6331910f,
93 
		-0.3915155f,  0.0495129f, -0.1361100f, -0.4175417f,
94 
		 0.2237058f, -0.7210123f, -0.4157107f, -0.2672763f,
95 
		-0.0878600f,  0.5703519f, -0.6006525f,  0.0284689f,
96 
		 0.1848291f, -0.9070992f,  0.2150897f, -0.6589518f
97 
	};
98
99 
	std::array<float, FFT_LEN> output, tmpBuf;
100 
	auto result = FFTReal<FFT_LEN_L2>::execute(input, output, tmpBuf);
101
102 
	// Expected output from numpy.fft.rfft() - stored in FFTReal format
103 
	// NOTE: FFTReal uses exp(+2πi) convention, so imaginary parts are negated from numpy
104 
	// Real parts: indices 0 to N/2 (17 values)
105 
	// Imaginary parts: indices N/2+1 to N-1 (15 values, negated from numpy)
106 
	std::array<float, FFT_LEN> expected = {
107 
		// Real parts [0 .. 16]
108 
		-4.1280258f,  1.2382417f, -0.9404192f,  2.6338796f,
109 
		 1.3872274f,  1.4370065f, -2.2326224f, -4.0005569f,
110 
		 2.8098016f, -0.2843922f, -3.1922000f, -0.5017270f,
111 
		-3.2142730f,  0.9296157f,  3.3380206f, -0.3273013f,
112 
		-2.0620086f,
113 
		// Imaginary parts [17 .. 31] (negated from numpy values)
114 
		 1.4368339f, -0.9718257f,  0.3860649f,  5.0175196f,
115 
		 1.1372134f,  4.0448431f, -0.4110472f, -0.8750027f,
116 
		 1.2779141f, -0.0101518f, -0.2213625f,  4.0391229f,
117 
		 1.4086333f,  0.2290648f,  3.1862191f
118 
	};
119
120 
	float error = compareArrays(result.data(), expected.data(), FFT_LEN);
121 
	CHECK(error < tolerance);
122 
}
123
124 
TEST_CASE("FFTReal: HalfZeroPadded mode")
125 
{
126 
	// Test that HalfZeroPadded mode (N/2 inputs) produces same result as Full mode
127 
	// with second half explicitly zero-filled
128 
	constexpr float tolerance = 1e-5f;
129 
	constexpr unsigned FFT_LEN_L2 = 8;
130 
	constexpr unsigned FFT_LEN = 1 << FFT_LEN_L2; // 256
131
132 
	// First half: some arbitrary pattern
133 
	std::array<float, FFT_LEN / 2> halfInput = {
134 
		0.5f, -0.3f, 0.8f, -0.2f, 0.1f, 0.9f, -0.7f, 0.4f,
135 
		0.2f, -0.6f, 0.3f, -0.1f, 0.7f, -0.5f, 0.6f, -0.4f,
136 
		0.1f,  0.3f, -0.2f, 0.5f, -0.3f, 0.4f, -0.6f, 0.2f,
137 
		-0.1f, 0.7f, -0.4f, 0.6f, -0.5f, 0.3f, -0.7f, 0.8f,
138 
		0.4f, -0.2f, 0.9f, -0.1f, 0.5f, -0.8f, 0.3f, -0.6f,
139 
		0.7f, -0.4f, 0.2f, -0.9f, 0.1f, -0.3f, 0.8f, -0.5f,
140 
		-0.2f, 0.6f, -0.7f, 0.4f, -0.1f, 0.9f, -0.3f, 0.5f,
141 
		0.8f, -0.6f, 0.2f, -0.4f, 0.7f, -0.9f, 0.3f, -0.1f,
142 
		0.5f, -0.7f, 0.4f, -0.2f, 0.9f, -0.3f, 0.6f, -0.8f,
143 
		0.1f, -0.5f, 0.3f, -0.6f, 0.8f, -0.2f, 0.7f, -0.4f,
144 
		-0.3f, 0.9f, -0.1f, 0.4f, -0.7f, 0.2f, -0.5f, 0.6f,
145 
		0.3f, -0.8f, 0.5f, -0.1f, 0.4f, -0.6f, 0.9f, -0.2f,
146 
		0.7f, -0.4f, 0.1f, -0.9f, 0.3f, -0.5f, 0.8f, -0.6f,
147 
		-0.2f, 0.4f, -0.7f, 0.9f, -0.3f, 0.6f, -0.1f, 0.5f,
148 
		0.8f, -0.2f, 0.4f, -0.6f, 0.7f, -0.9f, 0.1f, -0.3f,
149 
		0.5f, -0.8f, 0.6f, -0.4f, 0.2f, -0.7f, 0.9f, -0.1f
150 
	};
151
152 
	// Full input: same pattern + zeros
153 
	std::array<float, FFT_LEN> fullInput = {};
154 
	for (unsigned i = 0; i < FFT_LEN / 2; ++i) {
155 
		fullInput[i] = halfInput[i];
156 
	}
157
158 
	std::array<float, FFT_LEN> outputFull, tmpBufFull;
159 
	std::array<float, FFT_LEN> outputHalf, tmpBufHalf;
160
161 
	auto resultFull = FFTReal<FFT_LEN_L2, FFTInputMode::Full>::execute(fullInput, outputFull, tmpBufFull);
162 
	auto resultHalf = FFTReal<FFT_LEN_L2, FFTInputMode::HalfZeroPadded>::execute(halfInput, outputHalf, tmpBufHalf);
163
164 
	float error = compareArrays(resultFull.data(), resultHalf.data(), FFT_LEN);
165 
	CHECK(error < tolerance);
166 
}
167
168 
TEST_CASE("FFTReal: Symmetric mode")
169 
{
170 
	// Test that Symmetric mode (N/2+1 inputs) produces same result as Full mode
171 
	// with mirrored second half. Output should be purely real (N/2+1 values).
172 
	constexpr float tolerance = 1e-4f;
173 
	constexpr unsigned FFT_LEN_L2 = 8;
174 
	constexpr unsigned FFT_LEN = 1 << FFT_LEN_L2; // 256
175 
	constexpr unsigned OUTPUT_LEN = FFT_LEN / 2 + 1; // 129
176
177 
	// Symmetric input (first half + center point)
178 
	std::array<float, FFT_LEN / 2 + 1> symInput = {
179 
		0.5f, -0.3f, 0.8f, -0.2f, 0.1f, 0.9f, -0.7f, 0.4f,
180 
		0.2f, -0.6f, 0.3f, -0.1f, 0.7f, -0.5f, 0.6f, -0.4f,
181 
		0.1f,  0.3f, -0.2f, 0.5f, -0.3f, 0.4f, -0.6f, 0.2f,
182 
		-0.1f, 0.7f, -0.4f, 0.6f, -0.5f, 0.3f, -0.7f, 0.8f,
183 
		0.4f, -0.2f, 0.9f, -0.1f, 0.5f, -0.8f, 0.3f, -0.6f,
184 
		0.7f, -0.4f, 0.2f, -0.9f, 0.1f, -0.3f, 0.8f, -0.5f,
185 
		-0.2f, 0.6f, -0.7f, 0.4f, -0.1f, 0.9f, -0.3f, 0.5f,
186 
		0.8f, -0.6f, 0.2f, -0.4f, 0.7f, -0.9f, 0.3f, -0.1f,
187 
		0.5f, -0.7f, 0.4f, -0.2f, 0.9f, -0.3f, 0.6f, -0.8f,
188 
		0.1f, -0.5f, 0.3f, -0.6f, 0.8f, -0.2f, 0.7f, -0.4f,
189 
		-0.3f, 0.9f, -0.1f, 0.4f, -0.7f, 0.2f, -0.5f, 0.6f,
190 
		0.3f, -0.8f, 0.5f, -0.1f, 0.4f, -0.6f, 0.9f, -0.2f,
191 
		0.7f, -0.4f, 0.1f, -0.9f, 0.3f, -0.5f, 0.8f, -0.6f,
192 
		-0.2f, 0.4f, -0.7f, 0.9f, -0.3f, 0.6f, -0.1f, 0.5f,
193 
		0.8f, -0.2f, 0.4f, -0.6f, 0.7f, -0.9f, 0.1f, -0.3f,
194 
		0.5f, -0.8f, 0.6f, -0.4f, 0.2f, -0.7f, 0.9f, -0.1f,
195 
		0.3f  // Center point (index 128)
196 
	};
197
198 
	// Full input: mirror to create symmetric array
199 
	std::array<float, FFT_LEN> fullInput;
200 
	for (unsigned i = 0; i <= FFT_LEN / 2; ++i) {
201 
		fullInput[i] = symInput[i];
202 
	}
203 
	for (unsigned i = FFT_LEN / 2 + 1; i < FFT_LEN; ++i) {
204 
		fullInput[i] = fullInput[FFT_LEN - i];
205 
	}
206
207 
	std::array<float, FFT_LEN> outputFull, tmpBufFull;
208 
	std::array<float, FFT_LEN> outputSym, tmpBufSym;
209
210 
	auto resultFull = FFTReal<FFT_LEN_L2, FFTInputMode::Full>::execute(fullInput, outputFull, tmpBufFull);
211 
	auto resultSym  = FFTReal<FFT_LEN_L2, FFTInputMode::Symmetric>::execute(symInput, outputSym, tmpBufSym);
212
213 
	float error = compareArrays(resultFull.data(), resultSym.data(), OUTPUT_LEN);
214 
	CHECK(error < tolerance);
215 
}