1+ #include < dsplib/fft.h>
2+ #include < dsplib/ifft.h>
3+
4+ #include " dsplib/array.h"
5+ #include " fftw3.h"
6+
7+ namespace dsplib {
8+
9+ namespace {
10+
11+ // ------------------------------------------------------------------------------------------------------------------------
12+ // complex FFT/IFFT
13+ class FFTWPlanC : public FftPlanC
14+ {
15+ public:
16+ explicit FFTWPlanC (int n, bool forward = true )
17+ : n_{n}
18+ , scale_(forward ? 1 : (1.0 / n_))
19+ , out_(n) {
20+ const int sign = forward ? FFTW_FORWARD : FFTW_BACKWARD;
21+ in_ = (fftw_complex*)fftw_malloc (sizeof (fftw_complex) * n);
22+ plan_ = fftw_plan_dft_1d (n, in_, reinterpret_cast <fftw_complex*>(out_.data ()), sign,
23+ FFTW_MEASURE | FFTW_DESTROY_INPUT);
24+ }
25+
26+ virtual ~FFTWPlanC () {
27+ fftw_destroy_plan (plan_);
28+ fftw_free (in_);
29+ }
30+
31+ dsplib::arr_cmplx solve (const dsplib::arr_cmplx& x) const final {
32+ DSPLIB_ASSERT (x.size () == n_, " input size must be equal `n`"
33+ std::memcpy (in_, x.data (), n_ * sizeof (x[0 ]));
34+ fftw_execute (plan_);
35+ out_ *= scale_;
36+ return out_;
37+ }
38+
39+ int size () const noexcept final {
40+ return n_;
41+ }
42+
43+ private:
44+ const int n_;
45+ double scale_{1.0 };
46+ fftw_complex* in_{nullptr };
47+ mutable dsplib::arr_cmplx out_;
48+ fftw_plan plan_;
49+ };
50+
51+ // ------------------------------------------------------------------------------------------------------------------------
52+ // real FFT
53+ class FFTWPlanR : public FftPlanR
54+ {
55+ public:
56+ explicit FFTWPlanR (int n)
57+ : n_{n}
58+ , out_(n) {
59+ in_ = (double *)fftw_malloc (sizeof (double ) * n);
60+ plan_ =
61+ fftw_plan_dft_r2c_1d (n, in_, reinterpret_cast <fftw_complex*>(out_.data ()), FFTW_MEASURE | FFTW_DESTROY_INPUT);
62+ }
63+
64+ virtual ~FFTWPlanR () {
65+ fftw_destroy_plan (plan_);
66+ fftw_free (in_);
67+ }
68+
69+ dsplib::arr_cmplx solve (const dsplib::arr_real& x) const final {
70+ DSPLIB_ASSERT (x.size () == n_, " input size must be equal `n`"
71+ std::memcpy (in_, x.data (), n_ * sizeof (x[0 ]));
72+ fftw_execute (plan_);
73+ for (size_t i = 1 ; i < n_ / 2 ; ++i) {
74+ out_[n_ - i] = out_[i].conj ();
75+ }
76+ return out_;
77+ }
78+
79+ int size () const noexcept final {
80+ return n_;
81+ }
82+
83+ private:
84+ const int n_;
85+ double * in_{nullptr };
86+ mutable dsplib::arr_cmplx out_;
87+ fftw_plan plan_;
88+ };
89+
90+ // ------------------------------------------------------------------------------------------------------------------------
91+ // real IFFT
92+ class IFFTWPlanR : public IfftPlanR
93+ {
94+ public:
95+ explicit IFFTWPlanR (int n)
96+ : n_{n}
97+ , scale_{1.0 / n_}
98+ , out_(n) {
99+ in_ = (fftw_complex*)fftw_malloc (sizeof (fftw_complex) * (n_ / 2 + 1 ));
100+ plan_ = fftw_plan_dft_c2r_1d (n, in_, out_.data (), FFTW_MEASURE | FFTW_DESTROY_INPUT);
101+ }
102+
103+ virtual ~IFFTWPlanR () {
104+ fftw_destroy_plan (plan_);
105+ fftw_free (in_);
106+ }
107+
108+ dsplib::arr_real solve (const dsplib::arr_cmplx& x) const final {
109+ const int n2 = n_ / 2 + 1 ;
110+ DSPLIB_ASSERT ((x.size () == n_) || (x.size () == n2), " input size must be equal `n` or `n/2+1`"
111+ std::memcpy (in_, x.data (), n2 * sizeof (x[0 ]));
112+ fftw_execute (plan_);
113+ out_ *= scale_;
114+ return out_;
115+ }
116+
117+ int size () const noexcept final {
118+ return n_;
119+ }
120+
121+ private:
122+ const int n_;
123+ const double scale_;
124+ fftw_complex* in_{nullptr };
125+ mutable dsplib::arr_real out_;
126+ fftw_plan plan_;
127+ };
128+
129+ } // namespace
130+
131+ std::shared_ptr<FftPlanC> fft_plan_c (int n) {
132+ return std::make_shared<FFTWPlanC>(n, true );
133+ }
134+
135+ std::shared_ptr<FftPlanR> fft_plan_r (int n) {
136+ return std::make_shared<FFTWPlanR>(n);
137+ }
138+
139+ std::shared_ptr<IfftPlanC> ifft_plan_c (int n) {
140+ return std::make_shared<FFTWPlanC>(n, false );
141+ }
142+
143+ std::shared_ptr<IfftPlanR> ifft_plan_r (int n) {
144+ return std::make_shared<IFFTWPlanR>(n);
145+ }
146+
147+ } // namespace dsplib
0 commit comments