论坛 › 算法及编程语言 › 快速傅里叶变换（转）

ljz 发表于 2012-10-14 00:56

快速傅里叶变换（转）

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#definePI3.14159265358979323846
structCOMPLEX
{
    floatre;
    floatim;
} cplx , * Hfield , * S , * R , * w;
intn , m;
intln , lm;
voidinitiate ();
voiddfft ();
voidrdfft ();
voidshowresult ();
voidfft (int l , int k);
intreverse (int t , int k);
voidW (int l);
intloop (int l);
voidconjugate ();
voidadd (structCOMPLEX* x , structCOMPLEX* y , structCOMPLEX* z);
voidsub (structCOMPLEX* x , structCOMPLEX* y , structCOMPLEX* z);
voidmul (structCOMPLEX* x , structCOMPLEX* y , structCOMPLEX* z);
structCOMPLEX* Hread(int i , int j);
voidHwrite (int i , int j , structCOMPLEX x);
voidmain ()
{
    initiate ();   
    printf("\n原始数据:\n");   
    showresult();   
    getchar ();   
    dfft ();   
    printf("\n快速复利叶变换后的结果:\n");   
    showresult ();   
    getchar ();   
    rdfft ();   
    printf("\n快速复利叶逆变换后的结果:\n");   
    showresult ();   
    getchar ();   
    free (Hfield);
}
voidinitiate ()
{//程序初始化操作，包括分配内存、读入要处理的数据、进行显示等   
    FILE* df;
    df = fopen ("data.txt" , "r");   
    fscanf (df , "%5d" , &n);   
    fscanf (df , "%5d" , &m);   
    if ((ln = loop (n)) == -1)   
    {   
      printf (" 列数不是2的整数次幂 ");      
      exit (1);
    }   
    if ((lm = loop (m)) == -1)   
    {   
      printf (" 行数不是2的整数次幂 ");      
      exit (1);
    }
    Hfield = (structCOMPLEX *) malloc (n * m * sizeof (cplx));
    if (fread (Hfield , sizeof (cplx) , m * n , df) != (unsigned) (m * n))
    {
      if (feof (df)) printf (" Premature end of file ");   
      else printf (" File read error ");
    }   
    fclose (df);
}
voiddfft ()
{//进行二维快速复利叶变换
    inti , j;   
    intl , k;
    l = n;
    k = ln;   
    w = (structCOMPLEX *) calloc (l , sizeof (cplx));
    R = (structCOMPLEX *) calloc (l , sizeof (cplx));
    S = (structCOMPLEX *) calloc (l , sizeof(cplx));
    W (l);
    for ( i = 0 ; i < m ; i++ )
    {//按行进行快速复利叶变换
      for (j = 0 ; j < n ; j++)
      {            
            S.re = Hread (i , j)->re;
            S.im = Hread (i , j)->im;
      }
      fft(l , k);
      for (j = 0 ; j < n ; j++)
            Hwrite (i , j , R);
    }
    free (R);
    free (S);
    free (w);   
    l = m;
    k = lm;
    w = (structCOMPLEX *) calloc (l , sizeof (cplx));
    R = (structCOMPLEX *) calloc (l , sizeof (cplx));
    S = (structCOMPLEX *) calloc (l , sizeof (cplx));
    W (l);
    for (i = 0 ; i < n ; i++)
    {//按列进行快速复利叶变换
      for(j = 0 ; j < m ; j++)
      {
            S.re = Hread(j , i)->re;   
            S.im = Hread(j , i)->im;
      }
      fft(l , k);
      for (j = 0 ; j < m ; j++)
            Hwrite (j , i , R);
    }
    free (R);
    free (S);
    free (w);
}
voidrdfft ()
{
    conjugate ();   
    dfft ();
    conjugate ();
}
voidshowresult ()
{
    inti , j;   
    for (i = 0 ; i < m ; i++)
    {
      printf ( " \n第%d行\n " , i);   
      for (j = 0 ; j < n ; j++)
      {
            if (j % 4 == 0) printf (" \n ");   
            printf(" (%5.2f,%5.2fi) " , Hread (i , j)->re , Hread (i , j)->im);
      }
    }
}
voidfft (int l , int k)
{
    inti , j , s , nv , t;   
    floatc;
    structCOMPLEXmp , r;
    nv = l;
    c = (float) l;
    c = pow (c , 0.5);
    for (i = 0 ; i < k ; i++)
    {
      for (t = 0 ; t < l ; t += nv)   
      {
            for (j = 0 ; j < nv / 2 ; j++)   
            {
                s = (t + j) >> (k - i -1);   
                s = reverse(s , k);
                r.re = S.re;
                r.im = S.im;
                mul (&w , &S , &mp);/////////讲解传递结构指针和结构本身的区别
                add (&r , &mp , &S);
                sub (&r , &mp , &S);               
            }
      }
      nv = nv >> 1;      
    }
    for (i = 0 ; i < l ; i++)
    {
      j = reverse(i , k);   
      R.re = S.re / c;
      R.im = S.im / c;
    }
}
intreverse (int t , int k)
{
    inti , x , y;   
    y = 0;
    for (i = 0 ; i < k ; i++)
    {
      x = t & 1;   
      t = t >> 1;
      y = (y << 1) + x;      
    }
    return y;
}
voidW (int l)
{
    int i;   
    float c , a;
    c = (float) l;
    c = 2 * PI / c;
    for (i = 0 ; i < l ; i++)
    {      
      a = (float) i;
      w.re = (float) cos(a * c);
      w.im = -(float) sin(a * c);
    }
}
intloop (int l)
{//检验输入数据是否为2的整数次幂，如果是返回用2进制表示时的位数
    inti , m;   
    if (l != 0)   
    {   
      for (i = 1 ; i < 32 ; i++)      
      {      
            m = l >> i;
            if (m == 0)            
                break;
      }      
      if (l == (1 << (i - 1)))      
            return i - 1;
    }   
    return -1;
}
voidconjugate ()
{//求复数矩阵的共轭矩阵
    inti , j;
    for (i = 0 ; i < m ; i++)   
    {
      for (j = 0 ; j < n ; j++)   
      {
            Hread (i , j)->im *= -1;
      }
    }
}
structCOMPLEX* Hread (int i , int j)
{//按读矩阵方式返回Hfield中指定位置的指针
    return (Hfield + i * n + j);
}
voidHwrite (int i , int j , structCOMPLEX x)
{//按写矩阵方式将复数结构x写到指定的Hfield位置上
    (Hfield + i * n + j)->re = x.re;   
    (Hfield + i * n + j)->im = x.im;
}
voidadd (structCOMPLEX* x , structCOMPLEX* y , structCOMPLEX* z)
{//定义复数加法
    z->re = x->re + y->re;
    z->im = x->im + y->im;   
}
voidsub (structCOMPLEX* x , structCOMPLEX* y , structCOMPLEX* z)
{//定义复数减法
    z->re = x->re - y->re;
    z->im = x->im - y->im;
}
voidmul (structCOMPLEX* x , structCOMPLEX* y , structCOMPLEX* z)
{//定义复数乘法
       z->re = (x->re) * (y->re) - (x->im) * (y->im);   
    z->im = (x->im) * (y->re) + (x->re) * (y->im);
}

Rainyboy 发表于 2012-10-18 16:15

用C/C++写出来挺不容易，要是使用PYTHON等脚本语言应该会更简单一些。

xiaojiong8 发表于 2013-7-20 14:09

完全看不懂。。。。

wanxiaoliu 发表于 2013-11-1 15:09

{:{26}:}膜拜一下

飞雪域 发表于 2013-11-4 15:20

感谢分享，好强悍

viomag 发表于 2013-11-5 20:29

感谢分享

江树空 发表于 2013-11-7 09:10

{:{39}:}这个好，终于在这里找到了

牛小贱 发表于 2013-12-15 21:08

感谢分享 ！但是看不大懂啊{:{26}:}

mayaview 发表于 2014-2-25 00:53

Rainyboy 发表于 2012-10-18 16:15
用C/C++写出来挺不容易，要是使用PYTHON等脚本语言应该会更简单一些。

不要尝试写FFT，实现一个经典算法非常容易，但是做一个真正的快速FT是很困难的，去看看fftpack和fftw之争就知道这个问题完全不是几行代码可以搞定的事情。这些成熟稳定的算法还是用现成库实在。说到数值分析，现在好些学校还在讲什么LU分解，最速下降法，哎，真的有用吗？

上完一个学期的数值分析，天天在讨论收敛性，收敛速率，SOR最佳因子。结果ANSYS用的是预处理共轭梯度法，感觉完全被骗了一样。

impulse 发表于 2014-2-25 08:49

mayaview 发表于 2014-2-25 00:53
不要尝试写FFT，实现一个经典算法非常容易，但是做一个真正的快速FT是很困难的，去看看fftpack和fftw之争 ...

你可以不去学习，但是作为一个国家，都不去学习，那么你永远只能跟在人家屁股后面，Ansys一年赚取国人多少银子，再说说仪器，每个领域所有高端仪器都是老外的，真准备子子孙孙也这样了？对了，中国发展北斗做什么，不是有GPS吗？

mayaview 发表于 2014-2-25 09:46

impulse 发表于 2014-2-25 08:49
你可以不去学习，但是作为一个国家，都不去学习，那么你永远只能跟在人家屁股后面，Ansys一年赚取国人多 ...

完全不是一个概念啊，你觉得单凭写个这样的C语言的FFT就能超过FFTW？像这样重复的造轮子，啥时候都赶不上现在的国外水平，还别说超过人家了。设备和北斗里或许多少都还有点前人没做过的东西吧，这种仅仅是实现一下经典算法的代码无非就是重复造轮子，而且还是质量很差的轮子。

impulse 发表于 2014-2-25 10:20

本帖最后由 impulse 于 2014-2-25 10:21 编辑

再好的工匠做的第一个轮子就能装在空客上？你可以指出这段代码哪些地方有不足，但是不能打击这种精神，创新也是要有基础的，有继承的。

mayaview 发表于 2014-2-25 12:07

impulse 发表于 2014-2-25 10:20
再好的工匠做的第一个轮子就能装在空客上？你可以指出这段代码哪些地方有不足，但是不能打击这种精神，创新 ...

这是转帖，谢谢！而且我已经提供了两个更好的库可供参考，那两个库都是开源的。做过点FFT的都应该知道（如果只是在Matlab里面用fft命令的就算了吧）FFT实现的经典算法基本是不可以拿来实际运用，了解了解基本算法还行，真正要用需要大量的优化和特殊情况的处理，不是这样的代码可以解决的。

查看完整版本: 快速傅里叶变换（转）