我试图了解如何利用 fft 来处理我用软件定义的无线电捕获的数据。我发现了 C++ 的 fftw3 库,并尝试查找一些文档或教程,发现没有太多可以参考的内容。
\n不过我确实找到了这个教程。执行 fft/ifft/等的代码。对我来说编译和工作得很好,但是当我尝试实现 fftShift 函数时,我收到以下编译器错误:
\nIn file included from /usr/include/c++/9/algorithm:62,\n from CppFFTW.cpp:13:\n/usr/include/c++/9/bits/stl_algo.h: In instantiation of \xe2\x80\x98_RandomAccessIterator std::_V2::__rotate(_RandomAccessIterator, _RandomAccessIterator, _RandomAccessIterator, std::random_access_iterator_tag) [with _RandomAccessIterator = double (*)[2]]\xe2\x80\x99:\n/usr/include/c++/9/bits/stl_algo.h:1449:27: required from \xe2\x80\x98_FIter std::_V2::rotate(_FIter, _FIter, _FIter) [with _FIter = double (*)[2]]\xe2\x80\x99\nCppFFTW.cpp:76:54: required from here\n/usr/include/c++/9/bits/stl_algo.h:1371:16: error: array must be initialized with a brace-enclosed initializer\n 1371 | _ValueType __t = _GLIBCXX_MOVE(*__p);\n | ^~~\n/usr/include/c++/9/bits/stl_algo.h:1373:22: error: invalid array assignment\n 1373 | *(__p + __n - 1) = _GLIBCXX_MOVE(__t);\n | ^\n/usr/include/c++/9/bits/stl_algo.h:1394:16: error: array must be initialized with a brace-enclosed initializer\n 1394 | _ValueType __t = _GLIBCXX_MOVE(*(__p + __n - 1));\n | ^~~\n/usr/include/c++/9/bits/stl_algo.h:1396:10: error: invalid array assignment\n 1396 | *__p = _GLIBCXX_MOVE(__t);\n | ^\nRun Code Online (Sandbox Code Playgroud)\n这是我用来编译代码的行:g++ CppFFTW.cpp -lfftw3 -lm
\ng++ --version 的输出:\n\'g++ (Ubuntu 9.4.0-1ubuntu1~20.04.1) 9.4.0\'
这个错误消息对我来说没有多大意义,但我认为这意味着算法的 std::rotate 函数由于某种原因不再与fftw_complex数据类型配合得很好。
\n这给我留下了几个问题
\nfftShift 应该做的是将零分量移动到中心,如本例所示
\nWhen # elements odd:\n{a, b, c, d, e, f, g} -> {e, f, g, a, b, c, d} \nRun Code Online (Sandbox Code Playgroud)\nWhen # elements even:\n{a, b, c, d, e, f, g, h} -> {e, f, g, h, a, b, c, d}\nRun Code Online (Sandbox Code Playgroud)\n这是 fftShift 定义:\n编辑:原始(损坏的)fftShift
\n// FFT shift for complex data \nvoid fftShift(fftw_complex *data)\n{\n // even number of elements\n if (N % 2 == 0) {\n std::rotate(&data[0], &data[N >> 1], &data[N]);\n }\n // odd number of elements\n else {\n std::rotate(&data[0], &data[(N >> 1) + 1], &data[N]);\n }\n}\nRun Code Online (Sandbox Code Playgroud)\n这是代码的其余部分(省略了与问题无关的函数/调用):\n编辑:根据 Ted\ 的贡献添加了 #include 和固定的 fftShift 函数。
\n* \n Example code for how to utilize the FFTW libs\n Adapted from damian-dz C++ Tutorial\n https://github.com/damian-dz/CppFFTW/tree/master/CppFFTW\n (https://www.youtube.com/watch?v=geYbCA137PU&t=0s)\n To compile, run: g++ CppFFTW.cpp -lfftw3 -lm\n -lfftw3 -lm links the code to the fftw3 library\n*/\n\n#include <fftw3.h>\n#include <iostream>\n#include <cmath>\n#include <algorithm>\n#include <complex> // Added post feedback from Ted\n\n// macros for real & imaginary parts\n#define REAL 0\n#define IMAG 1\n\n// length of the complex arrays\n#define N 8\n\n/* Computres the 1-D Fast Fourier Transform. */\nvoid fft(fftw_complex *in, fftw_complex *out)\n{\n // create a DFT plan\n fftw_plan plan = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);\n // execute the plan\n fftw_execute(plan);\n // clean up\n fftw_destroy_plan(plan);\n fftw_cleanup();\n}\n\n/* Displays complex numbers in the form a +/- bi */\nvoid displayComplex(fftw_complex *y)\n{\n for (int idx = 0; idx < N; ++idx) {\n if (y[idx][IMAG] < 0) {\n std::cout << y[idx][REAL] << " - " << abs(y[idx][IMAG]) << "i" << std::endl;\n } else {\n std::cout << y[idx][REAL] << " + " << abs(y[idx][IMAG]) << "i" << std::endl;\n }\n }\n}\n\n// Edit: Adjusted fftShift function per Ted\'s feedback\nvoid fftShift(std::complex<double>* data) {\n static_assert(sizeof(fftw_complex) == sizeof(std::complex<double>));\n\n // even number of elements\n if(N % 2 == 0) {\n std::rotate(&data[0], &data[N >> 1], &data[N]);\n }\n // odd number of elements\n else {\n std::rotate(&data[0], &data[(N >> 1) + 1], &data[N]);\n }\n}\n// Edit: Accompanying fftShift to re-cast data\nvoid fftShift(fftw_complex* data) {\n fftShift(reinterpret_cast<std::complex<double>*>(data));\n}\n\n// Original (broken) FFT shift for complex data \nvoid fftShift(fftw_complex *data)\n{\n // even number of elements\n if (N % 2 == 0) {\n std::rotate(&data[0], &data[N >> 1], &data[N]);\n }\n // odd number of elements\n else {\n std::rotate(&data[0], &data[(N >> 1) + 1], &data[N]);\n }\n}\n\n/* Test */\nint main()\n{\n // input array\n fftw_complex x[N];\n\n // output array\n fftw_complex y[N];\n\n // fill the first array with some numbers\n for (int idx = 0; idx < N; ++idx) {\n x[idx][REAL] = idx + 1;\n x[idx][IMAG] = 0;\n }\n \n // compute the FFT of x and store results in y\n fft(x, y);\n // display the results\n std::cout << "FFT =" << std::endl;\n displayComplex(y);\n\n // "shifted" results\n fftShift(y);\n std::cout << "\\nFFT shifted =" << std::endl;\n displayComplex(y);\n \n return 0;\n}\n\nRun Code Online (Sandbox Code Playgroud)\n没有看到任何明显的现有问题适用于我的案例(可能是错误的,我仍在学习 C++),所以我注册了一个帐户,这是我的第一个问题。请随时提供反馈,以便我可以使我的问题更容易理解。
\n对于任何
z类型的对象std::complex<T>,reinterpret_cast<T(&)[2]>(z)[0]是 的实部z,reinterpret_cast<T(&)[2]>(z)[1]是 的虚部z。
此要求的目的是保持 C++ 库复数类型和 C 语言复数类型(及其数组)之间的二进制兼容性,它们具有相同的对象表示要求。
现在,fftw_complex不一定是 C 语言中定义的复数类型 - 但您可能会通过执行类似的转换而逃脱惩罚。事实上,fftw3.h标题中的评论表明它很可能没问题:
/* 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
Run Code Online (Sandbox Code Playgroud)
C++ 的定义fftw_complex变成了double[2]这也是std::rotate失败的原因——你不能将数组分配给数组。
所以我会像这样定义我的函数fftw_complex*:
/* 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
Run Code Online (Sandbox Code Playgroud)
更好的替代方法是std::complex<double>在 C++ 代码中使用 s,并且仅fftw_complex*在调用fftw函数时进行强制转换。
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);
}
Run Code Online (Sandbox Code Playgroud)
- 为什么这个方法在教程视频中有效,但对我不起作用?
这是因为他们在教程中使用 Visual Studio,并且std::rotate该版本的标准库中的实现std::iter_swap实际上能够交换整个数组。我不知道标准是否iter_swap强制要求能够做到这一点,但iter_swap在 g++、clang++、icx 和 MSVC 实现中都能够做到这一点。iter_swap但是,g++、clang++ 和 icx在其实现中不使用自己的s rotate- 至少在本例中不使用。
我们可以从实际使用的rotate@cppreferencerotate借用一个实现,然后所有四个都可以使用它们自己的s 。演示std::iter_swaprotatefftw_complexiter_swap
- 有哪些方法可以帮助您理解此错误消息的含义?
这意味着std::rotate尝试将 a 分配double[2]给另一个double[2],就好像通过执行以下操作:
std::vector<std::complex<double>> wave;
fftw_function( reinterpret_cast<fftw_complex*>(wave.data()), wave.size() );
Run Code Online (Sandbox Code Playgroud)
- 如何实现工作
fftShift功能?
往上看。
- 我应该使用不同的编译器版本吗?
不,即使是较新版本的g++和clang++也会抱怨同样的问题。您可以使用 Visual Studio 或其他实现rotate(例如 cppreference.com 上的实现) - 但我建议在需要时调整并使用std::complex<double>您转换为fftw_complexs 的 s 。