commit 7c4bee81d8a99fb745ae7fe1ed8e4365a741c972 parent 4125bfd07f9ec955a5866e2039c47ffa942da465 Author: [email protected] <[email protected]> Date: Mon, 25 Nov 2019 05:59:00 +0000 4.0 (WIP) Diffstat:
34 files changed, 1632 insertions(+), 124 deletions(-)
diff --git a/dspplot/dspplot/dspplotting.py b/dspplot/dspplot/dspplotting.py @@ -12,7 +12,7 @@ from scipy import interpolate 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 plot(data, 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) diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt @@ -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) diff --git a/examples/iir.cpp b/examples/iir.cpp @@ -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; +} diff --git a/include/kfr/base/basic_expressions.hpp b/include/kfr/base/basic_expressions.hpp @@ -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; } diff --git a/include/kfr/base/expression.hpp b/include/kfr/base/expression.hpp @@ -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 diff --git a/include/kfr/cometa.hpp b/include/kfr/cometa.hpp @@ -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> diff --git a/include/kfr/dft/fft.hpp b/include/kfr/dft/fft.hpp @@ -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 { diff --git a/include/kfr/dft/impl/ft.hpp b/include/kfr/dft/impl/ft.hpp @@ -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>()); diff --git a/include/kfr/dsp.hpp b/include/kfr/dsp.hpp @@ -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" diff --git a/include/kfr/dsp/biquad.hpp b/include/kfr/dsp/biquad.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> { diff --git a/include/kfr/dsp/iir_design.hpp b/include/kfr/dsp/iir_design.hpp @@ -0,0 +1,1167 @@ +/** @addtogroup biquad + * @{ + */ +/* + Copyright (C) 2016 D Levin (https://www.kfrlib.com) + This file is part of KFR + + KFR is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + KFR is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with KFR. + + If GPL is not suitable for your project, you must purchase a commercial license to use KFR. + Buying a commercial license is mandatory as soon as you develop commercial activities without + disclosing the source code of your own applications. + See https://www.kfrlib.com for details. + */ +#pragma once + +#include "../base/filter.hpp" +#include "../base/pointer.hpp" +#include "../simd/impl/function.hpp" +#include "../simd/operators.hpp" +#include "../simd/vec.hpp" +#include "../testo/assert.hpp" +#include "biquad_design.hpp" + +namespace kfr +{ + +template <typename T> +struct zpk +{ + univector<complex<T>> z; + univector<complex<T>> p; + T k; +}; + +inline namespace CMT_ARCH_NAME +{ + +template <typename T> +KFR_FUNCTION zpk<T> butterworth(int N) +{ + switch (N) + { + case 0: + return { {}, {}, 1 }; + case 1: + return { {}, { complex<T>(-1., -0.) }, 1 }; + case 2: + return { {}, + { complex<T>(-0.7071067811865476, -0.7071067811865476), + complex<T>(-0.7071067811865476, +0.7071067811865476) }, + 1 }; + case 3: + return { {}, + { complex<T>(-1., -0.), complex<T>(-0.5000000000000001, -0.8660254037844386), + complex<T>(-0.5000000000000001, +0.8660254037844386) }, + 1 }; + case 4: + return { {}, + { complex<T>(-0.9238795325112867, -0.3826834323650898), + complex<T>(-0.9238795325112867, +0.3826834323650898), + complex<T>(-0.38268343236508984, -0.9238795325112867), + complex<T>(-0.38268343236508984, +0.9238795325112867) }, + 1 }; + case 5: + return { {}, + { complex<T>(-1., -0.), complex<T>(-0.8090169943749475, -0.5877852522924731), + complex<T>(-0.8090169943749475, +0.5877852522924731), + complex<T>(-0.30901699437494745, -0.9510565162951535), + complex<T>(-0.30901699437494745, +0.9510565162951535) }, + 1 }; + case 6: + return { {}, + { complex<T>(-0.9659258262890683, -0.25881904510252074), + complex<T>(-0.9659258262890683, +0.25881904510252074), + complex<T>(-0.7071067811865476, -0.7071067811865476), + complex<T>(-0.7071067811865476, +0.7071067811865476), + complex<T>(-0.25881904510252096, -0.9659258262890682), + complex<T>(-0.25881904510252096, +0.9659258262890682) }, + 1 }; + case 7: + return { {}, + { complex<T>(-1., -0.), complex<T>(-0.9009688679024191, -0.4338837391175581), + complex<T>(-0.9009688679024191, +0.4338837391175581), + complex<T>(-0.6234898018587336, -0.7818314824680298), + complex<T>(-0.6234898018587336, +0.7818314824680298), + complex<T>(-0.22252093395631445, -0.9749279121818236), + complex<T>(-0.22252093395631445, +0.9749279121818236) }, + 1 }; + case 8: + return { {}, + { complex<T>(-0.9807852804032304, -0.19509032201612825), + complex<T>(-0.9807852804032304, +0.19509032201612825), + complex<T>(-0.8314696123025452, -0.5555702330196022), + complex<T>(-0.8314696123025452, +0.5555702330196022), + complex<T>(-0.5555702330196023, -0.8314696123025452), + complex<T>(-0.5555702330196023, +0.8314696123025452), + complex<T>(-0.19509032201612833, -0.9807852804032304), + complex<T>(-0.19509032201612833, +0.9807852804032304) }, + 1 }; + case 9: + return { {}, + { complex<T>(-1., -0.), complex<T>(-0.9396926207859084, -0.3420201433256687), + complex<T>(-0.9396926207859084, +0.3420201433256687), + complex<T>(-0.766044443118978, -0.6427876096865393), + complex<T>(-0.766044443118978, +0.6427876096865393), + complex<T>(-0.5000000000000001, -0.8660254037844386), + complex<T>(-0.5000000000000001, +0.8660254037844386), + complex<T>(-0.17364817766693041, -0.984807753012208), + complex<T>(-0.17364817766693041, +0.984807753012208) }, + 1 }; + case 10: + return { {}, + { complex<T>(-0.9876883405951378, -0.15643446504023087), + complex<T>(-0.9876883405951378, +0.15643446504023087), + complex<T>(-0.8910065241883679, -0.45399049973954675), + complex<T>(-0.8910065241883679, +0.45399049973954675), + complex<T>(-0.7071067811865476, -0.7071067811865476), + complex<T>(-0.7071067811865476, +0.7071067811865476), + complex<T>(-0.4539904997395468, -0.8910065241883678), + complex<T>(-0.4539904997395468, +0.8910065241883678), + complex<T>(-0.15643446504023092, -0.9876883405951378), + complex<T>(-0.15643446504023092, +0.9876883405951378) }, + 1 }; + case 11: + return { {}, + { complex<T>(-1., -0.), complex<T>(-0.9594929736144974, -0.28173255684142967), + complex<T>(-0.9594929736144974, +0.28173255684142967), + complex<T>(-0.8412535328311812, -0.5406408174555976), + complex<T>(-0.8412535328311812, +0.5406408174555976), + complex<T>(-0.654860733945285, -0.7557495743542583), + complex<T>(-0.654860733945285, +0.7557495743542583), + complex<T>(-0.41541501300188644, -0.9096319953545183), + complex<T>(-0.41541501300188644, +0.9096319953545183), + complex<T>(-0.14231483827328512, -0.9898214418809327), + complex<T>(-0.14231483827328512, +0.9898214418809327) }, + 1 }; + case 12: + return { {}, + { complex<T>(-0.9914448613738104, -0.13052619222005157), + complex<T>(-0.9914448613738104, +0.13052619222005157), + complex<T>(-0.9238795325112867, -0.3826834323650898), + complex<T>(-0.9238795325112867, +0.3826834323650898), + complex<T>(-0.7933533402912353, -0.6087614290087205), + complex<T>(-0.7933533402912353, +0.6087614290087205), + complex<T>(-0.6087614290087207, -0.7933533402912352), + complex<T>(-0.6087614290087207, +0.7933533402912352), + complex<T>(-0.38268343236508984, -0.9238795325112867), + complex<T>(-0.38268343236508984, +0.9238795325112867), + complex<T>(-0.13052619222005193, -0.9914448613738104), + complex<T>(-0.13052619222005193, +0.9914448613738104) }, + 1 }; + case 13: + return { {}, + { complex<T>(-1., -0.), complex<T>(-0.970941817426052, -0.23931566428755777), + complex<T>(-0.970941817426052, +0.23931566428755777), + complex<T>(-0.8854560256532099, -0.46472317204376856), + complex<T>(-0.8854560256532099, +0.46472317204376856), + complex<T>(-0.7485107481711011, -0.6631226582407952), + complex<T>(-0.7485107481711011, +0.6631226582407952), + complex<T>(-0.5680647467311558, -0.8229838658936565), + complex<T>(-0.5680647467311558, +0.8229838658936565), + complex<T>(-0.35460488704253557, -0.9350162426854148), + complex<T>(-0.35460488704253557, +0.9350162426854148), + complex<T>(-0.120536680255323, -0.992708874098054), + complex<T>(-0.120536680255323, +0.992708874098054) }, + 1 }; + case 14: + return { {}, + { complex<T>(-0.9937122098932426, -0.11196447610330786), + complex<T>(-0.9937122098932426, +0.11196447610330786), + complex<T>(-0.9438833303083676, -0.3302790619551671), + complex<T>(-0.9438833303083676, +0.3302790619551671), + complex<T>(-0.8467241992282842, -0.5320320765153366), + complex<T>(-0.8467241992282842, +0.5320320765153366), + complex<T>(-0.7071067811865476, -0.7071067811865476), + complex<T>(-0.7071067811865476, +0.7071067811865476), + complex<T>(-0.5320320765153366, -0.8467241992282841), + complex<T>(-0.5320320765153366, +0.8467241992282841), + complex<T>(-0.3302790619551673, -0.9438833303083675), + complex<T>(-0.3302790619551673, +0.9438833303083675), + complex<T>(-0.11196447610330791, -0.9937122098932426), + complex<T>(-0.11196447610330791, +0.9937122098932426) }, + 1 }; + case 15: + return { {}, + { complex<T>(-1., -0.), complex<T>(-0.9781476007338057, -0.20791169081775931), + complex<T>(-0.9781476007338057, +0.20791169081775931), + complex<T>(-0.9135454576426009, -0.40673664307580015), + complex<T>(-0.9135454576426009, +0.40673664307580015), + complex<T>(-0.8090169943749475, -0.5877852522924731), + complex<T>(-0.8090169943749475, +0.5877852522924731), + complex<T>(-0.6691306063588582, -0.7431448254773941), + complex<T>(-0.6691306063588582, +0.7431448254773941), + complex<T>(-0.5000000000000001, -0.8660254037844386), + complex<T>(-0.5000000000000001, +0.8660254037844386), + complex<T>(-0.30901699437494745, -0.9510565162951535), + complex<T>(-0.30901699437494745, +0.9510565162951535), + complex<T>(-0.10452846326765346, -0.9945218953682733), + complex<T>(-0.10452846326765346, +0.9945218953682733) }, + 1 }; + case 16: + return { {}, + { complex<T>(-0.9951847266721969, -0.0980171403295606), + complex<T>(-0.9951847266721969, +0.0980171403295606), + complex<T>(-0.9569403357322088, -0.29028467725446233), + complex<T>(-0.9569403357322088, +0.29028467725446233), + complex<T>(-0.881921264348355, -0.47139673682599764), + complex<T>(-0.881921264348355, +0.47139673682599764), + complex<T>(-0.773010453362737, -0.6343932841636455), + complex<T>(-0.773010453362737, +0.6343932841636455), + complex<T>(-0.6343932841636455, -0.7730104533627369), + complex<T>(-0.6343932841636455, +0.7730104533627369), + complex<T>(-0.4713967368259978, -0.8819212643483549), + complex<T>(-0.4713967368259978, +0.8819212643483549), + complex<T>(-0.29028467725446233, -0.9569403357322089), + complex<T>(-0.29028467725446233, +0.9569403357322089), + complex<T>(-0.09801714032956077, -0.9951847266721968), + complex<T>(-0.09801714032956077, +0.9951847266721968) }, + 1 }; + case 17: + return { {}, + { complex<T>(-1., -0.), complex<T>(-0.9829730996839018, -0.18374951781657034), + complex<T>(-0.9829730996839018, +0.18374951781657034), + complex<T>(-0.9324722294043558, -0.3612416661871529), + complex<T>(-0.9324722294043558, +0.3612416661871529), + complex<T>(-0.8502171357296142, -0.5264321628773557), + complex<T>(-0.8502171357296142, +0.5264321628773557), + complex<T>(-0.7390089172206591, -0.6736956436465572), + complex<T>(-0.7390089172206591, +0.6736956436465572), + complex<T>(-0.6026346363792565, -0.7980172272802395), + complex<T>(-0.6026346363792565, +0.7980172272802395), + complex<T>(-0.4457383557765383, -0.8951632913550623), + complex<T>(-0.4457383557765383, +0.8951632913550623), + complex<T>(-0.273662990072083, -0.961825643172819), + complex<T>(-0.273662990072083, +0.961825643172819), + complex<T>(-0.09226835946330202, -0.9957341762950345), + complex<T>(-0.09226835946330202, +0.9957341762950345) }, + 1 }; + case 18: + return { {}, + { complex<T>(-0.9961946980917455, -0.08715574274765817), + complex<T>(-0.9961946980917455, +0.08715574274765817), + complex<T>(-0.9659258262890683, -0.25881904510252074), + complex<T>(-0.9659258262890683, +0.25881904510252074), + complex<T>(-0.9063077870366499, -0.4226182617406994), + complex<T>(-0.9063077870366499, +0.4226182617406994), + complex<T>(-0.8191520442889918, -0.573576436351046), + complex<T>(-0.8191520442889918, +0.573576436351046), + complex<T>(-0.7071067811865476, -0.7071067811865476), + complex<T>(-0.7071067811865476, +0.7071067811865476), + complex<T>(-0.5735764363510463, -0.8191520442889917), + complex<T>(-0.5735764363510463, +0.8191520442889917), + complex<T>(-0.42261826174069944, -0.9063077870366499), + complex<T>(-0.42261826174069944, +0.9063077870366499), + complex<T>(-0.25881904510252096, -0.9659258262890682), + complex<T>(-0.25881904510252096, +0.9659258262890682), + complex<T>(-0.08715574274765814, -0.9961946980917455), + complex<T>(-0.08715574274765814, +0.9961946980917455) }, + 1 }; + case 19: + return { {}, + { complex<T>(-1., -0.), complex<T>(-0.9863613034027223, -0.1645945902807339), + complex<T>(-0.9863613034027223, +0.1645945902807339), + complex<T>(-0.9458172417006346, -0.32469946920468346), + complex<T>(-0.9458172417006346, +0.32469946920468346), + complex<T>(-0.8794737512064891, -0.4759473930370735), + complex<T>(-0.8794737512064891, +0.4759473930370735), + complex<T>(-0.7891405093963936, -0.6142127126896678), + complex<T>(-0.7891405093963936, +0.6142127126896678), + complex<T>(-0.6772815716257411, -0.7357239106731316), + complex<T>(-0.6772815716257411, +0.7357239106731316), + complex<T>(-0.5469481581224269, -0.8371664782625285), + complex<T>(-0.5469481581224269, +0.8371664782625285), + complex<T>(-0.40169542465296953, -0.9157733266550574), + complex<T>(-0.40169542465296953, +0.9157733266550574), + complex<T>(-0.24548548714079924, -0.9694002659393304), + complex<T>(-0.24548548714079924, +0.9694002659393304), + complex<T>(-0.0825793454723324, -0.9965844930066698), + complex<T>(-0.0825793454723324, +0.9965844930066698) }, + 1 }; + case 20: + return { {}, + { complex<T>(-0.996917333733128, -0.07845909572784494), + complex<T>(-0.996917333733128, +0.07845909572784494), + complex<T>(-0.9723699203976766, -0.2334453638559054), + complex<T>(-0.9723699203976766, +0.2334453638559054), + complex<T>(-0.9238795325112867, -0.3826834323650898), + complex<T>(-0.9238795325112867, +0.3826834323650898), + complex<T>(-0.8526401643540923, -0.5224985647159488), + complex<T>(-0.8526401643540923, +0.5224985647159488), + complex<T>(-0.7604059656000309, -0.6494480483301837), + complex<T>(-0.7604059656000309, +0.6494480483301837), + complex<T>(-0.6494480483301837, -0.7604059656000309), + complex<T>(-0.6494480483301837, +0.7604059656000309), + complex<T>(-0.5224985647159487, -0.8526401643540923), + complex<T>(-0.5224985647159487, +0.8526401643540923), + complex<T>(-0.38268343236508984, -0.9238795325112867), + complex<T>(-0.38268343236508984, +0.9238795325112867), + complex<T>(-0.23344536385590525, -0.9723699203976767), + complex<T>(-0.23344536385590525, +0.9723699203976767), + complex<T>(-0.078459095727845, -0.996917333733128), + complex<T>(-0.078459095727845, +0.996917333733128) }, + 1 }; + case 21: + return { {}, + { complex<T>(-1., -0.), + complex<T>(-0.9888308262251285, -0.14904226617617441), + complex<T>(-0.9888308262251285, +0.14904226617617441), + complex<T>(-0.9555728057861408, -0.29475517441090415), + complex<T>(-0.9555728057861408, +0.29475517441090415), + complex<T>(-0.9009688679024191, -0.4338837391175581), + complex<T>(-0.9009688679024191, +0.4338837391175581), + complex<T>(-0.8262387743159949, -0.563320058063622), + complex<T>(-0.8262387743159949, +0.563320058063622), + complex<T>(-0.7330518718298263, -0.6801727377709194), + complex<T>(-0.7330518718298263, +0.6801727377709194), + complex<T>(-0.6234898018587336, -0.7818314824680298), + complex<T>(-0.6234898018587336, +0.7818314824680298), + complex<T>(-0.5000000000000001, -0.8660254037844386), + complex<T>(-0.5000000000000001, +0.8660254037844386), + complex<T>(-0.3653410243663952, -0.9308737486442041), + complex<T>(-0.3653410243663952, +0.9308737486442041), + complex<T>(-0.22252093395631445, -0.9749279121818236), + complex<T>(-0.22252093395631445, +0.9749279121818236), + complex<T>(-0.07473009358642439, -0.9972037971811801), + complex<T>(-0.07473009358642439, +0.9972037971811801) }, + 1 }; + case 22: + return { {}, + { complex<T>(-0.9974521146102535, -0.07133918319923234), + complex<T>(-0.9974521146102535, +0.07133918319923234), + complex<T>(-0.9771468659711595, -0.21256528955297668), + complex<T>(-0.9771468659711595, +0.21256528955297668), + complex<T>(-0.9369497249997617, -0.34946417959909837), + complex<T>(-0.9369497249997617, +0.34946417959909837), + complex<T>(-0.8776789895672557, -0.4792489867200568), + complex<T>(-0.8776789895672557, +0.4792489867200568), + complex<T>(-0.8005412409243604, -0.5992776665113468), + complex<T>(-0.8005412409243604, +0.5992776665113468), + complex<T>(-0.7071067811865477, -0.7071067811865475), + complex<T>(-0.7071067811865477, +0.7071067811865475), + complex<T>(-0.5992776665113468, -0.8005412409243605), + complex<T>(-0.5992776665113468, +0.8005412409243605), + complex<T>(-0.479248986720057, -0.8776789895672555), + complex<T>(-0.479248986720057, +0.8776789895672555), + complex<T>(-0.3494641795990982, -0.9369497249997618), + complex<T>(-0.3494641795990982, +0.9369497249997618), + complex<T>(-0.21256528955297682, -0.9771468659711595), + complex<T>(-0.21256528955297682, +0.9771468659711595), + complex<T>(-0.07133918319923235, -0.9974521146102535), + complex<T>(-0.07133918319923235, +0.9974521146102535) }, + 1 }; + case 23: + return { {}, + { complex<T>(-1., -0.), + complex<T>(-0.9906859460363308, -0.1361666490962466), + complex<T>(-0.9906859460363308, +0.1361666490962466), + complex<T>(-0.9629172873477994, -0.2697967711570243), + complex<T>(-0.9629172873477994, +0.2697967711570243), + complex<T>(-0.917211301505453, -0.39840108984624145), + complex<T>(-0.917211301505453, +0.39840108984624145), + complex<T>(-0.8544194045464886, -0.5195839500354336), + complex<T>(-0.8544194045464886, +0.5195839500354336), + complex<T>(-0.7757112907044199, -0.6310879443260528), + complex<T>(-0.7757112907044199, +0.6310879443260528), + complex<T>(-0.6825531432186541, -0.730835964278124), + complex<T>(-0.6825531432186541, +0.730835964278124), + complex<T>(-0.5766803221148672, -0.816969893010442), + complex<T>(-0.5766803221148672, +0.816969893010442), + complex<T>(-0.4600650377311522, -0.8878852184023752), + complex<T>(-0.4600650377311522, +0.8878852184023752), + complex<T>(-0.3348796121709863, -0.9422609221188204), + complex<T>(-0.3348796121709863, +0.9422609221188204), + complex<T>(-0.20345601305263397, -0.9790840876823228), + complex<T>(-0.20345601305263397, +0.9790840876823228), + complex<T>(-0.06824241336467123, -0.9976687691905392), + complex<T>(-0.06824241336467123, +0.9976687691905392) }, + 1 }; + case 24: + return { {}, + { complex<T>(-0.9978589232386035, -0.06540312923014306), + complex<T>(-0.9978589232386035, +0.06540312923014306), + complex<T>(-0.9807852804032304, -0.19509032201612825), + complex<T>(-0.9807852804032304, +0.19509032201612825), + complex<T>(-0.9469301294951057, -0.32143946530316153), + complex<T>(-0.9469301294951057, +0.32143946530316153), + complex<T>(-0.8968727415326884, -0.44228869021900125), + complex<T>(-0.8968727415326884, +0.44228869021900125), + complex<T>(-0.8314696123025452, -0.5555702330196022), + complex<T>(-0.8314696123025452, +0.5555702330196022), + complex<T>(-0.7518398074789775, -0.6593458151000687), + complex<T>(-0.7518398074789775, +0.6593458151000687), + complex<T>(-0.6593458151000688, -0.7518398074789774), + complex<T>(-0.6593458151000688, +0.7518398074789774), + complex<T>(-0.5555702330196024, -0.8314696123025451), + complex<T>(-0.5555702330196024, +0.8314696123025451), + complex<T>(-0.44228869021900125, -0.8968727415326884), + complex<T>(-0.44228869021900125, +0.8968727415326884), + complex<T>(-0.3214394653031617, -0.9469301294951056), + complex<T>(-0.3214394653031617, +0.9469301294951056), + complex<T>(-0.19509032201612833, -0.9807852804032304), + complex<T>(-0.19509032201612833, +0.9807852804032304), + complex<T>(-0.06540312923014327, -0.9978589232386035), + complex<T>(-0.06540312923014327, +0.9978589232386035) }, + 1 }; + default: + return { {}, {}, 1.0 }; + } +} + +template <typename T> +KFR_FUNCTION zpk<T> bessel(int N) +{ + switch (N) + { + case 0: + return { {}, {}, 1.0 }; + case 1: + return { {}, { complex<T>(-1., -0.) }, 1.0 }; + case 2: + return { {}, + { complex<T>(-0.8660254037844385, -0.4999999999999999), + complex<T>(-0.8660254037844385, +0.4999999999999999) }, + 1.0 }; + case 3: + return { {}, + { complex<T>(-0.9416000265332067, -0.), complex<T>(-0.7456403858480766, -0.7113666249728353), + complex<T>(-0.7456403858480766, +0.7113666249728353) }, + 1.0 }; + case 4: + return { {}, + { complex<T>(-0.904758796788245, -0.27091873300387465), + complex<T>(-0.904758796788245, +0.27091873300387465), + complex<T>(-0.6572111716718827, -0.830161435004873), + complex<T>(-0.6572111716718827, +0.830161435004873) }, + 1.0 }; + case 5: + return { {}, + { complex<T>(-0.92644207738776, -0.), complex<T>(-0.8515536193688396, -0.44271746394433276), + complex<T>(-0.8515536193688396, +0.44271746394433276), + complex<T>(-0.5905759446119192, -0.9072067564574549), + complex<T>(-0.5905759446119192, +0.9072067564574549) }, + 1.0 }; + case 6: + return { {}, + { complex<T>(-0.9093906830472273, -0.1856964396793047), + complex<T>(-0.9093906830472273, +0.1856964396793047), + complex<T>(-0.7996541858328287, -0.5621717346937316), + complex<T>(-0.7996541858328287, +0.5621717346937316), + complex<T>(-0.5385526816693108, -0.9616876881954276), + complex<T>(-0.5385526816693108, +0.9616876881954276) }, + 1.0 }; + case 7: + return { {}, + { complex<T>(-0.919487155649029, -0.), complex<T>(-0.8800029341523379, -0.32166527623077407), + complex<T>(-0.8800029341523379, +0.32166527623077407), + complex<T>(-0.7527355434093214, -0.6504696305522553), + complex<T>(-0.7527355434093214, +0.6504696305522553), + complex<T>(-0.4966917256672317, -1.0025085084544205), + complex<T>(-0.4966917256672317, +1.0025085084544205) }, + 1.0 }; + case 8: + return { {}, + { complex<T>(-0.9096831546652911, -0.1412437976671423), + complex<T>(-0.9096831546652911, +0.1412437976671423), + complex<T>(-0.8473250802359334, -0.42590175382729345), + complex<T>(-0.8473250802359334, +0.42590175382729345), + complex<T>(-0.7111381808485397, -0.71865173141084), + complex<T>(-0.7111381808485397, +0.71865173141084), + complex<T>(-0.46217404125321254, -1.0343886811269012), + complex<T>(-0.46217404125321254, +1.0343886811269012) }, + 1.0 }; + case 9: + return { {}, + { complex<T>(-0.9154957797499037, -0.), + complex<T>(-0.8911217017079759, -0.25265809345821644), + complex<T>(-0.8911217017079759, +0.25265809345821644), + complex<T>(-0.8148021112269014, -0.50858156896315), + complex<T>(-0.8148021112269014, +0.50858156896315), + complex<T>(-0.6743622686854763, -0.7730546212691185), + complex<T>(-0.6743622686854763, +0.7730546212691185), + complex<T>(-0.43314155615536226, -1.0600736701359301), + complex<T>(-0.43314155615536226, +1.0600736701359301) }, + 1.0 }; + case 10: + return { {}, + { complex<T>(-0.9091347320900505, -0.11395831373355114), + complex<T>(-0.9091347320900505, +0.11395831373355114), + complex<T>(-0.8688459641284763, -0.34300082337663096), + complex<T>(-0.8688459641284763, +0.34300082337663096), + complex<T>(-0.7837694413101445, -0.5759147538499948), + complex<T>(-0.7837694413101445, +0.5759147538499948), + complex<T>(-0.6417513866988322, -0.8175836167191022), + complex<T>(-0.6417513866988322, +0.8175836167191022), + complex<T>(-0.4083220732868863, -1.0812748428191246), + complex<T>(-0.4083220732868863, +1.0812748428191246) }, + 1.0 }; + case 11: + return { {}, + { complex<T>(-0.9129067244518985, -0.), + complex<T>(-0.8963656705721169, -0.20804803750710327), + complex<T>(-0.8963656705721169, +0.20804803750710327), + complex<T>(-0.8453044014712964, -0.4178696917801249), + complex<T>(-0.8453044014712964, +0.4178696917801249), + complex<T>(-0.7546938934722305, -0.6319150050721849), + complex<T>(-0.7546938934722305, +0.6319150050721849), + complex<T>(-0.6126871554915198, -0.8547813893314768), + complex<T>(-0.6126871554915198, +0.8547813893314768), + complex<T>(-0.3868149510055096, -1.0991174667631216), + complex<T>(-0.3868149510055096, +1.0991174667631216) }, + 1.0 }; + case 12: + return { {}, + { complex<T>(-0.9084478234140686, -0.09550636521345042), + complex<T>(-0.9084478234140686, +0.09550636521345042), + complex<T>(-0.880253434201683, -0.2871779503524227), + complex<T>(-0.880253434201683, +0.2871779503524227), + complex<T>(-0.8217296939939074, -0.48102121151006755), + complex<T>(-0.8217296939939074, +0.48102121151006755), + complex<T>(-0.7276681615395161, -0.6792961178764695), + complex<T>(-0.7276681615395161, +0.6792961178764695), + complex<T>(-0.5866369321861475, -0.8863772751320727), + complex<T>(-0.5866369321861475, +0.8863772751320727), + complex<T>(-0.3679640085526314, -1.1143735756415463), + complex<T>(-0.3679640085526314, +1.1143735756415463) }, + 1.0 }; + case 13: + return { {}, + { complex<T>(-0.9110914665984182, -0.), + complex<T>(-0.8991314665475196, -0.17683429561610436), + complex<T>(-0.8991314665475196, +0.17683429561610436), + complex<T>(-0.8625094198260553, -0.35474137311729914), + complex<T>(-0.8625094198260553, +0.35474137311729914), + complex<T>(-0.7987460692470972, -0.5350752120696802), + complex<T>(-0.7987460692470972, +0.5350752120696802), + complex<T>(-0.7026234675721276, -0.7199611890171304), + complex<T>(-0.7026234675721276, +0.7199611890171304), + complex<T>(-0.5631559842430193, -0.9135900338325103), + complex<T>(-0.5631559842430193, +0.9135900338325103), + complex<T>(-0.35127923233898156, -1.1275915483177048), + complex<T>(-0.35127923233898156, +1.1275915483177048) }, + 1.0 }; + case 14: + return { {}, + { complex<T>(-0.9077932138396491, -0.08219639941940154), + complex<T>(-0.9077932138396491, +0.08219639941940154), + complex<T>(-0.8869506674916446, -0.24700791787653334), + complex<T>(-0.8869506674916446, +0.24700791787653334), + complex<T>(-0.8441199160909851, -0.41316538251026935), + complex<T>(-0.8441199160909851, +0.41316538251026935), + complex<T>(-0.7766591387063624, -0.581917067737761), + complex<T>(-0.7766591387063624, +0.581917067737761), + complex<T>(-0.6794256425119227, -0.7552857305042031), + complex<T>(-0.6794256425119227, +0.7552857305042031), + complex<T>(-0.5418766775112306, -0.9373043683516926), + complex<T>(-0.5418766775112306, +0.9373043683516926), + complex<T>(-0.3363868224902031, -1.1391722978398593), + complex<T>(-0.3363868224902031, +1.1391722978398593) }, + 1.0 }; + case 15: + return { {}, + { complex<T>(-0.9097482363849062, -0.), complex<T>(-0.9006981694176978, -0.1537681197278439), + complex<T>(-0.9006981694176978, +0.1537681197278439), + complex<T>(-0.8731264620834984, -0.30823524705642674), + complex<T>(-0.8731264620834984, +0.30823524705642674), + complex<T>(-0.8256631452587148, -0.46423487527343266), + complex<T>(-0.8256631452587148, +0.46423487527343266), + complex<T>(-0.7556027168970723, -0.6229396358758262), + complex<T>(-0.7556027168970723, +0.6229396358758262), + complex<T>(-0.6579196593111002, -0.7862895503722519), + complex<T>(-0.6579196593111002, +0.7862895503722519), + complex<T>(-0.5224954069658334, -0.9581787261092526), + complex<T>(-0.5224954069658334, +0.9581787261092526), + complex<T>(-0.3229963059766445, -1.14941615458363), + complex<T>(-0.3229963059766445, +1.14941615458363) }, + 1.0 }; + case 16: + return { {}, + { complex<T>(-0.9072099595087003, -0.07214211304111734), + complex<T>(-0.9072099595087003, +0.07214211304111734), + complex<T>(-0.8911723070323642, -0.2167089659900575), + complex<T>(-0.8911723070323642, +0.2167089659900575), + complex<T>(-0.8584264231521322, -0.3621697271802063), + complex<T>(-0.8584264231521322, +0.3621697271802063), + complex<T>(-0.8074790293236005, -0.5092933751171799), + complex<T>(-0.8074790293236005, +0.5092933751171799), + complex<T>(-0.7356166304713119, -0.6591950877860395), + complex<T>(-0.7356166304713119, +0.6591950877860395), + complex<T>(-0.6379502514039066, -0.8137453537108762), + complex<T>(-0.6379502514039066, +0.8137453537108762), + complex<T>(-0.504760644442476, -0.9767137477799086), + complex<T>(-0.504760644442476, +0.9767137477799086), + complex<T>(-0.31087827556453856, -1.1585528411993304), + complex<T>(-0.31087827556453856, +1.1585528411993304) }, + 1.0 }; + case 17: + return { {}, + { complex<T>(-0.9087141161336388, -0.), + complex<T>(-0.9016273850787279, -0.13602679951730237), + complex<T>(-0.9016273850787279, +0.13602679951730237), + complex<T>(-0.8801100704438625, -0.2725347156478803), + complex<T>(-0.8801100704438625, +0.2725347156478803), + complex<T>(-0.8433414495836128, -0.41007592829100215), + complex<T>(-0.8433414495836128, +0.41007592829100215), + complex<T>(-0.7897644147799701, -0.5493724405281085), + complex<T>(-0.7897644147799701, +0.5493724405281085), + complex<T>(-0.7166893842372348, -0.6914936286393606), + complex<T>(-0.7166893842372348, +0.6914936286393606), + complex<T>(-0.6193710717342136, -0.8382497252826987), + complex<T>(-0.6193710717342136, +0.8382497252826987), + complex<T>(-0.48846293376727057, -0.9932971956316782), + complex<T>(-0.48846293376727057, +0.9932971956316782), + complex<T>(-0.29984894599900724, -1.1667612729256673), + complex<T>(-0.29984894599900724, +1.1667612729256673) }, + 1.0 }; + case 18: + return { {}, + { complex<T>(-0.9067004324162776, -0.06427924106393067), + complex<T>(-0.9067004324162776, +0.06427924106393067), + complex<T>(-0.8939764278132456, -0.19303746408947586), + complex<T>(-0.8939764278132456, +0.19303746408947586), + complex<T>(-0.8681095503628832, -0.32242049251632576), + complex<T>(-0.8681095503628832, +0.32242049251632576), + complex<T>(-0.8281885016242831, -0.45293856978159136), + complex<T>(-0.8281885016242831, +0.45293856978159136), + complex<T>(-0.7726285030739557, -0.5852778162086639), + complex<T>(-0.7726285030739557, +0.5852778162086639), + complex<T>(-0.698782144500527, -0.7204696509726628), + complex<T>(-0.698782144500527, +0.7204696509726628), + complex<T>(-0.6020482668090646, -0.8602708961893666), + complex<T>(-0.6020482668090646, +0.8602708961893666), + complex<T>(-0.4734268069916154, -1.0082343003148009), + complex<T>(-0.4734268069916154, +1.0082343003148009), + complex<T>(-0.28975920298804847, -1.1741830106000584), + complex<T>(-0.28975920298804847, +1.1741830106000584) }, + 1.0 }; + case 19: + return { {}, + { complex<T>(-0.9078934217899399, -0.), + complex<T>(-0.9021937639390656, -0.12195683818720263), + complex<T>(-0.9021937639390656, +0.12195683818720263), + complex<T>(-0.8849290585034385, -0.24425907575498182), + complex<T>(-0.8849290585034385, +0.24425907575498182), + complex<T>(-0.8555768765618422, -0.3672925896399872), + complex<T>(-0.8555768765618422, +0.3672925896399872), + complex<T>(-0.8131725551578203, -0.491536503556246), + complex<T>(-0.8131725551578203, +0.491536503556246), + complex<T>(-0.7561260971541627, -0.6176483917970176), + complex<T>(-0.7561260971541627, +0.6176483917970176), + complex<T>(-0.6818424412912442, -0.7466272357947761), + complex<T>(-0.6818424412912442, +0.7466272357947761), + complex<T>(-0.5858613321217832, -0.8801817131014564), + complex<T>(-0.5858613321217832, +0.8801817131014564), + complex<T>(-0.4595043449730983, -1.0217687769126707), + complex<T>(-0.4595043449730983, +1.0217687769126707), + complex<T>(-0.2804866851439361, -1.1809316284532905), + complex<T>(-0.2804866851439361, +1.1809316284532905) }, + 1.0 }; + case 20: + return { {}, + { complex<T>(-0.9062570115576768, -0.0579617802778495), + complex<T>(-0.9062570115576768, +0.0579617802778495), + complex<T>(-0.8959150941925766, -0.17403171759187044), + complex<T>(-0.8959150941925766, +0.17403171759187044), + complex<T>(-0.8749560316673335, -0.2905559296567909), + complex<T>(-0.8749560316673335, +0.2905559296567909), + complex<T>(-0.8427907479956664, -0.4078917326291931), + complex<T>(-0.8427907479956664, +0.4078917326291931), + complex<T>(-0.7984251191290602, -0.526494238881713), + complex<T>(-0.7984251191290602, +0.526494238881713), + complex<T>(-0.7402780309646764, -0.6469975237605227), + complex<T>(-0.7402780309646764, +0.6469975237605227), + complex<T>(-0.6658120544829932, -0.7703721701100759), + complex<T>(-0.6658120544829932, +0.7703721701100759), + complex<T>(-0.5707026806915716, -0.8982829066468254), + complex<T>(-0.5707026806915716, +0.8982829066468254), + complex<T>(-0.44657006982051484, -1.0340977025608422), + complex<T>(-0.44657006982051484, +1.0340977025608422), + complex<T>(-0.27192995802516506, -1.187099379810886), + complex<T>(-0.27192995802516506, +1.187099379810886) }, + 1.0 }; + case 21: + return { {}, + { complex<T>(-0.9072262653142963, -0.), + complex<T>(-0.9025428073192694, -0.11052525727898564), + complex<T>(-0.9025428073192694, +0.11052525727898564), + complex<T>(-0.8883808106664449, -0.221306921508435), + complex<T>(-0.8883808106664449, +0.221306921508435), + complex<T>(-0.8643915813643203, -0.33262585125221866), + complex<T>(-0.8643915813643203, +0.33262585125221866), + complex<T>(-0.8299435470674444, -0.44481777394079575), + complex<T>(-0.8299435470674444, +0.44481777394079575), + complex<T>(-0.7840287980408347, -0.5583186348022857), + complex<T>(-0.7840287980408347, +0.5583186348022857), + complex<T>(-0.7250839687106612, -0.6737426063024383), + complex<T>(-0.7250839687106612, +0.6737426063024383), + complex<T>(-0.6506315378609466, -0.7920349342629495), + complex<T>(-0.6506315378609466, +0.7920349342629495), + complex<T>(-0.5564766488918566, -0.9148198405846728), + complex<T>(-0.5564766488918566, +0.9148198405846728), + complex<T>(-0.4345168906815268, -1.045382255856986), + complex<T>(-0.4345168906815268, +1.045382255856986), + complex<T>(-0.2640041595834027, -1.192762031948052), + complex<T>(-0.2640041595834027, +1.192762031948052) }, + 1.0 }; + case 22: + return { {}, + { complex<T>(-0.9058702269930871, -0.05277490828999903), + complex<T>(-0.9058702269930871, +0.05277490828999903), + complex<T>(-0.8972983138153532, -0.15843519122898653), + complex<T>(-0.8972983138153532, +0.15843519122898653), + complex<T>(-0.8799661455640174, -0.26443630392015344), + complex<T>(-0.8799661455640174, +0.26443630392015344), + complex<T>(-0.8534754036851689, -0.37103893194823206), + complex<T>(-0.8534754036851689, +0.37103893194823206), + complex<T>(-0.8171682088462721, -0.47856194922027806), + complex<T>(-0.8171682088462721, +0.47856194922027806), + complex<T>(-0.7700332930556816, -0.5874255426351151), + complex<T>(-0.7700332930556816, +0.5874255426351151), + complex<T>(-0.7105305456418792, -0.6982266265924527), + complex<T>(-0.7105305456418792, +0.6982266265924527), + complex<T>(-0.6362427683267828, -0.811887504024635), + complex<T>(-0.6362427683267828, +0.811887504024635), + complex<T>(-0.5430983056306306, -0.9299947824439877), + complex<T>(-0.5430983056306306, +0.9299947824439877), + complex<T>(-0.4232528745642629, -1.055755605227546), + complex<T>(-0.4232528745642629, +1.055755605227546), + complex<T>(-0.2566376987939318, -1.1979824335552132), + complex<T>(-0.2566376987939318, +1.1979824335552132) }, + 1.0 }; + case 23: + return { {}, + { complex<T>(-0.9066732476324991, -0.), + complex<T>(-0.9027564979912508, -0.10105343353140452), + complex<T>(-0.9027564979912508, +0.10105343353140452), + complex<T>(-0.8909283242471254, -0.20230246993812237), + complex<T>(-0.8909283242471254, +0.20230246993812237), + complex<T>(-0.8709469395587415, -0.3039581993950041), + complex<T>(-0.8709469395587415, +0.3039581993950041), + complex<T>(-0.8423805948021129, -0.4062657948237603), + complex<T>(-0.8423805948021129, +0.4062657948237603), + complex<T>(-0.8045561642053178, -0.5095305912227259), + complex<T>(-0.8045561642053178, +0.5095305912227259), + complex<T>(-0.7564660146829886, -0.6141594859476034), + complex<T>(-0.7564660146829886, +0.6141594859476034), + complex<T>(-0.6965966033912708, -0.720734137475305), + complex<T>(-0.6965966033912708, +0.720734137475305), + complex<T>(-0.6225903228771341, -0.830155830281298), + complex<T>(-0.6225903228771341, +0.830155830281298), + complex<T>(-0.5304922463810198, -0.9439760364018306), + complex<T>(-0.5304922463810198, +0.9439760364018306), + complex<T>(-0.4126986617510148, -1.0653287944755134), + complex<T>(-0.4126986617510148, +1.0653287944755134), + complex<T>(-0.24976972022089572, -1.2028131878706978), + complex<T>(-0.24976972022089572, +1.2028131878706978) }, + 1.0 }; + case 24: + return { {}, + { complex<T>(-0.9055312363372773, -0.0484400665404787), + complex<T>(-0.9055312363372773, +0.0484400665404787), + complex<T>(-0.8983105104397872, -0.14540561338736102), + complex<T>(-0.8983105104397872, +0.14540561338736102), + complex<T>(-0.8837358034555707, -0.24263352344013836), + complex<T>(-0.8837358034555707, +0.24263352344013836), + complex<T>(-0.8615278304016355, -0.34032021126186246), + complex<T>(-0.8615278304016355, +0.34032021126186246), + complex<T>(-0.8312326466813242, -0.4386985933597306), + complex<T>(-0.8312326466813242, +0.4386985933597306), + complex<T>(-0.7921695462343489, -0.5380628490968016), + complex<T>(-0.7921695462343489, +0.5380628490968016), + complex<T>(-0.7433392285088533, -0.6388084216222569), + complex<T>(-0.7433392285088533, +0.6388084216222569), + complex<T>(-0.6832565803536519, -0.7415032695091649), + complex<T>(-0.6832565803536519, +0.7415032695091649), + complex<T>(-0.6096221567378332, -0.8470292433077199), + complex<T>(-0.6096221567378332, +0.8470292433077199), + complex<T>(-0.518591457482032, -0.9569048385259057), + complex<T>(-0.518591457482032, +0.9569048385259057), + complex<T>(-0.4027853855197519, -1.0741951965186751), + complex<T>(-0.4027853855197519, +1.0741951965186751), + complex<T>(-0.24334813375248746, -1.2072986837319728), + complex<T>(-0.24334813375248746, +1.2072986837319728) }, + 1.0 }; + + default: + return { {}, {}, 1.f }; + } +} + +// template <typename T> +// KFR_FUNCTION zpk<T> iirfilter(const zpk<T>& filter) +// { +// } + +namespace internal +{ + +template <typename T> +KFR_FUNCTION zpk<T> bilinear(const zpk<T>& filter, T fs) +{ + const T fs2 = 2.0 * fs; + zpk<T> result; + result.z = (fs2 + filter.z) / (fs2 - filter.z); + result.p = (fs2 + filter.p) / (fs2 - filter.p); + result.z.resize(result.p.size(), complex<T>(-1)); + result.k = filter.k * real(product(fs2 - filter.z) / product(fs2 - filter.p)); + return result; +} + +template <typename T> +struct zero_pole_pairs +{ + complex<T> p1, p2, z1, z2; +}; + +template <typename T> +KFR_FUNCTION vec<T, 3> zpk2tf_poly(const complex<T>& x, const complex<T>& y) +{ + return { T(1), -(x.re + y.re), x.re * y.re - x.im * y.im }; +} + +template <typename T> +KFR_FUNCTION biquad_params<T> zpk2tf(const zero_pole_pairs<T>& pairs, T k) +{ + // println("----zpk2tf"); + // println(fmt<'g', 24, 19>(pairs.z1), " | ", fmt<'g', 24, 19>(pairs.z2)); + // println(fmt<'g', 24, 19>(pairs.p1), " | ", fmt<'g', 24, 19>(pairs.p2)); + // println(k); + + vec<T, 3> zz = k * zpk2tf_poly(pairs.z1, pairs.z2); + vec<T, 3> pp = zpk2tf_poly(pairs.p1, pairs.p2); + // return { zz[0], zz[1], zz[2], pp[0], pp[1], pp[2] }; + return { pp[0], pp[1], pp[2], zz[0], zz[1], zz[2] }; +} + +template <typename T> +KFR_FUNCTION univector<complex<T>> cplxreal(const univector<complex<T>>& x) +{ + // conjugate pairs must be adjacent + T tol = std::numeric_limits<T>::epsilon() * 100; + univector<complex<T>> result = x; + for (size_t i = result.size(); i > 1; i--) + { + if (!isreal(result[i - 1]) && !isreal(result[i - 2])) + { + if (abs(result[i - 1].re - result[i - 2].re) < tol && + abs(result[i - 1].im + result[i - 2].im) < tol) + { + result.erase(result.begin() + i - 1); + result[i - 2].im = abs(result[i - 2].im); + } + } + } + return result; +} + +template <typename T> +KFR_FUNCTION size_t nearest_real_or_complex(const univector<complex<T>>& list, const complex<T>& val, + bool mustbereal = true) +{ + univector<complex<T>> filtered; + for (complex<T> v : list) + { + if (isreal(v) == mustbereal) + { + filtered.push_back(v); + } + } + TESTO_ASSERT(!filtered.empty()); + if (filtered.empty()) + return std::numeric_limits<size_t>::max(); + + size_t minidx = 0; + T minval = cabs(val - filtered[0]); + for (size_t i = 1; i < list.size(); i++) + { + T newminval = cabs(val - filtered[i]); + if (newminval < minval) + { + minval = newminval; + minidx = i; + } + } + return minidx; +} + +template <typename T> +KFR_FUNCTION int countreal(const univector<complex<T>>& list) +{ + int nreal = 0; + for (complex<T> c : list) + { + if (c.im == 0) + nreal++; + } + return nreal; +} + +template <typename T> +KFR_FUNCTION zpk<T> lp2lp_zpk(const zpk<T>& filter, T wo) +{ + zpk<T> result; + result.z = wo * filter.z; + result.p = wo * filter.p; + result.k = filter.k * pow(wo, filter.p.size() - filter.z.size()); + return result; +} + +template <typename T> +KFR_FUNCTION zpk<T> lp2hp_zpk(const zpk<T>& filter, identity<T> wo) +{ + zpk<T> result; + result.z = wo / filter.z; + result.p = wo / filter.p; + result.z.resize(result.p.size(), T(0)); + result.k = filter.k * real(product(-filter.z) / product(-filter.p)); + return result; +} + +template <typename T> +KFR_FUNCTION zpk<T> lp2bp_zpk(const zpk<T>& filter, identity<T> wo, identity<T> bw) +{ + zpk<T> lowpass; + lowpass.z = bw * 0.5 * filter.z; + lowpass.p = bw * 0.5 * filter.p; + + zpk<T> result; + result.z = concatenate(lowpass.z + csqrt(csqr(lowpass.z) - wo * wo), + lowpass.z - csqrt(csqr(lowpass.z) - wo * wo)); + result.p = concatenate(lowpass.p + csqrt(csqr(lowpass.p) - wo * wo), + lowpass.p - csqrt(csqr(lowpass.p) - wo * wo)); + + result.z.resize(result.z.size() + filter.p.size() - filter.z.size(), T(0)); + result.k = filter.k * pow(bw, filter.p.size() - filter.z.size()); + + return result; +} + +template <typename T> +KFR_FUNCTION zpk<T> lp2bs_zpk(const zpk<T>& filter, identity<T> wo, identity<T> bw) +{ + zpk<T> highpass; + highpass.z = (bw * 0.5) / filter.z; + highpass.p = (bw * 0.5) / filter.p; + + zpk<T> result; + result.z = concatenate(highpass.z + csqrt(csqr(highpass.z) - wo * wo), + highpass.z - csqrt(csqr(highpass.z) - wo * wo)); + result.p = concatenate(highpass.p + csqrt(csqr(highpass.p) - wo * wo), + highpass.p - csqrt(csqr(highpass.p) - wo * wo)); + + result.z.resize(result.z.size() + filter.p.size() - filter.z.size(), complex<T>(0, +wo)); + result.z.resize(result.z.size() + filter.p.size() - filter.z.size(), complex<T>(0, -wo)); + result.k = filter.k * real(product(-filter.z) / product(-filter.p)); + + return result; +} + +template <typename T> +KFR_FUNCTION T warp_freq(T frequency, T fs) +{ + frequency = 2 * frequency / fs; + fs = 2.0; + T warped = 2 * fs * tan(c_pi<T> * frequency / fs); + return warped; +} + +} // namespace internal + +template <typename T> +KFR_FUNCTION zpk<T> iir_lowpass(const zpk<T>& filter, identity<T> frequency, identity<T> fs = 2.0) +{ + T warped = internal::warp_freq(frequency, fs); + + zpk<T> result = filter; + result = internal::lp2lp_zpk(result, warped); + result = internal::bilinear(result, 2.0); // fs = 2.0 + return result; +} + +template <typename T> +KFR_FUNCTION zpk<T> iir_highpass(const zpk<T>& filter, identity<T> frequency, identity<T> fs = 2.0) +{ + T warped = internal::warp_freq(frequency, fs); + + zpk<T> result = filter; + result = internal::lp2hp_zpk(result, warped); + result = internal::bilinear(result, 2.0); // fs = 2.0 + return result; +} + +template <typename T> +KFR_FUNCTION zpk<T> iir_bandpass(const zpk<T>& filter, identity<T> lowfreq, identity<T> highfreq, + identity<T> fs = 2.0) +{ + T warpedlow = internal::warp_freq(lowfreq, fs); + T warpedhigh = internal::warp_freq(highfreq, fs); + + zpk<T> result = filter; + result = internal::lp2bp_zpk(result, std::sqrt(warpedlow * warpedhigh), warpedhigh - warpedlow); + result = internal::bilinear(result, 2.0); // fs = 2.0 + return result; +} + +template <typename T> +KFR_FUNCTION zpk<T> iir_bandstop(const zpk<T>& filter, identity<T> lowfreq, identity<T> highfreq, + identity<T> fs = 2.0) +{ + T warpedlow = internal::warp_freq(lowfreq, fs); + T warpedhigh = internal::warp_freq(highfreq, fs); + + zpk<T> result = filter; + result = internal::lp2bs_zpk(result, std::sqrt(warpedlow * warpedhigh), warpedhigh - warpedlow); + result = internal::bilinear(result, 2.0); // fs = 2.0 + return result; +} + +template <typename T> +KFR_FUNCTION std::vector<biquad_params<T>> to_sos(const zpk<T>& filter) +{ + if (filter.p.empty() && filter.z.empty()) + return { biquad_params<T>(filter.k, 0., 0., 1., 0., 0) }; + + zpk<T> filt = filter; + size_t length = std::max(filter.p.size(), filter.z.size()); + filt.p.resize(length, complex<T>(0)); + filt.z.resize(length, complex<T>(0)); + + size_t n_sections = (length + 1) / 2; + if (length & 1) + { + filt.z.push_back(complex<T>(0)); + filt.p.push_back(complex<T>(0)); + } + + filt.z = internal::cplxreal(filt.z); + filt.p = internal::cplxreal(filt.p); + std::vector<internal::zero_pole_pairs<T>> pairs(n_sections); + + for (size_t si = 0; si < n_sections; si++) + { + size_t worstidx = 0; + T worstval = abs(1 - cabs(filt.p[0])); + for (size_t i = 1; i < filt.p.size(); i++) + { + T val = abs(1 - cabs(filt.p[i])); + if (val < worstval) + { + worstidx = i; + worstval = val; + } + } + complex<T> p1 = filt.p[worstidx]; + filt.p.erase(filt.p.begin() + worstidx); + + complex<T> z1, p2, z2; + if (isreal(p1) && internal::countreal(filt.p) == 0) + { + size_t z1_idx = internal::nearest_real_or_complex(filt.z, p1, true); + z1 = filt.z[z1_idx]; + filt.z.erase(filt.z.begin() + z1_idx); + p2 = z2 = 0; + } + else + { + size_t z1_idx; + if (!isreal(p1) && internal::countreal(filt.z) == 1) + { + z1_idx = internal::nearest_real_or_complex(filt.z, p1, false); + } + else + { + size_t minidx = 0; + T minval = cabs(p1 - filt.z[0]); + for (size_t i = 1; i < filt.z.size(); i++) + { + T newminval = cabs(p1 - filt.z[i]); + if (newminval < minval) + { + minidx = i; + minval = newminval; + } + } + z1_idx = minidx; + } + z1 = filt.z[z1_idx]; + filt.z.erase(filt.z.begin() + z1_idx); + if (!isreal(p1)) + { + if (!isreal(z1)) + { + p2 = cconj(p1); + z2 = cconj(z1); + } + else + { + p2 = cconj(p1); + size_t z2_idx = internal::nearest_real_or_complex(filt.z, p1, true); + z2 = filt.z[z2_idx]; + TESTO_ASSERT(isreal(z2)); + filt.z.erase(filt.z.begin() + z2_idx); + } + } + else + { + size_t p2_idx; + size_t z2_idx; + if (!isreal(z1)) + { + z2 = cconj(z1); + p2_idx = internal::nearest_real_or_complex(filt.z, p1, true); + p2 = filt.p[p2_idx]; + TESTO_ASSERT(isreal(p2)); + } + else + { + size_t worstidx = 0; + T worstval = abs(cabs(filt.p[0]) - 1); + for (size_t i = 1; i < filt.p.size(); i++) + { + T val = abs(cabs(filt.p[i]) - 1); + if (val < worstval) + { + worstidx = i; + worstval = val; + } + } + p2_idx = worstidx; + p2 = filt.p[p2_idx]; + + TESTO_ASSERT(isreal(p2)); + z2_idx = internal::nearest_real_or_complex(filt.z, p2, true); + z2 = filt.z[z2_idx]; + TESTO_ASSERT(isreal(z2)); + filt.z.erase(filt.z.begin() + z2_idx); + } + filt.p.erase(filt.p.begin() + p2_idx); + } + } + pairs[si].p1 = p1; + pairs[si].p2 = p2; + pairs[si].z1 = z1; + pairs[si].z2 = z2; + } + + std::vector<biquad_params<T>> result(n_sections); + for (size_t si = 0; si < n_sections; si++) + { + result[si] = internal::zpk2tf(pairs[n_sections - 1 - si], si == 0 ? filt.k : T(1)); + } + return result; +} +} // namespace CMT_ARCH_NAME + +} // namespace kfr +\ No newline at end of file diff --git a/include/kfr/io/audiofile.hpp b/include/kfr/io/audiofile.hpp @@ -178,7 +178,7 @@ struct audio_writer_wav : audio_writer<T> (drwav_seek_proc)&internal_generic::drwav_writer_seek_proc, this->writer.get(), nullptr); } - ~audio_writer_wav() { close(); } + ~audio_writer_wav() override { close(); } using audio_writer<T>::write; @@ -243,7 +243,7 @@ struct audio_reader_wav : audio_reader<T> (drwav_seek_proc)&internal_generic::drwav_reader_seek_proc, this->reader.get(), nullptr); fmt.channels = f.channels; fmt.samplerate = f.sampleRate; - fmt.length = f.totalPCMFrameCount; + fmt.length = static_cast<imax>(f.totalPCMFrameCount); switch (f.translatedFormatTag) { case DR_WAVE_FORMAT_IEEE_FLOAT: @@ -288,7 +288,7 @@ struct audio_reader_wav : audio_reader<T> break; } } - ~audio_reader_wav() { drwav_uninit(&f); } + ~audio_reader_wav() override { drwav_uninit(&f); } /// @brief Returns audio format description const audio_format_and_length& format() const override { return fmt; } @@ -323,14 +323,13 @@ struct audio_reader_wav : audio_reader<T> switch (origin) { case seek_origin::current: - return drwav_seek_to_pcm_frame(&f, this->position + offset); + return drwav_seek_to_pcm_frame(&f, static_cast<drmp3_uint64>(this->position + offset)); case seek_origin::begin: - return drwav_seek_to_pcm_frame(&f, offset); + return drwav_seek_to_pcm_frame(&f, static_cast<drmp3_uint64>(offset)); case seek_origin::end: - return drwav_seek_to_pcm_frame(&f, fmt.length + offset); - default: - return false; + return drwav_seek_to_pcm_frame(&f, static_cast<drmp3_uint64>(fmt.length + offset)); } + return false; } private: @@ -355,10 +354,10 @@ struct audio_reader_flac : audio_reader<T> nullptr); fmt.channels = f->channels; fmt.samplerate = f->sampleRate; - fmt.length = f->totalPCMFrameCount; + fmt.length = static_cast<imax>(f->totalPCMFrameCount); fmt.type = audio_sample_type::i32; } - ~audio_reader_flac() { drflac_close(f); } + ~audio_reader_flac() override { drflac_close(f); } /// @brief Returns audio format description const audio_format_and_length& format() const override { return fmt; } @@ -394,14 +393,13 @@ struct audio_reader_flac : audio_reader<T> switch (origin) { case seek_origin::current: - return drflac_seek_to_pcm_frame(f, this->position + offset); + return drflac_seek_to_pcm_frame(f, static_cast<drmp3_uint64>(this->position + offset)); case seek_origin::begin: - return drflac_seek_to_pcm_frame(f, offset); + return drflac_seek_to_pcm_frame(f, static_cast<drmp3_uint64>(offset)); case seek_origin::end: - return drflac_seek_to_pcm_frame(f, fmt.length + offset); - default: - return false; + return drflac_seek_to_pcm_frame(f, static_cast<drmp3_uint64>(fmt.length + offset)); } + return false; } private: @@ -426,10 +424,10 @@ struct audio_reader_mp3 : audio_reader<T> nullptr); fmt.channels = f.channels; fmt.samplerate = f.sampleRate; - fmt.length = drmp3_get_pcm_frame_count(&f); + fmt.length = static_cast<imax>(drmp3_get_pcm_frame_count(&f)); fmt.type = audio_sample_type::i16; } - ~audio_reader_mp3() { drmp3_uninit(&f); } + ~audio_reader_mp3() override { drmp3_uninit(&f); } drmp3_config config{ 0, 0 }; @@ -467,14 +465,13 @@ struct audio_reader_mp3 : audio_reader<T> switch (origin) { case seek_origin::current: - return drmp3_seek_to_pcm_frame(&f, this->position + offset); + return drmp3_seek_to_pcm_frame(&f, static_cast<drmp3_uint64>(this->position + offset)); case seek_origin::begin: - return drmp3_seek_to_pcm_frame(&f, offset); + return drmp3_seek_to_pcm_frame(&f, static_cast<drmp3_uint64>(offset)); case seek_origin::end: - return drmp3_seek_to_pcm_frame(&f, fmt.length + offset); - default: - return false; + return drmp3_seek_to_pcm_frame(&f, static_cast<drmp3_uint64>(fmt.length + offset)); } + return false; } private: diff --git a/include/kfr/io/dr/dr_flac.h b/include/kfr/io/dr/dr_flac.h @@ -2264,20 +2264,20 @@ static DRFLAC_INLINE drflac_uint32 drflac__clz_lzcnt(drflac_cache_t x) #endif #else #if defined(__GNUC__) || defined(__clang__) - #if defined(DRFLAC_X64) + #if 0 /* defined(DRFLAC_X64) --- workaround for Clang bug */ { drflac_uint64 r; __asm__ __volatile__ ( - "lzcnt{ %1, %0| %0, %1}" : "=r"(r) : "r"(x) + "lzcnt{ %1, %0| %0, %1}" : "=r"(r) : "r"(x) : "cc" ); return (drflac_uint32)r; } - #elif defined(DRFLAC_X86) + #elif 0 /* defined(DRFLAC_X86) --- workaround for Clang bug */ { drflac_uint32 r; __asm__ __volatile__ ( - "lzcnt{l %1, %0| %0, %1}" : "=r"(r) : "r"(x) + "lzcnt{l %1, %0| %0, %1}" : "=r"(r) : "r"(x) : "cc" ); return r; diff --git a/include/kfr/io/python_plot.hpp b/include/kfr/io/python_plot.hpp @@ -41,9 +41,8 @@ namespace kfr namespace internal_generic { CMT_PRAGMA_GNU(GCC diagnostic push) -#if CMT_HAS_WARNING("-Wdeprecated-declarations") CMT_PRAGMA_GNU(GCC diagnostic ignored "-Wdeprecated-declarations") -#endif +CMT_PRAGMA_GNU(GCC diagnostic ignored "-Wunused-result") template <int = 0> void python(const std::string& name, const std::string& code) diff --git a/include/kfr/io/tostring.hpp b/include/kfr/io/tostring.hpp @@ -62,6 +62,10 @@ struct representation<cometa::special_value> return as_string(value.ll); case special_constant::floating_point: return as_string(value.d); + case special_constant::epsilon: + return "epsilon"; + case special_constant::random_bits: + return "random_bits"; default: return "unknown"; } @@ -152,6 +156,17 @@ struct representation<kfr::complex<T>> } }; +template <char t, int width, int prec, typename T> +struct representation<details::fmt_t<kfr::complex<T>, t, width, prec>> +{ + using type = std::string; + static std::string get(const details::fmt_t<kfr::complex<T>, t, width, prec>& value) + { + return as_string(fmt<t, width, prec>(value.value.real())) + " + " + + as_string(fmt<t, width, prec>(value.value.imag())) + "j"; + } +}; + template <> struct representation<kfr::cpu_t> { diff --git a/include/kfr/math/complex_math.hpp b/include/kfr/math/complex_math.hpp @@ -122,15 +122,24 @@ KFR_INTRINSIC vec<complex<T>, N> cartesian(const vec<complex<T>, N>& x) template <typename T, size_t N> KFR_INTRINSIC vec<T, N> cabsdup(const vec<T, N>& x) { - x = sqr(x); - return sqrt(x + swap<2>(x)); + vec<T, N> xx = sqr(x); + return sqrt(xx + swap<2>(xx)); } template <typename T, size_t N> KFR_INTRINSIC vec<complex<T>, N> csqrt(const vec<complex<T>, N>& x) { - const vec<T, N> t = (cabsdup(cdecom(x)) + cdecom(cnegimag(cdupreal(x)))) * T(0.5); - return ccomp(select(dupodd(x) < T(), cdecom(cnegimag(ccomp(t))), t)); + const vec<T, N> s = sqrt((abs(real(x)) + cabs(x)) * 0.5); + const vec<T, N> d = abs(imag(x)) * 0.5 / s; + const mask<T, N> posreal = real(x) >= T(0); + const vec<T, N> imagsign = imag(x) & special_constants<T>::highbitmask(); + return make_complex(select(posreal, s, d), select(posreal, d ^ imagsign, s ^ imagsign)); +} + +template <typename T, size_t N> +KFR_INTRINSIC vec<complex<T>, N> csqr(const vec<complex<T>, N>& x) +{ + return x * x; } KFR_HANDLE_SCALAR(cconj) @@ -147,6 +156,7 @@ KFR_HANDLE_SCALAR(cexp10) KFR_HANDLE_SCALAR(polar) KFR_HANDLE_SCALAR(cartesian) KFR_HANDLE_SCALAR(csqrt) +KFR_HANDLE_SCALAR(csqr) template <typename T, size_t N> KFR_INTRINSIC vec<T, N> cabs(const vec<T, N>& a) @@ -188,6 +198,7 @@ KFR_I_FN(cexp10) KFR_I_FN(polar) KFR_I_FN(cartesian) KFR_I_FN(csqrt) +KFR_I_FN(csqr) /// @brief Returns the sine of the complex number x template <typename T1, KFR_ENABLE_IF(is_numeric<T1>)> @@ -273,13 +284,6 @@ KFR_FUNCTION internal::expression_function<fn::carg, E1> carg(E1&& x) return { fn::carg(), std::forward<E1>(x) }; } -/// @brief Returns the complex conjugate of the complex number x -template <typename T1, KFR_ENABLE_IF(is_numeric<T1>)> -KFR_FUNCTION T1 cconj(const T1& x) -{ - return intrinsics::cconj(x); -} - /// @brief Returns template expression that returns the complex conjugate of the complex number x template <typename E1, KFR_ENABLE_IF(is_input_expression<E1>)> KFR_FUNCTION internal::expression_function<fn::cconj, E1> cconj(E1&& x) @@ -413,5 +417,19 @@ KFR_FUNCTION internal::expression_function<fn::csqrt, E1> csqrt(E1&& x) return { fn::csqrt(), std::forward<E1>(x) }; } +/// @brief Returns square of the complex number x +template <typename T1, KFR_ENABLE_IF(is_numeric<T1>)> +KFR_FUNCTION T1 csqr(const T1& x) +{ + return intrinsics::csqr(x); +} + +/// @brief Returns template expression that returns square of the complex number x +template <typename E1, KFR_ENABLE_IF(is_input_expression<E1>)> +KFR_FUNCTION internal::expression_function<fn::csqr, E1> csqr(E1&& x) +{ + return { fn::csqr(), std::forward<E1>(x) }; +} + } // namespace CMT_ARCH_NAME } // namespace kfr diff --git a/include/kfr/math/gamma.hpp b/include/kfr/math/gamma.hpp @@ -59,5 +59,37 @@ KFR_FUNCTION internal::expression_function<fn::factorial_approx, E1> factorial_a { return { fn::factorial_approx(), std::forward<E1>(x) }; } + +constexpr inline uint64_t factorial_table[21] = { + 0, + 1, + 2, + 6, + 24, + 120, + 720, + 5040, + 40320, + 362880, + 3628800, + 39916800, + 479001600, + 6227020800, + 87178291200, + 1307674368000, + 20922789888000, + 355687428096000, + 6402373705728000, + 121645100408832000, + 2432902008176640000, +}; + +/// @brief Returns the factorial of an argument. Returns max(uint64_t) if does not fit to uint64_t +constexpr uint64_t factorial(int n) +{ + if (n < 0 || n > 20) + return std::numeric_limits<uint64_t>::max(); + return factorial_table[n]; +} } // namespace CMT_ARCH_NAME } // namespace kfr diff --git a/include/kfr/simd/complex.hpp b/include/kfr/simd/complex.hpp @@ -307,6 +307,14 @@ KFR_INTRINSIC vec<complex<T>, N> cnegimag(const vec<complex<T>, N>& x) } KFR_FN(cnegimag) +/// @brief Returns mask with true for real elements +template <typename T> +KFR_INTRINSIC bool isreal(const complex<T>& x) +{ + return x.imag() == 0; +} +KFR_FN(isreal) + namespace internal { template <typename T> @@ -377,7 +385,7 @@ KFR_FN(real) template <typename E1, KFR_ENABLE_IF(is_input_expression<E1>)> KFR_INTRINSIC internal::expression_function<fn::real, E1> real(E1&& x) { - return { {}, std::forward<E1>(x) }; + return { fn::real{}, std::forward<E1>(x) }; } /// @brief Returns the imaginary part of the complex value @@ -399,7 +407,7 @@ KFR_FN(imag) template <typename E1, KFR_ENABLE_IF(is_input_expression<E1>)> KFR_INTRINSIC internal::expression_function<fn::imag, E1> imag(E1&& x) { - return { {}, std::forward<E1>(x) }; + return { fn::imag{}, std::forward<E1>(x) }; } /// @brief Constructs complex value from real and imaginary parts @@ -411,11 +419,20 @@ constexpr KFR_INTRINSIC vec<complex<T>, N> make_complex(const vec<T1, N>& real, } /// @brief Constructs complex value from real and imaginary parts -template <typename T1, typename T2 = T1, typename T = common_type<T1, T2>> +template <typename T1, typename T2 = T1, typename T = common_type<T1, T2>, + KFR_ENABLE_IF(is_numeric_args<T1, T2>)> constexpr KFR_INTRINSIC complex<T> make_complex(T1 real, T2 imag = T2(0)) { return complex<T>(innercast<T>(real), innercast<T>(imag)); } +KFR_FN(make_complex) + +/// @brief Constructs complex value from real and imaginary parts +template <typename E1, typename E2, KFR_ENABLE_IF(is_input_expressions<E1, E2>)> +KFR_INTRINSIC internal::expression_function<fn::make_complex, E1, E2> make_complex(E1&& re, E2&& im) +{ + return { fn::make_complex{}, std::forward<E1>(re), std::forward<E2>(im) }; +} namespace intrinsics { diff --git a/include/kfr/simd/impl/backend_generic.hpp b/include/kfr/simd/impl/backend_generic.hpp @@ -32,6 +32,41 @@ CMT_PRAGMA_GNU(GCC diagnostic ignored "-Wmaybe-uninitialized") namespace kfr { +template <size_t bits, size_t...> +struct shuffle_mask; + +template <size_t i0, size_t i1, size_t i2, size_t i3, size_t i4, size_t i5, size_t i6, size_t i7> +struct shuffle_mask<8, i0, i1, i2, i3, i4, i5, i6, i7> +{ + constexpr static inline size_t Nmax = 1; + constexpr static inline size_t value = (const_min(i7, Nmax) << 7) | (const_min(i6, Nmax) << 6) | + (const_min(i5, Nmax) << 5) | (const_min(i4, Nmax) << 4) | + (const_min(i3, Nmax) << 3) | (const_min(i2, Nmax) << 2) | + (const_min(i1, Nmax) << 1) | const_min(i0, Nmax); +}; + +template <size_t i0, size_t i1, size_t i2, size_t i3> +struct shuffle_mask<8, i0, i1, i2, i3> +{ + constexpr static inline size_t Nmax = 3; + constexpr static inline size_t value = (const_min(i3, Nmax) << 6) | (const_min(i2, Nmax) << 4) | + (const_min(i1, Nmax) << 2) | const_min(i0, Nmax); +}; + +template <size_t i0, size_t i1, size_t i2, size_t i3> +struct shuffle_mask<4, i0, i1, i2, i3> +{ + constexpr static inline size_t Nmax = 1; + constexpr static inline size_t value = (const_min(i3, Nmax) << 3) | (const_min(i2, Nmax) << 2) | + (const_min(i1, Nmax) << 1) | const_min(i0, Nmax); +}; + +template <size_t i0, size_t i1> +struct shuffle_mask<2, i0, i1> +{ + constexpr static inline size_t Nmax = 1; + constexpr static inline size_t value = (const_min(i1, Nmax) << 1) | const_min(i0, Nmax); +}; #if KFR_SHOW_NOT_OPTIMIZED CMT_PUBLIC_C CMT_DLL_EXPORT void not_optimized(const char* fn) CMT_NOEXCEPT; @@ -609,6 +644,22 @@ KFR_INTRIN_BITCAST(i32, f32, 8, _mm256_castps_si256(x)) KFR_INTRIN_BITCAST(f64, i64, 4, _mm256_castsi256_pd(x)) KFR_INTRIN_BITCAST(i64, f64, 4, _mm256_castpd_si256(x)) +KFR_INTRINSIC simd<float, 8> simd_shuffle(simd_t<float, 16>, const simd<float, 16>& x, + csizes_t<2, 3, 6, 7, 10, 11, 14, 15>, overload_priority<9>) +{ + const __m256 t1 = _mm256_permute2f128_ps(x.low, x.high, (0 << 0) | (2 << 4)); + const __m256 t2 = _mm256_permute2f128_ps(x.low, x.high, (1 << 0) | (3 << 4)); + return _mm256_shuffle_ps(t1, t2, shuffle_mask<8, 2, 3, 2, 3>::value); +} + +KFR_INTRINSIC simd<float, 8> simd_shuffle(simd_t<float, 16>, const simd<float, 16>& x, + csizes_t<0, 1, 4, 5, 8, 9, 12, 13>, overload_priority<9>) +{ + const __m256 t1 = _mm256_permute2f128_ps(x.low, x.high, (0 << 0) | (2 << 4)); + const __m256 t2 = _mm256_permute2f128_ps(x.low, x.high, (1 << 0) | (3 << 4)); + return _mm256_shuffle_ps(t1, t2, shuffle_mask<8, 0, 1, 0, 1>::value); +} + #ifndef CMT_ARCH_AVX2 KFR_INTRIN_SHUFFLE_DUPHALVES(i8, 16, _mm256_insertf128_si256(_mm256_castsi128_si256(x), x, 1)) KFR_INTRIN_SHUFFLE_DUPHALVES(u8, 16, _mm256_insertf128_si256(_mm256_castsi128_si256(x), x, 1)) @@ -1297,12 +1348,29 @@ KFR_INTRINSIC const simd<T, N1>& simd_concat(const simd<T, N1>& x) CMT_NOEXCEPT return x; } +template <typename T, size_t N1, size_t N2, size_t N3, size_t N4> +KFR_INTRINSIC simd<T, N1 + N2 + N3 + N4> simd_concat4(const simd<T, N1>& x, const simd<T, N2>& y, + const simd<T, N3>& z, const simd<T, N4>& w) CMT_NOEXCEPT +{ + return simd_shuffle(simd2_t<T, N1 + N2, N3 + N4>{}, + simd_shuffle(simd2_t<T, N1, N2>{}, x, y, csizeseq<N1 + N2>, overload_auto), + simd_shuffle(simd2_t<T, N3, N4>{}, z, w, csizeseq<N3 + N4>, overload_auto), + csizeseq<N1 + N2 + N3 + N4>, overload_auto); +} + template <typename T, size_t N1, size_t N2, size_t... Ns, size_t Nscount /*= csum(csizes<Ns...>)*/> KFR_INTRINSIC simd<T, N1 + N2 + Nscount> simd_concat(const simd<T, N1>& x, const simd<T, N2>& y, const simd<T, Ns>&... z) CMT_NOEXCEPT { - return simd_shuffle(simd2_t<T, N1, N2 + Nscount>{}, x, simd_concat<T, N2, Ns...>(y, z...), - csizeseq<N1 + N2 + Nscount>, overload_auto); + if constexpr (sizeof...(Ns) == 2) + { + return simd_concat4<T, N1, N2, Ns...>(x, y, z...); + } + else + { + return simd_shuffle(simd2_t<T, N1, N2 + Nscount>{}, x, simd_concat<T, N2, Ns...>(y, z...), + csizeseq<N1 + N2 + Nscount>, overload_auto); + } } template <typename Tout, typename Tin, size_t N, size_t... indices> @@ -1444,42 +1512,6 @@ SIMD_TYPE_INTRIN_EX(f64, 8, _mm512_cvtsd_f64(x), _mm512_castpd128_pd512(_mm_set_ KFR_mm512_setr_m256d(x, y)) #endif -template <size_t bits, size_t...> -struct shuffle_mask; - -template <size_t i0, size_t i1, size_t i2, size_t i3, size_t i4, size_t i5, size_t i6, size_t i7> -struct shuffle_mask<8, i0, i1, i2, i3, i4, i5, i6, i7> -{ - constexpr static inline size_t Nmax = 1; - constexpr static inline size_t value = (const_min(i7, Nmax) << 7) | (const_min(i6, Nmax) << 6) | - (const_min(i5, Nmax) << 5) | (const_min(i4, Nmax) << 4) | - (const_min(i3, Nmax) << 3) | (const_min(i2, Nmax) << 2) | - (const_min(i1, Nmax) << 1) | const_min(i0, Nmax); -}; - -template <size_t i0, size_t i1, size_t i2, size_t i3> -struct shuffle_mask<8, i0, i1, i2, i3> -{ - constexpr static inline size_t Nmax = 3; - constexpr static inline size_t value = (const_min(i3, Nmax) << 6) | (const_min(i2, Nmax) << 4) | - (const_min(i1, Nmax) << 2) | const_min(i0, Nmax); -}; - -template <size_t i0, size_t i1, size_t i2, size_t i3> -struct shuffle_mask<4, i0, i1, i2, i3> -{ - constexpr static inline size_t Nmax = 1; - constexpr static inline size_t value = (const_min(i3, Nmax) << 3) | (const_min(i2, Nmax) << 2) | - (const_min(i1, Nmax) << 1) | const_min(i0, Nmax); -}; - -template <size_t i0, size_t i1> -struct shuffle_mask<2, i0, i1> -{ - constexpr static inline size_t Nmax = 1; - constexpr static inline size_t value = (const_min(i1, Nmax) << 1) | const_min(i0, Nmax); -}; - #ifdef CMT_ARCH_SSE2 template <size_t I0, size_t I1, size_t I2, size_t I3> diff --git a/include/kfr/simd/impl/function.hpp b/include/kfr/simd/impl/function.hpp @@ -285,7 +285,7 @@ KFR_INTRINSIC void intrin(vec<T, N>& result, const T& a, const vec<T, N>& b, Fn& KFR_INTRINSIC vec<T, N> fn(const vec<T, N>& a, const vec<T, N>& b) CMT_NOEXCEPT \ { \ vec<T, N> r; \ - intrin(r, a, b, [](const auto& a, const auto& b) { return intrinsics::fn(a, b); }); \ + intrin(r, a, b, [](const auto& aa, const auto& bb) { return intrinsics::fn(aa, bb); }); \ return r; \ } \ template <typename T, size_t N, \ @@ -299,7 +299,7 @@ KFR_INTRINSIC void intrin(vec<T, N>& result, const T& a, const vec<T, N>& b, Fn& KFR_INTRINSIC vec<T, N> fn(const vec<T, N>& a, const T& b) CMT_NOEXCEPT \ { \ vec<T, N> r; \ - intrin(r, a, b, [](const auto& a, const auto& b) { return intrinsics::fn(a, b); }); \ + intrin(r, a, b, [](const auto& aa, const auto& bb) { return intrinsics::fn(aa, bb); }); \ return r; \ } \ template <typename T, size_t N, \ @@ -313,7 +313,7 @@ KFR_INTRINSIC void intrin(vec<T, N>& result, const T& a, const vec<T, N>& b, Fn& KFR_INTRINSIC vec<T, N> fn(const T& a, const vec<T, N>& b) CMT_NOEXCEPT \ { \ vec<T, N> r; \ - intrin(r, a, b, [](const auto& a, const auto& b) { return intrinsics::fn(a, b); }); \ + intrin(r, a, b, [](const auto& aa, const auto& bb) { return intrinsics::fn(aa, bb); }); \ return r; \ } diff --git a/include/kfr/simd/impl/simd.hpp b/include/kfr/simd/impl/simd.hpp @@ -95,7 +95,7 @@ struct simd_halves subtype low; subtype high; -#if KFR_DEFINE_CTORS_FOR_HALVES +#if defined KFR_DEFINE_CTORS_FOR_HALVES && KFR_DEFINE_CTORS_FOR_HALVES simd_halves() CMT_NOEXCEPT {} simd_halves(const subtype& l, const subtype& h) CMT_NOEXCEPT : low(l), high(h) {} simd_halves(const simd_halves& v) CMT_NOEXCEPT : low(v.low), high(v.high) {} diff --git a/include/kfr/simd/impl/specialconstants.hpp b/include/kfr/simd/impl/specialconstants.hpp @@ -77,7 +77,7 @@ template <> struct special_scalar_constants<float> { constexpr static float highbitmask() { return -0.f; } - KFR_CONSTEXPR_NON_INTEL static float allones() noexcept { return allones_f32(); }; + KFR_CONSTEXPR_NON_INTEL static float allones() noexcept { return allones_f32(); } constexpr static float allzeros() { return 0.f; } KFR_CONSTEXPR_NON_INTEL static float invhighbitmask() { return invhighbit_f32(); } }; @@ -86,7 +86,7 @@ template <> struct special_scalar_constants<double> { constexpr static double highbitmask() { return -0.; } - KFR_CONSTEXPR_NON_INTEL static double allones() noexcept { return allones_f64(); }; + KFR_CONSTEXPR_NON_INTEL static double allones() noexcept { return allones_f64(); } constexpr static double allzeros() { return 0.; } KFR_CONSTEXPR_NON_INTEL static double invhighbitmask() { return invhighbit_f64(); } }; diff --git a/include/kfr/simd/shuffle.hpp b/include/kfr/simd/shuffle.hpp @@ -501,12 +501,13 @@ KFR_INTRINSIC vec<T, N> splitpairs(const vec<T, N>& x) } KFR_FN(splitpairs) -template <size_t group = 1, typename T, size_t N, size_t size = N / group> +template <size_t group = 1, typename T, size_t N> KFR_INTRINSIC vec<T, N> reverse(const vec<T, N>& x) { + constexpr size_t size = N / group; return x.shuffle(scale<group>(csizeseq_t<size, size - 1, -1>())); } -template <typename T, size_t N1, size_t N2> +template <size_t group = 1, typename T, size_t N1, size_t N2> KFR_INTRINSIC vec<vec<T, N1>, N2> reverse(const vec<vec<T, N1>, N2>& x) { return swap<N1>(x.flatten()).v; diff --git a/include/kfr/simd/types.hpp b/include/kfr/simd/types.hpp @@ -309,7 +309,7 @@ template <typename T> struct special_scalar_constants<bit<T>> { constexpr static bit<T> highbitmask() { return true; } - constexpr static bit<T> allones() noexcept { return true; }; + constexpr static bit<T> allones() noexcept { return true; } constexpr static bit<T> allzeros() { return false; } constexpr static bit<T> invhighbitmask() { return false; } }; diff --git a/include/kfr/simd/vec.hpp b/include/kfr/simd/vec.hpp @@ -187,11 +187,8 @@ struct alignas(force_compiletime_size_t< using ST = typename compound_type_traits<T>::deep_subtype; using scalar_type = ST; - enum : size_t - { - SW = compound_type_traits<T>::deep_width, - SN = N * SW - }; + constexpr static size_t SW = compound_type_traits<T>::deep_width; + constexpr static size_t SN = N * SW; constexpr static size_t scalar_size() CMT_NOEXCEPT { return SN; } @@ -334,7 +331,11 @@ struct alignas(force_compiletime_size_t< KFR_ENABLE_IF(dummy == 0 && !compound_type_traits<T>::is_scalar)> KFR_MEM_INTRINSIC constexpr value_type get(size_t index) const CMT_NOEXCEPT { - return this->s[index]; + union { + simd_type v; + T s[N]; + } u{ this->v }; + return u.s[index]; } template <size_t index, KFR_ENABLE_IF(index < 1024 && compound_type_traits<T>::is_scalar)> @@ -366,7 +367,12 @@ struct alignas(force_compiletime_size_t< template <int dummy = 0, KFR_ENABLE_IF(dummy == 0 && !compound_type_traits<T>::is_scalar)> KFR_MEM_INTRINSIC constexpr void set(size_t index, const value_type& s) CMT_NOEXCEPT { - this->s[index] = s; + union { + simd_type v; + T s[N]; + } u{ this->v }; + u.s[index] = s; + this->v = u.v; } template <size_t index, KFR_ENABLE_IF(index < 1024 && compound_type_traits<T>::is_scalar)> @@ -437,8 +443,8 @@ public: union { simd_type v; vec_halves<T, N> h; - simd_element_type w[simd_element_count]; - T s[N]; + // simd_element_type w[simd_element_count]; + // T s[N]; }; }; @@ -1433,7 +1439,7 @@ namespace std { template <typename T, size_t N> -class tuple_size<kfr::vec<T, N>> : public integral_constant<size_t, N> +struct tuple_size<kfr::vec<T, N>> : public integral_constant<size_t, N> { }; @@ -1443,4 +1449,4 @@ struct tuple_element<I, kfr::vec<T, N>> using type = T; }; -} // namespace std -\ No newline at end of file +} // namespace std diff --git a/include/kfr/testo/comparison.hpp b/include/kfr/testo/comparison.hpp @@ -73,7 +73,9 @@ struct eplison_scope template <> struct eplison_scope<void> { - eplison_scope(float scale) : f(scale), d(scale), ld(scale) {} + eplison_scope(float scale) : f(scale), d(static_cast<double>(scale)), ld(static_cast<long double>(scale)) + { + } eplison_scope<float> f; eplison_scope<double> d; eplison_scope<long double> ld; diff --git a/include/kfr/testo/console_colors.hpp b/include/kfr/testo/console_colors.hpp @@ -15,9 +15,9 @@ namespace win32_lite typedef void* HANDLE; typedef uint32_t DWORD; -#define WIN32_LITE_STD_INPUT_HANDLE ((win32_lite::DWORD)-10) -#define WIN32_LITE_STD_OUTPUT_HANDLE ((win32_lite::DWORD)-11) -#define WIN32_LITE_STD_ERROR_HANDLE ((win32_lite::DWORD)-12) +#define WIN32_LITE_STD_INPUT_HANDLE (static_cast<win32_lite::DWORD>(-10)) +#define WIN32_LITE_STD_OUTPUT_HANDLE (static_cast<win32_lite::DWORD>(-11)) +#define WIN32_LITE_STD_ERROR_HANDLE (static_cast<win32_lite::DWORD>(-12)) #define WIN32_LITE_ENABLE_VIRTUAL_TERMINAL_PROCESSING (4) diff --git a/sources.cmake b/sources.cmake @@ -64,6 +64,7 @@ set( ${PROJECT_SOURCE_DIR}/include/kfr/dsp/fir_design.hpp ${PROJECT_SOURCE_DIR}/include/kfr/dsp/fracdelay.hpp ${PROJECT_SOURCE_DIR}/include/kfr/dsp/goertzel.hpp + ${PROJECT_SOURCE_DIR}/include/kfr/dsp/iir_design.hpp ${PROJECT_SOURCE_DIR}/include/kfr/dsp/mixdown.hpp ${PROJECT_SOURCE_DIR}/include/kfr/dsp/oscillators.hpp ${PROJECT_SOURCE_DIR}/include/kfr/dsp/sample_rate_conversion.hpp diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt @@ -24,6 +24,10 @@ if (CMAKE_CXX_COMPILER_ID MATCHES "Clang") add_compile_options(-fdiagnostics-absolute-paths) endif () +if (ENABLE_DFT) + add_definitions(-DHAVE_DFT) +endif () + if (MSVC) link_libraries(-DEBUG) else () diff --git a/tests/asm_test.cpp b/tests/asm_test.cpp @@ -7,7 +7,8 @@ #define KFR_EXTENDED_TESTS #include <kfr/base.hpp> -#ifdef KFR_ENABLE_DFT +#include <kfr/dsp.hpp> +#ifdef HAVE_DFT #include <kfr/dft/impl/fft-impl.hpp> #endif #include <kfr/io.hpp> @@ -245,7 +246,7 @@ TEST_ASM_UIF(read, TEST_READ) TEST_ASM_UIF(write, TEST_WRITE) -#ifdef KFR_ENABLE_DFT +#ifdef HAVE_DFT #define TEST_FFT_SPEC(ty, size) \ static intrinsics::fft_specialization<ty, size> fft__##ty##__##size(static_cast<size_t>(1 << size)); \ diff --git a/tests/complex_test.cpp b/tests/complex_test.cpp @@ -154,6 +154,19 @@ TEST(complex_basic_expressions) CHECK(uv3[14] == c32{ 14, 0 }); } +TEST(complex_functions) +{ + CHECK(csqr(complex<f32>(4.f, 0.f)) == c32{ 16.f, 0.f }); + CHECK(csqrt(complex<f32>(16.f, 0.f)) == c32{ 4.f, 0.f }); + + CHECK(csqr(complex<f32>(1.f, 4.f)) == c32{ -15.f, 8.f }); + + CHECK(csqrt(complex<f32>(15.f, 8.f)) == c32{ 4.f, 1.f }); + CHECK(csqrt(complex<f32>(-15.f, 8.f)) == c32{ 1.f, 4.f }); + CHECK(csqrt(complex<f32>(15.f, -8.f)) == c32{ 4.f, -1.f }); + CHECK(csqrt(complex<f32>(-15.f, -8.f)) == c32{ 1.f, -4.f }); +} + TEST(complex_function_expressions) { const univector<c32, 4> uv1 = sqr(counter()); diff --git a/tests/dft_test.cpp b/tests/dft_test.cpp @@ -143,6 +143,35 @@ TEST(fft_accuracy) } }); } + +TEST(dct) +{ + constexpr size_t size = 16; + dct_plan<float> plan(size); + univector<float, size> in = counter(); + univector<float, size> out; + univector<float, size> outinv; + univector<u8> tmp(plan.temp_size); + plan.execute(out, in, tmp, false); + + univector<float, size> refout = { 120., -51.79283109806667, 0., -5.6781471211595695, + 0., -1.9843883778092053, 0., -0.9603691873838152, + 0., -0.5308329190495176, 0., -0.3030379000702155, + 0., -0.1584982220313824, 0., -0.0494839805703826 }; + + CHECK(rms(refout - out) < 0.00001f); + + plan.execute(outinv, in, tmp, true); + + univector<float, size> refoutinv = { 59.00747544192212, -65.54341437693878, 27.70332758523579, + -24.56124678824279, 15.546989102481612, -14.293082621965974, + 10.08224348063459, -9.38097406470581, 6.795411054455922, + -6.320715753372687, 4.455202292297903, -4.0896421269390455, + 2.580439536964837, -2.2695816108369176, 0.9311870090070382, + -0.643618159997807 }; + + CHECK(rms(refoutinv - outinv) < 0.00001f); +} } // namespace CMT_ARCH_NAME #ifndef KFR_NO_MAIN diff --git a/tests/dsp_test.cpp b/tests/dsp_test.cpp @@ -275,7 +275,6 @@ TEST(mixdown) [](size_t i) { return i + i * 2 + 100; }); } -#ifdef CMT_COMPILER_CLANG__ TEST(mixdown_stereo) { const univector<double, 21> left = counter(); @@ -287,7 +286,6 @@ TEST(mixdown_stereo) CHECK_EXPRESSION(mid, 21, [](size_t i) { return i + i * 2.0 + 100.0; }); CHECK_EXPRESSION(side, 21, [](size_t i) { return i - (i * 2.0 + 100.0); }); } -#endif TEST(phasor) { diff --git a/tests/expression_test.cpp b/tests/expression_test.cpp @@ -18,6 +18,7 @@ namespace CMT_ARCH_NAME TEST(pack) { + static_assert(std::is_same_v<vec<f32x2, 1>, std::invoke_result_t<fn::reverse, vec<f32x2, 1>>>); const univector<float, 21> v1 = 1 + counter(); const univector<float, 21> v2 = v1 * 11;