GeographicLib: DST.cpp Source File (original) (raw)
1
2
3
4
5
6
7
8
9
11#include
12
13#include "kissfft.hh"
14
16
17 using namespace std;
18
20 : _nN(N < 0 ? 0 : N)
21 , _fft(make_shared<fft_t>(fft_t(2 * _nN, false)))
22 {}
23
26 if (N == _nN) return;
27 _nN = N;
28 _fft->assign(2 * _nN, false);
29 }
30
31 void DST::fft_transform(real data[], real F[], bool centerp) const {
32
33
34
35
36
37 if (_nN == 0) return;
38 if (centerp) {
39 for (int i = 0; i < _nN; ++i) {
40 data[_nN+i] = data[_nN-1-i];
41 data[2*_nN+i] = -data[i];
42 data[3*_nN+i] = -data[_nN-1-i];
43 }
44 } else {
45 data[0] = 0;
46 for (int i = 1; i < _nN; ++i)
47 data[_nN+i] = data[_nN-i];
48 for (int i = 0; i < 2*_nN; ++i)
49 data[2*_nN+i] = -data[i];
50 }
51 vector<complex> ctemp(2*_nN);
52 _fft->transform_real(data, ctemp.data());
53 if (centerp) {
55 for (int i = 0, j = 1; i < _nN; ++i, j+=2)
56 ctemp[j] *= exp(complex(0, j*d));
57 }
58 for (int i = 0, j = 1; i < _nN; ++i, j+=2) {
59 F[i] = -ctemp[j].imag() / (2*_nN);
60 }
61 }
62
63 void DST::fft_transform2(real data[], real F[]) const {
64
65
66
67
68 fft_transform(data, F+_nN, true);
69
70 for (int i = 0; i < _nN; ++i) data[i] = F[i+_nN];
71 for (int i = _nN; i < 2*_nN; ++i)
72
73 F[i] = (data[2*_nN-1-i] - F[2*_nN-1-i])/2;
74 for (int i = 0; i < _nN; ++i)
75
76 F[i] = (data[i] + F[i])/2;
77 }
78
80 vector data(4 * _nN);
81 real d = Math::pi()/(2 * _nN);
82 for (int i = 1; i <= _nN; ++i)
83 data[i] = f( i * d );
84 fft_transform(data.data(), F, false);
85 }
86
88 vector data(4 * _nN);
89 real d = Math::pi()/(4 * _nN);
90 for (int i = 0; i < _nN; ++i)
91 data[i] = f( (2*i + 1) * d );
92 fft_transform2(data.data(), F);
93 }
94
96
97
98
99
100 real
101 ar = 2 * (cosx - sinx) * (cosx + sinx),
102 y0 = N & 1 ? F[--N] : 0, y1 = 0;
103
104 while (N > 0) {
105
106 y1 = ar * y0 - y1 + F[--N];
107 y0 = ar * y1 - y0 + F[--N];
108 }
109 return sinx * (y0 + y1);
110 }
111
113
114
115
116
117 real
118 ar = 2 * (cosx - sinx) * (cosx + sinx),
119 y0 = 0, y1 = 0;
121 real t = ar * y0 - y1 + F[N]/(2*N+1);
122 y1 = y0; y0 = t;
123 }
124 return cosx * (y1 - y0);
125 }
126
128 const real F[], int N) {
129
130 real
131
132 ac = +2 * (cosy * cosx + siny * sinx) * (cosy * cosx - siny * sinx),
133
134 as = -2 * (siny * cosx - cosy * sinx) * (siny * cosx + cosy * sinx),
135 y0 = 0, y1 = 0, z0 = 0, z1 = 0;
137 real
138 ty = ac * y0 + as * z0 - y1 + F[N]/(2*N+1),
139 tz = as * y0 + ac * z0 - z1;
140 y1 = y0; y0 = ty;
141 z1 = z0; z0 = tz;
142 }
143
144
145
146
147
148
149 return (y1 - y0) * (cosy - cosx) + (z1 - z0) * (cosy + cosx);
150 }
151
152}
Header for GeographicLib::DST class.
GeographicLib::Math::real real
void reset(int N)
Definition DST.cpp:24
void transform(std::function< real(real)> f, real F[]) const
Definition DST.cpp:79
DST(int N=0)
Definition DST.cpp:19
static real eval(real sinx, real cosx, const real F[], int N)
Definition DST.cpp:95
void refine(std::function< real(real)> f, real F[]) const
Definition DST.cpp:87
static real integral(real sinx, real cosx, const real F[], int N)
Definition DST.cpp:112
Namespace for GeographicLib.