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

Commits on Oct 6, 2020

  1. Verified

    This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
    Copy the full SHA
    dfa859d View commit details
  2. fix a test

    eggrobin committed Oct 6, 2020
    Copy the full SHA
    9ff61bb View commit details

Commits on Oct 7, 2020

  1. After pleroy’s review

    eggrobin committed Oct 7, 2020
    Copy the full SHA
    ca72454 View commit details
  2. Merge pull request #2752 from eggrobin/no-time

    Allow arguments other than Time in FFT
    eggrobin authored Oct 7, 2020
    Copy the full SHA
    a127da7 View commit details
39 changes: 27 additions & 12 deletions numerics/fast_fourier_transform.hpp
Original file line number Diff line number Diff line change
@@ -21,46 +21,61 @@ using base::FloorLog2;
using geometry::Complexification;
using geometry::Hilbert;
using geometry::Interval;
using quantities::AngularFrequency;
using quantities::Time;

// This class computes Fourier[{...}, FourierParameters -> {1, -1}] in
// Mathematica notation (the "signal processing" Fourier transform).
template<typename Value, std::size_t size_>
using quantities::Angle;
using quantities::Derivative;
using quantities::Difference;

// Given (u₀, ..., uₙ₋₁), this class computes the discrete Fourier transform
// Uₛ = ∑ᵣ uᵣ exp(-2πirs/n),
// corresponding to
// Fourier[{...}, FourierParameters -> {1, -1}]
// in Mathematica notation (the "signal processing" Fourier transform).
template<typename Value, typename Argument, std::size_t size_>
class FastFourierTransform {
public:
// This is only an actual angular frequency if |Argument| is time-like.
// If |Argument| is an angular frequency, this is a time.
using AngularFrequency = Derivative<Angle, Argument>;

// 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. The samples
// are assumed to be separated by Δt.
// In the constructors, the container must have |size| elements. For the
// purpose of expressing the frequencies, the values are assumed to be
// sampled at intervals of Δt.

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

template<typename Iterator,
typename = std::enable_if_t<std::is_convertible_v<
typename std::iterator_traits<Iterator>::value_type,
Value>>>
FastFourierTransform(Iterator begin, Iterator end,
Time const& Δt);
Difference<Argument> const& Δt);

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

std::map<AngularFrequency, typename Hilbert<Value>::InnerProductType>
PowerSpectrum() const;

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

// Given s ∈ [0, size - 1] ∩ ℕ, returns the coefficient Uₛ.
Complexification<Value> const& operator[](int s) const;

// Given s ∈ [0, size - 1] ∩ ℕ, returns the frequency corresponding to Uₛ.
AngularFrequency frequency(int s) const;

private:
Time const Δt_;
Difference<Argument> const Δt_;
AngularFrequency const Δω_;

// The elements of transform_ are spaced in frequency by ω_.
44 changes: 30 additions & 14 deletions numerics/fast_fourier_transform_body.hpp
Original file line number Diff line number Diff line change
@@ -112,19 +112,19 @@ void DanielsonLánczos<Complex, array_size_, 4>::Transform(
}
}

template<typename Value, std::size_t size_>
template<typename Value, typename Argument, std::size_t size_>
template<typename Container, typename>
FastFourierTransform<Value, size_>::FastFourierTransform(
FastFourierTransform<Value, Argument, size_>::FastFourierTransform(
Container const& container,
Time const& Δt)
Difference<Argument> const& Δt)
: FastFourierTransform(container.cbegin(), container.cend(), Δt) {}

template<typename Value, std::size_t size_>
template<typename Value, typename Argument, std::size_t size_>
template<typename Iterator, typename>
FastFourierTransform<Value, size_>::FastFourierTransform(
FastFourierTransform<Value, Argument, size_>::FastFourierTransform(
Iterator const begin,
Iterator const end,
Time const& Δt)
Difference<Argument> const& Δt)
: Δt_(Δt),
Δω_(2 * π * Radian / (size * Δt_)) {
DCHECK_EQ(size, std::distance(begin, end));
@@ -143,15 +143,15 @@ FastFourierTransform<Value, size_>::FastFourierTransform(
transform_.begin());
}

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

template<typename Value, std::size_t size_>
std::map<AngularFrequency, typename Hilbert<Value>::InnerProductType>
FastFourierTransform<Value, size_>::PowerSpectrum() const {
template<typename Value, typename Argument, std::size_t size_>
auto FastFourierTransform<Value, Argument, size_>::PowerSpectrum() const
-> std::map<AngularFrequency, typename Hilbert<Value>::InnerProductType> {
std::map<AngularFrequency, typename Hilbert<Value>::InnerProductType>
spectrum;
int k = 0;
@@ -162,8 +162,9 @@ FastFourierTransform<Value, size_>::PowerSpectrum() const {
return spectrum;
}

template<typename Value, std::size_t size_>
Interval<AngularFrequency> FastFourierTransform<Value, size_>::Mode() const {
template<typename Value, typename Argument, std::size_t size_>
auto FastFourierTransform<Value, Argument, size_>::Mode() const
-> Interval<AngularFrequency> {
auto const spectrum = PowerSpectrum();
typename std::map<AngularFrequency,
typename Hilbert<Value>::InnerProductType>::const_iterator
@@ -187,6 +188,21 @@ Interval<AngularFrequency> FastFourierTransform<Value, size_>::Mode() const {
return result;
}

template<typename Value, typename Argument, std::size_t size_>
Complexification<Value> const&
FastFourierTransform<Value, Argument, size_>::operator[](
int const s) const {
return transform_[s];
}

template<typename Value, typename Argument, std::size_t size_>
typename FastFourierTransform<Value, Argument, size_>::AngularFrequency
FastFourierTransform<Value, Argument, size_>::frequency(int const s) const {
DCHECK_GE(s, 0);
DCHECK_LT(s, size);
return s * Δω_;
}

} // namespace internal_fast_fourier_transform
} // namespace numerics
} // namespace principia
68 changes: 46 additions & 22 deletions numerics/fast_fourier_transform_test.cpp
Original file line number Diff line number Diff line change
@@ -21,13 +21,19 @@ using geometry::Displacement;
using geometry::Frame;
using geometry::Handedness;
using geometry::Inertial;
using geometry::Instant;
using quantities::AngularFrequency;
using quantities::Cos;
using quantities::Length;
using quantities::Sqrt;
using quantities::Time;
using quantities::Voltage;
using quantities::si::Metre;
using quantities::si::Second;
using quantities::si::Volt;
using testing_utilities::AlmostEquals;
using ::testing::ElementsAre;
using ::testing::ElementsAreArray;
using ::testing::Lt;
using ::testing::Pair;

@@ -42,14 +48,14 @@ class FastFourierTransformTest : public ::testing::Test {

template<typename Scalar, std::size_t size_>
std::array<Complexification<double>, size_> Coefficients(
FastFourierTransform<Scalar, size_> const& fft) {
FastFourierTransform<Scalar, Instant, size_> const& fft) {
return fft.transform_;
}
};

TEST_F(FastFourierTransformTest, Square) {
FastFourierTransform<double, 8> const transform({1, 1, 1, 1, 0, 0, 0, 0},
1 * Second);
FastFourierTransform<double, Instant, 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)}, 0, 1),
@@ -63,23 +69,24 @@ TEST_F(FastFourierTransformTest, Square) {

TEST_F(FastFourierTransformTest, Sin) {
// Sin(x) on [0, 7].
FastFourierTransform<double, 16> const transform({+0,
+0.44991188055599964373,
+0.80360826369441117592,
+0.98544972998846018066,
+0.95654873748436662401,
+0.72308588173832461680,
+0.33498815015590491954,
-0.12474816864589884767,
-0.55780658091328209620,
-0.87157577241358806002,
-0.99895491709792831520,
-0.91270346343588987220,
-0.63126663787232131146,
-0.21483085764466499644,
+0.24754738092257664739,
+0.65698659871878909040},
1 * Second);
FastFourierTransform<double, Instant, 16> const transform(
{+0,
+0.44991188055599964373,
+0.80360826369441117592,
+0.98544972998846018066,
+0.95654873748436662401,
+0.72308588173832461680,
+0.33498815015590491954,
-0.12474816864589884767,
-0.55780658091328209620,
-0.87157577241358806002,
-0.99895491709792831520,
-0.91270346343588987220,
-0.63126663787232131146,
-0.21483085764466499644,
+0.24754738092257664739,
+0.65698659871878909040},
1 * Second);
EXPECT_THAT(
Coefficients(transform),
ElementsAre(
@@ -151,7 +158,7 @@ TEST_F(FastFourierTransformTest, Sin) {
}

TEST_F(FastFourierTransformTest, Mode) {
using FFT = FastFourierTransform<Length, 1 << 16>;
using FFT = FastFourierTransform<Length, Instant, 1 << 16>;
AngularFrequency const ω = 666 * π / FFT::size * Radian / Second;
Time const Δt = 1 * Second;
std::mt19937_64 random(42);
@@ -171,7 +178,7 @@ TEST_F(FastFourierTransformTest, Mode) {
}

TEST_F(FastFourierTransformTest, Vector) {
using FFT = FastFourierTransform<Displacement<World>, 1 << 16>;
using FFT = FastFourierTransform<Displacement<World>, Instant, 1 << 16>;
AngularFrequency const ω = 666 * π / FFT::size * Radian / Second;
Time const Δt = 1 * Second;
std::vector<Displacement<World>> signal;
@@ -190,6 +197,23 @@ TEST_F(FastFourierTransformTest, Vector) {
AlmostEquals(4 * π / FFT::size * Radian / Second, 24));
}

TEST_F(FastFourierTransformTest, Inverse) {
constexpr int n = 4;
constexpr Time Δt = 1 * Second;
std::array<Voltage, n> v{{1 * Volt, 0 * Volt, 1 * Volt, 0 * Volt}};
FastFourierTransform<Voltage, Instant, n> const V(v, Δt);

FastFourierTransform<Voltage, AngularFrequency, n> const nv(
{V[0].real_part(), V[1].real_part(), V[2].real_part(), V[3].real_part()},
V.frequency(1) - V.frequency(0));
EXPECT_THAT((std::array{nv[0].real_part() / n,
nv[1].real_part() / n,
nv[2].real_part() / n,
nv[3].real_part() / n}),
ElementsAreArray(v));
EXPECT_THAT(nv.frequency(1) - nv.frequency(0), AlmostEquals(Δt, 0));
}

} // namespace internal_fast_fourier_transform
} // namespace numerics
} // namespace principia
4 changes: 2 additions & 2 deletions numerics/frequency_analysis_test.cpp
Original file line number Diff line number Diff line change
@@ -106,7 +106,7 @@ class FrequencyAnalysisTest : public ::testing::Test {
};

TEST_F(FrequencyAnalysisTest, PreciseModeScalar) {
using FFT = FastFourierTransform<Length, 1 << 16>;
using FFT = FastFourierTransform<Length, Instant, 1 << 16>;
AngularFrequency const ω = 666.543 * π / FFT::size * Radian / Second;
Time const Δt = 1 * Second;
std::mt19937_64 random(42);
@@ -164,7 +164,7 @@ TEST_F(FrequencyAnalysisTest, PreciseModeScalar) {
}

TEST_F(FrequencyAnalysisTest, PreciseModeVector) {
using FFT = FastFourierTransform<Displacement<World>, 1 << 16>;
using FFT = FastFourierTransform<Displacement<World>, Instant, 1 << 16>;
AngularFrequency const ω = 666.543 * π / FFT::size * Radian / Second;
Time const Δt = 1 * Second;

1 change: 1 addition & 0 deletions physics/analytical_series_test.cpp
Original file line number Diff line number Diff line change
@@ -94,6 +94,7 @@ class AnalyticalSeriesTest : public ::testing::Test {
LOG(INFO) << "max_residual=" << max_residual;
auto fft =
std::make_unique<FastFourierTransform<Displacement<ICRS>,
Instant,
1 << log2_number_of_samples>>(
residuals, Δt);
auto const mode = fft->Mode();