我正在使用 fftw 库实现一个 2D fft 示例,以了解如何使用它,并且我遇到了就地 r2c 案例的问题。
据我了解,我正在按照应该完成的方式进行操作:i)分配正确大小的内存(
N_x * 2 * (N_y/2 + 1)
,我从here获得)和ii)创建计划将输入作为 double* 传递,并将其作为输出传递,将其转换为 (fftw_complex*) 类型。然而,结果与我实现的异地 r2c 和就地 c2c 都不匹配。
运行这三种情况并进行比较的代码如下:
#include <fftw3.h>
#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;
int main(){
// define parameters
int N_x = 20;
int N_y = 20;
int N_y_complex = N_y/2 + 1;
double a = 2;
double b = 2;
double range_x = 10;
double range_y = 10;
double dx = range_x/(N_x);
double dy = range_y/(N_y);
////////////////////////////////////////// OUT-OF-PLACE R2C TRANSFORM //////////////////////////////////////////
// declare functions
double *in = (double*)fftw_malloc(N_x * N_y * sizeof(double));
fftw_complex *out = (fftw_complex*)fftw_malloc(N_x * N_y_complex * sizeof(fftw_complex));
for(int i = 0; i < N_x; i++){
double x = i*dx;
for(int j = 0; j < N_y; j++){
double y = j*dy;
if(fabs(x) <= a/2 && fabs(y) <= b/2){
in[i*N_y + j] = 1;
}
else if(fabs(range_x - x) <= a/2 && fabs(range_y - y) <= b/2){
in[i*N_y + j] = 1;
}
else if(fabs(x) <= a/2 && fabs(range_y - y) <= b/2){
in[i*N_y + j] = 1;
}
else if(fabs(range_x - x) <= a/2 && fabs(y) <= b/2){
in[i*N_y + j] = 1;
}
else{
in[i*N_y + j] = 0;
}
}
}
// Check if the function is even
bool is_even = true;
for (int i = 0; i < N_x / 2; i++) {
for (int j = 0; j < N_y; j++) {
int i_neg = (N_x - i) % N_x;
int j_neg = (N_y - j) % N_y;
if (in[i * N_y + j] != in[i_neg * N_y + j_neg]) {
cout << "Symmetry error: f(" << i * dx << "," << j * dy << ") = " << in[i * N_y + j]
<< " but f(" << i_neg * dx << "," << j_neg * dy << ") = " << in[i_neg * N_y + j_neg] << endl;
is_even = false;
cout << i << ", " << j << " || " << i_neg << ", " << j_neg << endl;
}
}
}
cout << "Checking f(x,y) parity: " << endl;
if (is_even) {
cout << " The function is even." << endl;
} else {
cout << " The function is not even." << endl;
}
fftw_plan p = fftw_plan_dft_r2c_2d(N_x,N_y,in,out,FFTW_ESTIMATE);
fftw_execute(p);
////////////////////////////////////////// IN-PLACE C2C FORWARD TRANSFORM //////////////////////////////////////////
fftw_complex *in2 = (fftw_complex*)fftw_malloc(N_x * N_y * sizeof(fftw_complex));
for(int i = 0; i < N_x; i++){
for(int j = 0; j < N_y; j++){
in2[i*N_y + j][0] = in[i*N_y + j];
in2[i*N_y + j][1] = 0;
}
}
fftw_plan p_inPlace_c2c = fftw_plan_dft_2d(N_x, N_y, in2, in2, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(p_inPlace_c2c);
bool success = false;
double tol = 1e-10;
for(int i = 0; i < N_x; i++){
for(int j = 0; j < N_y/2; j++){
if(abs(out[i * N_y_complex + j][0] - in2[i * N_y_complex + j][0]) > tol && abs(out[i * N_y_complex + j][1] - in2[i * N_y_complex + j][1]) > tol){
success = false;
goto skiploop;
}
else{
success = true;
}
}
}
skiploop:
cout << "Checking if in-place (forward C2C) and out-of-place results are matching: " << endl;
if(success == true){
cout << " Matching results for in-place c2c and out-of-place" << endl;
}
else{
cout << " Something went wrong" << endl;
}
////////////////////////////////////////// IN-PLACE R2C TRANSFORM //////////////////////////////////////////
// Calculate padded size: N_x * 2 * (N_y/2 + 1) = N_x * 2 * N_y_complex
int padded_size = N_x * 2 * N_y_complex;
// Allocate enough memory to hold complex data
double *in3 = (double*)fftw_malloc(padded_size * sizeof(double));
// Assign values for the real function to be transformed
for(int i = 0; i < N_x; i++){
for(int j = 0; j < N_y; j++){
in3[i*N_y + j] = in[i*N_y + j];
}
}
// Create plan casting in3 to fftw_complex type
fftw_plan p_inPlace_r2c = fftw_plan_dft_r2c_2d(N_x, N_y, in3, (fftw_complex*)in3, FFTW_ESTIMATE);
// Execute plan
fftw_execute(p_inPlace_r2c);
// printing first real value for debugging:
cout << "out-of-place r2c: " << out[0][0] << " || in-place c2c: " << in2[0][0] << " || in-place r2c: " << ((fftw_complex*)in3)[0][0] << endl;
// Check final results compared to out-of-place transform
for(int i = 0; i < N_x; i++){
for(int j = 0; j < N_y/2; j++){
if(abs(out[i * N_y_complex + j][0] - ((fftw_complex*)in3)[i * N_y_complex + j][0]) > tol && abs(out[i * N_y_complex + j][1] - ((fftw_complex*)in3)[i * N_y_complex + j][1]) > tol){
success = false;
goto skiploop2;
}
else{
success = true;
}
}
}
skiploop2:
cout << "Checking if in-place (R2C) and out-of-place results are matching: " << endl;
if(success == true){
cout << " Matching results for in-place r2c and out-of-place" << endl;
}
else{
cout << " Something went wrong" << endl;
}
fftw_free(in);
fftw_free(out);
fftw_free(in2);
fftw_free(in3);
fftw_destroy_plan(p);
fftw_destroy_plan(p_inPlace_c2c);
fftw_destroy_plan(p_inPlace_r2c);
}
我为获取所有三种情况显示的第一个实际值而输入的打印语句:
out-of-place r2c: 25 || in-place c2c: 25 || in-place r2c: -2.58976e+284
在这种情况下我错过了什么?我将其编译为
g++ -I/path/to/include -o fft2d fft2d.cpp -L/path/to/lib -lfftw3
,因为我正在使用的机器在自定义目录中有 fftw 库。
感谢@TimRoberts 为我澄清了这一点!我原以为额外的填充仅在用新的复数值覆盖数组时才重要,所以我只是在
in3
的末尾留下了所有额外的空间。
因此,在将 2D 数组表示为 1D 行主数组时,我需要修改这部分代码,以确保所有额外的填充都放置在与每行末尾相对应的正确位置:
// Assign values for the real function to be transformed
int extra_row_padding;
for(int i = 0; i < N_x; i++){
for(int j = 0; j < N_y; j++){
if(N_y % 2 == 0){
extra_row_padding = 2;
}
else{
extra_row_padding = 1;
}
in3[i*(N_y + extra_row_padding) + j] = in[i * N_y + j];
}
}
这是完整的工作示例:
#include <fftw3.h>
#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;
int main(){
// define parameters
int N_x = 20;
int N_y = 20;
int N_y_complex = N_y/2 + 1;
double a = 2;
double b = 2;
double range_x = 10;
double range_y = 10;
double dx = range_x/(N_x);
double dy = range_y/(N_y);
////////////////////////////////////////// OUT-OF-PLACE R2C TRANSFORM //////////////////////////////////////////
// declare functions
double *in = (double*)fftw_malloc(N_x * N_y * sizeof(double));
fftw_complex *out = (fftw_complex*)fftw_malloc(N_x * N_y_complex * sizeof(fftw_complex));
for(int i = 0; i < N_x; i++){
double x = i*dx;
for(int j = 0; j < N_y; j++){
double y = j*dy;
if(fabs(x) <= a/2 && fabs(y) <= b/2){
in[i*N_y + j] = 1;
}
else if(fabs(range_x - x) <= a/2 && fabs(range_y - y) <= b/2){
in[i*N_y + j] = 1;
}
else if(fabs(x) <= a/2 && fabs(range_y - y) <= b/2){
in[i*N_y + j] = 1;
}
else if(fabs(range_x - x) <= a/2 && fabs(y) <= b/2){
in[i*N_y + j] = 1;
}
else{
in[i*N_y + j] = 0;
}
}
}
// Check if the function is even
bool is_even = true;
for (int i = 0; i < N_x / 2; i++) {
for (int j = 0; j < N_y; j++) {
int i_neg = (N_x - i) % N_x;
int j_neg = (N_y - j) % N_y;
if (in[i * N_y + j] != in[i_neg * N_y + j_neg]) {
cout << "Symmetry error: f(" << i * dx << "," << j * dy << ") = " << in[i * N_y + j]
<< " but f(" << i_neg * dx << "," << j_neg * dy << ") = " << in[i_neg * N_y + j_neg] << endl;
is_even = false;
cout << i << ", " << j << " || " << i_neg << ", " << j_neg << endl;
}
}
}
cout << "Checking f(x,y) parity: " << endl;
if (is_even) {
cout << " The function is even." << endl;
} else {
cout << " The function is not even." << endl;
}
fftw_plan p = fftw_plan_dft_r2c_2d(N_x,N_y,in,out,FFTW_ESTIMATE);
fftw_execute(p);
////////////////////////////////////////// IN-PLACE C2C FORWARD TRANSFORM //////////////////////////////////////////
fftw_complex *in2 = (fftw_complex*)fftw_malloc(N_x * N_y * sizeof(fftw_complex));
for(int i = 0; i < N_x; i++){
for(int j = 0; j < N_y; j++){
in2[i*N_y + j][0] = in[i*N_y + j];
in2[i*N_y + j][1] = 0;
}
}
fftw_plan p_inPlace_c2c = fftw_plan_dft_2d(N_x, N_y, in2, in2, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(p_inPlace_c2c);
bool success = false;
double tol = 1e-10;
for(int i = 0; i < N_x; i++){
for(int j = 0; j < N_y/2; j++){
if(abs(out[i * N_y_complex + j][0] - in2[i * N_y_complex + j][0]) > tol && abs(out[i * N_y_complex + j][1] - in2[i * N_y_complex + j][1]) > tol){
success = false;
goto skiploop;
}
else{
success = true;
}
}
}
skiploop:
cout << "Checking if in-place (forward C2C) and out-of-place results are matching: " << endl;
if(success == true){
cout << " Matching results for in-place c2c and out-of-place" << endl;
}
else{
cout << " Something went wrong" << endl;
}
////////////////////////////////////////// IN-PLACE R2C TRANSFORM //////////////////////////////////////////
// Calculate padded size: N_x * 2 * (N_y/2 + 1) = N_x * 2 * N_y_complex
int padded_size = N_x * 2 * N_y_complex;
// Allocate enough memory to hold complex data
double *in3 = (double*)fftw_malloc(padded_size * sizeof(double));
// Assign values for the real function to be transformed
int extra_row_padding;
for(int i = 0; i < N_x; i++){
for(int j = 0; j < N_y; j++){
if(N_y % 2 == 0){
extra_row_padding = 2;
}
else{
extra_row_padding = 1;
}
in3[i*(N_y + extra_row_padding) + j] = in[i * N_y + j];
}
}
// Create plan casting in3 to fftw_complex type
fftw_plan p_inPlace_r2c = fftw_plan_dft_r2c_2d(N_x, N_y, in3, (fftw_complex*)in3, FFTW_ESTIMATE);
// Execute plan
fftw_execute(p_inPlace_r2c);
// printing first real value for debugging:
cout << "out-of-place r2c: " << out[0][0] << " || in-place c2c: " << in2[0][0] << " || in-place r2c: " << ((fftw_complex*)in3)[0][0] << endl;
// Check final results compared to out-of-place transform
for(int i = 0; i < N_x; i++){
for(int j = 0; j < N_y/2; j++){
if(abs(out[i * N_y_complex + j][0] - ((fftw_complex*)in3)[i * N_y_complex + j][0]) > tol && abs(out[i * N_y_complex + j][1] - ((fftw_complex*)in3)[i * N_y_complex + j][1]) > tol){
success = false;
goto skiploop2;
}
else{
success = true;
}
}
}
skiploop2:
cout << "Checking if in-place (R2C) and out-of-place results are matching: " << endl;
if(success == true){
cout << " Matching results for in-place r2c and out-of-place" << endl;
}
else{
cout << " Something went wrong" << endl;
}
fftw_free(in);
fftw_free(out);
fftw_free(in2);
fftw_free(in3);
fftw_destroy_plan(p);
fftw_destroy_plan(p_inPlace_c2c);
fftw_destroy_plan(p_inPlace_r2c);
}