summaryrefslogtreecommitdiff
path: root/fourier.hpp
blob: 170866714681366b477f45094565bdd41a6e3976 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#pragma once
#include <cmath>
#include <fftw3.h>
#include <complex>
#include <vector>
#include <tuple>
#include <gsl/gsl_sf_gamma.h>

#ifndef FFTW_THREADS
#define FFTW_THREADS 1
#endif

using Real = double;
using Complex = std::complex<Real>;
using namespace std::complex_literals;

Real f(Real λ, unsigned p, unsigned s, Real q);
Real df(Real λ, unsigned p, unsigned s, Real q);
Real ddf(Real λ, unsigned p, unsigned s, Real q);

class FourierTransform {
private:
  Real* a;
  Complex* â;
  fftw_plan plan_r2c;
  fftw_plan plan_c2r;
  unsigned n;
  Real Δω;
  Real Δτ;
public:
  FourierTransform(unsigned n, Real Δω, Real Δτ, unsigned flags = 0);
  ~FourierTransform();
  std::vector<Complex> fourier(const std::vector<Real>& c);
  std::vector<Complex> fourier();
  std::vector<Real> inverse(const std::vector<Complex>& ĉ);
  void writeToA(unsigned i, Real ai);
  std::vector<Real> convolve(const std::vector<Real>& Γh, const std::vector<Real>& R);
};

class LogarithmicFourierTransform {
private:
  Complex* a;
  Complex* â;
  fftw_plan a_to_â;
  fftw_plan â_to_a;
  unsigned N;
  unsigned pad;
  Real k;
  Real Δτ;
  Real τₛ;
  Real ωₛ;
  Real sₛ;
public:
  LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad = 4);
  ~LogarithmicFourierTransform();
  Real τ(unsigned n) const;
  Real ω(unsigned n) const;
  Real t(unsigned n) const;
  Real ν(unsigned n) const;
  Real s(unsigned n) const;
  std::vector<Complex> fourier(const std::vector<Real>& c, bool symmetric);
  std::vector<Real> inverse(const std::vector<Complex>& ĉ);
};

Complex gamma(Complex z);

std::string fourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real y, unsigned log2n, Real τₘₐₓ);
Real energy(const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ, Real y, Real Δτ);
std::tuple<std::vector<Complex>, std::vector<Complex>> RddfCtdfCt(FourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ);
Real estimateZ(FourierTransform& fft, const std::vector<Real>& C, const std::vector<Complex>& Ct, const std::vector<Real>& R, const std::vector<Complex>& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real y);