6502/CLK
1
0
Fork 0
mirror of https://github.com/TomHarte/CLK.git synced 2026-04-20 10:17:05 +00:00
Code Issues Projects Releases Wiki Activity

Improve style.

Browse Source
This commit is contained in:
Thomas Harte
2025-02-11 21:48:56 -05:00
parent 9470775292
commit 94b972aaf4
2 changed files with 66 additions and 57 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SignalProcessing/FIRFilter.cpp
+61 -47
View File
@@ -9,9 +9,11 @@
#include "FIRFilter.hpp"
#include <cmath>
#include <numeric>
#ifndef M_PI
static constexpr float M_PI = 3.1415926f;
// TODO: use std::numbers::pi_v when switching to C++20.
#endif
using namespace SignalProcessing;
@@ -36,8 +38,10 @@ using namespace SignalProcessing;
"DIGITAL SIGNAL PROCESSING, II", IEEE Press, pages 123-126.
*/
namespace {
/*! Evaluates the 0th order Bessel function at @c a. */
float FIRFilter::ino(float a) {
constexpr float ino(float a) {
float d = 0.0f;
float ds = 1.0f;
float s = 1.0f;
@@ -51,73 +55,78 @@ float FIRFilter::ino(float a) {
return s;
}
void FIRFilter::coefficients_for_idealised_filter_response(short *filter_coefficients, float *A, float attenuation, std::size_t number_of_taps) {
/* calculate alpha, which is the Kaiser-Bessel window shape factor */
float a; // to take the place of alpha in the normal derivation
if(attenuation < 21.0f) {
a = 0.0f;
} else {
if(attenuation > 50.0f)
a = 0.1102f * (attenuation - 8.7f);
else
a = 0.5842f * powf(attenuation - 21.0f, 0.4f) + 0.7886f * (attenuation - 21.0f);
}
std::vector<short> coefficients_for_idealised_filter_response(
const std::vector<float> &A,
const float attenuation,
const std::size_t number_of_taps
) {
/* Calculate alpha, the Kaiser-Bessel window shape factor */
const float a = [&] {
if(attenuation < 21.0f) {
return 0.0f;
} else if(attenuation > 50.0f) {
return 0.1102f * (attenuation - 8.7f);
} else {
return 0.5842f * powf(attenuation - 21.0f, 0.4f) + 0.7886f * (attenuation - 21.0f);
}
} ();
std::vector<float> filter_coefficients_float(number_of_taps);
/* work out the right hand side of the filter coefficients */
std::size_t Np = (number_of_taps - 1) / 2;
float I0 = ino(a);
float Np_squared = float(Np * Np);
for(unsigned int i = 0; i <= Np; ++i) {
/* Work out the right hand side of the filter coefficients. */
const float I0 = ino(a);
const std::size_t Np = (number_of_taps - 1) / 2;
const float Np_squared = float(Np * Np);
for(std::size_t i = 0; i <= Np; ++i) {
filter_coefficients_float[Np + i] =
A[i] *
ino(a * sqrtf(1.0f - (float(i * i) / Np_squared) )) /
I0;
}
/* coefficients are symmetrical, so copy from right hand side to left side */
/* Coefficients are symmetrical, so copy from right hand side to left. */
for(std::size_t i = 0; i < Np; ++i) {
filter_coefficients_float[i] = filter_coefficients_float[number_of_taps - 1 - i];
}
/* scale back up so that we retain 100% of input volume */
float coefficientTotal = 0.0f;
for(std::size_t i = 0; i < number_of_taps; ++i) {
coefficientTotal += filter_coefficients_float[i];
}
/* Scale back up to retain 100% of input volume. */
const float coefficientTotal =
std::accumulate(filter_coefficients_float.begin(), filter_coefficients_float.end(), 0.0f);
/* we'll also need integer versions, potentially */
float coefficientMultiplier = 1.0f / coefficientTotal;
/* Hence produce integer versions. */
const float coefficientMultiplier = 1.0f / coefficientTotal;
std::vector<short> filter_coefficients;
filter_coefficients.reserve(number_of_taps);
for(std::size_t i = 0; i < number_of_taps; ++i) {
filter_coefficients[i] = short(filter_coefficients_float[i] * FixedMultiplier * coefficientMultiplier);
filter_coefficients.push_back(short(filter_coefficients_float[i] * FixedMultiplier * coefficientMultiplier));
}
return filter_coefficients;
}
}
std::vector<float> FIRFilter::get_coefficients() const {
std::vector<float> coefficients;
coefficients.reserve(filter_coefficients_.size());
for(const auto short_coefficient: filter_coefficients_) {
coefficients.push_back(float(short_coefficient) / FixedMultiplier);
}
return coefficients;
}
FIRFilter::FIRFilter(std::size_t number_of_taps, float input_sample_rate, float low_frequency, float high_frequency, float attenuation) {
// we must be asked to filter based on an odd number of
// taps, and at least three
if(number_of_taps < 3) number_of_taps = 3;
if(attenuation < 21.0f) attenuation = 21.0f;
// ensure we have an odd number of taps
number_of_taps |= 1;
// store instance variables
filter_coefficients_.resize(number_of_taps);
FIRFilter::FIRFilter(
std::size_t number_of_taps,
const float input_sample_rate,
const float low_frequency,
float high_frequency,
float attenuation
) {
// Ensure an odd number of taps greater than or equal to 3, with a minimum attenuation of 21.
number_of_taps = std::max<size_t>(3, number_of_taps) | 1;
attenuation = std::max(attenuation, 21.0f);
/* calculate idealised filter response */
std::size_t Np = (number_of_taps - 1) / 2;
float two_over_sample_rate = 2.0f / input_sample_rate;
const std::size_t Np = (number_of_taps - 1) / 2;
const float two_over_sample_rate = 2.0f / input_sample_rate;
// Clamp the high cutoff frequency.
high_frequency = std::min(high_frequency, input_sample_rate * 0.5f);
@@ -125,7 +134,7 @@ FIRFilter::FIRFilter(std::size_t number_of_taps, float input_sample_rate, float
std::vector<float> A(Np+1);
A[0] = 2.0f * (high_frequency - low_frequency) / input_sample_rate;
for(unsigned int i = 1; i <= Np; ++i) {
float i_pi = float(i) * float(M_PI);
const float i_pi = float(i) * float(M_PI);
A[i] =
(
sinf(two_over_sample_rate * i_pi * high_frequency) -
@@ -133,20 +142,22 @@ FIRFilter::FIRFilter(std::size_t number_of_taps, float input_sample_rate, float
) / i_pi;
}
FIRFilter::coefficients_for_idealised_filter_response(filter_coefficients_.data(), A.data(), attenuation, number_of_taps);
filter_coefficients_ = coefficients_for_idealised_filter_response(A, attenuation, number_of_taps);
}
FIRFilter::FIRFilter(const std::vector<float> &coefficients) {
filter_coefficients_.reserve(coefficients.size());
for(const auto coefficient: coefficients) {
filter_coefficients_.push_back(short(coefficient * FixedMultiplier));
}
}
FIRFilter FIRFilter::operator+(const FIRFilter &rhs) const {
std::vector<float> coefficients = get_coefficients();
std::vector<float> rhs_coefficients = rhs.get_coefficients();
const auto coefficients = get_coefficients();
const auto rhs_coefficients = rhs.get_coefficients();
std::vector<float> sum;
sum.reserve(coefficients.size());
for(std::size_t i = 0; i < coefficients.size(); ++i) {
sum.push_back((coefficients[i] + rhs_coefficients[i]) / 2.0f);
}
@@ -155,9 +166,11 @@ FIRFilter FIRFilter::operator+(const FIRFilter &rhs) const {
}
FIRFilter FIRFilter::operator-() const {
const auto coefficients = get_coefficients();
std::vector<float> negative_coefficients;
for(const auto coefficient: get_coefficients()) {
negative_coefficients.reserve(coefficients.size());
for(const auto coefficient: coefficients) {
negative_coefficients.push_back(1.0f - coefficient);
}
@@ -165,10 +178,11 @@ FIRFilter FIRFilter::operator-() const {
}
FIRFilter FIRFilter::operator*(const FIRFilter &rhs) const {
std::vector<float> coefficients = get_coefficients();
std::vector<float> rhs_coefficients = rhs.get_coefficients();
const std::vector<float> coefficients = get_coefficients();
const std::vector<float> rhs_coefficients = rhs.get_coefficients();
std::vector<float> sum;
sum.reserve(coefficients.size());
for(std::size_t i = 0; i < coefficients.size(); ++i) {
sum.push_back(coefficients[i] * rhs_coefficients[i]);
}