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: 9bb331004a00
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 8f3747110ab3
Choose a head ref
  • 8 commits
  • 3 files changed
  • 1 contributor

Commits on Jul 12, 2020

  1. Power spectrum.

    pleroy committed Jul 12, 2020
    Copy the full SHA
    94b61a8 View commit details
  2. Mode analysis.

    pleroy committed Jul 12, 2020
    Copy the full SHA
    366183a View commit details

Commits on Jul 13, 2020

  1. Merge.

    pleroy committed Jul 13, 2020
    Copy the full SHA
    d688177 View commit details
  2. The merge was a mess.

    pleroy committed Jul 13, 2020
    Copy the full SHA
    e4157b3 View commit details
  3. Tolerances.

    pleroy committed Jul 13, 2020
    Copy the full SHA
    f3809d6 View commit details
  4. Lint.

    pleroy committed Jul 13, 2020
    Copy the full SHA
    8e78e7e View commit details
  5. After egg's review.

    pleroy committed Jul 13, 2020
    Copy the full SHA
    0adcfe2 View commit details
  6. Merge pull request #2646 from pleroy/Spectrum

    Spectrum extraction using FFT
    pleroy authored Jul 13, 2020
    Copy the full SHA
    8f37471 View commit details
Showing with 145 additions and 15 deletions.
  1. +25 −5 numerics/fast_fourier_transform.hpp
  2. +51 −5 numerics/fast_fourier_transform_body.hpp
  3. +69 −5 numerics/fast_fourier_transform_test.cpp
30 changes: 25 additions & 5 deletions numerics/fast_fourier_transform.hpp
Original file line number Diff line number Diff line change
@@ -1,19 +1,25 @@


#pragma once

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

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

namespace principia {
namespace numerics {
namespace internal_fast_fourier_transform {

using base::FloorLog2;
using geometry::Interval;
using quantities::AngularFrequency;
using quantities::Square;
using quantities::Time;

// This class computes Fourier[{...}, FourierParameters -> {1, -1}] in
// Mathematica notation (the "signal processing" Fourier transform).
@@ -25,22 +31,36 @@ class FastFourierTransform {
static constexpr int log2_size = FloorLog2(size);
static_assert(size == 1 << log2_size);

// In the constructors, the container must have |size| elements.
// In the constructors, the container must have |size| elements. The samples
// are assumed to be separated by Δt.

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

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);
FastFourierTransform(Iterator begin, Iterator end,
Time const& Δt);

FastFourierTransform(std::array<Scalar, size> const& container,
Time const& Δt);

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

// Returns the interval that contains the largest peak of power.
Interval<AngularFrequency> Mode() const;

private:
Time const Δt_;
AngularFrequency const ω_;

// The elements of transform_ are in SI units of Scalar. They are spaced in
// frequency by ω_.
std::array<std::complex<double>, size> transform_;

friend class FastFourierTransformTest;
56 changes: 51 additions & 5 deletions numerics/fast_fourier_transform_body.hpp
Original file line number Diff line number Diff line change
@@ -2,6 +2,9 @@

#include "numerics/fast_fourier_transform.hpp"

#include <map>
#include <optional>

#include "base/bits.hpp"
#include "quantities/elementary_functions.hpp"
#include "quantities/numbers.hpp"
@@ -118,14 +121,18 @@ void DanielsonLánczos<array_size_, 4>::Transform(
template<typename Scalar, std::size_t size_>
template<typename Container, typename>
FastFourierTransform<Scalar, size_>::FastFourierTransform(
Container const& container)
: FastFourierTransform(container.cbegin(), container.cend()) {}
Container const& container,
Time const& Δt)
: FastFourierTransform(container.cbegin(), container.cend(), Δt) {}

template<typename Scalar, std::size_t size_>
template<typename Iterator, typename>
FastFourierTransform<Scalar, size_>::FastFourierTransform(
Iterator const begin,
Iterator const end) {
Iterator const end,
Time const& Δt)
: Δt_(Δt),
ω_(2 * π * Radian / (size * Δt_)) {
DCHECK_EQ(size, std::distance(begin, end));

// Type decay, reindexing, and promotion to complex.
@@ -143,8 +150,47 @@ FastFourierTransform<Scalar, size_>::FastFourierTransform(

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

template<typename Scalar, std::size_t size_>
std::map<AngularFrequency, Square<Scalar>>
FastFourierTransform<Scalar, size_>::PowerSpectrum() const {
std::map<AngularFrequency, Square<Scalar>> spectrum;
int k = 0;
for (auto const& coefficient : transform_) {
spectrum.emplace_hint(spectrum.end(),
k * ω_,
std::norm(coefficient) * si::Unit<Square<Scalar>>);
++k;
}
return spectrum;
}

template<typename Scalar, std::size_t size_>
Interval<AngularFrequency> FastFourierTransform<Scalar, size_>::Mode() const {
auto const spectrum = PowerSpectrum();
typename std::map<AngularFrequency, Square<Scalar>>::const_iterator max =
spectrum.end();

// Only look at the first size / 2 + 1 elements because the spectrum is
// symmetrical.
auto it = spectrum.begin();
for (int i = 0; i < size / 2 + 1; ++i, ++it) {
if (max == spectrum.end() || it->second > max->second) {
max = it;
}
}
Interval<AngularFrequency> result;
if (max == spectrum.begin()) {
result.Include(max->first);
} else {
result.Include(std::prev(max)->first);
}
result.Include(std::next(max)->first);
return result;
}

} // namespace internal_fast_fourier_transform
} // namespace numerics
74 changes: 69 additions & 5 deletions numerics/fast_fourier_transform_test.cpp
Original file line number Diff line number Diff line change
@@ -1,20 +1,27 @@
#include "numerics/fast_fourier_transform.hpp"
#include "numerics/fast_fourier_transform.hpp"

#include <algorithm>
#include <random>
#include <vector>

#include "gmock/gmock.h"
#include "gtest/gtest.h"
#include "quantities/elementary_functions.hpp"
#include "quantities/si.hpp"
#include "testing_utilities/almost_equals.hpp"

namespace principia {
namespace numerics {
namespace internal_fast_fourier_transform {

using quantities::Length;
using quantities::Sqrt;
using quantities::si::Metre;
using quantities::si::Second;
using testing_utilities::AlmostEquals;
using ::testing::ElementsAre;
using ::testing::Lt;
using ::testing::Pair;

class FastFourierTransformTest : public ::testing::Test {
protected:
@@ -28,10 +35,11 @@ class FastFourierTransformTest : public ::testing::Test {
};

TEST_F(FastFourierTransformTest, Square) {
FastFourierTransform<double, 8> const transform({1, 1, 1, 1, 0, 0, 0, 0});
FastFourierTransform<double, 8> const transform({1, 1, 1, 1, 0, 0, 0, 0},
1 * Second);
EXPECT_THAT(Coefficients(transform),
ElementsAre(AlmostEquals(Complex{4}, 0),
AlmostEquals(Complex{1, -1 - Sqrt(2)}, 1),
AlmostEquals(Complex{1, -1 - Sqrt(2)}, 0, 1),
AlmostEquals(Complex{0}, 0),
AlmostEquals(Complex{1, 1 - Sqrt(2)}, 4),
AlmostEquals(Complex{0}, 0),
@@ -57,7 +65,8 @@ TEST_F(FastFourierTransformTest, Sin) {
-0.63126663787232131146,
-0.21483085764466499644,
+0.24754738092257664739,
+0.65698659871878909040});
+0.65698659871878909040},
1 * Second);
EXPECT_THAT(
Coefficients(transform),
ElementsAre(
@@ -80,7 +89,7 @@ TEST_F(FastFourierTransformTest, Sin) {
AlmostEquals(Complex{-0.5400094598886346723188351187,
-0.0975085343055921838098913009}, 3),
AlmostEquals(Complex{-0.5504467338298529673542814016,
-0.2045798116680680999699630079}, 2),
-0.2045798116680680999699630079}, 2, 6),
AlmostEquals(Complex{-0.5726936754458237681916166204,
-0.3352695668656918501471058013}, 6),
AlmostEquals(Complex{-0.6197133589762012113067702090,
@@ -91,6 +100,61 @@ TEST_F(FastFourierTransformTest, Sin) {
-1.7603441806320655123815748876}, 2),
AlmostEquals(Complex{+4.0811719972720017737297705383,
+5.7509914354728044020475598497}, 1)));

EXPECT_THAT(
transform.PowerSpectrum(),
ElementsAre(Pair(AlmostEquals(0 * Radian / Second, 0),
AlmostEquals(0.7161225188062225589818894003, 0)),
Pair(AlmostEquals(π / 8 * Radian / Second, 0),
AlmostEquals(49.729867362198687411669694816, 0)),
Pair(AlmostEquals(π / 4 * Radian / Second, 0),
AlmostEquals(4.5768125922478623650817219388, 1)),
Pair(AlmostEquals(3 * π / 8 * Radian / Second, 0),
AlmostEquals(1.2458226823327073992355624995, 3)),
Pair(AlmostEquals(π / 2 * Radian / Second, 0),
AlmostEquals(0.6527765348939019854655626516, 1)),
Pair(AlmostEquals(5 * π / 8 * Radian / Second, 0),
AlmostEquals(0.4403837283619551481413056610, 2)),
Pair(AlmostEquals(3 * π / 4 * Radian / Second, 0),
AlmostEquals(0.3448445061260952118899789115, 1)),
Pair(AlmostEquals(7 * π / 8 * Radian / Second, 0),
AlmostEquals(0.3011181310316397868684530905, 18)),
Pair(AlmostEquals(π * Radian / Second, 0),
AlmostEquals(0.2882738863361058320489546747, 4)),
Pair(AlmostEquals(9 * π / 8 * Radian / Second, 0),
AlmostEquals(0.3011181310316397868684530905, 0)),
Pair(AlmostEquals(5 * π / 4 * Radian / Second, 0),
AlmostEquals(0.3448445061260952118899789115, 4)),
Pair(AlmostEquals(11 * π / 8 * Radian / Second, 0),
AlmostEquals(0.4403837283619551481413056610, 15)),
Pair(AlmostEquals(3 * π / 2 * Radian / Second, 0),
AlmostEquals(0.6527765348939019854655626516, 1)),
Pair(AlmostEquals(13 * π / 8 * Radian / Second, 0),
AlmostEquals(1.2458226823327073992355624995, 1)),
Pair(AlmostEquals(7 * π / 4 * Radian / Second, 0),
AlmostEquals(4.5768125922478623650817219388, 2)),
Pair(AlmostEquals(15 * π / 8 * Radian / Second, 0),
AlmostEquals(49.729867362198687411669694816, 0))));
}

TEST_F(FastFourierTransformTest, Mode) {
using FFT = FastFourierTransform<Length, 1 << 16>;
AngularFrequency const ω = 666 * π / FFT::size * Radian / Second;
Time const Δt = 1 * Second;
std::mt19937_64 random(42);
std::uniform_real_distribution<> noise(-0.5, 0.5);
std::vector<Length> signal;
for (int n = 0; n < FFT::size; ++n) {
signal.push_back((Sin(n * ω * Δt) + noise(random)) * Metre);
}

// Won't fit on the stack.
auto transform = std::make_unique<FFT>(signal, Δt);

auto const mode = transform->Mode();
EXPECT_THAT(mode.midpoint(), AlmostEquals(ω, 0));
EXPECT_THAT(mode.measure(),
AlmostEquals(4 * π / FFT::size * Radian / Second, 24));
}

} // namespace internal_fast_fourier_transform