3 #include "gsl/gsl_fft_real.h"
4 #include "gsl/gsl_fft_halfcomplex.h"
9 number_frequencies_m(N)
13 const int M = LineDensity.size();
14 gsl_fft_halfcomplex_wavetable *hc;
15 gsl_fft_real_wavetable *
real = gsl_fft_real_wavetable_alloc(M);
16 gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(M);
17 double *LD =
new double[M];
18 for(
int i = 0; i < M; ++ i) {
19 LD[i] = LineDensity[i];
21 gsl_fft_real_transform(LD, 1, M, real, work);
27 gsl_fft_real_wavetable_free(real);
28 hc = gsl_fft_halfcomplex_wavetable_alloc(M);
30 gsl_fft_halfcomplex_inverse(LD, 1, M, hc, work);
32 gsl_fft_halfcomplex_wavetable_free(hc);
33 gsl_fft_real_workspace_free(work);
35 for(
int i = 0; i < M; ++ i) {
36 LineDensity[i] = LD[i];
42 const int M = LineDensity.size();
43 const double gff = 2. *
Physics::pi / (h * (M - 1));
46 gsl_fft_halfcomplex_wavetable *hc;
47 gsl_fft_real_wavetable *
real = gsl_fft_real_wavetable_alloc(M);
48 gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(M);
49 double *LD =
new double[M];
51 for(
int i = 0; i < M; ++ i) {
52 LD[i] = LineDensity[i];
54 gsl_fft_real_transform(LD, 1, M, real, work);
56 gsl_fft_real_wavetable_free(real);
61 LD[i] = -LD[i+1] * gff * (i + 1) / 2;
62 LD[i+1] = temp * gff * (i + 1) / 2;
65 for(
int i = 2 * number_frequencies_m + 1; i < M; ++ i) {
69 hc = gsl_fft_halfcomplex_wavetable_alloc(M);
71 gsl_fft_halfcomplex_inverse(LD, 1, M, hc, work);
73 gsl_fft_halfcomplex_wavetable_free(hc);
74 gsl_fft_real_workspace_free(work);
76 for(
int i = 0; i < M; ++ i) {
77 LineDensity[i] = LD[i];
FixedFFTLowPassFilter(const int &N)
FLieGenerator< T, N > real(const FLieGenerator< std::complex< T >, N > &)
Take real part of a complex generator.
constexpr double pi
The value of .
void apply(std::vector< double > &LineDensity)
void calc_derivative(std::vector< double > &histogram, const double &h)