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