请问：这个复数FFT程序对吗？怎么验证？

只看该作者 · 2009-9-28 22:30

本帖最后由 void_c 于 2009-9-29 09:13 编辑

/*0001*/
/*0002*/  #include "util_complex.h"
/*0003*/
/*0004*/
/*0005*/  #include <math.h>
/*0006*/  #include "util_complex.h"
/*0007*/
/*0008*/  #define  PI    3.141592653589793238462643
/*0009*/
/*0010*/  #define FFT_N 8
/*0011*/  #define FFT_M 3
/*0012*/
/*0013*/  void fft_rev(complex_t f[])
/*0014*/  {
/*0015*/    complex_t temp;
/*0016*/    unsigned short ii,jj,kk,ll;
/*0017*/
/*0018*/    ll = FFT_N/2;
/*0019*/    jj = 0;
/*0020*/    for (ii=0;ii<FFT_N-1;ii++)
/*0021*/    {
/*0022*/       if (ii < jj)
/*0023*/       {
/*0024*/          temp=f[ii];
/*0025*/          f[ii]=f[jj];
/*0026*/          f[jj]=temp;
/*0027*/       }
/*0028*/       kk = ll;
/*0029*/       while (kk <= jj)
/*0030*/       {
/*0031*/          jj -= kk;
/*0032*/          kk >>= 1;
/*0033*/       }
/*0034*/       jj += kk;
/*0035*/    }
/*0036*/  }
/*0037*/
/*0038*/  void fft_c(complex_t f[])
/*0039*/  {
/*0040*/    unsigned short order;             //阶
/*0041*/    unsigned short group,group_count; //组
/*0042*/    unsigned short fly, fly_count;    //蝶形单元
/*0043*/    float arg;                         //旋转角度
/*0044*/    complex_t w;                   //旋转因子
/*0045*/
/*0046*/    complex_t p,q,r;             //临时变量，复数
/*0047*/
/*0048*/    unsigned short fly_i,fly_j;          //蝶形运算单元操作数序号
/*0049*/
/*0050*/          fft_rev(f);                            //倒序
/*0051*/
/*0052*/    for (order=0;order<FFT_M;order++)       //阶数，从0开始
/*0053*/    {
/*0054*/          group_count=FFT_N>>(order+1);    //第order阶组数, 第0阶组数FFT_N/2 ,第一阶组数 FFT_N/4，
/*0055*/
/*0056*/          for (group=0; group<group_count; group++)  //第order阶第group组
/*0057*/          {
/*0058*/                fly_count=FFT_N/group_count/2; //第order阶第group组蝶形单元个数，第0阶某组蝶形单元个数1
/*0059*/
/*0060*/             for (fly=0;fly<fly_count;fly++) //第order阶第group组第fly蝶形单元
/*0061*/             {
/*0062*/                fly_i=group*fly_count*2 +fly;    //第order阶第group组第fly蝶形单元,操作数序号fly_i
/*0063*/                fly_j=fly_i+fly_count;    //第order阶第group组第fly蝶形单元,操作数序号fly_j
/*0064*/
/*0065*/                arg = (float)fly/fly_count*PI; //第order阶第group组第fly蝶形单元,旋转角度
/*0066*/                arg = -(float)fly/fly_count*PI; //第order阶第group组第fly蝶形单元,旋转角度
/*0067*/                这里倒是是正还是负数，为什么？
/*0068*/
/*0069*/                w=(complex_t){cos(arg),sin(arg)}; //第order阶第group组第fly蝶形单元,旋转因子
/*0070*/
/*0071*/                p=f[fly_i];                //蝶形运算先保存操作数
/*0072*/                q=f[fly_j];                //蝶形运算先保存操作数
/*0073*/                                  r=complex_mul(q,w);                //中间变量
/*0074*/                f[fly_i]=complex_add(p,r);          //蝶形运算
/*0075*/                f[fly_j]=complex_sub(p,r);    //蝶形运算
/*0076*/             }
/*0077*/          }
/*0078*/    }
/*0079*/  }
/*0080*/
/*0081*/
/*0082*/

只看该作者 · 2009-9-28 23:57

试了下8点FFT，

输入：
(1,0) (2,0) (3,0) (4,0) (5,0) (6,0) (7,0) (8,0)

输出:
(36,0) (-4,-9.65685) (-4,-4) (-4,-1.65685) (-4,0) (-4,1.65685) (-4,4) (-4,9.65685)

http://zhidao.baidu.com/question/1246069.html
网上有人8点FFT结果是：
(36,0) (-4,9.65685) (-4,4) (-4,1.65685) (-4,0) (-4,-1.65685) (-4,-4) (-4,-9.65685)

和我的程序徐步完全相反，

程序如果把旋转角度相反，结果及一样了。
arg = (float)fly/fly_count*PI;
改成:
arg = -(float)fly/fly_count*PI;

为什么呢？

请问：这个复数FFT程序对吗？ 怎么验证？

相关下载

相关帖子

请问：这个复数FFT程序对吗？怎么验证？