论坛首页 › DSP 技术 › 任意数为基数的FFT 算法

icekoor 发表于 2014-8-7 11:18

任意数为基数的FFT 算法

本帖最后由 icekoor 于 2014-8-11 08:44 编辑

简单描述如下：
如 N 为复合数，可分解为两个整数P与Q的乘积，像前面以2为基数时一样，FFT的基本思想是将DFT的点数尽量分小，因此，在N=PQ情况下，也希望将N点的DFT分解为P个Q点DFT或Q个P点DFT，以减少计算量。



补零的方法不采用，因为影响频谱。

FFTW库有些复杂，实现起来比较消耗时间。
（如果有FFTW与DSPC6000配置过的程序，推荐给我哈，非常感谢）
//------------------------------------------------------------------------------------
参考网上找到的程序（吕娴娜的混合基FFt算法）完成了任意数为基数的FFT 算法。
吕娴娜的混合基FFt算法，有很多问题，运行结果不对，经过修改后，可以直接在VS C++环境下运行；
/*------------------------------------
作者：icekoor
时间：2014.08.08
功能：任意数为基数的FFT 算法
说明：直接添加到VS C++中就可以运行；
      FFT1.cpp : 定义控制台应用程序的入口点
------------------------------------*/
#include "stdafx.h"
#include "math.h"
#include "stdio.h"
#include "time.h"
#include "Windows.h"
#define Num 100       //定义采样点数量
#define Pi3.1416   //定义圆周率Pi
#define Fn10       //定义信号的有效频率
#define Fs1000   //定义信号的采样频率
//-----定义复数结构体-----------------------//
struct Compx
{
float real;
   float imag;
} Compx ;
//-----复数结构体清零函数-------------------//
void ZeroCompx(struct Compx *X,unsigned int N)
{
unsigned int i;
for(i=0;i<N;i++)
{
X.imag = 0;
X.real = 0;
}
return;
}
//-----复数乘法运算函数---------------------//
struct Compx EE(struct Compx b1,struct Compx b2)
{
struct Compx b3;
b3.real=b1.real*b2.real-b1.imag*b2.imag;
b3.imag=b1.real*b2.imag+b1.imag*b2.real;
return(b3);
}
//-----复数加法运算函数---------------------//
struct Compx add(struct Compx a,struct Compx b)
{
struct Compx c;
c.real=a.real+b.real;
c.imag=a.imag+b.imag;
return(c);
}
//-----采样点个数为N = r1*r2的FFT函数------//
/*------------------------------------------
n = n1*r2 + n0; n1,k0的范围0,1,2......r1-1;
k = k1*r1 + k0; k1,n0的范围0,1,2......r2-1;
-------------------------------------------*/
void FFT(struct Compx *Xin, struct Compx *Xout, int N, int r1, int r2)
{
   int k0,k1,n0,n1;            
   float p ,ps ;
   struct Compx Val,W,Temp;
   for (k1=0;k1<=r2-1;k1++)
   {
for(k0=0;k0<=r1-1;k0++)
      {   
   ZeroCompx(Val,Num);      //注意：每次计算FFT后的X(k)时，重新初始化Val
   for(n0=0;n0<=r2-1;n0++)
   {
    //计算X1(k0,n0)的过程
    for(n1=0;n1<=r1-1;n1++)
    {
                p=r2*n1*k0;   
                ps=2*Pi/N*p;
                W.real=cos(ps);
               W.imag=-sin(ps);
               Temp=EE(Xin,W);
               Val=add(Val,Temp);
    }
    //计算X1'(k0,n0)的过程
    p=n0*k0;
    ps=2*Pi/N*p;
    W.real=cos(ps);
    W.imag=-sin(ps);
    Val=EE(Val,W);
    //计算X2(k0,k0)的过程
    p=r1*n0*k1;
    ps=2*Pi/N*p;
    W.real=cos(ps);
    W.imag=-sin(ps);
    Val=EE(Val,W);
    //将X2(k0,k0)通过X（k1*r1+k0)恢复为X(k)的过程
    Xout=add(Xout,Val);
             } //end for(n0=0)
         } //end for(k0=0)
      } //end for(k1=0)
   return ;
}
float    ResultFFT;   //定义FFT输出的幅值
struct   Compx Source;//定义FFT的采样点存放数组
struct   Compx Result;//定义FFT的运算结果存放数组
int main()
{
unsigned int r1=10,r2=10;
//------------------
int i;
int milseconds1 =0;
int milseconds0 =0;
int icounter = 0;
//-------Window系统的时间变量-----------//
SYSTEMTIME sys;
while(1)
{
   GetLocalTime( &sys );
   milseconds1=sys.wMilliseconds;
      if(milseconds1 != milseconds0)    //以此为1kHz的采样频率；
   {
   Source.real= sin(2*Pi*Fn*icounter/1000)+0.5*sin(2*Pi*(5*Fn)*icounter/1000);   //采样信号的实部
      Source.imag=0;                               //采样信号的虚部为0
   printf("%.4f , ",Source.real);
   printf("%d\n",icounter);
   milseconds0 = milseconds1;
   icounter++;
   if(icounter >= Num)                     //每次采集到Num个点，进行一次FFT变换
   {
   FFT(Source,Result,Num,r1,r2);   //采用任意基数的FFT变换
   printf("\nFFT分析结果:\n");
   for(i=0;i<Num;i++)
   {
      printf("%.4f",      Result.real);
      printf("+%.4fj    ",Result.imag);
   }
   printf("\nFFT的原始幅值:\n");
   for(i=0;i<Num;i++)
   {
      ResultFFT=sqrt(pow(Result.real,2)+pow(Result.imag,2)); //计算FFT变换后的幅值，为原始的FFT模值
      printf("%.4f   ",ResultFFT);
   }
   printf("\nFFT的实际幅值:\n");
   ResultFFT = ResultFFT/Num;         //换成实际幅值时直流量除以Num
   printf("%.4f   ",ResultFFT);
   for(i=1;i<Num;i++)
   {
      ResultFFT = ResultFFT/(Num/2); //换成实际幅值时交流量除以Num/2
      printf("%.4f   ",ResultFFT);
   }
   printf("\n运行结束n");
   }
   }
}
}
运行结果如下：

icekoor 发表于 2014-8-7 11:22

自己顶一个

icekoor 发表于 2014-8-7 13:39

找到了一段靠谱的程序，推荐给大家：
/*****************fft programe*********************/
/*混合基FFT算法*/
//by SC07023029 吕娴娜
#include "math.h"
#include "stdio.h"

struct compx
{
        double real;
        double imag;
} compx ;


struct compx EE(struct compx b1,struct compx b2)
{
        struct compx b3;
        b3.real=b1.real*b2.real-b1.imag*b2.imag;
        b3.imag=b1.real*b2.imag+b1.imag*b2.real;
        return(b3);
}

struct compx add(struct compx a,struct compx b)
{
        struct compx c;
        c.real=a.real+b.real;
        c.imag=a.imag+b.imag;
        return(c);
}

void FFT(struct compx *xin,int N,int r1,int r2 )

{
        int k0,k1,n0,n1;
        double p ,ps ;
        float pi;
        struct compx v,w,t;
       {       
   for (k1=0;k1<=r1-1;k1++)
    {
        for(k0=0;k0<=r2-1;k0++)
      {
         pi=3.14159;
         for(n0=0;n0<=r2-1;n0++)
         {
             for(n1=0;n1<=r1-1;n1++)
             {
            p=r2*n1*k0;
            ps=2*pi/N*p;
            w.real=cos(ps);
            w.imag=-sin(ps);
            t=EE(xin,w);
            v=add(v,t);
            }
            p=n0*k0;
            ps=2*pi/N*p;
            w.real=cos(ps);
            w.imag=-sin(ps);
            v=EE(v,w);
            p=r1*n0*k1;
            ps=2*pi/N*p;
            w.real=cos(ps);
            w.imag=-sin(ps);
         v=EE(v,w);   
         xin=add(xin,v) ;
         }
       }
    }
}
return ;
}

icekoor 发表于 2014-8-7 13:40

以下是主程序：
/*****************main programe********************/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

floatresult;
structcompx s;
int   Num=16;
int r1=4,r2=4;
const float pp=3.14159;

main()
{

int i;
printf("原序列：\n");
for(i=0;i<16;i++)
{
s.real=sin(pp*i/32);
s.imag=0;
printf("%.4f    ",s.real);
}

FFT(s,Num,r1,r2);

for(i=0;i<16;i++)
{
printf("%.4f",s.real);
printf("+%.4fj    ",s.imag);
}
printf("\n功率密度:\n");
for(i=0;i<16;i++)
{
result=sqrt(pow(s.real,2)+pow(s.imag,2));
printf("+%.4f    ",result);
}
scanf("2");

}

zhangmangui 发表于 2014-8-7 23:15

学习啦

huangzj121 发表于 2014-8-8 08:06

基4 基16 基64 的都有
看复数乘法计算量 并不是越大越快

仔细看看这句
for(i=0;i<16;i++)
每次 都做16点 基二方法的快速DFT
这样 也行？

真扯淡

icekoor 发表于 2014-8-8 08:51

huangzj121 发表于 2014-8-8 08:06 static/image/common/back.gif
基4 基16 基64 的都有
看复数乘法计算量 并不是越大越快



没明白您的意思。虽说该程序是N=16，分为4x4的方式，但是也可以是N=15，分为3x5的方式，不是你说的基2基4的FFT。

icekoor 发表于 2014-8-8 13:44

程序已经调试完成，已结贴，谢谢！

huangzj121 发表于 2014-8-11 10:04

当别人 给你指出来的时候
不能够仔细 查找问题 自得洋洋
:shutup:

icekoor 发表于 2014-8-12 08:37

huangzj121 发表于 2014-8-11 10:04 static/image/common/back.gif
当别人 给你指出来的时候
不能够仔细 查找问题 自得洋洋

for(i=0;i<16;i++)
意思是就有16个点，分成4*4，计算一次就结束，不是每次都做16点的FFT。此程序的FFT是一种N=P*Q的方式，16个点也可以分成2*8，或者1*16，跟基2没什么关系，我讨论的是任意数为基的FFT分解方法，基2基4基16等等都仅仅是个个例，希望你理解。

ysszhk 发表于 2016-1-20 11:52

- -

小木欧尼 发表于 2016-1-24 17:10

好资料帮顶

edishen 发表于 2016-1-24 17:25

学习啦

查看完整版本: 任意数为基数的FFT 算法