GeographicLib::DST Class Reference

Discrete sine transforms. More...

#include <GeographicLib/DST.hpp>

Public Member Functions
	DST (int N=0)
void	reset (int N)
int	N () const
void	transform (std::function< real(real)> f, real F[]) const
void	refine (std::function< real(real)> f, real F[]) const

Static Public Member Functions
static real	eval (real sinx, real cosx, const real F[], int N)
static real	integral (real sinx, real cosx, const real F[], int N)
static real	integral (real sinx, real cosx, real siny, real cosy, const real F[], int N)

Detailed Description

Discrete sine transforms.

This decomposes periodic functions \( f(\sigma) \) (period \( 2\pi \)) which are odd about \( \sigma = 0 \) and even about \( \sigma = \frac12 \pi \) into a Fourier series

\[ f(\sigma) = \sum_{l=0}^\infty F_l \sin\bigl((2l+1)\sigma\bigr). \]

The first \( N \) components of \( F_l \), for \(0 \le l < N\) may be approximated by

\[ F_l = \frac2N \sum_{j=1}^{N} p_j f(\sigma_j) \sin\bigl((2l+1)\sigma_j\bigr), \]

where \( \sigma_j = j\pi/(2N) \) and \( p_j = \frac12 \) for \( j = N \) and \( 1 \) otherwise. \( F_l \) is a discrete sine transform of type DST-III and may be conveniently computed using the fast Fourier transform, FFT; this is implemented with the DST::transform method.

Having computed \( F_l \) based on \( N \) evaluations of \( f(\sigma) \) at \( \sigma_j = j\pi/(2N) \), it is possible to refine these transform values and add another \( N \) coefficients by evaluating \( f(\sigma) \) at \( (j-\frac12)\pi/(2N) \); this is implemented with the DST::refine method.

Here we compute FFTs using the kissfft package https://github.com/mborgerding/kissfft by Mark Borgerding.

Example of use:

// Example of using the GeographicLib::DST class
 
#include <iostream>
#include <exception>
#include <vector>
#include <GeographicLib/Math.hpp>
#include <GeographicLib/DST.hpp>
 
using namespace std;
using namespace GeographicLib;
 
class sawtooth {
private:
  double _a;
public:
  sawtooth(double a) : _a(a) {}
  // only called for x in (0, pi/2].  DST assumes function is periodic, period
  // 2*pi, is odd about 0, and is even about pi/2.
  double operator()(double x) const { return _a * x; }
};
 
int main() {
  try {
    sawtooth f(Math::pi()/4);
    DST dst;
    int N = 5, K = 2*N;
    vector<double> tx(N), txa(2*N);
    dst.reset(N);
    dst.transform(f, tx.data());
    cout << "Transform of sawtooth based on " << N << " points\n"
         << "approx 1, -1/9, 1/25, -1/49, ...\n";
    for (int i = 0; i < min(K,N); ++i) {
      int j = (2*i+1)*(2*i+1)*(1-((i&1)<<1));
      cout << i << " " << tx[i] << " " << tx[i]*j << "\n";
    }
    tx.resize(2*N);
    dst.refine(f, tx.data());
    cout << "Add another " << N << " points\n";
    for (int i = 0; i < min(K,2*N); ++i) {
      int j = (2*i+1)*(2*i+1)*(1-((i&1)<<1));
      cout << i << " " << tx[i] << " " << tx[i]*j << "\n";
    }
    dst.reset(2*N);
    dst.transform(f, txa.data());
    cout << "Retransform of sawtooth based on " << 2*N << " points\n";
    for (int i = 0; i < min(K,2*N); ++i) {
      int j = (2*i+1)*(2*i+1)*(1-((i&1)<<1));
      cout << i << " " << txa[i] << " " << txa[i]*j << "\n";
    }
    cout << "Table of values and integral\n";
    for (int i = 0; i <= K; ++i) {
      double x = i*Math::pi()/(2*K), sinx = sin(x), cosx = cos(x);
      cout << x << " " << f(x) << " "
           << DST::eval(sinx, cosx, txa.data(), 2*N) << " "
           << DST::integral(sinx, cosx, txa.data(), 2*N) << "\n";
    }
  }
  catch (const exception& e) {
    cerr << "Caught exception: " << e.what() << "\n";
    return 1;
  }
}

Note

The FFTW package https://www.fftw.org/ can also be used. However this is a more complicated dependency, its CMake support is broken, and it doesn't work with mpreals (GEOGRAPHICLIB_PRECISION = 5).

Definition at line 63 of file DST.hpp.

Constructor & Destructor Documentation

DST()

GeographicLib::DST::DST

(

int

N = 0

)

Constructor specifying the number of points to use.

Parameters

[in]

N

the number of points to use.

Definition at line 19 of file DST.cpp.

Member Function Documentation

reset()

void GeographicLib::DST::reset

(

int

N

)

Reset the given number of points.

Parameters

[in]

N

the number of points to use.

Definition at line 24 of file DST.cpp.

References N().

Referenced by GeographicLib::GeodesicExact::GeodesicExact().

N()

int GeographicLib::DST::N ( ) const

inline

Return the number of points.

Returns

the number of points to use.

Definition at line 93 of file DST.hpp.

Referenced by eval(), integral(), and reset().

transform()

void GeographicLib::DST::transform	(	std::function< real(real)>	f,
		real	F[]
	)		const

Determine first N terms in the Fourier series

Parameters

[in]	f	the function used for evaluation.
[out]	F	the first N coefficients of the Fourier series.

The evaluates \( f(\sigma) \) at \( \sigma = (j + 1) \pi / (2 N) \) for integer \( j \in [0, N) \). F should be an array of length at least N.

Definition at line 77 of file DST.cpp.

References GeographicLib::Math::pi().

refine()

void GeographicLib::DST::refine	(	std::function< real(real)>	f,
		real	F[]
	)		const

Refine the Fourier series by doubling the number of points sampled

Parameters

[in]	f	the function used for evaluation.
[in,out]	F	on input the first N coefficents of the Fourier series; on output the refined transform based on 2N points, i.e., the first 2N coefficents.

The evaluates \( f(\sigma) \) at additional points \( \sigma = (j + \frac12) \pi / (2 N) \) for integer \( j \in [0, N) \), computes the DST-IV transform of these, and combines this with the input F to compute the 2N term DST-III discrete sine transform. This is equivalent to calling transform with twice the value of N but is more efficient, given that the N term coefficients are already known. See the example code above.

Definition at line 85 of file DST.cpp.

References GeographicLib::Math::pi().

eval()

Math::real GeographicLib::DST::eval ( real  sinx, real  cosx, const real  F[], int  N  )

static

Evaluate the Fourier sum given the sine and cosine of the angle

Parameters

[in]	sinx	sinσ.
[in]	cosx	cosσ.
[in]	F	the array of Fourier coefficients.
[in]	N	the number of Fourier coefficients.

Returns

the value of the Fourier sum.

Definition at line 93 of file DST.cpp.

References N().

integral() [1/2]

Math::real GeographicLib::DST::integral ( real  sinx, real  cosx, const real  F[], int  N  )

static

Evaluate the integral of Fourier sum given the sine and cosine of the angle

Parameters

[in]	sinx	sinσ.
[in]	cosx	cosσ.
[in]	F	the array of Fourier coefficients.
[in]	N	the number of Fourier coefficients.

Returns

the value of the integral.

The constant of integration is chosen so that the integral is zero at \( \sigma = \frac12\pi \).

Definition at line 110 of file DST.cpp.

References N().

Referenced by GeographicLib::GeodesicLineExact::GenPosition().

integral() [2/2]

Math::real GeographicLib::DST::integral ( real  sinx, real  cosx, real  siny, real  cosy, const real  F[], int  N  )

static

Evaluate the definite integral of Fourier sum given the sines and cosines of the angles at the endpoints.

Parameters

[in]	sinx	sinσ1.
[in]	cosx	cosσ1.
[in]	siny	sinσ2.
[in]	cosy	cosσ2.
[in]	F	the array of Fourier coefficients.
[in]	N	the number of Fourier coefficients.

Returns

the value of the integral.

The integral is evaluated between limits σ₁ and σ₂.

Definition at line 125 of file DST.cpp.

References N().

---

The documentation for this class was generated from the following files:

• DST.hpp
• DST.cpp