Source browser
examples/fft_checks.hpp
Read the complete file below or download its original bytes. Line numbers provide stable links to a declaration or processing step.
File fingerprint
1/** @file fft_checks.hpp
2 * @brief FFT, window, and frequency-filter checks against independent signal references.
3 * Small direct DFTs provide an independent numerical reference for the FFT.
4 * Known tone mixtures check pass bands, stop bands, and reconstructed amplitudes.
5 */
6#ifndef SHDR_EXAMPLES_FFT_CHECKS_HPP // Include the reference helpers only once.
7#define SHDR_EXAMPLES_FFT_CHECKS_HPP // Mark this validation header as included.
8#include "runner.hpp" // Select the CPU or requested GPU backend.
9#include "fft_support.hpp" // Build reusable transform and filtering pipelines.
10#include <complex> // Accumulate independent reference transforms in double precision.
11#include <random> // Generate the same input on every backend.
12
13namespace shdr_examples { // Keep validation utilities with the other examples.
14namespace fft_checks_detail { // Separate reference calculations from shader implementations.
15typedef shdr::glsl::vec2 complex_value; // Store real and imaginary values in the shader layout.
16inline std::vector<complex_value> // Return complex spectra for all input channels.
17direct_dft( // Accept the original channel-major sample vector.
18 const std::vector<complex_value>& values, // Accept the original channel-major sample vector.
19 unsigned length) { // Calculate each frequency independently.
20 std::vector<complex_value> result(values.size()); // Keep the same channel-major shape.
21 const double tau = 6.283185307179586476925286766559; // Use a double-precision phase reference.
22 for (unsigned channel = 0; // Start one independent transform for each channel.
23 channel < values.size() / length; // Start one independent transform for each channel.
24 ++channel) { // Transform each channel as an independent signal.
25 for (unsigned bin = 0; bin < length; ++bin) { // Evaluate one unshifted frequency bin.
26 std::complex<double> total(0.0, // Begin the real part of the reference sum at zero.
27 0.0); // Accumulate with more precision than the shaders.
28 for (unsigned position = 0; // Visit every time-domain position for this bin.
29 position < length; // Visit every time-domain position for this bin.
30 ++position) { // Read every value in this reference block.
31 const complex_value value = // Read one original complex sample.
32 values[channel * length + position]; // Select the original complex input.
33 const double phase = // Calculate this sample’s phase in radians.
34 -tau * bin * position / length; // Follow the unscaled forward DFT convention.
35 total += std::complex<double>( // Multiply the complex sample by its unit-magnitude phase.
36 value.x, // Multiply the complex sample by its unit-magnitude phase.
37 value.y) * // Multiply the complex sample by its unit-magnitude phase.
38 std::polar(1.0, phase); // Add this input's complex frequency contribution.
39 } // This reference sums all complex Fourier terms directly.
40 result[channel * length + bin] = // Store this bin in the channel-major reference output.
41 complex_value( // Store this bin in the channel-major reference output.
42 float(total.real()), // Convert the accumulated real component to float.
43 float(total.imag())); // Compare the final reference with float shader outputs.
44 } // Continue with the next frequency in this channel.
45 } // Independent channels have separate reference spectra.
46 return result; // Supply a complete spectrum for verification.
47} // The reference has quadratic work and is used only on small test blocks.
48inline void compare(const shdr::buffer<complex_value>& // Read the completed transform through a const resource.
49 actual, // Read the completed transform through a const resource.
50 const std::vector<complex_value>& // Accept the independent reference values to compare.
51 expected, // Accept the independent reference values to compare.
52 const char* label) { // Check every complex value.
53 const complex_value* values = actual.data(); // Download GPU results only after processing.
54 for (unsigned i = 0; i < expected.size(); ++i) { // Compare both components at each position.
55 require_close( // Check the real component with an absolute and relative tolerance.
56 values[i].x, // Read the shader’s real output.
57 expected[i].x, // Use the direct reference’s real result.
58 label, // Keep the caller’s diagnostic label.
59 0.0003f + // Allow a small absolute float error.
60 0.00001f * // Add tolerance proportional to the expected amplitude.
61 std::fabs( // Use the magnitude so negative values receive the same tolerance.
62 expected[i].x)); // Allow float arithmetic error in the real component.
63 require_close( // Check the imaginary component separately.
64 values[i].y, // Read the shader’s imaginary output.
65 expected[i].y, // Use the direct reference’s imaginary result.
66 label, // Preserve the same diagnostic label.
67 0.0003f + 0.00001f * // Combine absolute and relative imaginary error limits.
68 std::fabs(expected[i].y)); // Use the same tolerance for imaginary values.
69 } // A mismatch identifies a failed transform or signal response.
70} // Const access preserves the input values for subsequent dispatches.
71inline double // Evaluate scalar filter gains with double precision.
72low_gain(double frequency, // Accept the positive frequency in hertz.
73 double cutoff, // Accept the transition center in hertz.
74 double width) { // Describe the expected block-filter response in double precision.
75 if (width == 0.0) // Select the hard-cutoff response for a zero transition width.
76 return frequency <= cutoff ? 1.0 : 0.0; // Match the documented hard cutoff convention.
77 const double position = // Normalize distance from the transition’s lower edge.
78 (frequency - cutoff + width * 0.5) / // Normalize distance from the transition’s lower edge.
79 width; // Express frequency within the cosine transition.
80 if (position <= 0.0) // Detect the full-pass region.
81 return 1.0; // Keep values below the transition.
82 if (position >= 1.0) // Detect the full-stop region.
83 return 0.0; // Reject values above the transition.
84 return 0.5 + // Use a raised cosine within the transition.
85 0.5 * // Use a raised cosine within the transition.
86 std::cos( // Use a raised cosine within the transition.
87 3.1415926535897932384626433832795 * // Use a raised cosine within the transition.
88 position); // Evaluate the smooth gain independently of the shader.
89} // Negative-frequency bins use this same real gain.
90} // namespace fft_checks_detail
91
92/** @brief Check small batched transforms, reconstruction, input aliasing, and invalid sizes. */
93inline void // Expose the transform validation routine.
94check_fft1d(runner& run) { // Use independent reference transforms on the selected backend.
95 const unsigned lengths[] = {// List the small block sizes that receive direct references.
96 1, // Include this documented transform length in the independent reference checks.
97 2, // Include this documented transform length in the independent reference checks.
98 4, // Include this documented transform length in the independent reference checks.
99 8, // Include this documented transform length in the independent reference checks.
100 16, // Include this documented transform length in the independent reference checks.
101 64, // Include this documented transform length in the independent reference checks.
102 128}; // Cover the copy case and several butterfly depths.
103 std::mt19937 random(12345u); // Reproduce inputs on different machines.
104 std::uniform_real_distribution<float> values( // Begin a uniform distribution at minus one.
105 -1.0f, // Begin a uniform distribution at minus one.
106 1.0f); // Keep reference amplitudes well bounded.
107 for (unsigned length : lengths) { // Validate each supported block size.
108 const unsigned channels = 3, // Test three separate channel transforms per length.
109 count = length * channels; // Include independently transformed channels.
110 std::vector<shdr::glsl::vec2> original( // Allocate the original complex fixture vector.
111 count); // Retain unmodified values for reference checks.
112 shdr::buffer<shdr::glsl::vec2> input(count), // Allocate input storage and forward-transform output.
113 spectrum(count), // Allocate input storage and forward-transform output.
114 restored(count); // Allocate reusable transform resources.
115 for (unsigned i = 0; i < count; // Visit every sample across the three channels.
116 ++i) { // Fill real and imaginary components with random input.
117 original[i] = shdr::glsl::vec2( // Assign both complex components of this input.
118 values(random), // Generate the real component independently.
119 values(random)); // Generate both real and imaginary sample components.
120 } // Different channels start with different values.
121 original[0] += shdr::glsl::vec2(3.0f, -2.0f); // Include a strong complex impulse.
122 std::copy(original.begin(), // Start copying at the original vector’s first value.
123 original.end(), // Stop at the end of the original vector.
124 input.data()); // Supply the same original values to either backend.
125 fft1d forward(input, spectrum, length, channels); // Plan an unscaled forward transform.
126 fft1d inverse(spectrum, // Read the forward spectrum as inverse-transform input.
127 restored, // Write reconstructed values into separate storage.
128 length, // Use the same length for both transforms.
129 channels, // Reconstruct all three channels independently.
130 true); // Plan inverse reconstruction with 1/N scaling.
131 shdr::pipeline pipeline; // Link both transforms before running.
132 forward.append_to(pipeline); // First produce the frequency spectrum.
133 inverse.append_to(pipeline); // Reconstruct from that spectrum on the same backend.
134 run.run(pipeline); // Execute every stage in dependency order.
135 fft_checks_detail::compare( // Compare the completed forward spectrum.
136 spectrum, // Compare the completed forward spectrum.
137 fft_checks_detail::direct_dft( // Compute a direct double-precision reference from original samples.
138 original, // Use the preserved original values as the reference.
139 length), // Compute a direct double-precision reference from original samples.
140 "1D FFT reference"); // Validate the actual forward FFT.
141 fft_checks_detail::compare( // Check the reconstructed signal after the inverse stage.
142 restored, // Read the inverse-transform output.
143 original, // Use the original signal as the reconstruction reference.
144 "1D inverse FFT reconstruction"); // Validate normalization and inverse phases.
145 double input_energy = 0.0, // Start the original signal’s energy sum at zero.
146 spectrum_energy = 0.0; // Check the transform's energy convention.
147 const shdr::glsl::vec2* bins = // Access forward bins through a read-only pointer.
148 static_cast<const shdr::buffer<shdr::glsl::vec2>&>( // Keep the spectrum resource const during readback.
149 spectrum) // Keep the spectrum resource const during readback.
150 .data(); // Read completed bins through const resource access.
151 for (unsigned i = 0; i < count; ++i) { // Sum squared magnitudes across all channels.
152 input_energy += // Accumulate each original sample’s squared magnitude.
153 double(original[i].x) * // Square the original real component in double precision.
154 original[i].x + // Square the original real component in double precision.
155 double(original[i].y) * original[i].y; // Measure original signal energy.
156 spectrum_energy += // Accumulate each spectrum bin’s squared magnitude.
157 double(bins[i].x) * // Square the spectrum’s real component in double precision.
158 bins[i].x + // Square the spectrum’s real component in double precision.
159 double(bins[i].y) * bins[i].y; // Measure the unscaled spectrum energy.
160 } // Parseval's identity predicts an N-fold increase for this forward convention.
161 if (std::fabs(spectrum_energy / length - // Compare normalized spectral energy with original energy.
162 input_energy) > // Compare normalized spectral energy with original energy.
163 0.00005 * input_energy) { // Allow the accumulated float error.
164 throw shdr::error( // Report a failed energy identity.
165 "FFT energy normalization differs from Parseval's identity"); // Report an incorrect transform convention.
166 } // The transform's scaling is now independently checked.
167 std::copy(original.begin(), // Begin restoring the original sample vector.
168 original.end(), // Copy its complete range before the aliasing check.
169 input.data()); // Restore the original values before testing aliases.
170 fft1d in_place( // Construct a transform with shared input and output storage.
171 input, // Read samples from the caller’s input buffer.
172 input, // Write results back into that same buffer.
173 length, // Keep the original transform length.
174 channels); // Use working buffers to support a shared input and output resource.
175 run.run(in_place.processing_pipeline()); // Perform the same FFT with caller storage reused.
176 fft_checks_detail::compare( // Compare the aliased transform with an independent result.
177 input, // Read the buffer overwritten by the transform.
178 fft_checks_detail::direct_dft( // Recompute the direct reference from preserved original values.
179 original, length), // Recompute the direct reference from preserved original values.
180 "1D in-place FFT"); // Verify that an in-place transform preserves the expected result.
181 } // All small sizes agree with an independent direct transform.
182 shdr::buffer<shdr::glsl::vec2> short_input( // Allocate a deliberately small input buffer for invalid plans.
183 4), // Allocate a deliberately small input buffer for invalid plans.
184 short_output(4); // Exercise dimension checks with small buffers.
185 for (unsigned length : {0u, 3u, 8u}) { // Reject zero, non-power-of-two, and undersized buffers.
186 bool rejected = false; // Record whether validation caught the invalid plan.
187 try { // Attempt construction so validation failures can be caught.
188 fft1d invalid(short_input, // Ask the public constructor to validate this plan.
189 short_output, // Ask the public constructor to validate this plan.
190 length); // Ask the public constructor to validate this plan.
191 } // Try to construct a plan with this invalid length or capacity.
192 catch (const shdr::error&) { // Catch the library’s dimension or capacity error.
193 rejected = true; // Remember that this unsafe plan was rejected.
194 } // Dimension failures should be named runtime errors.
195 if (!rejected) // Fail if an invalid plan escaped its checks.
196 throw shdr::error( // Report acceptance of an unsafe transform shape.
197 "invalid FFT dimensions were accepted"); // Keep unsafe plans out of application loops.
198 } // Valid plans have consistent lengths and resource capacities.
199 std::cout // Write the transform validation result to standard output.
200 << "1D FFT: batched complex spectra, inverse, energy, and input aliasing verified\n"; // Summarize the reference checks.
201} // The next examples use these transforms as signal-processing stages.
202
203/** @brief Validate all block-frequency filter modes with known tones and cosine transitions. */
204inline void check_frequency_filters( // Expose the frequency-response validation routine.
205 runner& run) { // Compare reconstructed signals with their known frequency response.
206 const unsigned length = 128, channels = 2, // Use two channels of 128 samples each.
207 count = length * channels; // Use bin-aligned tones in two independent channels.
208 const double tau = // Keep a double-precision full-turn phase constant.
209 6.283185307179586476925286766559; // Generate reference sinusoidal inputs in double precision.
210 const double frequencies[] = {// List tone frequencies used to probe the filter edges.
211 3.0, // Include pass, stop, and exact cutoff frequencies.
212 8.0, // Include pass, stop, and exact cutoff frequencies.
213 11.0, // Include pass, stop, and exact cutoff frequencies.
214 20.0, // Include pass, stop, and exact cutoff frequencies.
215 25.0}; // Include pass, stop, and exact cutoff frequencies.
216 const double amplitudes[] = {// Assign a distinguishable amplitude to every tone.
217 0.7, // Give each component an identifiable amplitude.
218 0.2, // Give each component an identifiable amplitude.
219 0.5, // Give each component an identifiable amplitude.
220 0.3, // Give each component an identifiable amplitude.
221 0.4}; // Give each component an identifiable amplitude.
222 shdr::buffer<shdr::glsl::vec2> input(count), // Allocate original channel-major complex samples.
223 output(count); // Keep original and filtered signals separate.
224 for (unsigned channel = 0; channel < channels; // Build each channel with its own phase offset.
225 ++channel) { // Generate a different phase for each channel.
226 for (unsigned position = 0; // Visit every sample in this channel’s block.
227 position < length; // Visit every sample in this channel’s block.
228 ++position) { // Form the complete real-valued tone mixture.
229 double value = 0.3; // Include DC so high-pass rejection is checked.
230 for (unsigned tone = 0; tone < 5; ++tone) { // Add all known components.
231 value += amplitudes[tone] * // Weight this tone by its known amplitude.
232 std::cos(tau * frequencies[tone] * position / // Generate its cosine with a bin-aligned phase.
233 length + // Generate its cosine with a bin-aligned phase.
234 0.31 * channel); // Keep every component aligned with FFT bins.
235 } // The frequencies share the sample rate of 128 hertz.
236 input.data()[channel * length + // Write the complete tone mixture into this sample.
237 position] = // Write the complete tone mixture into this sample.
238 shdr::glsl::vec2(float(value), 0.0f); // Store a real input in complex FFT storage.
239 } // Channel layout matches the transform helper.
240 } // Every filter receives the same original signal.
241 for (unsigned selected = 0; selected < 4; // Select each of the four public filter modes.
242 ++selected) { // Exercise low-pass, high-pass, band-pass, and band-stop modes.
243 for (float width : {0.0f, 6.0f}) { // Check both hard bin edges and smooth transitions.
244 const frequency_filter // Convert the test index into the public filter enumeration.
245 mode = // Convert the test index into the public filter enumeration.
246 static_cast<frequency_filter>(selected); // Use the public mode type.
247 fft_filter filter(input, // Read the original signal in the filtering pipeline.
248 output, // Write reconstructed filtered samples into separate storage.
249 length, // Use 128 samples per transform.
250 channels, // Process both channels with the same frequency response.
251 128.0f, // Set a sample rate of 128 hertz.
252 mode, // Choose the current low, high, band or notch response.
253 8.0f, // Place the lower cutoff at eight hertz.
254 20.0f, // Place the upper band edge at twenty hertz.
255 width); // Build the complete transform-filter-inverse pipeline.
256 run.run(filter.processing_pipeline()); // Keep intermediate spectra in their bound processing resources.
257 std::vector<shdr::glsl::vec2> expected( // Allocate a complete analytic output vector.
258 count); // Assemble the known output from tone amplitudes.
259 const auto response = // Build an independent scalar response function.
260 [&](double frequency) { // Evaluate the documented gain for any component.
261 const double lower = fft_checks_detail:: // Evaluate the low-pass response at the lower cutoff.
262 low_gain( // Evaluate the low-pass response at the lower cutoff.
263 frequency, // Describe the lower cutoff.
264 8.0, // Describe the lower cutoff.
265 width); // Describe the lower cutoff.
266 if (selected == 0) // Handle the low-pass mode.
267 return lower; // Low-pass keeps frequencies below this cutoff.
268 if (selected == 1) // Handle the complementary high-pass mode.
269 return 1.0 - lower; // High-pass uses its complementary response.
270 const double band = // Combine the two cutoff responses into a pass band.
271 (1.0 - lower) * // Reject frequencies retained by the lower low-pass edge.
272 fft_checks_detail::low_gain( // Retain frequencies below the upper edge.
273 frequency, // Combine both edges for a band.
274 20.0, // Combine both edges for a band.
275 width); // Combine both edges for a band.
276 return selected == 2 // Choose the pass band or its complement.
277 ? band // Return the pass-band response for band-pass mode.
278 : 1.0 - band; // Band-stop keeps the complement of band-pass.
279 }; // All responses are symmetric for positive and negative frequencies.
280 for (unsigned channel = 0; // Reconstruct each channel with its original phase.
281 channel < channels; // Reconstruct each channel with its original phase.
282 ++channel) { // Preserve each channel's phase in the expected output.
283 for (unsigned position = 0; // Generate each analytic time-domain sample.
284 position < length; // Generate each analytic time-domain sample.
285 ++position) { // Reconstruct the expected time-domain signal directly.
286 double value = 0.3 * response(0.0); // Apply the expected DC gain.
287 for (unsigned tone = 0; // Apply the filter response to every known tone.
288 tone < 5; // Apply the filter response to every known tone.
289 ++tone) { // Apply the expected gain to each tone.
290 value += // Scale the tone by both amplitude and expected gain.
291 amplitudes[tone] * // Scale the tone by both amplitude and expected gain.
292 response( // Scale the tone by both amplitude and expected gain.
293 frequencies // Scale the tone by both amplitude and expected gain.
294 [tone]) * // Scale the tone by both amplitude and expected gain.
295 std::cos( // Preserve the tone’s phase in the reconstructed sample.
296 tau * frequencies[tone] * position / // Preserve the tone’s phase in the reconstructed sample.
297 length + // Preserve the tone’s phase in the reconstructed sample.
298 0.31 * channel); // Keep each tone's original phase.
299 } // A real symmetric response has an expected imaginary component of zero.
300 expected[channel * length + // Store the analytic channel-major output value.
301 position] = // Store the analytic channel-major output value.
302 shdr::glsl::vec2(float(value), 0.0f); // Store the analytic filtered signal.
303 } // Continue through this channel's reconstructed block.
304 } // Apply the same response independently to both channels.
305 fft_checks_detail:: // Compare the actual filtered signal with its analytic response.
306 compare( // Compare the actual filtered signal with its analytic response.
307 output, // Read the reconstructed filter output.
308 expected, // Use the independently assembled tone mixture.
309 "FFT signal filter response"); // Compare every reconstructed value.
310 } // The response checks include the transition edges.
311 } // All four block-filter modes now have numerical coverage.
312 std::cout // Write the filter validation result to standard output.
313 << "FFT filters: low-pass, high-pass, band-pass, and band-stop tone responses verified\n"; // Describe the completed signal checks.
314} // Circular block filtering remains distinct from a streaming FIR or IIR filter.
315
316/** @brief Check rectangular, Hann, Hamming, and Blackman windows before a Fourier transform. */
317inline void check_signal_windows(runner& run) { // Verify window values and their linked spectra.
318 const unsigned length = 32; // Use a small block for an independent direct DFT.
319 shdr::buffer<shdr::glsl::vec2> input(length), // Allocate original and window-weighted complex values.
320 windowed(length), // Allocate original and window-weighted complex values.
321 spectrum(length); // Keep the original signal and both outputs.
322 for (unsigned i = 0; i < length; ++i) // Create every sample in the small window fixture.
323 input.data()[i] = // Assign the next complex fixture value.
324 shdr::glsl::vec2(float(i % 7) - 3.0f, 0.25f); // Include both complex components.
325 auto window = shdr::make_kernel("signal_window_kernel"); // Create the native window stage once.
326 window // Connect original samples to the window shader.
327 ->bind("input_values", input) // Connect original samples to the window shader.
328 .bind("output_values", windowed); // Link caller storage to windowed FFT input.
329 window->set("block_length", length).set("channels", 1u); // Set one complete block.
330 fft1d transform(windowed, spectrum, length); // Plan an FFT of the weighted values.
331 shdr::pipeline pipeline( // Start a pipeline with the window dispatch.
332 shdr::dispatch(*window, shdr::extent(1))); // Begin the chain with input weighting.
333 transform.append_to( // Append the FFT plan to the same linked pipeline.
334 pipeline); // Pass the windowed resource directly to the FFT.
335 for (unsigned mode = 0; mode < 4; // Select each supported amplitude window.
336 ++mode) { // Reuse the linked resources for every supported window.
337 window->set("mode", mode); // Select the next amplitude response.
338 run.run(pipeline); // Execute the window and transform on the same backend.
339 std::vector<shdr::glsl::vec2> // Allocate independent window-weighted reference values.
340 expected( // Allocate independent window-weighted reference values.
341 length); // Retain an independently weighted reference block.
342 for (unsigned i = 0; i < length; ++i) { // Compute the reference gain at each position.
343 const double phase = // Calculate phase from the sample’s position within the window.
344 6.283185307179586476925286766559 * // Calculate phase from the sample’s position within the window.
345 i / // Calculate phase from the sample’s position within the window.
346 (length - 1u); // Sweep the symmetric window endpoints.
347 double gain = 1.0; // Rectangular weighting preserves input amplitudes.
348 if (mode == 1) // Select the Hann formula.
349 gain = 0.5 - 0.5 * std::cos(phase); // Evaluate a Hann window.
350 if (mode == 2) // Select the Hamming formula.
351 gain = 0.54 - 0.46 * std::cos(phase); // Evaluate a Hamming window.
352 if (mode == 3) // Select the Blackman formula.
353 gain = 0.42 - 0.5 * std::cos(phase) + // Combine the first two Blackman terms.
354 0.08 * std::cos(2.0 * phase); // Evaluate a Blackman window.
355 expected[i] = // Scale both complex components by the analytic window gain.
356 static_cast<const shdr::buffer<shdr::glsl::vec2>&>(input) // Scale both complex components by the analytic window gain.
357 .data()[i] * // Scale both complex components by the analytic window gain.
358 float(gain); // Apply the reference to both complex components.
359 } // The windowed signal intentionally differs from the original input.
360 fft_checks_detail::compare( // Compare window output with independent weighted samples.
361 windowed, // Validate the weighting stage itself.
362 expected, // Use the independently calculated expected values.
363 "signal window"); // Validate the weighting stage itself.
364 fft_checks_detail::compare( // Compare the linked FFT output after windowing.
365 spectrum, // Compare the linked FFT output after windowing.
366 fft_checks_detail::direct_dft( // Transform the independent window reference by direct DFT.
367 expected, length), // Transform the independent window reference by direct DFT.
368 "windowed FFT spectrum"); // Validate the complete linked chain.
369 } // Every window can be selected dynamically after pipeline construction.
370 std::cout // Write the window validation result to standard output.
371 << "Signal windows: rectangular, Hann, Hamming, and Blackman stages verified\n"; // Report the spectral-analysis checks.
372} // Applications can bind these windowed spectra to further GPU consumers.
373} // namespace shdr_examples
374#endif // Finish the FFT validation include guard.