GeographicLib 2.1.1
DST.hpp
Go to the documentation of this file.
1/**
2 * \file DST.hpp
3 * \brief Header for GeographicLib::DST class
4 *
5 * Copyright (c) Charles Karney (2022) <charles@karney.com> and licensed under
6 * the MIT/X11 License. For more information, see
7 * https://geographiclib.sourceforge.io/
8 **********************************************************************/
9
10#if !defined(GEOGRAPHICLIB_DST_HPP)
11#define GEOGRAPHICLIB_DST_HPP 1
12
14
15#include <functional>
16#include <memory>
17
18template <typename scalar_t>
19class kissfft;
20
21namespace GeographicLib {
22
23 /**
24 * \brief Discrete sine transforms
25 *
26 * This decomposes periodic functions \f$ f(\sigma) \f$ (period \f$ 2\pi \f$)
27 * which are odd about \f$ \sigma = 0 \f$ and even about \f$ \sigma = \frac12
28 * \pi \f$ into a Fourier series
29 * \f[
30 * f(\sigma) = \sum_{l=0}^\infty F_l \sin\bigl((2l+1)\sigma\bigr).
31 * \f]
32 *
33 * The first \f$ N \f$ components of \f$ F_l \f$, for \f$0 \le l < N\f$ may
34 * be approximated by
35 * \f[
36 * F_l = \frac2N \sum_{j=1}^{N}
37 * p_j f(\sigma_j) \sin\bigl((2l+1)\sigma_j\bigr),
38 * \f]
39 * where \f$ \sigma_j = j\pi/(2N) \f$ and \f$ p_j = \frac12 \f$ for \f$ j = N
40 * \f$ and \f$ 1 \f$ otherwise. \f$ F_l \f$ is a discrete sine transform of
41 * type DST-III and may be conveniently computed using the fast Fourier
42 * transform, FFT; this is implemented with the DST::transform method.
43 *
44 * Having computed \f$ F_l \f$ based on \f$ N \f$ evaluations of \f$
45 * f(\sigma) \f$ at \f$ \sigma_j = j\pi/(2N) \f$, it is possible to
46 * refine these transform values and add another \f$ N \f$ coefficients by
47 * evaluating \f$ f(\sigma) \f$ at \f$ (j-\frac12)\pi/(2N) \f$; this is
48 * implemented with the DST::refine method.
49 *
50 * Here we compute FFTs using the kissfft package
51 * https://github.com/mborgerding/kissfft by Mark Borgerding.
52 *
53 * Example of use:
54 * \include example-DST.cpp
55 *
56 * \note The FFTW package https://www.fftw.org/ can also be used. However
57 * this is a more complicated dependency, its CMake support is broken, and it
58 * doesn't work with mpreals (GEOGRAPHICLIB_PRECISION = 5).
59 **********************************************************************/
60
61 class DST {
62 private:
63 typedef Math::real real;
64 int _N;
65 typedef kissfft<real> fft_t;
66 std::shared_ptr<fft_t> _fft;
67 // Implement DST-III (centerp = false) or DST-IV (centerp = true)
68 void fft_transform(real data[], real F[], bool centerp) const;
69 // Add another N terms to F
70 void fft_transform2(real data[], real F[]) const;
71 public:
72 /**
73 * Constructor specifying the number of points to use.
74 * @param[in] N the number of points to use.
75 **********************************************************************/
76 GEOGRAPHICLIB_EXPORT DST(int N = 0);
77
78 /**
79 * Reset the given number of points.
80 * @param[in] N the number of points to use.
81 **********************************************************************/
83
84 /**
85 * Return the number of points.
86 * @return the number of points to use.
87 **********************************************************************/
88 int N() const { return _N; }
89
90 /**
91 * Determine first \e N terms in the Fourier series
92 * @param[in] f the function used for evaluation.
93 * @param[out] F the first \e N coefficients of the Fourier series.
94 *
95 * The evaluates \f$ f(\sigma) \f$ at \f$ \sigma = (j + 1) \pi / (2 N) \f$
96 * for integer \f$ j \in [0, N) \f$. \e F should be an array of length at
97 * least \e N.
98 **********************************************************************/
99 void GEOGRAPHICLIB_EXPORT transform(std::function<real(real)> f, real F[])
100 const;
101
102 /**
103 * Refine the Fourier series by doubling the number of points sampled
104 * @param[in] f the function used for evaluation.
105 * @param[inout] F on input the first \e N coefficents of the Fourier
106 * series; on output the refined transform based on 2\e N points, i.e.,
107 * the first 2\e N coefficents.
108 *
109 * The evaluates \f$ f(\sigma) \f$ at additional points \f$ \sigma = (j +
110 * \frac12) \pi / (2 N) \f$ for integer \f$ j \in [0, N) \f$, computes the
111 * DST-IV transform of these, and combines this with the input \e F to
112 * compute the 2\e N term DST-III discrete sine transform. This is
113 * equivalent to calling transform with twice the value of \e N but is more
114 * efficient, given that the \e N term coefficients are already known. See
115 * the example code above.
116 **********************************************************************/
117 void GEOGRAPHICLIB_EXPORT refine(std::function<real(real)> f, real F[])
118 const;
119
120 /**
121 * Evaluate the Fourier sum given the sine and cosine of the angle
122 * @param[in] sinx sin&sigma;.
123 * @param[in] cosx cos&sigma;.
124 * @param[in] F the array of Fourier coefficients.
125 * @param[in] N the number of Fourier coefficients.
126 * @return the value of the Fourier sum.
127 **********************************************************************/
128 static real GEOGRAPHICLIB_EXPORT eval(real sinx, real cosx,
129 const real F[], int N);
130
131 /**
132 * Evaluate the integral of Fourier sum given the sine and cosine of the
133 * angle
134 * @param[in] sinx sin&sigma;.
135 * @param[in] cosx cos&sigma;.
136 * @param[in] F the array of Fourier coefficients.
137 * @param[in] N the number of Fourier coefficients.
138 * @return the value of the integral.
139 *
140 * The constant of integration is chosen so that the integral is zero at
141 * \f$ \sigma = \frac12\pi \f$.
142 **********************************************************************/
143 static real GEOGRAPHICLIB_EXPORT integral(real sinx, real cosx,
144 const real F[], int N);
145 };
146
147} // namespace GeographicLib
148
149#endif // GEOGRAPHICLIB_DST_HPP
Header for GeographicLib::Constants class.
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:67
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
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
DST(int N=0)
Definition: DST.cpp:19
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
Definition: DST.hpp:19
Namespace for GeographicLib.
Definition: Accumulator.cpp:12