首页 诗词 字典 板报 句子 名言 友答 励志 学校 网站地图
当前位置: 首页 > 教程频道 > 开发语言 > C++ >

FFT赋值有关问题,拜请高手来解决

2012-03-24 
求助FFT赋值问题,拜请高手来解决我参与一个项目,涉及到FFT算法(快速傅里叶变换算法),参考网上的代码,如下:

求助FFT赋值问题,拜请高手来解决
我参与一个项目,涉及到FFT算法(快速傅里叶变换算法),参考网上的代码,如下:
void CFFT::FFT(float xreal [], float ximag [], int n)
{
  // 快速傅立叶变换,将复数 x 变换后仍保存在 x 中,xreal, ximag 分别是 x 的实部和虚部
  float wreal [N / 2], wimag [N / 2], treal, timag, ureal, uimag, arg;
  int m, k, j, t, index1, index2;
   
  bitrp (xreal, ximag, n);//位反转置换

  // 计算 1 的前 n / 2 个 n 次方根的共轭复数 W'j = wreal [j] + i * wimag [j] , j = 0, 1, ... , n / 2 - 1
  arg = - 2 * PI / n;
  treal = cos (arg);
  timag = sin (arg);
  wreal [0] = 1.0;
  wimag [0] = 0.0;
  for (j = 1; j < n/2 ; j ++)
  {
  wreal [j] = wreal [j - 1] * treal - wimag [j - 1] * timag;
  wimag [j] = wreal [j - 1] * timag + wimag [j - 1] * treal;
  }

  for (m = 2; m <= n; m *= 2)
  {
  for (k = 0; k < n; k += m)
  {
  for (j = 0; j < m / 2; j ++)
  {
  index1 = k + j;
  index2 = index1 + m / 2;
  t = n * j / m; // 旋转因子 w 的实部在 wreal [] 中的下标为 t
  treal = wreal [t] * xreal [index2] - wimag [t] * ximag [index2];
  timag = wreal [t] * ximag [index2] + wimag [t] * xreal [index2];
  ureal = xreal [index1];
  uimag = ximag [index1];
  xreal [index1] = ureal + treal;
  ximag [index1] = uimag + timag;
  xreal [index2] = ureal - treal;
  ximag [index2] = uimag - timag;
  }
  }
  }
}
经过与MATLAB的测试,该算法结果正确无误,与MATLAB一致。
但是现在的问题出现在赋值给算法的过程中,原来的程序,是从txt文本里把数值转到数组再执行FFT函数,代码如下:
void CFFT::FFT_test1 ()
{
  char inputfile [] = "F:\\txtdata\\input.txt"; // 从文件 input.txt 中读入原始数据
  char outputfile [] = "F:\\txtdata\\output.txt"; // 将结果输出到文件 output.txt 中
  float xreal [N/2] = {}, ximag [N/2] = {};
  int n, i;
  FILE *input, *output;
//测试文本是否能够找到
  if (!(input = fopen (inputfile, "r")))
  {
printf ("Cannot open file. ");
exit (1);
  }
  if (!(output = fopen (outputfile, "w")))
  {
printf ("Cannot open file. ");
exit (1);
  }
   
  i = 0;
  while ((fscanf (input, "%f%f", xreal + i, ximag + i)) != EOF)
  {
i ++;
  }
  n = i; // 要求 n 为 2 的整数幂
  while (i > 1)
  {
  if (i % 2)
  {
  fprintf (output, "%d is not a power of 2! ", n);
  exit (1);
  }
  i /= 2;
  }
  FFT (xreal, ximag, n);
  fprintf (output, "FFT: i real imag ");
  for (i = 0; i < n; i ++)
  {
  fprintf (output, "%4d %8.4f %8.4f ", i, xreal [i], ximag [i]);
  }
  fprintf (output, "================================= ");

  IFFT (xreal, ximag, n);
  fprintf (output, "IFFT: i real imag ");
  for (i = 0; i < n; i ++)
  {
  fprintf (output, "%4d %8.4f %8.4f ", i, xreal [i], ximag [i]);
  }

  if ( fclose (input))
  {
  printf ("File close error. ");
  exit (1);


  }
  if ( fclose (output))
  {
  printf ("File close error. ");
  exit (1);
  }
}
这段代码经过测试,没有问题,赋值正确,得的结果也正确。
我实际需求,需要从SQL数据库取出数据,经过FFT运算以后,把值(复数的实部、虚部和模)再存到SQL数据库里。我写的代码如下:
void CFFT::FFT_test (LPCTSTR c,LPCTSTR d)
{
//CString strAllData="-1.11131700 0 -1.23425700 0 -0.86663300 0 -0.54973000 0";
char input[20480];
int i,n,r;
float xreal[N/2],ximag[N/2];
char *p;
float modulo;
sprintf(input, "%s" ,c);
p = input;
i = 0;
while (1)
{
if ( p== NULL) 
break;
r=sscanf( p, "%f%f%n", xreal + i, ximag + i,&n);
//if ( EOF== r) 
//break;
if ( 2== r) 
{
i++;
if (i >= N/2-1) 
break;
p+=n;
}
}
n = i;
printf("n=%d\n",n); // 要求 n 为 2 的整数幂
for (i=0;i <= n;i++) 
{
printf("xreal[%d]==%lg,ximag[%d]==%lg\n",i,xreal[i],i,ximag[i]);
}
//执行傅里叶变换
FFT (xreal, ximag, n);
///////////////以下代码为把数据写入数据库/////////////////////
//连接数据库
ADOConn m_AdoConn;
m_AdoConn.OnInitADOConn();
//设置INSERT语句
for (i=0;i <= n;i++) 
{
printf("xreal[%d]==%lg,ximag[%d]==%lg\n",i,xreal[i],i,ximag[i]);
modulo=sqrt(xreal[i]*xreal[i] + ximag[i]*ximag[i]);
CString strSQL;
strSQL.Format( "INSERT INTO tb_pdatavalue(p_no, p_type, p_value_s, p_value_x,p_value) VALUES ('%s', 'Channel 1 Timewave Amplitude', '%f', '%f','%f')",d, xreal[i], ximag[i],modulo);
m_AdoConn.ExecuteSQL( (_bstr_t)strSQL );
}
//断开与数据库的连接
m_AdoConn.ExitConnect();
//执行傅里叶逆变换
  //IFFT (xreal, ximag, n);
}
我的问题是,我通过测试,传值那里貌似没有什么问题,可是运算的结果却不正确。一共传了512个数进去,前面256个数运算后的结果(比如实部)是正确结果的两倍左右(不完全等于2倍,但是只差了一点),而后面的256个数则非常大(已经达到几亿大了),请高手指点我,谢谢了!!!!在线急等!!!!!

[解决办法]

C/C++ code
FILE *fp = fopen( "output.txt", "w" );for( int i = 0; i < n; i++ )    fprintf( fp, "%10f %10f\n", xreal[i], ximag[i] );fclose( fp ); 

热点排行