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

Commits on Nov 28, 2020

  1. Copy the full SHA
    23a1430 View commit details
  2. Change the test to construct origin-compatible series. Fix a bug in t…

    …he construction of the base. Start with more points in Clenshaw-Curtis.
    pleroy committed Nov 28, 2020
    Copy the full SHA
    8b93ba5 View commit details
  3. Adjust the tolerances.

    pleroy committed Nov 28, 2020
    Copy the full SHA
    eac6d73 View commit details
  4. More tolerances.

    pleroy committed Nov 28, 2020
    Copy the full SHA
    0d80ee7 View commit details
  5. Even more tolerances.

    pleroy committed Nov 28, 2020
    Copy the full SHA
    c5249f4 View commit details

Commits on Nov 29, 2020

  1. Fix the tests.

    pleroy committed Nov 29, 2020
    Copy the full SHA
    3227419 View commit details
  2. Merge pull request #2800 from pleroy/CentredApodization

    Change the apodization functions to be centred
    pleroy authored Nov 29, 2020
    Copy the full SHA
    5377688 View commit details
Showing with 181 additions and 140 deletions.
  1. +69 −50 numerics/apodization_body.hpp
  2. +8 −8 numerics/apodization_test.cpp
  3. +87 −72 numerics/frequency_analysis_test.cpp
  4. +2 −2 numerics/poisson_series_basis_body.hpp
  5. +14 −7 numerics/poisson_series_test.cpp
  6. +1 −1 numerics/quadrature.hpp
119 changes: 69 additions & 50 deletions numerics/apodization_body.hpp
Original file line number Diff line number Diff line change
@@ -8,6 +8,15 @@
#include "quantities/numbers.hpp"
#include "quantities/si.hpp"

// The formulæ below differs from those in
// https://en.wikipedia.org/wiki/Window_function because our functions have as
// their origin the centre of the interval, while those in Wikipedia have as
// their origin the lower bound of the interval. Specifically consider:
// f(t) = a₀ - a₁ cos ωt - a₂ cos 2ωt - a₃ cos 3ωt - a₄ cos 4ωt
// for t ∈ [0, 2π/ω]. Rewriting this formula in terms of τ = t - π/ω yields:
// f(τ) = a₀ + a₁ cos ωτ - a₂ cos 2ωτ + a₃ cos 3ωτ - a₄ cos 4ωτ
// for τ ∈ [-π/ω, π/ω].

namespace principia {
namespace numerics {
namespace apodization {
@@ -22,7 +31,8 @@ PoissonSeries<double, 0, 0, Evaluator> Dirichlet(Instant const& t_min,
Instant const& t_max) {
using Result = PoissonSeries<double, 0, 0, Evaluator>;
using AperiodicPolynomial = typename Result::AperiodicPolynomial;
return Result(AperiodicPolynomial({1}, t_min), {});
Instant const t_mid = Barycentre<Instant, int>({t_min, t_max}, {1, 1});
return Result(AperiodicPolynomial({1}, t_mid), {});
}

template<template<typename, typename, int> class Evaluator>
@@ -32,10 +42,11 @@ PoissonSeries<double, 0, 0, Evaluator> Sine(Instant const& t_min,
using AperiodicPolynomial = typename Result::AperiodicPolynomial;
using PeriodicPolynomial = typename Result::PeriodicPolynomial;
AngularFrequency const ω = π * Radian / (t_max - t_min);
return Result(AperiodicPolynomial({0}, t_min),
Instant const t_mid = Barycentre<Instant, int>({t_min, t_max}, {1, 1});
return Result(AperiodicPolynomial({0}, t_mid),
{{ω,
{/*sin=*/PeriodicPolynomial({1}, t_min),
/*cos=*/PeriodicPolynomial({0}, t_min)}}});
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial({1}, t_mid)}}});
}

template<template<typename, typename, int> class Evaluator>
@@ -45,10 +56,11 @@ PoissonSeries<double, 0, 0, Evaluator> Hann(Instant const& t_min,
using AperiodicPolynomial = typename Result::AperiodicPolynomial;
using PeriodicPolynomial = typename Result::PeriodicPolynomial;
AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
return Result(AperiodicPolynomial({0.5}, t_min),
Instant const t_mid = Barycentre<Instant, int>({t_min, t_max}, {1, 1});
return Result(AperiodicPolynomial({0.5}, t_mid),
{{ω,
{/*sin=*/PeriodicPolynomial({0}, t_min),
/*cos=*/PeriodicPolynomial({-0.5}, t_min)}}});
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial({0.5}, t_mid)}}});
}

template<template<typename, typename, int> class Evaluator>
@@ -58,10 +70,11 @@ PoissonSeries<double, 0, 0, Evaluator> Hamming(Instant const& t_min,
AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
using AperiodicPolynomial = typename Result::AperiodicPolynomial;
using PeriodicPolynomial = typename Result::PeriodicPolynomial;
return Result(AperiodicPolynomial({25.0 / 46.0}, t_min),
Instant const t_mid = Barycentre<Instant, int>({t_min, t_max}, {1, 1});
return Result(AperiodicPolynomial({25.0 / 46.0}, t_mid),
{{ω,
{/*sin=*/PeriodicPolynomial({0}, t_min),
/*cos=*/PeriodicPolynomial({-21.0 / 46.0}, t_min)}}});
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial({21.0 / 46.0}, t_mid)}}});
}

template<template<typename, typename, int> class Evaluator>
@@ -71,13 +84,14 @@ PoissonSeries<double, 0, 0, Evaluator> Blackman(Instant const& t_min,
using AperiodicPolynomial = typename Result::AperiodicPolynomial;
using PeriodicPolynomial = typename Result::PeriodicPolynomial;
AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
return Result(AperiodicPolynomial({0.42}, t_min),
Instant const t_mid = Barycentre<Instant, int>({t_min, t_max}, {1, 1});
return Result(AperiodicPolynomial({0.42}, t_mid),
{{ω,
{/*sin=*/PeriodicPolynomial({0}, t_min),
/*cos=*/PeriodicPolynomial({-0.5}, t_min)}},
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial({0.5}, t_mid)}},
{2 * ω,
{/*sin=*/PeriodicPolynomial({0}, t_min),
/*cos=*/PeriodicPolynomial({0.08}, t_min)}}});
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial({0.08}, t_mid)}}});
}

template<template<typename, typename, int> class Evaluator>
@@ -87,13 +101,14 @@ PoissonSeries<double, 0, 0, Evaluator> ExactBlackman(Instant const& t_min,
using AperiodicPolynomial = typename Result::AperiodicPolynomial;
using PeriodicPolynomial = typename Result::PeriodicPolynomial;
AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
return Result(AperiodicPolynomial({3969.0 / 9304.0}, t_min),
Instant const t_mid = Barycentre<Instant, int>({t_min, t_max}, {1, 1});
return Result(AperiodicPolynomial({3969.0 / 9304.0}, t_mid),
{{ω,
{/*sin=*/PeriodicPolynomial({0}, t_min),
/*cos=*/PeriodicPolynomial({-1155.0 / 2326.0}, t_min)}},
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial({1155.0 / 2326.0}, t_mid)}},
{2 * ω,
{/*sin=*/PeriodicPolynomial({0}, t_min),
/*cos=*/PeriodicPolynomial({715.0 / 9304.0}, t_min)}}});
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial({715.0 / 9304.0}, t_mid)}}});
}

template<template<typename, typename, int> class Evaluator>
@@ -103,16 +118,17 @@ PoissonSeries<double, 0, 0, Evaluator> Nuttall(Instant const& t_min,
using AperiodicPolynomial = typename Result::AperiodicPolynomial;
using PeriodicPolynomial = typename Result::PeriodicPolynomial;
AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
return Result(AperiodicPolynomial({0.355768}, t_min),
Instant const t_mid = Barycentre<Instant, int>({t_min, t_max}, {1, 1});
return Result(AperiodicPolynomial({0.355768}, t_mid),
{{ω,
{/*sin=*/PeriodicPolynomial({0}, t_min),
/*cos=*/PeriodicPolynomial({-0.487396}, t_min)}},
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial({0.487396}, t_mid)}},
{2 * ω,
{/*sin=*/PeriodicPolynomial({0}, t_min),
/*cos=*/PeriodicPolynomial({0.144232}, t_min)}},
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial({0.144232}, t_mid)}},
{3 * ω,
{/*sin=*/PeriodicPolynomial({0}, t_min),
/*cos=*/PeriodicPolynomial({-0.012604}, t_min)}}});
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial({0.012604}, t_mid)}}});
}

template<template<typename, typename, int> class Evaluator>
@@ -122,16 +138,17 @@ PoissonSeries<double, 0, 0, Evaluator> BlackmanNuttall(Instant const& t_min,
using AperiodicPolynomial = typename Result::AperiodicPolynomial;
using PeriodicPolynomial = typename Result::PeriodicPolynomial;
AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
return Result(AperiodicPolynomial({0.3635819}, t_min),
Instant const t_mid = Barycentre<Instant, int>({t_min, t_max}, {1, 1});
return Result(AperiodicPolynomial({0.3635819}, t_mid),
{{ω,
{/*sin=*/PeriodicPolynomial({0}, t_min),
/*cos=*/PeriodicPolynomial({-0.4891775}, t_min)}},
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial({0.4891775}, t_mid)}},
{2 * ω,
{/*sin=*/PeriodicPolynomial({0}, t_min),
/*cos=*/PeriodicPolynomial({0.1365995}, t_min)}},
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial({0.1365995}, t_mid)}},
{3 * ω,
{/*sin=*/PeriodicPolynomial({0}, t_min),
/*cos=*/PeriodicPolynomial({-0.0106411}, t_min)}}});
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial({0.0106411}, t_mid)}}});
}

template<template<typename, typename, int> class Evaluator>
@@ -141,16 +158,17 @@ PoissonSeries<double, 0, 0, Evaluator> BlackmanHarris(Instant const& t_min,
using AperiodicPolynomial = typename Result::AperiodicPolynomial;
using PeriodicPolynomial = typename Result::PeriodicPolynomial;
AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
return Result(AperiodicPolynomial({0.35875}, t_min),
Instant const t_mid = Barycentre<Instant, int>({t_min, t_max}, {1, 1});
return Result(AperiodicPolynomial({0.35875}, t_mid),
{{ω,
{/*sin=*/PeriodicPolynomial({0}, t_min),
/*cos=*/PeriodicPolynomial({-0.48829}, t_min)}},
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial({0.48829}, t_mid)}},
{2 * ω,
{/*sin=*/PeriodicPolynomial({0}, t_min),
/*cos=*/PeriodicPolynomial({0.14128}, t_min)}},
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial({0.14128}, t_mid)}},
{3 * ω,
{/*sin=*/PeriodicPolynomial({0}, t_min),
/*cos=*/PeriodicPolynomial({-0.01168}, t_min)}}});
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial({0.01168}, t_mid)}}});
}

template<template<typename, typename, int> class Evaluator>
@@ -160,23 +178,24 @@ PoissonSeries<double, 0, 0, Evaluator> ISO18431_2(Instant const& t_min,
using AperiodicPolynomial = typename Result::AperiodicPolynomial;
using PeriodicPolynomial = typename Result::PeriodicPolynomial;
AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
return Result(AperiodicPolynomial({1.0 / 4.63867187}, t_min),
Instant const t_mid = Barycentre<Instant, int>({t_min, t_max}, {1, 1});
return Result(AperiodicPolynomial({1.0 / 4.63867187}, t_mid),
{{ω,
{/*sin=*/PeriodicPolynomial({0}, t_min),
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial(
{-1.93261719 / 4.63867187}, t_min)}},
{1.93261719 / 4.63867187}, t_mid)}},
{2 * ω,
{/*sin=*/PeriodicPolynomial({0}, t_min),
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial(
{1.28613281 / 4.63867187}, t_min)}},
{1.28613281 / 4.63867187}, t_mid)}},
{3 * ω,
{/*sin=*/PeriodicPolynomial({0}, t_min),
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial(
{-0.38769531 / 4.63867187}, t_min)}},
{0.38769531 / 4.63867187}, t_mid)}},
{4 * ω,
{/*sin=*/PeriodicPolynomial({0}, t_min),
{/*sin=*/PeriodicPolynomial({0}, t_mid),
/*cos=*/PeriodicPolynomial(
{0.03222656 / 4.63867187}, t_min)}}});
{0.03222656 / 4.63867187}, t_mid)}}});
}

} // namespace internal_apodization
16 changes: 8 additions & 8 deletions numerics/apodization_test.cpp
Original file line number Diff line number Diff line change
@@ -42,24 +42,24 @@ TEST_F(ApodizationTest, Dirichlet) {

TEST_F(ApodizationTest, Sine) {
auto a = apodization::Sine<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a(t1_), AlmostEquals(0, 0));
EXPECT_THAT(a(t0_), AlmostEquals(Sqrt(3) / 2, 0));
EXPECT_THAT(a(t1_), VanishesBefore(1, 0));
EXPECT_THAT(a(t0_), AlmostEquals(Sqrt(3) / 2, 1));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), VanishesBefore(1, 1));
EXPECT_THAT(a(t2_), VanishesBefore(1, 0));
}

TEST_F(ApodizationTest, Hann) {
auto a = apodization::Hann<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a(t1_), AlmostEquals(0, 0));
EXPECT_THAT(a(t0_), AlmostEquals(0.75, 1));
EXPECT_THAT(a(t0_), AlmostEquals(0.75, 0));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), VanishesBefore(1, 0));
}

TEST_F(ApodizationTest, Hamming) {
auto a = apodization::Hamming<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a(t1_), AlmostEquals(2.0 / 23.0, 0));
EXPECT_THAT(a(t0_), AlmostEquals(71.0 / 92.0, 1));
EXPECT_THAT(a(t0_), AlmostEquals(71.0 / 92.0, 0));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), AlmostEquals(2.0 / 23.0, 0));
}
@@ -75,15 +75,15 @@ TEST_F(ApodizationTest, Blackman) {
TEST_F(ApodizationTest, ExactBlackman) {
auto a = apodization::ExactBlackman<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a(t1_), AlmostEquals(8.0 / 1163.0, 37));
EXPECT_THAT(a(t0_), AlmostEquals(11843.0 / 18608.0, 1));
EXPECT_THAT(a(t0_), AlmostEquals(11843.0 / 18608.0, 2));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), AlmostEquals(8.0 / 1163.0, 37));
}

TEST_F(ApodizationTest, Nuttall) {
auto a = apodization::Nuttall<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a(t1_), VanishesBefore(1, 0));
EXPECT_THAT(a(t0_), AlmostEquals(0.514746, 2));
EXPECT_THAT(a(t0_), AlmostEquals(0.514746, 1));
EXPECT_THAT(a(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a(t2_), VanishesBefore(1, 0));
}
@@ -107,7 +107,7 @@ TEST_F(ApodizationTest, BlackmanHarris) {
TEST_F(ApodizationTest, ISO18431_2) {
auto a = apodization::ISO18431_2<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a(t1_), AlmostEquals(-0.00195313 / 4.63867187, 440));
EXPECT_THAT(a(t0_), AlmostEquals(0.9194336 / 4.63867187, 6));
EXPECT_THAT(a(t0_), AlmostEquals(0.9194336 / 4.63867187, 5));
EXPECT_THAT(a(mid_), AlmostEquals(1, 1));
EXPECT_THAT(a(t2_), AlmostEquals(-0.00195313 / 4.63867187, 440));
}
Loading