GeographicLib: GeographicLib::DST Class Reference (original) (raw)
Discrete sine transforms.
This decomposes periodic functions \( f(\sigma) \) (period \( 2\pi \)) which are odd about \( \sigma = 0 \) and even about \( \sigma = \frac12 \pi \) into a Fourier series
\[ f(\sigma) = \sum_{l=0}^\infty F_l \sin\bigl((2l+1)\sigma\bigr). \]
The first \( N \) components of \( F_l \), for \(0 \le l < N\) may be approximated by
\[ F_l = \frac2N \sum_{j=1}^{N} p_j f(\sigma_j) \sin\bigl((2l+1)\sigma_j\bigr), \]
where \( \sigma_j = j\pi/(2N) \) and \( p_j = \frac12 \) for \( j = N \) and \( 1 \) otherwise. \( F_l \) is a discrete sine transform of type DST-III and may be conveniently computed using the fast Fourier transform, FFT; this is implemented with the DST::transform method.
Having computed \( F_l \) based on \( N \) evaluations of \( f(\sigma) \) at \( \sigma_j = j\pi/(2N) \), it is possible to refine these transform values and add another \( N \) coefficients by evaluating \( f(\sigma) \) at \( (j-\frac12)\pi/(2N) \); this is implemented with the DST::refine method.
Here we compute FFTs using the kissfft package https://github.com/mborgerding/kissfft by Mark Borgerding.
Example of use:
#include
#include
#include
using namespace std;
class sawtooth {
private:
double _a;
public:
sawtooth(double a) : _a(a) {}
double operator()(double x) const { return _a * x; }
};
try {
sawtooth f(Math::pi()/4);
int N = 5, K = 2*N;
vector tx(N), txa(2*N);
dst.reset(N);
dst.transform(f, tx.data());
cout << "Transform of sawtooth based on " << N << " points\n"
<< "approx 1, -1/9, 1/25, -1/49, ...\n";
for (int i = 0; i < min(K,N); ++i) {
int j = (2*i+1)*(2*i+1)*(1-((i&1)<<1));
cout << i << " " << tx[i] << " " << tx[i]*j << "\n";
}
tx.resize(2*N);
dst.refine(f, tx.data());
cout << "Add another " << N << " points\n";
for (int i = 0; i < min(K,2*N); ++i) {
int j = (2*i+1)*(2*i+1)*(1-((i&1)<<1));
cout << i << " " << tx[i] << " " << tx[i]*j << "\n";
}
dst.reset(2*N);
dst.transform(f, txa.data());
cout << "Retransform of sawtooth based on " << 2*N << " points\n";
for (int i = 0; i < min(K,2*N); ++i) {
int j = (2*i+1)*(2*i+1)*(1-((i&1)<<1));
cout << i << " " << txa[i] << " " << txa[i]*j << "\n";
}
cout << "Table of values and integral\n";
for (int i = 0; i <= K; ++i) {
double x = i*Math::pi()/(2*K), sinx = sin(x), cosx = cos(x);
cout << x << " " << f(x) << " "
<< DST::eval(sinx, cosx, txa.data(), 2*N) << " "
<< DST::integral(sinx, cosx, txa.data(), 2*N) << "\n";
}
}
catch (const exception& e) {
cerr << "Caught exception: " << e.what() << "\n";
return 1;
}
}
int main(int argc, const char *const argv[])
Header for GeographicLib::DST class.
Header for GeographicLib::Math class.
Discrete sine transforms.
void transform(std::function< real(real)> f, real F[]) const
void refine(std::function< real(real)> f, real F[]) const
Namespace for GeographicLib.
Note
The FFTW package https://www.fftw.org/ can also be used. However this is a more complicated dependency, its CMake support is broken, and it doesn't work with mpreals (GEOGRAPHICLIB_PRECISION = 5).