GeographicLib 2.1.1
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)
 

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>
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;
}
}
int main(int argc, const char *const argv[])
Definition: CartConvert.cpp:29
Header for GeographicLib::DST class.
Header for GeographicLib::Math class.
Discrete sine transforms.
Definition: DST.hpp:61
void reset(int N)
Definition: DST.cpp:24
void transform(std::function< real(real)> f, real F[]) const
Definition: DST.cpp:77
static real eval(real sinx, real cosx, const real F[], int N)
Definition: DST.cpp:93
int N() const
Definition: DST.hpp:88
void refine(std::function< real(real)> f, real F[]) const
Definition: DST.cpp:85
static real integral(real sinx, real cosx, const real F[], int N)
Definition: DST.cpp:110
static T pi()
Definition: Math.hpp:190
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
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 61 of file DST.hpp.

Constructor & Destructor Documentation

◆ DST()

GeographicLib::DST::DST ( int  N = 0)

Constructor specifying the number of points to use.

Parameters
[in]Nthe 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]Nthe 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 88 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]fthe function used for evaluation.
[out]Fthe 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]fthe function used for evaluation.
[in,out]Fon 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]sinxsinσ.
[in]cosxcosσ.
[in]Fthe array of Fourier coefficients.
[in]Nthe number of Fourier coefficients.
Returns
the value of the Fourier sum.

Definition at line 93 of file DST.cpp.

References N().

◆ integral()

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]sinxsinσ.
[in]cosxcosσ.
[in]Fthe array of Fourier coefficients.
[in]Nthe 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().


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