I am trying to understand how to utilize the fft to process data I capture with a software defined radio. I discovered C++'s fftw3 library and attempted to look up some documentation or tutorials, and found there wasn't a whole ton to go off of.
I did however find this tutorial. The code to do the fft/ifft/etc. compiles and works just fine for me, but when I attempted to implement the fftShift function I get the following compiler error:
In file included from /usr/include/c++/9/algorithm:62,
from CppFFTW.cpp:13:
/usr/include/c++/9/bits/stl_algo.h: In instantiation of ‘_RandomAccessIterator std::_V2::__rotate(_RandomAccessIterator, _RandomAccessIterator, _RandomAccessIterator, std::random_access_iterator_tag) [with _RandomAccessIterator = double (*)[2]]’:
/usr/include/c++/9/bits/stl_algo.h:1449:27: required from ‘_FIter std::_V2::rotate(_FIter, _FIter, _FIter) [with _FIter = double (*)[2]]’
CppFFTW.cpp:76:54: required from here
/usr/include/c++/9/bits/stl_algo.h:1371:16: error: array must be initialized with a brace-enclosed initializer
1371 | _ValueType __t = _GLIBCXX_MOVE(*__p);
| ^~~
/usr/include/c++/9/bits/stl_algo.h:1373:22: error: invalid array assignment
1373 | *(__p + __n - 1) = _GLIBCXX_MOVE(__t);
| ^
/usr/include/c++/9/bits/stl_algo.h:1394:16: error: array must be initialized with a brace-enclosed initializer
1394 | _ValueType __t = _GLIBCXX_MOVE(*(__p + __n - 1));
| ^~~
/usr/include/c++/9/bits/stl_algo.h:1396:10: error: invalid array assignment
1396 | *__p = _GLIBCXX_MOVE(__t);
| ^
Here is the line I use to compile the code: g++ CppFFTW.cpp -lfftw3 -lm
Output for g++ --version:
'g++ (Ubuntu 9.4.0-1ubuntu1~20.04.1) 9.4.0'
This error message didn't make much sense to me, but I figured it meant that algorithm's std::rotate function didn't play nice with the fftw_complex data type anymore for some reason.
This left me with a few questions
What fftShift is supposed to do is shift the zero-components to the center like in this example
When # elements odd:
{a, b, c, d, e, f, g} -> {e, f, g, a, b, c, d}
When # elements even:
{a, b, c, d, e, f, g, h} -> {e, f, g, h, a, b, c, d}
Here is the fftShift definition definition: Edit: The original (broken) fftShift
// FFT shift for complex data
void fftShift(fftw_complex *data)
{
// even number of elements
if (N % 2 == 0) {
std::rotate(&data[0], &data[N >> 1], &data[N]);
}
// odd number of elements
else {
std::rotate(&data[0], &data[(N >> 1) + 1], &data[N]);
}
}
Here is the rest of the code (functions/calls not relevant to the question were omitted): Edit: Added #include and fixed fftShift functions as per Ted's contribution.
*
Example code for how to utilize the FFTW libs
Adapted from damian-dz C++ Tutorial
https://github.com/damian-dz/CppFFTW/tree/master/CppFFTW
(https://www.youtube.com/watch?v=geYbCA137PU&t=0s)
To compile, run: g++ CppFFTW.cpp -lfftw3 -lm
-lfftw3 -lm links the code to the fftw3 library
*/
#include <fftw3.h>
#include <iostream>
#include <cmath>
#include <algorithm>
#include <complex> // Added post feedback from Ted
// macros for real & imaginary parts
#define REAL 0
#define IMAG 1
// length of the complex arrays
#define N 8
/* Computres the 1-D Fast Fourier Transform. */
void fft(fftw_complex *in, fftw_complex *out)
{
// create a DFT plan
fftw_plan plan = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
// execute the plan
fftw_execute(plan);
// clean up
fftw_destroy_plan(plan);
fftw_cleanup();
}
/* Displays complex numbers in the form a +/- bi */
void displayComplex(fftw_complex *y)
{
for (int idx = 0; idx < N; ++idx) {
if (y[idx][IMAG] < 0) {
std::cout << y[idx][REAL] << " - " << abs(y[idx][IMAG]) << "i" << std::endl;
} else {
std::cout << y[idx][REAL] << " + " << abs(y[idx][IMAG]) << "i" << std::endl;
}
}
}
// Edit: Adjusted fftShift function per Ted's feedback
void fftShift(std::complex<double>* data) {
static_assert(sizeof(fftw_complex) == sizeof(std::complex<double>));
// even number of elements
if(N % 2 == 0) {
std::rotate(&data[0], &data[N >> 1], &data[N]);
}
// odd number of elements
else {
std::rotate(&data[0], &data[(N >> 1) + 1], &data[N]);
}
}
// Edit: Accompanying fftShift to re-cast data
void fftShift(fftw_complex* data) {
fftShift(reinterpret_cast<std::complex<double>*>(data));
}
// Original (broken) FFT shift for complex data
void fftShift(fftw_complex *data)
{
// even number of elements
if (N % 2 == 0) {
std::rotate(&data[0], &data[N >> 1], &data[N]);
}
// odd number of elements
else {
std::rotate(&data[0], &data[(N >> 1) + 1], &data[N]);
}
}
/* Test */
int main()
{
// input array
fftw_complex x[N];
// output array
fftw_complex y[N];
// fill the first array with some numbers
for (int idx = 0; idx < N; ++idx) {
x[idx][REAL] = idx + 1;
x[idx][IMAG] = 0;
}
// compute the FFT of x and store results in y
fft(x, y);
// display the results
std::cout << "FFT =" << std::endl;
displayComplex(y);
// "shifted" results
fftShift(y);
std::cout << "\nFFT shifted =" << std::endl;
displayComplex(y);
return 0;
}
Didn't see any obvious existing questions that I could see applying to my case (could be wrong, I am still learning C++), so I made an account and this is my first question here. Feel free to provide feedback so that I can make my question easier to understand.
For any object
zof typestd::complex<T>,reinterpret_cast<T(&)[2]>(z)[0]is the real part ofzandreinterpret_cast<T(&)[2]>(z)[1]is the imaginary part ofz.
The intent of this requirement is to preserve binary compatibility between the C++ library complex number types and the C language complex number types (and arrays thereof), which have an identical object representation requirement.
Now, fftw_complex isn't necessarily a complex number type defined in the C language - but you will likely get away with doing a similar cast. In fact, the comment in the fftw3.h header suggests it will most probably be ok:
/* If <complex.h> is included, use the C99 complex type. Otherwise
define a type bit-compatible with C99 complex */
#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I)
# define FFTW_DEFINE_COMPLEX(R, C) typedef R _Complex C
#else
# define FFTW_DEFINE_COMPLEX(R, C) typedef R C[2]
#endif
The C++ definition of fftw_complex becomes double[2] which is also why std::rotate fails - you can't assign arrays to arrays.
So I'd define my function taking a fftw_complex* like this:
void fftShift(std::complex<double>* data) {
static_assert(sizeof(fftw_complex) == sizeof(std::complex<double>));
// even number of elements
if(N % 2 == 0) {
std::rotate(&data[0], &data[N >> 1], &data[N]);
}
// odd number of elements
else {
std::rotate(&data[0], &data[(N >> 1) + 1], &data[N]);
}
}
void fftShift(fftw_complex* data) {
fftShift(reinterpret_cast<std::complex<double>*>(data));
}
int main() {
fftw_complex data[N];
fftShift(data);
}
A nicer alternative would be to use std::complex<double>s in your C++ code and only cast to fftw_complex* when calling fftw functions.
std::vector<std::complex<double>> wave;
fftw_function( reinterpret_cast<fftw_complex*>(wave.data()), wave.size() );
- How come this method worked in the tutorial video, but not for me?
It's because they use Visual Studio in the tutorial and the std::rotate implementation in that version of the standard library uses std::iter_swap that is actually capable of swapping whole arrays. If iter_swap is mandated by the standard to be able to do that, I don't know, but iter_swap in the g++, clang++, icx and MSVC implementations are all capable of doing that. However, g++, clang++ and icx do not use their own iter_swaps in their rotate implementations - at least not in this case.
We can borrow a rotate implementation from rotate @ cppreference that actually uses std::iter_swap and then all four can rotate your fftw_complexs - using their own iter_swaps. Demo
- What are some ways to go about understanding what this error message means?
It means that std::rotate tries to assign a double[2] to another double[2], as-if by doing:
double a[2]{};
double b[2]{};
a = b; // error: invalid array assignment
- How can I implement a working
fftShiftfunction?
See above.
- Is there a different compiler version I should be using?
No, even newer versions of g++ and clang++ would complain about the same issue. You could use Visual Studio or a different implementation of rotate (like the one on cppreference.com) - but I suggest just adapting and use std::complex<double>s that you cast to fftw_complexs when needed.