fftw3反变换不起作用 [英] fftw3 inverse transform not work

查看:329
本文介绍了fftw3反变换不起作用的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我在写一个简单的code复杂的复杂DFT在C与fftw3库。
我写的输入数组双数据的文件,所以我可以用MATLAB的FFT功能进行比较。
我试图执行逆变换,从变换数组,但结果和第一个输入数组是不同的。这是我的结果:

  FFTW3变换的欢迎及LT;<<<<输入采样的数量(整数)N(位:64)($ P $ 2 pferably功率):8输入样例
在[0] [0] =在-216448918.015237 [0] [1] = 0.000000
在[1] [0] =在948904790.062151 [1] [1] = 0.000000
在[2] [0] =在826811206.185300 [2] [1] = 0.000000
在[3] [0] =在1868763250.342451 [3] [1] = 0.000000
在[4] [0] =在703135606.077152 [4] [1] = 0.000000
在[5] [0] =在-1989016445.622210 [5] [1] = 0.000000
在[6] [0] =在1912963650.704585 [6] [1] = 0.000000
在[7] [0] =在811527262.805480 [7] [1] = 0.000000 按回车键继续......
 正向变换系数出[0] [0] = 4866640402.539672出[0] [1] = 0.000000
出[1] [0] = 410260768.150135出[1] [1] = -1738850319.926936
出[2] [0] = -2253088168.827970出[2] [1] = 3720402168.707990
出[3] [0] = -2249429816.334913出[3] [1] = -3911155208.965507
出[4] [0] = 1586282687.363928出[4] [1] = 0.000000
出[5] [0] = -2249429816.334913出[5] [1] = 3911155208.965507
出[6] [0] = -2253088168.827970出[6] [1] = -3720402168.707990
出[7] [0] = 410260768.150135出[7] [1] = 1738850319.926936
你要计算逆变换? (Y / N)
ÿ
逆变换系数
转[0] [0] = -1731591344.121896转[0] [1] = 0.000000
转[1] [0] = 7591238320.497208转[1] [1] = 0.000000
转[2] [0] = 6614489649.482399转[2] [1] = 0.000000
转[3] [0] = 14950106002.739609转[3] [1] = 0.000000
转[4] [0] = 5625084848.617215转[4] [1] = 0.000000
转[5] [0] = -15912131564.977680转[5] [1] = 0.000000
转[6] [0] = 15303709205.636681转[6] [1] = 0.000000
转[7] [0] = 6492218102.443840转[7] [1] = 0.000000

正如你看到'在'和'转'数组是不同的,但直接的变换是正确的。我用MATLAB进行了比较和结果都是一样的。
当我执行使用MATLAB逆变换我获得输入数组。
我该怎么办?

这是我的C code:

 的#include< fftw3.h>
#包括LT&;&math.h中GT;
#包括LT&;&stdio.h中GT;
#包括LT&;&stdlib.h中GT;
#定义PI 3.141592653589诠释的main()
{  fftw_complex *中,*出来,*转;
  INT I,F0,A,N;
  没有烧焦;
  FILE * FP;
  fftw_plan磷;  的printf(\\ n \\ n>>>>> FFTW3变换的欢迎及LT;<<<<);
  的printf(\\ n \\ n输入数字(整数)样品的N(位:%LD)($ P $ 2 pferably功率):(的sizeof(fftw_complex)/ 2)* 8);
  scanf函数(%d个,&安培; N);  // F0 = 50;
  // A = 1;
  //分配输入和放大器记忆;输出数组
  在=(fftw_complex *)fftw_malloc(的sizeof(fftw_complex)* N);
  OUT =(fftw_complex *)fftw_malloc(的sizeof(fftw_complex)* N);
  REV =(fftw_complex *)fftw_malloc(的sizeof(fftw_complex)* N);  //打开数据文件
  如果((FP = FOPEN(lista_numeri_double.dat,RB))== NULL)
  {
    的printf(\\ n错误读取档案\\ n);
    出口(1);
  }
  的printf(\\输入NSample个);
  //分配样品从文件中读取
  对于(i = 0; I< N;我++)
  {
    //在[I] [0] = A * COS(2 * PI * F0 * I / N);
    FREAD(安培;在[I] [0]的sizeof(双),1,FP);
    在[I] [1] = 0;    的printf(\\ n在[%D]。[0] =%F \\ t \\锡[%D]。[1] =%F,我在[我] [0],我在[我] [1] );
  }
  //计划和执行变换
  P = fftw_plan_dft_1d(N,IN,OUT,FFTW_FORWARD,FFTW_ESTIMATE);
  fftw_execute(P);  的printf(\\ n \\ n按Enter继续... \\ n);
  scanf函数(%C,&安培;无);
  //打印输出值
  的printf(\\ n \\ nFORWARD变换系数\\ n);
  对于(i = 0; I< N;我++)
  {
    的printf(\\ NOUT内容[%d] [0] =%F \\ t \\兜售[%D]。[1] =%F,我出[I] [0],我,出[I] [1] );
  }  fftw_destroy_plan(P);  的printf(\\ n您要计算逆变换(Y / N)\\ n吗?);
  scanf函数(%C,&安培;无);  如果(没有=='Y')
  {    //计划并执行反变换
    P = fftw_plan_dft_1d(N,出,转,FFTW_BACKWARD,FFTW_ESTIMATE);
    fftw_execute(P);    的printf(\\ n \\ nINVERSE变换系数\\ n);
    对于(i = 0; I< N;我++)
    {
      的printf(REV [%D]。[0] =%F \\ t \\崔佛[%D]。[1] =%F \\ N,我,转速[我] [0],我,转速[我] [1 ]);
    }    fftw_destroy_plan(P);
  }  返回0;
}


解决方案

Matlab和FFTW之间的差异自带应用的缩放因子转换。

Matlab的FFT 是标准化的,使用FFTW的算法在FFTW的文档描述,不归。换句话说,全圆周的一个因素变换使用FFTW(向前后向后)缩放结果 N

相应地,在比较阵列显示,由8(大小 N 在你的例子中,变换使用)。一致的因子按比例

I'm writing a simple code for complex to complex DFT in c with fftw3 library. i have written a file with input array double data so I can compare with matlab fft function. I try to execute the reverse transform from transform array but results and first input array are different. this is my results:

FFTW3 TRANSFORM WELCOME <<<<<

enter the number (integer) N of samples (Bit: 64) (preferably power of 2):8

SAMPLE INPUT
in[0][0] = -216448918.015237        in[0][1] = 0.000000 
in[1][0] = 948904790.062151         in[1][1] = 0.000000
in[2][0] = 826811206.185300         in[2][1] = 0.000000
in[3][0] = 1868763250.342451        in[3][1] = 0.000000
in[4][0] = 703135606.077152         in[4][1] = 0.000000
in[5][0] = -1989016445.622210       in[5][1] = 0.000000
in[6][0] = 1912963650.704585        in[6][1] = 0.000000
in[7][0] = 811527262.805480         in[7][1] = 0.000000

 Hit enter to continue ... 


 FORWARD TRANSFORM COEFFICIENTS

out[0][0] = 4866640402.539672       out[0][1] = 0.000000
out[1][0] = 410260768.150135        out[1][1] = -1738850319.926936
out[2][0] = -2253088168.827970      out[2][1] = 3720402168.707990
out[3][0] = -2249429816.334913      out[3][1] = -3911155208.965507
out[4][0] = 1586282687.363928       out[4][1] = 0.000000
out[5][0] = -2249429816.334913      out[5][1] = 3911155208.965507
out[6][0] = -2253088168.827970      out[6][1] = -3720402168.707990
out[7][0] = 410260768.150135        out[7][1] = 1738850319.926936
do you want to calculate the inverse-transform? (y/n) 
y


INVERSE TRANSFORM COEFFICIENTS
rev[0][0] = -1731591344.121896      rev[0][1] = 0.000000
rev[1][0] = 7591238320.497208       rev[1][1] = 0.000000
rev[2][0] = 6614489649.482399       rev[2][1] = 0.000000
rev[3][0] = 14950106002.739609      rev[3][1] = 0.000000
rev[4][0] = 5625084848.617215       rev[4][1] = 0.000000
rev[5][0] = -15912131564.977680        rev[5][1] = 0.000000
rev[6][0] = 15303709205.636681      rev[6][1] = 0.000000
rev[7][0] = 6492218102.443840       rev[7][1] = 0.000000

As you see 'in' and 'rev' arrays are different but direct transform is correct. I've compared it with matlab and results are the same. When I execute the inverse transform with matlab I obtain the input array. What can I do?

this is my c code:

#include <fftw3.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define PI 3.141592653589

int main()
{

  fftw_complex *in, *out, *rev;
  int i,f0,A,N;
  char no;  
  FILE *fp;
  fftw_plan p;

  printf("\n\n>>>>> FFTW3 TRANSFORM WELCOME <<<<<");
  printf("\n\n enter the number (integer) N of samples (bit: %ld) (preferably power of 2):",(sizeof(fftw_complex)/2)*8);
  scanf("%d",&N);

  //f0 = 50;
  //A = 1;


  //allocating memory for input & output arrays
  in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*N);
  out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*N);
  rev = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*N);



  //Opening the data file
  if((fp=fopen("lista_numeri_double.dat","rb"))==NULL)
  {
    printf("\nError reading file\n");
    exit(1);
  }


  printf("\nSAMPLE INPUT");
  //assigning samples read from the file
  for(i=0;i<N;i++)
  {
    //in[i][0] = A * cos(2*PI*f0*i/N);
    fread(&in[i][0],sizeof(double),1,fp);   
    in[i][1]=0;

    printf("\nin[%d][0] = %f \t\tin[%d][1] = %f",i,in[i][0],i,in[i][1]);
  }


  //plan and execute transform
  p = fftw_plan_dft_1d(N,in,out,FFTW_FORWARD,FFTW_ESTIMATE);
  fftw_execute(p);

  printf("\n\n Hit enter to continue ... \n");
  scanf("%c",&no);


  //print output values 
  printf("\n\nFORWARD TRANSFORM COEFFICIENTS\n");
  for(i=0;i<N;i++)
  {
    printf("\nout[%d][0] = %f \t\tout[%d][1] = %f",i,out[i][0],i,out[i][1]);
  }

  fftw_destroy_plan(p); 

  printf("\n do you want to calculate the inverse-transform? (y/n)  \n");
  scanf ("%c",&no);

  if(no=='y')
  {

    //plan and execute inverse transform
    p = fftw_plan_dft_1d(N,out,rev,FFTW_BACKWARD,FFTW_ESTIMATE);
    fftw_execute(p);

    printf("\n\nINVERSE TRANSFORM COEFFICIENTS\n");
    for(i=0;i<N;i++)
    {
      printf("rev[%d][0] = %f \t\trev[%d][1] = %f\n",i,rev[i][0],i,rev[i][1]);
    }

    fftw_destroy_plan(p);           
  }

  return 0;
}

解决方案

The difference between Matlab and FFTW comes with the scaling factor applied to the transform.

Whereas Matlab's FFT is normalized, the algorithm used by FFTW as described in FFTW's documentation, is not normalized. In other words, the full-circle transform using FFTW (forward followed by backward) scales the result by a factor N.

Correspondingly, comparing the in and rev array shows that rev is scaled by consistent factor of 8 (the size N of the transform used in your example).

这篇关于fftw3反变换不起作用的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆