Skip to content

Commit

Permalink
Merge pull request #62 from petiaccja/code-consolidation
Browse files Browse the repository at this point in the history
code consolidation #1
  • Loading branch information
petiaccja committed Jun 18, 2022
2 parents 2bcf452 + 96a7daf commit 04ac5a4
Show file tree
Hide file tree
Showing 20 changed files with 1,105 additions and 1,189 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/clang_format.yml
Original file line number Diff line number Diff line change
Expand Up @@ -23,5 +23,5 @@ jobs:
run: |
python ./support/run-clang-format.py -r --exclude "*pocketfft*" include
python ./support/run-clang-format.py -r test
python ./support/run-clang-format.py -r --exclude "*RtAudio*" examples
python ./support/run-clang-format.py -r --exclude "*RtAudio*" --exclude "*/out/*" examples
python ./support/run-clang-format.py -r benchmark
6 changes: 3 additions & 3 deletions benchmark/Bench_ApplyFilter.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -68,7 +68,7 @@ class FirFilterFixture : public celero::TestFixture {


template <class T, int64_t MaxOrder>
class IirFilterFixture : public celero::TestFixture {
class DesignFilterFixture : public celero::TestFixture {
public:
std::vector<ExperimentValue> getExperimentValues() const override {
std::vector<ExperimentValue> experimentValues;
Expand Down Expand Up @@ -118,8 +118,8 @@ class IirFilterFixture : public celero::TestFixture {
using BaselineFixture = FirFilterFixture<float, 1, 1, 8192>;
using ConvFixture = FirFilterFixture<float, 1, maxFirOrder, 16>;
using OlaFixture = FirFilterFixture<float, 32, maxFirOrder, 16>;
using TfFixture = IirFilterFixture<float, maxIirDirectOrder>;
using CascadeFixture = IirFilterFixture<float, maxIirCascadeOrder>;
using TfFixture = DesignFilterFixture<float, maxIirDirectOrder>;
using CascadeFixture = DesignFilterFixture<float, maxIirCascadeOrder>;

BASELINE_F(ApplyFilter, gain, BaselineFixture, 25, 1) {
Multiply(AsView(out).SubSignal(0, signal.Size()), signal, filter[0]);
Expand Down
20 changes: 10 additions & 10 deletions examples/03_IIR_Filtering.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -86,24 +86,24 @@ Signal<float> DialTone(char character) {
// - We set a loose pass-band ripple because the precise magnitude of the picked-up tone is not important
// - We set a strict stop-band ripple to heavily suppress noise outside the narrow band
constexpr int filterOrder = 6;
const auto filterDesc = Bandpass(ELLIPTIC).PassbandRipple(0.15f).StopbandRipple(0.02f);
const auto filterDesc = Iir.Bandpass.Elliptic.PassbandRipple(0.15f).StopbandRipple(0.02f);
float Normalize(float f) { return NormalizedFrequency(f, sampleRate); }; // A little helper since we need normalized frequencies.

// IirFilter returns the zero-pole representation of the designed filter. To apply it
// DesignFilter returns the zero-pole representation of the designed filter. To apply it
// to a signal, you have to convert it to a transfer function or a cascaded biquad.
// Unless you have a good reason, use cascaded biquads for their superior stability and accuracy.
const std::array<CascadedBiquad<float>, 4> filterBank1 = {
CascadedBiquad{ IirFilter<float>(filterOrder, filterDesc.Band(Normalize(frequencies1[0] - 10.f), Normalize(frequencies1[0] + 10.f))) },
CascadedBiquad{ IirFilter<float>(filterOrder, filterDesc.Band(Normalize(frequencies1[1] - 10.f), Normalize(frequencies1[1] + 10.f))) },
CascadedBiquad{ IirFilter<float>(filterOrder, filterDesc.Band(Normalize(frequencies1[2] - 10.f), Normalize(frequencies1[2] + 10.f))) },
CascadedBiquad{ IirFilter<float>(filterOrder, filterDesc.Band(Normalize(frequencies1[3] - 10.f), Normalize(frequencies1[3] + 10.f))) },
CascadedBiquad{ DesignFilter<float>(filterOrder, filterDesc.Band(Normalize(frequencies1[0] - 10.f), Normalize(frequencies1[0] + 10.f))) },
CascadedBiquad{ DesignFilter<float>(filterOrder, filterDesc.Band(Normalize(frequencies1[1] - 10.f), Normalize(frequencies1[1] + 10.f))) },
CascadedBiquad{ DesignFilter<float>(filterOrder, filterDesc.Band(Normalize(frequencies1[2] - 10.f), Normalize(frequencies1[2] + 10.f))) },
CascadedBiquad{ DesignFilter<float>(filterOrder, filterDesc.Band(Normalize(frequencies1[3] - 10.f), Normalize(frequencies1[3] + 10.f))) },
};

const std::array<CascadedBiquad<float>, 4> filterBank2 = {
CascadedBiquad{ IirFilter<float>(filterOrder, filterDesc.Band(Normalize(frequencies2[0] - 10.f), Normalize(frequencies2[0] + 10.f))) },
CascadedBiquad{ IirFilter<float>(filterOrder, filterDesc.Band(Normalize(frequencies2[1] - 10.f), Normalize(frequencies2[1] + 10.f))) },
CascadedBiquad{ IirFilter<float>(filterOrder, filterDesc.Band(Normalize(frequencies2[2] - 10.f), Normalize(frequencies2[2] + 10.f))) },
CascadedBiquad{ IirFilter<float>(filterOrder, filterDesc.Band(Normalize(frequencies2[3] - 10.f), Normalize(frequencies2[3] + 10.f))) },
CascadedBiquad{ DesignFilter<float>(filterOrder, filterDesc.Band(Normalize(frequencies2[0] - 10.f), Normalize(frequencies2[0] + 10.f))) },
CascadedBiquad{ DesignFilter<float>(filterOrder, filterDesc.Band(Normalize(frequencies2[1] - 10.f), Normalize(frequencies2[1] + 10.f))) },
CascadedBiquad{ DesignFilter<float>(filterOrder, filterDesc.Band(Normalize(frequencies2[2] - 10.f), Normalize(frequencies2[2] + 10.f))) },
CascadedBiquad{ DesignFilter<float>(filterOrder, filterDesc.Band(Normalize(frequencies2[3] - 10.f), Normalize(frequencies2[3] + 10.f))) },
};


Expand Down
4 changes: 2 additions & 2 deletions examples/04_FIR_Filtering.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
#include "Tools/LoadSound.hpp"
#include "Tools/PlaySound.hpp"

#include "dspbb/Filtering/FilterParameters.hpp"
#include "dspbb/Filtering/MeasureFilter.hpp"
#include <dspbb/Filtering/FIR.hpp>

#include <iostream>
Expand Down Expand Up @@ -51,7 +51,7 @@ class Equalizer {
void SetLevels(float bass, float mid, float treble) {
const auto normalizedResponse = [&](float nf) { return EqualizedResponse(nf * float(m_sampleRate) / 2.0f, bass, mid, treble); };
// We will use a least squares FIR filter with no weighting and default grid size.
FirFilter(m_filter, Arbitrary(LEAST_SQUARES).Response(normalizedResponse));
DesignFilter(m_filter, Fir.Arbitrary.LeastSquares.Response(normalizedResponse));
}

void Filter(SignalView<const float> leftIn, SignalView<const float> rightIn, SignalView<float> leftOut, SignalView<float> rightOut) {
Expand Down
2 changes: 1 addition & 1 deletion examples/Tools/LoadSound.cpp
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
#include "LoadSound.hpp"

#include <dspbb/Filtering/Interpolation.hpp>
#include <dspbb/Filtering/Resample.hpp>

#include <sndfile.hh>

Expand Down
177 changes: 86 additions & 91 deletions include/dspbb/Filtering/FIR.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -18,36 +18,36 @@ namespace dspbb {

// Lowpass
template <class SignalR, class ParamType, class WindowType, std::enable_if_t<is_mutable_signal_v<SignalR>, int> = 0>
void FirFilter(SignalR&& out, const impl::LowpassDesc<impl::FirMethodWindowed, ParamType, WindowType>& desc) {
void DesignFilter(SignalR&& out, const impl::windowed::LowpassDesc<ParamType, WindowType>& desc) {
fir::KernelWindowedLowpass(out, desc.cutoff, desc.window);
}

// Highpass
template <class SignalR, class ParamType, class WindowType, std::enable_if_t<is_mutable_signal_v<SignalR>, int> = 0>
void FirFilter(SignalR&& out, const impl::HighpassDesc<impl::FirMethodWindowed, ParamType, WindowType>& desc) {
FirFilter(out, Lowpass(WINDOWED).Cutoff(desc.cutoff).Window(desc.window));
void DesignFilter(SignalR&& out, const impl::windowed::HighpassDesc<ParamType, WindowType>& desc) {
DesignFilter(out, Fir.Lowpass.Windowed.Cutoff(desc.cutoff).Window(desc.window));
fir::ComplementaryResponse(out, out);
}

// Bandpass
template <class SignalR, class ParamType, class WindowType, std::enable_if_t<is_mutable_signal_v<SignalR>, int> = 0>
void FirFilter(SignalR&& out, const impl::BandpassDesc<impl::FirMethodWindowed, ParamType, WindowType>& desc) {
void DesignFilter(SignalR&& out, const impl::windowed::BandpassDesc<ParamType, WindowType>& desc) {
const ParamType bandWidth = desc.upper - desc.lower;
const ParamType bandCenter = (desc.upper + desc.lower) / ParamType(2);
FirFilter(out, Lowpass(WINDOWED).Cutoff(bandWidth / ParamType(2)).Window(desc.window));
DesignFilter(out, Fir.Lowpass.Windowed.Cutoff(bandWidth / ParamType(2)).Window(desc.window));
fir::ShiftResponse(out, out, bandCenter);
}

// Bandstop
template <class SignalR, class ParamType, class WindowType, std::enable_if_t<is_mutable_signal_v<SignalR>, int> = 0>
void FirFilter(SignalR&& out, const impl::BandstopDesc<impl::FirMethodWindowed, ParamType, WindowType>& desc) {
FirFilter(out, Bandpass(WINDOWED).Band(desc.lower, desc.upper).Window(desc.window));
void DesignFilter(SignalR&& out, const impl::windowed::BandstopDesc<ParamType, WindowType>& desc) {
DesignFilter(out, Fir.Bandpass.Windowed.Band(desc.lower, desc.upper).Window(desc.window));
fir::ComplementaryResponse(out, out);
}

// Arbitrary
template <class SignalR, class ResponseFunc, class WindowType, std::enable_if_t<is_mutable_signal_v<SignalR>, int> = 0>
void FirFilter(SignalR&& out, const impl::ArbitraryDesc<impl::FirMethodWindowed, ResponseFunc, WindowType>& desc) {
void DesignFilter(SignalR&& out, const impl::windowed::ArbitraryDesc<ResponseFunc, WindowType>& desc) {
fir::KernelWindowedArbitrary(out, desc.responseFunc, desc.window);
}

Expand Down Expand Up @@ -95,70 +95,64 @@ namespace impl {
}
return static_cast<F>(desc.weightHigh);
}
} // namespace impl

template <class ParamType>
auto TranslateLeastSquares(const LowpassDesc<FirMethodLeastSquares, ParamType>& desc) {
const auto response = [desc](auto f) {
using F = std::decay_t<decltype(f)>;
return Smoothstep(LerpParam(f, F(desc.cutoffEnd), F(desc.cutoffBegin)));
};
const auto weight = [desc](auto f) {
return LeastSquaresSplitWeight(f, desc);
};
return std::make_tuple(response, weight);
}

template <class ParamType>
auto TranslateLeastSquares(const HighpassDesc<FirMethodLeastSquares, ParamType>& desc) {
const auto response = [desc](auto f) {
using F = std::decay_t<decltype(f)>;
return Smoothstep(LerpParam(f, F(desc.cutoffBegin), F(desc.cutoffEnd)));
};
const auto weight = [desc](auto f) {
return LeastSquaresSplitWeight(f, desc);
};
return std::make_tuple(response, weight);
}

template <class ParamType>
auto TranslateLeastSquares(const BandpassDesc<FirMethodLeastSquares, ParamType>& desc) {
const ParamType fmid = (desc.lowerEnd + desc.upperBegin) / ParamType(2);
const auto response = [desc, fmid](auto f) {
using F = std::decay_t<decltype(f)>;
return f < fmid ? Smoothstep(LerpParam(f, F(desc.lowerBegin), F(desc.lowerEnd)))
: Smoothstep(LerpParam(f, F(desc.upperEnd), F(desc.upperBegin)));
};
const auto weight = [desc](auto f) {
return LeastSquaresBandWeight(f, desc);
};
return std::make_tuple(response, weight);
}
template <class SignalR, class ParamType, std::enable_if_t<is_mutable_signal_v<SignalR>, int> = 0>
auto DesignFilter(SignalR&& out, const impl::least_squares::LowpassDesc<ParamType>& desc) {
const auto response = [desc](auto f) {
using F = std::decay_t<decltype(f)>;
return impl::Smoothstep(impl::LerpParam(f, F(desc.cutoffEnd), F(desc.cutoffBegin)));
};
const auto weight = [desc](auto f) {
return impl::LeastSquaresSplitWeight(f, desc);
};
fir::KernelLeastSquares(out, response, weight, desc.grid);
}

template <class ParamType>
auto TranslateLeastSquares(const BandstopDesc<FirMethodLeastSquares, ParamType>& desc) {
const ParamType fmid = (desc.lowerEnd + desc.upperBegin) / ParamType(2);
const auto response = [desc, fmid](auto f) {
using F = std::decay_t<decltype(f)>;
return f < fmid ? Smoothstep(LerpParam(f, F(desc.lowerEnd), F(desc.lowerBegin)))
: Smoothstep(LerpParam(f, F(desc.upperBegin), F(desc.upperEnd)));
};
const auto weight = [desc](auto f) {
return LeastSquaresBandWeight(f, desc);
};
return std::make_tuple(response, weight);
}
template <class SignalR, class ParamType, std::enable_if_t<is_mutable_signal_v<SignalR>, int> = 0>
auto DesignFilter(SignalR&& out, const impl::least_squares::HighpassDesc<ParamType>& desc) {
const auto response = [desc](auto f) {
using F = std::decay_t<decltype(f)>;
return impl::Smoothstep(impl::LerpParam(f, F(desc.cutoffBegin), F(desc.cutoffEnd)));
};
const auto weight = [desc](auto f) {
return impl::LeastSquaresSplitWeight(f, desc);
};
fir::KernelLeastSquares(out, response, weight, desc.grid);
}

template <class ResponseFunc, class WeightFunc>
auto TranslateLeastSquares(const ArbitraryDesc<FirMethodLeastSquares, ResponseFunc, WeightFunc>& desc) {
return std::make_tuple(desc.responseFunc, desc.weightFunc);
}
template <class SignalR, class ParamType, std::enable_if_t<is_mutable_signal_v<SignalR>, int> = 0>
auto DesignFilter(SignalR&& out, const impl::least_squares::BandpassDesc<ParamType>& desc) {
const ParamType fmid = (desc.lowerEnd + desc.upperBegin) / ParamType(2);
const auto response = [desc, fmid](auto f) {
using F = std::decay_t<decltype(f)>;
return f < fmid ? impl::Smoothstep(impl::LerpParam(f, F(desc.lowerBegin), F(desc.lowerEnd)))
: impl::Smoothstep(impl::LerpParam(f, F(desc.upperEnd), F(desc.upperBegin)));
};
const auto weight = [desc](auto f) {
return impl::LeastSquaresBandWeight(f, desc);
};
fir::KernelLeastSquares(out, response, weight, desc.grid);
}

} // namespace impl
template <class SignalR, class ParamType, std::enable_if_t<is_mutable_signal_v<SignalR>, int> = 0>
auto DesignFilter(SignalR&& out, const impl::least_squares::BandstopDesc<ParamType>& desc) {
const ParamType fmid = (desc.lowerEnd + desc.upperBegin) / ParamType(2);
const auto response = [desc, fmid](auto f) {
using F = std::decay_t<decltype(f)>;
return f < fmid ? impl::Smoothstep(impl::LerpParam(f, F(desc.lowerEnd), F(desc.lowerBegin)))
: impl::Smoothstep(impl::LerpParam(f, F(desc.upperBegin), F(desc.upperEnd)));
};
const auto weight = [desc](auto f) {
return impl::LeastSquaresBandWeight(f, desc);
};
fir::KernelLeastSquares(out, response, weight, desc.grid);
}

template <class SignalR, template <typename, typename...> class Desc, class... Params, std::enable_if_t<is_mutable_signal_v<SignalR>, int> = 0>
auto FirFilter(SignalR&& out, const Desc<impl::FirMethodLeastSquares, Params...>& desc)
-> decltype(void(impl::TranslateLeastSquares(desc))) {
const auto [response, weight] = impl::TranslateLeastSquares(desc);
template <class SignalR, class ResponseFunc, class WeightFunc, std::enable_if_t<is_mutable_signal_v<SignalR>, int> = 0>
auto DesignFilter(SignalR&& out, const impl::least_squares::ArbitraryDesc<ResponseFunc, WeightFunc>& desc) {
const auto& response = desc.responseFunc;
const auto& weight = desc.weightFunc;
fir::KernelLeastSquares(out, response, weight, desc.grid);
}

Expand All @@ -169,43 +163,44 @@ auto FirFilter(SignalR&& out, const Desc<impl::FirMethodLeastSquares, Params...>

namespace impl {

template <class WindowType>
auto TranslateHilbert2HalfbandDesc(const HilbertDesc<FirMethodWindowed, WindowType>& desc) {
return Lowpass(WINDOWED).Cutoff(0.5f).Window(desc.window);
}

template <class ParamType>
auto TranslateHilbert2HalfbandDesc(const HilbertDesc<FirMethodLeastSquares, ParamType>& desc) {
const ParamType transitionBand = desc.transitionWidth;
return Lowpass(LEAST_SQUARES).Cutoff(ParamType(0.5) - transitionBand, ParamType(0.5) + transitionBand);
template <class SignalR, class HalfbandDesc, std::enable_if_t<is_mutable_signal_v<SignalR>, int> = 0>
void FirFilterHilbert(SignalR&& out, const HalfbandDesc& halfbandDesc) {
if (out.Size() % 2 == 0) {
const size_t halfbandSize = out.Size() * 2 - 1;
std::decay_t<SignalR> halfband(halfbandSize);
DesignFilter(halfband, halfbandDesc);
fir::HalfbandToHilbertEven(out, halfband);
}
else {
DesignFilter(out, halfbandDesc);
fir::HalfbandToHilbertOdd(out, out);
}
}

} // namespace impl

template <class SignalR, class Method, class... Params, std::enable_if_t<is_mutable_signal_v<SignalR>, int> = 0>
void FirFilter(SignalR&& out, const impl::HilbertDesc<Method, Params...>& desc) {
const auto halfbandDesc = impl::TranslateHilbert2HalfbandDesc(desc);
template <class SignalR, class WindowType, std::enable_if_t<is_mutable_signal_v<SignalR>, int> = 0>
auto DesignFilter(SignalR&& out, const impl::windowed::HilbertDesc<WindowType>& desc) {
const auto halfband = Fir.Lowpass.Windowed.Cutoff(0.5f).Window(desc.window);
impl::FirFilterHilbert(out, halfband);
}

if (out.Size() % 2 == 0) {
const size_t halfbandSize = out.Size() * 2 - 1;
std::decay_t<SignalR> halfband(halfbandSize);
FirFilter(halfband, halfbandDesc);
fir::HalfbandToHilbertEven(out, halfband);
}
else {
FirFilter(out, halfbandDesc);
fir::HalfbandToHilbertOdd(out, out);
}
template <class SignalR, class ParamType, std::enable_if_t<is_mutable_signal_v<SignalR>, int> = 0>
auto DesignFilter(SignalR&& out, const impl::least_squares::HilbertDesc<ParamType>& desc) {
const ParamType transitionBand = desc.transitionWidth;
const auto halfband = Fir.Lowpass.LeastSquares.Cutoff(ParamType(0.5) - transitionBand, ParamType(0.5) + transitionBand);
impl::FirFilterHilbert(out, halfband);
}


//------------------------------------------------------------------------------
// Out-of-place wrapper.
//------------------------------------------------------------------------------

template <class T, eSignalDomain Domain, class ResponseDesc>
auto FirFilter(size_t taps, ResponseDesc responseDesc) {
auto DesignFilter(size_t taps, ResponseDesc responseDesc) {
BasicSignal<T, Domain> out(taps);
FirFilter(out, responseDesc);
DesignFilter(out, responseDesc);
return out;
}

Expand Down
Loading

0 comments on commit 04ac5a4

Please sign in to comment.