RADLib
RADical C++ application framework
RADDsp.h
Go to the documentation of this file.
1 #ifndef RADDSP_H
2 #define RADDSP_H
3 
9 #pragma once
10 
11 #include <RADCore/RADLicenseGate.h>
12 
13 #include <array>
14 #include <cstddef>
15 #include <cstdint>
16 #include <initializer_list>
17 #include <vector>
18 
19 namespace RADDsp {
20 
22  struct ComplexSample {
24  float i = 0.0f;
26  float q = 0.0f;
27  };
28 
30  class FftPlan {
31  public:
33  FftPlan() = default;
35  explicit FftPlan(size_t size);
36 
38  bool reset(size_t size);
40  size_t size() const;
42  bool valid() const;
43 
45  bool transformInPlace(std::vector<ComplexSample>& samples, bool inverse = false) const;
47  bool transform(const std::vector<ComplexSample>& input, std::vector<ComplexSample>& output,
48  bool inverse = false) const;
50  std::vector<ComplexSample> transform(const std::vector<ComplexSample>& input,
51  bool inverse = false) const;
53  bool transformManyInPlace(std::vector<std::vector<ComplexSample>>& frames,
54  bool inverse = false) const;
56  bool transformMany(const std::vector<std::vector<ComplexSample>>& inputs,
57  std::vector<std::vector<ComplexSample>>& outputs, bool inverse = false) const;
59  std::vector<std::vector<ComplexSample>> transformMany(
60  const std::vector<std::vector<ComplexSample>>& inputs, bool inverse = false) const;
62  bool transformFramesInPlace(std::vector<ComplexSample>& frames, size_t frameCount,
63  bool inverse = false) const;
65  bool transformFrames(const std::vector<ComplexSample>& input, size_t frameCount,
66  std::vector<ComplexSample>& output, bool inverse = false) const;
68  std::vector<ComplexSample> transformFrames(const std::vector<ComplexSample>& input,
69  size_t frameCount, bool inverse = false) const;
70 
71  private:
72  size_t size_ = 0;
73  std::vector<size_t> bitReverse_;
74  std::vector<ComplexSample> twiddlesForward_;
75  std::vector<ComplexSample> twiddlesInverse_;
76  };
77 
79  class MatrixF {
80  public:
82  MatrixF() = default;
84  MatrixF(size_t rows, size_t columns, float value = 0.0f);
86  MatrixF(size_t rows, size_t columns, const std::vector<float>& rowMajor);
87 
89  size_t rows() const;
91  size_t columns() const;
93  bool empty() const;
95  const std::vector<float>& data() const;
97  std::vector<float>& data();
99  void resize(size_t rows, size_t columns, float value = 0.0f);
101  float& operator()(size_t row, size_t column);
103  float operator()(size_t row, size_t column) const;
105  float get(size_t row, size_t column) const;
107  bool set(size_t row, size_t column, float value);
108 
109  private:
110  size_t rows_ = 0;
111  size_t columns_ = 0;
112  std::vector<float> data_;
113  };
114 
118  float b0 = 1.0f;
120  float b1 = 0.0f;
122  float b2 = 0.0f;
124  float a1 = 0.0f;
126  float a2 = 0.0f;
127  };
128 
130  struct BiquadState {
132  float x1 = 0.0f;
134  float x2 = 0.0f;
136  float y1 = 0.0f;
138  float y2 = 0.0f;
139 
141  void reset();
143  float process(float input, const BiquadCoefficients& coeffs);
144  };
145 
147  struct EqualizerBand {
149  float frequencyHz = 1000.0f;
151  float gainDb = 0.0f;
153  float q = 1.0f;
154  };
155 
158  public:
160  void setSections(std::vector<BiquadCoefficients> sections);
164  void clear();
166  void reset();
168  const std::vector<BiquadCoefficients>& sections() const;
170  float processSample(float input);
172  std::vector<float> process(const std::vector<float>& samples);
173 
174  private:
175  std::vector<BiquadCoefficients> sections_;
176  std::vector<BiquadState> states_;
177  };
178 
181  public:
183  explicit Equalizer16Band(uint32_t sampleRate = 48000);
185  void setSampleRate(uint32_t sampleRate);
187  uint32_t sampleRate() const;
189  const std::array<EqualizerBand, 16>& bands() const;
191  bool setBand(size_t index, float frequencyHz, float gainDb, float q = 1.0f);
193  bool setBandGain(size_t index, float gainDb);
195  float bandGain(size_t index) const;
197  void reset();
199  std::vector<float> process(const std::vector<float>& samples);
200 
201  private:
202  void rebuild();
203 
204  uint32_t sampleRate_ = 48000;
205  std::array<EqualizerBand, 16> bands_{};
206  BiquadCascade cascade_;
207  };
208 
210  void applyGain(std::vector<float>& samples, float gain);
212  std::vector<float> mix(const std::vector<float>& a, const std::vector<float>& b);
214  float peak(const std::vector<float>& samples);
216  float rms(const std::vector<float>& samples);
218  void normalizePeak(std::vector<float>& samples, float targetPeak = 1.0f);
219 
221  std::vector<float> sine(float frequencyHz, float seconds, uint32_t sampleRate, float amplitude = 1.0f);
223  std::vector<float> whiteNoise(size_t count, float amplitude = 1.0f, uint32_t seed = 0x524144u);
224 
226  std::vector<float> linearResample(const std::vector<float>& interleaved, uint16_t channels,
227  uint32_t sourceRate, uint32_t targetRate);
229  std::vector<std::vector<float>> deinterleave(const std::vector<float>& interleaved, uint16_t channels);
231  std::vector<float> interleave(const std::vector<std::vector<float>>& channels);
233  std::vector<float> magnitude(const std::vector<ComplexSample>& samples);
235  std::vector<float> phase(const std::vector<ComplexSample>& samples);
237  void applyGain(std::vector<ComplexSample>& samples, float gain);
239  std::vector<ComplexSample> frequencyShift(const std::vector<ComplexSample>& samples,
240  float frequencyHz, uint32_t sampleRate);
241 
243  std::vector<float> hannWindow(size_t count);
245  std::vector<float> hammingWindow(size_t count);
247  std::vector<float> blackmanWindow(size_t count);
249  std::vector<float> rectangularWindow(size_t count);
251  void applyWindow(std::vector<float>& samples, const std::vector<float>& window);
253  std::vector<float> firFilter(const std::vector<float>& samples, const std::vector<float>& taps);
255  std::vector<float> convolve(const std::vector<float>& a, const std::vector<float>& b);
257  std::vector<float> crossCorrelate(const std::vector<float>& a, const std::vector<float>& b);
259  std::vector<float> normalizedCrossCorrelate(const std::vector<float>& a, const std::vector<float>& b);
261  std::vector<float> onePoleLowPass(const std::vector<float>& samples, float alpha);
263  FftPlan createFftPlan(size_t size);
265  std::vector<ComplexSample> fft(const std::vector<ComplexSample>& samples, bool inverse = false);
267  std::vector<std::vector<ComplexSample>> fftMany(
268  const std::vector<std::vector<ComplexSample>>& frames, bool inverse = false);
270  std::vector<ComplexSample> fftFrames(const std::vector<ComplexSample>& frames,
271  size_t frameSize, bool inverse = false);
273  std::vector<float> magnitudeSpectrum(const std::vector<ComplexSample>& spectrum);
275  std::vector<float> powerSpectrum(const std::vector<ComplexSample>& spectrum);
277  std::vector<float> designLowPassFIR(size_t taps, float cutoffHz, uint32_t sampleRate);
279  BiquadCoefficients designBiquadLowPass(float cutoffHz, uint32_t sampleRate, float q = 0.70710678f);
281  BiquadCoefficients designBiquadHighPass(float cutoffHz, uint32_t sampleRate, float q = 0.70710678f);
283  BiquadCoefficients designBiquadBandPass(float centerHz, uint32_t sampleRate, float q = 1.0f);
285  BiquadCoefficients designBiquadNotch(float centerHz, uint32_t sampleRate, float q = 1.0f);
287  BiquadCoefficients designBiquadPeakingEQ(float centerHz, uint32_t sampleRate, float gainDb, float q = 1.0f);
289  std::vector<float> decimate(const std::vector<float>& samples, size_t factor);
291  std::vector<float> interpolateZeros(const std::vector<float>& samples, size_t factor);
293  std::vector<float> mixCosine(const std::vector<float>& samples, float frequencyHz, uint32_t sampleRate);
295  void automaticGainControl(std::vector<float>& samples, float targetRms = 0.25f, float maxGain = 32.0f);
296 
298  MatrixF identityMatrix(size_t size);
300  MatrixF diagonalMatrix(const std::vector<float>& diagonal);
302  MatrixF transpose(const MatrixF& matrix);
304  MatrixF add(const MatrixF& a, const MatrixF& b);
306  MatrixF subtract(const MatrixF& a, const MatrixF& b);
308  MatrixF scale(const MatrixF& matrix, float scalar);
310  MatrixF hadamard(const MatrixF& a, const MatrixF& b);
312  MatrixF multiply(const MatrixF& a, const MatrixF& b);
314  std::vector<float> multiply(const MatrixF& matrix, const std::vector<float>& vector);
316  std::vector<float> multiply(const MatrixF& matrix, std::initializer_list<float> vector);
318  float dot(const std::vector<float>& a, const std::vector<float>& b);
320  float l2Norm(const std::vector<float>& vector);
322  std::vector<float> normalizeVector(const std::vector<float>& vector);
324  MatrixF outerProduct(const std::vector<float>& a, const std::vector<float>& b);
326  MatrixF covarianceMatrix(const MatrixF& observations);
328  bool invert(const MatrixF& matrix, MatrixF& inverse, float epsilon = 1.0e-6f);
330  bool solveLinearSystem(const MatrixF& matrix, const std::vector<float>& rhs,
331  std::vector<float>& solution, float epsilon = 1.0e-6f);
332 
333 } // namespace RADDsp
334 
335 #endif
Stateful cascade of biquad sections.
Definition: RADDsp.h:157
float processSample(float input)
Processes one sample through all sections.
void reset()
Resets all section states.
std::vector< float > process(const std::vector< float > &samples)
Processes a block through all sections and updates states.
void addSection(BiquadCoefficients section)
Appends one section and creates matching state.
void clear()
Removes all sections and states.
const std::vector< BiquadCoefficients > & sections() const
Returns current section coefficients.
void setSections(std::vector< BiquadCoefficients > sections)
Replaces all sections and resets section states.
16-band peaking-EQ built from a biquad cascade.
Definition: RADDsp.h:180
float bandGain(size_t index) const
Returns one band gain or 0 for invalid index.
void reset()
Resets all filter states.
Equalizer16Band(uint32_t sampleRate=48000)
Creates an equalizer for sampleRate with default band centers.
bool setBand(size_t index, float frequencyHz, float gainDb, float q=1.0f)
Sets one band center, gain, and Q.
bool setBandGain(size_t index, float gainDb)
Sets one band gain while preserving center/Q.
void setSampleRate(uint32_t sampleRate)
Sets sample rate and rebuilds coefficients.
const std::array< EqualizerBand, 16 > & bands() const
Returns all band settings.
uint32_t sampleRate() const
Returns configured sample rate.
std::vector< float > process(const std::vector< float > &samples)
Processes a block through the equalizer.
Reusable power-of-two complex FFT plan with cached bit-reversal and twiddle factors.
Definition: RADDsp.h:30
std::vector< ComplexSample > transform(const std::vector< ComplexSample > &input, bool inverse=false) const
Transforms input and returns an empty vector if the plan or sample count is invalid.
std::vector< ComplexSample > transformFrames(const std::vector< ComplexSample > &input, size_t frameCount, bool inverse=false) const
Transforms frameCount contiguous frames, returning an empty vector on invalid input.
bool transformManyInPlace(std::vector< std::vector< ComplexSample >> &frames, bool inverse=false) const
Transforms many equal-sized frames in place. Returns false if any frame has the wrong size.
bool valid() const
Returns true when the plan has a non-zero power-of-two size.
FftPlan()=default
Creates an empty invalid plan.
bool transformInPlace(std::vector< ComplexSample > &samples, bool inverse=false) const
Transforms samples in place. Returns false if the plan or sample count is invalid.
bool transformFramesInPlace(std::vector< ComplexSample > &frames, size_t frameCount, bool inverse=false) const
Transforms frameCount contiguous frames of size() in place.
bool transform(const std::vector< ComplexSample > &input, std::vector< ComplexSample > &output, bool inverse=false) const
Transforms input into output. Returns false if the plan or sample count is invalid.
bool reset(size_t size)
Rebuilds the plan and returns false for zero or non-power-of-two sizes.
bool transformMany(const std::vector< std::vector< ComplexSample >> &inputs, std::vector< std::vector< ComplexSample >> &outputs, bool inverse=false) const
Transforms many equal-sized frames. Returns false if any frame has the wrong size.
FftPlan(size_t size)
Creates a plan for size when size is a non-zero power of two.
size_t size() const
Returns the FFT length owned by this plan.
bool transformFrames(const std::vector< ComplexSample > &input, size_t frameCount, std::vector< ComplexSample > &output, bool inverse=false) const
Transforms frameCount contiguous frames of size(). Returns false on invalid sizes.
std::vector< std::vector< ComplexSample > > transformMany(const std::vector< std::vector< ComplexSample >> &inputs, bool inverse=false) const
Transforms many equal-sized frames, returning an empty vector on invalid input.
Row-major dense floating-point matrix for DSP transforms and small linear algebra.
Definition: RADDsp.h:79
std::vector< float > & data()
Returns the mutable matrix backing storage in row-major order.
MatrixF(size_t rows, size_t columns, const std::vector< float > &rowMajor)
Creates a rows x columns matrix from row-major values; missing values are zero-filled.
MatrixF()=default
Creates an empty matrix.
float & operator()(size_t row, size_t column)
Fast unchecked element access for hot DSP loops.
size_t rows() const
Returns the number of rows.
bool empty() const
Returns true when either dimension is zero.
const std::vector< float > & data() const
Returns the matrix backing storage in row-major order.
float operator()(size_t row, size_t column) const
Fast unchecked element access for hot DSP loops.
size_t columns() const
Returns the number of columns.
bool set(size_t row, size_t column, float value)
Safely writes one value and returns false for invalid indices.
void resize(size_t rows, size_t columns, float value=0.0f)
Resizes the matrix and fills all values with value.
float get(size_t row, size_t column) const
Safely reads one value, returning 0 for invalid indices.
MatrixF(size_t rows, size_t columns, float value=0.0f)
Creates a rows x columns matrix initialized to value.
Definition: RADDsp.h:19
float peak(const std::vector< float > &samples)
Returns the absolute peak amplitude in a scalar sample buffer.
std::vector< float > mixCosine(const std::vector< float > &samples, float frequencyHz, uint32_t sampleRate)
Mixes scalar samples with a cosine oscillator.
std::vector< float > interpolateZeros(const std::vector< float > &samples, size_t factor)
Inserts factor-1 zeros between samples.
MatrixF outerProduct(const std::vector< float > &a, const std::vector< float > &b)
Builds the outer product a * b^T.
BiquadCoefficients designBiquadPeakingEQ(float centerHz, uint32_t sampleRate, float gainDb, float q=1.0f)
Designs an RBJ peaking equalizer biquad.
MatrixF hadamard(const MatrixF &a, const MatrixF &b)
Multiplies matching matrix elements, or returns an empty matrix.
std::vector< float > designLowPassFIR(size_t taps, float cutoffHz, uint32_t sampleRate)
Designs a windowed-sinc low-pass FIR filter.
std::vector< float > crossCorrelate(const std::vector< float > &a, const std::vector< float > &b)
Cross-correlates a and b using FFT-based convolution with reversed b.
std::vector< float > whiteNoise(size_t count, float amplitude=1.0f, uint32_t seed=0x524144u)
Generates deterministic white noise using seed for repeatable tests and examples.
void applyGain(std::vector< float > &samples, float gain)
Multiplies every scalar sample by gain in place.
std::vector< std::vector< float > > deinterleave(const std::vector< float > &interleaved, uint16_t channels)
Splits interleaved channel data into one vector per channel.
std::vector< float > rectangularWindow(size_t count)
Builds a rectangular window with count coefficients.
BiquadCoefficients designBiquadNotch(float centerHz, uint32_t sampleRate, float q=1.0f)
Designs an RBJ notch biquad.
float rms(const std::vector< float > &samples)
Returns root-mean-square amplitude for a scalar sample buffer.
MatrixF diagonalMatrix(const std::vector< float > &diagonal)
Builds a diagonal matrix from diagonal values.
bool invert(const MatrixF &matrix, MatrixF &inverse, float epsilon=1.0e-6f)
Computes the inverse of a square matrix with partial pivoting.
std::vector< float > convolve(const std::vector< float > &a, const std::vector< float > &b)
Convolves two scalar buffers using an FFT-based implementation.
void normalizePeak(std::vector< float > &samples, float targetPeak=1.0f)
Scales samples so their absolute peak equals targetPeak when possible.
std::vector< float > firFilter(const std::vector< float > &samples, const std::vector< float > &taps)
Applies an FIR filter using taps and returns filtered samples.
std::vector< float > phase(const std::vector< ComplexSample > &samples)
Converts complex I/Q samples to phase values in radians.
std::vector< float > powerSpectrum(const std::vector< ComplexSample > &spectrum)
Converts complex spectrum bins to power values.
MatrixF subtract(const MatrixF &a, const MatrixF &b)
Subtracts matrices with matching dimensions, or returns an empty matrix.
MatrixF scale(const MatrixF &matrix, float scalar)
Multiplies every matrix value by scalar.
FftPlan createFftPlan(size_t size)
Creates a reusable FFT plan for power-of-two sizes.
BiquadCoefficients designBiquadBandPass(float centerHz, uint32_t sampleRate, float q=1.0f)
Designs an RBJ band-pass biquad.
std::vector< float > blackmanWindow(size_t count)
Builds a Blackman window with count coefficients.
std::vector< float > magnitude(const std::vector< ComplexSample > &samples)
Converts complex I/Q samples to magnitude values.
void automaticGainControl(std::vector< float > &samples, float targetRms=0.25f, float maxGain=32.0f)
Applies simple automatic gain control toward targetRms.
void applyWindow(std::vector< float > &samples, const std::vector< float > &window)
Applies a window to samples in place up to the smaller buffer size.
std::vector< float > decimate(const std::vector< float > &samples, size_t factor)
Keeps every factor-th sample after low-pass filtering should be applied by caller when needed.
BiquadCoefficients designBiquadHighPass(float cutoffHz, uint32_t sampleRate, float q=0.70710678f)
Designs an RBJ high-pass biquad.
std::vector< float > mix(const std::vector< float > &a, const std::vector< float > &b)
Adds two scalar buffers sample-by-sample; extra tail samples are preserved from the longer input.
MatrixF covarianceMatrix(const MatrixF &observations)
Computes a covariance matrix from row-major observations where each row is one observation.
float dot(const std::vector< float > &a, const std::vector< float > &b)
Returns the dot product of equal-size vectors, or 0 for mismatched sizes.
std::vector< float > normalizedCrossCorrelate(const std::vector< float > &a, const std::vector< float > &b)
Cross-correlates and normalizes by signal energy when possible.
float l2Norm(const std::vector< float > &vector)
Returns the Euclidean norm of vector.
std::vector< float > onePoleLowPass(const std::vector< float > &samples, float alpha)
Applies a simple one-pole low-pass filter with alpha in the 0..1 smoothing range.
std::vector< float > hammingWindow(size_t count)
Builds a Hamming window with count coefficients.
bool solveLinearSystem(const MatrixF &matrix, const std::vector< float > &rhs, std::vector< float > &solution, float epsilon=1.0e-6f)
Solves matrix * x = rhs using Gaussian elimination with partial pivoting.
BiquadCoefficients designBiquadLowPass(float cutoffHz, uint32_t sampleRate, float q=0.70710678f)
Designs an RBJ low-pass biquad.
MatrixF transpose(const MatrixF &matrix)
Returns the transpose of matrix.
std::vector< ComplexSample > fftFrames(const std::vector< ComplexSample > &frames, size_t frameSize, bool inverse=false)
Computes many contiguous equal-sized power-of-two FFT frames using one reusable plan.
std::vector< ComplexSample > fft(const std::vector< ComplexSample > &samples, bool inverse=false)
Computes a complex FFT for power-of-two sizes and DFT fallback otherwise.
std::vector< float > normalizeVector(const std::vector< float > &vector)
Returns a normalized copy of vector, or zeros for zero-length input energy.
std::vector< float > linearResample(const std::vector< float > &interleaved, uint16_t channels, uint32_t sourceRate, uint32_t targetRate)
Resamples interleaved scalar audio using linear interpolation.
std::vector< std::vector< ComplexSample > > fftMany(const std::vector< std::vector< ComplexSample >> &frames, bool inverse=false)
Computes many equal-sized power-of-two FFTs using one reusable plan.
std::vector< ComplexSample > frequencyShift(const std::vector< ComplexSample > &samples, float frequencyHz, uint32_t sampleRate)
Frequency-shifts complex I/Q samples by frequencyHz at sampleRate.
std::vector< float > sine(float frequencyHz, float seconds, uint32_t sampleRate, float amplitude=1.0f)
Generates a sine wave at frequencyHz for seconds at sampleRate.
std::vector< float > magnitudeSpectrum(const std::vector< ComplexSample > &spectrum)
Converts complex spectrum bins to magnitudes.
std::vector< float > hannWindow(size_t count)
Builds a Hann window with count coefficients.
std::vector< float > interleave(const std::vector< std::vector< float >> &channels)
Combines per-channel buffers into interleaved channel order.
MatrixF identityMatrix(size_t size)
Builds a square identity matrix.
MatrixF multiply(const MatrixF &a, const MatrixF &b)
Multiplies two dense matrices, or returns an empty matrix for incompatible dimensions.
MatrixF add(const MatrixF &a, const MatrixF &b)
Adds matrices with matching dimensions, or returns an empty matrix.
Normalized direct-form biquad coefficients with a0 folded into all terms.
Definition: RADDsp.h:116
float b1
Feed-forward coefficient for x[n-1].
Definition: RADDsp.h:120
float a2
Feedback coefficient for y[n-2], assuming y = b*x - a*y.
Definition: RADDsp.h:126
float a1
Feedback coefficient for y[n-1], assuming y = b*x - a*y.
Definition: RADDsp.h:124
float b2
Feed-forward coefficient for x[n-2].
Definition: RADDsp.h:122
float b0
Feed-forward coefficient for x[n].
Definition: RADDsp.h:118
Stateful direct-form biquad delay elements.
Definition: RADDsp.h:130
float x2
Input sample before x1.
Definition: RADDsp.h:134
float process(float input, const BiquadCoefficients &coeffs)
Processes one sample through coeffs and updates state.
float x1
Previous input sample.
Definition: RADDsp.h:132
void reset()
Clears all delay elements.
float y2
Output sample before y1.
Definition: RADDsp.h:138
float y1
Previous output sample.
Definition: RADDsp.h:136
Complex I/Q sample used by SDR and constellation-oriented helpers.
Definition: RADDsp.h:22
float q
Quadrature component.
Definition: RADDsp.h:26
float i
In-phase component.
Definition: RADDsp.h:24
One 16-band equalizer band setting.
Definition: RADDsp.h:147
float frequencyHz
Center frequency in hertz.
Definition: RADDsp.h:149
float q
Peaking-EQ Q factor.
Definition: RADDsp.h:153
float gainDb
Gain in decibels.
Definition: RADDsp.h:151