GeographicLib  2.1.1
DST.cpp
Go to the documentation of this file.
1 /**
2  * \file DST.cpp
3  * \brief Implementation 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 #include <GeographicLib/DST.hpp>
11 #include <vector>
12 
13 #include "kissfft.hh"
14 
15 namespace GeographicLib {
16 
17  using namespace std;
18 
19  DST::DST(int N)
20  : _N(N < 0 ? 0 : N)
21  , _fft(make_shared<fft_t>(fft_t(2 * _N, false)))
22  {}
23 
24  void DST::reset(int N) {
25  N = N < 0 ? 0 : N;
26  if (N == _N) return;
27  _N = N;
28  _fft->assign(2 * _N, false);
29  }
30 
31  void DST::fft_transform(real data[], real F[], bool centerp) const {
32  // Implement DST-III (centerp = false) or DST-IV (centerp = true).
33 
34  // Elements (0,N], resp. [0,N), of data should be set on input for centerp
35  // = false, resp. true. F must have a size of at least N and on output
36  // elements [0,N) of F contain the transform.
37  if (_N == 0) return;
38  if (centerp) {
39  for (int i = 0; i < _N; ++i) {
40  data[_N+i] = data[_N-1-i];
41  data[2*_N+i] = -data[i];
42  data[3*_N+i] = -data[_N-1-i];
43  }
44  } else {
45  data[0] = 0; // set [0]
46  for (int i = 1; i < _N; ++i) data[_N+i] = data[_N-i]; // set [N+1,2*N-1]
47  for (int i = 0; i < 2*_N; ++i) data[2*_N+i] = -data[i]; // [2*N, 4*N-1]
48  }
49  vector<complex<real>> ctemp(2*_N);
50  _fft->transform_real(data, ctemp.data());
51  if (centerp) {
52  real d = -Math::pi()/(4*_N);
53  for (int i = 0, j = 1; i < _N; ++i, j+=2)
54  ctemp[j] *= exp(complex<real>(0, j*d));
55  }
56  for (int i = 0, j = 1; i < _N; ++i, j+=2) {
57  F[i] = -ctemp[j].imag() / (2*_N);
58  }
59  }
60 
61  void DST::fft_transform2(real data[], real F[]) const {
62  // Elements [0,N), of data should be set to the N grid center values and F
63  // should have size of at least 2*N. On input elements [0,N) of F contain
64  // the size N transform; on output elements [0,2*N) of F contain the size
65  // 2*N transform.
66  fft_transform(data, F+_N, true);
67  // Copy DST-IV order N tx to [0,N) elements of data
68  for (int i = 0; i < _N; ++i) data[i] = F[i+_N];
69  for (int i = _N; i < 2*_N; ++i)
70  // (DST-IV order N - DST-III order N) / 2
71  F[i] = (data[2*_N-1-i] - F[2*_N-1-i])/2;
72  for (int i = 0; i < _N; ++i)
73  // (DST-IV order N + DST-III order N) / 2
74  F[i] = (data[i] + F[i])/2;
75  }
76 
77  void DST::transform(function<real(real)> f, real F[]) const {
78  vector<real> data(4 * _N);
79  real d = Math::pi()/(2 * _N);
80  for (int i = 1; i <= _N; ++i)
81  data[i] = f( i * d );
82  fft_transform(data.data(), F, false);
83  }
84 
85  void DST::refine(function<real(real)> f, real F[]) const {
86  vector<real> data(4 * _N);
87  real d = Math::pi()/(4 * _N);
88  for (int i = 0; i < _N; ++i)
89  data[i] = f( (2*i + 1) * d );
90  fft_transform2(data.data(), F);
91  }
92 
93  Math::real DST::eval(real sinx, real cosx, const real F[], int N) {
94  // Evaluate
95  // y = sum(F[i] * sin((2*i+1) * x), i, 0, N-1)
96  // using Clenshaw summation.
97  // Approx operation count = (N + 5) mult and (2 * N + 2) add
98  real
99  ar = 2 * (cosx - sinx) * (cosx + sinx), // 2 * cos(2 * x)
100  y0 = N & 1 ? F[--N] : 0, y1 = 0; // accumulators for sum
101  // Now N is even
102  while (N > 0) {
103  // Unroll loop x 2, so accumulators return to their original role
104  y1 = ar * y0 - y1 + F[--N];
105  y0 = ar * y1 - y0 + F[--N];
106  }
107  return sinx * (y0 + y1); // sin(x) * (y0 + y1)
108  }
109 
110  Math::real DST::integral(real sinx, real cosx, const real F[], int N) {
111  // Evaluate
112  // y = -sum(F[i]/(2*i+1) * cos((2*i+1) * x), i, 0, N-1)
113  // using Clenshaw summation.
114  // Approx operation count = (N + 5) mult and (2 * N + 2) add
115  int l = N;
116  real
117  ar = 2 * (cosx - sinx) * (cosx + sinx), // 2 * cos(2 * x)
118  y0 = N & 1 ? F[--N]/(2*(--l)+1) : 0, y1 = 0; // accumulators for sum
119  // Now N is even
120  while (N > 0) {
121  // Unroll loop x 2, so accumulators return to their original role
122  y1 = ar * y0 - y1 + F[--N]/(2*(--l)+1);
123  y0 = ar * y1 - y0 + F[--N]/(2*(--l)+1);
124  }
125  return cosx * (y1 - y0); // cos(x) * (y1 - y0)
126  }
127 
128 } // namespace GeographicLib
Header for GeographicLib::DST class.
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
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
static T pi()
Definition: Math.hpp:190
Definition: DST.hpp:19
Namespace for GeographicLib.
Definition: Accumulator.cpp:12