您的位置 首页 动态

快速傅里叶变换FFT的C程序代码完成

  一、完全了解傅里叶变换   快速傅里叶变换(Fast Fourier Transform)是离散傅里叶变换的一种快速算法,简称FFT,经过FFT能够将一个信号从时域变换到频域。…

  一、完全了解傅里叶变换

  快速傅里叶变换(Fast Fourier Transform)是离散傅里叶变换的一种快速算法,简称FFT,经过FFT能够将一个信号从时域变换到频域。
  模拟信号经过A/D转化变为数字信号的进程称为采样。为确保采样后信号的频谱形状不失真,采样频率有必要大于信号中最高频率成分的2倍,这称之为采样定理。
  假定采样频率为fs,采样点数为N,那么FFT成果便是一个N点的复数,每一个点就对应着一个频率点,某一点n(n从1开端)表明的频率为:fn=(n-1)*fs/N。
  举例阐明:用1kHz的采样频率采样128点,则FFT成果的128个数据即对应的频率点分别是0,1k/128,2k/128,3k/128,…,127k/128 Hz。
  这个频率点的幅值为:该点复数的模值除以N/2(n=1时是直流重量,其幅值是该点的模值除以N)。

  二、傅里叶变换的C言语编程

  1、关于快速傅里叶变换FFT,第一个要处理的问题便是码位倒序。

  假定一个N点的输入序列,那么它的序号二进制数位数便是t=log2N.
  码位倒序要处理两个问题:①将t位二进制数倒序;②将倒序后的两个存储单元进行交流。
  假如输入序列的天然顺序号i用二进制数表明,例如若最大序号为15,即用4位就可表明n3n2n1n0,则其倒序后j对应的二进制数便是n0n1n2n3,那么怎样才能完结倒序呢?使用C言语的移位功用!
  程序如下,我不多说,看不懂者智商一定在180以下!

  复数类型界说及其运算
  #define N 64 //64点
  #define log2N 6 //log2N=6
  /*复数类型*/
  typedef struct
  {
  float real;
  float img;
  }complex;
  complex xdata x[N]; //输入序列
  /*复数加法*/
  complex add(complex a,complex b)
  {
  complex c;
  c.real=a.real+b.real;
  c.img=a.img+b.img;
  return c;
  }
  /*复数减法*/
  complex sub(complex a,complex b)
  {
  complex c;
  c.real=a.real-b.real;
  c.img=a.img-b.img;
  return c;
  }
  /*复数乘法*/
  complex mul(complex a,complex b)
  {
  complex c;
  c.real=a.real*b.real – a.img*b.img;
  c.img=a.real*b.img + a.img*b.real;
  return c;
  }
  /***码位倒序函数***/
  void Reverse(void)
  {
  unsigned int i,j,k;
  unsigned int t;
  complex temp;//暂时交流变量
  for(i=0;iN;i++)//从第0个序号到第N-1个序号
  {
  k=i;//当时第i个序号
  j=0;//存储倒序后的序号,先初始化为0
  for(t=0;tlog2N;t++)//共移位t次,其间log2N是事前宏界说算好的
  {
  j=1;
  j|=(k1);//j左移一位然后加上k的最低位
  k>>=1;//k右移一位,次低位变为最低位
  }
  if(j>i)//假如倒序后大于原序数,就将两个存储单元进行交流(判别j>i是为了避免重复交流)
  {
  temp=x;
  x=x[j];
  x[j]=temp;
  }
  }
  }

  2、第二个要处理的问题便是蝶形运算

  ①第1级(第1列)每个蝶形的两节点“间隔”为1,第2级每个蝶形的两节点“间隔”为2,第3级每个蝶形的两节点“间隔”为4,第4级每个蝶形的两节点“间隔”为8。由此推得,
  第m级蝶形运算,每个蝶形的两节点“间隔”L=2m-1。
  ②关于16点的FFT,第1级有16组蝶形,每组有1个蝶形;第2级有4组蝶形,每组有2个蝶形;第3级有2组蝶形,每组有4个蝶形;第4级有1组蝶形,每组有8个蝶形。由此可推出,
  关于N点的FFT,第m级有N/2L组蝶形,每组有L=2m-1个蝶形。
  ③旋转因子的确认
  以16点FFT为例,第m级第k个旋转因子为,其间k=0~2m-1-1,即第m级共有2m-1个旋转因子,依据旋转因子的可约性,,所以第m级第k个旋转因子为,其间k=0~2m-1-1

  为进步FFT的运算速度,咱们能够事前树立一个旋转因子数组,然后经过查表法来完结。

  complex code WN[N]=//旋转因子数组
  { //为节约CPU核算时刻,旋转因子选用查表处理
  //★依据实践FFT的点数N,该表数据需自行修正
  //以下成果经过Excel主动生成
  // WN[k].real=cos(2*PI/N*k);
  // WN[k].img=-sin(2*PI/N*k);
  {1.00000,0.00000},{0.99518,-0.09802},{0.98079,-0.19509},{0.95694,-0.29028},
  {0.92388,-0.38268},{0.88192,-0.47140},{0.83147,-0.55557},{0.77301,-0.63439},
  {0.70711,-0.70711},{0.63439,-0.77301},{0.55557,-0.83147},{0.47140,-0.88192},
  {0.38268,-0.92388},{0.29028,-0.95694},{0.19509,-0.98079},{0.09802,-0.99518},
  {0.00000,-1.00000},{-0.09802,-0.99518},{-0.19509,-0.98079},{-0.29028,-0.95694},
  {-0.38268,-0.92388},{-0.47140,-0.88192},{-0.55557,-0.83147},{-0.63439,-0.77301},
  {-0.70711,-0.70711},{-0.77301,-0.63439},{-0.83147,-0.55557},{-0.88192,-0.47140},
  {-0.92388,-0.38268},{-0.95694,-0.29028},{-0.98079,-0.19509},{-0.99518,-0.09802},
  {-1.00000,0.00000},{-0.99518,0.09802},{-0.98079,0.19509},{-0.95694,0.29028},
  {-0.92388,0.38268},{-0.88192,0.47140},{-0.83147,0.55557},{-0.77301,0.63439},
  {-0.70711,0.70711},{-0.63439,0.77301},{-0.55557,0.83147},{-0.47140,0.88192},
  {-0.38268,0.92388},{-0.29028,0.95694},{-0.19509,0.98079},{-0.09802,0.99518},
  {0.00000,1.00000},{0.09802,0.99518},{0.19509,0.98079},{0.29028,0.95694},
  {0.38268,0.92388},{0.47140,0.88192},{0.55557,0.83147},{0.63439,0.77301},
  {0.70711,0.70711},{0.77301,0.63439},{0.83147,0.55557},{0.88192,0.47140},
  {0.92388,0.38268},{0.95694,0.29028},{0.98079,0.19509},{0.99518,0.09802}
  };

  3、算法完结

  咱们现已知道,N点FFT从左到右共有log2N级蝶形,每级有N/2L组,每组有L个。所以FFT的C言语编程只需用3层循环即可完结:最外层循环完结每一级的蝶形运算(整个FFT共log2N级),中间层循环完结每一组的蝶形运算(每一级有N/2L组),最内层循环完结独自1个蝶形运算(每一组有L个)。
  /***【快速傅里叶变换】***/
  void FFT(void)
  {
  unsigned int i,j,k,l;
  complex top,bottom,xW;
  Reverse(); //码位倒序
  for(i=0;ilog2N;i++) /*共log2N级*/
  { //一级蝶形运算
  l=1i;//l等于2的i次方
  for(j=0;jN;j+=2*l) /*每L个蝶形是一组,每级有N/2L组*/
  { //一组蝶形运算
  for(k=0;kl;k++) /*每组有L个*/
  { //一个蝶形运算
  xW=mul(x[j+k+l],WN[N/(2*l)*k]); //碟距离为l
  top=add(x[j+k],xW); //每组的第k个蝶形
  bottom=sub(x[j+k],xW);
  x[j+k]=top;
  x[j+k+l]=bottom;
  }
  }
  }
  }

  三、FFT核算成果验证

  随意输入一个64点序列,例如
  x[N]={{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0}};
  在Keil中Debug检查数组变量x的FFT核算成果并与MATLAB核算成果进行比对,能够发现十分精确,阐明程序编写正确!

声明:本文内容来自网络转载或用户投稿,文章版权归原作者和原出处所有。文中观点,不代表本站立场。若有侵权请联系本站删除(kf@86ic.com)https://www.86ic.net/news/dongtai/166554.html

为您推荐

联系我们

联系我们

在线咨询: QQ交谈

邮箱: kf@86ic.com

关注微信
微信扫一扫关注我们

微信扫一扫关注我们

返回顶部