Source browser
examples/fft_support.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_support.hpp
2 * @brief Batched complex FFT and block signal-processing examples.
3 * A vec2 stores real and imaginary parts. Channels occupy consecutive blocks.
4 * Dispatch boundaries order the dependent FFT stages.
5 */
6#ifndef SHDR_EXAMPLES_FFT_SUPPORT_HPP // Prevent duplicate helper definitions.
7#define SHDR_EXAMPLES_FFT_SUPPORT_HPP // Mark this helper as included.
8#include "shdr.h" // Use the single-header compute runtime.
9#include <cmath> // Compute shared phase tables once on the host.
10#include <memory> // Keep each stage's kernel alive with its pipeline.
11#include <vector> // Store the kernels for all butterfly stages.
12
13namespace shdr_examples { // Keep example helpers separate from the runtime API.
14namespace fft_detail { // Share dimension checks between signal helpers.
15inline unsigned power_of_two(unsigned length) { // Require a positive power-of-two transform length.
16 if (!length || (length & (length - 1u))) { // Test that exactly one length bit is set.
17 throw shdr::error( // Report the failed format, numerical or resource precondition.
18 "FFT length must be a positive power of two"); // Report the unsupported length before allocation.
19 } // A length of one uses the permutation copy directly.
20 return length; // Use the checked length in member initialization.
21} // Use the requested block length exactly.
22inline unsigned channel_count(unsigned channels) { // Require at least one independent signal.
23 if (!channels) { // Require at least one signal channel before creating dispatches.
24 throw shdr::error( // Report the failed format, numerical or resource precondition.
25 "signal channel count must be positive"); // Explain the invalid dimension.
26 } // Positive channel counts are checked for overflow with their length next.
27 return channels; // Preserve the requested channel count.
28} // Every helper uses the same channel-major layout.
29inline unsigned elements(unsigned length, // Check shader index arithmetic before allocation.
30 unsigned channels) { // Check shader index arithmetic before allocation.
31 if (length > 0xffffffffu / channels) { // Validate the product capacity within shader unsigned indexing.
32 throw shdr::error( // Report the failed format, numerical or resource precondition.
33 "signal dimensions exceed GLSL uint indexing"); // Keep host and shader indexing consistent.
34 } // The product now fits the shader's unsigned type.
35 return length * channels; // Return the total number of complex values.
36} // Buffer constructors also check allocation byte counts.
37inline void // Inspect the input's allocated size.
38buffers(const shdr::buffer<shdr::glsl::vec2>& input, // Inspect the input's allocated size.
39 const shdr::buffer<shdr::glsl::vec2>& output, // Read the output capacity from host resource metadata.
40 unsigned count) { // Read the output capacity from host resource metadata.
41 if (input.size() < count || // Require storage for every requested channel.
42 output.size() < count) { // Require storage for every requested channel.
43 throw shdr::error( // Report the failed format, numerical or resource precondition.
44 "signal buffers do not contain every requested channel"); // Fail before a dispatch could access missing values.
45 } // Extra capacity remains available to other application code.
46} // Dimension checks use the resource capacity metadata.
47inline void append(shdr::pipeline& destination, // Link existing stages in their execution order.
48 const shdr::pipeline& source) { // Link existing stages in their execution order.
49 for (const shdr::stage& stage : // Retain the references already stored by each helper.
50 source.stages()) { // Retain the references already stored by each helper.
51 destination += stage; // Place this dispatch after all previously appended work.
52 } // The helper objects must remain alive while their stages execute.
53} // This avoids copying any input or intermediate values.
54} // namespace fft_detail
55
56/** @brief An out-of-place radix-two complex FFT with O(N log N) work.
57 * Forward transforms are unscaled. Inverse transforms include 1/N.
58 * Input and output may share storage because working buffers separate stages.
59 * Keep this object alive while using its pipeline or an appended dispatch.
60 */
61class fft1d { // Own one fixed-size transform plan and its working resources.
62 unsigned length_, channels_, // Keep dimensions checked and consistent across all stages.
63 count_; // Keep dimensions checked and consistent across all stages.
64 shdr::buffer<shdr::glsl::vec2> first_, second_, // Alternate working buffers and share one phase table.
65 twiddles_; // Alternate working buffers and share one phase table.
66 std::vector<std::unique_ptr<shdr::kernel>> // Retain each permutation and butterfly kernel.
67 kernels_; // Retain each permutation and butterfly kernel.
68 shdr::pipeline pipeline_; // Preserve the dependency order established at construction.
69 public: // Expose execution while keeping working resources private.
70 /** @brief Plan a batched forward or inverse transform over channel-major buffers. */
71 fft1d(const shdr::buffer<shdr::glsl::vec2>& input, // Bind the original complex values.
72 const shdr::buffer<shdr::glsl::vec2>& output, // Bind the final transform result.
73 unsigned length, // Select block dimensions and transform direction.
74 unsigned channels = 1, // Select block dimensions and transform direction.
75 bool inverse = false) // Select block dimensions and transform direction.
76 : length_(fft_detail::power_of_two(length)), // Check the radix-two requirement first.
77 channels_(fft_detail::channel_count( // Check the channel count before computing products.
78 channels)), // Check the channel count before computing products.
79 count_( // Keep every shader index representable.
80 fft_detail::elements(length_, channels_)), // Keep every shader index representable.
81 first_(count_), second_(count_), // Allocate working storage once.
82 twiddles_(std::max(1u, length_ / 2u)) { // Allocate working storage once.
83 fft_detail::buffers(input, output, count_); // Validate caller-owned resource capacities.
84 const double tau = // Use one full turn for complex phases.
85 6.283185307179586476925286766559; // Use one full turn for complex phases.
86 shdr::glsl::vec2* phases = twiddles_.data(); // Fill the host table before GPU preparation.
87 for (unsigned k = 0; k < twiddles_.size(); // Compute the phases needed by every butterfly stage.
88 ++k) { // Compute the phases needed by every butterfly stage.
89 const double angle = -tau * double(k) / length_; // Store forward-transform phases.
90 phases[k] = shdr::glsl::vec2( // Store matching float values for both backends.
91 float(std::cos(angle)), // Store matching float values for both backends.
92 float(std::sin(angle))); // Store matching float values for both backends.
93 } // Butterflies read these precomputed phase values.
94 unsigned bits = 0; // Count the binary digits in a block position.
95 for (unsigned size = length_; size > 1u; size /= 2u) { // Repeated halving measures log2(N).
96 ++bits; // Add one bit for this power of two.
97 } // Length one has zero permutation bits.
98 std::unique_ptr<shdr::kernel> permutation = // Create the index-ordering stage.
99 shdr::make_kernel("fft_bit_reverse_kernel"); // Create the index-ordering stage.
100 permutation->bind("input_values", input); // Read the application input.
101 permutation->bind("output_values", // Length one needs only a copy.
102 length_ == 1u ? output : first_); // Length one needs only a copy.
103 permutation->set("block_length", length_) // Supply the per-channel block length.
104 .set("channels", channels_) // Supply the channels processing parameter.
105 .set("bits", bits); // Supply the permutation dimensions.
106 pipeline_ += shdr::dispatch( // Dispatch this configured shader with the following work group geometry.
107 *permutation,
108 shdr::extent( // Use quotient and remainder to cover every value.
109 count_ / 64u + // Use quotient and remainder to cover every value.
110 (count_ % 64u != 0u))); // Use quotient and remainder to cover every value.
111 kernels_.push_back( // Keep the stage alive after this constructor returns.
112 std::move(permutation)); // Keep the stage alive after this constructor returns.
113 bool read_first = true; // The permutation initially fills the first working buffer.
114 for (unsigned span = 2u; // Grow butterfly groups from two values to the full block.
115 span <= length_ && // Grow butterfly groups from two values to the full block.
116 length_ > 1u;) { // Grow butterfly groups from two values to the full block.
117 std::unique_ptr<shdr::kernel> butterfly = // Give each stage independent uniforms.
118 shdr::make_kernel("fft_butterfly_kernel"); // Give each stage independent uniforms.
119 butterfly->bind("input_values", // Read the preceding stage's output.
120 read_first ? first_ : second_); // Read the preceding stage's output.
121 butterfly->bind( // Send the last stage to caller storage.
122 "output_values", // Use the exact diagnostic or comparison label for this case.
123 span == length_ // Send the last stage to caller storage.
124 ? output // Send the last stage to caller storage.
125 : (read_first ? second_ : first_)); // Send the last stage to caller storage.
126 butterfly->bind("twiddles", twiddles_); // Reuse the same phase table at every span.
127 butterfly->set("block_length", length_) // Supply the per-channel block length.
128 .set("channels", channels_) // Supply the channels processing parameter.
129 .set("span", span); // Set this stage's butterfly width.
130 butterfly->set("direction", // Conjugate the phase table for an inverse transform.
131 inverse ? -1.0f // Conjugate the phase table for an inverse transform.
132 : 1.0f); // Conjugate the phase table for an inverse transform.
133 butterfly->set("normalization", // Scale exactly once for the inverse.
134 inverse && span == length_ // Scale exactly once for the inverse.
135 ? 1.0f / length_ // Scale exactly once for the inverse.
136 : 1.0f); // Scale exactly once for the inverse.
137 pipeline_ += shdr::dispatch( // Dispatch this configured shader with the following work group geometry.
138 *butterfly,
139 shdr::extent( // Use quotient and remainder to cover every butterfly output.
140 count_ / 64u + // Use quotient and remainder to cover every butterfly output.
141 (count_ % 64u != // Use quotient and remainder to cover every butterfly output.
142 0u))); // Use quotient and remainder to cover every butterfly output.
143 kernels_.push_back( // Retain the stage's source, bindings, and uniforms.
144 std::move(butterfly)); // Retain the stage's source, bindings, and uniforms.
145 read_first = !read_first; // Exchange the working-buffer roles for the next stage.
146 if (span == length_) // Stop before the final span could overflow on doubling.
147 break; // Stop before the final span could overflow on doubling.
148 span *= 2u; // Double the butterfly group size.
149 } // The last dispatch has written the complete transformed block.
150 } // Preparation and execution are separate from plan construction.
151 /** @brief Return the fixed transform stages for CPU, GPU, or asynchronous execution. */
152 const shdr::pipeline& // Return the existing plan and its bound resources.
153 processing_pipeline() const { // Return the existing plan and its bound resources.
154 return pipeline_; // Run this pipeline on the CPU or a GPU device.
155 } // GPU preparation can compile and upload the plan before a hot loop.
156 /** @brief Append all FFT stages to a larger ordered pipeline. */
157 void append_to( // Keep original kernel ownership in this FFT plan.
158 shdr::pipeline& destination) const { // Keep original kernel ownership in this FFT plan.
159 fft_detail::append(destination, // Link the transform after earlier application work.
160 pipeline_); // Link the transform after earlier application work.
161 } // Values remain in their bound CPU or GPU buffers.
162}; // End the reusable complex-transform helper.
163
164/** @brief Frequency responses available for independent FFT signal blocks. */
165enum class frequency_filter { // Match the native mask shader's mode values.
166 low_pass, // Match the native mask shader's mode values.
167 high_pass, // Match the native mask shader's mode values.
168 band_pass, // Match the native mask shader's mode values.
169 band_stop // Match the native mask shader's mode values.
170}; // Match the native mask shader's mode values.
171
172/** @brief Forward FFT followed by a symmetric spectral mask and an inverse FFT.
173 * This is circular block filtering. Use padding and overlap for a linear stream.
174 * A positive transition width joins pass and stop bands with a cosine response.
175 */
176class fft_filter { // Own all resources needed to filter one fixed-size batch.
177 unsigned length_, channels_, count_; // Store checked FFT dimensions.
178 shdr::buffer<shdr::glsl::vec2> spectrum_, // Retain spectra for later inspection or linked processing.
179 masked_; // Retain spectra for later inspection or linked processing.
180 fft1d forward_, inverse_; // Reuse the same transform implementation in both directions.
181 std::unique_ptr<shdr::kernel> mask_; // Keep the frequency-selection dispatch alive.
182 shdr::pipeline pipeline_; // Link forward transform, mask, and inverse transform.
183 public: // Expose the complete processing chain and optional spectral results.
184 /** @brief Construct a frequency-selective filter using hertz for rate and cutoffs. */
185 fft_filter(const shdr::buffer<shdr::glsl::vec2>& // Bind one or more real or complex input blocks.
186 input, // Bind one or more real or complex input blocks.
187 const shdr::buffer<shdr::glsl::vec2>& output, // Bind reconstructed filtered blocks.
188 unsigned length, // Define block geometry and frequency units.
189 unsigned channels, // Define block geometry and frequency units.
190 float sample_rate, // Define block geometry and frequency units.
191 frequency_filter mode, // Select one cutoff or two band edges.
192 float lower_hz, // Select one cutoff or two band edges.
193 float upper_hz = 0.0f, // Select one cutoff or two band edges.
194 float transition_hz = 0.0f) // Select one cutoff or two band edges.
195 : length_( // Reject a block length unsupported by the FFT.
196 fft_detail::power_of_two(length)), // Reject a block length unsupported by the FFT.
197 channels_(fft_detail::channel_count(channels)), // Require at least one channel.
198 count_( // Check index products before allocation.
199 fft_detail::elements(length_, channels_)), // Check index products before allocation.
200 spectrum_(count_), masked_(count_), // Allocate reusable frequency-domain buffers.
201 forward_(input, spectrum_, length_, channels_), // Plan the forward transform.
202 inverse_( // Plan normalized inverse reconstruction.
203 masked_, output, length_, channels_, true), // Plan normalized inverse reconstruction.
204 mask_(shdr::make_kernel("signal_frequency_mask_kernel")) { // Create the response stage.
205 if (!std::isfinite(sample_rate) || // Frequency units need a valid positive rate.
206 sample_rate <= 0.0f) { // Frequency units need a valid positive rate.
207 throw shdr::error( // Report the failed format, numerical or resource precondition.
208 "FFT filter sample rate must be finite and positive"); // Reject invalid conversion from bins to hertz.
209 } // Nyquist is half this sample rate.
210 if (!std::isfinite(lower_hz) || lower_hz <= 0.0f || // Validate the first cutoff.
211 lower_hz >= sample_rate * 0.5f) { // Validate the first cutoff.
212 throw shdr::error( // Report the failed format, numerical or resource precondition.
213 "FFT filter cutoff must lie strictly between zero and Nyquist"); // Keep the response inside representable frequencies.
214 } // Low-pass and high-pass use this single cutoff.
215 if (mode >= frequency_filter::band_pass && // Validate both edges of a band.
216 (!std::isfinite(upper_hz) || upper_hz <= lower_hz || // Validate both edges of a band.
217 upper_hz >= sample_rate * 0.5f)) { // Validate both edges of a band.
218 throw shdr::error( // Report the failed format, numerical or resource precondition.
219 "FFT filter band edges must increase below Nyquist"); // Reject reversed or missing edges.
220 } // The band can now be evaluated symmetrically.
221 if (!std::isfinite(transition_hz) || // Require a finite nonnegative transition width.
222 transition_hz < 0.0f) { // Require a finite nonnegative transition width.
223 throw shdr::error( // Report the failed format, numerical or resource precondition.
224 "FFT filter transition must be finite and nonnegative"); // Zero selects a hard bin response.
225 } // Wider transitions deliberately soften the filter response.
226 const unsigned selected = // Match the host enumeration to the GLSL uniform.
227 static_cast<unsigned>(mode); // Match the host enumeration to the GLSL uniform.
228 if (selected > 3u) // Reject values outside the documented responses.
229 throw shdr::error( // Report the failed format, numerical or resource precondition.
230 "unknown FFT filter mode"); // Reject values outside the documented responses.
231 mask_->bind("input_values", spectrum_) // Connect the input_values shader port to spectrum_ storage.
232 .bind("output_values", masked_); // Link the transform result to the mask.
233 mask_->set("block_length", length_) // Supply the per-channel block length.
234 .set("channels", channels_) // Supply the channels processing parameter.
235 .set("mode", selected); // Supply channel and bin dimensions.
236 mask_->set("sample_rate", sample_rate) // Supply the sample_rate processing parameter.
237 .set("lower_hz", lower_hz) // Supply the lower_hz processing parameter.
238 .set("upper_hz", upper_hz); // Describe the response in hertz.
239 mask_->set("transition_hz", transition_hz); // Select hard or smooth cutoff behavior.
240 forward_.append_to(pipeline_); // First compute the unshifted spectrum.
241 pipeline_ += shdr::dispatch( // Dispatch this configured shader with the following work group geometry.
242 *mask_,
243 shdr::extent(count_ / 64u + // Dispatch the mask over every resident spectral value.
244 (count_ % 64u != // Dispatch the mask over every resident spectral value.
245 0u))); // Dispatch the mask over every resident spectral value.
246 inverse_.append_to(pipeline_); // Reconstruct the filtered signal on the same backend.
247 } // The plan can be prepared once and executed repeatedly.
248 /** @brief Return the full transform, mask, and reconstruction pipeline. */
249 const shdr::pipeline& // Expose all linked stages as one operation.
250 processing_pipeline() const { // Expose all linked stages as one operation.
251 return pipeline_; // CPU and GPU execution use the same plan.
252 } // The final transform writes to the chosen output, including aliased input storage.
253 /** @brief Append the complete block filter to another processing chain. */
254 void append_to( // Preserve dependencies when linking further stages.
255 shdr::pipeline& destination) const { // Preserve dependencies when linking further stages.
256 fft_detail::append(destination, pipeline_); // Add every dispatch in processing order.
257 } // The filter object owns its kernels and must remain alive.
258 /** @brief Inspect the original unshifted FFT spectrum when needed. */
259 const shdr::buffer<shdr::glsl::vec2>& // Prefer a const reference for optional downloads.
260 spectrum() const { // Prefer a const reference for optional downloads.
261 return spectrum_; // GPU data is read back only when its data pointer is requested.
262 } // Downstream GPU kernels can bind this resource directly.
263 /** @brief Inspect the spectrum after applying the frequency response. */
264 const shdr::buffer<shdr::glsl::vec2>& // Keep the masked values available as a shared resource.
265 filtered_spectrum() const { // Keep the masked values available as a shared resource.
266 return masked_; // Return the shared masked spectrum for another consumer.
267 } // Reading on the host remains an explicit synchronization point.
268}; // End the frequency-domain block-filter helper.
269} // namespace shdr_examples
270#endif // Finish the FFT helper include guard.