diff --git a/CMakeLists.txt b/CMakeLists.txt index 2ea3ef0..2b3d765 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -38,6 +38,7 @@ option(EAR_EXAMPLES "Build examples" ${IS_ROOT_PROJECT}) option(EAR_PACKAGE_AND_INSTALL "Package and install libear" ${IS_ROOT_PROJECT}) option(EAR_USE_INTERNAL_EIGEN "should we use our own version of Eigen, or find one with find_package?" TRUE) option(EAR_USE_INTERNAL_XSIMD "should we use our own version of xsimd, or find one with find_package?" TRUE) +option(EAR_SIMD "try to use SIMD extensions" TRUE) set(INSTALL_LIB_DIR lib CACHE PATH "Installation directory for libraries") set(INSTALL_BIN_DIR bin CACHE PATH "Installation directory for executables") set(INSTALL_INCLUDE_DIR include CACHE PATH "Installation directory for header files") @@ -100,6 +101,7 @@ add_feature_info(EAR_EXAMPLES ${EAR_EXAMPLES} "Build examples") add_feature_info(EAR_PACKAGE_AND_INSTALL ${EAR_PACKAGE_AND_INSTALL} "Package and install libear") add_feature_info(EAR_USE_INTERNAL_EIGEN ${EAR_USE_INTERNAL_EIGEN} "use internal version of Eigen") add_feature_info(EAR_USE_INTERNAL_XSIMD ${EAR_USE_INTERNAL_XSIMD} "use internal version of xsimd") +add_feature_info(EAR_SIMD ${EAR_SIMD} "try to use SIMD extensions") feature_summary(WHAT ALL) if(EAR_PACKAGE_AND_INSTALL) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d4db70b..2afa40f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -27,7 +27,9 @@ add_library(ear hoa/hoa.cpp layout.cpp screen.cpp - object_based/extent.cpp + object_based/polar_extent.cpp + object_based/polar_extent_scalar.cpp + object_based/polar_extent_simd.cpp object_based/gain_calculator_objects.cpp warnings.cpp ${CMAKE_CURRENT_BINARY_DIR}/resources.hpp @@ -81,6 +83,62 @@ else() endif() set_target_properties(ear PROPERTIES CXX_EXTENSIONS OFF) +include(CheckCompilerFlag) + +set(PER_ARCH_SOURCES object_based/polar_extent_simd_instance.cpp) +function(add_arch_lib TARGET ARCH FLAGS) + check_compiler_flag(CXX "${FLAGS}" ear_has_flag_${ARCH}) # var is cached and must be unique + if (ear_has_flag_${ARCH}) + list(APPEND XSIMD_ARCHS xsimd::${ARCH}) + add_library(${TARGET}_${ARCH} OBJECT ${PER_ARCH_SOURCES}) + target_sources(${TARGET} PRIVATE $) + + # copy relevant target properties to feature-specific objects + foreach(PROPERTY COMPILE_DEFINITIONS COMPILE_FEATURES INCLUDE_DIRECTORIES LINK_LIBRARIES + POSITION_INDEPENDENT_CODE CXX_VISIBILITY_PRESET C_VISIBILITY_PRESET VISIBILITY_INLINES_HIDDEN) + get_target_property(VALUE ${TARGET} ${PROPERTY}) + if (NOT VALUE STREQUAL "VALUE-NOTFOUND") + set_property(TARGET ${TARGET}_${ARCH} PROPERTY ${PROPERTY} ${VALUE}) + endif() + endforeach() + + target_compile_definitions(${TARGET}_${ARCH} PRIVATE XSIMD_ARCH=xsimd::${ARCH}) + + target_compile_options(${TARGET}_${ARCH} PRIVATE ${FLAGS}) + endif() + + set(XSIMD_ARCHS "${XSIMD_ARCHS}" PARENT_SCOPE) +endfunction() + +if (EAR_SIMD) + # add_arch_lib must be called from best to worst instruction set + if (CMAKE_SYSTEM_PROCESSOR MATCHES "x86_64|AMD64|X86|i686") + if (MSVC) + add_arch_lib(ear avx512bw "/arch:AVX512") + add_arch_lib(ear avx2_fma "/arch:AVX2") + add_arch_lib(ear avx "/arch:AVX") + else() + add_arch_lib(ear avx512bw "-mavx512bw;-mavx512dq;-mavx512cd;-mavx512f;-mfma;-mavx2;-mavx") + add_arch_lib(ear avx2_fma "-mfma;-mavx2;-mavx") + add_arch_lib(ear avx "-mavx") + add_arch_lib(ear sse4_2 "-msse4.2") + endif() + add_arch_lib(ear default_arch "") + elseif(CMAKE_SYSTEM_PROCESSOR MATCHES "ARM.*|arm.*|aarch64.*") + # for arm, flags are not generally needed as the FP extension support is + # closely tied to the ISA, so just build the default and fallback + # + # this can go wrong if the compiler doesn't actually support neon, but it's + # tricky to detect that + add_arch_lib(ear default_arch "") + endif() +endif() + +list(APPEND XSIMD_ARCHS xsimd::generic_for_dispatch) + +list(JOIN XSIMD_ARCHS "," XSIMD_ARCHS_STR) +target_compile_definitions(ear PRIVATE "XSIMD_ARCHS=${XSIMD_ARCHS_STR}") + include(GenerateExportHeader) generate_export_header(ear EXPORT_FILE_NAME ${PROJECT_BINARY_DIR}/generated/export.hpp diff --git a/src/common/helpers/xsimd_extension.hpp b/src/common/helpers/xsimd_extension.hpp new file mode 100644 index 0000000..92c7ab2 --- /dev/null +++ b/src/common/helpers/xsimd_extension.hpp @@ -0,0 +1,22 @@ +#pragma once +#include + +namespace xsimd { + // aliases for things which will be used via preprocessor defines, to avoid + // quoting issues + using avx2_fma = fma3; + + // a version of generic which works with xsimd::dispatch + struct generic_for_dispatch : generic { + static constexpr bool supported() noexcept { return true; } + static constexpr bool available() noexcept { return true; } + static constexpr bool requires_alignment() noexcept { return false; } + static constexpr unsigned version() noexcept { + return generic::version(0, 0, 0); + } + static constexpr std::size_t alignment() noexcept { return 0; } + static constexpr char const* name() noexcept { + return "generic_for_dispatch"; + } + }; +} // namespace xsimd diff --git a/src/object_based/extent.cpp b/src/object_based/extent.cpp deleted file mode 100644 index f5d0fc1..0000000 --- a/src/object_based/extent.cpp +++ /dev/null @@ -1,242 +0,0 @@ -#include "extent.hpp" - -#include -#include -#include -#include -#include -#include "../common/geom.hpp" -#include "../common/helpers/eigen_helpers.hpp" - -const double PI = boost::math::constants::pi(); - -namespace ear { - - /** @brief Normalise position or return {0,1,0}. - * - * @param position Position to normalise. - * - * @returns normalised position - */ - Eigen::Vector3d safeNormPosition(Eigen::Vector3d position) { - double norm = position.norm(); - if (norm < 1e-10) { - return Eigen::Vector3d{0.0, 1.0, 0.0}; - } else { - return position / norm; - } - } - - double extentMod(double extent, double distance) { - double minSize = 0.2; - double size = interp(extent, Eigen::Vector2d(0.0, 360.0), - Eigen::Vector2d(minSize, 1.0)); - double extent1 = 4.0 * degrees(atan2(size, 1.0)); - return interp(4.0 * degrees(atan2(size, distance)), - Eigen::Vector3d(0.0, extent1, 360.0), - Eigen::Vector3d(0.0, extent, 360.0)); - } - - Eigen::Matrix3d calcBasis(Eigen::Vector3d position) { - position = safeNormPosition(position); - double az = azimuth(position); - double el = elevation(position); - - // points near the poles have indeterminate azimuth; assume 0 - if (std::abs(el) > (90.0 - 1e-5)) { - az = 0.0; - } - return localCoordinateSystem(az, el); - } - - Eigen::Vector3d cartOnBasis(Eigen::Matrix3d basis, double azimuth, - double elevation) { - Eigen::RowVector3d cartPosRel{sin(azimuth) * cos(elevation), - cos(azimuth) * cos(elevation), - sin(elevation)}; - return cartPosRel * basis; - } - - std::pair azimuthElevationOnBasis( - Eigen::Matrix3d basis, Eigen::RowVector3d position) { - // project onto each basis, and clip components to keep asin happy - Eigen::Vector3d components = - (position * basis.transpose()).cwiseMin(1.0).cwiseMax(-1.0); - - double azimuth = atan2(components(0), components(1)); - double elevation = asin(components(2)); - - return std::make_pair(azimuth, elevation); - } - - WeightingFunction::WeightingFunction(Eigen::Vector3d position, double width, - double height) { - _width = radians(width) / 2; - _height = radians(height) / 2; - - // basis vectors to rotate the vsource positions towards position - Eigen::Matrix3d basises = calcBasis(position); - - _circleRadius = std::min(_width, _height); - - // Flip the width and the height such that it is always wider than it is - // high from here in. - if (_height > _width) { - std::swap(_height, _width); - _flippedBasis = basises.colwise().reverse(); - } else { - _flippedBasis = basises; - } - - // modify the width to make it meet at the back. - double widthFull = PI + _height; - // interpolate to this from a width of pi/2 to pi - double widthMod = interp(_width, Eigen::Vector3d{0.0, PI / 2.0, PI}, - Eigen::Vector3d{0.0, PI / 2.0, widthFull}); - // apply this fully for a height of less than pi/4; tail off until pi/2 - _width = interp(_height, Eigen::Vector4d{0, PI / 4.0, PI / 2.0, PI}, // - Eigen::Vector4d{widthMod, widthMod, _width, _width}); // - - // angle of the circle centres from the source position; width is to the - // end of the rectangle. - _circlePos = _width - _circleRadius; - - // Cartesian circle centres - _circlePositions << cartOnBasis(_flippedBasis, -_circlePos, 0.0), - cartOnBasis(_flippedBasis, _circlePos, 0.0); - } - - double WeightingFunction::operator()(Eigen::Vector3d position) const { - // Flipped azimuths and elevations; the straight edges are always along - // azimuth lines. - double azimuth, elevation; - std::tie(azimuth, elevation) = - azimuthElevationOnBasis(_flippedBasis, position); - - // The distance is the angle away from the defined shape; 0 or negative is - // inside. - double distance = 0.0; - - // for the straight lines - if (std::abs(azimuth) <= _circlePos) { - distance = std::abs(elevation) - _circleRadius; - } else { - // distance from the closest circle centre - size_t nearest_circle = azimuth < 0 ? 0 : 1; - double angle = - position.transpose() * _circlePositions.col(nearest_circle); - double circleDistance = acos(boost::algorithm::clamp(angle, -1.0, 1.0)); - distance = circleDistance - _circleRadius; - } - // fade the weight from one to zero over fadeWidth - return interp(distance, Eigen::Vector2d{0.0, radians(_fadeWidth)}, - Eigen::Vector2d{1.0, 0.0}); - } - - SpreadingPanner::SpreadingPanner(std::shared_ptr psp, - int nRows) - : _psp(psp), _nRows(nRows) { - _panningPositions = _generatePanningPositionsEven(); - _panningPositionsResults = _generatePanningPositionsResults(); - } - - Eigen::VectorXd SpreadingPanner::panningValuesForWeight( - const WeightingFunction& weightFunc) { - Eigen::VectorXd weights(_panningPositions.rows()); - for (int i = 0; i < _panningPositions.rows(); ++i) { - weights(i) = weightFunc(_panningPositions.row(i)); - } - Eigen::VectorXd totalPv = weights.transpose() * _panningPositionsResults; - return totalPv / totalPv.norm(); - } - - Eigen::MatrixXd SpreadingPanner::_generatePanningPositionsEven() { - Eigen::VectorXd elevations = - Eigen::VectorXd::LinSpaced(_nRows, -90.0, 90.0); - Eigen::MatrixXd positions(0, 3); - - for (double el : elevations) { - double radius = cos(radians(el)); - double perimiter = 2 * PI * radius; - double perimiter_centre = 2 * PI; - - int nPoints = static_cast( - std::round((perimiter / perimiter_centre) * 2 * (_nRows - 1))); - if (nPoints == 0) { - nPoints = 1; - } - Eigen::VectorXd azimuths = - Eigen::VectorXd::LinSpaced(nPoints + 1, 0.0, 360.0); - for (int i = 0; i < azimuths.size() - 1; ++i) { - double az = azimuths(i); - positions.conservativeResize(positions.rows() + 1, Eigen::NoChange); - positions.row(positions.rows() - 1) = cart(az, el, 1.0); - } - } - return positions; - } - - Eigen::MatrixXd SpreadingPanner::_generatePanningPositionsResults() { - Eigen::MatrixXd results(_panningPositions.rows(), - _psp->numberOfOutputChannels()); - for (int i = 0; i < _panningPositions.rows(); ++i) { - results.row(i) = _psp->handle(_panningPositions.row(i)).get(); - } - return results; - } - - PolarExtentPanner::PolarExtentPanner(std::shared_ptr psp) - : _psp(psp), _spreadingPanner(SpreadingPanner(psp, _nRows)){}; - - Eigen::VectorXd PolarExtentPanner::calcPvSpread(Eigen::Vector3d position, - double width, double height) { - // When calculating the spread panning values the width and height are - // set to at least fade_width. For sizes where any of the dimensions is - // less than this, interpolate linearly between the point and spread - // panning values. - double ammount_spread = - interp(std::max(width, height), Eigen::Vector2d(0.0, _fadeWidth), - Eigen::Vector2d(0.0, 1.0)); - double ammountPoint = 1.0 - ammount_spread; - Eigen::ArrayXd pv = Eigen::ArrayXd::Zero(_psp->numberOfOutputChannels()); - if (ammountPoint > 1e-10) { - pv += ammountPoint * _psp->handle(position).get().array().square(); - } - if (ammount_spread > 1e-10) { - // minimum width and height as above - width = std::max(width, _fadeWidth / 2.0); - height = std::max(height, _fadeWidth / 2.0); - - WeightingFunction weightingFunction(position, width, height); - Eigen::VectorXd panning_values = - _spreadingPanner.panningValuesForWeight(weightingFunction); - pv += ammount_spread * panning_values.array().square(); - } - return pv.sqrt().matrix(); - } - - Eigen::VectorXd PolarExtentPanner::handle(Eigen::Vector3d position, - double width, double height, - double depth) { - double distance = position.norm(); - - if (depth != 0.0) { - double distanceMin = distance - depth / 2.0; - double distanceMax = distance + depth / 2.0; - distanceMin = (distanceMin < 0) ? 0.0 : distanceMin; - distanceMax = (distanceMax < 0) ? 0.0 : distanceMax; - Eigen::VectorXd pvsMin = - calcPvSpread(position, extentMod(width, distanceMin), - extentMod(height, distanceMin)); - Eigen::VectorXd pvsMax = - calcPvSpread(position, extentMod(width, distanceMax), - extentMod(height, distanceMax)); - return ((pvsMin.array().square() + pvsMax.array().square()) / 2.0).sqrt(); - } else { - Eigen::VectorXd pvs = calcPvSpread(position, extentMod(width, distance), - extentMod(height, distance)); - return pvs; - } - } - -} // namespace ear diff --git a/src/object_based/extent.hpp b/src/object_based/extent.hpp deleted file mode 100644 index 469d1f5..0000000 --- a/src/object_based/extent.hpp +++ /dev/null @@ -1,143 +0,0 @@ -#pragma once -#include -#include -#include "../common/point_source_panner.hpp" - -namespace ear { - - class ExtentPanner {}; - - /** @brief Modify an extent parameter given a distance. - * - * A right triangle if formed, with the adjacent edge being the distance, and - * the opposite edge being determined from the extent. The angle formed is - * then used to determine the new extent. - * - * - at distance=0, the extent is always 360 - * - at distance=1, the original extent is used - * - at distance>1, the extent decreases - * - in 0 < distance < 1, the extent changes more steeply around 0 for smaller - * extents - */ - double extentMod(double extent, double distance); - - /** @brief Calculate basis vectors that rotate (0, 1, 0) - onto source_pos. */ - Eigen::Matrix3d calcBasis(Eigen::Vector3d position); - - /** @brief Polar to Cartesian in radians with no distance, in a given basis.*/ - Eigen::Vector3d cartOnBasis(Eigen::Matrix3d basis, double azimuth, - double elevation); - - std::pair azimuthElevationOnBasis( - Eigen::Matrix3d basis, Eigen::RowVector3d position); - - /** @brief Weighting function for spread sources. - * - * The weighting function is one inside a region approximately determined by a - * width x height rectangle in azimuth-elevation space, with maximally-sized - * rounded corners; the shape of the corners is calculated using the vector - * angle from their centres (always directly above or below the source - * position) so as to avoid issues at the poles. - * - * The two straight edges of the rectangle are always parallel in Cartesian - * space; this is achieved by following azimuth lines; for tall sources, the - * whole coordinate system is rotated 90 degrees about the source position to - * achieve this. - * - * Note that for sources where width == height, this degrades to a circular - * region relative to the source position. - * - * To make the two ends meet, the width is adjusted such that a width of 180 - * degrees is mapped to width + height. - * - */ - class WeightingFunction { - public: - EIGEN_MAKE_ALIGNED_OPERATOR_NEW - /** @brief Ctor - * - * @param position Centre of the extent. - * @param width Width of the extent in degrees from one edge to the other. - * @param height Height of the extent in degrees from one edge to the other. - */ - WeightingFunction(Eigen::Vector3d position, double width, double height); - - /** @brief Calculate weight for position */ - double operator()(Eigen::Vector3d position) const; - - private: - const double _fadeWidth = 10.0; - double _width; - double _height; - double _circleRadius; - Eigen::Matrix3d _flippedBasis; - double _circlePos; - Eigen::Matrix _circlePositions; - }; - - class SpreadingPanner { - public: - SpreadingPanner(std::shared_ptr psp, int nRows); - - /** @brief Panning values for a given weighting function. - * - * @param weightFunc function from Cartesian position to weight in range - * (0, 1) - * - * @return panning value for each speaker. - */ - Eigen::VectorXd panningValuesForWeight(const WeightingFunction& weightFunc); - - private: - /** @brief Generate points spread evenly on the sphere. - * - * Based on - * http://web.archive.org/web/20150108040043/http://www.math.niu.edu/~rusin/known-math/95/equispace.elect - * - * @param nRows number of rows to place on sphere, e.g. 37 for 5 degree - * spacing - * - * @returns cartesian array. - */ - Eigen::MatrixXd _generatePanningPositionsEven(); - Eigen::MatrixXd _generatePanningPositionsResults(); - - std::shared_ptr _psp; - int _nRows; - Eigen::MatrixXd _panningPositions; - Eigen::MatrixXd _panningPositionsResults; - }; - - class PolarExtentPanner : public ExtentPanner { - public: - PolarExtentPanner(std::shared_ptr psp); - - /** @brief Calculate loudspeaker gains given position and extent parameters. - * - * @param position Cartesian source position - * @param width block format width parameter - * @param height block format height parameter - * @param depth block format depth parameter - * - * @returns loudspeaker gains for each channel - */ - Eigen::VectorXd handle(Eigen::Vector3d position, double width, - double height, double depth); - - /** @brief Calculate the speaker panning values for the position, width, and - * height of a source; this just deals with the positioning and spreading. - */ - Eigen::VectorXd calcPvSpread(Eigen::Vector3d position, double width, - double height); - - private: - const int _nRows = 37; // 5 degrees per row - const double _fadeWidth = 10.0; // degrees - std::shared_ptr _psp; - SpreadingPanner _spreadingPanner; - }; - - class CartesianExtentPanner : public ExtentPanner {}; - -} // namespace ear diff --git a/src/object_based/gain_calculator_objects.cpp b/src/object_based/gain_calculator_objects.cpp index 7074d26..4132d99 100644 --- a/src/object_based/gain_calculator_objects.cpp +++ b/src/object_based/gain_calculator_objects.cpp @@ -29,9 +29,10 @@ namespace ear { _screenEdgeLockHandler(ScreenEdgeLockHandler(_layout.screen())), _screenScaleHandler(ScreenScaleHandler(_layout.screen())), _channelLockHandler(ChannelLockHandler(_layout.withoutLfe())), - _polarExtentPanner(PolarExtentPanner(_pointSourcePanner)), + _polarExtentPanner(_pointSourcePanner), _zoneExclusionHandler(ZoneExclusionHandler(_layout.withoutLfe())), - _isLfe(copy_vector(layout.isLfe())){}; + _isLfe(copy_vector(layout.isLfe())), + _pvTmp(_pointSourcePanner->numberOfOutputChannels()){}; void GainCalculatorObjectsImpl::calculate(const ObjectsTypeMetadata& metadata, OutputGains& direct, @@ -47,12 +48,12 @@ namespace ear { Eigen::Vector3d position = toCartesianVector3d(boost::get(metadata.position)); - auto pv = _polarExtentPanner.handle(position, metadata.width, - metadata.height, metadata.depth); - pv *= metadata.gain; + _polarExtentPanner.handle(position, metadata.width, metadata.height, + metadata.depth, _pvTmp); + _pvTmp *= metadata.gain; Eigen::VectorXd pv_full = Eigen::VectorXd::Zero(_isLfe.size()); - mask_write(pv_full, !_isLfe, pv); + mask_write(pv_full, !_isLfe, _pvTmp); // apply diffuse split direct.write_vector(pv_full * std::sqrt(1.0 - metadata.diffuse)); diff --git a/src/object_based/gain_calculator_objects.hpp b/src/object_based/gain_calculator_objects.hpp index 8b0bae4..cdb5922 100644 --- a/src/object_based/gain_calculator_objects.hpp +++ b/src/object_based/gain_calculator_objects.hpp @@ -9,7 +9,7 @@ #include "ear/layout.hpp" #include "ear/metadata.hpp" #include "ear/warnings.hpp" -#include "extent.hpp" +#include "polar_extent.hpp" #include "zone_exclusion.hpp" namespace ear { @@ -36,9 +36,10 @@ namespace ear { ScreenEdgeLockHandler _screenEdgeLockHandler; ScreenScaleHandler _screenScaleHandler; ChannelLockHandler _channelLockHandler; - PolarExtentPanner _polarExtentPanner; + PolarExtent _polarExtentPanner; ZoneExclusionHandler _zoneExclusionHandler; Eigen::Array _isLfe; + Eigen::VectorXd _pvTmp; }; } // namespace ear diff --git a/src/object_based/polar_extent.cpp b/src/object_based/polar_extent.cpp new file mode 100644 index 0000000..c236fd8 --- /dev/null +++ b/src/object_based/polar_extent.cpp @@ -0,0 +1,303 @@ +#include "polar_extent.hpp" +#include +#include +#include "../common/geom.hpp" +#include "../common/helpers/eigen_helpers.hpp" + +const double PI = boost::math::constants::pi(); + +namespace ear { + + namespace { + constexpr double fadeWidth = 10.0; + constexpr int nRows = 37; // 5 degrees per row + + Eigen::MatrixXd generatePanningPositionsEven(int nRows) { + Eigen::VectorXd elevations = + Eigen::VectorXd::LinSpaced(nRows, -90.0, 90.0); + Eigen::MatrixXd positions(0, 3); + + for (double el : elevations) { + double radius = cos(radians(el)); + double perimiter = 2 * PI * radius; + double perimiter_centre = 2 * PI; + + int nPoints = static_cast( + std::round((perimiter / perimiter_centre) * 2 * (nRows - 1))); + if (nPoints == 0) { + nPoints = 1; + } + Eigen::VectorXd azimuths = + Eigen::VectorXd::LinSpaced(nPoints + 1, 0.0, 360.0); + for (int i = 0; i < azimuths.size() - 1; ++i) { + double az = azimuths(i); + positions.conservativeResize(positions.rows() + 1, Eigen::NoChange); + positions.row(positions.rows() - 1) = cart(az, el, 1.0); + } + } + return positions; + } + + Eigen::MatrixXd generatePanningPositionsResults( + const std::shared_ptr &psp, + const Eigen::Ref &positions) { + Eigen::MatrixXd results(psp->numberOfOutputChannels(), positions.rows()); + for (int i = 0; i < positions.rows(); ++i) { + results.col(i) = psp->handle(positions.row(i)).get(); + } + return results; + } + + size_t round_up(size_t n, size_t multiple) { + return ((n + multiple - 1) / multiple) * multiple; + } + + /** @brief Normalise position or return {0,1,0}. + * + * @param position Position to normalise. + * + * @returns normalised position + */ + static Eigen::Vector3d safeNormPosition(Eigen::Vector3d position) { + double norm = position.norm(); + if (norm < 1e-10) { + return Eigen::Vector3d{0.0, 1.0, 0.0}; + } else { + return position / norm; + } + } + + double extentMod(double extent, double distance) { + double minSize = 0.2; + double size = interp(extent, Eigen::Vector2d(0.0, 360.0), + Eigen::Vector2d(minSize, 1.0)); + double extent1 = 4.0 * degrees(atan2(size, 1.0)); + return interp(4.0 * degrees(atan2(size, distance)), + Eigen::Vector3d(0.0, extent1, 360.0), + Eigen::Vector3d(0.0, extent, 360.0)); + } + } // namespace + + Eigen::Matrix3d calcBasis(Eigen::Vector3d position) { + position = safeNormPosition(position); + double az = azimuth(position); + double el = elevation(position); + + // points near the poles have indeterminate azimuth; assume 0 + if (std::abs(el) > (90.0 - 1e-5)) { + az = 0.0; + } + return localCoordinateSystem(az, el); + } + + PolarExtentCore::~PolarExtentCore() {} + + PolarExtent::PolarExtent(const std::shared_ptr &psp_) + : PolarExtent(psp_, get_polar_extent_core()) {} + + PolarExtent::PolarExtent(const std::shared_ptr &psp_, + std::unique_ptr core_impl_) + : core_impl(std::move(core_impl_)), + psp(psp_), + pvsMin(psp->numberOfOutputChannels()), + pvsMax(psp->numberOfOutputChannels()) { + auto panning_positions = generatePanningPositionsEven(nRows); + auto panning_results = + generatePanningPositionsResults(psp, panning_positions); + + size_t num_speakers = panning_results.rows(); + size_t num_points = panning_results.cols(); + size_t batch_size = core_impl->batch_size(); + ctx.real_num_speakers = num_speakers; + ctx.num_speakers = round_up(num_speakers, batch_size); + ctx.num_points = round_up(num_points, batch_size); + + // use regular malloc and free for memory management. reasoning: + // - the parts compiled with different architectures should not use any + // in-line functions which could get merged; therefore using any kind of + // smart pointer in the interface is a bad plan + // - we could use a smart pointer and extract a float pointer, but that + // doesn't make it much cleaner, would still break copying, and would + // ultimately be more code + // - using smart pointers would be awkward anyway, as the alignment is + // known at run-time + auto malloc_float = [&](size_t n) { + auto mem = (extent_float_t *)xsimd::aligned_malloc( + sizeof(extent_float_t) * n, core_impl->alignment()); + for (size_t i = 0; i < n; i++) mem[i] = 0.0; + return mem; + }; + + ctx.xs = malloc_float(ctx.num_points); + ctx.ys = malloc_float(ctx.num_points); + ctx.zs = malloc_float(ctx.num_points); + + ctx.panning_results = malloc_float(ctx.num_points * ctx.num_speakers); + size_t num_batches = ctx.num_points / batch_size; + ctx.summed_panning_results = malloc_float(num_batches * ctx.num_speakers); + + ctx.results = malloc_float(ctx.num_speakers); + + for (size_t p = 0; p < num_points; p++) { + ctx.xs[p] = static_cast(panning_positions(p, 0)); + ctx.ys[p] = static_cast(panning_positions(p, 1)); + ctx.zs[p] = static_cast(panning_positions(p, 2)); + + for (size_t s = 0; s < num_speakers; s++) { + ctx.panning_results[s * ctx.num_points + p] = + static_cast(panning_results(s, p)); + size_t batch = p / batch_size; + ctx.summed_panning_results[batch * ctx.num_speakers + s] += + static_cast(panning_results(s, p)); + } + } + + for (size_t p = num_points; p < ctx.num_points; p++) { + // the results of these will be ignored; make them consistent with the + // rest of the batch + ctx.xs[p] = + static_cast(panning_positions(num_points - 1, 0)); + ctx.ys[p] = + static_cast(panning_positions(num_points - 1, 1)); + ctx.zs[p] = + static_cast(panning_positions(num_points - 1, 2)); + } + } + + PolarExtent::~PolarExtent() { + xsimd::aligned_free(ctx.xs); + xsimd::aligned_free(ctx.ys); + xsimd::aligned_free(ctx.zs); + xsimd::aligned_free(ctx.panning_results); + xsimd::aligned_free(ctx.summed_panning_results); + xsimd::aligned_free(ctx.results); + } + + void PolarExtent::setup_angle_to_weight(double start_angle, + double end_angle) const { + ear_assert(start_angle >= 0, "start angle should be +ve"); + ear_assert(end_angle >= 0, "end angle should be +ve"); + ear_assert(end_angle > start_angle, "end angle should be > start angle"); + + ctx.cos_start_angle = static_cast( + start_angle < PI ? std::cos(start_angle) : -1.0); + ctx.cos_end_angle = static_cast( + end_angle < PI ? std::cos(end_angle) : -(1.0 + 1e-6)); + // sin is only used for less than PI/2 + ctx.sin_start_angle = static_cast( + start_angle < PI / 2 ? std::sin(start_angle) : 1.0); + ctx.sin_end_angle = static_cast( + end_angle < PI / 2 ? std::sin(end_angle) : 1.0 + 1e-6); + // between start and end angle, we want: + // out = (angle - end_angle) / (start_angle - end_angle) + // therefore the slope is: + ctx.m = static_cast(1.0 / (start_angle - end_angle)); + // out = m * (angle - end_angle) + // out = m * angle - m * end_angle + // intercept is: + ctx.c = static_cast(-ctx.m * end_angle); + } + + void PolarExtent::setup_weighting_function(Eigen::Vector3d position, + double width, + double height) const { + width = radians(width) / 2; + height = radians(height) / 2; + + // basis vectors to rotate the vsource positions towards position + Eigen::Matrix3d m = calcBasis(position); + + // Flip the width and the height such that it is always wider than it is + // high from here in. + if (height > width) { + std::swap(height, width); + // rotate rather than flipping x and y to preserve triangle winding + Eigen::Matrix3d flip; + flip << 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, -1.0, 0.0, 0.0; + m = flip * m; + } + + for (size_t i = 0; i < 3; i++) + for (size_t j = 0; j < 3; j++) + ctx.flippedBasis[i * 3 + j] = static_cast(m(i, j)); + + // modify the width to make it meet at the back. + double widthFull = PI + height; + // interpolate to this from a width of pi/2 to pi + double widthMod = interp(width, Eigen::Vector3d{0.0, PI / 2.0, PI}, + Eigen::Vector3d{0.0, PI / 2.0, widthFull}); + // apply this fully for a height of less than pi/4; tail off until pi/2 + width = interp(height, Eigen::Vector4d{0, PI / 4.0, PI / 2.0, PI}, + Eigen::Vector4d{widthMod, widthMod, width, width}); + + ctx.is_circular = (width - height) < 1e-6; + + double circlePos = width - height; + + ctx.right_circle_centre[0] = + static_cast(std::sin(circlePos)); + ctx.right_circle_centre[1] = + static_cast(std::cos(circlePos)); + ctx.circle_test[0] = static_cast(-std::cos(circlePos)); + ctx.circle_test[1] = static_cast(std::sin(circlePos)); + + setup_angle_to_weight(height, height + radians(fadeWidth)); + } + + void PolarExtent::calc_pv_spread(Eigen::Vector3d position, double width, + double height, + Eigen::Ref out) const { + // When calculating the spread panning values the width and height are + // set to at least fade_width. For sizes where any of the dimensions is + // less than this, interpolate linearly between the point and spread + // panning values. + double ammount_spread = + interp(std::max(width, height), Eigen::Vector2d(0.0, fadeWidth), + Eigen::Vector2d(0.0, 1.0)); + double ammountPoint = 1.0 - ammount_spread; + out.setZero(); + if (ammountPoint > 1e-10) { + out.array() += + ammountPoint * psp->handle(position).get().array().square(); + } + if (ammount_spread > 1e-10) { + // minimum width and height as above + width = std::max(width, fadeWidth / 2.0); + height = std::max(height, fadeWidth / 2.0); + + setup_weighting_function(position, width, height); + core_impl->run(ctx); + + auto results_map = Eigen::Map>( + ctx.results, psp->numberOfOutputChannels()); + results_map *= 1.0 / results_map.norm(); + + out.array() += + ammount_spread * results_map.array().square().cast(); + } + out.array() = out.array().sqrt(); + } + + void PolarExtent::handle(Eigen::Vector3d position, double width, + double height, double depth, + Eigen::Ref out) const { + double distance = position.norm(); + + if (depth != 0.0) { + double distanceMin = distance - depth / 2.0; + double distanceMax = distance + depth / 2.0; + distanceMin = (distanceMin < 0) ? 0.0 : distanceMin; + distanceMax = (distanceMax < 0) ? 0.0 : distanceMax; + calc_pv_spread(position, extentMod(width, distanceMin), + extentMod(height, distanceMin), pvsMin); + calc_pv_spread(position, extentMod(width, distanceMax), + extentMod(height, distanceMax), pvsMax); + out = ((pvsMin.array().square() + pvsMax.array().square()) / 2.0) + .sqrt() + .matrix(); + } else { + calc_pv_spread(position, extentMod(width, distance), + extentMod(height, distance), out); + } + } +} // namespace ear diff --git a/src/object_based/polar_extent.hpp b/src/object_based/polar_extent.hpp new file mode 100644 index 0000000..570058f --- /dev/null +++ b/src/object_based/polar_extent.hpp @@ -0,0 +1,78 @@ +#pragma once +#include "../common/point_source_panner.hpp" +#include "polar_extent_core.hpp" + +namespace ear { + + Eigen::Matrix3d calcBasis(Eigen::Vector3d position); + + /// polar extent implementation + /// + /// this has a different structure from the reference implementation for + /// efficiency: rather than the SpreadingPanner, which sums pre-calculated + /// gains corresponding to points on the sphere, with weights calculated by a + /// WeightingFunction instance, we have instead: + /// + /// - a core, which calculates the weight of each point according to some + /// some pre-processed parameters, and sums the weighted gains for each + /// point (without normalisation) + /// + /// - the PolarExtent class, which sets up the conditions for the core + /// (generating positions and corresponding gains, and pre-processing the + /// parameters), runs the core, and post-processes the results + /// + /// this makes it possible to have core implementations that process multiple + /// points at once + /// + /// the various moving parts of this are as follows: + /// + /// - PolarExtentCore in polar_extent_core.hpp is the interface that core + /// implementations implement + /// + /// - PolarExtentCoreContext in polar_extent_core.hpp is the structure passed + /// between PolarExtent and the core implementations + /// + /// - get_polar_extent_core in polar_extent_core.hpp constructs the best core + /// implementation for this platform + /// + /// - PolarExtentCoreScalar in polar_extent_scalar.cpp is a scalar core + /// implementation to fall back on + /// + /// - PolarExtentCoreSimd in polar_extent_simd.hpp is a SIMD core + /// implementation using xsimd, parametrised by the instruction set to use + /// + /// - polar_extent_simd_instance.cpp is compiled once for each supported + /// architecture (see src/CMakeLists.txt), only instantiating + /// PolarExtentCoreSimd and an instance of get_polar_extent_core_arch for + /// that architecture + class PolarExtent { + public: + PolarExtent(const std::shared_ptr &psp); + // for testing + PolarExtent(const std::shared_ptr &psp, + std::unique_ptr core_impl); + ~PolarExtent(); + + PolarExtent(PolarExtent const &) = delete; + PolarExtent &operator=(PolarExtent const &) = delete; + + void handle(Eigen::Vector3d position, double width, double height, + double depth, Eigen::Ref out) const; + + private: + void setup_angle_to_weight(double start_angle, double end_angle) const; + void setup_weighting_function(Eigen::Vector3d position, double width, + double height) const; + + void calc_pv_spread(Eigen::Vector3d position, double width, double height, + Eigen::Ref out) const; + + mutable PolarExtentCoreContext ctx; + std::unique_ptr core_impl; + std::shared_ptr psp; + + mutable Eigen::VectorXd pvsMin; + mutable Eigen::VectorXd pvsMax; + }; + +} // namespace ear diff --git a/src/object_based/polar_extent_core.hpp b/src/object_based/polar_extent_core.hpp new file mode 100644 index 0000000..90d8090 --- /dev/null +++ b/src/object_based/polar_extent_core.hpp @@ -0,0 +1,66 @@ +#pragma once +#include +#include + +namespace ear { + + /// the type used for internal extent calculations (calculating the weights + /// and summing the gains) + using extent_float_t = float; + + /// data shared between the core and PolarExtent; see PolarExtent + struct PolarExtentCoreContext { + // sizes rounded up to a multiple of the batch size + size_t num_points; + size_t num_speakers; + + size_t real_num_speakers; + + // position of each point + extent_float_t *xs, *ys, *zs; + // panning result for point and speaker is panning_results[speaker * + // num_points + point] + extent_float_t *panning_results; + // total panning results in each batch is summed_panning_results[batch * + // num_speakers + speaker] + extent_float_t *summed_panning_results; + + // use the circular weighting function? + bool is_circular; + + // for determining the distance from the centre + extent_float_t flippedBasis[9]; + extent_float_t circle_test[2]; + extent_float_t right_circle_centre[2]; + + // for mapping from the distance from the centre to a weight + extent_float_t cos_start_angle, cos_end_angle; + extent_float_t sin_start_angle, sin_end_angle; + extent_float_t m, c; + + extent_float_t *results; + }; + + /// interface for polar extent cores; see PolarExtent + class PolarExtentCore { + public: + virtual ~PolarExtentCore(); + + /// the required alignment of arrays in bytes + virtual size_t alignment() const = 0; + /// the number of points processed at once; affects the size and shape of + /// arrays + virtual size_t batch_size() const = 0; + /// calculate the weight for each point, and sum the gains into ctx.results + virtual void run(PolarExtentCoreContext &ctx) const = 0; + }; + + /// get a PolarExtentCore for a given xsimd architecture + template + std::unique_ptr get_polar_extent_core_arch(); + /// get a scalar implementation of PolarExtentCore + std::unique_ptr get_polar_extent_core_scalar(); + + /// get the default implementation of PolarExtentCore for this platform + std::unique_ptr get_polar_extent_core(); +} // namespace ear diff --git a/src/object_based/polar_extent_scalar.cpp b/src/object_based/polar_extent_scalar.cpp new file mode 100644 index 0000000..cd5036f --- /dev/null +++ b/src/object_based/polar_extent_scalar.cpp @@ -0,0 +1,115 @@ +#include +#include +#include "polar_extent_core.hpp" + +namespace ear { + + struct ExtentPosition { + extent_float_t x; + extent_float_t y; + extent_float_t z; + + extent_float_t dot(const extent_float_t *vec) const { + return x * vec[0] + y * vec[1] + z * vec[2]; + } + + ExtentPosition transform(const extent_float_t *m) const { + ExtentPosition res; + res.x = dot(m); + res.y = dot(m + 3); + res.z = dot(m + 6); + return res; + } + }; + + /// fallback scalar implementation of PolarExtentCore; see PolarExtent + class PolarExtentCoreScalar : public PolarExtentCore { + public: + virtual size_t alignment() const override { + // required for posix_memalign + return sizeof(void *); + } + virtual size_t batch_size() const override { return 1; } + + extent_float_t weight_from_cos(PolarExtentCoreContext &ctx, + extent_float_t cos_angle) const { + if (cos_angle >= ctx.cos_start_angle) return 1.0; + if (cos_angle <= ctx.cos_end_angle) return 0.0; + + return ctx.m * std::acos(cos_angle) + ctx.c; + } + + extent_float_t weight_from_sin(PolarExtentCoreContext &ctx, + extent_float_t sin_angle) const { + if (sin_angle <= ctx.sin_start_angle) return 1.0; + if (sin_angle >= ctx.sin_end_angle) return 0.0; + + return ctx.m * std::asin(sin_angle) + ctx.c; + } + + extent_float_t weight_circle(PolarExtentCoreContext &ctx, + const ExtentPosition &position) const { + // simplified dot product assuming that circle_centre is {0, 1, 0} + auto dot = position.dot(ctx.flippedBasis + 3); + return weight_from_cos(ctx, dot); + } + + extent_float_t weight_stadium(PolarExtentCoreContext &ctx, + const ExtentPosition &position) const { + ExtentPosition position_t = position.transform(ctx.flippedBasis); + + ExtentPosition position_t_right = position_t; + position_t_right.x = std::abs(position_t.x); + + auto circle_test_dot = position_t_right.x * ctx.circle_test[0] + + position_t_right.y * ctx.circle_test[1]; + bool in_straight_line_part = circle_test_dot >= 0.0; + + if (in_straight_line_part) + return weight_from_sin(ctx, std::abs(position_t.z)); + else { + auto circle_dot = position_t_right.x * ctx.right_circle_centre[0] + + position_t_right.y * ctx.right_circle_centre[1]; + return weight_from_cos(ctx, circle_dot); + } + } + + void run(PolarExtentCoreContext &ctx) const override { + for (size_t speaker_idx = 0; speaker_idx < ctx.num_speakers; + speaker_idx++) + ctx.results[speaker_idx] = 0.0; + + for (size_t i = 0; i < ctx.num_points; i++) { + ExtentPosition pos; + pos.x = ctx.xs[i]; + pos.y = ctx.ys[i]; + pos.z = ctx.zs[i]; + + extent_float_t weight = ctx.is_circular ? weight_circle(ctx, pos) + : weight_stadium(ctx, pos); + + // add to results, using summed_panning_results for both as this is in + // a more efficient order + if (weight == 1.0) { + for (size_t speaker_idx = 0; speaker_idx < ctx.num_speakers; + speaker_idx++) { + extent_float_t gain = + ctx.summed_panning_results[i * ctx.num_speakers + speaker_idx]; + ctx.results[speaker_idx] += gain; + } + } else if (weight != 0) { + for (size_t speaker_idx = 0; speaker_idx < ctx.num_speakers; + speaker_idx++) { + extent_float_t gain = + ctx.summed_panning_results[i * ctx.num_speakers + speaker_idx]; + ctx.results[speaker_idx] += weight * gain; + } + } + } + } + }; + + std::unique_ptr get_polar_extent_core_scalar() { + return boost::make_unique(); + } +} // namespace ear diff --git a/src/object_based/polar_extent_simd.cpp b/src/object_based/polar_extent_simd.cpp new file mode 100644 index 0000000..90a6729 --- /dev/null +++ b/src/object_based/polar_extent_simd.cpp @@ -0,0 +1,23 @@ +#include +#include "../common/helpers/xsimd_extension.hpp" +#include "polar_extent.hpp" + +namespace ear { + using arch_list = xsimd::arch_list; + + struct GetPolarExtentCore { + template + std::unique_ptr operator()(Arch) { + return get_polar_extent_core_arch(); + } + + std::unique_ptr operator()(xsimd::generic_for_dispatch) { + return get_polar_extent_core_scalar(); + } + }; + + std::unique_ptr get_polar_extent_core() { + auto dispatcher = xsimd::dispatch(GetPolarExtentCore{}); + return dispatcher(); + } +} // namespace ear diff --git a/src/object_based/polar_extent_simd.hpp b/src/object_based/polar_extent_simd.hpp new file mode 100644 index 0000000..3ce3af2 --- /dev/null +++ b/src/object_based/polar_extent_simd.hpp @@ -0,0 +1,138 @@ +#pragma once +#include +#include "polar_extent_core.hpp" +#include "xsimd/xsimd.hpp" + +namespace ear { + + template + struct PositionBatch { + T x; + T y; + T z; + + T dot(const extent_float_t *vec) const { + return xsimd::fma(x, T{vec[0]}, xsimd::fma(y, T{vec[1]}, z * T{vec[2]})); + } + + PositionBatch transform(const extent_float_t *m) const { + PositionBatch res; + res.x = dot(m); + res.y = dot(m + 3); + res.z = dot(m + 6); + return res; + } + }; + + /// SIMD implementation of PolarExtentCore for a given architecture; see + /// PolarExtent + /// + /// these are to be instantiated in polar_extent_simd_instance.cpp + template + class PolarExtentCoreSimd : public PolarExtentCore { + public: + using batch = xsimd::batch; + + virtual size_t alignment() const override { return arch::alignment(); } + virtual size_t batch_size() const override { return batch::size; } + + batch weight_from_cos(PolarExtentCoreContext &ctx, batch cos_angle) const { + auto start = cos_angle >= ctx.cos_start_angle; + auto end = cos_angle <= ctx.cos_end_angle; + if (xsimd::all(start)) return batch{1.0}; + if (xsimd::all(end)) return batch{0.0}; + + auto ramp = ctx.m * xsimd::acos(cos_angle) + ctx.c; + return xsimd::select(start, batch{1.0}, + xsimd::select(end, batch{0.0}, ramp)); + } + + batch weight_from_sin(PolarExtentCoreContext &ctx, batch sin_angle) const { + auto start = sin_angle <= ctx.sin_start_angle; + auto end = sin_angle >= ctx.sin_end_angle; + if (xsimd::all(start)) return batch{1.0}; + if (xsimd::all(end)) return batch{0.0}; + + auto ramp = ctx.m * xsimd::asin(sin_angle) + ctx.c; + return xsimd::select(start, batch{1.0}, + xsimd::select(end, batch{0.0}, ramp)); + } + + batch weight_circle(PolarExtentCoreContext &ctx, + const PositionBatch &position) const { + // simplified dot product assuming that circle_centre is {0, 1, 0} + auto dot = position.dot(ctx.flippedBasis + 3); + return weight_from_cos(ctx, dot); + } + + batch weight_stadium(PolarExtentCoreContext &ctx, + const PositionBatch &position) const { + PositionBatch position_t = position.transform(ctx.flippedBasis); + + PositionBatch position_t_right = position_t; + position_t_right.x = xsimd::abs(position_t.x); + + auto circle_test_dot = + xsimd::fma(position_t_right.x, batch{ctx.circle_test[0]}, + position_t_right.y * batch{ctx.circle_test[1]}); + auto in_straight_line_part = circle_test_dot >= batch{0.0}; + + // for the straight line part + auto straight_val = weight_from_sin(ctx, xsimd::abs(position_t.z)); + + // for the circle part + auto circle_dot = + xsimd::fma(position_t_right.x, batch{ctx.right_circle_centre[0]}, + position_t_right.y * batch{ctx.right_circle_centre[1]}); + auto circle_val = weight_from_cos(ctx, circle_dot); + + if (xsimd::all(in_straight_line_part)) + return straight_val; + else if (xsimd::none(in_straight_line_part)) + return circle_val; + else + return xsimd::select(in_straight_line_part, straight_val, circle_val); + } + + void run(PolarExtentCoreContext &ctx) const override { + for (size_t speaker_idx = 0; speaker_idx < ctx.num_speakers; + speaker_idx++) + ctx.results[speaker_idx] = 0.0; + + for (size_t i = 0; i < ctx.num_points; i += batch::size) { + PositionBatch pos; + pos.x = batch::load_aligned(ctx.xs + i); + pos.y = batch::load_aligned(ctx.ys + i); + pos.z = batch::load_aligned(ctx.zs + i); + + batch weight = ctx.is_circular ? weight_circle(ctx, pos) + : weight_stadium(ctx, pos); + + bool all_zeros = xsimd::all(weight == 0.0); + bool all_ones = xsimd::all(weight == 1.0); + if (all_ones) { + // add pre-computed sum for this batch + size_t batch_i = i / batch::size; + for (size_t speaker_idx = 0; speaker_idx < ctx.num_speakers; + speaker_idx += batch::size) { + batch gain = + batch::load_aligned(ctx.summed_panning_results + + (batch_i * ctx.num_speakers) + speaker_idx); + + extent_float_t *result_ptr = ctx.results + speaker_idx; + xsimd::store_aligned(result_ptr, + batch::load_aligned(result_ptr) + gain); + } + } else if (!all_zeros) { + for (size_t speaker_idx = 0; speaker_idx < ctx.real_num_speakers; + speaker_idx++) { + batch gain = batch::load_aligned(ctx.panning_results + + speaker_idx * ctx.num_points + i); + ctx.results[speaker_idx] += xsimd::reduce_add(gain * weight); + } + } + } + } + }; + +} // namespace ear diff --git a/src/object_based/polar_extent_simd_instance.cpp b/src/object_based/polar_extent_simd_instance.cpp new file mode 100644 index 0000000..7c0b34e --- /dev/null +++ b/src/object_based/polar_extent_simd_instance.cpp @@ -0,0 +1,14 @@ +#include +#include "../common/helpers/xsimd_extension.hpp" +#include "polar_extent_simd.hpp" + +namespace ear { + + template + std::unique_ptr get_polar_extent_core_arch() { + return boost::make_unique>(); + } + + template std::unique_ptr + get_polar_extent_core_arch(); +} // namespace ear diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index a77e9c6..57a0059 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -27,6 +27,7 @@ file(COPY "${CMAKE_CURRENT_SOURCE_DIR}/test_data" DESTINATION ${PROJECT_BINARY_D if(EAR_INTERNAL_TESTS) add_ear_test("decorrelate_tests") add_ear_test("extent_tests") + target_sources(extent_tests PRIVATE "reference/extent.cpp") add_ear_test("gain_calculator_direct_speakers_tests") add_ear_test("geom_tests") add_ear_test("hoa_tests") diff --git a/tests/extent_tests.cpp b/tests/extent_tests.cpp index 0acd35d..72bedba 100644 --- a/tests/extent_tests.cpp +++ b/tests/extent_tests.cpp @@ -1,4 +1,5 @@ #include +#include #include #include #include @@ -6,7 +7,9 @@ #include "common/geom.hpp" #include "common/helpers/eigen_helpers.hpp" #include "ear/bs2051.hpp" -#include "object_based/extent.hpp" +#include "eigen_utils.hpp" +#include "object_based/polar_extent.hpp" +#include "reference/extent.hpp" using namespace ear; @@ -49,38 +52,9 @@ TEST_CASE("test_basis") { REQUIRE(calcBasis(cart(90.0, -90.0 + 1e-6, 1.0)).isApprox(expected, eps)); } -TEST_CASE("test_azimuth_elevation_on_basis") { - double eps = 1e-6; - - double azimuth, elevation; - Eigen::Matrix3d basis; - basis = calcBasis(cart(0.0, 10.0, 1.0)); - std::tie(azimuth, elevation) = - azimuthElevationOnBasis(basis, cart(0.0, 10.0, 1.0)); - REQUIRE(azimuth == Approx(0.0).margin(eps)); - REQUIRE(elevation == Approx(0.0).margin(eps)); - std::tie(azimuth, elevation) = - azimuthElevationOnBasis(basis, cart(0.0, 20.0, 1.0)); - REQUIRE(azimuth == Approx(0.0).margin(eps)); - REQUIRE(elevation == Approx(radians(10.0)).margin(eps)); - basis = calcBasis(cart(-10.0, 0.0, 1.0)); - std::tie(azimuth, elevation) = - azimuthElevationOnBasis(basis, cart(-20.0, 0.0, 1.0)); - REQUIRE(azimuth == Approx(radians(10.0)).margin(eps)); - REQUIRE(elevation == Approx(0.0).margin(eps)); -} - -TEST_CASE("test_cart_on_basis") { - Eigen::Matrix3d basis; - basis = calcBasis(cart(0.0, 10.0, 1.0)); - REQUIRE( - cartOnBasis(basis, 0.0, radians(10.0)).isApprox(cart(0.0, 20.0, 1.0))); - basis = calcBasis(cart(-10.0, 0.0, 1.0)); - REQUIRE( - cartOnBasis(basis, radians(10.0), 0.0).isApprox(cart(-20.0, 0.0, 1.0))); -} - TEST_CASE("test_weight_func") { + // this isn't exposed by the real implementation; test the reference version, + // which is checked against the real implementation in same_as_reference double fade = 10.0; double height = 10.0; @@ -100,13 +74,15 @@ TEST_CASE("test_weight_func") { Eigen::Vector4d{0, 1, 1, 0}); point = cart(azimuth, elevation, 1.0); - WeightingFunction weightingFunc(cart(0.0, 0.0, 1.0), width, height); - actual = weightingFunc(point); + reference::WeightingFunction weightingFunc(cart(0.0, 0.0, 1.0), width, + height); + actual = weightingFunc(point.transpose()); REQUIRE(actual == Approx(expected)); // Swapped point = cart(elevation, azimuth, 1.0); - WeightingFunction weightingFuncSwap(cart(0.0, 0.0, 1.0), height, width); - actual = weightingFuncSwap(point); + reference::WeightingFunction weightingFuncSwap(cart(0.0, 0.0, 1.0), + height, width); + actual = weightingFuncSwap(point.transpose()); REQUIRE(actual == Approx(expected)); } } @@ -118,13 +94,15 @@ TEST_CASE("test_weight_func") { Eigen::Vector4d{0, 1, 1, 0}); point = cart(azimuth, 0.0, 1.0); - WeightingFunction weightingFunc(cart(0.0, 0.0, 1.0), width, height); - actual = weightingFunc(point); + reference::WeightingFunction weightingFunc(cart(0.0, 0.0, 1.0), width, + height); + actual = weightingFunc(point.transpose()); REQUIRE(actual == Approx(expected)); // Swapped point = cart(0.0, azimuth, 1.0); - WeightingFunction weightingFuncSwap(cart(0.0, 0.0, 1.0), height, width); - actual = weightingFuncSwap(point); + reference::WeightingFunction weightingFuncSwap(cart(0.0, 0.0, 1.0), height, + width); + actual = weightingFuncSwap(point.transpose()); REQUIRE(actual == Approx(expected)); } } @@ -132,24 +110,95 @@ TEST_CASE("test_weight_func") { TEST_CASE("test_pv") { Layout layout = getLayout("9+10+3").withoutLfe(); std::shared_ptr psp = configurePolarPanner(layout); - PolarExtentPanner extentPanner(psp); + PolarExtent extentPanner(psp); + Eigen::VectorXd pvs{psp->numberOfOutputChannels()}; - REQUIRE(extentPanner.calcPvSpread(cart(0.0, 0.0, 1.0), 0.0, 0.0) == - psp->handle(cart(0.0, 0.0, 1.0)).get()); - REQUIRE(extentPanner.calcPvSpread(cart(10.0, 20.0, 1.0), 0.0, 0.0) == - psp->handle(cart(10.0, 20.0, 1.0)).get()); + extentPanner.handle(cart(0.0, 0.0, 1.0), 0.0, 0.0, 0.0, pvs); + REQUIRE(pvs == psp->handle(cart(0.0, 0.0, 1.0)).get()); + extentPanner.handle(cart(10.0, 20.0, 1.0), 0.0, 0.0, 0.0, pvs); + REQUIRE(pvs == psp->handle(cart(10.0, 20.0, 1.0)).get()); std::vector> positions{ - std::make_pair(cart(0.0, 0.0, 1.0), 1e-10), + std::make_pair(cart(0.0, 0.0, 1.0), 1e-5), std::make_pair(cart(30.0, 10.0, 1.0), 1e-2)}; for (auto position : positions) { Eigen::Vector3d pos = position.first; double tol = position.second; - Eigen::VectorXd spread_pv = extentPanner.calcPvSpread(pos, 20.0, 10.0); - REQUIRE(spread_pv.norm() == Approx(1.0)); + extentPanner.handle(pos, 20.0, 10.0, 0.0, pvs); + REQUIRE(pvs.norm() == Approx(1.0)); Eigen::VectorXd vv = - spread_pv.transpose() * toPositionsMatrix(layout.positions()); + pvs.transpose() * toPositionsMatrix(layout.positions()); vv /= vv.norm(); REQUIRE(vv.isApprox(pos, tol)); } } + +TEST_CASE("same_as_reference") { + Layout layout = getLayout("9+10+3").withoutLfe(); + std::shared_ptr psp = configurePolarPanner(layout); + + PolarExtent extentPanner(psp); + PolarExtent extentPannerScalar(psp, get_polar_extent_core_scalar()); + Eigen::VectorXd spread_pv{psp->numberOfOutputChannels()}; + + reference::PolarExtentPanner extentPannerReference(psp); + + for (size_t i = 0; i < 1000; i++) { + Eigen::Vector3d pos = Eigen::Vector3d::Random(); + pos.normalize(); + Eigen::Vector2d size = (Eigen::Vector2d::Random().array() + 1) * 180; + double width = size(0); + double height = size(1); + Eigen::VectorXd spread_pv_reference = + extentPannerReference.handle(pos, width, height, 0.0); + + INFO("size " << width << " " << height); + INFO("spread_pv_reference " << spread_pv_reference.transpose()); + + extentPanner.handle(pos, width, height, 0.0, spread_pv); + INFO("spread_pv " << spread_pv.transpose()); + CHECK(spread_pv.isApprox(spread_pv_reference, 1e-5)); + + extentPannerScalar.handle(pos, width, height, 0.0, spread_pv); + INFO("spread_pv scalar " << spread_pv.transpose()); + CHECK(spread_pv.isApprox(spread_pv_reference, 1e-5)); + } +} + +#ifdef CATCH_CONFIG_ENABLE_BENCHMARKING +TEST_CASE("bench_panner", "[.benchmark]") { + Layout layout = getLayout("9+10+3").withoutLfe(); + std::shared_ptr psp = configurePolarPanner(layout); + reference::PolarExtentPanner extentPannerRef(psp); + + PolarExtent extentPanner(psp); + PolarExtent extentPannerScalar(psp, get_polar_extent_core_scalar()); + Eigen::VectorXd pvs{psp->numberOfOutputChannels()}; + + Eigen::Vector3d pos(1, 1, 1); + pos.normalize(); + + std::vector> sizes = { + {5, 5}, {30, 30}, {30, 50}, {90, 90}, {135, 135}, {360, 360}, + }; + + for (auto &width_height : sizes) { + double width, height; + std::tie(width, height) = width_height; + std::string size_name = + std::to_string((int)width) + " " + std::to_string((int)height); + + BENCHMARK("ref " + size_name) { + return extentPannerRef.handle(pos, width, height, 0.0); + }; + + BENCHMARK("impl scalar " + size_name) { + return extentPannerScalar.handle(pos, width, height, 0.0, pvs); + }; + + BENCHMARK("impl " + size_name) { + return extentPanner.handle(pos, width, height, 0.0, pvs); + }; + } +} +#endif diff --git a/tests/reference/extent.cpp b/tests/reference/extent.cpp new file mode 100644 index 0000000..ee66223 --- /dev/null +++ b/tests/reference/extent.cpp @@ -0,0 +1,246 @@ +#include "extent.hpp" + +#include +#include +#include +#include +#include +#include "common/geom.hpp" +#include "common/helpers/eigen_helpers.hpp" + +const double PI = boost::math::constants::pi(); + +namespace ear { + namespace reference { + + /** @brief Normalise position or return {0,1,0}. + * + * @param position Position to normalise. + * + * @returns normalised position + */ + Eigen::Vector3d safeNormPosition(Eigen::Vector3d position) { + double norm = position.norm(); + if (norm < 1e-10) { + return Eigen::Vector3d{0.0, 1.0, 0.0}; + } else { + return position / norm; + } + } + + double extentMod(double extent, double distance) { + double minSize = 0.2; + double size = interp(extent, Eigen::Vector2d(0.0, 360.0), + Eigen::Vector2d(minSize, 1.0)); + double extent1 = 4.0 * degrees(atan2(size, 1.0)); + return interp(4.0 * degrees(atan2(size, distance)), + Eigen::Vector3d(0.0, extent1, 360.0), + Eigen::Vector3d(0.0, extent, 360.0)); + } + + Eigen::Matrix3d calcBasis(Eigen::Vector3d position) { + position = safeNormPosition(position); + double az = azimuth(position); + double el = elevation(position); + + // points near the poles have indeterminate azimuth; assume 0 + if (abs(el) > (90.0 - 1e-5)) { + az = 0.0; + } + return localCoordinateSystem(az, el); + } + + Eigen::Vector3d cartOnBasis(Eigen::Matrix3d basis, double azimuth, + double elevation) { + Eigen::RowVector3d cartPosRel{sin(azimuth) * cos(elevation), + cos(azimuth) * cos(elevation), + sin(elevation)}; + return cartPosRel * basis; + } + + std::pair azimuthElevationOnBasis( + Eigen::Matrix3d basis, Eigen::RowVector3d position) { + // project onto each basis, and clip components to keep asin happy + Eigen::Vector3d components = + (position * basis.transpose()).cwiseMin(1.0).cwiseMax(-1.0); + + double azimuth = atan2(components(0), components(1)); + double elevation = asin(components(2)); + + return std::make_pair(azimuth, elevation); + } + + WeightingFunction::WeightingFunction(Eigen::Vector3d position, double width, + double height) { + _width = radians(width) / 2; + _height = radians(height) / 2; + + // basis vectors to rotate the vsource positions towards position + Eigen::Matrix3d basises = calcBasis(position); + + _circleRadius = std::min(_width, _height); + + // Flip the width and the height such that it is always wider than it is + // high from here in. + if (_height > _width) { + std::swap(_height, _width); + _flippedBasis = basises.colwise().reverse(); + } else { + _flippedBasis = basises; + } + + // modify the width to make it meet at the back. + double widthFull = PI + _height; + // interpolate to this from a width of pi/2 to pi + double widthMod = interp(_width, Eigen::Vector3d{0.0, PI / 2.0, PI}, + Eigen::Vector3d{0.0, PI / 2.0, widthFull}); + // apply this fully for a height of less than pi/4; tail off until pi/2 + _width = interp(_height, Eigen::Vector4d{0, PI / 4.0, PI / 2.0, PI}, // + Eigen::Vector4d{widthMod, widthMod, _width, _width}); // + + // angle of the circle centres from the source position; width is to the + // end of the rectangle. + _circlePos = _width - _circleRadius; + + // Cartesian circle centres + _circlePositions << cartOnBasis(_flippedBasis, -_circlePos, 0.0), + cartOnBasis(_flippedBasis, _circlePos, 0.0); + } + + double WeightingFunction::operator()(Eigen::Vector3d position) const { + // Flipped azimuths and elevations; the straight edges are always along + // azimuth lines. + double azimuth, elevation; + std::tie(azimuth, elevation) = + azimuthElevationOnBasis(_flippedBasis, position); + + // The distance is the angle away from the defined shape; 0 or negative is + // inside. + double distance = 0.0; + + // for the straight lines + if (abs(azimuth) <= _circlePos) { + distance = abs(elevation) - _circleRadius; + } else { + // distance from the closest circle centre + size_t nearest_circle = azimuth < 0 ? 0 : 1; + double angle = + position.transpose() * _circlePositions.col(nearest_circle); + double circleDistance = acos(boost::algorithm::clamp(angle, -1.0, 1.0)); + distance = circleDistance - _circleRadius; + } + // fade the weight from one to zero over fadeWidth + return interp(distance, Eigen::Vector2d{0.0, radians(_fadeWidth)}, + Eigen::Vector2d{1.0, 0.0}); + } + + SpreadingPanner::SpreadingPanner(std::shared_ptr psp, + int nRows) + : _psp(psp), _nRows(nRows) { + _panningPositions = _generatePanningPositionsEven(); + _panningPositionsResults = _generatePanningPositionsResults(); + } + + Eigen::VectorXd SpreadingPanner::panningValuesForWeight( + const WeightingFunction& weightFunc) { + Eigen::VectorXd weights(_panningPositions.rows()); + for (int i = 0; i < _panningPositions.rows(); ++i) { + weights(i) = weightFunc(_panningPositions.row(i)); + } + Eigen::VectorXd totalPv = weights.transpose() * _panningPositionsResults; + return totalPv / totalPv.norm(); + } + + Eigen::MatrixXd SpreadingPanner::_generatePanningPositionsEven() { + Eigen::VectorXd elevations = + Eigen::VectorXd::LinSpaced(_nRows, -90.0, 90.0); + Eigen::MatrixXd positions(0, 3); + + for (double el : elevations) { + double radius = cos(radians(el)); + double perimiter = 2 * PI * radius; + double perimiter_centre = 2 * PI; + + int nPoints = static_cast( + std::round((perimiter / perimiter_centre) * 2 * (_nRows - 1))); + if (nPoints == 0) { + nPoints = 1; + } + Eigen::VectorXd azimuths = + Eigen::VectorXd::LinSpaced(nPoints + 1, 0.0, 360.0); + for (int i = 0; i < azimuths.size() - 1; ++i) { + double az = azimuths(i); + positions.conservativeResize(positions.rows() + 1, Eigen::NoChange); + positions.row(positions.rows() - 1) = cart(az, el, 1.0); + } + } + return positions; + } + + Eigen::MatrixXd SpreadingPanner::_generatePanningPositionsResults() { + Eigen::MatrixXd results(_panningPositions.rows(), + _psp->numberOfOutputChannels()); + for (int i = 0; i < _panningPositions.rows(); ++i) { + results.row(i) = _psp->handle(_panningPositions.row(i)).get(); + } + return results; + } + + PolarExtentPanner::PolarExtentPanner(std::shared_ptr psp) + : _psp(psp), _spreadingPanner(SpreadingPanner(psp, _nRows)){}; + + Eigen::VectorXd PolarExtentPanner::calcPvSpread(Eigen::Vector3d position, + double width, + double height) { + // When calculating the spread panning values the width and height are + // set to at least fade_width. For sizes where any of the dimensions is + // less than this, interpolate linearly between the point and spread + // panning values. + double ammount_spread = + interp(std::max(width, height), Eigen::Vector2d(0.0, _fadeWidth), + Eigen::Vector2d(0.0, 1.0)); + double ammountPoint = 1.0 - ammount_spread; + Eigen::ArrayXd pv = Eigen::ArrayXd::Zero(_psp->numberOfOutputChannels()); + if (ammountPoint > 1e-10) { + pv += ammountPoint * _psp->handle(position).get().array().square(); + } + if (ammount_spread > 1e-10) { + // minimum width and height as above + width = std::max(width, _fadeWidth / 2.0); + height = std::max(height, _fadeWidth / 2.0); + + WeightingFunction weightingFunction(position, width, height); + Eigen::VectorXd panning_values = + _spreadingPanner.panningValuesForWeight(weightingFunction); + pv += ammount_spread * panning_values.array().square(); + } + return pv.sqrt().matrix(); + } + + Eigen::VectorXd PolarExtentPanner::handle(Eigen::Vector3d position, + double width, double height, + double depth) { + double distance = position.norm(); + + if (depth != 0.0) { + double distanceMin = distance - depth / 2.0; + double distanceMax = distance + depth / 2.0; + distanceMin = (distanceMin < 0) ? 0.0 : distanceMin; + distanceMax = (distanceMax < 0) ? 0.0 : distanceMax; + Eigen::VectorXd pvsMin = + calcPvSpread(position, extentMod(width, distanceMin), + extentMod(height, distanceMin)); + Eigen::VectorXd pvsMax = + calcPvSpread(position, extentMod(width, distanceMax), + extentMod(height, distanceMax)); + return ((pvsMin.array().square() + pvsMax.array().square()) / 2.0) + .sqrt(); + } else { + Eigen::VectorXd pvs = calcPvSpread(position, extentMod(width, distance), + extentMod(height, distance)); + return pvs; + } + } + + } // namespace reference +} // namespace ear diff --git a/tests/reference/extent.hpp b/tests/reference/extent.hpp new file mode 100644 index 0000000..7dbb282 --- /dev/null +++ b/tests/reference/extent.hpp @@ -0,0 +1,151 @@ +#pragma once +#include +#include +#include "common/point_source_panner.hpp" + +namespace ear { + namespace reference { + + class ExtentPanner {}; + + /** @brief Modify an extent parameter given a distance. + * + * A right triangle if formed, with the adjacent edge being the distance, + * and the opposite edge being determined from the extent. The angle formed + * is then used to determine the new extent. + * + * - at distance=0, the extent is always 360 + * - at distance=1, the original extent is used + * - at distance>1, the extent decreases + * - in 0 < distance < 1, the extent changes more steeply around 0 for + * smaller extents + */ + double extentMod(double extent, double distance); + + /** @brief Calculate basis vectors that rotate (0, 1, 0) + onto source_pos. */ + Eigen::Matrix3d calcBasis(Eigen::Vector3d position); + + /** @brief Polar to Cartesian in radians with no distance, in a given + * basis.*/ + Eigen::Vector3d cartOnBasis(Eigen::Matrix3d basis, double azimuth, + double elevation); + + std::pair azimuthElevationOnBasis( + Eigen::Matrix3d basis, Eigen::RowVector3d position); + + /** @brief Weighting function for spread sources. + * + * The weighting function is one inside a region approximately determined by + * a width x height rectangle in azimuth-elevation space, with + * maximally-sized rounded corners; the shape of the corners is calculated + * using the vector angle from their centres (always directly above or below + * the source position) so as to avoid issues at the poles. + * + * The two straight edges of the rectangle are always parallel in Cartesian + * space; this is achieved by following azimuth lines; for tall sources, the + * whole coordinate system is rotated 90 degrees about the source position + * to achieve this. + * + * Note that for sources where width == height, this degrades to a circular + * region relative to the source position. + * + * To make the two ends meet, the width is adjusted such that a width of 180 + * degrees is mapped to width + height. + * + */ + class WeightingFunction { + public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + /** @brief Ctor + * + * @param position Centre of the extent. + * @param width Width of the extent in degrees from one edge to the other. + * @param height Height of the extent in degrees from one edge to the + * other. + */ + WeightingFunction(Eigen::Vector3d position, double width, double height); + + /** @brief Calculate weight for position */ + double operator()(Eigen::Vector3d position) const; + + private: + const double _fadeWidth = 10.0; + double _width; + double _height; + double _circleRadius; + Eigen::Matrix3d _flippedBasis; + double _circlePos; + Eigen::Matrix _circlePositions; + }; + + class SpreadingPanner { + public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + SpreadingPanner(std::shared_ptr psp, int nRows); + + /** @brief Panning values for a given weighting function. + * + * @param weightFunc function from Cartesian position to weight in range + * (0, 1) + * + * @return panning value for each speaker. + */ + Eigen::VectorXd panningValuesForWeight( + const WeightingFunction& weightFunc); + + private: + /** @brief Generate points spread evenly on the sphere. + * + * Based on + * http://web.archive.org/web/20150108040043/http://www.math.niu.edu/~rusin/known-math/95/equispace.elect + * + * @param nRows number of rows to place on sphere, e.g. 37 for 5 degree + * spacing + * + * @returns cartesian array. + */ + Eigen::MatrixXd _generatePanningPositionsEven(); + Eigen::MatrixXd _generatePanningPositionsResults(); + + std::shared_ptr _psp; + int _nRows; + Eigen::MatrixXd _panningPositions; + Eigen::MatrixXd _panningPositionsResults; + }; + + class PolarExtentPanner : public ExtentPanner { + public: + PolarExtentPanner(std::shared_ptr psp); + + /** @brief Calculate loudspeaker gains given position and extent + * parameters. + * + * @param position Cartesian source position + * @param width block format width parameter + * @param height block format height parameter + * @param depth block format depth parameter + * + * @returns loudspeaker gains for each channel + */ + Eigen::VectorXd handle(Eigen::Vector3d position, double width, + double height, double depth); + + /** @brief Calculate the speaker panning values for the position, width, + * and height of a source; this just deals with the positioning and + * spreading. + */ + Eigen::VectorXd calcPvSpread(Eigen::Vector3d position, double width, + double height); + + private: + const int _nRows = 37; // 5 degrees per row + const double _fadeWidth = 10.0; // degrees + std::shared_ptr _psp; + SpreadingPanner _spreadingPanner; + }; + + class CartesianExtentPanner : public ExtentPanner {}; + + } // namespace reference +} // namespace ear