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 4708af551a1ccb5f58221f7f3bb50fb0650931c0
parent 46020a812cfcb5f8022a389f683c5f4f96f85c05
Author: samuriddle@gmail.com <samuriddle@gmail.com>
Date:   Sun,  7 Aug 2016 23:28:12 +0300

rename sample rate converter

Diffstat:
Mexamples/CMakeLists.txt | 2+-
Dexamples/resampling.cpp | 105-------------------------------------------------------------------------------
Aexamples/sample_rate_conversion.cpp | 105+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Minclude/kfr/dsp/resample.hpp | 188+------------------------------------------------------------------------------
Ainclude/kfr/dsp/sample_rate_conversion.hpp | 223+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Msources.cmake | 1+
6 files changed, 331 insertions(+), 293 deletions(-)

diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt @@ -34,5 +34,5 @@ include_directories(../include) add_executable(biquads biquads.cpp ${KFR_SRC}) add_executable(window window.cpp ${KFR_SRC}) add_executable(fir fir.cpp ${KFR_SRC}) -add_executable(resampling resampling.cpp ${KFR_SRC}) +add_executable(sample_rate_conversion sample_rate_conversion.cpp ${KFR_SRC}) add_executable(dft dft.cpp ${KFR_SRC} ${DFT_SRC}) diff --git a/examples/resampling.cpp b/examples/resampling.cpp @@ -1,105 +0,0 @@ -/** - * KFR (http://kfrlib.com) - * Copyright (C) 2016 D Levin - * See LICENSE.txt for details - */ - -// library_version() -#include <kfr/version.hpp> - -// print(), format() -#include <kfr/cometa/string.hpp> - -#include <kfr/math.hpp> - -// resample* -#include <kfr/dsp/resample.hpp> - -// file* -#include <kfr/io/audiofile.hpp> - -// swept -#include <kfr/dsp/oscillators.hpp> - -// plot_save() -#include <kfr/io/python_plot.hpp> - -#include <iostream> - -using namespace kfr; - -constexpr size_t input_sr = 96000; -constexpr size_t output_sr = 44100; -constexpr size_t len = 96000 * 6; -constexpr f64 i32max = 2147483647.0; - -int main(int argc, char** argv) -{ - println(library_version()); - - const std::string options = "phaseresp=False"; - - univector<f64> swept_sine = swept(0.5, len); - - { - auto r = resampler(resample_quality::high, output_sr, input_sr, 1.0, 0.496); - univector<f64> resampled(len * output_sr / input_sr); - - const size_t destsize = r(resampled.data(), swept_sine); - - univector<i32> i32data = clamp(resampled.slice(0, destsize) * i32max, -i32max, +i32max); - univector2d<i32> data = { i32data }; - - auto wr = sequential_file_writer("audio_high_quality.wav"); - audio_encode(wr, data, audioformat(data, output_sr)); - - plot_save("audio_high_quality", "audio_high_quality.wav", ""); - } - - { - auto r = resampler(resample_quality::normal, output_sr, input_sr, 1.0, 0.496); - univector<f64> resampled(len * output_sr / input_sr); - - const size_t destsize = r(resampled.data(), swept_sine); - - univector<i32> i32data = clamp(resampled.slice(0, destsize) * i32max, -i32max, +i32max); - univector2d<i32> data = { i32data }; - - auto wr = sequential_file_writer("audio_normal_quality.wav"); - audio_encode(wr, data, audioformat(data, output_sr)); - - plot_save("audio_normal_quality", "audio_normal_quality.wav", ""); - } - - { - auto r = resampler(resample_quality::low, output_sr, input_sr, 1.0, 0.496); - univector<f64> resampled(len * output_sr / input_sr); - - const size_t destsize = r(resampled.data(), swept_sine); - - univector<i32> i32data = clamp(resampled.slice(0, destsize) * i32max, -i32max, +i32max); - univector2d<i32> data = { i32data }; - - auto wr = sequential_file_writer("audio_low_quality.wav"); - audio_encode(wr, data, audioformat(data, output_sr)); - - plot_save("audio_low_quality", "audio_low_quality.wav", ""); - } - - { - auto r = resampler(resample_quality::draft, output_sr, input_sr, 1.0, 0.496); - univector<f64> resampled(len * output_sr / input_sr); - - const size_t destsize = r(resampled.data(), swept_sine); - - univector<i32> i32data = clamp(resampled.slice(0, destsize) * i32max, -i32max, +i32max); - univector2d<i32> data = { i32data }; - - auto wr = sequential_file_writer("audio_draft_quality.wav"); - audio_encode(wr, data, audioformat(data, output_sr)); - - plot_save("audio_draft_quality", "audio_draft_quality.wav", ""); - } - - return 0; -} diff --git a/examples/sample_rate_conversion.cpp b/examples/sample_rate_conversion.cpp @@ -0,0 +1,105 @@ +/** + * KFR (http://kfrlib.com) + * Copyright (C) 2016 D Levin + * See LICENSE.txt for details + */ + +// library_version() +#include <kfr/version.hpp> + +// print(), format() +#include <kfr/cometa/string.hpp> + +#include <kfr/math.hpp> + +// resample* +#include <kfr/dsp/sample_rate_conversion.hpp> + +// file* +#include <kfr/io/audiofile.hpp> + +// swept +#include <kfr/dsp/oscillators.hpp> + +// plot_save() +#include <kfr/io/python_plot.hpp> + +#include <iostream> + +using namespace kfr; + +constexpr size_t input_sr = 96000; +constexpr size_t output_sr = 44100; +constexpr size_t len = 96000 * 6; +constexpr f64 i32max = 2147483647.0; + +int main(int argc, char** argv) +{ + println(library_version()); + + const std::string options = "phaseresp=False"; + + univector<f64> swept_sine = swept(0.5, len); + + { + auto r = resampler(resample_quality::high, output_sr, input_sr, 1.0, 0.496); + univector<f64> resampled(len * output_sr / input_sr); + + const size_t destsize = r(resampled.data(), swept_sine); + + univector<i32> i32data = clamp(resampled.slice(0, destsize) * i32max, -i32max, +i32max); + univector2d<i32> data = { i32data }; + + auto wr = sequential_file_writer("audio_high_quality.wav"); + audio_encode(wr, data, audioformat(data, output_sr)); + + plot_save("audio_high_quality", "audio_high_quality.wav", ""); + } + + { + auto r = resampler(resample_quality::normal, output_sr, input_sr, 1.0, 0.496); + univector<f64> resampled(len * output_sr / input_sr); + + const size_t destsize = r(resampled.data(), swept_sine); + + univector<i32> i32data = clamp(resampled.slice(0, destsize) * i32max, -i32max, +i32max); + univector2d<i32> data = { i32data }; + + auto wr = sequential_file_writer("audio_normal_quality.wav"); + audio_encode(wr, data, audioformat(data, output_sr)); + + plot_save("audio_normal_quality", "audio_normal_quality.wav", ""); + } + + { + auto r = resampler(resample_quality::low, output_sr, input_sr, 1.0, 0.496); + univector<f64> resampled(len * output_sr / input_sr); + + const size_t destsize = r(resampled.data(), swept_sine); + + univector<i32> i32data = clamp(resampled.slice(0, destsize) * i32max, -i32max, +i32max); + univector2d<i32> data = { i32data }; + + auto wr = sequential_file_writer("audio_low_quality.wav"); + audio_encode(wr, data, audioformat(data, output_sr)); + + plot_save("audio_low_quality", "audio_low_quality.wav", ""); + } + + { + auto r = resampler(resample_quality::draft, output_sr, input_sr, 1.0, 0.496); + univector<f64> resampled(len * output_sr / input_sr); + + const size_t destsize = r(resampled.data(), swept_sine); + + univector<i32> i32data = clamp(resampled.slice(0, destsize) * i32max, -i32max, +i32max); + univector2d<i32> data = { i32data }; + + auto wr = sequential_file_writer("audio_draft_quality.wav"); + audio_encode(wr, data, audioformat(data, output_sr)); + + plot_save("audio_draft_quality", "audio_draft_quality.wav", ""); + } + + return 0; +} diff --git a/include/kfr/dsp/resample.hpp b/include/kfr/dsp/resample.hpp @@ -22,190 +22,4 @@ */ #pragma once -#include "../base/function.hpp" -#include "../base/memory.hpp" -#include "../base/reduce.hpp" -#include "../base/vec.hpp" -#include "window.hpp" - -namespace kfr -{ -namespace resample_quality -{ -constexpr csize_t<4> draft{}; -constexpr csize_t<6> low{}; -constexpr csize_t<8> normal{}; -constexpr csize_t<10> high{}; -} - -namespace internal -{ -template <typename T1, typename T2> -KFR_SINTRIN T1 resample_blackman(T1 n, T2 a) -{ - const T1 a0 = (1 - a) * 0.5; - const T1 a1 = 0.5; - const T1 a2 = a * 0.5; - n = n * c_pi<T1, 2>; - return a0 - a1 * cos(n) + a2 * cos(2 * n); -} - -template <typename T, size_t quality, KFR_ARCH_DEP> -struct resampler -{ - using itype = i64; - - constexpr static itype depth = static_cast<itype>(1 << (quality + 1)); - - resampler(itype interpolation_factor, itype decimation_factor, T scale = T(1), T cutoff = 0.49) - : input_position(0), output_position(0) - { - const i64 gcf = gcd(interpolation_factor, decimation_factor); - interpolation_factor /= gcf; - decimation_factor /= gcf; - - taps = depth * interpolation_factor; - order = size_t(depth * interpolation_factor - 1); - - this->interpolation_factor = interpolation_factor; - this->decimation_factor = decimation_factor; - - const itype halftaps = taps / 2; - filter = univector<T>(size_t(taps), T()); - delay = univector<T>(size_t(depth), T()); - - cutoff = cutoff / std::max(decimation_factor, interpolation_factor); - - for (itype j = 0, jj = 0; j < taps; j++) - { - filter[size_t(j)] = scale * 2 * interpolation_factor * cutoff * - sinc((jj - halftaps) * cutoff * c_pi<T, 2>) * - resample_blackman(T(jj) / T(taps - 1), T(0.16)); - jj += size_t(interpolation_factor); - if (jj >= taps) - jj = jj - taps + 1; - } - - const T s = reciprocal(sum(filter)) * interpolation_factor; - filter = filter * s; - } - KFR_INLINE size_t operator()(T* dest, size_t zerosize) - { - size_t outputsize = 0; - const itype srcsize = itype(zerosize); - - 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) - { - delay = zeros(); - } - else - { - delay.slice(0, size_t(depth - srcsize)) = delay.slice(size_t(srcsize)); - delay.slice(size_t(depth - srcsize)) = zeros(); - } - - input_position += srcsize; - output_position += outputsize; - return outputsize; - } - KFR_INLINE size_t operator()(T* dest, univector_ref<const T> src) - { - size_t outputsize = 0; - const itype srcsize = itype(src.size()); - - 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] = dotproduct(src.slice(size_t(srcindex - input_position), size_t(depth)), - tap_ptr /*, depth*/); - } - else - { - const itype prev_count = input_position - srcindex; - dest[i] = - dotproduct(delay.slice(size_t(depth - prev_count)), - tap_ptr /*, size_t(prev_count)*/) + - dotproduct( - src, tap_ptr.slice(size_t(prev_count), - size_t(depth - prev_count)) /*, size_t(depth - prev_count)*/); - } - } - } - if (srcsize >= depth) - { - delay = src.slice(size_t(srcsize - depth)); - } - else - { - delay.slice(0, size_t(depth - srcsize)) = delay.slice(size_t(srcsize)); - delay.slice(size_t(depth - srcsize)) = src; - } - - input_position += srcsize; - output_position += outputsize; - return outputsize; - } - itype taps; - size_t order; - itype interpolation_factor; - itype decimation_factor; - univector<T> filter; - univector<T> delay; - itype input_position; - itype output_position; -}; -} - -template <typename T, size_t quality> -inline internal::resampler<T, quality> resampler(csize_t<quality>, size_t interpolation_factor, - size_t decimation_factor, T scale = T(1), T cutoff = 0.49) -{ - using itype = typename internal::resampler<T, quality>::itype; - return internal::resampler<T, quality>(itype(interpolation_factor), itype(decimation_factor), scale, - cutoff); -} -} +#include "sample_rate_conversion.hpp" diff --git a/include/kfr/dsp/sample_rate_conversion.hpp b/include/kfr/dsp/sample_rate_conversion.hpp @@ -0,0 +1,223 @@ +/** + * Copyright (C) 2016 D Levin (http://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 http://www.kfrlib.com for details. + */ +#pragma once + +#include "../base/function.hpp" +#include "../base/memory.hpp" +#include "../base/reduce.hpp" +#include "../base/vec.hpp" +#include "window.hpp" + +namespace kfr +{ +namespace sample_rate_conversion_quality +{ +constexpr csize_t<4> draft{}; +constexpr csize_t<6> low{}; +constexpr csize_t<8> normal{}; +constexpr csize_t<10> high{}; +} + +namespace resample_quality = sample_rate_conversion_quality; + +namespace internal +{ +template <typename T1, typename T2> +KFR_SINTRIN T1 sample_rate_converter_blackman(T1 n, T2 a) +{ + const T1 a0 = (1 - a) * 0.5; + const T1 a1 = 0.5; + const T1 a2 = a * 0.5; + n = n * c_pi<T1, 2>; + return a0 - a1 * cos(n) + a2 * cos(2 * n); +} + +template <typename T, size_t quality, KFR_ARCH_DEP> +struct sample_rate_converter +{ + using itype = i64; + + constexpr static itype depth = static_cast<itype>(1 << (quality + 1)); + + sample_rate_converter(itype interpolation_factor, itype decimation_factor, T scale = T(1), T cutoff = 0.49) + : input_position(0), output_position(0) + { + const i64 gcf = gcd(interpolation_factor, decimation_factor); + interpolation_factor /= gcf; + decimation_factor /= gcf; + + taps = depth * interpolation_factor; + order = size_t(depth * interpolation_factor - 1); + + this->interpolation_factor = interpolation_factor; + this->decimation_factor = decimation_factor; + + const itype halftaps = taps / 2; + filter = univector<T>(size_t(taps), T()); + delay = univector<T>(size_t(depth), T()); + + cutoff = cutoff / std::max(decimation_factor, interpolation_factor); + + for (itype j = 0, jj = 0; j < taps; j++) + { + filter[size_t(j)] = scale * 2 * interpolation_factor * cutoff * + sinc((jj - halftaps) * cutoff * c_pi<T, 2>) * + sample_rate_converter_blackman(T(jj) / T(taps - 1), T(0.16)); + jj += size_t(interpolation_factor); + if (jj >= taps) + jj = jj - taps + 1; + } + + const T s = reciprocal(sum(filter)) * interpolation_factor; + filter = filter * s; + } + KFR_INLINE size_t operator()(T* dest, size_t zerosize) + { + size_t outputsize = 0; + const itype srcsize = itype(zerosize); + + 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) + { + delay = zeros(); + } + else + { + delay.slice(0, size_t(depth - srcsize)) = delay.slice(size_t(srcsize)); + delay.slice(size_t(depth - srcsize)) = zeros(); + } + + input_position += srcsize; + output_position += outputsize; + return outputsize; + } + KFR_INLINE size_t operator()(T* dest, univector_ref<const T> src) + { + size_t outputsize = 0; + const itype srcsize = itype(src.size()); + + 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] = dotproduct(src.slice(size_t(srcindex - input_position), size_t(depth)), + tap_ptr /*, depth*/); + } + else + { + const itype prev_count = input_position - srcindex; + dest[i] = + dotproduct(delay.slice(size_t(depth - prev_count)), + tap_ptr /*, size_t(prev_count)*/) + + dotproduct( + src, tap_ptr.slice(size_t(prev_count), + size_t(depth - prev_count)) /*, size_t(depth - prev_count)*/); + } + } + } + if (srcsize >= depth) + { + delay = src.slice(size_t(srcsize - depth)); + } + else + { + delay.slice(0, size_t(depth - srcsize)) = delay.slice(size_t(srcsize)); + delay.slice(size_t(depth - srcsize)) = src; + } + + input_position += srcsize; + output_position += outputsize; + return outputsize; + } + itype taps; + size_t order; + itype interpolation_factor; + itype decimation_factor; + univector<T> filter; + univector<T> delay; + itype input_position; + itype output_position; +}; +} + +template <typename T, size_t quality> +inline internal::sample_rate_converter<T, quality> sample_rate_converter(csize_t<quality>, size_t interpolation_factor, + size_t decimation_factor, T scale = T(1), T cutoff = 0.49) +{ + using itype = typename internal::sample_rate_converter<T, quality>::itype; + return internal::sample_rate_converter<T, quality>(itype(interpolation_factor), itype(decimation_factor), scale, + cutoff); +} + +// Deprecated in 0.9.2 +template <typename T, size_t quality> +inline internal::sample_rate_converter<T, quality> resampler(csize_t<quality>, size_t interpolation_factor, + size_t decimation_factor, T scale = T(1), T cutoff = 0.49) +{ + using itype = typename internal::sample_rate_converter<T, quality>::itype; + return internal::sample_rate_converter<T, quality>(itype(interpolation_factor), itype(decimation_factor), scale, + cutoff); +} +} diff --git a/sources.cmake b/sources.cmake @@ -75,6 +75,7 @@ set( ${PROJECT_SOURCE_DIR}/include/kfr/dsp/mixdown.hpp ${PROJECT_SOURCE_DIR}/include/kfr/dsp/oscillators.hpp ${PROJECT_SOURCE_DIR}/include/kfr/dsp/resample.hpp + ${PROJECT_SOURCE_DIR}/include/kfr/dsp/sample_rate_conversion.hpp ${PROJECT_SOURCE_DIR}/include/kfr/dsp/speaker.hpp ${PROJECT_SOURCE_DIR}/include/kfr/dsp/units.hpp ${PROJECT_SOURCE_DIR}/include/kfr/dsp/waveshaper.hpp