Skip to content
Permalink
Browse files

4.0 (WIP)

  • Loading branch information
dlevin256 committed Nov 25, 2019
1 parent 4125bfd commit 7c4bee81d8a99fb745ae7fe1ed8e4365a741c972
@@ -12,7 +12,7 @@

def gen_ticks(stop, start=10):
yield start
for s in range(0, 10):
for s in range(1, 10):
if start * s > stop:
yield stop
raise StopIteration
@@ -21,8 +21,8 @@ def gen_ticks(stop, start=10):
yield t

def gen_tick_labels(stop, start=10):
yield (str(start) + 'Hz').replace('000Hz', 'kHz')
for s in range(0, 10):
yield (str(int(start)) + 'Hz').replace('000Hz', 'kHz')
for s in range(1, 10):
if start * s > stop:
yield (str(int(stop)) + 'Hz').replace('000Hz', 'kHz')
raise StopIteration
@@ -153,7 +153,7 @@ def plot(data,
style = {'linewidth': 1.4, 'color': '#0072bd'}
grid_style = {'color': '#777777'}

dataplot = a[0] if freqresp or phaseresp else a
dataplot = a[0] if freqresp or phaseresp else a

dataplot.plot(np.linspace(0, n, n, False), data, marker='.' if dots else None, **style)
dataplot.set_xlabel('Samples')
@@ -180,11 +180,16 @@ def set_freq(a):
if normalized_freq:
a.set_xlabel(freq_label[0])
X = np.linspace(0, 1, len(Y), False)
a.set_xlim([0, 1])
if log_freq:
a.set_xscale('log')
a.set_xlim([0.01, 1])
else:
a.set_xlim([0, 1])
else:
a.set_xlabel(freq_label[1])
if log_freq:
a.set_xscale('log')

a.set_xticks(list(gen_ticks(Fs / 2)))
a.set_xticklabels(list(gen_tick_labels(Fs / 2)))
X = np.linspace(0, Fs / 2, len(Y), False)
@@ -26,6 +26,9 @@ file(MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/svg)
add_executable(biquads biquads.cpp)
target_link_libraries(biquads kfr use_arch)

add_executable(iir iir.cpp)
target_link_libraries(iir kfr use_arch)

add_executable(window window.cpp)
target_link_libraries(window kfr use_arch)

@@ -0,0 +1,93 @@
/**
* KFR (http://kfrlib.com)
* Copyright (C) 2016 D Levin
* See LICENSE.txt for details
*/

#include <kfr/base.hpp>
#include <kfr/dsp.hpp>
#include <kfr/io.hpp>

using namespace kfr;

int main()
{
println(library_version());

constexpr size_t maxorder = 32;

const std::string options = "phaseresp=True, log_freq=True, freq_dB_lim=(-160, 10), padwidth=8192";

univector<fbase, 1024> output;
{
zpk<fbase> filt = iir_lowpass(bessel<fbase>(24), 1000, 48000);
std::vector<biquad_params<fbase>> bqs = to_sos(filt);
output = biquad<maxorder>(bqs, unitimpulse());
}
plot_save("bessel_lowpass24", output, options + ", title='24th-order Bessel filter, lowpass 1khz'");

{
zpk<fbase> filt = iir_lowpass(bessel<fbase>(12), 1000, 48000);
std::vector<biquad_params<fbase>> bqs = to_sos(filt);
output = biquad<maxorder>(bqs, unitimpulse());
}
plot_save("bessel_lowpass12", output, options + ", title='12th-order Bessel filter, lowpass 1khz'");

{
zpk<fbase> filt = iir_lowpass(bessel<fbase>(6), 1000, 48000);
std::vector<biquad_params<fbase>> bqs = to_sos(filt);
output = biquad<maxorder>(bqs, unitimpulse());
}
plot_save("bessel_lowpass6", output, options + ", title='6th-order Bessel filter, lowpass 1khz'");

{
zpk<fbase> filt = iir_lowpass(butterworth<fbase>(24), 1000, 48000);
std::vector<biquad_params<fbase>> bqs = to_sos(filt);
output = biquad<maxorder>(bqs, unitimpulse());
}
plot_save("butterworth_lowpass24", output,
options + ", title='24th-order Butterworth filter, lowpass 1khz'");

{
zpk<fbase> filt = iir_lowpass(butterworth<fbase>(12), 1000, 48000);
std::vector<biquad_params<fbase>> bqs = to_sos(filt);
output = biquad<maxorder>(bqs, unitimpulse());
}
plot_save("butterworth_lowpass12", output,
options + ", title='12th-order Butterworth filter, lowpass 1khz'");

{
zpk<fbase> filt = iir_highpass(butterworth<fbase>(12), 1000, 48000);
std::vector<biquad_params<fbase>> bqs = to_sos(filt);
output = biquad<maxorder>(bqs, unitimpulse());
}
plot_save("butterworth_highpass12", output,
options + ", title='12th-order Butterworth filter, highpass 1khz'");

{
zpk<fbase> filt = iir_bandpass(butterworth<fbase>(12), 0.1, 0.2);
std::vector<biquad_params<fbase>> bqs = to_sos(filt);
output = biquad<maxorder>(bqs, unitimpulse());
}
plot_save("butterworth_bandpass12", output,
options + ", title='12th-order Butterworth filter, bandpass'");

{
zpk<fbase> filt = iir_bandstop(butterworth<fbase>(12), 0.1, 0.2);
std::vector<biquad_params<fbase>> bqs = to_sos(filt);
output = biquad<maxorder>(bqs, unitimpulse());
}
plot_save("butterworth_bandstop12", output,
options + ", title='12th-order Butterworth filter, bandstop'");

{
zpk<fbase> filt = iir_bandpass(butterworth<fbase>(4), 0.005, 0.9);
std::vector<biquad_params<fbase>> bqs = to_sos(filt);
output = biquad<maxorder>(bqs, unitimpulse());
}
plot_save("butterworth_bandpass4", output, options + ", title='4th-order Butterworth filter, bandpass'");

println("SVG plots have been saved to svg directory");

return 0;
}
@@ -689,7 +689,7 @@ struct concatenate_expression : expression_with_arguments<E1, E2>
}
for (size_t i = size0 - index; i < N; ++i)
{
result[i] = self.argument(cinput, csize<1>, index + i, vec_shape<T, 1>{})[0];
result[i] = self.argument(cinput, csize<1>, index + i - size0, vec_shape<T, 1>{})[0];
}
return result;
}
@@ -333,26 +333,37 @@ struct expression_scalar : input_expression
}
};

template <typename, typename T, typename = void>
template <typename T1, typename T2, typename = void>
struct arg_impl
{
using type = T;
using type = T2;
using value_type = typename T1::value_type;
};

template <typename T1, typename T2>
struct arg_impl<T1, T2, void_t<enable_if<is_vec_element<T1>>>>
{
using type = expression_scalar<T1>;
using type = expression_scalar<T1>;
using value_type = T1;
};

template <typename T>
using arg = typename internal::arg_impl<decay<T>, T>::type;

template <typename T>
using arg_type = typename internal::arg_impl<decay<T>, T>::value_type;

template <typename Fn, typename... Args>
struct function_value_type
{
using type = typename invoke_result<Fn, vec<arg_type<Args>, 1>...>::value_type;
};

template <typename Fn, typename... Args>
struct expression_function : expression_with_arguments<arg<Args>...>
{
using value_type =
subtype<decltype(std::declval<Fn>()(std::declval<vec<value_type_of<arg<Args>>, 1>>()...))>;
using value_type = typename function_value_type<Fn, Args...>::type;
// subtype<decltype(std::declval<Fn>()(std::declval<vec<value_type_of<arg<Args>>, 1>>()...))>;
using T = value_type;

expression_function(Fn&& fn, arg<Args>&&... args) CMT_NOEXCEPT
@@ -1584,7 +1584,7 @@ inline size_t cfind(cvals_t<T, values...>, identity<T> value)
}

template <typename Fn, typename... Args>
CMT_NOINLINE static invoke_result<Fn, Args...> noinline(Fn&& fn, Args&&... args)
CMT_UNUSED CMT_NOINLINE static invoke_result<Fn, Args...> noinline(Fn&& fn, Args&&... args)
{
return fn(std::forward<Args>(args)...);
}
@@ -2083,9 +2083,12 @@ struct special_value
return static_cast<T>(d);
case special_constant::random_bits:
return random_bits<T>();
default:
return T{};
case special_constant::epsilon:
return std::numeric_limits<subtype<T>>::epsilon();
// default:
// return T{};
}
return T();
}

template <typename T>
@@ -422,6 +422,59 @@ struct dft_plan_real : dft_plan<T>
dft_pack_format) const = delete;
};

/// @brief DCT type 2 (unscaled)
template <typename T>
struct dct_plan : dft_plan<T>
{
dct_plan(size_t size) : dft_plan<T>(size) { this->temp_size += sizeof(complex<T>) * size * 2; }

KFR_MEM_INTRINSIC void execute(T* out, const T* in, u8* temp, bool inverse = false) const
{
const size_t size = this->size;
const size_t halfSize = size / 2;
univector_ref<complex<T>> mirrored = make_univector(
ptr_cast<complex<T>>(temp + this->temp_size - sizeof(complex<T>) * size * 2), size);
univector_ref<complex<T>> mirrored_dft =
make_univector(ptr_cast<complex<T>>(temp + this->temp_size - sizeof(complex<T>) * size), size);
auto t = counter() * c_pi<T> / (size * 2);
if (!inverse)
{
for (size_t i = 0; i < halfSize; i++)
{
mirrored[i] = in[i * 2];
mirrored[size - 1 - i] = in[i * 2 + 1];
}
if (size % 2)
{
mirrored[halfSize] = in[size - 1];
}
dft_plan<T>::execute(mirrored_dft.data(), mirrored.data(), temp, cfalse);
make_univector(out, size) = real(mirrored_dft) * cos(t) + imag(mirrored_dft) * sin(t);
}
else
{
mirrored = make_complex(make_univector(in, size) * cos(t), make_univector(in, size) * -sin(t));
dft_plan<T>::execute(mirrored_dft.data(), mirrored.data(), temp, cfalse);
for (size_t i = 0; i < halfSize; i++)
{
out[i * 2 + 0] = mirrored_dft[i].re;
out[i * 2 + 1] = mirrored_dft[size - 1 - i].re;
}
if (size % 2)
{
out[size - 1] = mirrored_dft[halfSize].re;
}
}
}

template <univector_tag Tag1, univector_tag Tag2, univector_tag Tag3>
KFR_MEM_INTRINSIC void execute(univector<T, Tag1>& out, const univector<T, Tag2>& in,
univector<u8, Tag3>& temp, bool inverse = false) const
{
execute(out.data(), in.data(), temp.data(), inverse);
}
};

inline namespace CMT_ARCH_NAME
{

@@ -685,7 +685,7 @@ template <bool inverse = false, typename T>
KFR_INTRINSIC void butterfly32(cvec<T, 32>& v32)
{
cvec<T, 4> w0, w1, w2, w3, w4, w5, w6, w7;
split<T, 64>(v32, w0, w1, w2, w3, w4, w5, w6, w7);
split(v32, w0, w1, w2, w3, w4, w5, w6, w7);
butterfly8<4, inverse>(w0, w1, w2, w3, w4, w5, w6, w7);

w1 = cmul(w1, fixed_twiddle<T, 4, 32, 0, 1, inverse>());
@@ -33,6 +33,7 @@
#include "dsp/fir_design.hpp"
#include "dsp/fracdelay.hpp"
#include "dsp/goertzel.hpp"
#include "dsp/iir_design.hpp"
#include "dsp/mixdown.hpp"
#include "dsp/oscillators.hpp"
#include "dsp/sample_rate_conversion.hpp"
@@ -321,6 +321,12 @@ KFR_FUNCTION expression_pointer<T> biquad(const biquad_params<T>* bq, size_t cou
[&] { return to_pointer(zeros<T>()); });
}

template <size_t maxfiltercount = 4, typename T, typename E1>
KFR_FUNCTION expression_pointer<T> biquad(const std::vector<biquad_params<T>>& bq, E1&& e1)
{
return biquad<maxfiltercount>(bq.data(), bq.size(), std::forward<E1>(e1));
}

template <typename T, size_t maxfiltercount = 4>
class biquad_filter : public expression_filter<T>
{

0 comments on commit 7c4bee8

Please sign in to comment.
You can’t perform that action at this time.