3 #include "gsl/gsl_fft_real.h"
4 #include "gsl/gsl_fft_halfcomplex.h"
13 const int M = LineDensity.size();
14 double max_four_coef = 0.0;
16 gsl_fft_real_wavetable *
real = gsl_fft_real_wavetable_alloc(M);
17 gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(M);
18 double *LD =
new double[M];
20 for(
int i = 0; i < M; ++ i) {
21 LD[i] = LineDensity[i];
23 gsl_fft_real_transform(LD, 1, M, real, work);
25 gsl_fft_real_wavetable_free(real);
27 for(
int i = 0; i < M; ++ i) {
28 if(
fabs(LD[i]) > max_four_coef) {
29 max_four_coef =
fabs(LD[i]);
34 for(
int i = 0; i < M; ++ i) {
35 if(
fabs(LD[i]) < max_four_coef) {
40 gsl_fft_halfcomplex_wavetable *hc = gsl_fft_halfcomplex_wavetable_alloc(M);
42 gsl_fft_halfcomplex_inverse(LD, 1, M, hc, work);
44 gsl_fft_halfcomplex_wavetable_free(hc);
45 gsl_fft_real_workspace_free(work);
47 for(
int i = 0; i < M; ++ i) {
48 LineDensity[i] = LD[i];
55 const int M = LineDensity.size();
57 double max_four_coef = 0.0;
60 gsl_fft_real_wavetable *
real = gsl_fft_real_wavetable_alloc(M);
61 gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(M);
62 double *LD =
new double[M];
64 for(
int i = 0; i < M; ++ i) {
65 LD[i] = LineDensity[i];
67 gsl_fft_real_transform(LD, 1, M, real, work);
69 gsl_fft_real_wavetable_free(real);
71 for(
int i = 1; i < M; ++ i) {
72 if(
fabs(LD[i]) > max_four_coef) {
73 max_four_coef =
fabs(LD[i]);
79 for(
int i = 1; i < M; i += 2) {
81 if(
fabs(LD[i+1]) > max_four_coef) {
83 LD[i] = -LD[i+1] * gff * (i + 1) / 2;
87 if(
fabs(temp) > max_four_coef) {
88 LD[i+1] = temp * gff * (i + 1) / 2;
94 gsl_fft_halfcomplex_wavetable *hc = gsl_fft_halfcomplex_wavetable_alloc(M);
96 gsl_fft_halfcomplex_inverse(LD, 1, M, hc, work);
98 gsl_fft_halfcomplex_wavetable_free(hc);
99 gsl_fft_real_workspace_free(work);
101 for(
int i = 0; i < M; ++ i) {
102 LineDensity[i] = LD[i];
PETE_TUTree< FnFabs, typename T::PETE_Expr_t > fabs(const PETE_Expr< T > &l)
FLieGenerator< T, N > real(const FLieGenerator< std::complex< T >, N > &)
Take real part of a complex generator.
void calc_derivative(std::vector< double > &histogram, const double &h)
constexpr double pi
The value of .
RelativeFFTLowPassFilter(const double &threshold)
void apply(std::vector< double > &LineDensity)