diff options
47 files changed, 4756 insertions, 437 deletions
diff --git a/Makefile.am b/Makefile.am index d29b530..47d6f43 100644 --- a/Makefile.am +++ b/Makefile.am @@ -35,10 +35,11 @@ endif bin_PROGRAMS = odr-dabmod -odr_dabmod_CFLAGS = -Wall -Isrc -Ilib \ - $(GITVERSION_FLAGS) -odr_dabmod_CXXFLAGS = -Wall -Isrc -Ilib \ - $(GITVERSION_FLAGS) $(BOOST_CPPFLAGS) +KISS_FLAGS=-DFIXED_POINT=32 +odr_dabmod_CFLAGS = -Wall -Isrc -Ilib -Ikiss \ + $(GITVERSION_FLAGS) $(KISS_FLAGS) +odr_dabmod_CXXFLAGS = -Wall -Isrc -Ilib -Ikiss \ + $(GITVERSION_FLAGS) $(BOOST_CPPFLAGS) $(KISS_FLAGS) odr_dabmod_LDADD = $(BOOST_LDFLAGS) $(BOOST_THREAD_LIB) $(UHD_LIBS) $(LIMESDR_LIBS) $(ADDITIONAL_UHD_LIBS) odr_dabmod_SOURCES = src/DabMod.cpp \ src/PcDebug.h \ @@ -175,7 +176,17 @@ odr_dabmod_SOURCES = src/DabMod.cpp \ src/PAPRStats.cpp \ src/PAPRStats.h \ src/TII.cpp \ - src/TII.h + src/TII.h \ + kiss/kfc.h \ + kiss/kfc.c \ + kiss/kiss_fft.c \ + kiss/kiss_fft.h \ + kiss/kiss_fftnd.c \ + kiss/kiss_fftnd.h \ + kiss/kiss_fftndr.c \ + kiss/kiss_fftndr.h \ + kiss/kiss_fftr.c \ + kiss/kiss_fftr.h man_MANS = man/odr-dabmod.1 diff --git a/doc/example.ini b/doc/example.ini index eda50a5..0d0f8e3 100644 --- a/doc/example.ini +++ b/doc/example.ini @@ -103,6 +103,8 @@ gainmode=var ; If not defined, use Transmission Mode 1 ;mode=1 +fixed_point=1 + ; The digital gain is a value that is multiplied to each sample. It is used ; to tune the chain to make sure that no non-linearities appear up to the ; USRP daughterboard programmable gain amplifier (PGA). diff --git a/fpm/LICENSE b/fpm/LICENSE new file mode 100644 index 0000000..bb86b71 --- /dev/null +++ b/fpm/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2019 Mike Lankamp + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/fpm/README.md b/fpm/README.md new file mode 100644 index 0000000..38ee444 --- /dev/null +++ b/fpm/README.md @@ -0,0 +1,48 @@ +# fpm +A C++ header-only fixed-point math library. "fpm" stands for "fixed-point math". + +It is designed to serve as a drop-in replacement for floating-point types and aims to provide as much of the standard library's functionality as possible with exclusively integers. `fpm` requires C++11 or higher. + +[](https://travis-ci.org/MikeLankamp/fpm) +[](https://ci.appveyor.com/project/MikeLankamp/fpm) + +`fpm` is designed to guard against accidental conversion to and from floats and supports many of the standard C++ maths functions, including trigonometry, power and logarithmic functions, with performance and accuracy generally comparable to alternative libraries. + +## Why use fixed-point math? +There are several reasons why you can not or choose not to use floating-point math, but still want a similar type: +* Your target platform lacks an FPU, does not support floating-point operations or its floating-point operations are + considerably slower than fixed-point integer operations. +* You require deterministic calculations. + +If any of these reasons apply for you, and your problem domain has a clearly outlined range and required resolution, +then fixed-point numbers might be a solution for you. + +## Quick Start +To use `fpm`, include its header `<fpm/fixed.hpp>` and use the `fpm::fixed_16_16`, `fpm::fixed_24_8` or `fpm::fixed_8_24` +types as if they were native floating-pointer types: +```c++ +#include <fpm/fixed.hpp> // For fpm::fixed_16_16 +#include <fpm/math.hpp> // For fpm::cos +#include <fpm/ios.hpp> // For fpm::operator<< +#include <iostream> // For std::cin, std::cout + +int main() { + std::cout << "Please input a number: "; + fpm::fixed_16_16 x; + std::cin >> x; + std::cout << "The cosine of " << x << " radians is: " << cos(x) << std::endl; + return 0; +} +``` + +To use the fixed-point equivalents of the `<math.h>` functions such as `sqrt`, `sin` and `log`, include the header `<fpm/math.hpp>`. +To stream fixed-point values to or from streams, include the header `<fpm/ios.hpp>`. + +## Documentation +Please refer to the [documentation](docs/index.md) for detailed information how to use `fpm`, or skip straight to the [performance](docs/performance.md) or [accuracy](docs/accuracy.md) results. + +## Contributions +This library is a work-in-progress. We welcome any contributions that improve the functional coverage or the performance or accuracy of the mathematical functions. + +## License +See the [LICENSE](LICENSE) file diff --git a/fpm/fixed.hpp b/fpm/fixed.hpp new file mode 100644 index 0000000..e2e71bf --- /dev/null +++ b/fpm/fixed.hpp @@ -0,0 +1,490 @@ +#ifndef FPM_FIXED_HPP +#define FPM_FIXED_HPP + +#include <cassert> +#include <cmath> +#include <cstdint> +#include <functional> +#include <limits> +#include <type_traits> + +namespace fpm +{ + +//! Fixed-point number type +//! \tparam BaseType the base integer type used to store the fixed-point number. This can be a signed or unsigned type. +//! \tparam IntermediateType the integer type used to store intermediate results during calculations. +//! \tparam FractionBits the number of bits of the BaseType used to store the fraction +//! \tparam EnableRounding enable rounding of LSB for multiplication, division, and type conversion +template <typename BaseType, typename IntermediateType, unsigned int FractionBits, bool EnableRounding = true> +class fixed +{ + static_assert(std::is_integral<BaseType>::value, "BaseType must be an integral type"); + static_assert(FractionBits > 0, "FractionBits must be greater than zero"); + static_assert(FractionBits <= sizeof(BaseType) * 8 - 1, "BaseType must at least be able to contain entire fraction, with space for at least one integral bit"); + static_assert(sizeof(IntermediateType) > sizeof(BaseType), "IntermediateType must be larger than BaseType"); + static_assert(std::is_signed<IntermediateType>::value == std::is_signed<BaseType>::value, "IntermediateType must have same signedness as BaseType"); + + // Although this value fits in the BaseType in terms of bits, if there's only one integral bit, this value + // is incorrect (flips from positive to negative), so we must extend the size to IntermediateType. + static constexpr IntermediateType FRACTION_MULT = IntermediateType(1) << FractionBits; + + struct raw_construct_tag {}; + constexpr inline fixed(BaseType val, raw_construct_tag) noexcept : m_value(val) {} + +public: + inline fixed() noexcept = default; + + // Converts an integral number to the fixed-point type. + // Like static_cast, this truncates bits that don't fit. + template <typename T, typename std::enable_if<std::is_integral<T>::value>::type* = nullptr> + constexpr inline explicit fixed(T val) noexcept + : m_value(static_cast<BaseType>(val * FRACTION_MULT)) + {} + + // Converts an floating-point number to the fixed-point type. + // Like static_cast, this truncates bits that don't fit. + template <typename T, typename std::enable_if<std::is_floating_point<T>::value>::type* = nullptr> + constexpr inline explicit fixed(T val) noexcept + : m_value(static_cast<BaseType>((EnableRounding) ? + (val >= 0.0) ? (val * FRACTION_MULT + T{0.5}) : (val * FRACTION_MULT - T{0.5}) + : (val * FRACTION_MULT))) + {} + + // Constructs from another fixed-point type with possibly different underlying representation. + // Like static_cast, this truncates bits that don't fit. + template <typename B, typename I, unsigned int F, bool R> + constexpr inline explicit fixed(fixed<B,I,F,R> val) noexcept + : m_value(from_fixed_point<F>(val.raw_value()).raw_value()) + {} + + // Explicit conversion to a floating-point type + template <typename T, typename std::enable_if<std::is_floating_point<T>::value>::type* = nullptr> + constexpr inline explicit operator T() const noexcept + { + return static_cast<T>(m_value) / FRACTION_MULT; + } + + // Explicit conversion to an integral type + template <typename T, typename std::enable_if<std::is_integral<T>::value>::type* = nullptr> + constexpr inline explicit operator T() const noexcept + { + return static_cast<T>(m_value / FRACTION_MULT); + } + + // Returns the raw underlying value of this type. + // Do not use this unless you know what you're doing. + constexpr inline BaseType raw_value() const noexcept + { + return m_value; + } + + //! Constructs a fixed-point number from another fixed-point number. + //! \tparam NumFractionBits the number of bits used by the fraction in \a value. + //! \param value the integer fixed-point number + template <unsigned int NumFractionBits, typename T, typename std::enable_if<(NumFractionBits > FractionBits)>::type* = nullptr> + static constexpr inline fixed from_fixed_point(T value) noexcept + { + // To correctly round the last bit in the result, we need one more bit of information. + // We do this by multiplying by two before dividing and adding the LSB to the real result. + return (EnableRounding) ? fixed(static_cast<BaseType>( + value / (T(1) << (NumFractionBits - FractionBits)) + + (value / (T(1) << (NumFractionBits - FractionBits - 1)) % 2)), + raw_construct_tag{}) : + fixed(static_cast<BaseType>(value / (T(1) << (NumFractionBits - FractionBits))), + raw_construct_tag{}); + } + + template <unsigned int NumFractionBits, typename T, typename std::enable_if<(NumFractionBits <= FractionBits)>::type* = nullptr> + static constexpr inline fixed from_fixed_point(T value) noexcept + { + return fixed(static_cast<BaseType>( + value * (T(1) << (FractionBits - NumFractionBits))), + raw_construct_tag{}); + } + + // Constructs a fixed-point number from its raw underlying value. + // Do not use this unless you know what you're doing. + static constexpr inline fixed from_raw_value(BaseType value) noexcept + { + return fixed(value, raw_construct_tag{}); + } + + // + // Constants + // + static constexpr fixed e() { return from_fixed_point<61>(6267931151224907085ll); } + static constexpr fixed pi() { return from_fixed_point<61>(7244019458077122842ll); } + static constexpr fixed half_pi() { return from_fixed_point<62>(7244019458077122842ll); } + static constexpr fixed two_pi() { return from_fixed_point<60>(7244019458077122842ll); } + + // + // Arithmetic member operators + // + + constexpr inline fixed operator-() const noexcept + { + return fixed::from_raw_value(-m_value); + } + + inline fixed& operator+=(const fixed& y) noexcept + { + m_value += y.m_value; + return *this; + } + + template <typename I, typename std::enable_if<std::is_integral<I>::value>::type* = nullptr> + inline fixed& operator+=(I y) noexcept + { + m_value += y * FRACTION_MULT; + return *this; + } + + inline fixed& operator-=(const fixed& y) noexcept + { + m_value -= y.m_value; + return *this; + } + + template <typename I, typename std::enable_if<std::is_integral<I>::value>::type* = nullptr> + inline fixed& operator-=(I y) noexcept + { + m_value -= y * FRACTION_MULT; + return *this; + } + + inline fixed& operator*=(const fixed& y) noexcept + { + if (EnableRounding){ + // Normal fixed-point multiplication is: x * y / 2**FractionBits. + // To correctly round the last bit in the result, we need one more bit of information. + // We do this by multiplying by two before dividing and adding the LSB to the real result. + auto value = (static_cast<IntermediateType>(m_value) * y.m_value) / (FRACTION_MULT / 2); + m_value = static_cast<BaseType>((value / 2) + (value % 2)); + } else { + auto value = (static_cast<IntermediateType>(m_value) * y.m_value) / FRACTION_MULT; + m_value = static_cast<BaseType>(value); + } + return *this; + } + + template <typename I, typename std::enable_if<std::is_integral<I>::value>::type* = nullptr> + inline fixed& operator*=(I y) noexcept + { + m_value *= y; + return *this; + } + + inline fixed& operator/=(const fixed& y) noexcept + { + assert(y.m_value != 0); + if (EnableRounding){ + // Normal fixed-point division is: x * 2**FractionBits / y. + // To correctly round the last bit in the result, we need one more bit of information. + // We do this by multiplying by two before dividing and adding the LSB to the real result. + auto value = (static_cast<IntermediateType>(m_value) * FRACTION_MULT * 2) / y.m_value; + m_value = static_cast<BaseType>((value / 2) + (value % 2)); + } else { + auto value = (static_cast<IntermediateType>(m_value) * FRACTION_MULT) / y.m_value; + m_value = static_cast<BaseType>(value); + } + return *this; + } + + template <typename I, typename std::enable_if<std::is_integral<I>::value>::type* = nullptr> + inline fixed& operator/=(I y) noexcept + { + m_value /= y; + return *this; + } + +private: + BaseType m_value; +}; + +// +// Convenience typedefs +// + +using fixed_16_16 = fixed<std::int32_t, std::int64_t, 16>; +using fixed_24_8 = fixed<std::int32_t, std::int64_t, 8>; +using fixed_8_24 = fixed<std::int32_t, std::int64_t, 24>; + +// +// Addition +// + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline fixed<B, I, F, R> operator+(const fixed<B, I, F, R>& x, const fixed<B, I, F, R>& y) noexcept +{ + return fixed<B, I, F, R>(x) += y; +} + +template <typename B, typename I, unsigned int F, bool R, typename T, typename std::enable_if<std::is_integral<T>::value>::type* = nullptr> +constexpr inline fixed<B, I, F, R> operator+(const fixed<B, I, F, R>& x, T y) noexcept +{ + return fixed<B, I, F, R>(x) += y; +} + +template <typename B, typename I, unsigned int F, bool R, typename T, typename std::enable_if<std::is_integral<T>::value>::type* = nullptr> +constexpr inline fixed<B, I, F, R> operator+(T x, const fixed<B, I, F, R>& y) noexcept +{ + return fixed<B, I, F, R>(y) += x; +} + +// +// Subtraction +// + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline fixed<B, I, F, R> operator-(const fixed<B, I, F, R>& x, const fixed<B, I, F, R>& y) noexcept +{ + return fixed<B, I, F, R>(x) -= y; +} + +template <typename B, typename I, unsigned int F, bool R, typename T, typename std::enable_if<std::is_integral<T>::value>::type* = nullptr> +constexpr inline fixed<B, I, F, R> operator-(const fixed<B, I, F, R>& x, T y) noexcept +{ + return fixed<B, I, F, R>(x) -= y; +} + +template <typename B, typename I, unsigned int F, bool R, typename T, typename std::enable_if<std::is_integral<T>::value>::type* = nullptr> +constexpr inline fixed<B, I, F, R> operator-(T x, const fixed<B, I, F, R>& y) noexcept +{ + return fixed<B, I, F, R>(x) -= y; +} + +// +// Multiplication +// + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline fixed<B, I, F, R> operator*(const fixed<B, I, F, R>& x, const fixed<B, I, F, R>& y) noexcept +{ + return fixed<B, I, F, R>(x) *= y; +} + +template <typename B, typename I, unsigned int F, bool R, typename T, typename std::enable_if<std::is_integral<T>::value>::type* = nullptr> +constexpr inline fixed<B, I, F, R> operator*(const fixed<B, I, F, R>& x, T y) noexcept +{ + return fixed<B, I, F, R>(x) *= y; +} + +template <typename B, typename I, unsigned int F, bool R, typename T, typename std::enable_if<std::is_integral<T>::value>::type* = nullptr> +constexpr inline fixed<B, I, F, R> operator*(T x, const fixed<B, I, F, R>& y) noexcept +{ + return fixed<B, I, F, R>(y) *= x; +} + +// +// Division +// + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline fixed<B, I, F, R> operator/(const fixed<B, I, F, R>& x, const fixed<B, I, F, R>& y) noexcept +{ + return fixed<B, I, F, R>(x) /= y; +} + +template <typename B, typename I, unsigned int F, typename T, bool R, typename std::enable_if<std::is_integral<T>::value>::type* = nullptr> +constexpr inline fixed<B, I, F, R> operator/(const fixed<B, I, F, R>& x, T y) noexcept +{ + return fixed<B, I, F, R>(x) /= y; +} + +template <typename B, typename I, unsigned int F, typename T, bool R, typename std::enable_if<std::is_integral<T>::value>::type* = nullptr> +constexpr inline fixed<B, I, F, R> operator/(T x, const fixed<B, I, F, R>& y) noexcept +{ + return fixed<B, I, F, R>(x) /= y; +} + +// +// Comparison operators +// + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline bool operator==(const fixed<B, I, F, R>& x, const fixed<B, I, F, R>& y) noexcept +{ + return x.raw_value() == y.raw_value(); +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline bool operator!=(const fixed<B, I, F, R>& x, const fixed<B, I, F, R>& y) noexcept +{ + return x.raw_value() != y.raw_value(); +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline bool operator<(const fixed<B, I, F, R>& x, const fixed<B, I, F, R>& y) noexcept +{ + return x.raw_value() < y.raw_value(); +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline bool operator>(const fixed<B, I, F, R>& x, const fixed<B, I, F, R>& y) noexcept +{ + return x.raw_value() > y.raw_value(); +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline bool operator<=(const fixed<B, I, F, R>& x, const fixed<B, I, F, R>& y) noexcept +{ + return x.raw_value() <= y.raw_value(); +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline bool operator>=(const fixed<B, I, F, R>& x, const fixed<B, I, F, R>& y) noexcept +{ + return x.raw_value() >= y.raw_value(); +} + +namespace detail +{ +// Number of base-10 digits required to fully represent a number of bits +static constexpr int max_digits10(int bits) +{ + // 8.24 fixed-point equivalent of (int)ceil(bits * std::log10(2)); + using T = long long; + return static_cast<int>((T{bits} * 5050445 + (T{1} << 24) - 1) >> 24); +} + +// Number of base-10 digits that can be fully represented by a number of bits +static constexpr int digits10(int bits) +{ + // 8.24 fixed-point equivalent of (int)(bits * std::log10(2)); + using T = long long; + return static_cast<int>((T{bits} * 5050445) >> 24); +} + +} // namespace detail +} // namespace fpm + +// Specializations for customization points +namespace std +{ + +template <typename B, typename I, unsigned int F, bool R> +struct hash<fpm::fixed<B,I,F,R>> +{ + using argument_type = fpm::fixed<B, I, F, R>; + using result_type = std::size_t; + + result_type operator()(argument_type arg) const noexcept(noexcept(std::declval<std::hash<B>>()(arg.raw_value()))) { + return m_hash(arg.raw_value()); + } + +private: + std::hash<B> m_hash; +}; + +template <typename B, typename I, unsigned int F, bool R> +struct numeric_limits<fpm::fixed<B,I,F,R>> +{ + static constexpr bool is_specialized = true; + static constexpr bool is_signed = std::numeric_limits<B>::is_signed; + static constexpr bool is_integer = false; + static constexpr bool is_exact = true; + static constexpr bool has_infinity = false; + static constexpr bool has_quiet_NaN = false; + static constexpr bool has_signaling_NaN = false; + static constexpr std::float_denorm_style has_denorm = std::denorm_absent; + static constexpr bool has_denorm_loss = false; + static constexpr std::float_round_style round_style = std::round_to_nearest; + static constexpr bool is_iec_559 = false; + static constexpr bool is_bounded = true; + static constexpr bool is_modulo = std::numeric_limits<B>::is_modulo; + static constexpr int digits = std::numeric_limits<B>::digits; + + // Any number with `digits10` significant base-10 digits (that fits in + // the range of the type) is guaranteed to be convertible from text and + // back without change. Worst case, this is 0.000...001, so we can only + // guarantee this case. Nothing more. + static constexpr int digits10 = 1; + + // This is equal to max_digits10 for the integer and fractional part together. + static constexpr int max_digits10 = + fpm::detail::max_digits10(std::numeric_limits<B>::digits - F) + fpm::detail::max_digits10(F); + + static constexpr int radix = 2; + static constexpr int min_exponent = 1 - F; + static constexpr int min_exponent10 = -fpm::detail::digits10(F); + static constexpr int max_exponent = std::numeric_limits<B>::digits - F; + static constexpr int max_exponent10 = fpm::detail::digits10(std::numeric_limits<B>::digits - F); + static constexpr bool traps = true; + static constexpr bool tinyness_before = false; + + static constexpr fpm::fixed<B,I,F,R> lowest() noexcept { + return fpm::fixed<B,I,F,R>::from_raw_value(std::numeric_limits<B>::lowest()); + }; + + static constexpr fpm::fixed<B,I,F,R> min() noexcept { + return lowest(); + } + + static constexpr fpm::fixed<B,I,F,R> max() noexcept { + return fpm::fixed<B,I,F,R>::from_raw_value(std::numeric_limits<B>::max()); + }; + + static constexpr fpm::fixed<B,I,F,R> epsilon() noexcept { + return fpm::fixed<B,I,F,R>::from_raw_value(1); + }; + + static constexpr fpm::fixed<B,I,F,R> round_error() noexcept { + return fpm::fixed<B,I,F,R>(1) / 2; + }; + + static constexpr fpm::fixed<B,I,F,R> denorm_min() noexcept { + return min(); + } +}; + +template <typename B, typename I, unsigned int F, bool R> +constexpr bool numeric_limits<fpm::fixed<B,I,F,R>>::is_specialized; +template <typename B, typename I, unsigned int F, bool R> +constexpr bool numeric_limits<fpm::fixed<B,I,F,R>>::is_signed; +template <typename B, typename I, unsigned int F, bool R> +constexpr bool numeric_limits<fpm::fixed<B,I,F,R>>::is_integer; +template <typename B, typename I, unsigned int F, bool R> +constexpr bool numeric_limits<fpm::fixed<B,I,F,R>>::is_exact; +template <typename B, typename I, unsigned int F, bool R> +constexpr bool numeric_limits<fpm::fixed<B,I,F,R>>::has_infinity; +template <typename B, typename I, unsigned int F, bool R> +constexpr bool numeric_limits<fpm::fixed<B,I,F,R>>::has_quiet_NaN; +template <typename B, typename I, unsigned int F, bool R> +constexpr bool numeric_limits<fpm::fixed<B,I,F,R>>::has_signaling_NaN; +template <typename B, typename I, unsigned int F, bool R> +constexpr std::float_denorm_style numeric_limits<fpm::fixed<B,I,F,R>>::has_denorm; +template <typename B, typename I, unsigned int F, bool R> +constexpr bool numeric_limits<fpm::fixed<B,I,F,R>>::has_denorm_loss; +template <typename B, typename I, unsigned int F, bool R> +constexpr std::float_round_style numeric_limits<fpm::fixed<B,I,F,R>>::round_style; +template <typename B, typename I, unsigned int F, bool R> +constexpr bool numeric_limits<fpm::fixed<B,I,F,R>>::is_iec_559; +template <typename B, typename I, unsigned int F, bool R> +constexpr bool numeric_limits<fpm::fixed<B,I,F,R>>::is_bounded; +template <typename B, typename I, unsigned int F, bool R> +constexpr bool numeric_limits<fpm::fixed<B,I,F,R>>::is_modulo; +template <typename B, typename I, unsigned int F, bool R> +constexpr int numeric_limits<fpm::fixed<B,I,F,R>>::digits; +template <typename B, typename I, unsigned int F, bool R> +constexpr int numeric_limits<fpm::fixed<B,I,F,R>>::digits10; +template <typename B, typename I, unsigned int F, bool R> +constexpr int numeric_limits<fpm::fixed<B,I,F,R>>::max_digits10; +template <typename B, typename I, unsigned int F, bool R> +constexpr int numeric_limits<fpm::fixed<B,I,F,R>>::radix; +template <typename B, typename I, unsigned int F, bool R> +constexpr int numeric_limits<fpm::fixed<B,I,F,R>>::min_exponent; +template <typename B, typename I, unsigned int F, bool R> +constexpr int numeric_limits<fpm::fixed<B,I,F,R>>::min_exponent10; +template <typename B, typename I, unsigned int F, bool R> +constexpr int numeric_limits<fpm::fixed<B,I,F,R>>::max_exponent; +template <typename B, typename I, unsigned int F, bool R> +constexpr int numeric_limits<fpm::fixed<B,I,F,R>>::max_exponent10; +template <typename B, typename I, unsigned int F, bool R> +constexpr bool numeric_limits<fpm::fixed<B,I,F,R>>::traps; +template <typename B, typename I, unsigned int F, bool R> +constexpr bool numeric_limits<fpm::fixed<B,I,F,R>>::tinyness_before; + +} + +#endif diff --git a/fpm/ios.hpp b/fpm/ios.hpp new file mode 100644 index 0000000..69581fb --- /dev/null +++ b/fpm/ios.hpp @@ -0,0 +1,740 @@ +#ifndef FPM_IOS_HPP +#define FPM_IOS_HPP + +#include "fixed.hpp" +#include "math.hpp" +#include <array> +#include <algorithm> +#include <cctype> +#include <climits> +#include <limits> +#include <ios> +#include <vector> + +namespace fpm +{ + +template <typename CharT, typename B, typename I, unsigned int F, bool R> +std::basic_ostream<CharT>& operator<<(std::basic_ostream<CharT>& os, fixed<B, I, F, R> x) noexcept +{ + const auto uppercase = ((os.flags() & std::ios_base::uppercase) != 0); + const auto showpoint = ((os.flags() & std::ios_base::showpoint) != 0); + const auto adjustfield = (os.flags() & std::ios_base::adjustfield); + const auto width = os.width(); + const auto& ctype = std::use_facet<std::ctype<CharT>>(os.getloc()); + const auto& numpunct = std::use_facet<std::numpunct<CharT>>(os.getloc()); + + auto floatfield = (os.flags() & std::ios_base::floatfield); + auto precision = os.precision(); + auto show_trailing_zeros = true; + auto use_significant_digits = false; + + // Invalid precision? Reset to the default + if (precision < 0) + { + precision = 6; + } + + // Output buffer. Needs to be big enough for the formatted number without padding. + // Optional prefixes (i.e. "+"/"-", decimal separator, exponent "e+/-" and/or "0x"). + constexpr auto worst_case_constant_size = 6; + // Maximum number of digits from the base type (covers integral + fractional digits) + constexpr auto worst_case_digit_count = std::numeric_limits<B>::digits10 + 2; + // Exponent suffixes (i.e. maximum digits based on log of the base type size). + // Needs a log10, but that isn't constexpr, so we're over-allocating on the stack. Can't hurt. + constexpr auto worst_case_suffix_size = std::numeric_limits<B>::digits; + // Double the digit count: in the worst case the thousands grouping add a character per digit. + using buffer_t = std::array<CharT, worst_case_constant_size + worst_case_digit_count * 2 + worst_case_suffix_size>; + buffer_t buffer; + + // Output cursor + auto end = buffer.begin(); + + // Keep track of the start of "internal" padding + typename buffer_t::iterator internal_pad = buffer.end(); + + // Representation of a number. + // The value of the number is: raw / divisor * (10|2) ^ exponent + // The base of the exponent is 2 in hexfloat mode, or 10 otherwise. + struct number_t { + I raw; // raw fixed-point value + I divisor; // the divisor indicating the place of the decimal point + int exponent; // the exponent applied + }; + + // Convert a value without exponent to scientific representation + // where the part before the decimal point is less than 10. + const auto as_scientific = [](number_t value) { + assert(value.exponent == 0); + if (value.raw > 0) + { + while (value.raw / 10 >= value.divisor) { + value.divisor *= 10; + ++value.exponent; + } + while (value.raw < value.divisor) { + value.raw *= 10; + --value.exponent; + } + } + return value; + }; + + number_t value = { x.raw_value(), I{1} << F, 0}; + + auto base = B{10}; + + // First write the sign + if (value.raw < 0) + { + *end++ = ctype.widen('-'); + value.raw = -value.raw; + internal_pad = end; + } + else if (os.flags() & std::ios_base::showpos) + { + *end++ = ctype.widen('+'); + internal_pad = end; + } + assert(value.raw >= 0); + + switch (floatfield) + { + case std::ios_base::fixed | std::ios_base::scientific: + // Hexadecimal mode: figure out the hexadecimal exponent and write "0x" + if (value.raw > 0) + { + auto bit = detail::find_highest_bit(value.raw); + value.exponent = bit - F; // exponent is applied to base 2 + value.divisor = I{1} << bit; // divisor is at the highest bit, ensuring it starts with "1." + precision = (bit + 3) / 4; // precision is number of nibbles, so we show all of them + } + base = 16; + show_trailing_zeros = false; // Always strip trailing zeros in hexfloat mode + + *end++ = ctype.widen('0'); + *end++ = ctype.widen(uppercase ? 'X' : 'x'); + break; + + case std::ios_base::scientific: + // Scientific mode, normalize value to scientific notation + value = as_scientific(value); + break; + + case std::ios_base::fixed: + // Fixed mode. Nothing to do. + break; + + default: + { + // "auto" mode: figure out the exponent + const number_t sci_value = as_scientific(value); + + // Now `precision` indicates the number of *significant digits* (not fractional digits). + use_significant_digits = true; + precision = std::max<std::streamsize>(precision, 1); + + if (sci_value.exponent >= precision || sci_value.exponent < -4) { + // Display as scientific format + floatfield = std::ios_base::scientific; + value = sci_value; + } else { + // Display as fixed format. + // "showpoint" indicates whether or not we show trailing zeros + floatfield = std::ios_base::fixed; + show_trailing_zeros = showpoint; + } + break; + } + }; + + // If we didn't write a sign, any internal padding starts here + // (after a potential "0x" for hexfloats). + if (internal_pad == buffer.end()) { + internal_pad = end; + } + + // Separate out the integral part of the number + I integral = value.raw / value.divisor; + value.raw %= value.divisor; + + // Here we start printing the number itself + const char* const digits = uppercase ? "0123456789ABCDEF" : "0123456789abcdef"; + const auto digits_start = end; + + // Are we already printing significant digits? (yes if we're not counting significant digits) + bool significant_digits = !use_significant_digits; + + // Print the integral part + int last_digit = 0; + if (integral == 0) { + *end++ = ctype.widen('0'); + if (value.raw == 0) { + // If the fraction is zero too, all zeros including the integral count + // as significant digits. + significant_digits = true; + } + } else { + while (integral > 0) { + last_digit = integral % base; + *end++ = ctype.widen(digits[last_digit]); + integral /= base; + } + std::reverse(digits_start, end); + significant_digits = true; + } + + if (use_significant_digits && significant_digits) + { + // Apparently the integral part was significant; subtract its + // length from the remaining significant digits. + precision -= (end - digits_start); + } + + // At this point, `value` contains only the fraction and + // `precision` holds the number of digits to print. + assert(value.raw < value.divisor); + assert(precision >= 0); + + // Location of decimal point + typename buffer_t::iterator point = buffer.end(); + + // Start (and length) of the trailing zeros to insert while printing + // By tracking this to print them later instead of actually printing them now, + // we can support large precisions with a small printing buffer. + typename buffer_t::iterator trailing_zeros_start = buffer.end(); + std::streamsize trailing_zeros_count = 0; + + if (precision > 0) + { + // Print the fractional part + *(point = end++) = numpunct.decimal_point(); + + for (int i = 0; i < precision; ++i) + { + if (value.raw == 0) + { + // The rest of the digits are all zeros, mark them + // to be printed in this spot. + trailing_zeros_start = end; + trailing_zeros_count = precision - i; + break; + } + + // Shift the divisor if we can to avoid overflow on the value + if (value.divisor % base == 0) { + value.divisor /= base; + } else { + value.raw *= base; + } + assert(value.divisor > 0); + assert(value.raw >= 0); + last_digit = (value.raw / value.divisor) % base; + value.raw %= value.divisor; + *end++ = ctype.widen(digits[last_digit]); + + if (!significant_digits) { + // We're still finding the first significant digit + if (last_digit != 0) { + // Found it + significant_digits = true; + } else { + // Not yet; increment number of digits to print + ++precision; + } + } + } + } + else if (showpoint) + { + // No fractional part to print, but we still want the point + *(point = end++) = numpunct.decimal_point(); + } + + // Insert `ch` into the output at `position`, updating all references accordingly + const auto insert_character = [&](typename buffer_t::iterator position, CharT ch) { + assert(position >= buffer.begin() && position < end); + std::move_backward(position, end, end + 1); + if (point != buffer.end() && position < point) { + ++point; + } + if (trailing_zeros_start != buffer.end() && position < trailing_zeros_start) { + ++trailing_zeros_start; + } + ++end; + *position = ch; + }; + + // Round the number: round to nearest + bool increment = false; + if (value.raw > value.divisor / 2) { + // Round up + increment = true; + } else if (value.raw == value.divisor / 2) { + // It's a tie (i.e. "xyzw.5"): round to even + increment = ((last_digit % 2) == 1); + } + + if (increment) + { + auto p = end - 1; + // Increment all digits backwards while we see "9" + while (p >= digits_start) { + if (p == point) { + // Skip over the decimal point + --p; + } + if ((*p)++ != ctype.widen('9')) { + break; + } + *p-- = ctype.widen('0'); + } + + if (p < digits_start) { + // We've incremented all the way to the start (all 9's), we need to insert the + // carried-over 1 from incrementing the last 9. + assert(p == digits_start - 1); + insert_character(++p, ctype.widen('1')); + + if (floatfield == std::ios::scientific) + { + // We just made the integral part equal to 10, so we shift the decimal point + // back one place (if any) and tweak the exponent, so that we keep the integer part + // less than 10. + if (point != buffer.end()) { + assert(p + 2 == point); + std::swap(*(point - 1), *point); + --point; + } + ++value.exponent; + + // We've introduced an extra digit so we need to strip the last digit + // to maintain the same precision + --end; + } + } + + if (use_significant_digits && *p == ctype.widen('1') && point != buffer.end()) { + // We've converted a leading zero to a 1 so we need to strip the last digit + // (behind the decimal point) to maintain the same significant digit count. + --end; + } + } + + if (point != buffer.end()) + { + if (!show_trailing_zeros) + { + // Remove trailing zeros + while (*(end - 1) == ctype.widen('0')) { + --end; + } + + // Also clear the "trailing zeros to append during printing" range + trailing_zeros_start = buffer.end(); + trailing_zeros_count = 0; + } + + if (end - 1 == point && trailing_zeros_count == 0 && !showpoint) { + // Remove the decimal point, too + --end; + } + } + + // Apply thousands grouping + const auto& grouping = numpunct.grouping(); + if (!grouping.empty()) + { + // Step backwards from the end or decimal point, inserting the + // thousands separator at every group interval. + const CharT thousands_sep = ctype.widen(numpunct.thousands_sep()); + std::size_t group = 0; + auto p = point != buffer.end() ? point : end; + auto size = static_cast<int>(grouping[group]); + while (size > 0 && size < CHAR_MAX && p - digits_start > size) { + p -= size; + insert_character(p, thousands_sep); + if (group < grouping.size() - 1) { + size = static_cast<int>(grouping[++group]); + } + } + } + + // Print the exponent if required + assert(floatfield != 0); + if (floatfield & std::ios_base::scientific) + { + // Hexadecimal (%a/%A) or decimal (%e/%E) scientific notation + if (floatfield & std::ios_base::fixed) { + *end++ = ctype.widen(uppercase ? 'P' : 'p'); + } else { + *end++ = ctype.widen(uppercase ? 'E' : 'e'); + } + + if (value.exponent < 0) { + *end++ = ctype.widen('-'); + value.exponent = -value.exponent; + } else { + *end++ = ctype.widen('+'); + } + + if (floatfield == std::ios_base::scientific) { + // In decimal scientific notation (%e/%E), the exponent is at least two digits + if (value.exponent < 10) { + *end++ = ctype.widen('0'); + } + } + + const auto exponent_start = end; + if (value.exponent == 0) { + *end++ = ctype.widen('0'); + } else while (value.exponent > 0) { + *end++ = ctype.widen(digits[value.exponent % 10]); + value.exponent /= 10; + } + std::reverse(exponent_start, end); + } + + // Write character `ch` `count` times to the stream + const auto sputcn = [&](CharT ch, std::streamsize count){ + // Fill a buffer to output larger chunks + constexpr std::streamsize chunk_size = 64; + std::array<CharT, chunk_size> fill_buffer; + std::fill_n(fill_buffer.begin(), std::min(count, chunk_size), ch); + + for (std::streamsize size, left = count; left > 0; left -= size) { + size = std::min(chunk_size, left); + os.rdbuf()->sputn(&fill_buffer[0], size); + } + }; + + // Outputs a range of characters, making sure to output the trailing zeros range + // if it lies in the specified range + const auto put_range = [&](typename buffer_t::const_iterator begin, typename buffer_t::const_iterator end) { + assert(end >= begin); + if (trailing_zeros_start >= begin && trailing_zeros_start <= end) { + // Print range with trailing zeros range in the middle + assert(trailing_zeros_count > 0); + os.rdbuf()->sputn(&*begin, trailing_zeros_start - begin); + sputcn(ctype.widen('0'), trailing_zeros_count); + os.rdbuf()->sputn(&*trailing_zeros_start, end - trailing_zeros_start); + } else { + // Print range as-is + os.rdbuf()->sputn(&*begin, end - begin); + } + }; + + // Pad the buffer if necessary. + // Note that the length of trailing zeros is counted towards the length of the content. + const auto content_size = end - buffer.begin() + trailing_zeros_count; + if (content_size >= width) + { + // Buffer needs no padding, output as-is + put_range(buffer.begin(), end); + } + else + { + const auto pad_size = width - content_size; + switch (adjustfield) + { + case std::ios_base::left: + // Content is left-aligned, so output the buffer, followed by the padding + put_range(buffer.begin(), end); + sputcn(os.fill(), pad_size); + break; + case std::ios_base::internal: + // Content is internally aligned, so output the buffer up to the "internal pad" + // point, followed by the padding, followed by the remainder of the buffer. + put_range(buffer.begin(), internal_pad); + sputcn(os.fill(), pad_size); + put_range(internal_pad, end); + break; + default: + // Content is right-aligned, so output the padding, followed by the buffer + sputcn(os.fill(), pad_size); + put_range(buffer.begin(), end); + break; + } + } + + // Width is reset after every write + os.width(0); + + return os; +} + + +template <typename CharT, class Traits, typename B, typename I, unsigned int F, bool R> +std::basic_istream<CharT, Traits>& operator>>(std::basic_istream<CharT, Traits>& is, fixed<B, I, F, R>& x) +{ + typename std::basic_istream<CharT, Traits>::sentry sentry(is); + if (!sentry) + { + return is; + } + + const auto& ctype = std::use_facet<std::ctype<CharT>>(is.getloc()); + const auto& numpunct = std::use_facet<std::numpunct<CharT>>(is.getloc()); + + bool thousands_separator_allowed = false; + const bool supports_thousands_separators = !numpunct.grouping().empty(); + + const auto& is_valid_character = [](char ch) { + // Note: allowing ['p', 'i', 'n', 't', 'y'] is technically in violation of the spec (we are emulating std::num_get), + // but otherwise we cannot parse hexfloats and "infinity". This is a known issue with the spec (LWG #2381). + return std::isxdigit(ch) || + ch == 'x' || ch == 'X' || ch == 'p' || ch == 'P' || + ch == 'i' || ch == 'I' || ch == 'n' || ch == 'N' || + ch == 't' || ch == 'T' || ch == 'y' || ch == 'Y' || + ch == '-' || ch == '+'; + }; + + const auto& peek = [&]() { + for(;;) { + auto ch = is.rdbuf()->sgetc(); + if (ch == Traits::eof()) { + is.setstate(std::ios::eofbit); + return '\0'; + } + if (ch == numpunct.decimal_point()) { + return '.'; + } + if (ch == numpunct.thousands_sep()) + { + if (!supports_thousands_separators || !thousands_separator_allowed) { + return '\0'; + } + // Ignore valid thousands separators + is.rdbuf()->sbumpc(); + continue; + } + auto res = ctype.narrow(ch, 0); + if (!is_valid_character(res)) { + // Invalid character: end input + return '\0'; + } + return res; + } + }; + + const auto& bump = [&]() { + is.rdbuf()->sbumpc(); + }; + + const auto& next = [&]() { + bump(); + return peek(); + }; + + bool negate = false; + auto ch = peek(); + if (ch == '-') { + negate = true; + ch = next(); + } else if (ch == '+') { + ch = next(); + } + + const char infinity[] = "infinity"; + // Must be "inf" or "infinity" + int i = 0; + while (i < 8 && ch == infinity[i]) { + ++i; + ch = next(); + } + + if (i > 0) { + if (i == 3 || i == 8) { + x = negate ? std::numeric_limits<fixed<B, I, F, R>>::min() : std::numeric_limits<fixed<B, I, F, R>>::max(); + } else { + is.setstate(std::ios::failbit); + } + return is; + } + + char exponent_char = 'e'; + int base = 10; + + constexpr auto NoFraction = std::numeric_limits<std::size_t>::max(); + std::size_t fraction_start = NoFraction; + std::vector<unsigned char> significand; + + if (ch == '0') { + ch = next(); + if (ch == 'x' || ch == 'X') { + // Hexfloat + exponent_char = 'p'; + base = 16; + ch = next(); + } else { + significand.push_back(0); + } + } + + // Parse the significand + thousands_separator_allowed = true; + for (;; ch = next()) { + if (ch == '.') { + if (fraction_start != NoFraction) { + // Double decimal point. Stop parsing. + break; + } + fraction_start = significand.size(); + thousands_separator_allowed = false; + } else { + unsigned char val = base; + if (ch >= '0' && ch <= '9') { + val = ch - '0'; + } else if (ch >= 'a' && ch <= 'f') { + val = ch - 'a' + 10; + } else if (ch >= 'A' && ch <= 'F') { + val = ch - 'A' + 10; + } + if (val < 0 || val >= base) { + break; + } + significand.push_back(val); + } + } + if (significand.empty()) { + // We need a significand + is.setstate(std::ios::failbit); + return is; + } + thousands_separator_allowed = false; + + if (fraction_start == NoFraction) { + // If we haven't seen a fraction yet, place it at the end of the significand + fraction_start = significand.size(); + } + + // Parse the exponent + bool exponent_overflow = false; + std::size_t exponent = 0; + bool exponent_negate = false; + if (std::tolower(ch) == exponent_char) + { + ch = next(); + if (ch == '-') { + exponent_negate = true; + ch = next(); + } else if (ch == '+') { + ch = next(); + } + + bool parsed = false; + while (std::isdigit(ch)) { + if (exponent <= std::numeric_limits<int>::max() / 10) { + exponent = exponent * 10 + (ch - '0'); + } else { + exponent_overflow = true; + } + parsed = true; + ch = next(); + } + if (!parsed) { + // If the exponent character is given, the exponent value may not be empty + is.setstate(std::ios::failbit); + return is; + } + } + + // We've parsed all we need. Construct the value. + if (exponent_overflow) { + // Absolute exponent is too large + if (std::all_of(significand.begin(), significand.end(), [](unsigned char x){ return x == 0; })) { + // Significand is zero. Exponent doesn't matter. + x = fixed<B, I, F, R>(0); + } else if (exponent_negate) { + // A huge negative exponent approaches 0. + x = fixed<B, I, F, R>::from_raw_value(0); + } else { + // A huge positive exponent approaches infinity. + x = std::numeric_limits<fixed<B, I, F, R>>::max(); + } + return is; + } + + // Shift the fraction offset according to exponent + { + const auto exponent_mult = (base == 10) ? 1: 4; + if (exponent_negate) { + const auto adjust = std::min(exponent / exponent_mult, fraction_start); + fraction_start -= adjust; + exponent -= adjust * exponent_mult; + } else { + const auto adjust = std::min(exponent / exponent_mult, significand.size() - fraction_start); + fraction_start += adjust; + exponent -= adjust * exponent_mult; + } + } + + constexpr auto IsSigned = std::is_signed<B>::value; + constexpr auto IntBits = sizeof(B) * 8 - F - (IsSigned ? 1 : 0); + constexpr auto MaxInt = (I{1} << IntBits) - 1; + constexpr auto MaxFraction = (I{1} << F) - 1; + constexpr auto MaxValue = (I{1} << sizeof(B) * 8) - 1; + + // Parse the integer part + I integer = 0; + for (std::size_t i = 0; i < fraction_start; ++i) { + if (integer > MaxInt / base) { + // Overflow + x = negate ? std::numeric_limits<fixed<B, I, F, R>>::min() : std::numeric_limits<fixed<B, I, F, R>>::max(); + return is; + } + assert(significand[i] < base); + integer = integer * base + significand[i]; + } + + // Parse the fractional part + I fraction = 0; + I divisor = 1; + for (std::size_t i = fraction_start; i < significand.size(); ++i) { + assert(significand[i] < base); + if (divisor > MaxFraction / base) { + // We're done + break; + } + fraction = fraction * base + significand[i]; + divisor *= base; + } + + // Construct the value from the parsed parts + I raw_value = (integer << F) + (fraction << F) / divisor; + + // Apply remaining exponent + if (exponent_char == 'p') { + // Base-2 exponent + if (exponent_negate) { + raw_value >>= exponent; + } else { + raw_value <<= exponent; + } + } else { + // Base-10 exponent + if (exponent_negate) { + I remainder = 0; + for (std::size_t e = 0; e < exponent; ++e) { + remainder = raw_value % 10; + raw_value /= 10; + } + raw_value += remainder / 5; + } else { + for (std::size_t e = 0; e < exponent; ++e) { + if (raw_value > MaxValue / 10) { + // Overflow + x = negate ? std::numeric_limits<fixed<B, I, F, R>>::min() : std::numeric_limits<fixed<B, I, F, R>>::max(); + return is; + } + raw_value *= 10; + } + } + } + x = fixed<B, I, F, R>::from_raw_value(static_cast<B>(negate ? -raw_value : raw_value)); + return is; +} + +} + +#endif diff --git a/fpm/math.hpp b/fpm/math.hpp new file mode 100644 index 0000000..7a76349 --- /dev/null +++ b/fpm/math.hpp @@ -0,0 +1,684 @@ +#ifndef FPM_MATH_HPP +#define FPM_MATH_HPP + +#include "fixed.hpp" +#include <cmath> + +#ifdef _MSC_VER +#include <intrin.h> +#endif + +namespace fpm +{ + +// +// Helper functions +// +namespace detail +{ + +// Returns the index of the most-signifcant set bit +inline long find_highest_bit(unsigned long long value) noexcept +{ + assert(value != 0); +#if defined(_MSC_VER) + unsigned long index; +#if defined(_WIN64) + _BitScanReverse64(&index, value); +#else + if (_BitScanReverse(&index, static_cast<unsigned long>(value >> 32)) != 0) { + index += 32; + } else { + _BitScanReverse(&index, static_cast<unsigned long>(value & 0xfffffffflu)); + } +#endif + return index; +#elif defined(__GNUC__) || defined(__clang__) + return sizeof(value) * 8 - 1 - __builtin_clzll(value); +#else +# error "your platform does not support find_highest_bit()" +#endif +} + +} + +// +// Classification methods +// + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline int fpclassify(fixed<B, I, F, R> x) noexcept +{ + return (x.raw_value() == 0) ? FP_ZERO : FP_NORMAL; +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline bool isfinite(fixed<B, I, F, R>) noexcept +{ + return true; +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline bool isinf(fixed<B, I, F, R>) noexcept +{ + return false; +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline bool isnan(fixed<B, I, F, R>) noexcept +{ + return false; +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline bool isnormal(fixed<B, I, F, R> x) noexcept +{ + return x.raw_value() != 0; +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline bool signbit(fixed<B, I, F, R> x) noexcept +{ + return x.raw_value() < 0; +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline bool isgreater(fixed<B, I, F, R> x, fixed<B, I, F, R> y) noexcept +{ + return x > y; +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline bool isgreaterequal(fixed<B, I, F, R> x, fixed<B, I, F, R> y) noexcept +{ + return x >= y; +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline bool isless(fixed<B, I, F, R> x, fixed<B, I, F, R> y) noexcept +{ + return x < y; +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline bool islessequal(fixed<B, I, F, R> x, fixed<B, I, F, R> y) noexcept +{ + return x <= y; +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline bool islessgreater(fixed<B, I, F, R> x, fixed<B, I, F, R> y) noexcept +{ + return x != y; +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline bool isunordered(fixed<B, I, F, R> x, fixed<B, I, F, R> y) noexcept +{ + return false; +} + +// +// Nearest integer operations +// +template <typename B, typename I, unsigned int F, bool R> +inline fixed<B, I, F, R> ceil(fixed<B, I, F, R> x) noexcept +{ + constexpr auto FRAC = B(1) << F; + auto value = x.raw_value(); + if (value > 0) value += FRAC - 1; + return fixed<B, I, F, R>::from_raw_value(value / FRAC * FRAC); +} + +template <typename B, typename I, unsigned int F, bool R> +inline fixed<B, I, F, R> floor(fixed<B, I, F, R> x) noexcept +{ + constexpr auto FRAC = B(1) << F; + auto value = x.raw_value(); + if (value < 0) value -= FRAC - 1; + return fixed<B, I, F, R>::from_raw_value(value / FRAC * FRAC); +} + +template <typename B, typename I, unsigned int F, bool R> +inline fixed<B, I, F, R> trunc(fixed<B, I, F, R> x) noexcept +{ + constexpr auto FRAC = B(1) << F; + return fixed<B, I, F, R>::from_raw_value(x.raw_value() / FRAC * FRAC); +} + +template <typename B, typename I, unsigned int F, bool R> +inline fixed<B, I, F, R> round(fixed<B, I, F, R> x) noexcept +{ + constexpr auto FRAC = B(1) << F; + auto value = x.raw_value() / (FRAC / 2); + return fixed<B, I, F, R>::from_raw_value(((value / 2) + (value % 2)) * FRAC); +} + +template <typename B, typename I, unsigned int F, bool R> +fixed<B, I, F, R> nearbyint(fixed<B, I, F, R> x) noexcept +{ + // Rounding mode is assumed to be FE_TONEAREST + constexpr auto FRAC = B(1) << F; + auto value = x.raw_value(); + const bool is_half = std::abs(value % FRAC) == FRAC / 2; + value /= FRAC / 2; + value = (value / 2) + (value % 2); + value -= (value % 2) * is_half; + return fixed<B, I, F, R>::from_raw_value(value * FRAC); +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline fixed<B, I, F, R> rint(fixed<B, I, F, R> x) noexcept +{ + // Rounding mode is assumed to be FE_TONEAREST + return nearbyint(x); +} + +// +// Mathematical functions +// +template <typename B, typename I, unsigned int F, bool R> +constexpr inline fixed<B, I, F, R> abs(fixed<B, I, F, R> x) noexcept +{ + return (x >= fixed<B, I, F, R>{0}) ? x : -x; +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline fixed<B, I, F, R> fmod(fixed<B, I, F, R> x, fixed<B, I, F, R> y) noexcept +{ + return + assert(y.raw_value() != 0), + fixed<B, I, F, R>::from_raw_value(x.raw_value() % y.raw_value()); +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline fixed<B, I, F, R> remainder(fixed<B, I, F, R> x, fixed<B, I, F, R> y) noexcept +{ + return + assert(y.raw_value() != 0), + x - nearbyint(x / y) * y; +} + +template <typename B, typename I, unsigned int F, bool R> +inline fixed<B, I, F, R> remquo(fixed<B, I, F, R> x, fixed<B, I, F, R> y, int* quo) noexcept +{ + assert(y.raw_value() != 0); + assert(quo != nullptr); + *quo = x.raw_value() / y.raw_value(); + return fixed<B, I, F, R>::from_raw_value(x.raw_value() % y.raw_value()); +} + +// +// Manipulation functions +// + +template <typename B, typename I, unsigned int F, bool R, typename C, typename J, unsigned int G, bool S> +constexpr inline fixed<B, I, F, R> copysign(fixed<B, I, F, R> x, fixed<C, J, G, S> y) noexcept +{ + return + x = abs(x), + (y >= fixed<C, J, G, S>{0}) ? x : -x; +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline fixed<B, I, F, R> nextafter(fixed<B, I, F, R> from, fixed<B, I, F, R> to) noexcept +{ + return from == to ? to : + to > from ? fixed<B, I, F, R>::from_raw_value(from.raw_value() + 1) + : fixed<B, I, F, R>::from_raw_value(from.raw_value() - 1); +} + +template <typename B, typename I, unsigned int F, bool R> +constexpr inline fixed<B, I, F, R> nexttoward(fixed<B, I, F, R> from, fixed<B, I, F, R> to) noexcept +{ + return nextafter(from, to); +} + +template <typename B, typename I, unsigned int F, bool R> +inline fixed<B, I, F, R> modf(fixed<B, I, F, R> x, fixed<B, I, F, R>* iptr) noexcept +{ + const auto raw = x.raw_value(); + constexpr auto FRAC = B{1} << F; + *iptr = fixed<B, I, F, R>::from_raw_value(raw / FRAC * FRAC); + return fixed<B, I, F, R>::from_raw_value(raw % FRAC); +} + + +// +// Power functions +// + +template <typename B, typename I, unsigned int F, bool R, typename T, typename std::enable_if<std::is_integral<T>::value>::type* = nullptr> +fixed<B, I, F, R> pow(fixed<B, I, F, R> base, T exp) noexcept +{ + using Fixed = fixed<B, I, F, R>; + + if (base == Fixed(0)) { + assert(exp > 0); + return Fixed(0); + } + + Fixed result {1}; + if (exp < 0) + { + for (Fixed intermediate = base; exp != 0; exp /= 2, intermediate *= intermediate) + { + if ((exp % 2) != 0) + { + result /= intermediate; + } + } + } + else + { + for (Fixed intermediate = base; exp != 0; exp /= 2, intermediate *= intermediate) + { + if ((exp % 2) != 0) + { + result *= intermediate; + } + } + } + return result; +} + +template <typename B, typename I, unsigned int F, bool R> +fixed<B, I, F, R> pow(fixed<B, I, F, R> base, fixed<B, I, F, R> exp) noexcept +{ + using Fixed = fixed<B, I, F, R>; + + if (base == Fixed(0)) { + assert(exp > Fixed(0)); + return Fixed(0); + } + + if (exp < Fixed(0)) + { + return 1 / pow(base, -exp); + } + + constexpr auto FRAC = B(1) << F; + if (exp.raw_value() % FRAC == 0) + { + // Non-fractional exponents are easier to calculate + return pow(base, exp.raw_value() / FRAC); + } + + // For negative bases we do not support fractional exponents. + // Technically fractions with odd denominators could work, + // but that's too much work to figure out. + assert(base > Fixed(0)); + return exp2(log2(base) * exp); +} + +template <typename B, typename I, unsigned int F, bool R> +fixed<B, I, F, R> exp(fixed<B, I, F, R> x) noexcept +{ + using Fixed = fixed<B, I, F, R>; + if (x < Fixed(0)) { + return 1 / exp(-x); + } + constexpr auto FRAC = B(1) << F; + const B x_int = x.raw_value() / FRAC; + x -= x_int; + assert(x >= Fixed(0) && x < Fixed(1)); + + constexpr auto fA = Fixed::template from_fixed_point<63>( 128239257017632854ll); // 1.3903728105644451e-2 + constexpr auto fB = Fixed::template from_fixed_point<63>( 320978614890280666ll); // 3.4800571158543038e-2 + constexpr auto fC = Fixed::template from_fixed_point<63>(1571680799599592947ll); // 1.7040197373796334e-1 + constexpr auto fD = Fixed::template from_fixed_point<63>(4603349000587966862ll); // 4.9909609871464493e-1 + constexpr auto fE = Fixed::template from_fixed_point<62>(4612052447974689712ll); // 1.0000794567422495 + constexpr auto fF = Fixed::template from_fixed_point<63>(9223361618412247875ll); // 9.9999887043019773e-1 + return pow(Fixed::e(), x_int) * (((((fA * x + fB) * x + fC) * x + fD) * x + fE) * x + fF); +} + +template <typename B, typename I, unsigned int F, bool R> +fixed<B, I, F, R> exp2(fixed<B, I, F, R> x) noexcept +{ + using Fixed = fixed<B, I, F, R>; + if (x < Fixed(0)) { + return 1 / exp2(-x); + } + constexpr auto FRAC = B(1) << F; + const B x_int = x.raw_value() / FRAC; + x -= x_int; + assert(x >= Fixed(0) && x < Fixed(1)); + + constexpr auto fA = Fixed::template from_fixed_point<63>( 17491766697771214ll); // 1.8964611454333148e-3 + constexpr auto fB = Fixed::template from_fixed_point<63>( 82483038782406547ll); // 8.9428289841091295e-3 + constexpr auto fC = Fixed::template from_fixed_point<63>( 515275173969157690ll); // 5.5866246304520701e-2 + constexpr auto fD = Fixed::template from_fixed_point<63>(2214897896212987987ll); // 2.4013971109076949e-1 + constexpr auto fE = Fixed::template from_fixed_point<63>(6393224161192452326ll); // 6.9315475247516736e-1 + constexpr auto fF = Fixed::template from_fixed_point<63>(9223371050976163566ll); // 9.9999989311082668e-1 + return Fixed(1 << x_int) * (((((fA * x + fB) * x + fC) * x + fD) * x + fE) * x + fF); +} + +template <typename B, typename I, unsigned int F, bool R> +fixed<B, I, F, R> expm1(fixed<B, I, F, R> x) noexcept +{ + return exp(x) - 1; +} + +template <typename B, typename I, unsigned int F, bool R> +fixed<B, I, F, R> log2(fixed<B, I, F, R> x) noexcept +{ + using Fixed = fixed<B, I, F, R>; + assert(x > Fixed(0)); + + // Normalize input to the [1:2] domain + B value = x.raw_value(); + const long highest = detail::find_highest_bit(value); + if (highest >= F) { + value >>= (highest - F); + } else { + value <<= (F - highest); + } + x = Fixed::from_raw_value(value); + assert(x >= Fixed(1) && x < Fixed(2)); + + constexpr auto fA = Fixed::template from_fixed_point<63>( 413886001457275979ll); // 4.4873610194131727e-2 + constexpr auto fB = Fixed::template from_fixed_point<63>(-3842121857793256941ll); // -4.1656368651734915e-1 + constexpr auto fC = Fixed::template from_fixed_point<62>( 7522345947206307744ll); // 1.6311487636297217 + constexpr auto fD = Fixed::template from_fixed_point<61>(-8187571043052183818ll); // -3.5507929249026341 + constexpr auto fE = Fixed::template from_fixed_point<60>( 5870342889289496598ll); // 5.0917108110420042 + constexpr auto fF = Fixed::template from_fixed_point<61>(-6457199832668582866ll); // -2.8003640347009253 + return Fixed(highest - F) + (((((fA * x + fB) * x + fC) * x + fD) * x + fE) * x + fF); +} + +template <typename B, typename I, unsigned int F, bool R> +fixed<B, I, F, R> log(fixed<B, I, F, R> x) noexcept +{ + using Fixed = fixed<B, I, F, R>; + return log2(x) / log2(Fixed::e()); +} + +template <typename B, typename I, unsigned int F, bool R> +fixed<B, I, F, R> log10(fixed<B, I, F, R> x) noexcept +{ + using Fixed = fixed<B, I, F, R>; + return log2(x) / log2(Fixed(10)); +} + +template <typename B, typename I, unsigned int F, bool R> +fixed<B, I, F, R> log1p(fixed<B, I, F, R> x) noexcept +{ + return log(1 + x); +} + +template <typename B, typename I, unsigned int F, bool R> +fixed<B, I, F, R> cbrt(fixed<B, I, F, R> x) noexcept +{ + using Fixed = fixed<B, I, F, R>; + + if (x == Fixed(0)) + { + return x; + } + if (x < Fixed(0)) + { + return -cbrt(-x); + } + assert(x >= Fixed(0)); + + // Finding the cube root of an integer, taken from Hacker's Delight, + // based on the square root algorithm. + + // We start at the greatest power of eight that's less than the argument. + int ofs = ((detail::find_highest_bit(x.raw_value()) + 2*F) / 3 * 3); + I num = I{x.raw_value()}; + I res = 0; + + const auto do_round = [&] + { + for (; ofs >= 0; ofs -= 3) + { + res += res; + const I val = (3*res*(res + 1) + 1) << ofs; + if (num >= val) + { + num -= val; + res++; + } + } + }; + + // We should shift by 2*F (since there are two multiplications), but that + // could overflow even the intermediate type, so we have to split the + // algorithm up in two rounds of F bits each. Each round will deplete + // 'num' digit by digit, so after a round we can shift it again. + num <<= F; + ofs -= F; + do_round(); + + num <<= F; + ofs += F; + do_round(); + + return Fixed::from_raw_value(static_cast<B>(res)); +} + +template <typename B, typename I, unsigned int F, bool R> +fixed<B, I, F, R> sqrt(fixed<B, I, F, R> x) noexcept +{ + using Fixed = fixed<B, I, F, R>; + + assert(x >= Fixed(0)); + if (x == Fixed(0)) + { + return x; + } + + // Finding the square root of an integer in base-2, from: + // https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_.28base_2.29 + + // Shift by F first because it's fixed-point. + I num = I{x.raw_value()} << F; + I res = 0; + + // "bit" starts at the greatest power of four that's less than the argument. + for (I bit = I{1} << ((detail::find_highest_bit(x.raw_value()) + F) / 2 * 2); bit != 0; bit >>= 2) + { + const I val = res + bit; + res >>= 1; + if (num >= val) + { + num -= val; + res += bit; + } + } + + // Round the last digit up if necessary + if (num > res) + { + res++; + } + + return Fixed::from_raw_value(static_cast<B>(res)); +} + +template <typename B, typename I, unsigned int F, bool R> +fixed<B, I, F, R> hypot(fixed<B, I, F, R> x, fixed<B, I, F, R> y) noexcept +{ + assert(x != 0 || y != 0); + return sqrt(x*x + y*y); +} + +// +// Trigonometry functions +// + +template <typename B, typename I, unsigned int F, bool R> +fixed<B, I, F, R> sin(fixed<B, I, F, R> x) noexcept +{ + // This sine uses a fifth-order curve-fitting approximation originally + // described by Jasper Vijn on coranac.com which has a worst-case + // relative error of 0.07% (over [-pi:pi]). + using Fixed = fixed<B, I, F, R>; + + // Turn x from [0..2*PI] domain into [0..4] domain + x = fmod(x, Fixed::two_pi()); + x = x / Fixed::half_pi(); + + // Take x modulo one rotation, so [-4..+4]. + if (x < Fixed(0)) { + x += Fixed(4); + } + + int sign = +1; + if (x > Fixed(2)) { + // Reduce domain to [0..2]. + sign = -1; + x -= Fixed(2); + } + + if (x > Fixed(1)) { + // Reduce domain to [0..1]. + x = Fixed(2) - x; + } + + const Fixed x2 = x*x; + return sign * x * (Fixed::pi() - x2*(Fixed::two_pi() - 5 - x2*(Fixed::pi() - 3)))/2; +} + +template <typename B, typename I, unsigned int F, bool R> +inline fixed<B, I, F, R> cos(fixed<B, I, F, R> x) noexcept +{ + using Fixed = fixed<B, I, F, R>; + if (x > Fixed(0)) { // Prevent an overflow due to the addition of π/2 + return sin(x - (Fixed::two_pi() - Fixed::half_pi())); + } else { + return sin(Fixed::half_pi() + x); + } +} + +template <typename B, typename I, unsigned int F, bool R> +inline fixed<B, I, F, R> tan(fixed<B, I, F, R> x) noexcept +{ + auto cx = cos(x); + + // Tangent goes to infinity at 90 and -90 degrees. + // We can't represent that with fixed-point maths. + assert(abs(cx).raw_value() > 1); + + return sin(x) / cx; +} + +namespace detail { + +// Calculates atan(x) assuming that x is in the range [0,1] +template <typename B, typename I, unsigned int F, bool R> +fixed<B, I, F, R> atan_sanitized(fixed<B, I, F, R> x) noexcept +{ + using Fixed = fixed<B, I, F, R>; + assert(x >= Fixed(0) && x <= Fixed(1)); + + constexpr auto fA = Fixed::template from_fixed_point<63>( 716203666280654660ll); // 0.0776509570923569 + constexpr auto fB = Fixed::template from_fixed_point<63>(-2651115102768076601ll); // -0.287434475393028 + constexpr auto fC = Fixed::template from_fixed_point<63>( 9178930894564541004ll); // 0.995181681698119 (PI/4 - A - B) + + const auto xx = x * x; + return ((fA*xx + fB)*xx + fC)*x; +} + +// Calculate atan(y / x), assuming x != 0. +// +// If x is very, very small, y/x can easily overflow the fixed-point range. +// If q = y/x and q > 1, atan(q) would calculate atan(1/q) as intermediate step +// anyway. We can shortcut that here and avoid the loss of information, thus +// improving the accuracy of atan(y/x) for very small x. +template <typename B, typename I, unsigned int F, bool R> +fixed<B, I, F, R> atan_div(fixed<B, I, F, R> y, fixed<B, I, F, R> x) noexcept +{ + using Fixed = fixed<B, I, F, R>; + assert(x != Fixed(0)); + + // Make sure y and x are positive. + // If y / x is negative (when y or x, but not both, are negative), negate the result to + // keep the correct outcome. + if (y < Fixed(0)) { + if (x < Fixed(0)) { + return atan_div(-y, -x); + } + return -atan_div(-y, x); + } + if (x < Fixed(0)) { + return -atan_div(y, -x); + } + assert(y >= Fixed(0)); + assert(x > Fixed(0)); + + if (y > x) { + return Fixed::half_pi() - detail::atan_sanitized(x / y); + } + return detail::atan_sanitized(y / x); +} + +} + +template <typename B, typename I, unsigned int F, bool R> +fixed<B, I, F, R> atan(fixed<B, I, F, R> x) noexcept +{ + using Fixed = fixed<B, I, F, R>; + if (x < Fixed(0)) + { + return -atan(-x); + } + + if (x > Fixed(1)) + { + return Fixed::half_pi() - detail::atan_sanitized(Fixed(1) / x); + } + + return detail::atan_sanitized(x); +} + +template <typename B, typename I, unsigned int F, bool R> +fixed<B, I, F, R> asin(fixed<B, I, F, R> x) noexcept +{ + using Fixed = fixed<B, I, F, R>; + assert(x >= Fixed(-1) && x <= Fixed(+1)); + + const auto yy = Fixed(1) - x * x; + if (yy == Fixed(0)) + { + return copysign(Fixed::half_pi(), x); + } + return detail::atan_div(x, sqrt(yy)); +} + +template <typename B, typename I, unsigned int F, bool R> +fixed<B, I, F, R> acos(fixed<B, I, F, R> x) noexcept +{ + using Fixed = fixed<B, I, F, R>; + assert(x >= Fixed(-1) && x <= Fixed(+1)); + + if (x == Fixed(-1)) + { + return Fixed::pi(); + } + const auto yy = Fixed(1) - x * x; + return Fixed(2)*detail::atan_div(sqrt(yy), Fixed(1) + x); +} + +template <typename B, typename I, unsigned int F, bool R> +fixed<B, I, F, R> atan2(fixed<B, I, F, R> y, fixed<B, I, F, R> x) noexcept +{ + using Fixed = fixed<B, I, F, R>; + if (x == Fixed(0)) + { + assert(y != Fixed(0)); + return (y > Fixed(0)) ? Fixed::half_pi() : -Fixed::half_pi(); + } + + auto ret = detail::atan_div(y, x); + + if (x < Fixed(0)) + { + return (y >= Fixed(0)) ? ret + Fixed::pi() : ret - Fixed::pi(); + } + return ret; +} + +} + +#endif diff --git a/kiss/CHANGELOG b/kiss/CHANGELOG new file mode 100644 index 0000000..2dd3603 --- /dev/null +++ b/kiss/CHANGELOG @@ -0,0 +1,123 @@ +1.3.0 2012-07-18 + removed non-standard malloc.h from kiss_fft.h + + moved -lm to end of link line + + checked various return values + + converted python Numeric code to NumPy + + fixed test of int32_t on 64 bit OS + + added padding in a couple of places to allow SIMD alignment of structs + +1.2.9 2010-05-27 + threadsafe ( including OpenMP ) + + first edition of kissfft.hh the C++ template fft engine + +1.2.8 + Changed memory.h to string.h -- apparently more standard + + Added openmp extensions. This can have fairly linear speedups for larger FFT sizes. + +1.2.7 + Shrank the real-fft memory footprint. Thanks to Galen Seitz. + +1.2.6 (Nov 14, 2006) The "thanks to GenArts" release. + Added multi-dimensional real-optimized FFT, see tools/kiss_fftndr + Thanks go to GenArts, Inc. for sponsoring the development. + +1.2.5 (June 27, 2006) The "release for no good reason" release. + Changed some harmless code to make some compilers' warnings go away. + Added some more digits to pi -- why not. + Added kiss_fft_next_fast_size() function to help people decide how much to pad. + Changed multidimensional test from 8 dimensions to only 3 to avoid testing + problems with fixed point (sorry Buckaroo Banzai). + +1.2.4 (Oct 27, 2005) The "oops, inverse fixed point real fft was borked" release. + Fixed scaling bug for inverse fixed point real fft -- also fixed test code that should've been failing. + Thanks to Jean-Marc Valin for bug report. + + Use sys/types.h for more portable types than short,int,long => int16_t,int32_t,int64_t + If your system does not have these, you may need to define them -- but at least it breaks in a + loud and easily fixable way -- unlike silently using the wrong size type. + + Hopefully tools/psdpng.c is fixed -- thanks to Steve Kellog for pointing out the weirdness. + +1.2.3 (June 25, 2005) The "you want to use WHAT as a sample" release. + Added ability to use 32 bit fixed point samples -- requires a 64 bit intermediate result, a la 'long long' + + Added ability to do 4 FFTs in parallel by using SSE SIMD instructions. This is accomplished by + using the __m128 (vector of 4 floats) as kiss_fft_scalar. Define USE_SIMD to use this. + + I know, I know ... this is drifting a bit from the "kiss" principle, but the speed advantages + make it worth it for some. Also recent gcc makes it SOO easy to use vectors of 4 floats like a POD type. + +1.2.2 (May 6, 2005) The Matthew release + Replaced fixed point division with multiply&shift. Thanks to Jean-Marc Valin for + discussions regarding. Considerable speedup for fixed-point. + + Corrected overflow protection in real fft routines when using fixed point. + Finder's Credit goes to Robert Oschler of robodance for pointing me at the bug. + This also led to the CHECK_OVERFLOW_OP macro. + +1.2.1 (April 4, 2004) + compiles cleanly with just about every -W warning flag under the sun + + reorganized kiss_fft_state so it could be read-only/const. This may be useful for embedded systems + that are willing to predeclare twiddle factors, factorization. + + Fixed C_MUL,S_MUL on 16-bit platforms. + + tmpbuf will only be allocated if input & output buffers are same + scratchbuf will only be allocated for ffts that are not multiples of 2,3,5 + + NOTE: The tmpbuf,scratchbuf changes may require synchronization code for multi-threaded apps. + + +1.2 (Feb 23, 2004) + interface change -- cfg object is forward declaration of struct instead of void* + This maintains type saftey and lets the compiler warn/error about stupid mistakes. + (prompted by suggestion from Erik de Castro Lopo) + + small speed improvements + + added psdpng.c -- sample utility that will create png spectrum "waterfalls" from an input file + ( not terribly useful yet) + +1.1.1 (Feb 1, 2004 ) + minor bug fix -- only affects odd rank, in-place, multi-dimensional FFTs + +1.1 : (Jan 30,2004) + split sample_code/ into test/ and tools/ + + Removed 2-D fft and added N-D fft (arbitrary) + + modified fftutil.c to allow multi-d FFTs + + Modified core fft routine to allow an input stride via kiss_fft_stride() + (eased support of multi-D ffts) + + Added fast convolution filtering (FIR filtering using overlap-scrap method, with tail scrap) + + Add kfc.[ch]: the KISS FFT Cache. It takes care of allocs for you ( suggested by Oscar Lesta ). + +1.0.1 (Dec 15, 2003) + fixed bug that occurred when nfft==1. Thanks to Steven Johnson. + +1.0 : (Dec 14, 2003) + changed kiss_fft function from using a single buffer, to two buffers. + If the same buffer pointer is supplied for both in and out, kiss will + manage the buffer copies. + + added kiss_fft2d and kiss_fftr as separate source files (declarations in kiss_fft.h ) + +0.4 :(Nov 4,2003) optimized for radix 2,3,4,5 + +0.3 :(Oct 28, 2003) woops, version 2 didn't actually factor out any radices other than 2. + Thanks to Steven Johnson for finding this one. + +0.2 :(Oct 27, 2003) added mixed radix, only radix 2,4 optimized versions + +0.1 :(May 19 2003) initial release, radix 2 only diff --git a/kiss/COPYING b/kiss/COPYING new file mode 100644 index 0000000..6b4b622 --- /dev/null +++ b/kiss/COPYING @@ -0,0 +1,11 @@ +Copyright (c) 2003-2010 Mark Borgerding . All rights reserved. + +KISS FFT is provided under: + + SPDX-License-Identifier: BSD-3-Clause + +Being under the terms of the BSD 3-clause "New" or "Revised" License, +according with: + + LICENSES/BSD-3-Clause + diff --git a/kiss/README.md b/kiss/README.md new file mode 100644 index 0000000..1138a0c --- /dev/null +++ b/kiss/README.md @@ -0,0 +1,245 @@ +# KISS FFT [](https://travis-ci.com/mborgerding/kissfft) + +KISS FFT - A mixed-radix Fast Fourier Transform based up on the principle, +"Keep It Simple, Stupid." + +There are many great fft libraries already around. Kiss FFT is not trying +to be better than any of them. It only attempts to be a reasonably efficient, +moderately useful FFT that can use fixed or floating data types and can be +incorporated into someone's C program in a few minutes with trivial licensing. + +## USAGE: + +The basic usage for 1-d complex FFT is: + +```c + #include "kiss_fft.h" + kiss_fft_cfg cfg = kiss_fft_alloc( nfft ,is_inverse_fft ,0,0 ); + while ... + + ... // put kth sample in cx_in[k].r and cx_in[k].i + + kiss_fft( cfg , cx_in , cx_out ); + + ... // transformed. DC is in cx_out[0].r and cx_out[0].i + + kiss_fft_free(cfg); +``` + - **Note**: frequency-domain data is stored from dc up to 2pi. + so cx_out[0] is the dc bin of the FFT + and cx_out[nfft/2] is the Nyquist bin (if exists) + +Declarations are in "kiss_fft.h", along with a brief description of the +functions you'll need to use. + +Code definitions for 1d complex FFTs are in kiss_fft.c. + +You can do other cool stuff with the extras you'll find in tools/ +> - multi-dimensional FFTs +> - real-optimized FFTs (returns the positive half-spectrum: + (nfft/2+1) complex frequency bins) +> - fast convolution FIR filtering (not available for fixed point) +> - spectrum image creation + +The core fft and most tools/ code can be compiled to use float, double, + Q15 short or Q31 samples. The default is float. + +## BUILDING: + +There are two functionally-equivalent build systems supported by kissfft: + + - Make (traditional Makefiles for Unix / Linux systems) + - CMake (more modern and feature-rich build system developed by Kitware) + +To build kissfft, the following build environment can be used: + + - GNU build environment with GCC, Clang and GNU Make or CMake (>= 3.6) + - Microsoft Visual C++ (MSVC) with CMake (>= 3.6) + +Additional libraries required to build and test kissfft include: + + - libpng for psdpng tool, + - libfftw3 to validate kissfft results against it, + - python 2/3 with Numpy to validate kissfft results against it. + - OpenMP supported by GCC, Clang or MSVC for multi-core FFT transformations + +Environments like Cygwin and MinGW can be highly likely used to build kissfft +targeting Windows platform, but no tests were performed to the date. + +Both Make and CMake builds are easily configurable: + + - `KISSFFT_DATATYPE=<datatype>` (for Make) or `-DKISSFFT_DATATYPE=<datatype>` + (for CMake) denote the principal datatype used by kissfft. It can be one + of the following: + + - float (default) + - double + - int16_t + - int32_t + - SIMD (requires SSE instruction set support on target CPU) + + - `KISSFFT_OPENMP=1` (for Make) or `-DKISSFFT_OPENMP=ON` (for CMake) builds kissfft + with OpenMP support. Please note that a supported compiler is required and this + option is turned off by default. + + - `KISSFFT_STATIC=1` (for Make) or `-DKISSFFT_STATIC=ON` (for CMake) instructs + the builder to create static library ('.lib' for Windows / '.a' for Unix or Linux). + By default, this option is turned off and the shared library is created + ('.dll' for Windows, '.so' for Linux or Unix, '.dylib' for Mac OSX) + + - `-DKISSFFT_TEST=OFF` (for CMake) disables building tests for kissfft. On Make, + building tests is done separately by 'make testall' or 'make testsingle', so + no specific setting is required. + + - `KISSFFT_TOOLS=0` (for Make) or `-DKISSFFT_TOOLS=OFF` (for CMake) builds kissfft + without command-line tools like 'fastconv'. By default the tools are built. + + - `KISSFFT_USE_ALLOCA=1` (for Make) or `-DKISSFFT_USE_ALLOCA=ON` (for CMake) + build kissfft with 'alloca' usage instead of 'malloc' / 'free'. + + - `PREFIX=/full/path/to/installation/prefix/directory` (for Make) or + `-DCMAKE_INSTALL_PREFIX=/full/path/to/installation/prefix/directory` (for CMake) + specifies the prefix directory to install kissfft into. + +For example, to build kissfft as a static library with 'int16_t' datatype and +OpenMP support using Make, run the command from kissfft source tree: + +``` +make KISSFFT_DATATYPE=int16_t KISSFFT_STATIC=1 KISSFFT_OPENMP=1 all +``` + +The same configuration for CMake is: + +``` +mkdir build && cd build +cmake -DKISSFFT_DATATYPE=int16_t -DKISSFFT_STATIC=ON -DKISSFFT_OPENMP=ON .. +make all +``` + +To specify '/tmp/1234' as installation prefix directory, run: + + +``` +make PREFIX=/tmp/1234 KISSFFT_DATATYPE=int16_t KISSFFT_STATIC=1 KISSFFT_OPENMP=1 install +``` + +or + +``` +mkdir build && cd build +cmake -DCMAKE_INSTALL_PREFIX=/tmp/1234 -DKISSFFT_DATATYPE=int16_t -DKISSFFT_STATIC=ON -DKISSFFT_OPENMP=ON .. +make all +make install +``` + +## TESTING: + +To validate the build configured as an example above, run the following command from +kissfft source tree: + +``` +make KISSFFT_DATATYPE=int16_t KISSFFT_STATIC=1 KISSFFT_OPENMP=1 testsingle +``` + +if using Make, or: + +``` +make test +``` + +if using CMake. + +To test all possible build configurations, please run an extended testsuite from +kissfft source tree: + +``` +sh test/kissfft-testsuite.sh +``` + +Please note that the extended testsuite takes around 20-40 minutes depending on device +it runs on. This testsuite is useful for reporting bugs or testing the pull requests. + +## BACKGROUND + +I started coding this because I couldn't find a fixed point FFT that didn't +use assembly code. I started with floating point numbers so I could get the +theory straight before working on fixed point issues. In the end, I had a +little bit of code that could be recompiled easily to do ffts with short, float +or double (other types should be easy too). + +Once I got my FFT working, I was curious about the speed compared to +a well respected and highly optimized fft library. I don't want to criticize +this great library, so let's call it FFT_BRANDX. +During this process, I learned: + +> 1. FFT_BRANDX has more than 100K lines of code. The core of kiss_fft is about 500 lines (cpx 1-d). +> 2. It took me an embarrassingly long time to get FFT_BRANDX working. +> 3. A simple program using FFT_BRANDX is 522KB. A similar program using kiss_fft is 18KB (without optimizing for size). +> 4. FFT_BRANDX is roughly twice as fast as KISS FFT in default mode. + +It is wonderful that free, highly optimized libraries like FFT_BRANDX exist. +But such libraries carry a huge burden of complexity necessary to extract every +last bit of performance. + +**Sometimes simpler is better, even if it's not better.** + +## FREQUENTLY ASKED QUESTIONS: +> Q: Can I use kissfft in a project with a ___ license?</br> +> A: Yes. See LICENSE below. + +> Q: Why don't I get the output I expect?</br> +> A: The two most common causes of this are +> 1) scaling : is there a constant multiplier between what you got and what you want? +> 2) mixed build environment -- all code must be compiled with same preprocessor +> definitions for FIXED_POINT and kiss_fft_scalar + +> Q: Will you write/debug my code for me?</br> +> A: Probably not unless you pay me. I am happy to answer pointed and topical questions, but +> I may refer you to a book, a forum, or some other resource. + + +## PERFORMANCE + (on Athlon XP 2100+, with gcc 2.96, float data type) + +Kiss performed 10000 1024-pt cpx ffts in .63 s of cpu time. +For comparison, it took md5sum twice as long to process the same amount of data. +Transforming 5 minutes of CD quality audio takes less than a second (nfft=1024). + +**DO NOT:** +- use Kiss if you need the Fastest Fourier Transform in the World +- ask me to add features that will bloat the code + +## UNDER THE HOOD + +Kiss FFT uses a time decimation, mixed-radix, out-of-place FFT. If you give it an input buffer +and output buffer that are the same, a temporary buffer will be created to hold the data. + +No static data is used. The core routines of kiss_fft are thread-safe (but not all of the tools directory).[ + +No scaling is done for the floating point version (for speed). +Scaling is done both ways for the fixed-point version (for overflow prevention). + +Optimized butterflies are used for factors 2,3,4, and 5. + +The real (i.e. not complex) optimization code only works for even length ffts. It does two half-length +FFTs in parallel (packed into real&imag), and then combines them via twiddling. The result is +nfft/2+1 complex frequency bins from DC to Nyquist. If you don't know what this means, search the web. + +The fast convolution filtering uses the overlap-scrap method, slightly +modified to put the scrap at the tail. + +## LICENSE + Revised BSD License, see COPYING for verbiage. + Basically, "free to use&change, give credit where due, no guarantees" + Note this license is compatible with GPL at one end of the spectrum and closed, commercial software at + the other end. See http://www.fsf.org/licensing/licenses + +## TODO + - Add real optimization for odd length FFTs + - Document/revisit the input/output fft scaling + - Make doc describing the overlap (tail) scrap fast convolution filtering in kiss_fastfir.c + - Test all the ./tools/ code with fixed point (kiss_fastfir.c doesn't work, maybe others) + +## AUTHOR + Mark Borgerding + Mark@Borgerding.net diff --git a/kiss/_kiss_fft_guts.h b/kiss/_kiss_fft_guts.h new file mode 100644 index 0000000..4bd8d1c --- /dev/null +++ b/kiss/_kiss_fft_guts.h @@ -0,0 +1,167 @@ +/* + * Copyright (c) 2003-2010, Mark Borgerding. All rights reserved. + * This file is part of KISS FFT - https://github.com/mborgerding/kissfft + * + * SPDX-License-Identifier: BSD-3-Clause + * See COPYING file for more information. + */ + +/* kiss_fft.h + defines kiss_fft_scalar as either short or a float type + and defines + typedef struct { kiss_fft_scalar r; kiss_fft_scalar i; }kiss_fft_cpx; */ + +#ifndef _kiss_fft_guts_h +#define _kiss_fft_guts_h + +#include "kiss_fft.h" +#include "kiss_fft_log.h" +#include <limits.h> + +#define MAXFACTORS 32 +/* e.g. an fft of length 128 has 4 factors + as far as kissfft is concerned + 4*4*4*2 + */ + +struct kiss_fft_state{ + int nfft; + int inverse; + int factors[2*MAXFACTORS]; + kiss_fft_cpx twiddles[1]; +}; + +/* + Explanation of macros dealing with complex math: + + C_MUL(m,a,b) : m = a*b + C_FIXDIV( c , div ) : if a fixed point impl., c /= div. noop otherwise + C_SUB( res, a,b) : res = a - b + C_SUBFROM( res , a) : res -= a + C_ADDTO( res , a) : res += a + * */ +#ifdef FIXED_POINT +#include <stdint.h> +#if (FIXED_POINT==32) +# define FRACBITS 31 +# define SAMPPROD int64_t +#define SAMP_MAX INT32_MAX +#define SAMP_MIN INT32_MIN +#else +# define FRACBITS 15 +# define SAMPPROD int32_t +#define SAMP_MAX INT16_MAX +#define SAMP_MIN INT16_MIN +#endif + +#if defined(CHECK_OVERFLOW) +# define CHECK_OVERFLOW_OP(a,op,b) \ + if ( (SAMPPROD)(a) op (SAMPPROD)(b) > SAMP_MAX || (SAMPPROD)(a) op (SAMPPROD)(b) < SAMP_MIN ) { \ + KISS_FFT_WARNING("overflow (%d " #op" %d) = %ld", (a),(b),(SAMPPROD)(a) op (SAMPPROD)(b)); } +#endif + + +# define smul(a,b) ( (SAMPPROD)(a)*(b) ) +# define sround( x ) (kiss_fft_scalar)( ( (x) + (1<<(FRACBITS-1)) ) >> FRACBITS ) + +# define S_MUL(a,b) sround( smul(a,b) ) + +# define C_MUL(m,a,b) \ + do{ (m).r = sround( smul((a).r,(b).r) - smul((a).i,(b).i) ); \ + (m).i = sround( smul((a).r,(b).i) + smul((a).i,(b).r) ); }while(0) + +# define DIVSCALAR(x,k) \ + (x) = sround( smul( x, SAMP_MAX/k ) ) + +# define C_FIXDIV(c,div) \ + do { DIVSCALAR( (c).r , div); \ + DIVSCALAR( (c).i , div); }while (0) + +# define C_MULBYSCALAR( c, s ) \ + do{ (c).r = sround( smul( (c).r , s ) ) ;\ + (c).i = sround( smul( (c).i , s ) ) ; }while(0) + +#else /* not FIXED_POINT*/ + +# define S_MUL(a,b) ( (a)*(b) ) +#define C_MUL(m,a,b) \ + do{ (m).r = (a).r*(b).r - (a).i*(b).i;\ + (m).i = (a).r*(b).i + (a).i*(b).r; }while(0) +# define C_FIXDIV(c,div) /* NOOP */ +# define C_MULBYSCALAR( c, s ) \ + do{ (c).r *= (s);\ + (c).i *= (s); }while(0) +#endif + +#ifndef CHECK_OVERFLOW_OP +# define CHECK_OVERFLOW_OP(a,op,b) /* noop */ +#endif + +#define C_ADD( res, a,b)\ + do { \ + CHECK_OVERFLOW_OP((a).r,+,(b).r)\ + CHECK_OVERFLOW_OP((a).i,+,(b).i)\ + (res).r=(a).r+(b).r; (res).i=(a).i+(b).i; \ + }while(0) +#define C_SUB( res, a,b)\ + do { \ + CHECK_OVERFLOW_OP((a).r,-,(b).r)\ + CHECK_OVERFLOW_OP((a).i,-,(b).i)\ + (res).r=(a).r-(b).r; (res).i=(a).i-(b).i; \ + }while(0) +#define C_ADDTO( res , a)\ + do { \ + CHECK_OVERFLOW_OP((res).r,+,(a).r)\ + CHECK_OVERFLOW_OP((res).i,+,(a).i)\ + (res).r += (a).r; (res).i += (a).i;\ + }while(0) + +#define C_SUBFROM( res , a)\ + do {\ + CHECK_OVERFLOW_OP((res).r,-,(a).r)\ + CHECK_OVERFLOW_OP((res).i,-,(a).i)\ + (res).r -= (a).r; (res).i -= (a).i; \ + }while(0) + + +#ifdef FIXED_POINT +# define KISS_FFT_COS(phase) floor(.5+SAMP_MAX * cos (phase)) +# define KISS_FFT_SIN(phase) floor(.5+SAMP_MAX * sin (phase)) +# define HALF_OF(x) ((x)>>1) +#elif defined(USE_SIMD) +# define KISS_FFT_COS(phase) _mm_set1_ps( cos(phase) ) +# define KISS_FFT_SIN(phase) _mm_set1_ps( sin(phase) ) +# define HALF_OF(x) ((x)*_mm_set1_ps(.5)) +#else +# define KISS_FFT_COS(phase) (kiss_fft_scalar) cos(phase) +# define KISS_FFT_SIN(phase) (kiss_fft_scalar) sin(phase) +# define HALF_OF(x) ((x)*((kiss_fft_scalar).5)) +#endif + +#define kf_cexp(x,phase) \ + do{ \ + (x)->r = KISS_FFT_COS(phase);\ + (x)->i = KISS_FFT_SIN(phase);\ + }while(0) + + +/* a debugging function */ +#define pcpx(c)\ + KISS_FFT_DEBUG("%g + %gi\n",(double)((c)->r),(double)((c)->i)) + + +#ifdef KISS_FFT_USE_ALLOCA +// define this to allow use of alloca instead of malloc for temporary buffers +// Temporary buffers are used in two case: +// 1. FFT sizes that have "bad" factors. i.e. not 2,3 and 5 +// 2. "in-place" FFTs. Notice the quotes, since kissfft does not really do an in-place transform. +#include <alloca.h> +#define KISS_FFT_TMP_ALLOC(nbytes) alloca(nbytes) +#define KISS_FFT_TMP_FREE(ptr) +#else +#define KISS_FFT_TMP_ALLOC(nbytes) KISS_FFT_MALLOC(nbytes) +#define KISS_FFT_TMP_FREE(ptr) KISS_FFT_FREE(ptr) +#endif + +#endif /* _kiss_fft_guts_h */ + diff --git a/kiss/kfc.c b/kiss/kfc.c new file mode 100644 index 0000000..a405d9b --- /dev/null +++ b/kiss/kfc.c @@ -0,0 +1,109 @@ +/* + * Copyright (c) 2003-2004, Mark Borgerding. All rights reserved. + * This file is part of KISS FFT - https://github.com/mborgerding/kissfft + * + * SPDX-License-Identifier: BSD-3-Clause + * See COPYING file for more information. + */ + +#include "kfc.h" + +typedef struct cached_fft *kfc_cfg; + +struct cached_fft +{ + int nfft; + int inverse; + kiss_fft_cfg cfg; + kfc_cfg next; +}; + +static kfc_cfg cache_root=NULL; +static int ncached=0; + +static kiss_fft_cfg find_cached_fft(int nfft,int inverse) +{ + size_t len; + kfc_cfg cur=cache_root; + kfc_cfg prev=NULL; + while ( cur ) { + if ( cur->nfft == nfft && inverse == cur->inverse ) + break;/*found the right node*/ + prev = cur; + cur = prev->next; + } + if (cur== NULL) { + /* no cached node found, need to create a new one*/ + kiss_fft_alloc(nfft,inverse,0,&len); +#ifdef USE_SIMD + int padding = (16-sizeof(struct cached_fft)) & 15; + // make sure the cfg aligns on a 16 byte boundary + len += padding; +#endif + cur = (kfc_cfg)KISS_FFT_MALLOC((sizeof(struct cached_fft) + len )); + if (cur == NULL) + return NULL; + cur->cfg = (kiss_fft_cfg)(cur+1); +#ifdef USE_SIMD + cur->cfg = (kiss_fft_cfg) ((char*)(cur+1)+padding); +#endif + kiss_fft_alloc(nfft,inverse,cur->cfg,&len); + cur->nfft=nfft; + cur->inverse=inverse; + cur->next = NULL; + if ( prev ) + prev->next = cur; + else + cache_root = cur; + ++ncached; + } + return cur->cfg; +} + +void kfc_cleanup(void) +{ + kfc_cfg cur=cache_root; + kfc_cfg next=NULL; + while (cur){ + next = cur->next; + free(cur); + cur=next; + } + ncached=0; + cache_root = NULL; +} +void kfc_fft(int nfft, const kiss_fft_cpx * fin,kiss_fft_cpx * fout) +{ + kiss_fft( find_cached_fft(nfft,0),fin,fout ); +} + +void kfc_ifft(int nfft, const kiss_fft_cpx * fin,kiss_fft_cpx * fout) +{ + kiss_fft( find_cached_fft(nfft,1),fin,fout ); +} + +#ifdef KFC_TEST +static void check(int nc) +{ + if (ncached != nc) { + fprintf(stderr,"ncached should be %d,but it is %d\n",nc,ncached); + exit(1); + } +} + +int main(void) +{ + kiss_fft_cpx buf1[1024],buf2[1024]; + memset(buf1,0,sizeof(buf1)); + check(0); + kfc_fft(512,buf1,buf2); + check(1); + kfc_fft(512,buf1,buf2); + check(1); + kfc_ifft(512,buf1,buf2); + check(2); + kfc_cleanup(); + check(0); + return 0; +} +#endif diff --git a/kiss/kfc.h b/kiss/kfc.h new file mode 100644 index 0000000..d7d8c1b --- /dev/null +++ b/kiss/kfc.h @@ -0,0 +1,54 @@ +/* + * Copyright (c) 2003-2004, Mark Borgerding. All rights reserved. + * This file is part of KISS FFT - https://github.com/mborgerding/kissfft + * + * SPDX-License-Identifier: BSD-3-Clause + * See COPYING file for more information. + */ + +#ifndef KFC_H +#define KFC_H +#include "kiss_fft.h" + +#ifdef __cplusplus +extern "C" { +#endif + +/* +KFC -- Kiss FFT Cache + +Not needing to deal with kiss_fft_alloc and a config +object may be handy for a lot of programs. + +KFC uses the underlying KISS FFT functions, but caches the config object. +The first time kfc_fft or kfc_ifft for a given FFT size, the cfg +object is created for it. All subsequent calls use the cached +configuration object. + +NOTE: +You should probably not use this if your program will be using a lot +of various sizes of FFTs. There is a linear search through the +cached objects. If you are only using one or two FFT sizes, this +will be negligible. Otherwise, you may want to use another method +of managing the cfg objects. + + There is no automated cleanup of the cached objects. This could lead +to large memory usage in a program that uses a lot of *DIFFERENT* +sized FFTs. If you want to force all cached cfg objects to be freed, +call kfc_cleanup. + + */ + +/*forward complex FFT */ +void KISS_FFT_API kfc_fft(int nfft, const kiss_fft_cpx * fin,kiss_fft_cpx * fout); +/*reverse complex FFT */ +void KISS_FFT_API kfc_ifft(int nfft, const kiss_fft_cpx * fin,kiss_fft_cpx * fout); + +/*free all cached objects*/ +void KISS_FFT_API kfc_cleanup(void); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/kiss/kiss_fft.c b/kiss/kiss_fft.c new file mode 100644 index 0000000..58c24a0 --- /dev/null +++ b/kiss/kiss_fft.c @@ -0,0 +1,420 @@ +/* + * Copyright (c) 2003-2010, Mark Borgerding. All rights reserved. + * This file is part of KISS FFT - https://github.com/mborgerding/kissfft + * + * SPDX-License-Identifier: BSD-3-Clause + * See COPYING file for more information. + */ + + +#include "_kiss_fft_guts.h" +/* The guts header contains all the multiplication and addition macros that are defined for + fixed or floating point complex numbers. It also delares the kf_ internal functions. + */ + +static void kf_bfly2( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + int m + ) +{ + kiss_fft_cpx * Fout2; + kiss_fft_cpx * tw1 = st->twiddles; + kiss_fft_cpx t; + Fout2 = Fout + m; + do{ + C_FIXDIV(*Fout,2); C_FIXDIV(*Fout2,2); + + C_MUL (t, *Fout2 , *tw1); + tw1 += fstride; + C_SUB( *Fout2 , *Fout , t ); + C_ADDTO( *Fout , t ); + ++Fout2; + ++Fout; + }while (--m); +} + +static void kf_bfly4( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + const size_t m + ) +{ + kiss_fft_cpx *tw1,*tw2,*tw3; + kiss_fft_cpx scratch[6]; + size_t k=m; + const size_t m2=2*m; + const size_t m3=3*m; + + + tw3 = tw2 = tw1 = st->twiddles; + + do { + C_FIXDIV(*Fout,4); C_FIXDIV(Fout[m],4); C_FIXDIV(Fout[m2],4); C_FIXDIV(Fout[m3],4); + + C_MUL(scratch[0],Fout[m] , *tw1 ); + C_MUL(scratch[1],Fout[m2] , *tw2 ); + C_MUL(scratch[2],Fout[m3] , *tw3 ); + + C_SUB( scratch[5] , *Fout, scratch[1] ); + C_ADDTO(*Fout, scratch[1]); + C_ADD( scratch[3] , scratch[0] , scratch[2] ); + C_SUB( scratch[4] , scratch[0] , scratch[2] ); + C_SUB( Fout[m2], *Fout, scratch[3] ); + tw1 += fstride; + tw2 += fstride*2; + tw3 += fstride*3; + C_ADDTO( *Fout , scratch[3] ); + + if(st->inverse) { + Fout[m].r = scratch[5].r - scratch[4].i; + Fout[m].i = scratch[5].i + scratch[4].r; + Fout[m3].r = scratch[5].r + scratch[4].i; + Fout[m3].i = scratch[5].i - scratch[4].r; + }else{ + Fout[m].r = scratch[5].r + scratch[4].i; + Fout[m].i = scratch[5].i - scratch[4].r; + Fout[m3].r = scratch[5].r - scratch[4].i; + Fout[m3].i = scratch[5].i + scratch[4].r; + } + ++Fout; + }while(--k); +} + +static void kf_bfly3( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + size_t m + ) +{ + size_t k=m; + const size_t m2 = 2*m; + kiss_fft_cpx *tw1,*tw2; + kiss_fft_cpx scratch[5]; + kiss_fft_cpx epi3; + epi3 = st->twiddles[fstride*m]; + + tw1=tw2=st->twiddles; + + do{ + C_FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3); + + C_MUL(scratch[1],Fout[m] , *tw1); + C_MUL(scratch[2],Fout[m2] , *tw2); + + C_ADD(scratch[3],scratch[1],scratch[2]); + C_SUB(scratch[0],scratch[1],scratch[2]); + tw1 += fstride; + tw2 += fstride*2; + + Fout[m].r = Fout->r - HALF_OF(scratch[3].r); + Fout[m].i = Fout->i - HALF_OF(scratch[3].i); + + C_MULBYSCALAR( scratch[0] , epi3.i ); + + C_ADDTO(*Fout,scratch[3]); + + Fout[m2].r = Fout[m].r + scratch[0].i; + Fout[m2].i = Fout[m].i - scratch[0].r; + + Fout[m].r -= scratch[0].i; + Fout[m].i += scratch[0].r; + + ++Fout; + }while(--k); +} + +static void kf_bfly5( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + int m + ) +{ + kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4; + int u; + kiss_fft_cpx scratch[13]; + kiss_fft_cpx * twiddles = st->twiddles; + kiss_fft_cpx *tw; + kiss_fft_cpx ya,yb; + ya = twiddles[fstride*m]; + yb = twiddles[fstride*2*m]; + + Fout0=Fout; + Fout1=Fout0+m; + Fout2=Fout0+2*m; + Fout3=Fout0+3*m; + Fout4=Fout0+4*m; + + tw=st->twiddles; + for ( u=0; u<m; ++u ) { + C_FIXDIV( *Fout0,5); C_FIXDIV( *Fout1,5); C_FIXDIV( *Fout2,5); C_FIXDIV( *Fout3,5); C_FIXDIV( *Fout4,5); + scratch[0] = *Fout0; + + C_MUL(scratch[1] ,*Fout1, tw[u*fstride]); + C_MUL(scratch[2] ,*Fout2, tw[2*u*fstride]); + C_MUL(scratch[3] ,*Fout3, tw[3*u*fstride]); + C_MUL(scratch[4] ,*Fout4, tw[4*u*fstride]); + + C_ADD( scratch[7],scratch[1],scratch[4]); + C_SUB( scratch[10],scratch[1],scratch[4]); + C_ADD( scratch[8],scratch[2],scratch[3]); + C_SUB( scratch[9],scratch[2],scratch[3]); + + Fout0->r += scratch[7].r + scratch[8].r; + Fout0->i += scratch[7].i + scratch[8].i; + + scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r); + scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r); + + scratch[6].r = S_MUL(scratch[10].i,ya.i) + S_MUL(scratch[9].i,yb.i); + scratch[6].i = -S_MUL(scratch[10].r,ya.i) - S_MUL(scratch[9].r,yb.i); + + C_SUB(*Fout1,scratch[5],scratch[6]); + C_ADD(*Fout4,scratch[5],scratch[6]); + + scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r); + scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r); + scratch[12].r = - S_MUL(scratch[10].i,yb.i) + S_MUL(scratch[9].i,ya.i); + scratch[12].i = S_MUL(scratch[10].r,yb.i) - S_MUL(scratch[9].r,ya.i); + + C_ADD(*Fout2,scratch[11],scratch[12]); + C_SUB(*Fout3,scratch[11],scratch[12]); + + ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4; + } +} + +/* perform the butterfly for one stage of a mixed radix FFT */ +static void kf_bfly_generic( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + int m, + int p + ) +{ + int u,k,q1,q; + kiss_fft_cpx * twiddles = st->twiddles; + kiss_fft_cpx t; + int Norig = st->nfft; + + kiss_fft_cpx * scratch = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC(sizeof(kiss_fft_cpx)*p); + if (scratch == NULL){ + KISS_FFT_ERROR("Memory allocation failed."); + return; + } + + for ( u=0; u<m; ++u ) { + k=u; + for ( q1=0 ; q1<p ; ++q1 ) { + scratch[q1] = Fout[ k ]; + C_FIXDIV(scratch[q1],p); + k += m; + } + + k=u; + for ( q1=0 ; q1<p ; ++q1 ) { + int twidx=0; + Fout[ k ] = scratch[0]; + for (q=1;q<p;++q ) { + twidx += fstride * k; + if (twidx>=Norig) twidx-=Norig; + C_MUL(t,scratch[q] , twiddles[twidx] ); + C_ADDTO( Fout[ k ] ,t); + } + k += m; + } + } + KISS_FFT_TMP_FREE(scratch); +} + +static +void kf_work( + kiss_fft_cpx * Fout, + const kiss_fft_cpx * f, + const size_t fstride, + int in_stride, + int * factors, + const kiss_fft_cfg st + ) +{ + kiss_fft_cpx * Fout_beg=Fout; + const int p=*factors++; /* the radix */ + const int m=*factors++; /* stage's fft length/p */ + const kiss_fft_cpx * Fout_end = Fout + p*m; + +#ifdef _OPENMP + // use openmp extensions at the + // top-level (not recursive) + if (fstride==1 && p<=5 && m!=1) + { + int k; + + // execute the p different work units in different threads +# pragma omp parallel for + for (k=0;k<p;++k) + kf_work( Fout +k*m, f+ fstride*in_stride*k,fstride*p,in_stride,factors,st); + // all threads have joined by this point + + switch (p) { + case 2: kf_bfly2(Fout,fstride,st,m); break; + case 3: kf_bfly3(Fout,fstride,st,m); break; + case 4: kf_bfly4(Fout,fstride,st,m); break; + case 5: kf_bfly5(Fout,fstride,st,m); break; + default: kf_bfly_generic(Fout,fstride,st,m,p); break; + } + return; + } +#endif + + if (m==1) { + do{ + *Fout = *f; + f += fstride*in_stride; + }while(++Fout != Fout_end ); + }else{ + do{ + // recursive call: + // DFT of size m*p performed by doing + // p instances of smaller DFTs of size m, + // each one takes a decimated version of the input + kf_work( Fout , f, fstride*p, in_stride, factors,st); + f += fstride*in_stride; + }while( (Fout += m) != Fout_end ); + } + + Fout=Fout_beg; + + // recombine the p smaller DFTs + switch (p) { + case 2: kf_bfly2(Fout,fstride,st,m); break; + case 3: kf_bfly3(Fout,fstride,st,m); break; + case 4: kf_bfly4(Fout,fstride,st,m); break; + case 5: kf_bfly5(Fout,fstride,st,m); break; + default: kf_bfly_generic(Fout,fstride,st,m,p); break; + } +} + +/* facbuf is populated by p1,m1,p2,m2, ... + where + p[i] * m[i] = m[i-1] + m0 = n */ +static +void kf_factor(int n,int * facbuf) +{ + int p=4; + double floor_sqrt; + floor_sqrt = floor( sqrt((double)n) ); + + /*factor out powers of 4, powers of 2, then any remaining primes */ + do { + while (n % p) { + switch (p) { + case 4: p = 2; break; + case 2: p = 3; break; + default: p += 2; break; + } + if (p > floor_sqrt) + p = n; /* no more factors, skip to end */ + } + n /= p; + *facbuf++ = p; + *facbuf++ = n; + } while (n > 1); +} + +/* + * + * User-callable function to allocate all necessary storage space for the fft. + * + * The return value is a contiguous block of memory, allocated with malloc. As such, + * It can be freed with free(), rather than a kiss_fft-specific function. + * */ +kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem ) +{ + KISS_FFT_ALIGN_CHECK(mem) + + kiss_fft_cfg st=NULL; + size_t memneeded = KISS_FFT_ALIGN_SIZE_UP(sizeof(struct kiss_fft_state) + + sizeof(kiss_fft_cpx)*(nfft-1)); /* twiddle factors*/ + + if ( lenmem==NULL ) { + st = ( kiss_fft_cfg)KISS_FFT_MALLOC( memneeded ); + }else{ + if (mem != NULL && *lenmem >= memneeded) + st = (kiss_fft_cfg)mem; + *lenmem = memneeded; + } + if (st) { + int i; + st->nfft=nfft; + st->inverse = inverse_fft; + + for (i=0;i<nfft;++i) { + const double pi=3.141592653589793238462643383279502884197169399375105820974944; + double phase = -2*pi*i / nfft; + if (st->inverse) + phase *= -1; + kf_cexp(st->twiddles+i, phase ); + } + + kf_factor(nfft,st->factors); + } + return st; +} + + +void kiss_fft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride) +{ + if (fin == fout) { + //NOTE: this is not really an in-place FFT algorithm. + //It just performs an out-of-place FFT into a temp buffer + if (fout == NULL){ + KISS_FFT_ERROR("fout buffer NULL."); + return; + } + + kiss_fft_cpx * tmpbuf = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC( sizeof(kiss_fft_cpx)*st->nfft); + if (tmpbuf == NULL){ + KISS_FFT_ERROR("Memory allocation error."); + return; + } + + + + kf_work(tmpbuf,fin,1,in_stride, st->factors,st); + memcpy(fout,tmpbuf,sizeof(kiss_fft_cpx)*st->nfft); + KISS_FFT_TMP_FREE(tmpbuf); + }else{ + kf_work( fout, fin, 1,in_stride, st->factors,st ); + } +} + +void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout) +{ + kiss_fft_stride(cfg,fin,fout,1); +} + + +void kiss_fft_cleanup(void) +{ + // nothing needed any more +} + +int kiss_fft_next_fast_size(int n) +{ + while(1) { + int m=n; + while ( (m%2) == 0 ) m/=2; + while ( (m%3) == 0 ) m/=3; + while ( (m%5) == 0 ) m/=5; + if (m<=1) + break; /* n is completely factorable by twos, threes, and fives */ + n++; + } + return n; +} diff --git a/kiss/kiss_fft.h b/kiss/kiss_fft.h new file mode 100644 index 0000000..dce1034 --- /dev/null +++ b/kiss/kiss_fft.h @@ -0,0 +1,160 @@ +/* + * Copyright (c) 2003-2010, Mark Borgerding. All rights reserved. + * This file is part of KISS FFT - https://github.com/mborgerding/kissfft + * + * SPDX-License-Identifier: BSD-3-Clause + * See COPYING file for more information. + */ + +#ifndef KISS_FFT_H +#define KISS_FFT_H + +#include <stdlib.h> +#include <stdio.h> +#include <math.h> +#include <string.h> + +// Define KISS_FFT_SHARED macro to properly export symbols +#ifdef KISS_FFT_SHARED +# ifdef _WIN32 +# ifdef KISS_FFT_BUILD +# define KISS_FFT_API __declspec(dllexport) +# else +# define KISS_FFT_API __declspec(dllimport) +# endif +# else +# define KISS_FFT_API __attribute__ ((visibility ("default"))) +# endif +#else +# define KISS_FFT_API +#endif + +#ifdef __cplusplus +extern "C" { +#endif + +/* + ATTENTION! + If you would like a : + -- a utility that will handle the caching of fft objects + -- real-only (no imaginary time component ) FFT + -- a multi-dimensional FFT + -- a command-line utility to perform ffts + -- a command-line utility to perform fast-convolution filtering + + Then see kfc.h kiss_fftr.h kiss_fftnd.h fftutil.c kiss_fastfir.c + in the tools/ directory. +*/ + +/* User may override KISS_FFT_MALLOC and/or KISS_FFT_FREE. */ +#ifdef USE_SIMD +# include <xmmintrin.h> +# define kiss_fft_scalar __m128 +# ifndef KISS_FFT_MALLOC +# define KISS_FFT_MALLOC(nbytes) _mm_malloc(nbytes,16) +# define KISS_FFT_ALIGN_CHECK(ptr) +# define KISS_FFT_ALIGN_SIZE_UP(size) ((size + 15UL) & ~0xFUL) +# endif +# ifndef KISS_FFT_FREE +# define KISS_FFT_FREE _mm_free +# endif +#else +# define KISS_FFT_ALIGN_CHECK(ptr) +# define KISS_FFT_ALIGN_SIZE_UP(size) (size) +# ifndef KISS_FFT_MALLOC +# define KISS_FFT_MALLOC malloc +# endif +# ifndef KISS_FFT_FREE +# define KISS_FFT_FREE free +# endif +#endif + + +#ifdef FIXED_POINT +#include <stdint.h> +# if (FIXED_POINT == 32) +# define kiss_fft_scalar int32_t +# else +# define kiss_fft_scalar int16_t +# endif +#else +# ifndef kiss_fft_scalar +/* default is float */ +# define kiss_fft_scalar float +# endif +#endif + +typedef struct { + kiss_fft_scalar r; + kiss_fft_scalar i; +}kiss_fft_cpx; + +typedef struct kiss_fft_state* kiss_fft_cfg; + +/* + * kiss_fft_alloc + * + * Initialize a FFT (or IFFT) algorithm's cfg/state buffer. + * + * typical usage: kiss_fft_cfg mycfg=kiss_fft_alloc(1024,0,NULL,NULL); + * + * The return value from fft_alloc is a cfg buffer used internally + * by the fft routine or NULL. + * + * If lenmem is NULL, then kiss_fft_alloc will allocate a cfg buffer using malloc. + * The returned value should be free()d when done to avoid memory leaks. + * + * The state can be placed in a user supplied buffer 'mem': + * If lenmem is not NULL and mem is not NULL and *lenmem is large enough, + * then the function places the cfg in mem and the size used in *lenmem + * and returns mem. + * + * If lenmem is not NULL and ( mem is NULL or *lenmem is not large enough), + * then the function returns NULL and places the minimum cfg + * buffer size in *lenmem. + * */ + +kiss_fft_cfg KISS_FFT_API kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem); + +/* + * kiss_fft(cfg,in_out_buf) + * + * Perform an FFT on a complex input buffer. + * for a forward FFT, + * fin should be f[0] , f[1] , ... ,f[nfft-1] + * fout will be F[0] , F[1] , ... ,F[nfft-1] + * Note that each element is complex and can be accessed like + f[k].r and f[k].i + * */ +void KISS_FFT_API kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout); + +/* + A more generic version of the above function. It reads its input from every Nth sample. + * */ +void KISS_FFT_API kiss_fft_stride(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int fin_stride); + +/* If kiss_fft_alloc allocated a buffer, it is one contiguous + buffer and can be simply free()d when no longer needed*/ +#define kiss_fft_free KISS_FFT_FREE + +/* + Cleans up some memory that gets managed internally. Not necessary to call, but it might clean up + your compiler output to call this before you exit. +*/ +void KISS_FFT_API kiss_fft_cleanup(void); + + +/* + * Returns the smallest integer k, such that k>=n and k has only "fast" factors (2,3,5) + */ +int KISS_FFT_API kiss_fft_next_fast_size(int n); + +/* for real ffts, we need an even size */ +#define kiss_fftr_next_fast_size_real(n) \ + (kiss_fft_next_fast_size( ((n)+1)>>1)<<1) + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/kiss/kiss_fft_log.h b/kiss/kiss_fft_log.h new file mode 100644 index 0000000..b5b631a --- /dev/null +++ b/kiss/kiss_fft_log.h @@ -0,0 +1,36 @@ +/* + * Copyright (c) 2003-2010, Mark Borgerding. All rights reserved. + * This file is part of KISS FFT - https://github.com/mborgerding/kissfft + * + * SPDX-License-Identifier: BSD-3-Clause + * See COPYING file for more information. + */ + +#ifndef kiss_fft_log_h +#define kiss_fft_log_h + +#define ERROR 1 +#define WARNING 2 +#define INFO 3 +#define DEBUG 4 + +#define STRINGIFY(x) #x +#define TOSTRING(x) STRINGIFY(x) + +#if defined(NDEBUG) +# define KISS_FFT_LOG_MSG(severity, ...) ((void)0) +#else +# define KISS_FFT_LOG_MSG(severity, ...) \ + fprintf(stderr, "[" #severity "] " __FILE__ ":" TOSTRING(__LINE__) " "); \ + fprintf(stderr, __VA_ARGS__); \ + fprintf(stderr, "\n") +#endif + +#define KISS_FFT_ERROR(...) KISS_FFT_LOG_MSG(ERROR, __VA_ARGS__) +#define KISS_FFT_WARNING(...) KISS_FFT_LOG_MSG(WARNING, __VA_ARGS__) +#define KISS_FFT_INFO(...) KISS_FFT_LOG_MSG(INFO, __VA_ARGS__) +#define KISS_FFT_DEBUG(...) KISS_FFT_LOG_MSG(DEBUG, __VA_ARGS__) + + + +#endif /* kiss_fft_log_h */
\ No newline at end of file diff --git a/kiss/kiss_fftnd.c b/kiss/kiss_fftnd.c new file mode 100644 index 0000000..5d5b089 --- /dev/null +++ b/kiss/kiss_fftnd.c @@ -0,0 +1,188 @@ +/* + * Copyright (c) 2003-2004, Mark Borgerding. All rights reserved. + * This file is part of KISS FFT - https://github.com/mborgerding/kissfft + * + * SPDX-License-Identifier: BSD-3-Clause + * See COPYING file for more information. + */ + +#include "kiss_fftnd.h" +#include "_kiss_fft_guts.h" + +struct kiss_fftnd_state{ + int dimprod; /* dimsum would be mighty tasty right now */ + int ndims; + int *dims; + kiss_fft_cfg *states; /* cfg states for each dimension */ + kiss_fft_cpx * tmpbuf; /*buffer capable of hold the entire input */ +}; + +kiss_fftnd_cfg kiss_fftnd_alloc(const int *dims,int ndims,int inverse_fft,void*mem,size_t*lenmem) +{ + KISS_FFT_ALIGN_CHECK(mem) + + kiss_fftnd_cfg st = NULL; + int i; + int dimprod=1; + size_t memneeded = KISS_FFT_ALIGN_SIZE_UP(sizeof(struct kiss_fftnd_state)); + char * ptr = NULL; + + for (i=0;i<ndims;++i) { + size_t sublen=0; + kiss_fft_alloc (dims[i], inverse_fft, NULL, &sublen); + memneeded += sublen; /* st->states[i] */ + dimprod *= dims[i]; + } + memneeded += KISS_FFT_ALIGN_SIZE_UP(sizeof(int) * ndims);/* st->dims */ + memneeded += KISS_FFT_ALIGN_SIZE_UP(sizeof(void*) * ndims);/* st->states */ + memneeded += KISS_FFT_ALIGN_SIZE_UP(sizeof(kiss_fft_cpx) * dimprod); /* st->tmpbuf */ + + if (lenmem == NULL) {/* allocate for the caller*/ + ptr = (char *) malloc (memneeded); + } else { /* initialize supplied buffer if big enough */ + if (*lenmem >= memneeded) + ptr = (char *) mem; + *lenmem = memneeded; /*tell caller how big struct is (or would be) */ + } + if (!ptr) + return NULL; /*malloc failed or buffer too small */ + + st = (kiss_fftnd_cfg) ptr; + st->dimprod = dimprod; + st->ndims = ndims; + ptr += KISS_FFT_ALIGN_SIZE_UP(sizeof(struct kiss_fftnd_state)); + + st->states = (kiss_fft_cfg *)ptr; + ptr += KISS_FFT_ALIGN_SIZE_UP(sizeof(void*) * ndims); + + st->dims = (int*)ptr; + ptr += KISS_FFT_ALIGN_SIZE_UP(sizeof(int) * ndims); + + st->tmpbuf = (kiss_fft_cpx*)ptr; + ptr += KISS_FFT_ALIGN_SIZE_UP(sizeof(kiss_fft_cpx) * dimprod); + + for (i=0;i<ndims;++i) { + size_t len; + st->dims[i] = dims[i]; + kiss_fft_alloc (st->dims[i], inverse_fft, NULL, &len); + st->states[i] = kiss_fft_alloc (st->dims[i], inverse_fft, ptr,&len); + ptr += len; + } + /* +Hi there! + +If you're looking at this particular code, it probably means you've got a brain-dead bounds checker +that thinks the above code overwrites the end of the array. + +It doesn't. + +-- Mark + +P.S. +The below code might give you some warm fuzzies and help convince you. + */ + if ( ptr - (char*)st != (int)memneeded ) { + fprintf(stderr, + "################################################################################\n" + "Internal error! Memory allocation miscalculation\n" + "################################################################################\n" + ); + } + return st; +} + +/* + This works by tackling one dimension at a time. + + In effect, + Each stage starts out by reshaping the matrix into a DixSi 2d matrix. + A Di-sized fft is taken of each column, transposing the matrix as it goes. + +Here's a 3-d example: +Take a 2x3x4 matrix, laid out in memory as a contiguous buffer + [ [ [ a b c d ] [ e f g h ] [ i j k l ] ] + [ [ m n o p ] [ q r s t ] [ u v w x ] ] ] + +Stage 0 ( D=2): treat the buffer as a 2x12 matrix + [ [a b ... k l] + [m n ... w x] ] + + FFT each column with size 2. + Transpose the matrix at the same time using kiss_fft_stride. + + [ [ a+m a-m ] + [ b+n b-n] + ... + [ k+w k-w ] + [ l+x l-x ] ] + + Note fft([x y]) == [x+y x-y] + +Stage 1 ( D=3) treats the buffer (the output of stage D=2) as an 3x8 matrix, + [ [ a+m a-m b+n b-n c+o c-o d+p d-p ] + [ e+q e-q f+r f-r g+s g-s h+t h-t ] + [ i+u i-u j+v j-v k+w k-w l+x l-x ] ] + + And perform FFTs (size=3) on each of the columns as above, transposing + the matrix as it goes. The output of stage 1 is + (Legend: ap = [ a+m e+q i+u ] + am = [ a-m e-q i-u ] ) + + [ [ sum(ap) fft(ap)[0] fft(ap)[1] ] + [ sum(am) fft(am)[0] fft(am)[1] ] + [ sum(bp) fft(bp)[0] fft(bp)[1] ] + [ sum(bm) fft(bm)[0] fft(bm)[1] ] + [ sum(cp) fft(cp)[0] fft(cp)[1] ] + [ sum(cm) fft(cm)[0] fft(cm)[1] ] + [ sum(dp) fft(dp)[0] fft(dp)[1] ] + [ sum(dm) fft(dm)[0] fft(dm)[1] ] ] + +Stage 2 ( D=4) treats this buffer as a 4*6 matrix, + [ [ sum(ap) fft(ap)[0] fft(ap)[1] sum(am) fft(am)[0] fft(am)[1] ] + [ sum(bp) fft(bp)[0] fft(bp)[1] sum(bm) fft(bm)[0] fft(bm)[1] ] + [ sum(cp) fft(cp)[0] fft(cp)[1] sum(cm) fft(cm)[0] fft(cm)[1] ] + [ sum(dp) fft(dp)[0] fft(dp)[1] sum(dm) fft(dm)[0] fft(dm)[1] ] ] + + Then FFTs each column, transposing as it goes. + + The resulting matrix is the 3d FFT of the 2x3x4 input matrix. + + Note as a sanity check that the first element of the final + stage's output (DC term) is + sum( [ sum(ap) sum(bp) sum(cp) sum(dp) ] ) + , i.e. the summation of all 24 input elements. + +*/ +void kiss_fftnd(kiss_fftnd_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout) +{ + int i,k; + const kiss_fft_cpx * bufin=fin; + kiss_fft_cpx * bufout; + + /*arrange it so the last bufout == fout*/ + if ( st->ndims & 1 ) { + bufout = fout; + if (fin==fout) { + memcpy( st->tmpbuf, fin, sizeof(kiss_fft_cpx) * st->dimprod ); + bufin = st->tmpbuf; + } + }else + bufout = st->tmpbuf; + + for ( k=0; k < st->ndims; ++k) { + int curdim = st->dims[k]; + int stride = st->dimprod / curdim; + + for ( i=0 ; i<stride ; ++i ) + kiss_fft_stride( st->states[k], bufin+i , bufout+i*curdim, stride ); + + /*toggle back and forth between the two buffers*/ + if (bufout == st->tmpbuf){ + bufout = fout; + bufin = st->tmpbuf; + }else{ + bufout = st->tmpbuf; + bufin = fout; + } + } +} diff --git a/kiss/kiss_fftnd.h b/kiss/kiss_fftnd.h new file mode 100644 index 0000000..956ba94 --- /dev/null +++ b/kiss/kiss_fftnd.h @@ -0,0 +1,26 @@ +/* + * Copyright (c) 2003-2004, Mark Borgerding. All rights reserved. + * This file is part of KISS FFT - https://github.com/mborgerding/kissfft + * + * SPDX-License-Identifier: BSD-3-Clause + * See COPYING file for more information. + */ + +#ifndef KISS_FFTND_H +#define KISS_FFTND_H + +#include "kiss_fft.h" + +#ifdef __cplusplus +extern "C" { +#endif + +typedef struct kiss_fftnd_state * kiss_fftnd_cfg; + +kiss_fftnd_cfg KISS_FFT_API kiss_fftnd_alloc(const int *dims,int ndims,int inverse_fft,void*mem,size_t*lenmem); +void KISS_FFT_API kiss_fftnd(kiss_fftnd_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout); + +#ifdef __cplusplus +} +#endif +#endif diff --git a/kiss/kiss_fftndr.c b/kiss/kiss_fftndr.c new file mode 100644 index 0000000..e979d03 --- /dev/null +++ b/kiss/kiss_fftndr.c @@ -0,0 +1,120 @@ +/* + * Copyright (c) 2003-2004, Mark Borgerding. All rights reserved. + * This file is part of KISS FFT - https://github.com/mborgerding/kissfft + * + * SPDX-License-Identifier: BSD-3-Clause + * See COPYING file for more information. + */ + +#include "kiss_fftndr.h" +#include "_kiss_fft_guts.h" +#define MAX(x,y) ( ( (x)<(y) )?(y):(x) ) + +struct kiss_fftndr_state +{ + int dimReal; + int dimOther; + kiss_fftr_cfg cfg_r; + kiss_fftnd_cfg cfg_nd; + void * tmpbuf; +}; + +static int prod(const int *dims, int ndims) +{ + int x=1; + while (ndims--) + x *= *dims++; + return x; +} + +kiss_fftndr_cfg kiss_fftndr_alloc(const int *dims,int ndims,int inverse_fft,void*mem,size_t*lenmem) +{ + KISS_FFT_ALIGN_CHECK(mem) + + kiss_fftndr_cfg st = NULL; + size_t nr=0 , nd=0,ntmp=0; + int dimReal = dims[ndims-1]; + int dimOther = prod(dims,ndims-1); + size_t memneeded; + char * ptr = NULL; + + (void)kiss_fftr_alloc(dimReal,inverse_fft,NULL,&nr); + (void)kiss_fftnd_alloc(dims,ndims-1,inverse_fft,NULL,&nd); + ntmp = + MAX( 2*dimOther , dimReal+2) * sizeof(kiss_fft_scalar) // freq buffer for one pass + + dimOther*(dimReal+2) * sizeof(kiss_fft_scalar); // large enough to hold entire input in case of in-place + + memneeded = KISS_FFT_ALIGN_SIZE_UP(sizeof( struct kiss_fftndr_state )) + KISS_FFT_ALIGN_SIZE_UP(nr) + KISS_FFT_ALIGN_SIZE_UP(nd) + KISS_FFT_ALIGN_SIZE_UP(ntmp); + + if (lenmem==NULL) { + ptr = (char*) malloc(memneeded); + }else{ + if (*lenmem >= memneeded) + ptr = (char *)mem; + *lenmem = memneeded; + } + if (ptr==NULL) + return NULL; + + st = (kiss_fftndr_cfg) ptr; + memset( st , 0 , memneeded); + ptr += KISS_FFT_ALIGN_SIZE_UP(sizeof(struct kiss_fftndr_state)); + + st->dimReal = dimReal; + st->dimOther = dimOther; + st->cfg_r = kiss_fftr_alloc( dimReal,inverse_fft,ptr,&nr); + ptr += KISS_FFT_ALIGN_SIZE_UP(nr); + st->cfg_nd = kiss_fftnd_alloc(dims,ndims-1,inverse_fft, ptr,&nd); + ptr += KISS_FFT_ALIGN_SIZE_UP(nd); + st->tmpbuf = ptr; + + return st; +} + +void kiss_fftndr(kiss_fftndr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata) +{ + int k1,k2; + int dimReal = st->dimReal; + int dimOther = st->dimOther; + int nrbins = dimReal/2+1; + + kiss_fft_cpx * tmp1 = (kiss_fft_cpx*)st->tmpbuf; + kiss_fft_cpx * tmp2 = tmp1 + MAX(nrbins,dimOther); + + // timedata is N0 x N1 x ... x Nk real + + // take a real chunk of data, fft it and place the output at correct intervals + for (k1=0;k1<dimOther;++k1) { + kiss_fftr( st->cfg_r, timedata + k1*dimReal , tmp1 ); // tmp1 now holds nrbins complex points + for (k2=0;k2<nrbins;++k2) + tmp2[ k2*dimOther+k1 ] = tmp1[k2]; + } + + for (k2=0;k2<nrbins;++k2) { + kiss_fftnd(st->cfg_nd, tmp2+k2*dimOther, tmp1); // tmp1 now holds dimOther complex points + for (k1=0;k1<dimOther;++k1) + freqdata[ k1*(nrbins) + k2] = tmp1[k1]; + } +} + +void kiss_fftndri(kiss_fftndr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata) +{ + int k1,k2; + int dimReal = st->dimReal; + int dimOther = st->dimOther; + int nrbins = dimReal/2+1; + kiss_fft_cpx * tmp1 = (kiss_fft_cpx*)st->tmpbuf; + kiss_fft_cpx * tmp2 = tmp1 + MAX(nrbins,dimOther); + + for (k2=0;k2<nrbins;++k2) { + for (k1=0;k1<dimOther;++k1) + tmp1[k1] = freqdata[ k1*(nrbins) + k2 ]; + kiss_fftnd(st->cfg_nd, tmp1, tmp2+k2*dimOther); + } + + for (k1=0;k1<dimOther;++k1) { + for (k2=0;k2<nrbins;++k2) + tmp1[k2] = tmp2[ k2*dimOther+k1 ]; + kiss_fftri( st->cfg_r,tmp1,timedata + k1*dimReal); + } +} diff --git a/kiss/kiss_fftndr.h b/kiss/kiss_fftndr.h new file mode 100644 index 0000000..0d56a1f --- /dev/null +++ b/kiss/kiss_fftndr.h @@ -0,0 +1,55 @@ +/* + * Copyright (c) 2003-2004, Mark Borgerding. All rights reserved. + * This file is part of KISS FFT - https://github.com/mborgerding/kissfft + * + * SPDX-License-Identifier: BSD-3-Clause + * See COPYING file for more information. + */ + +#ifndef KISS_NDR_H +#define KISS_NDR_H + +#include "kiss_fft.h" +#include "kiss_fftr.h" +#include "kiss_fftnd.h" + +#ifdef __cplusplus +extern "C" { +#endif + +typedef struct kiss_fftndr_state *kiss_fftndr_cfg; + + +kiss_fftndr_cfg KISS_FFT_API kiss_fftndr_alloc(const int *dims,int ndims,int inverse_fft,void*mem,size_t*lenmem); +/* + dims[0] must be even + + If you don't care to allocate space, use mem = lenmem = NULL +*/ + + +void KISS_FFT_API kiss_fftndr( + kiss_fftndr_cfg cfg, + const kiss_fft_scalar *timedata, + kiss_fft_cpx *freqdata); +/* + input timedata has dims[0] X dims[1] X ... X dims[ndims-1] scalar points + output freqdata has dims[0] X dims[1] X ... X dims[ndims-1]/2+1 complex points +*/ + +void KISS_FFT_API kiss_fftndri( + kiss_fftndr_cfg cfg, + const kiss_fft_cpx *freqdata, + kiss_fft_scalar *timedata); +/* + input and output dimensions are the exact opposite of kiss_fftndr +*/ + + +#define kiss_fftndr_free free + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/kiss/kiss_fftr.c b/kiss/kiss_fftr.c new file mode 100644 index 0000000..778a9a6 --- /dev/null +++ b/kiss/kiss_fftr.c @@ -0,0 +1,155 @@ +/* + * Copyright (c) 2003-2004, Mark Borgerding. All rights reserved. + * This file is part of KISS FFT - https://github.com/mborgerding/kissfft + * + * SPDX-License-Identifier: BSD-3-Clause + * See COPYING file for more information. + */ + +#include "kiss_fftr.h" +#include "_kiss_fft_guts.h" + +struct kiss_fftr_state{ + kiss_fft_cfg substate; + kiss_fft_cpx * tmpbuf; + kiss_fft_cpx * super_twiddles; +#ifdef USE_SIMD + void * pad; +#endif +}; + +kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem) +{ + KISS_FFT_ALIGN_CHECK(mem) + + int i; + kiss_fftr_cfg st = NULL; + size_t subsize = 0, memneeded; + + if (nfft & 1) { + KISS_FFT_ERROR("Real FFT optimization must be even."); + return NULL; + } + nfft >>= 1; + + kiss_fft_alloc (nfft, inverse_fft, NULL, &subsize); + memneeded = sizeof(struct kiss_fftr_state) + subsize + sizeof(kiss_fft_cpx) * ( nfft * 3 / 2); + + if (lenmem == NULL) { + st = (kiss_fftr_cfg) KISS_FFT_MALLOC (memneeded); + } else { + if (*lenmem >= memneeded) + st = (kiss_fftr_cfg) mem; + *lenmem = memneeded; + } + if (!st) + return NULL; + + st->substate = (kiss_fft_cfg) (st + 1); /*just beyond kiss_fftr_state struct */ + st->tmpbuf = (kiss_fft_cpx *) (((char *) st->substate) + subsize); + st->super_twiddles = st->tmpbuf + nfft; + kiss_fft_alloc(nfft, inverse_fft, st->substate, &subsize); + + for (i = 0; i < nfft/2; ++i) { + double phase = + -3.14159265358979323846264338327 * ((double) (i+1) / nfft + .5); + if (inverse_fft) + phase *= -1; + kf_cexp (st->super_twiddles+i,phase); + } + return st; +} + +void kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata) +{ + /* input buffer timedata is stored row-wise */ + int k,ncfft; + kiss_fft_cpx fpnk,fpk,f1k,f2k,tw,tdc; + + if ( st->substate->inverse) { + KISS_FFT_ERROR("kiss fft usage error: improper alloc"); + return;/* The caller did not call the correct function */ + } + + ncfft = st->substate->nfft; + + /*perform the parallel fft of two real signals packed in real,imag*/ + kiss_fft( st->substate , (const kiss_fft_cpx*)timedata, st->tmpbuf ); + /* The real part of the DC element of the frequency spectrum in st->tmpbuf + * contains the sum of the even-numbered elements of the input time sequence + * The imag part is the sum of the odd-numbered elements + * + * The sum of tdc.r and tdc.i is the sum of the input time sequence. + * yielding DC of input time sequence + * The difference of tdc.r - tdc.i is the sum of the input (dot product) [1,-1,1,-1... + * yielding Nyquist bin of input time sequence + */ + + tdc.r = st->tmpbuf[0].r; + tdc.i = st->tmpbuf[0].i; + C_FIXDIV(tdc,2); + CHECK_OVERFLOW_OP(tdc.r ,+, tdc.i); + CHECK_OVERFLOW_OP(tdc.r ,-, tdc.i); + freqdata[0].r = tdc.r + tdc.i; + freqdata[ncfft].r = tdc.r - tdc.i; +#ifdef USE_SIMD + freqdata[ncfft].i = freqdata[0].i = _mm_set1_ps(0); +#else + freqdata[ncfft].i = freqdata[0].i = 0; +#endif + + for ( k=1;k <= ncfft/2 ; ++k ) { + fpk = st->tmpbuf[k]; + fpnk.r = st->tmpbuf[ncfft-k].r; + fpnk.i = - st->tmpbuf[ncfft-k].i; + C_FIXDIV(fpk,2); + C_FIXDIV(fpnk,2); + + C_ADD( f1k, fpk , fpnk ); + C_SUB( f2k, fpk , fpnk ); + C_MUL( tw , f2k , st->super_twiddles[k-1]); + + freqdata[k].r = HALF_OF(f1k.r + tw.r); + freqdata[k].i = HALF_OF(f1k.i + tw.i); + freqdata[ncfft-k].r = HALF_OF(f1k.r - tw.r); + freqdata[ncfft-k].i = HALF_OF(tw.i - f1k.i); + } +} + +void kiss_fftri(kiss_fftr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata) +{ + /* input buffer timedata is stored row-wise */ + int k, ncfft; + + if (st->substate->inverse == 0) { + KISS_FFT_ERROR("kiss fft usage error: improper alloc"); + return;/* The caller did not call the correct function */ + } + + ncfft = st->substate->nfft; + + st->tmpbuf[0].r = freqdata[0].r + freqdata[ncfft].r; + st->tmpbuf[0].i = freqdata[0].r - freqdata[ncfft].r; + C_FIXDIV(st->tmpbuf[0],2); + + for (k = 1; k <= ncfft / 2; ++k) { + kiss_fft_cpx fk, fnkc, fek, fok, tmp; + fk = freqdata[k]; + fnkc.r = freqdata[ncfft - k].r; + fnkc.i = -freqdata[ncfft - k].i; + C_FIXDIV( fk , 2 ); + C_FIXDIV( fnkc , 2 ); + + C_ADD (fek, fk, fnkc); + C_SUB (tmp, fk, fnkc); + C_MUL (fok, tmp, st->super_twiddles[k-1]); + C_ADD (st->tmpbuf[k], fek, fok); + C_SUB (st->tmpbuf[ncfft - k], fek, fok); +#ifdef USE_SIMD + st->tmpbuf[ncfft - k].i *= _mm_set1_ps(-1.0); +#else + st->tmpbuf[ncfft - k].i *= -1; +#endif + } + kiss_fft (st->substate, st->tmpbuf, (kiss_fft_cpx *) timedata); +} diff --git a/kiss/kiss_fftr.h b/kiss/kiss_fftr.h new file mode 100644 index 0000000..7fd73d2 --- /dev/null +++ b/kiss/kiss_fftr.h @@ -0,0 +1,54 @@ +/* + * Copyright (c) 2003-2004, Mark Borgerding. All rights reserved. + * This file is part of KISS FFT - https://github.com/mborgerding/kissfft + * + * SPDX-License-Identifier: BSD-3-Clause + * See COPYING file for more information. + */ + +#ifndef KISS_FTR_H +#define KISS_FTR_H + +#include "kiss_fft.h" +#ifdef __cplusplus +extern "C" { +#endif + + +/* + + Real optimized version can save about 45% cpu time vs. complex fft of a real seq. + + + + */ + +typedef struct kiss_fftr_state *kiss_fftr_cfg; + + +kiss_fftr_cfg KISS_FFT_API kiss_fftr_alloc(int nfft,int inverse_fft,void * mem, size_t * lenmem); +/* + nfft must be even + + If you don't care to allocate space, use mem = lenmem = NULL +*/ + + +void KISS_FFT_API kiss_fftr(kiss_fftr_cfg cfg,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata); +/* + input timedata has nfft scalar points + output freqdata has nfft/2+1 complex points +*/ + +void KISS_FFT_API kiss_fftri(kiss_fftr_cfg cfg,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata); +/* + input freqdata has nfft/2+1 complex points + output timedata has nfft scalar points +*/ + +#define kiss_fftr_free KISS_FFT_FREE + +#ifdef __cplusplus +} +#endif +#endif diff --git a/m4/ax_cxx_compile_stdcxx.m4 b/m4/ax_cxx_compile_stdcxx.m4 index 43087b2..8edf515 100644 --- a/m4/ax_cxx_compile_stdcxx.m4 +++ b/m4/ax_cxx_compile_stdcxx.m4 @@ -10,13 +10,13 @@ # # Check for baseline language coverage in the compiler for the specified # version of the C++ standard. If necessary, add switches to CXX and -# CXXCPP to enable support. VERSION may be '11' (for the C++11 standard) -# or '14' (for the C++14 standard). +# CXXCPP to enable support. VERSION may be '11', '14', '17', or '20' for +# the respective C++ standard version. # # The second argument, if specified, indicates whether you insist on an # extended mode (e.g. -std=gnu++11) or a strict conformance mode (e.g. # -std=c++11). If neither is specified, you get whatever works, with -# preference for an extended mode. +# preference for no added switch, and then for an extended mode. # # The third argument, if specified 'mandatory' or if left unspecified, # indicates that baseline support for the specified C++ standard is @@ -35,13 +35,15 @@ # Copyright (c) 2015 Moritz Klammler <moritz@klammler.eu> # Copyright (c) 2016, 2018 Krzesimir Nowak <qdlacz@gmail.com> # Copyright (c) 2019 Enji Cooper <yaneurabeya@gmail.com> +# Copyright (c) 2020 Jason Merrill <jason@redhat.com> +# Copyright (c) 2021 Jörn Heusipp <osmanx@problemloesungsmaschine.de> # # Copying and distribution of this file, with or without modification, are # permitted in any medium without royalty provided the copyright notice # and this notice are preserved. This file is offered as-is, without any # warranty. -#serial 11 +#serial 18 dnl This macro is based on the code from the AX_CXX_COMPILE_STDCXX_11 macro dnl (serial version number 13). @@ -50,6 +52,7 @@ AC_DEFUN([AX_CXX_COMPILE_STDCXX], [dnl m4_if([$1], [11], [ax_cxx_compile_alternatives="11 0x"], [$1], [14], [ax_cxx_compile_alternatives="14 1y"], [$1], [17], [ax_cxx_compile_alternatives="17 1z"], + [$1], [20], [ax_cxx_compile_alternatives="20"], [m4_fatal([invalid first argument `$1' to AX_CXX_COMPILE_STDCXX])])dnl m4_if([$2], [], [], [$2], [ext], [], @@ -62,6 +65,16 @@ AC_DEFUN([AX_CXX_COMPILE_STDCXX], [dnl AC_LANG_PUSH([C++])dnl ac_success=no + m4_if([$2], [], [dnl + AC_CACHE_CHECK(whether $CXX supports C++$1 features by default, + ax_cv_cxx_compile_cxx$1, + [AC_COMPILE_IFELSE([AC_LANG_SOURCE([_AX_CXX_COMPILE_STDCXX_testbody_$1])], + [ax_cv_cxx_compile_cxx$1=yes], + [ax_cv_cxx_compile_cxx$1=no])]) + if test x$ax_cv_cxx_compile_cxx$1 = xyes; then + ac_success=yes + fi]) + m4_if([$2], [noext], [], [dnl if test x$ac_success = xno; then for alternative in ${ax_cxx_compile_alternatives}; do @@ -91,9 +104,18 @@ AC_DEFUN([AX_CXX_COMPILE_STDCXX], [dnl dnl HP's aCC needs +std=c++11 according to: dnl http://h21007.www2.hp.com/portal/download/files/unprot/aCxx/PDF_Release_Notes/769149-001.pdf dnl Cray's crayCC needs "-h std=c++11" + dnl MSVC needs -std:c++NN for C++17 and later (default is C++14) for alternative in ${ax_cxx_compile_alternatives}; do - for switch in -std=c++${alternative} +std=c++${alternative} "-h std=c++${alternative}"; do - cachevar=AS_TR_SH([ax_cv_cxx_compile_cxx$1_$switch]) + for switch in -std=c++${alternative} +std=c++${alternative} "-h std=c++${alternative}" MSVC; do + if test x"$switch" = xMSVC; then + dnl AS_TR_SH maps both `:` and `=` to `_` so -std:c++17 would collide + dnl with -std=c++17. We suffix the cache variable name with _MSVC to + dnl avoid this. + switch=-std:c++${alternative} + cachevar=AS_TR_SH([ax_cv_cxx_compile_cxx$1_${switch}_MSVC]) + else + cachevar=AS_TR_SH([ax_cv_cxx_compile_cxx$1_$switch]) + fi AC_CACHE_CHECK(whether $CXX supports C++$1 features with $switch, $cachevar, [ac_save_CXX="$CXX" @@ -140,7 +162,6 @@ m4_define([_AX_CXX_COMPILE_STDCXX_testbody_11], _AX_CXX_COMPILE_STDCXX_testbody_new_in_11 ) - dnl Test body for checking C++14 support m4_define([_AX_CXX_COMPILE_STDCXX_testbody_14], @@ -148,12 +169,24 @@ m4_define([_AX_CXX_COMPILE_STDCXX_testbody_14], _AX_CXX_COMPILE_STDCXX_testbody_new_in_14 ) +dnl Test body for checking C++17 support + m4_define([_AX_CXX_COMPILE_STDCXX_testbody_17], _AX_CXX_COMPILE_STDCXX_testbody_new_in_11 _AX_CXX_COMPILE_STDCXX_testbody_new_in_14 _AX_CXX_COMPILE_STDCXX_testbody_new_in_17 ) +dnl Test body for checking C++20 support + +m4_define([_AX_CXX_COMPILE_STDCXX_testbody_20], + _AX_CXX_COMPILE_STDCXX_testbody_new_in_11 + _AX_CXX_COMPILE_STDCXX_testbody_new_in_14 + _AX_CXX_COMPILE_STDCXX_testbody_new_in_17 + _AX_CXX_COMPILE_STDCXX_testbody_new_in_20 +) + + dnl Tests for new features in C++11 m4_define([_AX_CXX_COMPILE_STDCXX_testbody_new_in_11], [[ @@ -165,7 +198,11 @@ m4_define([_AX_CXX_COMPILE_STDCXX_testbody_new_in_11], [[ #error "This is not a C++ compiler" -#elif __cplusplus < 201103L +// MSVC always sets __cplusplus to 199711L in older versions; newer versions +// only set it correctly if /Zc:__cplusplus is specified as well as a +// /std:c++NN switch: +// https://devblogs.microsoft.com/cppblog/msvc-now-correctly-reports-__cplusplus/ +#elif __cplusplus < 201103L && !defined _MSC_VER #error "This is not a C++11 compiler" @@ -456,7 +493,7 @@ m4_define([_AX_CXX_COMPILE_STDCXX_testbody_new_in_14], [[ #error "This is not a C++ compiler" -#elif __cplusplus < 201402L +#elif __cplusplus < 201402L && !defined _MSC_VER #error "This is not a C++14 compiler" @@ -580,7 +617,7 @@ m4_define([_AX_CXX_COMPILE_STDCXX_testbody_new_in_17], [[ #error "This is not a C++ compiler" -#elif __cplusplus < 201703L +#elif __cplusplus < 201703L && !defined _MSC_VER #error "This is not a C++17 compiler" @@ -946,6 +983,36 @@ namespace cxx17 } // namespace cxx17 -#endif // __cplusplus < 201703L +#endif // __cplusplus < 201703L && !defined _MSC_VER + +]]) + + +dnl Tests for new features in C++20 + +m4_define([_AX_CXX_COMPILE_STDCXX_testbody_new_in_20], [[ + +#ifndef __cplusplus + +#error "This is not a C++ compiler" + +#elif __cplusplus < 202002L && !defined _MSC_VER + +#error "This is not a C++20 compiler" + +#else + +#include <version> + +namespace cxx20 +{ + +// As C++20 supports feature test macros in the standard, there is no +// immediate need to actually test for feature availability on the +// Autoconf side. + +} // namespace cxx20 + +#endif // __cplusplus < 202002L && !defined _MSC_VER ]]) diff --git a/src/Buffer.h b/src/Buffer.h index 077d654..711b804 100644 --- a/src/Buffer.h +++ b/src/Buffer.h @@ -34,7 +34,10 @@ #include <vector> #include <memory> #include <complex> +#include "fpm/fixed.hpp" + typedef std::complex<float> complexf; +typedef std::complex<fpm::fixed_16_16> complexfix; /* Buffer is a container for a byte array, which is memory-aligned * to 32 bytes for SSE performance. diff --git a/src/ConfigParser.cpp b/src/ConfigParser.cpp index 4a1e362..a48f7e7 100644 --- a/src/ConfigParser.cpp +++ b/src/ConfigParser.cpp @@ -154,6 +154,8 @@ static void parse_configfile( mod_settings.showProcessTime); // modulator parameters: + mod_settings.fixedPoint = pt.GetInteger("modulator.fixed_point", mod_settings.fixedPoint); + const string gainMode_setting = pt.Get("modulator.gainmode", "var"); mod_settings.gainMode = parse_gainmode(gainMode_setting); mod_settings.gainmodeVariance = pt.GetReal("modulator.normalise_variance", diff --git a/src/ConfigParser.h b/src/ConfigParser.h index ae76dee..f3a2af9 100644 --- a/src/ConfigParser.h +++ b/src/ConfigParser.h @@ -51,6 +51,8 @@ struct mod_settings_t { bool useLimeOutput = false; bool useBladeRFOutput = false; + bool fixedPoint = false; + size_t outputRate = 2048000; size_t clockRate = 0; unsigned dabMode = 1; diff --git a/src/DabMod.cpp b/src/DabMod.cpp index 4e338f9..5f8412b 100644 --- a/src/DabMod.cpp +++ b/src/DabMod.cpp @@ -249,7 +249,15 @@ static shared_ptr<ModOutput> prepare_output(mod_settings_t& s) { shared_ptr<ModOutput> output; - if (s.useFileOutput) { + if (s.fixedPoint) { + if (s.useFileOutput) { + output = make_shared<OutputFile>(s.outputName, s.fileOutputShowMetadata); + } + else { + throw runtime_error("Fixed point only works with file output"); + } + } + else if (s.useFileOutput) { if (s.fileOutputFormat == "complexf") { output = make_shared<OutputFile>(s.outputName, s.fileOutputShowMetadata); } @@ -413,7 +421,7 @@ int launch_modulator(int argc, char* argv[]) ModulatorData m; rcs.enrol(&m); - { + if (not mod_settings.fixedPoint) { // This is mostly useful on ARM systems where FFTW planning takes some time. If we do it here // it will be done before the modulator starts up etiLog.level(debug) << "Running FFTW planning..."; diff --git a/src/DabModulator.cpp b/src/DabModulator.cpp index 4a29132..d48f1a2 100644 --- a/src/DabModulator.cpp +++ b/src/DabModulator.cpp @@ -142,14 +142,14 @@ int DabModulator::process(Buffer* dataOut) auto cifMux = make_shared<FrameMultiplexer>(m_etiSource); auto cifPart = make_shared<BlockPartitioner>(mode); - auto cifMap = make_shared<QpskSymbolMapper>(m_nbCarriers); - auto cifRef = make_shared<PhaseReference>(mode); - auto cifFreq = make_shared<FrequencyInterleaver>(mode); - auto cifDiff = make_shared<DifferentialModulator>(m_nbCarriers); + auto cifMap = make_shared<QpskSymbolMapper>(m_nbCarriers, m_settings.fixedPoint); + auto cifRef = make_shared<PhaseReference>(mode, m_settings.fixedPoint); + auto cifFreq = make_shared<FrequencyInterleaver>(mode, m_settings.fixedPoint); + auto cifDiff = make_shared<DifferentialModulator>(m_nbCarriers, m_settings.fixedPoint); - auto cifNull = make_shared<NullSymbol>(m_nbCarriers); - auto cifSig = make_shared<SignalMultiplexer>( - (1 + m_nbSymbols) * m_nbCarriers * sizeof(complexf)); + auto cifNull = make_shared<NullSymbol>(m_nbCarriers, + m_settings.fixedPoint ? sizeof(complexfix) : sizeof(complexf)); + auto cifSig = make_shared<SignalMultiplexer>(); // TODO this needs a review bool useCicEq = false; @@ -182,44 +182,66 @@ int DabModulator::process(Buffer* dataOut) m_settings.dabMode, m_settings.tiiConfig); rcs.enrol(tii.get()); - tiiRef = make_shared<PhaseReference>(mode); + tiiRef = make_shared<PhaseReference>(mode, m_settings.fixedPoint); } catch (const TIIError& e) { etiLog.level(error) << "Could not initialise TII: " << e.what(); } - auto cifOfdm = make_shared<OfdmGenerator>( - (1 + m_nbSymbols), - m_nbCarriers, - m_spacing, - m_settings.enableCfr, - m_settings.cfrClip, - m_settings.cfrErrorClip); + shared_ptr<ModPlugin> cifOfdm; + + if (m_settings.fixedPoint) { + cifOfdm = make_shared<OfdmGeneratorFixed>( + (1 + m_nbSymbols), + m_nbCarriers, + m_spacing, + m_settings.enableCfr, + m_settings.cfrClip, + m_settings.cfrErrorClip); + } + else { + auto ofdm = make_shared<OfdmGeneratorCF32>( + (1 + m_nbSymbols), + m_nbCarriers, + m_spacing, + m_settings.enableCfr, + m_settings.cfrClip, + m_settings.cfrErrorClip); + + rcs.enrol(ofdm.get()); + cifOfdm = ofdm; + } - rcs.enrol(cifOfdm.get()); + shared_ptr<GainControl> cifGain; - auto cifGain = make_shared<GainControl>( - m_spacing, - m_settings.gainMode, - m_settings.digitalgain, - m_settings.normalise, - m_settings.gainmodeVariance); + if (not m_settings.fixedPoint) { + cifGain = make_shared<GainControl>( + m_spacing, + m_settings.gainMode, + m_settings.digitalgain, + m_settings.normalise, + m_settings.gainmodeVariance); - rcs.enrol(cifGain.get()); + rcs.enrol(cifGain.get()); + } auto cifGuard = make_shared<GuardIntervalInserter>( m_nbSymbols, m_spacing, m_nullSize, m_symSize, - m_settings.ofdmWindowOverlap); + m_settings.ofdmWindowOverlap, m_settings.fixedPoint); rcs.enrol(cifGuard.get()); shared_ptr<FIRFilter> cifFilter; if (not m_settings.filterTapsFilename.empty()) { + if (m_settings.fixedPoint) throw std::runtime_error("fixed point doesn't support fir filter"); + cifFilter = make_shared<FIRFilter>(m_settings.filterTapsFilename); rcs.enrol(cifFilter.get()); } shared_ptr<MemlessPoly> cifPoly; if (not m_settings.polyCoefFilename.empty()) { + if (m_settings.fixedPoint) throw std::runtime_error("fixed point doesn't support predistortion"); + cifPoly = make_shared<MemlessPoly>(m_settings.polyCoefFilename, m_settings.polyNumThreads); rcs.enrol(cifPoly.get()); @@ -227,6 +249,8 @@ int DabModulator::process(Buffer* dataOut) shared_ptr<Resampler> cifRes; if (m_settings.outputRate != 2048000) { + if (m_settings.fixedPoint) throw std::runtime_error("fixed point doesn't support resampler"); + cifRes = make_shared<Resampler>( 2048000, m_settings.outputRate, @@ -234,6 +258,8 @@ int DabModulator::process(Buffer* dataOut) } if (not m_format.empty()) { + if (m_settings.fixedPoint) throw std::runtime_error("fixed point doesn't support format converter"); + m_formatConverter = make_shared<FormatConverter>(m_format); } diff --git a/src/DifferentialModulator.cpp b/src/DifferentialModulator.cpp index 54e97e2..cfebf65 100644 --- a/src/DifferentialModulator.cpp +++ b/src/DifferentialModulator.cpp @@ -26,9 +26,10 @@ #include <stdexcept> #include <cstring> -DifferentialModulator::DifferentialModulator(size_t carriers) : +DifferentialModulator::DifferentialModulator(size_t carriers, bool fixedPoint) : ModMux(), - d_carriers(carriers) + m_carriers(carriers), + m_fixedPoint(fixedPoint) { PDEBUG("DifferentialModulator::DifferentialModulator(%zu)\n", carriers); @@ -38,10 +39,42 @@ DifferentialModulator::DifferentialModulator(size_t carriers) : DifferentialModulator::~DifferentialModulator() { PDEBUG("DifferentialModulator::~DifferentialModulator()\n"); - } +template<typename T> +void do_process(size_t carriers, std::vector<Buffer*> dataIn, Buffer* dataOut) +{ + size_t phaseSize = dataIn[0]->getLength() / sizeof(T); + size_t dataSize = dataIn[1]->getLength() / sizeof(T); + dataOut->setLength((phaseSize + dataSize) * sizeof(T)); + + const T* phase = reinterpret_cast<const T*>(dataIn[0]->getData()); + const T* in = reinterpret_cast<const T*>(dataIn[1]->getData()); + T* out = reinterpret_cast<T*>(dataOut->getData()); + + if (phaseSize != carriers) { + throw std::runtime_error( + "DifferentialModulator::process input phase size not valid!"); + } + if (dataSize % carriers != 0) { + throw std::runtime_error( + "DifferentialModulator::process input data size not valid!"); + } + + memcpy(dataOut->getData(), phase, phaseSize * sizeof(T)); + for (size_t i = 0; i < dataSize; i += carriers) { + for (size_t j = 0; j < carriers; j += 4) { + out[carriers + j] = out[j] * in[j]; + out[carriers + j + 1] = out[j + 1] * in[j + 1]; + out[carriers + j + 2] = out[j + 2] * in[j + 2]; + out[carriers + j + 3] = out[j + 3] * in[j + 3]; + } + in += carriers; + out += carriers; + } +} + // dataIn[0] -> phase reference // dataIn[1] -> data symbols int DifferentialModulator::process(std::vector<Buffer*> dataIn, Buffer* dataOut) @@ -63,33 +96,11 @@ int DifferentialModulator::process(std::vector<Buffer*> dataIn, Buffer* dataOut) "DifferentialModulator::process nb of input streams not 2!"); } - size_t phaseSize = dataIn[0]->getLength() / sizeof(complexf); - size_t dataSize = dataIn[1]->getLength() / sizeof(complexf); - dataOut->setLength((phaseSize + dataSize) * sizeof(complexf)); - - const complexf* phase = reinterpret_cast<const complexf*>(dataIn[0]->getData()); - const complexf* in = reinterpret_cast<const complexf*>(dataIn[1]->getData()); - complexf* out = reinterpret_cast<complexf*>(dataOut->getData()); - - if (phaseSize != d_carriers) { - throw std::runtime_error( - "DifferentialModulator::process input phase size not valid!"); - } - if (dataSize % d_carriers != 0) { - throw std::runtime_error( - "DifferentialModulator::process input data size not valid!"); + if (m_fixedPoint) { + do_process<complexfix>(m_carriers, dataIn, dataOut); } - - memcpy(dataOut->getData(), phase, phaseSize * sizeof(complexf)); - for (size_t i = 0; i < dataSize; i += d_carriers) { - for (size_t j = 0; j < d_carriers; j += 4) { - out[d_carriers + j] = out[j] * in[j]; - out[d_carriers + j + 1] = out[j + 1] * in[j + 1]; - out[d_carriers + j + 2] = out[j + 2] * in[j + 2]; - out[d_carriers + j + 3] = out[j + 3] * in[j + 3]; - } - in += d_carriers; - out += d_carriers; + else { + do_process<complexf>(m_carriers, dataIn, dataOut); } return dataOut->getLength(); diff --git a/src/DifferentialModulator.h b/src/DifferentialModulator.h index b26ea8b..9cc5081 100644 --- a/src/DifferentialModulator.h +++ b/src/DifferentialModulator.h @@ -35,7 +35,7 @@ class DifferentialModulator : public ModMux { public: - DifferentialModulator(size_t carriers); + DifferentialModulator(size_t carriers, bool fixedPoint); virtual ~DifferentialModulator(); DifferentialModulator(const DifferentialModulator&); DifferentialModulator& operator=(const DifferentialModulator&); @@ -45,6 +45,7 @@ public: const char* name() { return "DifferentialModulator"; } protected: - size_t d_carriers; + size_t m_carriers; + size_t m_fixedPoint; }; diff --git a/src/FrequencyInterleaver.cpp b/src/FrequencyInterleaver.cpp index a62e9f4..856e8d0 100644 --- a/src/FrequencyInterleaver.cpp +++ b/src/FrequencyInterleaver.cpp @@ -28,8 +28,9 @@ #include <cstdlib> -FrequencyInterleaver::FrequencyInterleaver(size_t mode) : - ModCodec() +FrequencyInterleaver::FrequencyInterleaver(size_t mode, bool fixedPoint) : + ModCodec(), + m_fixedPoint(fixedPoint) { PDEBUG("FrequencyInterleaver::FrequencyInterleaver(%zu) @ %p\n", mode, this); @@ -39,45 +40,43 @@ FrequencyInterleaver::FrequencyInterleaver(size_t mode) : size_t beta; switch (mode) { case 1: - d_carriers = 1536; + m_carriers = 1536; num = 2048; beta = 511; break; case 2: - d_carriers = 384; + m_carriers = 384; num = 512; beta = 127; break; case 3: - d_carriers = 192; + m_carriers = 192; num = 256; beta = 63; break; case 0: case 4: - d_carriers = 768; + m_carriers = 768; num = 1024; beta = 255; break; default: PDEBUG("Carriers: %zu\n", (d_carriers >> 1) << 1); - throw std::runtime_error("FrequencyInterleaver::FrequencyInterleaver " - "nb of carriers invalid!"); - break; + throw std::runtime_error("FrequencyInterleaver: invalid dab mode"); } - const int ret = posix_memalign((void**)(&d_indexes), 16, d_carriers * sizeof(size_t)); + const int ret = posix_memalign((void**)(&m_indices), 16, m_carriers * sizeof(size_t)); if (ret != 0) { throw std::runtime_error("memory allocation failed: " + std::to_string(ret)); } - size_t* index = d_indexes; + size_t *index = m_indices; size_t perm = 0; PDEBUG("i: %4u, R: %4u\n", 0, 0); for (size_t j = 1; j < num; ++j) { perm = (alpha * perm + beta) & (num - 1); - if (perm >= ((num - d_carriers) / 2) - && perm <= (num - (num - d_carriers) / 2) + if (perm >= ((num - m_carriers) / 2) + && perm <= (num - (num - m_carriers) / 2) && perm != (num / 2)) { PDEBUG("i: %4zu, R: %4zu, d: %4zu, n: %4zu, k: %5zi, index: %zu\n", j, perm, perm, index - d_indexes, perm - num / 2, @@ -85,8 +84,9 @@ FrequencyInterleaver::FrequencyInterleaver(size_t mode) : ? perm - (1 + (num / 2)) : perm + (d_carriers - (num / 2))); *(index++) = perm > num / 2 ? - perm - (1 + (num / 2)) : perm + (d_carriers - (num / 2)); - } else { + perm - (1 + (num / 2)) : perm + (m_carriers - (num / 2)); + } + else { PDEBUG("i: %4zu, R: %4zu\n", j, perm); } } @@ -97,9 +97,33 @@ FrequencyInterleaver::~FrequencyInterleaver() { PDEBUG("FrequencyInterleaver::~FrequencyInterleaver() @ %p\n", this); - free(d_indexes); + free(m_indices); } +template<typename T> +void do_process(Buffer* const dataIn, Buffer* dataOut, + size_t carriers, const size_t * const indices) +{ + const T* in = reinterpret_cast<const T*>(dataIn->getData()); + T* out = reinterpret_cast<T*>(dataOut->getData()); + size_t sizeIn = dataIn->getLength() / sizeof(T); + + if (sizeIn % carriers != 0) { + throw std::runtime_error( + "FrequencyInterleaver::process input size not valid!"); + } + + for (size_t i = 0; i < sizeIn;) { +// memset(out, 0, d_carriers * sizeof(T)); + for (size_t j = 0; j < carriers; i += 4, j += 4) { + out[indices[j]] = in[i]; + out[indices[j + 1]] = in[i + 1]; + out[indices[j + 2]] = in[i + 2]; + out[indices[j + 3]] = in[i + 3]; + } + out += carriers; + } +} int FrequencyInterleaver::process(Buffer* const dataIn, Buffer* dataOut) { @@ -109,24 +133,11 @@ int FrequencyInterleaver::process(Buffer* const dataIn, Buffer* dataOut) dataOut->setLength(dataIn->getLength()); - const complexf* in = reinterpret_cast<const complexf*>(dataIn->getData()); - complexf* out = reinterpret_cast<complexf*>(dataOut->getData()); - size_t sizeIn = dataIn->getLength() / sizeof(complexf); - - if (sizeIn % d_carriers != 0) { - throw std::runtime_error( - "FrequencyInterleaver::process input size not valid!"); + if (m_fixedPoint) { + do_process<complexfix>(dataIn, dataOut, m_carriers, m_indices); } - - for (size_t i = 0; i < sizeIn;) { -// memset(out, 0, d_carriers * sizeof(complexf)); - for (size_t j = 0; j < d_carriers; i += 4, j += 4) { - out[d_indexes[j]] = in[i]; - out[d_indexes[j + 1]] = in[i + 1]; - out[d_indexes[j + 2]] = in[i + 2]; - out[d_indexes[j + 3]] = in[i + 3]; - } - out += d_carriers; + else { + do_process<complexf>(dataIn, dataOut, m_carriers, m_indices); } return 1; diff --git a/src/FrequencyInterleaver.h b/src/FrequencyInterleaver.h index 43ca21a..b31b968 100644 --- a/src/FrequencyInterleaver.h +++ b/src/FrequencyInterleaver.h @@ -25,16 +25,14 @@ # include <config.h> #endif - #include "ModPlugin.h" #include <sys/types.h> - class FrequencyInterleaver : public ModCodec { public: - FrequencyInterleaver(size_t mode); + FrequencyInterleaver(size_t mode, bool fixedPoint); virtual ~FrequencyInterleaver(); FrequencyInterleaver(const FrequencyInterleaver&) = delete; FrequencyInterleaver& operator=(const FrequencyInterleaver&) = delete; @@ -43,7 +41,8 @@ public: const char* name() override { return "FrequencyInterleaver"; } protected: - size_t d_carriers; - size_t* d_indexes; + bool m_fixedPoint; + size_t m_carriers; + size_t *m_indices; }; diff --git a/src/GuardIntervalInserter.cpp b/src/GuardIntervalInserter.cpp index e200b84..0b39de8 100644 --- a/src/GuardIntervalInserter.cpp +++ b/src/GuardIntervalInserter.cpp @@ -31,35 +31,45 @@ #include <stdexcept> #include <mutex> +GuardIntervalInserter::Params::Params( + size_t nbSymbols, + size_t spacing, + size_t nullSize, + size_t symSize, + size_t& windowOverlap) : + nbSymbols(nbSymbols), + spacing(spacing), + nullSize(nullSize), + symSize(symSize), + windowOverlap(windowOverlap) {} GuardIntervalInserter::GuardIntervalInserter( size_t nbSymbols, size_t spacing, size_t nullSize, size_t symSize, - size_t& windowOverlap) : + size_t& windowOverlap, + bool fixedPoint) : ModCodec(), RemoteControllable("guardinterval"), - d_nbSymbols(nbSymbols), - d_spacing(spacing), - d_nullSize(nullSize), - d_symSize(symSize), - d_windowOverlap(windowOverlap) + m_fixedPoint(fixedPoint), + m_params(nbSymbols, spacing, nullSize, symSize, windowOverlap) { - if (d_nullSize == 0) { + if (nullSize == 0) { throw std::logic_error("NULL symbol must be present"); } + RC_ADD_PARAMETER(windowlen, "Window length for OFDM windowng [0 to disable]"); /* We use a raised-cosine window for the OFDM windowing. - * Each symbol is extended on both sides by d_windowOverlap samples. + * Each symbol is extended on both sides by windowOverlap samples. * * * Sym n |####################| * Sym n+1 |####################| * - * We now extend the symbols by d_windowOverlap (one dash) + * We now extend the symbols by windowOverlap (one dash) * * Sym n extended -|####################|- * Sym n+1 extended -|####################|- @@ -73,7 +83,7 @@ GuardIntervalInserter::GuardIntervalInserter( * / \ * ... ________________/ \__ ... * - * The window length is 2*d_windowOverlap. + * The window length is 2*windowOverlap. */ update_window(windowOverlap); @@ -85,44 +95,45 @@ GuardIntervalInserter::GuardIntervalInserter( void GuardIntervalInserter::update_window(size_t new_window_overlap) { - std::lock_guard<std::mutex> lock(d_windowMutex); + std::lock_guard<std::mutex> lock(m_params.windowMutex); - d_windowOverlap = new_window_overlap; + m_params.windowOverlap = new_window_overlap; - // d_window only contains the rising window edge. - d_window.resize(2*d_windowOverlap); - for (size_t i = 0; i < 2*d_windowOverlap; i++) { - d_window[i] = (float)(0.5 * (1.0 - cos(M_PI * i / (2*d_windowOverlap - 1)))); + // m_params.window only contains the rising window edge. + m_params.window.resize(2*m_params.windowOverlap); + for (size_t i = 0; i < 2*m_params.windowOverlap; i++) { + m_params.window[i] = (float)(0.5 * (1.0 - cos(M_PI * i / (2*m_params.windowOverlap - 1)))); } } -int GuardIntervalInserter::process(Buffer* const dataIn, Buffer* dataOut) +template<typename T> +int do_process(const GuardIntervalInserter::Params& p, Buffer* const dataIn, Buffer* dataOut) { - PDEBUG("GuardIntervalInserter::process(dataIn: %p, dataOut: %p)\n", + PDEBUG("GuardIntervalInserter do_process(dataIn: %p, dataOut: %p)\n", dataIn, dataOut); - std::lock_guard<std::mutex> lock(d_windowMutex); + std::lock_guard<std::mutex> lock(p.windowMutex); - // Every symbol overlaps over a length of d_windowOverlap with + // Every symbol overlaps over a length of windowOverlap with // the previous symbol, and with the next symbol. First symbol // receives no prefix window, because we don't remember the // last symbol from the previous TF (yet). Last symbol also // receives no suffix window, for the same reason. // Overall output buffer length must stay independent of the windowing. - dataOut->setLength((d_nullSize + (d_nbSymbols * d_symSize)) * sizeof(complexf)); + dataOut->setLength((p.nullSize + (p.nbSymbols * p.symSize)) * sizeof(T)); - const complexf* in = reinterpret_cast<const complexf*>(dataIn->getData()); - complexf* out = reinterpret_cast<complexf*>(dataOut->getData()); - size_t sizeIn = dataIn->getLength() / sizeof(complexf); + const T* in = reinterpret_cast<const T*>(dataIn->getData()); + T* out = reinterpret_cast<T*>(dataOut->getData()); + size_t sizeIn = dataIn->getLength() / sizeof(T); - const size_t num_symbols = d_nbSymbols + 1; - if (sizeIn != num_symbols * d_spacing) + const size_t num_symbols = p.nbSymbols + 1; + if (sizeIn != num_symbols * p.spacing) { - PDEBUG("Nb symbols: %zu\n", d_nbSymbols); - PDEBUG("Spacing: %zu\n", d_spacing); - PDEBUG("Null size: %zu\n", d_nullSize); - PDEBUG("Sym size: %zu\n", d_symSize); - PDEBUG("\n%zu != %zu\n", sizeIn, (d_nbSymbols + 1) * d_spacing); + PDEBUG("Nb symbols: %zu\n", p.nbSymbols); + PDEBUG("Spacing: %zu\n", p.spacing); + PDEBUG("Null size: %zu\n", p.nullSize); + PDEBUG("Sym size: %zu\n", p.symSize); + PDEBUG("\n%zu != %zu\n", sizeIn, (p.nbSymbols + 1) * p.spacing); throw std::runtime_error( "GuardIntervalInserter::process input size not valid!"); } @@ -130,65 +141,66 @@ int GuardIntervalInserter::process(Buffer* const dataIn, Buffer* dataOut) // TODO remember the end of the last TF so that we can do some // windowing too. - if (d_windowOverlap) { + + if (p.windowOverlap) { if constexpr (std::is_same_v<complexf, T>) { { // Handle Null symbol separately because it is longer - const size_t prefixlength = d_nullSize - d_spacing; + const size_t prefixlength = p.nullSize - p.spacing; // end = spacing - memcpy(out, &in[d_spacing - prefixlength], - prefixlength * sizeof(complexf)); + memcpy(out, &in[p.spacing - prefixlength], + prefixlength * sizeof(T)); - memcpy(&out[prefixlength], in, (d_spacing - d_windowOverlap) * sizeof(complexf)); + memcpy(&out[prefixlength], in, (p.spacing - p.windowOverlap) * sizeof(T)); // The remaining part of the symbol must have half of the window applied, // sloping down from 1 to 0.5 - for (size_t i = 0; i < d_windowOverlap; i++) { - const size_t out_ix = prefixlength + d_spacing - d_windowOverlap + i; - const size_t in_ix = d_spacing - d_windowOverlap + i; - out[out_ix] = in[in_ix] * d_window[2*d_windowOverlap - (i+1)]; + for (size_t i = 0; i < p.windowOverlap; i++) { + const size_t out_ix = prefixlength + p.spacing - p.windowOverlap + i; + const size_t in_ix = p.spacing - p.windowOverlap + i; + out[out_ix] = in[in_ix] * p.window[2*p.windowOverlap - (i+1)]; } // Suffix is taken from the beginning of the symbol, and sees the other // half of the window applied. - for (size_t i = 0; i < d_windowOverlap; i++) { - const size_t out_ix = prefixlength + d_spacing + i; - out[out_ix] = in[i] * d_window[d_windowOverlap - (i+1)]; + for (size_t i = 0; i < p.windowOverlap; i++) { + const size_t out_ix = prefixlength + p.spacing + i; + out[out_ix] = in[i] * p.window[p.windowOverlap - (i+1)]; } - in += d_spacing; - out += d_nullSize; + in += p.spacing; + out += p.nullSize; // out is now pointing to the proper end of symbol. There are - // d_windowOverlap samples ahead that were already written. + // windowOverlap samples ahead that were already written. } // Data symbols - for (size_t sym_ix = 0; sym_ix < d_nbSymbols; sym_ix++) { + for (size_t sym_ix = 0; sym_ix < p.nbSymbols; sym_ix++) { /* _ix variables are indices into in[], _ox variables are * indices for out[] */ - const ssize_t start_rise_ox = -d_windowOverlap; - const size_t start_rise_ix = 2 * d_spacing - d_symSize - d_windowOverlap; + const ssize_t start_rise_ox = -p.windowOverlap; + const size_t start_rise_ix = 2 * p.spacing - p.symSize - p.windowOverlap; /* const size_t start_real_symbol_ox = 0; - const size_t start_real_symbol_ix = 2 * d_spacing - d_symSize; + const size_t start_real_symbol_ix = 2 * p.spacing - p.symSize; */ - const ssize_t end_rise_ox = d_windowOverlap; - const size_t end_rise_ix = 2 * d_spacing - d_symSize + d_windowOverlap; - const ssize_t end_cyclic_prefix_ox = d_symSize - d_spacing; + const ssize_t end_rise_ox = p.windowOverlap; + const size_t end_rise_ix = 2 * p.spacing - p.symSize + p.windowOverlap; + const ssize_t end_cyclic_prefix_ox = p.symSize - p.spacing; /* end_cyclic_prefix_ix = end of symbol - const size_t begin_fall_ox = d_symSize - d_windowOverlap; - const size_t begin_fall_ix = d_spacing - d_windowOverlap; - const size_t end_real_symbol_ox = d_symSize; + const size_t begin_fall_ox = p.symSize - p.windowOverlap; + const size_t begin_fall_ix = p.spacing - p.windowOverlap; + const size_t end_real_symbol_ox = p.symSize; end_real_symbol_ix = end of symbol - const size_t end_fall_ox = d_symSize + d_windowOverlap; - const size_t end_fall_ix = d_spacing + d_windowOverlap; + const size_t end_fall_ox = p.symSize + p.windowOverlap; + const size_t end_fall_ix = p.spacing + p.windowOverlap; */ ssize_t ox = start_rise_ox; size_t ix = start_rise_ix; for (size_t i = 0; ix < end_rise_ix; i++) { - out[ox] += in[ix] * d_window.at(i); + out[ox] += in[ix] * p.window.at(i); ix++; ox++; } @@ -196,75 +208,88 @@ int GuardIntervalInserter::process(Buffer* const dataIn, Buffer* dataOut) const size_t remaining_prefix_length = end_cyclic_prefix_ox - end_rise_ox; memcpy( &out[ox], &in[ix], - remaining_prefix_length * sizeof(complexf)); + remaining_prefix_length * sizeof(T)); ox += remaining_prefix_length; assert(ox == end_cyclic_prefix_ox); ix = 0; - const bool last_symbol = (sym_ix + 1 >= d_nbSymbols); + const bool last_symbol = (sym_ix + 1 >= p.nbSymbols); if (last_symbol) { // No windowing at all at end - memcpy(&out[ox], &in[ix], d_spacing * sizeof(complexf)); - ox += d_spacing; + memcpy(&out[ox], &in[ix], p.spacing * sizeof(T)); + ox += p.spacing; } else { - // Copy the middle part of the symbol, d_windowOverlap samples + // Copy the middle part of the symbol, p.windowOverlap samples // short of the end. memcpy( &out[ox], &in[ix], - (d_spacing - d_windowOverlap) * sizeof(complexf)); - ox += d_spacing - d_windowOverlap; - ix += d_spacing - d_windowOverlap; - assert(ox == (ssize_t)(d_symSize - d_windowOverlap)); + (p.spacing - p.windowOverlap) * sizeof(T)); + ox += p.spacing - p.windowOverlap; + ix += p.spacing - p.windowOverlap; + assert(ox == (ssize_t)(p.symSize - p.windowOverlap)); // Apply window from 1 to 0.5 for the end of the symbol - for (size_t i = 0; ox < (ssize_t)d_symSize; i++) { - out[ox] = in[ix] * d_window[2*d_windowOverlap - (i+1)]; + for (size_t i = 0; ox < (ssize_t)p.symSize; i++) { + out[ox] = in[ix] * p.window[2*p.windowOverlap - (i+1)]; ox++; ix++; } - assert(ix == d_spacing); + assert(ix == p.spacing); ix = 0; // Cyclic suffix, with window from 0.5 to 0 - for (size_t i = 0; ox < (ssize_t)(d_symSize + d_windowOverlap); i++) { - out[ox] = in[ix] * d_window[d_windowOverlap - (i+1)]; + for (size_t i = 0; ox < (ssize_t)(p.symSize + p.windowOverlap); i++) { + out[ox] = in[ix] * p.window[p.windowOverlap - (i+1)]; ox++; ix++; } - assert(ix == d_windowOverlap); + assert(ix == p.windowOverlap); } - out += d_symSize; - in += d_spacing; + out += p.symSize; + in += p.spacing; // out is now pointing to the proper end of symbol. There are - // d_windowOverlap samples ahead that were already written. + // windowOverlap samples ahead that were already written. } - } + } } else { // Handle Null symbol separately because it is longer // end - (nullSize - spacing) = 2 * spacing - nullSize - memcpy(out, &in[2 * d_spacing - d_nullSize], - (d_nullSize - d_spacing) * sizeof(complexf)); - memcpy(&out[d_nullSize - d_spacing], in, d_spacing * sizeof(complexf)); - in += d_spacing; - out += d_nullSize; + memcpy(out, &in[2 * p.spacing - p.nullSize], + (p.nullSize - p.spacing) * sizeof(T)); + memcpy(&out[p.nullSize - p.spacing], in, p.spacing * sizeof(T)); + in += p.spacing; + out += p.nullSize; // Data symbols - for (size_t i = 0; i < d_nbSymbols; ++i) { + for (size_t i = 0; i < p.nbSymbols; ++i) { // end - (symSize - spacing) = 2 * spacing - symSize - memcpy(out, &in[2 * d_spacing - d_symSize], - (d_symSize - d_spacing) * sizeof(complexf)); - memcpy(&out[d_symSize - d_spacing], in, d_spacing * sizeof(complexf)); - in += d_spacing; - out += d_symSize; + memcpy(out, &in[2 * p.spacing - p.symSize], + (p.symSize - p.spacing) * sizeof(T)); + memcpy(&out[p.symSize - p.spacing], in, p.spacing * sizeof(T)); + in += p.spacing; + out += p.symSize; } } return sizeIn; } +int GuardIntervalInserter::process(Buffer* const dataIn, Buffer* dataOut) +{ + if (m_fixedPoint) { + if (m_params.windowOverlap) { + throw std::runtime_error("fixed point and ofdm windowing not supported"); + } + return do_process<complexfix>(m_params, dataIn, dataOut); + } + else { + return do_process<complexf>(m_params, dataIn, dataOut); + } +} + void GuardIntervalInserter::set_parameter( const std::string& parameter, const std::string& value) @@ -291,7 +316,7 @@ const std::string GuardIntervalInserter::get_parameter(const std::string& parame using namespace std; stringstream ss; if (parameter == "windowlen") { - ss << d_windowOverlap; + ss << m_params.windowOverlap; } else { ss << "Parameter '" << parameter << @@ -304,6 +329,6 @@ const std::string GuardIntervalInserter::get_parameter(const std::string& parame const json::map_t GuardIntervalInserter::get_all_values() const { json::map_t map; - map["windowlen"].v = d_windowOverlap; + map["windowlen"].v = m_params.windowOverlap; return map; } diff --git a/src/GuardIntervalInserter.h b/src/GuardIntervalInserter.h index f78ac91..380142e 100644 --- a/src/GuardIntervalInserter.h +++ b/src/GuardIntervalInserter.h @@ -50,7 +50,8 @@ class GuardIntervalInserter : public ModCodec, public RemoteControllable size_t spacing, size_t nullSize, size_t symSize, - size_t& windowOverlap); + size_t& windowOverlap, + bool fixedPoint); virtual ~GuardIntervalInserter() {} @@ -62,16 +63,30 @@ class GuardIntervalInserter : public ModCodec, public RemoteControllable virtual const std::string get_parameter(const std::string& parameter) const override; virtual const json::map_t get_all_values() const override; + struct Params { + Params( + size_t nbSymbols, + size_t spacing, + size_t nullSize, + size_t symSize, + size_t& windowOverlap); + + size_t nbSymbols; + size_t spacing; + size_t nullSize; + size_t symSize; + size_t& windowOverlap; + + mutable std::mutex windowMutex; + std::vector<float> window; + }; + protected: void update_window(size_t new_window_overlap); - size_t d_nbSymbols; - size_t d_spacing; - size_t d_nullSize; - size_t d_symSize; + bool m_fixedPoint; + + Params m_params; - mutable std::mutex d_windowMutex; - size_t& d_windowOverlap; - std::vector<float> d_window; }; diff --git a/src/NullSymbol.cpp b/src/NullSymbol.cpp index 023ae42..526e662 100644 --- a/src/NullSymbol.cpp +++ b/src/NullSymbol.cpp @@ -31,11 +31,12 @@ #include <cstdlib> #include <cstring> -NullSymbol::NullSymbol(size_t nbCarriers) : +NullSymbol::NullSymbol(size_t numCarriers, size_t typeSize) : ModInput(), - myNbCarriers(nbCarriers) + m_numCarriers(numCarriers), + m_typeSize(typeSize) { - PDEBUG("NullSymbol::NullSymbol(%zu) @ %p\n", nbCarriers, this); + PDEBUG("NullSymbol::NullSymbol(%zu) @ %p\n", numCarriers, this); } @@ -49,7 +50,7 @@ int NullSymbol::process(Buffer* dataOut) { PDEBUG("NullSymbol::process(dataOut: %p)\n", dataOut); - dataOut->setLength(myNbCarriers * 2 * sizeof(float)); + dataOut->setLength(m_numCarriers * m_typeSize); memset(dataOut->getData(), 0, dataOut->getLength()); return dataOut->getLength(); diff --git a/src/NullSymbol.h b/src/NullSymbol.h index 814e434..6ba9e63 100644 --- a/src/NullSymbol.h +++ b/src/NullSymbol.h @@ -39,14 +39,14 @@ class NullSymbol : public ModInput { public: - NullSymbol(size_t nbCarriers); + NullSymbol(size_t nunCarriers, size_t typeSize); virtual ~NullSymbol(); int process(Buffer* dataOut); const char* name() { return "NullSymbol"; } private: - size_t myNbCarriers; - + size_t m_numCarriers; + size_t m_typeSize; }; diff --git a/src/OfdmGenerator.cpp b/src/OfdmGenerator.cpp index cb799d3..e679694 100644 --- a/src/OfdmGenerator.cpp +++ b/src/OfdmGenerator.cpp @@ -27,17 +27,19 @@ #include "OfdmGenerator.h" #include "PcDebug.h" -#define FFT_TYPE fftwf_complex - -#include <string.h> #include <stdexcept> #include <assert.h> #include <string> #include <numeric> +#include <vector> +#include <cstring> +#include <complex> static const size_t MAX_CLIP_STATS = 10; -OfdmGenerator::OfdmGenerator(size_t nbSymbols, +using FFTW_TYPE = fftwf_complex; + +OfdmGeneratorCF32::OfdmGeneratorCF32(size_t nbSymbols, size_t nbCarriers, size_t spacing, bool& enableCfr, @@ -102,29 +104,29 @@ OfdmGenerator::OfdmGenerator(size_t nbSymbols, PDEBUG(" myZeroSize: %u\n", myZeroSize); const int N = mySpacing; // The size of the FFT - myFftIn = (FFT_TYPE*)fftwf_malloc(sizeof(FFT_TYPE) * N); - myFftOut = (FFT_TYPE*)fftwf_malloc(sizeof(FFT_TYPE) * N); + myFftIn = (FFTW_TYPE*)fftwf_malloc(sizeof(FFTW_TYPE) * N); + myFftOut = (FFTW_TYPE*)fftwf_malloc(sizeof(FFTW_TYPE) * N); fftwf_set_timelimit(2); myFftPlan = fftwf_plan_dft_1d(N, myFftIn, myFftOut, FFTW_BACKWARD, FFTW_MEASURE); - myCfrPostClip = (FFT_TYPE*)fftwf_malloc(sizeof(FFT_TYPE) * N); - myCfrPostFft = (FFT_TYPE*)fftwf_malloc(sizeof(FFT_TYPE) * N); + myCfrPostClip = (FFTW_TYPE*)fftwf_malloc(sizeof(FFTW_TYPE) * N); + myCfrPostFft = (FFTW_TYPE*)fftwf_malloc(sizeof(FFTW_TYPE) * N); myCfrFft = fftwf_plan_dft_1d(N, myCfrPostClip, myCfrPostFft, FFTW_FORWARD, FFTW_MEASURE); - if (sizeof(complexf) != sizeof(FFT_TYPE)) { + if (sizeof(complexf) != sizeof(FFTW_TYPE)) { printf("sizeof(complexf) %zu\n", sizeof(complexf)); - printf("sizeof(FFT_TYPE) %zu\n", sizeof(FFT_TYPE)); + printf("sizeof(FFT_TYPE) %zu\n", sizeof(FFTW_TYPE)); throw std::runtime_error( "OfdmGenerator::process complexf size is not FFT_TYPE size!"); } } -OfdmGenerator::~OfdmGenerator() +OfdmGeneratorCF32::~OfdmGeneratorCF32() { PDEBUG("OfdmGenerator::~OfdmGenerator() @ %p\n", this); @@ -153,15 +155,15 @@ OfdmGenerator::~OfdmGenerator() } } -int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut) +int OfdmGeneratorCF32::process(Buffer* const dataIn, Buffer* dataOut) { PDEBUG("OfdmGenerator::process(dataIn: %p, dataOut: %p)\n", dataIn, dataOut); dataOut->setLength(myNbSymbols * mySpacing * sizeof(complexf)); - FFT_TYPE* in = reinterpret_cast<FFT_TYPE*>(dataIn->getData()); - FFT_TYPE* out = reinterpret_cast<FFT_TYPE*>(dataOut->getData()); + FFTW_TYPE* in = reinterpret_cast<FFTW_TYPE*>(dataIn->getData()); + FFTW_TYPE* out = reinterpret_cast<FFTW_TYPE*>(dataOut->getData()); size_t sizeIn = dataIn->getLength() / sizeof(complexf); size_t sizeOut = dataOut->getLength() / sizeof(complexf); @@ -212,17 +214,17 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut) * PosSrc=0 PosDst=1 PosSize=768 * NegSrc=768 NegDst=1280 NegSize=768 */ - memset(&myFftIn[myZeroDst], 0, myZeroSize * sizeof(FFT_TYPE)); + memset(&myFftIn[myZeroDst], 0, myZeroSize * sizeof(FFTW_TYPE)); memcpy(&myFftIn[myPosDst], &in[myPosSrc], - myPosSize * sizeof(FFT_TYPE)); + myPosSize * sizeof(FFTW_TYPE)); memcpy(&myFftIn[myNegDst], &in[myNegSrc], - myNegSize * sizeof(FFT_TYPE)); + myNegSize * sizeof(FFTW_TYPE)); if (myCfr) { reference.resize(mySpacing); memcpy(reinterpret_cast<fftwf_complex*>(reference.data()), - myFftIn, mySpacing * sizeof(FFT_TYPE)); + myFftIn, mySpacing * sizeof(FFTW_TYPE)); } fftwf_execute(myFftPlan); // IFFT from myFftIn to myFftOut @@ -235,7 +237,7 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut) if (myMERCalcIndex == i) { before_cfr.resize(mySpacing); memcpy(reinterpret_cast<fftwf_complex*>(before_cfr.data()), - myFftOut, mySpacing * sizeof(FFT_TYPE)); + myFftOut, mySpacing * sizeof(FFTW_TYPE)); } /* cfr_one_iteration runs the myFftPlan again at the end, and @@ -277,7 +279,7 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut) num_error_clip += stat.errclip_count; } - memcpy(out, myFftOut, mySpacing * sizeof(FFT_TYPE)); + memcpy(out, myFftOut, mySpacing * sizeof(FFTW_TYPE)); in += myNbCarriers; out += mySpacing; @@ -308,14 +310,14 @@ int OfdmGenerator::process(Buffer* const dataIn, Buffer* dataOut) return sizeOut; } -OfdmGenerator::cfr_iter_stat_t OfdmGenerator::cfr_one_iteration( +OfdmGeneratorCF32::cfr_iter_stat_t OfdmGeneratorCF32::cfr_one_iteration( complexf *symbol, const complexf *reference) { // use std::norm instead of std::abs to avoid calculating the // square roots const float clip_squared = myCfrClip * myCfrClip; - OfdmGenerator::cfr_iter_stat_t ret; + OfdmGeneratorCF32::cfr_iter_stat_t ret; // Clip for (size_t i = 0; i < mySpacing; i++) { @@ -331,7 +333,7 @@ OfdmGenerator::cfr_iter_stat_t OfdmGenerator::cfr_one_iteration( } // Take FFT of our clipped signal - memcpy(myCfrPostClip, symbol, mySpacing * sizeof(FFT_TYPE)); + memcpy(myCfrPostClip, symbol, mySpacing * sizeof(FFTW_TYPE)); fftwf_execute(myCfrFft); // FFT from myCfrPostClip to myCfrPostFft // Calculate the error in frequency domain by subtracting our reference @@ -374,7 +376,7 @@ OfdmGenerator::cfr_iter_stat_t OfdmGenerator::cfr_one_iteration( } -void OfdmGenerator::set_parameter(const std::string& parameter, +void OfdmGeneratorCF32::set_parameter(const std::string& parameter, const std::string& value) { using namespace std; @@ -404,7 +406,7 @@ void OfdmGenerator::set_parameter(const std::string& parameter, } } -const std::string OfdmGenerator::get_parameter(const std::string& parameter) const +const std::string OfdmGeneratorCF32::get_parameter(const std::string& parameter) const { using namespace std; stringstream ss; @@ -458,9 +460,127 @@ const std::string OfdmGenerator::get_parameter(const std::string& parameter) con return ss.str(); } -const json::map_t OfdmGenerator::get_all_values() const +const json::map_t OfdmGeneratorCF32::get_all_values() const { json::map_t map; // TODO needs rework of the values return map; } + +OfdmGeneratorFixed::OfdmGeneratorFixed(size_t nbSymbols, + size_t nbCarriers, + size_t spacing, + bool& enableCfr, + float& cfrClip, + float& cfrErrorClip, + bool inverse) : + ModCodec(), + myNbSymbols(nbSymbols), + myNbCarriers(nbCarriers), + mySpacing(spacing) +{ + PDEBUG("OfdmGenerator::OfdmGenerator(%zu, %zu, %zu, %s) @ %p\n", + nbSymbols, nbCarriers, spacing, inverse ? "true" : "false", this); + + etiLog.level(info) << "Using KISS FFT by Mark Borgerding for fixed-point transform"; + + if (nbCarriers > spacing) { + throw std::runtime_error( + "OfdmGenerator::OfdmGenerator nbCarriers > spacing!"); + } + + if (inverse) { + myPosDst = (nbCarriers & 1 ? 0 : 1); + myPosSrc = 0; + myPosSize = (nbCarriers + 1) / 2; + myNegDst = spacing - (nbCarriers / 2); + myNegSrc = (nbCarriers + 1) / 2; + myNegSize = nbCarriers / 2; + } + else { + myPosDst = (nbCarriers & 1 ? 0 : 1); + myPosSrc = nbCarriers / 2; + myPosSize = (nbCarriers + 1) / 2; + myNegDst = spacing - (nbCarriers / 2); + myNegSrc = 0; + myNegSize = nbCarriers / 2; + } + myZeroDst = myPosDst + myPosSize; + myZeroSize = myNegDst - myZeroDst; + + PDEBUG(" myPosDst: %u\n", myPosDst); + PDEBUG(" myPosSrc: %u\n", myPosSrc); + PDEBUG(" myPosSize: %u\n", myPosSize); + PDEBUG(" myNegDst: %u\n", myNegDst); + PDEBUG(" myNegSrc: %u\n", myNegSrc); + PDEBUG(" myNegSize: %u\n", myNegSize); + PDEBUG(" myZeroDst: %u\n", myZeroDst); + PDEBUG(" myZeroSize: %u\n", myZeroSize); + + const int N = mySpacing; // The size of the FFT + + const size_t nbytes = N * sizeof(kiss_fft_cpx); + myFftIn = (kiss_fft_cpx*)KISS_FFT_MALLOC(nbytes); + myFftOut = (kiss_fft_cpx*)KISS_FFT_MALLOC(nbytes); + memset(myFftIn, 0, nbytes); + + myKissCfg = kiss_fft_alloc(N, inverse, nullptr, nullptr); +} + +OfdmGeneratorFixed::~OfdmGeneratorFixed() +{ + if (myKissCfg) KISS_FFT_FREE(myKissCfg); + if (myFftIn) KISS_FFT_FREE(myFftIn); + if (myFftOut) KISS_FFT_FREE(myFftOut); +} + +int OfdmGeneratorFixed::process(Buffer* const dataIn, Buffer* dataOut) +{ + dataOut->setLength(myNbSymbols * mySpacing * sizeof(kiss_fft_cpx)); + + kiss_fft_cpx* in = reinterpret_cast<kiss_fft_cpx*>(dataIn->getData()); + kiss_fft_cpx* out = reinterpret_cast<kiss_fft_cpx*>(dataOut->getData()); + + size_t sizeIn = dataIn->getLength() / sizeof(kiss_fft_cpx); + size_t sizeOut = dataOut->getLength() / sizeof(kiss_fft_cpx); + + if (sizeIn != myNbSymbols * myNbCarriers) { + PDEBUG("Nb symbols: %zu\n", myNbSymbols); + PDEBUG("Nb carriers: %zu\n", myNbCarriers); + PDEBUG("Spacing: %zu\n", mySpacing); + PDEBUG("\n%zu != %zu\n", sizeIn, myNbSymbols * myNbCarriers); + throw std::runtime_error( + "OfdmGenerator::process input size not valid!"); + } + if (sizeOut != myNbSymbols * mySpacing) { + PDEBUG("Nb symbols: %zu\n", myNbSymbols); + PDEBUG("Nb carriers: %zu\n", myNbCarriers); + PDEBUG("Spacing: %zu\n", mySpacing); + PDEBUG("\n%zu != %zu\n", sizeIn, myNbSymbols * mySpacing); + throw std::runtime_error( + "OfdmGenerator::process output size not valid!"); + } + + for (size_t i = 0; i < myNbSymbols; ++i) { + myFftIn[0].r = 0; + myFftIn[0].i = 0; + + /* For TM I this is: + * ZeroDst=769 ZeroSize=511 + * PosSrc=0 PosDst=1 PosSize=768 + * NegSrc=768 NegDst=1280 NegSize=768 + */ + memset(&myFftIn[myZeroDst], 0, myZeroSize * sizeof(kiss_fft_cpx)); + memcpy(&myFftIn[myPosDst], &in[myPosSrc], myPosSize * sizeof(kiss_fft_cpx)); + memcpy(&myFftIn[myNegDst], &in[myNegSrc], myNegSize * sizeof(kiss_fft_cpx)); + + kiss_fft(myKissCfg, myFftIn, myFftOut); + + memcpy(out, myFftOut, mySpacing * sizeof(kiss_fft_cpx)); + + in += myNbCarriers; + out += mySpacing; + } + + return sizeOut; +} diff --git a/src/OfdmGenerator.h b/src/OfdmGenerator.h index cab936e..2e1aa63 100644 --- a/src/OfdmGenerator.h +++ b/src/OfdmGenerator.h @@ -33,25 +33,26 @@ #include "ModPlugin.h" #include "RemoteControl.h" #include "PAPRStats.h" -#include "fftw3.h" +#include "kiss_fft.h" + +#include <fftw3.h> #include <cstddef> -#include <vector> -#include <complex> #include <atomic> -class OfdmGenerator : public ModCodec, public RemoteControllable +// Complex Float uses FFTW +class OfdmGeneratorCF32 : public ModCodec, public RemoteControllable { public: - OfdmGenerator(size_t nbSymbols, + OfdmGeneratorCF32(size_t nbSymbols, size_t nbCarriers, size_t spacing, bool& enableCfr, float& cfrClip, float& cfrErrorClip, bool inverse = true); - virtual ~OfdmGenerator(); - OfdmGenerator(const OfdmGenerator&) = delete; - OfdmGenerator& operator=(const OfdmGenerator&) = delete; + virtual ~OfdmGeneratorCF32(); + OfdmGeneratorCF32(const OfdmGeneratorCF32&) = delete; + OfdmGeneratorCF32& operator=(const OfdmGeneratorCF32&) = delete; int process(Buffer* const dataIn, Buffer* dataOut) override; const char* name() override { return "OfdmGenerator"; } @@ -105,4 +106,37 @@ class OfdmGenerator : public ModCodec, public RemoteControllable std::deque<double> myMERs; }; +// Fixed point implementation uses KISS FFT with -DFIXED_POINT=32 +class OfdmGeneratorFixed : public ModCodec +{ + public: + OfdmGeneratorFixed(size_t nbSymbols, + size_t nbCarriers, + size_t spacing, + bool& enableCfr, + float& cfrClip, + float& cfrErrorClip, + bool inverse = true); + virtual ~OfdmGeneratorFixed(); + OfdmGeneratorFixed(const OfdmGeneratorFixed&) = delete; + OfdmGeneratorFixed& operator=(const OfdmGeneratorFixed&) = delete; + + int process(Buffer* const dataIn, Buffer* dataOut) override; + const char* name() override { return "OfdmGenerator"; } + + private: + kiss_fft_cfg myKissCfg = nullptr; + kiss_fft_cpx *myFftIn, *myFftOut; + const size_t myNbSymbols; + const size_t myNbCarriers; + const size_t mySpacing; + unsigned myPosSrc; + unsigned myPosDst; + unsigned myPosSize; + unsigned myNegSrc; + unsigned myNegDst; + unsigned myNegSize; + unsigned myZeroDst; + unsigned myZeroSize; +}; diff --git a/src/OutputMemory.cpp b/src/OutputMemory.cpp index d6ef917..ac8a67b 100644 --- a/src/OutputMemory.cpp +++ b/src/OutputMemory.cpp @@ -39,7 +39,7 @@ OutputMemory::OutputMemory(Buffer* dataOut) { PDEBUG("OutputMemory::OutputMemory(%p) @ %p\n", dataOut, this); - setOutput(dataOut); + myDataOut = dataOut; #if OUTPUT_MEM_HISTOGRAM myMax = 0.0f; @@ -67,12 +67,6 @@ OutputMemory::~OutputMemory() } -void OutputMemory::setOutput(Buffer* dataOut) -{ - myDataOut = dataOut; -} - - int OutputMemory::process(Buffer* dataIn) { PDEBUG("OutputMemory::process(dataIn: %p)\n", diff --git a/src/OutputMemory.h b/src/OutputMemory.h index f0a5fbb..e7252d3 100644 --- a/src/OutputMemory.h +++ b/src/OutputMemory.h @@ -61,8 +61,6 @@ public: meta_vec_t get_latest_metadata(void); - void setOutput(Buffer* dataOut); - protected: Buffer* myDataOut; meta_vec_t myMetadata; diff --git a/src/PhaseReference.cpp b/src/PhaseReference.cpp index 568e15e..d7b89bf 100644 --- a/src/PhaseReference.cpp +++ b/src/PhaseReference.cpp @@ -29,12 +29,10 @@ #include <stdexcept> -using complexf = std::complex<float>; - /* ETSI EN 300 401 Table 43 (Clause 14.3.2) * Contains h_{i,k} values */ -const uint8_t PhaseReference::d_h[4][32] = { +static const uint8_t d_h[4][32] = { /* h0 */ { 0, 2, 0, 0, 0, 0, 1, 1, 2, 0, 0, 0, 2, 2, 1, 1, 0, 2, 0, 0, 0, 0, 1, 1, 2, 0, 0, 0, 2, 2, 1, 1 }, /* h1 */ { 0, 3, 2, 3, 0, 1, 3, 0, 2, 1, 2, 3, 2, 3, 3, 0, @@ -54,41 +52,80 @@ const uint8_t PhaseReference::d_h[4][32] = { * Tables 44 to 47 describe the frequency interleaving done in * FrequencyInterleaver. */ -PhaseReference::PhaseReference(unsigned int dabmode) : +PhaseReference::PhaseReference(unsigned int dabmode, bool fixedPoint) : ModInput(), - d_dabmode(dabmode) + d_dabmode(dabmode), + d_fixedPoint(fixedPoint) { PDEBUG("PhaseReference::PhaseReference(%u) @ %p\n", dabmode, this); switch (d_dabmode) { case 1: d_carriers = 1536; - d_num = 2048; break; case 2: d_carriers = 384; - d_num = 512; break; case 3: d_carriers = 192; - d_num = 256; break; case 4: d_dabmode = 0; case 0: d_carriers = 768; - d_num = 1024; break; default: throw std::runtime_error( "PhaseReference::PhaseReference DAB mode not valid!"); } - d_dataIn.resize(d_carriers); - fillData(); + + if (d_fixedPoint) { + d_phaseRefFixed.fillData(d_dabmode, d_carriers); + } + else { + d_phaseRefCF32.fillData(d_dabmode, d_carriers); + } } -complexf convert(uint8_t data) { +static const int table[][48][2] = { + { // Mode 0/4 + // Positive part + { 0, 0 }, { 3, 1 }, { 2, 0 }, { 1, 2 }, { 0, 0 }, { 3, 1 }, + { 2, 2 }, { 1, 2 }, { 0, 2 }, { 3, 1 }, { 2, 3 }, { 1, 0 }, + // Negative part + { 0, 0 }, { 1, 1 }, { 2, 1 }, { 3, 2 }, { 0, 2 }, { 1, 2 }, + { 2, 0 }, { 3, 3 }, { 0, 3 }, { 1, 1 }, { 2, 3 }, { 3, 2 }, + }, + { // Mode 1 + // Positive part + { 0, 3 }, { 3, 1 }, { 2, 1 }, { 1, 1 }, { 0, 2 }, { 3, 2 }, + { 2, 1 }, { 1, 0 }, { 0, 2 }, { 3, 2 }, { 2, 3 }, { 1, 3 }, + { 0, 0 }, { 3, 2 }, { 2, 1 }, { 1, 3 }, { 0, 3 }, { 3, 3 }, + { 2, 3 }, { 1, 0 }, { 0, 3 }, { 3, 0 }, { 2, 1 }, { 1, 1 }, + // Negative part + { 0, 1 }, { 1, 2 }, { 2, 0 }, { 3, 1 }, { 0, 3 }, { 1, 2 }, + { 2, 2 }, { 3, 3 }, { 0, 2 }, { 1, 1 }, { 2, 2 }, { 3, 3 }, + { 0, 1 }, { 1, 2 }, { 2, 3 }, { 3, 3 }, { 0, 2 }, { 1, 2 }, + { 2, 2 }, { 3, 1 }, { 0, 1 }, { 1, 3 }, { 2, 1 }, { 3, 2 }, + }, + { // Mode 2 + // Positive part + { 2, 0 }, { 1, 2 }, { 0, 2 }, { 3, 1 }, { 2, 0 }, { 1, 3 }, + // Negative part + { 0, 2 }, { 1, 3 }, { 2, 2 }, { 3, 2 }, { 0, 1 }, { 1, 2 }, + }, + { // Mode 3 + // Positive part + { 3, 2 }, { 2, 2 }, { 1, 2 }, + // Negative part + { 0, 2 }, { 1, 3 }, { 2, 0 }, + }, +}; + + +template <> +complexf PhaseRefGen<complexf>::convert(uint8_t data) { const complexf value[] = { complexf(1, 0), complexf(0, 1), @@ -98,62 +135,37 @@ complexf convert(uint8_t data) { return value[data % 4]; } +template <> +complexfix PhaseRefGen<complexfix>::convert(uint8_t data) { + constexpr auto one = fpm::fixed_16_16{1}; + constexpr auto zero = fpm::fixed_16_16{0}; -void PhaseReference::fillData() -{ - const int table[][48][2] = { - { // Mode 0/4 - // Positive part - { 0, 0 }, { 3, 1 }, { 2, 0 }, { 1, 2 }, { 0, 0 }, { 3, 1 }, - { 2, 2 }, { 1, 2 }, { 0, 2 }, { 3, 1 }, { 2, 3 }, { 1, 0 }, - // Negative part - { 0, 0 }, { 1, 1 }, { 2, 1 }, { 3, 2 }, { 0, 2 }, { 1, 2 }, - { 2, 0 }, { 3, 3 }, { 0, 3 }, { 1, 1 }, { 2, 3 }, { 3, 2 }, - }, - { // Mode 1 - // Positive part - { 0, 3 }, { 3, 1 }, { 2, 1 }, { 1, 1 }, { 0, 2 }, { 3, 2 }, - { 2, 1 }, { 1, 0 }, { 0, 2 }, { 3, 2 }, { 2, 3 }, { 1, 3 }, - { 0, 0 }, { 3, 2 }, { 2, 1 }, { 1, 3 }, { 0, 3 }, { 3, 3 }, - { 2, 3 }, { 1, 0 }, { 0, 3 }, { 3, 0 }, { 2, 1 }, { 1, 1 }, - // Negative part - { 0, 1 }, { 1, 2 }, { 2, 0 }, { 3, 1 }, { 0, 3 }, { 1, 2 }, - { 2, 2 }, { 3, 3 }, { 0, 2 }, { 1, 1 }, { 2, 2 }, { 3, 3 }, - { 0, 1 }, { 1, 2 }, { 2, 3 }, { 3, 3 }, { 0, 2 }, { 1, 2 }, - { 2, 2 }, { 3, 1 }, { 0, 1 }, { 1, 3 }, { 2, 1 }, { 3, 2 }, - }, - { // Mode 2 - // Positive part - { 2, 0 }, { 1, 2 }, { 0, 2 }, { 3, 1 }, { 2, 0 }, { 1, 3 }, - // Negative part - { 0, 2 }, { 1, 3 }, { 2, 2 }, { 3, 2 }, { 0, 1 }, { 1, 2 }, - }, - { // Mode 3 - // Positive part - { 3, 2 }, { 2, 2 }, { 1, 2 }, - // Negative part - { 0, 2 }, { 1, 3 }, { 2, 0 }, - }, + const complexfix value[] = { + complexfix(one, zero), + complexfix(zero, one), + complexfix(-one, zero), + complexfix(zero, -one), }; + return value[data % 4]; +} - if (d_dabmode > 3) { - throw std::runtime_error( - "PhaseReference::fillData invalid DAB mode!"); - } - - if (d_dataIn.size() != d_carriers) { +template <typename T> +void PhaseRefGen<T>::fillData(unsigned int dabmode, size_t carriers) +{ + dataIn.resize(carriers); + if (dataIn.size() != carriers) { throw std::runtime_error( - "PhaseReference::fillData d_dataIn has incorrect size!"); + "PhaseReference::fillData dataIn has incorrect size!"); } for (size_t index = 0, offset = 0; - index < d_dataIn.size(); + index < dataIn.size(); ++offset) { for (size_t k = 0; k < 32; ++k) { - d_dataIn[index++] = convert( - d_h[ table[d_dabmode][offset][0] ][k] + - table[d_dabmode][offset][1] ); + dataIn[index++] = convert( + d_h[ table[dabmode][offset][0] ][k] + + table[dabmode][offset][1] ); } } } @@ -163,7 +175,12 @@ int PhaseReference::process(Buffer* dataOut) { PDEBUG("PhaseReference::process(dataOut: %p)\n", dataOut); - dataOut->setData(&d_dataIn[0], d_carriers * sizeof(complexf)); + if (d_fixedPoint) { + dataOut->setData(&d_phaseRefFixed.dataIn[0], d_carriers * sizeof(complexfix)); + } + else { + dataOut->setData(&d_phaseRefCF32.dataIn[0], d_carriers * sizeof(complexf)); + } return 1; } diff --git a/src/PhaseReference.h b/src/PhaseReference.h index 6ecdc4e..735009c 100644 --- a/src/PhaseReference.h +++ b/src/PhaseReference.h @@ -32,25 +32,33 @@ #include "ModPlugin.h" -#include <cstddef> -#include <complex> #include <vector> +#include <cstddef> + +template <typename T> +struct PhaseRefGen { + std::vector<T> dataIn; + void fillData(unsigned int dabmode, size_t carriers); + + private: + T convert(uint8_t data); +}; + class PhaseReference : public ModInput { public: - PhaseReference(unsigned int dabmode); + PhaseReference(unsigned int dabmode, bool fixedPoint); int process(Buffer* dataOut) override; const char* name() override { return "PhaseReference"; } protected: unsigned int d_dabmode; + bool d_fixedPoint; size_t d_carriers; - size_t d_num; - const static uint8_t d_h[4][32]; - std::vector<std::complex<float> > d_dataIn; - void fillData(); + PhaseRefGen<complexf> d_phaseRefCF32; + PhaseRefGen<complexfix> d_phaseRefFixed; }; diff --git a/src/QpskSymbolMapper.cpp b/src/QpskSymbolMapper.cpp index c9b01fb..c12ad80 100644 --- a/src/QpskSymbolMapper.cpp +++ b/src/QpskSymbolMapper.cpp @@ -31,9 +31,10 @@ #include "QpskSymbolMapper.h" #include "PcDebug.h" -QpskSymbolMapper::QpskSymbolMapper(size_t carriers) : +QpskSymbolMapper::QpskSymbolMapper(size_t carriers, bool fixedPoint) : ModCodec(), - d_carriers(carriers) { } + m_fixedPoint(fixedPoint), + m_carriers(carriers) { } int QpskSymbolMapper::process(Buffer* const dataIn, Buffer* dataOut) { @@ -41,112 +42,172 @@ int QpskSymbolMapper::process(Buffer* const dataIn, Buffer* dataOut) "(dataIn: %p, dataOut: %p)\n", dataIn, dataOut); - dataOut->setLength(dataIn->getLength() * 4 * 2 * sizeof(float)); // 4 output complex symbols per input byte -#ifdef __SSE__ - const uint8_t* in = reinterpret_cast<const uint8_t*>(dataIn->getData()); - __m128* out = reinterpret_cast<__m128*>(dataOut->getData()); - - if (dataIn->getLength() % (d_carriers / 4) != 0) { - throw std::runtime_error( - "QpskSymbolMapper::process input size not valid: " + - std::to_string(dataIn->getLength()) + - "(input size) % (" + std::to_string(d_carriers) + - " (carriers) / 4) != 0"); - } + // 4 output complex symbols per input byte + + if (m_fixedPoint) { + dataOut->setLength(dataIn->getLength() * 4 * sizeof(complexfix)); + + using fixed_t = complexfix::value_type; + + const uint8_t* in = reinterpret_cast<const uint8_t*>(dataIn->getData()); + fixed_t* out = reinterpret_cast<fixed_t*>(dataOut->getData()); - const static __m128 symbols[16] = { - _mm_setr_ps( M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, M_SQRT1_2), - _mm_setr_ps( M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2), - _mm_setr_ps( M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2, M_SQRT1_2), - _mm_setr_ps( M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2), - _mm_setr_ps( M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2), - _mm_setr_ps( M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2), - _mm_setr_ps( M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2), - _mm_setr_ps( M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2), - _mm_setr_ps(-M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, M_SQRT1_2), - _mm_setr_ps(-M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2), - _mm_setr_ps(-M_SQRT1_2,- M_SQRT1_2, M_SQRT1_2, M_SQRT1_2), - _mm_setr_ps(-M_SQRT1_2,- M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2), - _mm_setr_ps(-M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2), - _mm_setr_ps(-M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2), - _mm_setr_ps(-M_SQRT1_2,- M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2), - _mm_setr_ps(-M_SQRT1_2,- M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2) - }; - size_t inOffset = 0; - size_t outOffset = 0; - uint8_t tmp = 0; - for (size_t i = 0; i < dataIn->getLength(); i += d_carriers / 4) { - for (size_t j = 0; j < d_carriers / 8; ++j) { - tmp = (in[inOffset] & 0xc0) >> 4; - tmp |= (in[inOffset + (d_carriers / 8)] & 0xc0) >> 6; - out[outOffset] = symbols[tmp]; - tmp = (in[inOffset] & 0x30) >> 2; - tmp |= (in[inOffset + (d_carriers / 8)] & 0x30) >> 4; - out[outOffset + 1] = symbols[tmp]; - tmp = (in[inOffset] & 0x0c); - tmp |= (in[inOffset + (d_carriers / 8)] & 0x0c) >> 2; - out[outOffset + 2] = symbols[tmp]; - tmp = (in[inOffset] & 0x03) << 2; - tmp |= (in[inOffset + (d_carriers / 8)] & 0x03); - out[outOffset + 3] = symbols[tmp]; - ++inOffset; - outOffset += 4; + if (dataIn->getLength() % (m_carriers / 4) != 0) { + throw std::runtime_error( + "QpskSymbolMapper::process input size not valid!"); + } + + constexpr fixed_t v = static_cast<fixed_t>(M_SQRT1_2); + + const static fixed_t symbols[16][4] = { + { v, v, v, v}, + { v, v, v, -v}, + { v, -v, v, v}, + { v, -v, v, -v}, + { v, v, -v, v}, + { v, v, -v, -v}, + { v, -v, -v, v}, + { v, -v, -v, -v}, + {-v, v, v, v}, + {-v, v, v, -v}, + {-v, -v, v, v}, + {-v, -v, v, -v}, + {-v, v, -v, v}, + {-v, v, -v, -v}, + {-v, -v, -v, v}, + {-v, -v, -v, -v} + }; + size_t inOffset = 0; + size_t outOffset = 0; + uint8_t tmp; + for (size_t i = 0; i < dataIn->getLength(); i += m_carriers / 4) { + for (size_t j = 0; j < m_carriers / 8; ++j) { + tmp = (in[inOffset] & 0xc0) >> 4; + tmp |= (in[inOffset + (m_carriers / 8)] & 0xc0) >> 6; + memcpy(&out[outOffset], symbols[tmp], sizeof(fixed_t) * 4); + tmp = (in[inOffset] & 0x30) >> 2; + tmp |= (in[inOffset + (m_carriers / 8)] & 0x30) >> 4; + memcpy(&out[outOffset + 4], symbols[tmp], sizeof(fixed_t) * 4); + tmp = (in[inOffset] & 0x0c); + tmp |= (in[inOffset + (m_carriers / 8)] & 0x0c) >> 2; + memcpy(&out[outOffset + 8], symbols[tmp], sizeof(fixed_t) * 4); + tmp = (in[inOffset] & 0x03) << 2; + tmp |= (in[inOffset + (m_carriers / 8)] & 0x03); + memcpy(&out[outOffset + 12], symbols[tmp], sizeof(fixed_t) * 4); + ++inOffset; + outOffset += 4*4; + } + inOffset += m_carriers / 8; } - inOffset += d_carriers / 8; } + else { + dataOut->setLength(dataIn->getLength() * 4 * sizeof(complexf)); +#ifdef __SSE__ + const uint8_t* in = reinterpret_cast<const uint8_t*>(dataIn->getData()); + __m128* out = reinterpret_cast<__m128*>(dataOut->getData()); + + if (dataIn->getLength() % (m_carriers / 4) != 0) { + throw std::runtime_error( + "QpskSymbolMapper::process input size not valid: " + + std::to_string(dataIn->getLength()) + + "(input size) % (" + std::to_string(m_carriers) + + " (carriers) / 4) != 0"); + } + + const static __m128 symbols[16] = { + _mm_setr_ps( M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, M_SQRT1_2), + _mm_setr_ps( M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2), + _mm_setr_ps( M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2, M_SQRT1_2), + _mm_setr_ps( M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2), + _mm_setr_ps( M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2), + _mm_setr_ps( M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2), + _mm_setr_ps( M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2), + _mm_setr_ps( M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2), + _mm_setr_ps(-M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, M_SQRT1_2), + _mm_setr_ps(-M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2), + _mm_setr_ps(-M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2, M_SQRT1_2), + _mm_setr_ps(-M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2), + _mm_setr_ps(-M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2), + _mm_setr_ps(-M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2), + _mm_setr_ps(-M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2), + _mm_setr_ps(-M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2) + }; + size_t inOffset = 0; + size_t outOffset = 0; + uint8_t tmp = 0; + for (size_t i = 0; i < dataIn->getLength(); i += m_carriers / 4) { + for (size_t j = 0; j < m_carriers / 8; ++j) { + tmp = (in[inOffset] & 0xc0) >> 4; + tmp |= (in[inOffset + (m_carriers / 8)] & 0xc0) >> 6; + out[outOffset] = symbols[tmp]; + tmp = (in[inOffset] & 0x30) >> 2; + tmp |= (in[inOffset + (m_carriers / 8)] & 0x30) >> 4; + out[outOffset + 1] = symbols[tmp]; + tmp = (in[inOffset] & 0x0c); + tmp |= (in[inOffset + (m_carriers / 8)] & 0x0c) >> 2; + out[outOffset + 2] = symbols[tmp]; + tmp = (in[inOffset] & 0x03) << 2; + tmp |= (in[inOffset + (m_carriers / 8)] & 0x03); + out[outOffset + 3] = symbols[tmp]; + ++inOffset; + outOffset += 4; + } + inOffset += m_carriers / 8; + } #else // !__SSE__ - const uint8_t* in = reinterpret_cast<const uint8_t*>(dataIn->getData()); - float* out = reinterpret_cast<float*>(dataOut->getData()); - if (dataIn->getLength() % (d_carriers / 4) != 0) { - throw std::runtime_error( - "QpskSymbolMapper::process input size not valid!"); - } - if (dataOut->getLength() / sizeof(float) != dataIn->getLength() * 4 * 2) { // 4 output complex symbols per input byte - throw std::runtime_error( - "QpskSymbolMapper::process output size not valid!"); - } + const uint8_t* in = reinterpret_cast<const uint8_t*>(dataIn->getData()); + float* out = reinterpret_cast<float*>(dataOut->getData()); + if (dataIn->getLength() % (m_carriers / 4) != 0) { + throw std::runtime_error( + "QpskSymbolMapper::process input size not valid!"); + } + if (dataOut->getLength() / sizeof(float) != dataIn->getLength() * 4 * 2) { // 4 output complex symbols per input byte + throw std::runtime_error( + "QpskSymbolMapper::process output size not valid!"); + } - const static float symbols[16][4] = { - { M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, M_SQRT1_2}, - { M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2}, - { M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2, M_SQRT1_2}, - { M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2}, - { M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2}, - { M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2}, - { M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2}, - { M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2}, - {-M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, M_SQRT1_2}, - {-M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2}, - {-M_SQRT1_2,- M_SQRT1_2, M_SQRT1_2, M_SQRT1_2}, - {-M_SQRT1_2,- M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2}, - {-M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2}, - {-M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2}, - {-M_SQRT1_2,- M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2}, - {-M_SQRT1_2,- M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2} - }; - size_t inOffset = 0; - size_t outOffset = 0; - uint8_t tmp; - for (size_t i = 0; i < dataIn->getLength(); i += d_carriers / 4) { - for (size_t j = 0; j < d_carriers / 8; ++j) { - tmp = (in[inOffset] & 0xc0) >> 4; - tmp |= (in[inOffset + (d_carriers / 8)] & 0xc0) >> 6; - memcpy(&out[outOffset], symbols[tmp], sizeof(float) * 4); - tmp = (in[inOffset] & 0x30) >> 2; - tmp |= (in[inOffset + (d_carriers / 8)] & 0x30) >> 4; - memcpy(&out[outOffset + 4], symbols[tmp], sizeof(float) * 4); - tmp = (in[inOffset] & 0x0c); - tmp |= (in[inOffset + (d_carriers / 8)] & 0x0c) >> 2; - memcpy(&out[outOffset + 8], symbols[tmp], sizeof(float) * 4); - tmp = (in[inOffset] & 0x03) << 2; - tmp |= (in[inOffset + (d_carriers / 8)] & 0x03); - memcpy(&out[outOffset + 12], symbols[tmp], sizeof(float) * 4); - ++inOffset; - outOffset += 4*4; + const static float symbols[16][4] = { + { M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, M_SQRT1_2}, + { M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2}, + { M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2, M_SQRT1_2}, + { M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2}, + { M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2}, + { M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2}, + { M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2}, + { M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2}, + {-M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, M_SQRT1_2}, + {-M_SQRT1_2, M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2}, + {-M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2, M_SQRT1_2}, + {-M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2}, + {-M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2}, + {-M_SQRT1_2, M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2}, + {-M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2, M_SQRT1_2}, + {-M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2, -M_SQRT1_2} + }; + size_t inOffset = 0; + size_t outOffset = 0; + uint8_t tmp; + for (size_t i = 0; i < dataIn->getLength(); i += m_carriers / 4) { + for (size_t j = 0; j < m_carriers / 8; ++j) { + tmp = (in[inOffset] & 0xc0) >> 4; + tmp |= (in[inOffset + (m_carriers / 8)] & 0xc0) >> 6; + memcpy(&out[outOffset], symbols[tmp], sizeof(float) * 4); + tmp = (in[inOffset] & 0x30) >> 2; + tmp |= (in[inOffset + (m_carriers / 8)] & 0x30) >> 4; + memcpy(&out[outOffset + 4], symbols[tmp], sizeof(float) * 4); + tmp = (in[inOffset] & 0x0c); + tmp |= (in[inOffset + (m_carriers / 8)] & 0x0c) >> 2; + memcpy(&out[outOffset + 8], symbols[tmp], sizeof(float) * 4); + tmp = (in[inOffset] & 0x03) << 2; + tmp |= (in[inOffset + (m_carriers / 8)] & 0x03); + memcpy(&out[outOffset + 12], symbols[tmp], sizeof(float) * 4); + ++inOffset; + outOffset += 4*4; + } + inOffset += m_carriers / 8; } - inOffset += d_carriers / 8; - } #endif // __SSE__ + } return 1; } diff --git a/src/QpskSymbolMapper.h b/src/QpskSymbolMapper.h index dbcf4dd..6cf7a2e 100644 --- a/src/QpskSymbolMapper.h +++ b/src/QpskSymbolMapper.h @@ -31,12 +31,13 @@ class QpskSymbolMapper : public ModCodec { public: - QpskSymbolMapper(size_t carriers); + QpskSymbolMapper(size_t carriers, bool fixedPoint); int process(Buffer* const dataIn, Buffer* dataOut); const char* name() { return "QpskSymbolMapper"; } protected: - size_t d_carriers; + bool m_fixedPoint; + size_t m_carriers; }; diff --git a/src/SignalMultiplexer.cpp b/src/SignalMultiplexer.cpp index 1d95bdd..8ecbe78 100644 --- a/src/SignalMultiplexer.cpp +++ b/src/SignalMultiplexer.cpp @@ -28,9 +28,8 @@ #include <string.h> -SignalMultiplexer::SignalMultiplexer(size_t framesize) : - ModMux(), - d_frameSize(framesize) +SignalMultiplexer::SignalMultiplexer() : + ModMux() { PDEBUG("SignalMultiplexer::SignalMultiplexer(%zu) @ %p\n", framesize, this); diff --git a/src/SignalMultiplexer.h b/src/SignalMultiplexer.h index 5186a8d..1f6bc12 100644 --- a/src/SignalMultiplexer.h +++ b/src/SignalMultiplexer.h @@ -36,7 +36,7 @@ class SignalMultiplexer : public ModMux { public: - SignalMultiplexer(size_t frameSize); + SignalMultiplexer(); virtual ~SignalMultiplexer(); SignalMultiplexer(const SignalMultiplexer&); SignalMultiplexer& operator=(const SignalMultiplexer&); @@ -44,8 +44,5 @@ public: int process(std::vector<Buffer*> dataIn, Buffer* dataOut); const char* name() { return "SignalMultiplexer"; } - -protected: - size_t d_frameSize; }; diff --git a/src/Utils.cpp b/src/Utils.cpp index 788d125..0065bc1 100644 --- a/src/Utils.cpp +++ b/src/Utils.cpp @@ -26,9 +26,9 @@ */ #include <ctime> +#include <cstring> #include <sstream> #include "Utils.h" -#include "GainControl.h" #if defined(HAVE_PRCTL) # include <sys/prctl.h> #endif |