kfr

Fast, modern C++ DSP framework, FFT, Sample Rate Conversion, FIR/IIR/Biquad Filters (SSE, AVX, AVX-512, ARM NEON)
Log | Files | Refs | README

commit 9d8fcebe092f6eded923b43d565ae36c10cf198f
parent 28c6ce8d1f98c198f5403f7139dee5597f8e32da
Author: d.levin256@gmail.com <d.levin256@gmail.com>
Date:   Fri,  7 Dec 2018 05:11:32 +0000

Sample Rate Convertor refactoring and fixes

Diffstat:
Mexamples/CMakeLists.txt | 3+++
Mexamples/sample_rate_conversion.cpp | 61++++++++++++++++++++++++-------------------------------------
Aexamples/sample_rate_converter.cpp | 71+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Minclude/kfr/base/conversion.hpp | 36++++++++++++++++++++++++++++++++++++
Minclude/kfr/base/univector.hpp | 12+++++++-----
Minclude/kfr/cometa.hpp | 67+++++++++++++++++++++++++++++++++++++++++++++----------------------
Minclude/kfr/cometa/array.hpp | 58++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Minclude/kfr/dsp/sample_rate_conversion.hpp | 157++++++++++++++++++++++++++++++++++++++++---------------------------------------
Minclude/kfr/io/audiofile.hpp | 10+++++++++-
Minclude/kfr/io/file.hpp | 20++++++++++++++++++++
10 files changed, 353 insertions(+), 142 deletions(-)

diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt @@ -33,5 +33,8 @@ target_link_libraries(fir kfr) add_executable(sample_rate_conversion sample_rate_conversion.cpp) target_link_libraries(sample_rate_conversion kfr kfr_io) +add_executable(sample_rate_converter sample_rate_converter.cpp) +target_link_libraries(sample_rate_converter kfr kfr_io) + add_executable(dft dft.cpp) target_link_libraries(dft kfr kfr_dft) diff --git a/examples/sample_rate_conversion.cpp b/examples/sample_rate_conversion.cpp @@ -13,7 +13,6 @@ using namespace kfr; constexpr size_t input_sr = 96000; constexpr size_t output_sr = 44100; constexpr size_t len = 96000 * 6; -constexpr fbase i32max = 2147483647.0; int main() { @@ -25,64 +24,52 @@ int main() { auto r = resampler<fbase>(resample_quality::high, output_sr, input_sr, 1.0, 0.496); - univector<fbase> resampled(len * output_sr / input_sr); + univector<fbase> resampled(len * output_sr / input_sr + r.get_delay()); + r.process(resampled, swept_sine); - const size_t destsize = r(resampled.data(), swept_sine); - - univector<i32> i32data = clamp(resampled.truncate(destsize) * i32max, -i32max, +i32max); - { - audio_writer_wav<i32> writer(open_file_for_writing(KFR_FILEPATH("audio_high_quality.wav")), - audio_format{ 2, audio_sample_type::i32, output_sr }); - writer.write(i32data.data(), i32data.size()); - } + audio_writer_wav<fbase> writer(open_file_for_writing(KFR_FILEPATH("audio_high_quality.wav")), + audio_format{ 1, audio_sample_type::i32, output_sr }); + writer.write(resampled.data(), resampled.size()); + writer.close(); plot_save("audio_high_quality", "audio_high_quality.wav", ""); } { auto r = resampler<fbase>(resample_quality::normal, output_sr, input_sr, 1.0, 0.496); - univector<fbase> resampled(len * output_sr / input_sr); - - const size_t destsize = r(resampled.data(), swept_sine); + univector<fbase> resampled(len * output_sr / input_sr + r.get_delay()); + r.process(resampled, swept_sine); - univector<i32> i32data = clamp(resampled.truncate(destsize) * i32max, -i32max, +i32max); - { - audio_writer_wav<i32> writer(open_file_for_writing(KFR_FILEPATH("audio_normal_quality.wav")), - audio_format{ 2, audio_sample_type::i32, output_sr }); - writer.write(i32data.data(), i32data.size()); - } + audio_writer_wav<fbase> writer(open_file_for_writing(KFR_FILEPATH("audio_normal_quality.wav")), + audio_format{ 1, audio_sample_type::i32, output_sr }); + writer.write(resampled.data(), resampled.size()); + writer.close(); plot_save("audio_normal_quality", "audio_normal_quality.wav", ""); } { auto r = resampler<fbase>(resample_quality::low, output_sr, input_sr, 1.0, 0.496); - univector<fbase> resampled(len * output_sr / input_sr); + univector<fbase> resampled(len * output_sr / input_sr + r.get_delay()); + r.process(resampled, swept_sine); - const size_t destsize = r(resampled.data(), swept_sine); - - univector<i32> i32data = clamp(resampled.truncate(destsize) * i32max, -i32max, +i32max); - { - audio_writer_wav<i32> writer(open_file_for_writing(KFR_FILEPATH("audio_low_quality.wav")), - audio_format{ 2, audio_sample_type::i32, output_sr }); - writer.write(i32data.data(), i32data.size()); - } + audio_writer_wav<fbase> writer(open_file_for_writing(KFR_FILEPATH("audio_low_quality.wav")), + audio_format{ 1, audio_sample_type::i32, output_sr }); + writer.write(resampled.data(), resampled.size()); + writer.close(); plot_save("audio_low_quality", "audio_low_quality.wav", ""); } { auto r = resampler<fbase>(resample_quality::draft, output_sr, input_sr, 1.0, 0.496); - univector<fbase> resampled(len * output_sr / input_sr); - - const size_t destsize = r(resampled.data(), swept_sine); + univector<fbase> resampled(len * output_sr / input_sr + r.get_delay()); + r.process(resampled, swept_sine); - univector<i32> i32data = clamp(resampled.truncate(destsize) * i32max, -i32max, +i32max); - { - audio_writer_wav<i32> writer(open_file_for_writing(KFR_FILEPATH("audio_draft_quality.wav")), - audio_format{ 2, audio_sample_type::i32, output_sr }); - writer.write(i32data.data(), i32data.size()); - } + audio_writer_wav<fbase> writer(open_file_for_writing(KFR_FILEPATH("audio_draft_quality.wav")), + audio_format{ 1, audio_sample_type::i32, output_sr }); + writer.write(resampled.data(), resampled.size()); + writer.close(); plot_save("audio_draft_quality", "audio_draft_quality.wav", ""); } diff --git a/examples/sample_rate_converter.cpp b/examples/sample_rate_converter.cpp @@ -0,0 +1,71 @@ +/** + * 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(int argc, char** argv) +{ + println(library_version()); + if (argc < 4) + { + println("Usage: sample_rate_converter <INPUT_FILE> <OUTPUT_FILE> <TARGET_SAMPLE_RATE>"); + println("Supported formats: WAV/W64, 16, 24, 32-bit PCM, 32, 64-bit IEEE"); + return 1; + } + + const size_t output_sr = std::atol(argv[3]); + + audio_reader_wav<double> reader(open_file_for_reading(argv[1])); + const size_t input_sr = reader.format().samplerate; + univector<double> input_interleaved(reader.format().length * reader.format().channels); + reader.read(input_interleaved.data(), input_interleaved.size()); + univector2d<double> input_channels( + reader.format().channels, univector<double>(input_interleaved.size() / reader.format().channels)); + deinterleave(input_channels, input_interleaved); + + univector2d<double> output_channels; + println("Input channels: ", reader.format().channels); + println("Input sample rate: ", reader.format().samplerate); + println("Input bit depth: ", audio_sample_bit_depth(reader.format().type)); + + for (size_t ch = 0; ch < input_channels.size(); ++ch) + { + println("Processing ", ch, " of ", reader.format().channels); + const univector<double>& input = input_channels[ch]; + auto r = resampler<double>(resample_quality::high, output_sr, input_sr, 1.0, 0.492); + const size_t output_size = input.size() * output_sr / input_sr; + univector<double> output(output_size); + const size_t input_step = r.skip(r.get_delay(), input.slice()); + + size_t input_pos = input_step; + size_t output_pos = 0; + for (;;) + { + const size_t block_size = std::min(size_t(4096), output.size() - output_pos); + if (block_size == 0) + break; + const size_t input_step = + r.process(output.slice(output_pos, block_size).ref(), input.slice(input_pos)); + input_pos += input_step; + output_pos += block_size; + } + + output_channels.push_back(std::move(output)); + } + + univector<double> output_interleved = interleave(output_channels); + + audio_writer_wav<double> writer( + open_file_for_writing(argv[2]), + audio_format{ reader.format().channels, reader.format().type, kfr::fmax(output_sr) }); + writer.write(output_interleved.data(), output_interleved.size()); + + return 0; +} diff --git a/include/kfr/base/conversion.hpp b/include/kfr/base/conversion.hpp @@ -190,6 +190,19 @@ void deinterleave(Tout* out[], const Tin* in, size_t channels, size_t size) } } +template <typename Tout, size_t Tag1, size_t Tag2, typename Tin, size_t Tag3> +void deinterleave(univector2d<Tout, Tag1, Tag2>& out, const univector<Tin, Tag3>& in) +{ + if (in.empty() || out.empty()) + return; + std::vector<Tout*> ptrs(out.size()); + for (size_t i = 0; i < out.size(); ++i) + { + ptrs[i] = out[i].data(); + } + return deinterleave(ptrs.data(), in.data(), out.size(), in.size() / out.size()); +} + template <typename Tout, typename Tin, typename Tout_traits = audio_sample_traits<Tout>, typename Tin_traits = audio_sample_traits<Tin>> void interleave(Tout* out, const Tin* in[], size_t channels, size_t size) @@ -201,6 +214,29 @@ void interleave(Tout* out, const Tin* in[], size_t channels, size_t size) } } +template <typename Tout, size_t Tag1, typename Tin, size_t Tag2, size_t Tag3> +void interleave(univector<Tout, Tag1>& out, const univector2d<Tin, Tag2, Tag3>& in) +{ + if (in.empty() || out.empty()) + return; + std::vector<const Tin*> ptrs(in.size()); + for (size_t i = 0; i < in.size(); ++i) + { + ptrs[i] = in[i].data(); + } + return interleave(out.data(), ptrs.data(), in.size(), out.size() / in.size()); +} + +template <typename Tin, size_t Tag1, size_t Tag2> +univector<Tin> interleave(const univector2d<Tin, Tag1, Tag2>& in) +{ + if (in.empty()) + return {}; + univector<Tin> result(in.size() * in[0].size()); + interleave(result, in); + return result; +} + template <typename Tout, typename Tin, typename Tout_traits = audio_sample_traits<Tout>, typename Tin_traits = audio_sample_traits<Tin>> void convert(Tout* out, const Tin* in, size_t size) diff --git a/include/kfr/base/univector.hpp b/include/kfr/base/univector.hpp @@ -92,13 +92,13 @@ struct univector_base : input_expression, output_expression { T* data = derived_cast<Class>(this)->data(); const size_t this_size = derived_cast<Class>(this)->size(); - return array_ref<T>(data + start, std::min(size, this_size - start)); + return array_ref<T>(data + start, std::min(size, start < this_size ? this_size - start : 0)); } univector<const T, 0> slice(size_t start = 0, size_t size = max_size_t) const { const T* data = derived_cast<Class>(this)->data(); const size_t this_size = derived_cast<Class>(this)->size(); - return array_ref<const T>(data + start, std::min(size, this_size - start)); + return array_ref<const T>(data + start, std::min(size, start < this_size ? this_size - start : 0)); } univector<T, 0> truncate(size_t size = max_size_t) { @@ -296,6 +296,8 @@ struct univector<T, tag_array_ref> : array_ref<T>, univector_base<T, univector<T return index < this->size() ? this->operator[](index) : fallback_value; } using univector_base<T, univector>::operator=; + + univector<T, tag_array_ref>& ref() && { return *this; } }; template <typename T> @@ -335,7 +337,7 @@ struct univector<T, tag_dynamic_vector> : std::vector<T, allocator<T>>, template <typename Input, KFR_ENABLE_IF(is_input_expression<Input>::value)> CMT_INLINE univector& operator=(Input&& input) { - if (this->empty()) + if (input.size() != infinite_size) this->resize(input.size()); this->assign_expr(std::forward<Input>(input)); return *this; @@ -353,7 +355,7 @@ using univector2d = univector<univector<T, Size2>, Size1>; template <typename T, size_t Size1 = tag_dynamic_vector, size_t Size2 = tag_dynamic_vector, size_t Size3 = tag_dynamic_vector> -using univector3d = univector<univector<univector<T, Size3>, Size2>, Size1>; +using univector3d = univector<univector<univector<T, Size3>, Size2>, Size1>; /// @brief Creates univector from data and size template <typename T> @@ -455,6 +457,6 @@ private: char cacheline_filler[64 - sizeof(std::atomic<size_t>)]; std::atomic<size_t> tail; }; -} +} // namespace kfr CMT_PRAGMA_MSVC(warning(pop)) diff --git a/include/kfr/cometa.hpp b/include/kfr/cometa.hpp @@ -19,8 +19,8 @@ CMT_PRAGMA_MSVC(warning(disable : 4814)) namespace cometa { -using std::size_t; using std::ptrdiff_t; +using std::size_t; #if __cplusplus >= 201103L || CMT_MSC_VER >= 1900 || CMT_HAS_FEATURE(cxx_constexpr) @@ -100,7 +100,7 @@ template <typename T, typename... Ts> struct and_t_impl<T, Ts...> : std::integral_constant<bool, T::value && and_t_impl<Ts...>::value> { }; -} +} // namespace details template <typename... Ts> using or_t = details::or_t_impl<Ts...>; @@ -253,7 +253,7 @@ namespace ops struct empty { }; -} +} // namespace ops template <typename T, T val> struct cval_t : ops::empty @@ -319,7 +319,7 @@ template <typename T, T val> struct is_val_impl<cval_t<T, val>> : std::true_type { }; -} +} // namespace details template <typename T> using is_inheritable = typename details::is_inheritable_impl<T>::type; @@ -376,7 +376,7 @@ template <size_t index> struct get_nth_type<index> { }; -} +} // namespace details template <typename T, T... values> struct cvals_t : ops::empty @@ -509,7 +509,7 @@ struct concat_impl<ctypes_t<types1...>, ctypes_t<types2...>> { using type = ctypes_t<types1..., types2...>; }; -} +} // namespace details template <typename T1, typename T2> using concat_lists = typename details::concat_impl<T1, T2>::type; @@ -561,7 +561,7 @@ struct filter_impl<cvals_t<T, value, values...>, cvals_t<bool, flag, flags...>> using filtered = typename filter_impl<cvals_t<T, values...>, cvals_t<bool, flags...>>::type; using type = conditional<flag, concat_lists<cvals_t<T, value>, filtered>, filtered>; }; -} +} // namespace details template <typename Fn> using function_arguments = typename details::function_arguments_impl<decltype(&Fn::operator())>::args; @@ -641,7 +641,7 @@ CMT_BIN_OP(&) CMT_BIN_OP(|) CMT_BIN_OP(^) // clang-format on -} +} // namespace ops namespace details { @@ -665,7 +665,7 @@ template <typename T, T Nstart, ptrdiff_t Nstep> struct cvalseq_impl<T, 1, Nstart, Nstep> : cvals_t<T, static_cast<T>(Nstart)> { }; -} +} // namespace details template <typename T, size_t size, T start = T(), ptrdiff_t step = 1> using cvalseq_t = typename details::cvalseq_impl<T, size, start, step>::type; @@ -738,7 +738,7 @@ struct unique_enum_impl #endif #define CMT_ENABLE_IF(...) CMT_ENABLE_IF_IMPL(__LINE__, __VA_ARGS__) -} +} // namespace details template <typename T> struct is_enabled : details::is_enabled_impl<T> @@ -768,7 +768,7 @@ inline auto call_if_callable(Fn&& fn) { return std::forward<Fn>(fn); } -} +} // namespace details template <typename Fn, typename... Args> inline auto bind_func(Fn&& fn, Args&&... args) @@ -854,6 +854,28 @@ constexpr inline T lcm(T a, T b, T c, Ts... rest) return lcm(a, lcm(b, c, rest...)); } +inline std::div_t floor_div(int a, int b) +{ + std::div_t d = std::div(a, b); + if (d.rem < 0) + { + d.rem += b; + --d.quot; + } + return d; +} + +inline std::lldiv_t floor_div(long long a, long long b) +{ + std::lldiv_t d = std::lldiv(a, b); + if (d.rem < 0) + { + d.rem += b; + --d.quot; + } + return d; +} + namespace details { @@ -939,7 +961,7 @@ template <typename T> using is_number_impl = std::integral_constant<bool, ((std::is_integral<T>::value) || (std::is_floating_point<T>::value)) && !std::is_same<T, bool>::value>; -} +} // namespace details template <size_t bits> using float_type = typename details::float_type_impl<bits>::type; @@ -977,7 +999,7 @@ constexpr size_t elementsize<void>() { return 1; }; -} +} // namespace details /// @brief Utility typedef used to disable type deduction template <typename T> @@ -1100,6 +1122,7 @@ struct carray : carray<T, N - 1> CMT_INTRIN constexpr const T* data() const noexcept { return begin(); } CMT_INTRIN constexpr T* data() noexcept { return begin(); } CMT_INTRIN constexpr bool empty() const noexcept { return false; } + private: T val; }; @@ -1258,7 +1281,7 @@ struct value_type_impl<T, Fallback, void_t<typename T::value_type>> { using type = typename T::value_type; }; -} +} // namespace details template <typename T> using has_begin_end = details::has_begin_end_impl<decay<T>>; @@ -1277,7 +1300,7 @@ void cforeach_impl(Fn&& fn) { fn(cval_t<T, value>()); } -} +} // namespace details #endif template <typename T, T... values, typename Fn> @@ -1322,7 +1345,7 @@ CMT_INTRIN void cforeach_types_impl(ctypes_t<T0, types...> type_list, Fn&& fn, c { swallow{ (fn(get_type_arg<indices>(type_list)), void(), 0)... }; } -} +} // namespace details template <typename... Ts, typename Fn> CMT_INTRIN void cforeach(ctypes_t<Ts...> types, Fn&& fn) @@ -1428,7 +1451,7 @@ inline decltype(auto) cmatch_impl(T&& value, Fn&& last) { return last(std::forward<T>(value)); } -} +} // namespace details template <typename T, typename Fn, typename... Args> inline decltype(auto) cmatch(T&& value, Fn&& fn, Args... args) @@ -1523,7 +1546,7 @@ struct signed_type_impl<T, void_t<enable_if<std::is_unsigned<T>::value>>> { using type = findinttype<std::numeric_limits<T>::min(), std::numeric_limits<T>::max()>; }; -} +} // namespace details template <typename T> using signed_type = typename details::signed_type_impl<T>::type; @@ -1631,7 +1654,7 @@ constexpr std::false_type test_sequence(...) { return {}; } -} +} // namespace details template <size_t number, size_t... numbers> constexpr bool is_sequence(csizes_t<number, numbers...>) @@ -1712,18 +1735,18 @@ constexpr indicesfor_t<List...> indicesfor{}; // Workaround for GCC 4.8 template <typename T> -constexpr conditional<std::is_scalar<T>::value, T, const T &> const_max(const T& x, const T& y) +constexpr conditional<std::is_scalar<T>::value, T, const T&> const_max(const T& x, const T& y) { return x > y ? x : y; } template <typename T> -constexpr conditional<std::is_scalar<T>::value, T, const T &> const_min(const T& x, const T& y) +constexpr conditional<std::is_scalar<T>::value, T, const T&> const_min(const T& x, const T& y) { return x < y ? x : y; } CMT_PRAGMA_GNU(GCC diagnostic pop) -} +} // namespace cometa CMT_PRAGMA_GNU(GCC diagnostic pop) diff --git a/include/kfr/cometa/array.hpp b/include/kfr/cometa/array.hpp @@ -126,4 +126,62 @@ inline array_ref<const T> make_array_ref(const std::vector<T>& cont) { return array_ref<const T>(cont.data(), cont.size()); } + +template <typename C> +constexpr auto datatype(C& c) +{ + return c[0]; +} +template <typename C> +constexpr auto datatype(const C& c) +{ + return c[0]; +} +template <typename E> +constexpr E datatype(const std::initializer_list<E>& il) +{ + return {}; +} +template <typename T, std::size_t N> +constexpr T datatype(T (&array)[N]) +{ + return {}; +} + +template <typename C> +constexpr auto data(C& c) -> decltype(c.data()) +{ + return c.data(); +} +template <typename C> +constexpr auto data(const C& c) -> decltype(c.data()) +{ + return c.data(); +} +template <typename T, std::size_t N> +constexpr T* data(T (&array)[N]) noexcept +{ + return array; +} +template <typename T> +constexpr T* data(T* array) noexcept +{ + return array; +} +template <typename E> +constexpr const E* data(const std::initializer_list<E>& il) noexcept +{ + return il.begin(); +} + +template <typename C> +constexpr auto size(const C& c) -> decltype(c.size()) +{ + return c.size(); +} +template <typename T, std::size_t N> +constexpr std::size_t size(const T (&array)[N]) noexcept +{ + return N; +} } // namespace cometa diff --git a/include/kfr/dsp/sample_rate_conversion.hpp b/include/kfr/dsp/sample_rate_conversion.hpp @@ -39,7 +39,7 @@ constexpr csize_t<4> draft{}; constexpr csize_t<6> low{}; constexpr csize_t<8> normal{}; constexpr csize_t<10> high{}; -} +} // namespace sample_rate_conversion_quality namespace resample_quality = sample_rate_conversion_quality; @@ -95,102 +95,105 @@ struct sample_rate_converter const T s = reciprocal(sum(filter)) * interpolation_factor; filter = filter * s; } - CMT_INLINE size_t operator()(T* dest, size_t zerosize) + + itype input_position_to_intermediate(itype in_pos) const { return in_pos * interpolation_factor; } + itype output_position_to_intermediate(itype out_pos) const { return out_pos * decimation_factor; } + + itype input_position_to_output(itype in_pos) const { - size_t outputsize = 0; - const itype srcsize = itype(zerosize); + return floor_div(input_position_to_intermediate(in_pos), decimation_factor).quot; + } + itype output_position_to_input(itype out_pos) const + { + return floor_div(output_position_to_intermediate(out_pos), interpolation_factor).quot; + } - for (size_t i = 0;; i++) - { - const itype ii = itype(i) + output_position; - const itype workindex = ii * (decimation_factor); - const itype workindex_rem = workindex % (interpolation_factor); - const itype start = workindex_rem ? (interpolation_factor)-workindex_rem : 0; - itype srcindex = workindex / (interpolation_factor); - srcindex = workindex_rem ? srcindex + 1 : srcindex; - const univector_ref<T> tap_ptr = filter.slice(static_cast<size_t>(start * depth)); - srcindex = srcindex - (depth - 1); - - if (srcindex + depth >= input_position + srcsize) - break; - outputsize++; - - if (dest) - { - if (srcindex >= input_position) - { - dest[i] = T(0); - } - else - { - const itype prev_count = input_position - srcindex; - dest[i] = dotproduct(delay.slice(size_t(depth - prev_count)), tap_ptr); - } - } - } - if (srcsize >= depth) + itype output_size_for_input(itype input_size) const + { + return input_position_to_output(input_position + input_size - 1) - + input_position_to_output(input_position - 1); + } + + itype input_size_for_output(itype output_size) const + { + return output_position_to_input(output_position + output_size - 1) - + output_position_to_input(output_position - 1); + } + + size_t skip(size_t output_size, univector_ref<const T> input) + { + const itype required_input_size = input_size_for_output(output_size); + + if (required_input_size >= depth) { - delay = zeros(); + delay.slice(0, delay.size()) = padded(input.slice(size_t(required_input_size - depth))); } else { - delay.truncate(size_t(depth - srcsize)) = delay.slice(size_t(srcsize)); - delay.slice(size_t(depth - srcsize)) = zeros(); + delay.truncate(size_t(depth - required_input_size)) = delay.slice(size_t(required_input_size)); + delay.slice(size_t(depth - required_input_size)) = padded(input); } - input_position += srcsize; - output_position += outputsize; - return outputsize; + input_position += required_input_size; + output_position += output_size; + + return required_input_size; } - CMT_INLINE size_t operator()(T* dest, univector_ref<const T> src) + + /// @brief Writes output.size() samples to output reading at most input.size(), then consuming zeros as + /// input. + /// @returns Number of processed input samples (may be less than input.size()). + template <size_t Tag> + size_t process(univector<T, Tag>& output, univector_ref<const T> input) { - size_t outputsize = 0; - const itype srcsize = itype(src.size()); + const itype required_input_size = input_size_for_output(output.size()); - for (size_t i = 0;; i++) + const itype input_size = input.size(); + for (size_t i = 0; i < output.size(); i++) { - const itype ii = itype(i) + output_position; - const itype workindex = ii * (decimation_factor); - const itype workindex_rem = workindex % (interpolation_factor); - const itype start = workindex_rem ? (interpolation_factor)-workindex_rem : 0; - itype srcindex = workindex / (interpolation_factor); - srcindex = workindex_rem ? srcindex + 1 : srcindex; - const univector_ref<T> tap_ptr = filter.slice(static_cast<size_t>(start * depth)); - srcindex = srcindex - (depth - 1); - - if (srcindex + depth >= input_position + srcsize) - break; - outputsize++; - - if (dest) + const itype intermediate_index = + output_position_to_intermediate(static_cast<itype>(i) + output_position); + const itype intermediate_start = intermediate_index - taps + 1; + const std::lldiv_t input_pos = + floor_div(intermediate_start + interpolation_factor - 1, interpolation_factor); + const itype input_start = input_pos.quot; // first input sample + const itype input_end = input_start + depth; + const itype tap_start = interpolation_factor - 1 - input_pos.rem; + const univector_ref<T> tap_ptr = filter.slice(static_cast<size_t>(tap_start * depth)); + + if (input_start >= input_position + input_size) + { + output[i] = T(0); + } + else if (input_start >= input_position) + { + output[i] = dotproduct(input.slice(input_start - input_position, depth), tap_ptr); + } + else { - if (srcindex >= input_position) - { - dest[i] = - dotproduct(src.slice(size_t(srcindex - input_position), size_t(depth)), tap_ptr); - } - else - { - const itype prev_count = input_position - srcindex; - dest[i] = dotproduct(delay.slice(size_t(depth - prev_count)), tap_ptr) + - dotproduct(src, tap_ptr.slice(size_t(prev_count), size_t(depth - prev_count))); - } + const itype prev_count = input_position - input_start; + output[i] = dotproduct(delay.slice(size_t(depth - prev_count)), tap_ptr) + + dotproduct(input.slice(0, size_t(depth - prev_count)), + tap_ptr.slice(size_t(prev_count), size_t(depth - prev_count))); } } - if (srcsize >= depth) + + if (required_input_size >= depth) { - delay = src.slice(size_t(srcsize - depth)); + delay.slice(0, delay.size()) = padded(input.slice(size_t(required_input_size - depth))); } else { - delay.truncate(size_t(depth - srcsize)) = delay.slice(size_t(srcsize)); - delay.slice(size_t(depth - srcsize)) = src; + delay.truncate(size_t(depth - required_input_size)) = delay.slice(size_t(required_input_size)); + delay.slice(size_t(depth - required_input_size)) = padded(input); } - input_position += srcsize; - output_position += outputsize; - return outputsize; + input_position += required_input_size; + output_position += output.size(); + + return required_input_size; } + size_t get_delay() const { return depth / 2 * interpolation_factor / decimation_factor; } itype taps; size_t order; itype interpolation_factor; @@ -301,7 +304,7 @@ struct expression_downsample<4, offset, E> : expression_base<E> return x.shuffle(csizeseq_t<N, offset, 4>()); } }; -} +} // namespace internal template <typename E1, size_t offset = 0> CMT_INLINE internal::expression_downsample<2, offset, E1> downsample2(E1&& e1, csize_t<offset> = csize_t<0>()) @@ -348,4 +351,4 @@ inline internal::sample_rate_converter<T, quality> resampler(csize_t<quality>, s return internal::sample_rate_converter<T, quality>(itype(interpolation_factor), itype(decimation_factor), scale, cutoff); } -} +} // namespace kfr diff --git a/include/kfr/io/audiofile.hpp b/include/kfr/io/audiofile.hpp @@ -125,10 +125,12 @@ struct audio_writer_wav : audio_writer<T> f = drwav_open_write(&wav_fmt, (drwav_write_proc)&internal::drwav_writer_write_proc, (drwav_seek_proc)&internal::drwav_writer_seek_proc, this->writer.get()); } - ~audio_writer_wav() { drwav_close(f); } + ~audio_writer_wav() { close(); } size_t write(const T* data, size_t size) { + if (!f) + return 0; if (fmt.type == audio_sample_type::unknown) return 0; if (fmt.type == audio_sample_traits<T>::type) @@ -146,6 +148,12 @@ struct audio_writer_wav : audio_writer<T> return sz; } } + void close() + { + drwav_close(f); + f = nullptr; + writer.reset(); + } const audio_format_and_length& format() const { return fmt; } imax tell() const { return fmt.length; } bool seek(imax position, seek_origin origin) { return false; } diff --git a/include/kfr/io/file.hpp b/include/kfr/io/file.hpp @@ -211,6 +211,26 @@ inline std::shared_ptr<file_writer<T>> open_file_for_appending(const filepath& p return std::make_shared<file_writer<T>>(fopen_portable(path.c_str(), KFR_FILEPATH("ab"))); } +#ifdef CMT_OS_WIN +template <typename T = void> +inline std::shared_ptr<file_reader<T>> open_file_for_reading(const std::string& path) +{ + return std::make_shared<file_reader<T>>(fopen(path.c_str(), "rb")); +} + +template <typename T = void> +inline std::shared_ptr<file_writer<T>> open_file_for_writing(const std::string& path) +{ + return std::make_shared<file_writer<T>>(fopen(path.c_str(), "wb")); +} + +template <typename T = void> +inline std::shared_ptr<file_writer<T>> open_file_for_appending(const std::string& path) +{ + return std::make_shared<file_writer<T>>(fopen(path.c_str(), "ab")); +} +#endif + namespace internal { struct expression_file_base