FFT results of c++ vs MATLAB

Pages: 12
mbozzi (3932)
I just I am trying to reproduce it using FFTW library and sort of failing


See the following. The important part is just the examples of how to make the plans.

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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
 
#include <cmath>
#include <cassert>
#include <iostream> 
#include <iomanip>
#include <utility>
#include <complex>

#include "fftw3.h"

int constexpr ny = 11;
int constexpr nx = 10;
int constexpr nxk = nx / 2 + 1;
int constexpr nyk = ny / 2 + 1;

unsigned constexpr flags = FFTW_ESTIMATE;

double const pi = 3.1415926535897932384626433;

struct cplx_buffer
{
  fftw_complex* a;
  int rows;
  int cols;

  fftw_complex& operator()(int i, int j) const { return a[i * cols + j]; }
};

struct real_buffer
{
  double* a;
  int rows;
  int cols;

  double& operator()(int i, int j) const { return a[i * cols + j]; }
};

real_buffer my_fftw_allocate_real(int x, int y) { return real_buffer{ fftw_alloc_real(x * y), x, y }; }
cplx_buffer my_fftw_allocate_cplx(int x, int y) { return cplx_buffer{ fftw_alloc_complex(x * y), x, y }; }


void print_matrix(real_buffer const& m)
{
  std::cout << std::fixed << std::setprecision(2);

  for (int i = 0; i < m.rows; i++)
  {
    for (int j = 0; j < m.cols; j++)
    {
      std::cout << std::setw(16) << m(i, j);
    }
    std::cout << '\n';
  }
  std::cout << '\n';
}

void print_matrix(cplx_buffer const& m)
{
  std::cout << std::fixed << std::setprecision(2);

  for (int i = 0; i < m.rows; i++)
  {
    for (int j = 0; j < m.cols; j++)
    {
      std::cout << std::setw(16) << std::complex<double> { m(i, j)[0], m(i, j)[1] };
    }
    std::cout << '\n';
  }
  std::cout << '\n';
}

int main()
{
  real_buffer in = my_fftw_allocate_real(ny, nx);
  for (int i = 0; i < ny; ++i)
    for (int j = 0; j < nx; ++j)
    {
      double const XXij = +std::sin(j * 2.0 * pi / nx);
      double const YYij = -std::cos(i * pi / ny) - 0.5;
      in(i, j) = YYij * YYij * XXij;
    }
  std::cout << "input matrix\n"; print_matrix(in);

  { // Transform the first row
    cplx_buffer out = my_fftw_allocate_cplx(1, nxk);
    fftw_execute(fftw_plan_dft_r2c_1d(nx, in.a, out.a, flags));
    std::cout << "DFT of first row of input\n"; print_matrix(out);
  }

  { // Transform all the rows
    cplx_buffer out = my_fftw_allocate_cplx(ny, nxk);
    fftw_execute(fftw_plan_many_dft_r2c(1, &nx, in.rows, in.a, &in.cols, 1, in.cols, out.a, &in.cols, 1, out.cols, flags));
    std::cout << "DFT of input, row-wise\n"; print_matrix(out);
  }

  { // Transform the first column 
    cplx_buffer out = my_fftw_allocate_cplx(nyk, 1);
    fftw_execute(fftw_plan_many_dft_r2c(1, &in.rows, 1, in.a, &in.rows, in.cols, 1, out.a, &in.rows, out.cols, 1, flags));
    std::cout << "DFT of first column of input\n";  print_matrix(out);
  }

  { // Transform all the columns: 
    cplx_buffer out = my_fftw_allocate_cplx(nyk, nx);
    fftw_execute(fftw_plan_many_dft_r2c(1, &in.rows, in.cols, in.a, &in.rows, in.cols, 1, out.a, &in.rows, out.cols, 1, flags));
    std::cout << "DFT of input, column-wise\n"; print_matrix(out);
  }

  { // Do a 2D FFT the hard way, first rows, then columns
    cplx_buffer tmp = my_fftw_allocate_cplx(ny, nxk);
    cplx_buffer out = my_fftw_allocate_cplx(ny, nxk);
    fftw_execute(fftw_plan_many_dft_r2c(1, &in.cols, in.rows, in.a, &in.cols, 1, in.cols, tmp.a, &in.cols, 1, tmp.cols, flags));
    fftw_execute(fftw_plan_many_dft(1, &in.rows, tmp.cols, tmp.a, &in.rows, tmp.cols, 1, out.a, &in.rows, out.cols, 1, FFTW_FORWARD, flags));
    std::cout << "2D transform, the hard way\n"; print_matrix(out);
  }

  { // Do a 2D FFT the easy way:
    cplx_buffer out = my_fftw_allocate_cplx(ny, nxk);
    fftw_execute(fftw_plan_dft_r2c_2d(in.rows, in.cols, in.a, out.a, flags));
    std::cout << "2D transform, the easy way\n"; print_matrix(out);
  }

  return 0;
}

Last edited on
mbozzi (3932)
Here's a version of the program that doesn't leak memory or fftw_plans, is a little more ergonomic, and marginally less error prone. As before the important part is the main function. 

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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
 
#include <cmath>
#include <cassert>
#include <iostream> 
#include <iomanip>
#include <utility>
#include <complex>

#include "fftw3.h"

template <typename T>
  class my_fftw_buffer_base
  {
    T* a;
    int m; 
    int n;
  
  protected:
    ~my_fftw_buffer_base() { if (a) fftw_free(a); };
  
  public:
    using value_type = T;
  
    my_fftw_buffer_base(int m, int n)
      : a{(T*)fftw_malloc(m * n * sizeof T)}, m{ m }, n{ n }
    {       
      if (a == nullptr) abort(); 
    }
  
    my_fftw_buffer_base& operator=(my_fftw_buffer_base const&) = delete;
    my_fftw_buffer_base(my_fftw_buffer_base const&) = delete;
  
    my_fftw_buffer_base& operator=(my_fftw_buffer_base&& that) noexcept
    {
      if (this == &that) return;
      this->a = std::exchange(that.a, nullptr);
      this->m = std::exchange(that.m, 0);
      this->n = std::exchange(that.n, 0);
      return *this;
    }
    
    my_fftw_buffer_base(my_fftw_buffer_base&& that) noexcept
    {
      *this = std::move(that);
    }
  
    int rows() const noexcept { return m; }
    int cols() const noexcept { return n; }
  
    T&       operator[](int n)              { return checked_access(n); }
    T const& operator[](int n) const        { return checked_access(n); }
    T&       operator()(int n)              { return checked_access(n); }
    T const& operator()(int n) const        { return checked_access(n); }
    T&       operator()(int i, int j)       { return checked_access(i, j); }
    T const& operator()(int i, int j) const { return checked_access(i, j); }
  
    T*       data()      noexcept       { return a; }
    T const* data()      const noexcept { return a; }
    operator T* ()       noexcept       { return a; }
    operator T const* () const noexcept { return a; }
  
    int size() const noexcept { return m * n; }
  
  private:
    T& checked_access(int i, int j) const
    { 
      assert(i >= 0 && i < m);
      assert(j >= 0 && j < n);
      return a[i * n + j];
    }
  
    T& checked_access(int i) const
    {
      assert(i >= 0 && i < size());
      return a[i];
    }
  
  protected: 
    T* unsafe_get_ptr() const noexcept { return a; };
  };

template <typename T>
  struct my_fftw_buffer_impl : my_fftw_buffer_base<T>
  { 
    using base_type = my_fftw_buffer_base<T>;
    using base_type::base_type;
    operator bool() = delete; 
    operator bool() const = delete;
  };

struct my_fftw_buffer_real
  : public my_fftw_buffer_impl<double>
{
  using base_type = my_fftw_buffer_impl<double>;
  using base_type::base_type;
};

struct my_fftw_buffer_complex 
  : public my_fftw_buffer_impl<std::complex<double>>
{
  using base_type = my_fftw_buffer_impl<std::complex<double>>;
  using base_type::base_type;

  fftw_complex* data()            noexcept       { return array_oriented(); }
  fftw_complex const* data()      const noexcept { return array_oriented(); }
  operator fftw_complex* ()       noexcept       { return array_oriented(); }
  operator fftw_complex const* () const noexcept { return array_oriented(); }

private:
  fftw_complex* array_oriented() const 
  { // This conversion is likely okay, see [complex.numbers.general]/4
    return reinterpret_cast<fftw_complex*>(unsafe_get_ptr()); 
  }
};

class my_fftw_plan
{
  fftw_plan p;

public:
  ~my_fftw_plan() { if (p) fftw_destroy_plan(p); }

  explicit my_fftw_plan(fftw_plan p)
        : p(p) { if (!p) abort(); }

  my_fftw_plan& operator=(my_fftw_plan&& that) noexcept
  { if (this != &that) this->p = std::exchange(that.p, nullptr); }
  my_fftw_plan(my_fftw_plan&& that) noexcept { *this = std::move(that); }
  my_fftw_plan& operator=(my_fftw_plan const&) = delete;
  my_fftw_plan(my_fftw_plan const&) = delete;
  operator fftw_plan() const noexcept { return p; }
};

template <typename T>
  void print_matrix(T const& m, typename T::value_type q = 1.0)
  {
    std::cout << std::fixed << std::setprecision(2); 
  
    for (int i = 0; i < m.rows(); i++)
    {
      for (int j = 0; j < m.cols(); j++) 
      {
        std::cout << std::setw(16) << (m(i, j) / q);
      }
      std::cout << '\n';
    }
    std::cout << '\n';
  }

int constexpr ny = 11;
int constexpr nx = 10;
int constexpr nxk = nx / 2 + 1;
int constexpr nyk = ny / 2 + 1;

unsigned constexpr flags = FFTW_ESTIMATE;

double const pi = 3.1415926535897932384626433;

int main()
{
  my_fftw_buffer_real in(ny, nx);
  for (int i = 0; i < ny; ++i)
    for (int j = 0; j < nx; ++j)
    {
      double const XXij = +std::sin(j * 2.0 * pi / nx);
      double const YYij = -std::cos(i * pi / ny) - 0.5;
      in(i, j) = YYij * YYij * XXij;
    }
  std::cout << "input matrix\n"; print_matrix(in);

  { // Transform the first row
    my_fftw_buffer_complex out{ 1, nxk };
    fftw_execute(my_fftw_plan{ fftw_plan_dft_r2c_1d(nx, in, out, flags) });
    std::cout << "DFT of first row of input\n"; print_matrix(out);
  }

  { // Transform all the rows
    my_fftw_buffer_complex out{ ny, nxk };
    fftw_execute(my_fftw_plan{ fftw_plan_many_dft_r2c(1, &nx, in.rows(), in, &nx, 1, nx, out, &nx, 1, out.cols(), flags) });
    std::cout << "DFT of input, row-wise\n"; print_matrix(out);
  }

  { // Transform the first column 
    my_fftw_buffer_complex out{ nyk, 1 };
    fftw_execute(my_fftw_plan{ fftw_plan_many_dft_r2c(1, &ny, 1, in, &ny, in.cols(), 1, out, &ny, out.cols(), 1, flags) });
    std::cout << "DFT of first column of input\n";  print_matrix(out);
  }

  { // Transform all the columns: 
    my_fftw_buffer_complex out{ nyk, nx };
    fftw_execute(my_fftw_plan{ fftw_plan_many_dft_r2c(1, &ny, in.cols(), in, &ny, in.cols(), 1, out, &ny, out.cols(), 1, flags)});
    std::cout << "DFT of input, column-wise\n"; print_matrix(out);
  }

  { // Do a 2D FFT the hard way, first rows, then columns
    my_fftw_buffer_complex tmp{ ny, nxk };
    my_fftw_buffer_complex out{ ny, nxk };
    fftw_execute(my_fftw_plan{ fftw_plan_many_dft_r2c(1, &nx, in.rows(), in, &nx, 1, nx, tmp, &nx, 1, tmp.cols(), flags) });
    fftw_execute(my_fftw_plan{ fftw_plan_many_dft(1, &ny, tmp.cols(), tmp, &ny, tmp.cols(), 1, out, &ny, out.cols(), 1, FFTW_FORWARD, flags) });
    std::cout << "2D transform, the hard way\n"; print_matrix(out);
  }

  { // Do a 2D FFT the easy way:
    my_fftw_buffer_complex out{ ny, nxk };
    fftw_execute(my_fftw_plan{ fftw_plan_dft_r2c_2d(in.rows(), in.cols(), in, out, flags) });
    std::cout << "2D transform, the easy way\n"; print_matrix(out);
  }
}

Last edited on
JamieAl (236)
This is pretty HELPFUL, thanks a lot everyone!
Topic archived. No new replies allowed.
Pages: 12