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

Commits on Jul 5, 2020

  1. Sine.

    pleroy committed Jul 5, 2020
    Copy the full SHA
    3e852fa View commit details
  2. Hahn window.

    pleroy committed Jul 5, 2020
    Copy the full SHA
    fe1188c View commit details
  3. Hamming and Blackman.

    pleroy committed Jul 5, 2020
    Copy the full SHA
    b7c1da2 View commit details
  4. A zoo of windows.

    pleroy committed Jul 5, 2020
    Copy the full SHA
    6f78e43 View commit details
  5. ISO flat top.

    pleroy committed Jul 5, 2020
    Copy the full SHA
    aa7511a View commit details
  6. After egg's review.

    pleroy committed Jul 5, 2020
    Copy the full SHA
    b40598c View commit details
  7. Merge pull request #2628 from pleroy/Apodization

    Apodization functions
    pleroy authored Jul 5, 2020
    Copy the full SHA
    35898da View commit details
Showing with 280 additions and 7 deletions.
  1. +53 −0 numerics/apodization.hpp
  2. +142 −3 numerics/apodization_body.hpp
  3. +85 −4 numerics/apodization_test.cpp
53 changes: 53 additions & 0 deletions numerics/apodization.hpp
Original file line number Diff line number Diff line change
@@ -4,20 +4,73 @@
#include "geometry/named_quantities.hpp"
#include "numerics/poisson_series.hpp"

// The order and terminology in this file follows
// https://en.wikipedia.org/wiki/Window_function.

namespace principia {
namespace numerics {
namespace apodization {
namespace internal_apodization {

using geometry::Instant;

// ISO 18431-2:2004, section 5.4.
template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> Dirichlet(Instant const& t_min,
Instant const& t_max);

template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> Sine(Instant const& t_min,
Instant const& t_max);

// ISO 18431-2:2004, section 5.2.
template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> Hann(Instant const& t_min,
Instant const& t_max);

template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> Hamming(Instant const& t_min,
Instant const& t_max);

template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> Blackman(Instant const& t_min,
Instant const& t_max);

template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> ExactBlackman(Instant const& t_min,
Instant const& t_max);

template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> Nuttall(Instant const& t_min,
Instant const& t_max);

template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> BlackmanNuttall(Instant const& t_min,
Instant const& t_max);

template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> BlackmanHarris(Instant const& t_min,
Instant const& t_max);

// The flat-top window in Wikipedia is not normalized and comes from Matlab (?).
// We use the normalized ISO 18431-2:2004, section 5.3 instead, which is very
// close.
template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> ISO18431_2(Instant const& t_min,
Instant const& t_max);

} // namespace internal_apodization

using internal_apodization::Blackman;
using internal_apodization::BlackmanHarris;
using internal_apodization::BlackmanNuttall;
using internal_apodization::Dirichlet;
using internal_apodization::ExactBlackman;
using internal_apodization::Hann;
using internal_apodization::Hamming;
using internal_apodization::ISO18431_2;
using internal_apodization::Nuttall;
using internal_apodization::Sine;

} // namespace apodization
} // namespace numerics
145 changes: 142 additions & 3 deletions numerics/apodization_body.hpp
Original file line number Diff line number Diff line change
@@ -1,23 +1,162 @@


#pragma once

#include "numerics/apodization.hpp"

#include "geometry/barycentre_calculator.hpp"
#include "quantities/named_quantities.hpp"
#include "quantities/numbers.hpp"
#include "quantities/si.hpp"

namespace principia {
namespace numerics {
namespace apodization {
namespace internal_apodization {

using geometry::Barycentre;
using quantities::AngularFrequency;
using quantities::si::Radian;

template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> Dirichlet(Instant const& t_min,
Instant const& t_max) {
using Result = PoissonSeries<double, 0, Evaluator>;
auto const midpoint = Barycentre<Instant, double>({t_min, t_max}, {0.5, 0.5});
return Result(typename Result::Polynomial({1}, midpoint), {});
return Result(typename Result::Polynomial({1}, t_min), {});
}

template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> Sine(Instant const& t_min,
Instant const& t_max) {
using Result = PoissonSeries<double, 0, Evaluator>;
AngularFrequency const ω = π * Radian / (t_max - t_min);
return Result(typename Result::Polynomial({0}, t_min),
{{ω,
{/*sin=*/typename Result::Polynomial({1}, t_min),
/*cos=*/typename Result::Polynomial({0}, t_min)}}});
}

template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> Hann(Instant const& t_min,
Instant const& t_max) {
using Result = PoissonSeries<double, 0, Evaluator>;
AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
return Result(typename Result::Polynomial({0.5}, t_min),
{{ω,
{/*sin=*/typename Result::Polynomial({0}, t_min),
/*cos=*/typename Result::Polynomial({-0.5}, t_min)}}});
}

template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> Hamming(Instant const& t_min,
Instant const& t_max) {
using Result = PoissonSeries<double, 0, Evaluator>;
AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
return Result(
typename Result::Polynomial({25.0 / 46.0}, t_min),
{{ω,
{/*sin=*/typename Result::Polynomial({0}, t_min),
/*cos=*/typename Result::Polynomial({-21.0 / 46.0}, t_min)}}});
}

template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> Blackman(Instant const& t_min,
Instant const& t_max) {
using Result = PoissonSeries<double, 0, Evaluator>;
AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
return Result(typename Result::Polynomial({0.42}, t_min),
{{ω,
{/*sin=*/typename Result::Polynomial({0}, t_min),
/*cos=*/typename Result::Polynomial({-0.5}, t_min)}},
{2 * ω,
{/*sin=*/typename Result::Polynomial({0}, t_min),
/*cos=*/typename Result::Polynomial({0.08}, t_min)}}});
}

template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> ExactBlackman(Instant const& t_min,
Instant const& t_max) {
using Result = PoissonSeries<double, 0, Evaluator>;
AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
return Result(
typename Result::Polynomial({3969.0 / 9304.0}, t_min),
{{ω,
{/*sin=*/typename Result::Polynomial({0}, t_min),
/*cos=*/typename Result::Polynomial({-1155.0 / 2326.0}, t_min)}},
{2 * ω,
{/*sin=*/typename Result::Polynomial({0}, t_min),
/*cos=*/typename Result::Polynomial({715.0 / 9304.0}, t_min)}}});
}

template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> Nuttall(Instant const& t_min,
Instant const& t_max) {
using Result = PoissonSeries<double, 0, Evaluator>;
AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
return Result(typename Result::Polynomial({0.355768}, t_min),
{{ω,
{/*sin=*/typename Result::Polynomial({0}, t_min),
/*cos=*/typename Result::Polynomial({-0.487396}, t_min)}},
{2 * ω,
{/*sin=*/typename Result::Polynomial({0}, t_min),
/*cos=*/typename Result::Polynomial({0.144232}, t_min)}},
{3 * ω,
{/*sin=*/typename Result::Polynomial({0}, t_min),
/*cos=*/typename Result::Polynomial({-0.012604}, t_min)}}});
}

template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> BlackmanNuttall(Instant const& t_min,
Instant const& t_max) {
using Result = PoissonSeries<double, 0, Evaluator>;
AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
return Result(typename Result::Polynomial({0.3635819}, t_min),
{{ω,
{/*sin=*/typename Result::Polynomial({0}, t_min),
/*cos=*/typename Result::Polynomial({-0.4891775}, t_min)}},
{2 * ω,
{/*sin=*/typename Result::Polynomial({0}, t_min),
/*cos=*/typename Result::Polynomial({0.1365995}, t_min)}},
{3 * ω,
{/*sin=*/typename Result::Polynomial({0}, t_min),
/*cos=*/typename Result::Polynomial({-0.0106411}, t_min)}}});
}

template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> BlackmanHarris(Instant const& t_min,
Instant const& t_max) {
using Result = PoissonSeries<double, 0, Evaluator>;
AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
return Result(typename Result::Polynomial({0.35875}, t_min),
{{ω,
{/*sin=*/typename Result::Polynomial({0}, t_min),
/*cos=*/typename Result::Polynomial({-0.48829}, t_min)}},
{2 * ω,
{/*sin=*/typename Result::Polynomial({0}, t_min),
/*cos=*/typename Result::Polynomial({0.14128}, t_min)}},
{3 * ω,
{/*sin=*/typename Result::Polynomial({0}, t_min),
/*cos=*/typename Result::Polynomial({-0.01168}, t_min)}}});
}

template<template<typename, typename, int> class Evaluator>
PoissonSeries<double, 0, Evaluator> ISO18431_2(Instant const& t_min,
Instant const& t_max) {
using Result = PoissonSeries<double, 0, Evaluator>;
AngularFrequency const ω = 2 * π * Radian / (t_max - t_min);
return Result(
typename Result::Polynomial({1.0 / 4.6392}, t_min),
{{ω,
{/*sin=*/typename Result::Polynomial({0}, t_min),
/*cos=*/typename Result::Polynomial({-1.933 / 4.6392}, t_min)}},
{2 * ω,
{/*sin=*/typename Result::Polynomial({0}, t_min),
/*cos=*/typename Result::Polynomial({1.286 / 4.6392}, t_min)}},
{3 * ω,
{/*sin=*/typename Result::Polynomial({0}, t_min),
/*cos=*/typename Result::Polynomial({-0.388 / 4.6392}, t_min)}},
{4 * ω,
{/*sin=*/typename Result::Polynomial({0}, t_min),
/*cos=*/typename Result::Polynomial({0.0322 / 4.6392}, t_min)}}});
}

} // namespace internal_apodization
89 changes: 85 additions & 4 deletions numerics/apodization_test.cpp
Original file line number Diff line number Diff line change
@@ -5,30 +5,111 @@
#include "gmock/gmock.h"
#include "gtest/gtest.h"
#include "numerics/polynomial_evaluators.hpp"
#include "quantities/elementary_functions.hpp"
#include "quantities/si.hpp"
#include "testing_utilities/almost_equals.hpp"
#include "testing_utilities/vanishes_before.hpp"

namespace principia {
namespace numerics {

using geometry::Instant;
using quantities::Sqrt;
using quantities::si::Second;
using testing_utilities::AlmostEquals;
using testing_utilities::VanishesBefore;

class ApodizationTest : public ::testing::Test {
protected:
ApodizationTest() : t1_(t0_ - 1 * Second), t2_(t0_ + 2 * Second) {}
ApodizationTest()
: t1_(t0_ - 1 * Second),
t2_(t0_ + 2 * Second),
mid_(t0_ + 0.5 * Second) {}

Instant const t0_;
Instant const t1_;
Instant const t2_;
Instant const mid_;
};

TEST_F(ApodizationTest, Dirichlet) {
auto a = apodization::Dirichlet<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(1, AlmostEquals(a.Evaluate(t0_), 0));
EXPECT_THAT(1, AlmostEquals(a.Evaluate(t0_ + 0.5 * Second), 0));
EXPECT_THAT(1, AlmostEquals(a.Evaluate(t0_ + 1 * Second), 0));
EXPECT_THAT(a.Evaluate(t1_), AlmostEquals(1, 0));
EXPECT_THAT(a.Evaluate(t0_), AlmostEquals(1, 0));
EXPECT_THAT(a.Evaluate(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a.Evaluate(t2_), AlmostEquals(1, 0));
}

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

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

TEST_F(ApodizationTest, Hamming) {
auto a = apodization::Hamming<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a.Evaluate(t1_), AlmostEquals(2.0 / 23.0, 0));
EXPECT_THAT(a.Evaluate(t0_), AlmostEquals(71.0 / 92.0, 1));
EXPECT_THAT(a.Evaluate(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a.Evaluate(t2_), AlmostEquals(2.0 / 23.0, 0));
}

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

TEST_F(ApodizationTest, ExactBlackman) {
auto a = apodization::ExactBlackman<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a.Evaluate(t1_), AlmostEquals(8.0 / 1163.0, 37));
EXPECT_THAT(a.Evaluate(t0_), AlmostEquals(11843.0 / 18608.0, 1));
EXPECT_THAT(a.Evaluate(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a.Evaluate(t2_), AlmostEquals(8.0 / 1163.0, 37));
}

TEST_F(ApodizationTest, Nuttall) {
auto a = apodization::Nuttall<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a.Evaluate(t1_), VanishesBefore(1, 0));
EXPECT_THAT(a.Evaluate(t0_), AlmostEquals(0.514746, 2));
EXPECT_THAT(a.Evaluate(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a.Evaluate(t2_), VanishesBefore(1, 0));
}

TEST_F(ApodizationTest, BlackmanNuttall) {
auto a = apodization::BlackmanNuttall<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a.Evaluate(t1_), AlmostEquals(0.0003628, 703));
EXPECT_THAT(a.Evaluate(t0_), AlmostEquals(0.5292298, 1));
EXPECT_THAT(a.Evaluate(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a.Evaluate(t2_), AlmostEquals(0.0003628, 703));
}

TEST_F(ApodizationTest, BlackmanHarris) {
auto a = apodization::BlackmanHarris<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a.Evaluate(t1_), AlmostEquals(0.00006, 151));
EXPECT_THAT(a.Evaluate(t0_), AlmostEquals(0.520575, 1));
EXPECT_THAT(a.Evaluate(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a.Evaluate(t2_), AlmostEquals(0.00006, 151));
}

TEST_F(ApodizationTest, ISO18431_2) {
auto a = apodization::ISO18431_2<HornerEvaluator>(t1_, t2_);
EXPECT_THAT(a.Evaluate(t1_), AlmostEquals(-0.0028 / 4.6392, 224));
EXPECT_THAT(a.Evaluate(t0_), AlmostEquals(0.9194 / 4.6392, 6));
EXPECT_THAT(a.Evaluate(mid_), AlmostEquals(1, 0));
EXPECT_THAT(a.Evaluate(t2_), AlmostEquals(-0.0028 / 4.6392, 224));
}

} // namespace numerics