Discrete Fourier Transform C++

Question

I'm trying to write simple DFT and IDFT functions which will be my core for future projects. Trouble means that IDFT returns different value from input value, and i can't understand, where is the mistake. Below my source code:

vector<double> input; vector<double> result; vector<complex<double>> output; double IDFT(int n) { double a = 0; double b = 0; int N = output.size(); for(int k = 0; k < N; k++) { double value = abs(output[k]); a+= cos((2 * M_PI * k * n) / N) * value; b+= sin((2 * M_PI * k * n) / N) * value; } complex<double> temp(a, b); double result = abs(temp); result /= N; return result; } complex<double> DFT(double in, int k) { double a = 0; double b = 0; int N = input.size(); for(int n = 0; n < N; n++) { a+= cos((2 * M_PI * k * n) / N) * input[n]; b+= -sin((2 * M_PI * k * n) / N) * input[n]; } complex<double> temp(a, b); return temp; } int main() { input.push_back(55); input.push_back(15); input.push_back(86); input.push_back(24); input.push_back(66); input.push_back(245); input.push_back(76); for(int k = 0; k < input.size(); k++) { output.push_back(DFT(input[k], k)); cout << "#" << k << ":\t" << input[k] << " \t>> abs: " << abs(output[k]) << " >> phase: " << arg(output[k]) << endl; } for(int n = 0; n < output.size(); n++) { result.push_back(IDFT(n)); cout << result[n] << endl; } return 0; }

If this is not just an educational exercise, why reinvent the wheel? — Ben Jones
– Ben Jones, Commented Aug 3, 2018 at 20:21
Maybe, if i decide to compile it for some microcontroller, i won't have problems with any additional libraries, and also, i must understand how it works in low level. — Roman Kulaha
– Roman Kulaha, Commented Aug 3, 2018 at 20:26
Strongly recommend sourceforge.net/projects/kissfft - small, self contained, fast. — mtrw
– mtrw, Commented Aug 3, 2018 at 20:31

Evg · Accepted Answer · 2018-08-03 21:17:05Z

Your inverse Fourier transform is obviously broken: you ignore the arguments of the complex numbers output[k].

It should look like this:

double IDFT(size_t n) { const auto ci = std::complex<double>(0, 1); std::complex<double> result; size_t N = output.size(); for (size_t k = 0; k < N; k++) result += std::exp((1. / N) * 2 * M_PI * k * n * ci) * output[k]; result /= N; return std::abs(result); }

Edit.

If you want to separate real and imaginary parts explicitly, you can use:

double IDFT(size_t n) { double a = 0; size_t N = output.size(); for (size_t k = 0; k < N; k++) { auto phase = (2 * M_PI * k * n) / N; a += cos(phase) * output[k].real() - sin(phase) * output[k].imag(); } a /= N; return a; }

In case the exp is confusing, this is equivalent to changing value to output[k].real() and output[k].imag() in OP's IDFT() accordingly.

Coral Kashri · Accepted Answer · 2020-03-07 11:45:35Z

For Intel cores computers:

There is a library of Intel- IPP. It offers many functions with a very high performance. It is very hard to write something that is faster than their functions, due to the vectors operations they use. Try it: https://software.intel.com/en-us/intel-ipp

https://software.intel.com/en-us/articles/how-to-use-intel-ipp-s-1d-fourier-transform-functions

Collectives™ on Stack Overflow

Discrete Fourier Transform C++

2 Answers 2

1 Comment

Comments

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

1 Comment

Comments

Related