21 , _fft(make_shared<
fft_t>(
fft_t(2 * _N, false)))
28 _fft->assign(2 * _N,
false);
31 void DST::fft_transform(
real data[],
real F[],
bool centerp)
const {
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];
46 for (
int i = 1; i < _N; ++i) data[_N+i] = data[_N-i];
47 for (
int i = 0; i < 2*_N; ++i) data[2*_N+i] = -data[i];
49 vector<complex<real>> ctemp(2*_N);
50 _fft->transform_real(data, ctemp.data());
53 for (
int i = 0, j = 1; i < _N; ++i, j+=2)
54 ctemp[j] *= exp(complex<real>(0, j*d));
56 for (
int i = 0, j = 1; i < _N; ++i, j+=2) {
57 F[i] = -ctemp[j].imag() / (2*_N);
61 void DST::fft_transform2(
real data[],
real F[])
const {
67 fft_transform(data, F+_N,
true);
68 for (
int i = 0; i < _N; ++i) data[i] = F[i+_N];
69 for (
int i = _N; i < 2*_N; ++i)
70 F[i] = (-data[2*_N-1-i] + F[2*_N-1-i])/2;
71 for (
int i = 0; i < _N; ++i)
72 F[i] = (data[i] + F[i])/2;
76 vector<real> data(4 * _N);
78 for (
int i = 1; i <= _N; ++i)
80 fft_transform(data.data(), F,
false);
84 vector<real> data(4 * _N);
86 for (
int i = 0; i < _N; ++i)
87 data[i] = f( (2*i + 1) * d );
88 fft_transform2(data.data(), F);
97 ar = 2 * (cosx - sinx) * (cosx + sinx),
98 y0 =
N & 1 ? F[--
N] : 0, y1 = 0;
102 y1 = ar * y0 - y1 + F[--
N];
103 y0 = ar * y1 - y0 + F[--
N];
105 return sinx * (y0 + y1);
115 ar = 2 * (cosx - sinx) * (cosx + sinx),
116 y0 =
N & 1 ? F[--
N]/(2*(--l)+1) : 0, y1 = 0;
120 y1 = ar * y0 - y1 + F[--
N]/(2*(--l)+1);
121 y0 = ar * y1 - y0 + F[--
N]/(2*(--l)+1);
123 return cosx * (y1 - y0);
Header for GeographicLib::DST class.
GeographicLib::Math::real real
void transform(std::function< real(real)> f, real F[]) const
static real eval(real sinx, real cosx, const real F[], int N)
void refine(std::function< real(real)> f, real F[]) const
static real integral(real sinx, real cosx, const real F[], int N)
Namespace for GeographicLib.