请教FFTW3的应用
各位,我写了一段FFTW3做二维DFT的代码,但是结果不正确,请帮我找找原因何在?#include "stdafx.h"
#include <iostream>
#include <string>
#include <iomanip>
#include <fstream>
#include <conio.h>
#include <math.h>
#include <complex>
#include <fftw3.h>
#include <stdio.h>
#define M 64
#define N 64
#define Pi 3.14159265359
#define complex_d complex<double>
#define ELEM(r,c) (r*M+c)
double data[M][2*N];
using namespace std;
void dataread( ) //读取文件数据 经过检验读取的数据是正确的
{
using namespace std;
int a[2];
double b[N+1];
double c[1];
streampos weizhi;
streampos WZhi;
ifstream infile("ISAR2D_1.dat",ios::in);
if(!infile)
{
cerr<<"open error!"<<endl;
exit(1);
}
for(int i=0;i<2;i++)
{
infile>>a[i];
}
cout<<endl;
streampos weizhi1 =infile.tellg(); //weizhi1 是 64 64
infile.seekg(weizhi1);
for(int i=0;i<N+1;i++)
{
infile>>b[i];
}
streampos weizhi2 =infile.tellg(); //weizhi2 第二行结束
infile.seekg(weizhi2);
infile>>c[0];
streampos weizhi3=infile.tellg(); //weizhi3 第三行第二个数开始位置
weizhi=weizhi3-weizhi2; //weizhi 求出第一个数所占的空间
;
WZhi= weizhi3;
for(int i=0 ,k=1;i<M;i++,k++)
{
infile.seekg(WZhi);
for(int j=0;j<2*N;j++)
{
infile>>data[i][j];
}
streampos weizhi4=infile.tellg();
infile.seekg(weizhi4);
infile>>c[0];
//cout<<c[k]<<endl;
streampos weizhi5=infile.tellg();
WZhi=weizhi5;
}
//c[65]=0.0;
}
int showresult(fftw_complex* in, fftw_complex* out)
{
FILE* shibu = NULL;
fopen_s(&shibu, "shibu.txt", "w");
if(shibu == NULL)
{
throw "创建数据文件失败。";
}
for(int i=0;i<M;i++)
{
for(int j=0;j<N;j++)
{
fprintf_s(shibu, "%f\t",out[ELEM(i, j)][0]);
}
fprintf_s(shibu, "\n");
}
FILE* xubu = NULL;
fopen_s(&xubu, "xubu.txt", "w");
if(xubu == NULL)
{
throw "创建数据文件失败。";
}
for(int i=0;i<M;i++)
{
for(int j=0;j<N;j++)
{
fprintf_s(xubu, "%f\t", out[ELEM(i, j)][1]);
}
fprintf_s(xubu, "\n");
}
fclose(shibu);
fclose(xubu);
return 1;
}
int main()
{
fftw_complex *in, *out;
fftw_plan p;
int i, j;
//complex_d G[M][N];
//读取数据
dataread( );
// 分配存储空间
in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * M* N);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * M *N);
// 设置变换计划
p = fftw_plan_dft_2d(M, N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
//加窗函数
for (i=0; i<M; i++)
{
for (j=0; j<2*N; j=j+2)
{
data[i][j]=data[i][j];//*(0.5-0.5*cos(2*Pi*28/64));
}
}
// 设置测试数据
for (i=0; i<M; i++)
{
int k=0;
for (j=0; j<2*N; j=j+2)
{
in[ELEM(i,k)][0] =data[i][j];
in[ELEM(i,k)][1] =data[i][j+1];
//cout<<in[ELEM(i, k)][0]<<endl;
k=k+1;
}
}
// 执行变换
fftw_execute(p);
// 显示测试结果
showresult(in, out);
// 释放内存
fftw_destroy_plan(p);
fftw_free(in);
fftw_free(out);
return 1;
}