Skip to content
Permalink

Comparing changes

Choose two branches to see what’s changed or to start a new pull request. If you need to, you can also or learn more about diff comparisons.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also . Learn more about diff comparisons here.
base repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected base ref is valid, then try again.
Loading
base: 8248da04898e
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 9bb331004a00
Choose a head ref
  • 16 commits
  • 13 files changed
  • 1 contributor

Commits on Jul 8, 2020

  1. A skeleton.

    pleroy committed Jul 8, 2020
    Copy the full SHA
    57d80f5 View commit details

Commits on Jul 9, 2020

  1. Some implementation code.

    pleroy committed Jul 9, 2020
    Copy the full SHA
    c0dd007 View commit details
  2. Copy the full SHA
    d11f438 View commit details

Commits on Jul 11, 2020

  1. Bit reversal.

    pleroy committed Jul 11, 2020
    Copy the full SHA
    a596d75 View commit details
  2. Danielson-Lánczos algorithm.

    pleroy committed Jul 11, 2020
    Copy the full SHA
    2fc732a View commit details
  3. A test and compilation errors.

    pleroy committed Jul 11, 2020
    Copy the full SHA
    5648468 View commit details
  4. Matchers and tests.

    pleroy committed Jul 11, 2020
    Copy the full SHA
    e9ff0fe View commit details
  5. Cleanup.

    pleroy committed Jul 11, 2020
    Copy the full SHA
    88b6dd4 View commit details
  6. Short FFTs specializations.

    pleroy committed Jul 11, 2020
    Copy the full SHA
    e734c4d View commit details
  7. Cleanup.

    pleroy committed Jul 11, 2020
    Copy the full SHA
    69af03b View commit details
  8. Lint.

    pleroy committed Jul 11, 2020
    Copy the full SHA
    f8dbf7c View commit details

Commits on Jul 12, 2020

  1. After egg's review.

    pleroy committed Jul 12, 2020
    Copy the full SHA
    d8d6c5c View commit details
  2. Formatting.

    pleroy committed Jul 12, 2020
    Copy the full SHA
    9162bec View commit details
  3. Cleanup.

    pleroy committed Jul 12, 2020
    Copy the full SHA
    6840c62 View commit details

Commits on Jul 13, 2020

  1. After egg's review.

    pleroy committed Jul 13, 2020
    Copy the full SHA
    33b2bbe View commit details
  2. Merge pull request #2637 from pleroy/FFT

    A rough implementation of FFT
    pleroy authored Jul 13, 2020
    Copy the full SHA
    9bb3310 View commit details
5 changes: 5 additions & 0 deletions base/bits.hpp
Original file line number Diff line number Diff line change
@@ -10,6 +10,11 @@ constexpr int FloorLog2(int n);
// Greatest power of 2 less than or equal to n. 8 ↦ 8, 7 ↦ 4.
constexpr int PowerOf2Le(int n);

// Computes bitreversed(bitreversed(n) + 1) assuming that n is represented on
// the given number of bits. For 4 bits:
// 0 ↦ 8 ↦ 4 ↦ C ↦ 2 ↦ A ↦ 6 ↦ E ↦ 1 ↦ 9 ↦ 5 ↦ D ↦ 3 ↦ B ↦ 7 ↦ F
constexpr int BitReversedIncrement(int n, int bits);

} // namespace base
} // namespace principia

21 changes: 21 additions & 0 deletions base/bits_body.hpp
Original file line number Diff line number Diff line change
@@ -2,6 +2,11 @@

#include "base/bits.hpp"

#include <cstdint>

#include "base/macros.hpp"
#include "glog/logging.h"

namespace principia {
namespace base {

@@ -13,5 +18,21 @@ constexpr int PowerOf2Le(int const n) {
return n == 0 ? 0 : n == 1 ? 1 : PowerOf2Le(n >> 1) << 1;
}

constexpr int BitReversedIncrement(int const n, int const bits) {
CONSTEXPR_DCHECK(n >= 0 && n < 1 << bits);
CONSTEXPR_DCHECK(bits > 0 && bits < 32);
// [War03], chapter 7.1 page 105.
std::uint32_t mask = 0x8000'0000;
std::uint32_t x = n << (32 - bits);
x ^= mask;
if (static_cast<std::int32_t>(x) >= 0) {
do {
mask >>= 1;
x ^= mask;
} while (x < mask);
}
return x >> (32 - bits);
}

} // namespace base
} // namespace principia
18 changes: 18 additions & 0 deletions base/bits_test.cpp
Original file line number Diff line number Diff line change
@@ -25,5 +25,23 @@ TEST(BitsTest, PowerOf2Le) {
EXPECT_EQ(8, PowerOf2Le(8));
}

TEST(BitsTest, BitReversedIncrement) {
EXPECT_EQ(0x8, BitReversedIncrement(0x0, 4));
EXPECT_EQ(0x9, BitReversedIncrement(0x1, 4));
EXPECT_EQ(0xA, BitReversedIncrement(0x2, 4));
EXPECT_EQ(0xB, BitReversedIncrement(0x3, 4));
EXPECT_EQ(0xC, BitReversedIncrement(0x4, 4));
EXPECT_EQ(0xD, BitReversedIncrement(0x5, 4));
EXPECT_EQ(0xE, BitReversedIncrement(0x6, 4));
EXPECT_EQ(0xF, BitReversedIncrement(0x7, 4));
EXPECT_EQ(0x4, BitReversedIncrement(0x8, 4));
EXPECT_EQ(0x5, BitReversedIncrement(0x9, 4));
EXPECT_EQ(0x6, BitReversedIncrement(0xA, 4));
EXPECT_EQ(0x7, BitReversedIncrement(0xB, 4));
EXPECT_EQ(0x2, BitReversedIncrement(0xC, 4));
EXPECT_EQ(0x3, BitReversedIncrement(0xD, 4));
EXPECT_EQ(0x1, BitReversedIncrement(0xE, 4));
}

} // namespace base
} // namespace principia
1 change: 1 addition & 0 deletions base/macros.hpp
Original file line number Diff line number Diff line change
@@ -153,6 +153,7 @@ inline void noreturn() { std::exit(0); }
// Needed to circumvent lint warnings in constexpr functions where CHECK_LT and
// friends cannot be used.
#define CONSTEXPR_CHECK(condition) CHECK(condition)
#define CONSTEXPR_DCHECK(condition) DCHECK(condition)

// Clang for some reason doesn't like FP arithmetic that yields infinities in
// constexpr code (MSVC and GCC are fine with that). This will be fixed in
20 changes: 20 additions & 0 deletions documentation/bibliography.bib
Original file line number Diff line number Diff line change
@@ -247,6 +247,18 @@ @article{ChinKidwell2000
volume = {62},
}

@article{DanielsonLánczos1942,
author = {Danielson, G. C. and Lánczos, C.},
url = {https://www.sciencedirect.com/science/article/abs/pii/S0016003242907671},
date = {1942},
doi = {10.1016/S0016-0032(42)90767-1},
journaltitle = {Journal of the Franklin Institute},
number = {4},
pages = {365--380},
title = {Some improvements in practical Fourier analysis and their application to x-ray scattering from liquids},
volume = {233},
}

@article{DieleMarangi2011,
author = {Diele, F. and Marangi, C.},
url = {http://www.sciencedirect.com/science/article/pii/S0168927411000353},
@@ -594,6 +606,14 @@ @article{Montenbruck1992
volume = {53},
}

@article{Myrnyy2007,
author = {Myrnyy, Vlodymyr},
date = {2007},
eprint = {https://www.drdobbs.com/cpp/a-simple-and-efficient-fft-implementatio/199500857},
journaltitle = {Dr. Dobb's},
title = {A Simple and Efficient FFT Implementation in C++},
}

@article{Newhall1989,
author = {Newhall, X. X.},
date = {1989},
Binary file modified documentation/bibliography.pdf
Binary file not shown.
56 changes: 56 additions & 0 deletions numerics/fast_fourier_transform.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,56 @@

#pragma once

#include <array>
#include <complex>
#include <type_traits>
#include <vector>

#include "base/bits.hpp"
#include "quantities/named_quantities.hpp"

namespace principia {
namespace numerics {
namespace internal_fast_fourier_transform {

using base::FloorLog2;

// This class computes Fourier[{...}, FourierParameters -> {1, -1}] in
// Mathematica notation (the "signal processing" Fourier transform).
template<typename Scalar, std::size_t size_>
class FastFourierTransform {
public:
// The size must be a power of 2.
static constexpr int size = size_;
static constexpr int log2_size = FloorLog2(size);
static_assert(size == 1 << log2_size);

// In the constructors, the container must have |size| elements.

template<typename Container,
typename = std::enable_if_t<
std::is_convertible_v<typename Container::value_type, Scalar>>>
explicit FastFourierTransform(Container const& container);

template<typename Iterator,
typename = std::enable_if_t<std::is_convertible_v<
typename std::iterator_traits<Iterator>::value_type,
Scalar>>>
FastFourierTransform(Iterator begin, Iterator end);

explicit FastFourierTransform(std::array<Scalar, size> const& container);

private:
std::array<std::complex<double>, size> transform_;

friend class FastFourierTransformTest;
};

} // namespace internal_fast_fourier_transform

using internal_fast_fourier_transform::FastFourierTransform;

} // namespace numerics
} // namespace principia

#include "numerics/fast_fourier_transform_body.hpp"
151 changes: 151 additions & 0 deletions numerics/fast_fourier_transform_body.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,151 @@
#pragma once

#include "numerics/fast_fourier_transform.hpp"

#include "base/bits.hpp"
#include "quantities/elementary_functions.hpp"
#include "quantities/numbers.hpp"
#include "quantities/quantities.hpp"
#include "quantities/si.hpp"

namespace principia {
namespace numerics {
namespace internal_fast_fourier_transform {

using base::BitReversedIncrement;
using quantities::Angle;
using quantities::Sin;
using quantities::si::Radian;
using quantities::si::Unit;
namespace si = quantities::si;

// Implementation of the Danielson-Lánczos algorithm using templates for
// recursion and template specializations for short FFTs [DL42, Myr07].
template<int array_size_, int chunk_size_ = array_size_>
class DanielsonLánczos {
public:
static void Transform(typename std::array<std::complex<double>,
array_size_>::iterator const begin);
};

template<int array_size_>
class DanielsonLánczos<array_size_, 1> {
public:
static void Transform(typename std::array<std::complex<double>,
array_size_>::iterator const begin);
};

template<int array_size_>
class DanielsonLánczos<array_size_, 2> {
public:
static void Transform(typename std::array<std::complex<double>,
array_size_>::iterator const begin);
};

template<int array_size_>
class DanielsonLánczos<array_size_, 4> {
public:
static void Transform(typename std::array<std::complex<double>,
array_size_>::iterator const begin);
};

template<int array_size_, int chunk_size_>
void DanielsonLánczos<array_size_, chunk_size_>::Transform(
typename std::array<std::complex<double>, array_size_>::iterator const
begin) {
constexpr int N = chunk_size_;

DanielsonLánczos<array_size_, N / 2>::Transform(begin);
DanielsonLánczos<array_size_, N / 2>::Transform(begin + N / 2);

Angle const θ = π * Radian / N;
double const sin_θ = Sin(θ);
double const cos_2θ_minus_1 = -2 * sin_θ * sin_θ;
double const sin_2θ = Sin(2 * θ);
// Computing e⁻²ⁱ⁽ᵏ⁺¹⁾ᶿ as e⁻²ⁱᵏᶿ + e⁻²ⁱᵏᶿ (e⁻²ⁱᶿ - 1) rather than e⁻²ⁱᵏᶿe⁻²ⁱᶿ
// improves accuracy [Myr07].
std::complex<double> const e⁻²ⁱᶿ_minus_1(cos_2θ_minus_1, -sin_2θ);
std::complex<double> e⁻²ⁱᵏᶿ = 1;
auto it = begin;
for (int k = 0; k < N / 2; ++it, ++k, e⁻²ⁱᵏᶿ += e⁻²ⁱᵏᶿ * (e⁻²ⁱᶿ_minus_1)) {
auto const t = *(it + N / 2) * e⁻²ⁱᵏᶿ;
*(it + N / 2) = *it - t;
*it += t;
}
}

template<int array_size_>
void DanielsonLánczos<array_size_, 1>::Transform(
typename std::array<std::complex<double>, array_size_>::iterator const
begin) {}

template<int array_size_>
void DanielsonLánczos<array_size_, 2>::Transform(
typename std::array<std::complex<double>, array_size_>::iterator const
begin) {
auto const t = *(begin + 1);
*(begin + 1) = *begin - t;
*begin += t;
}

template<int array_size_>
void DanielsonLánczos<array_size_, 4>::Transform(
typename std::array<std::complex<double>, array_size_>::iterator const
begin) {
{
auto const t = *(begin + 1);
*(begin + 1) = *begin - t;
*begin += t;
}
{
auto const t = *(begin + 3);
*(begin + 3) = {(begin + 2)->imag() - t.imag(),
t.real() - (begin + 2)->real()};
*(begin + 2) += t;
}
{
auto const t = *(begin + 2);
*(begin + 2) = *begin - t;
*begin += t;
}
{
auto const t = *(begin + 3);
*(begin + 3) = *(begin + 1) - t;
*(begin + 1) += t;
}
}

template<typename Scalar, std::size_t size_>
template<typename Container, typename>
FastFourierTransform<Scalar, size_>::FastFourierTransform(
Container const& container)
: FastFourierTransform(container.cbegin(), container.cend()) {}

template<typename Scalar, std::size_t size_>
template<typename Iterator, typename>
FastFourierTransform<Scalar, size_>::FastFourierTransform(
Iterator const begin,
Iterator const end) {
DCHECK_EQ(size, std::distance(begin, end));

// Type decay, reindexing, and promotion to complex.
int bit_reversed_index = 0;
for (auto it = begin;
it != end;
++it,
bit_reversed_index = BitReversedIncrement(bit_reversed_index,
log2_size)) {
transform_[bit_reversed_index] = *it / si::Unit<Scalar>;
}

DanielsonLánczos<size>::Transform(transform_.begin());
}

template<typename Scalar, std::size_t size_>
FastFourierTransform<Scalar, size_>::FastFourierTransform(
std::array<Scalar, size> const& container)
: FastFourierTransform(container.cbegin(), container.cend()) {}

} // namespace internal_fast_fourier_transform
} // namespace numerics
} // namespace principia
Loading