本帖最后由 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 Pi 3.1416 //定义圆周率Pi
#define Fn 10 //定义信号的有效频率
#define Fs 1000 //定义信号的采样频率
//-----定义复数结构体-----------------------//
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[Num],W,Temp;
for (k1=0;k1<=r2-1;k1++)
{
for(k0=0;k0<=r1-1;k0++)
{
ZeroCompx(Val,Num); //注意：每次计算FFT后的X(k)时，重新初始化Val[Num]
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[n1*r2+n0],W);
Val[n0]=add(Val[n0],Temp);
}
//计算X1'(k0,n0)的过程
p=n0*k0;
ps=2*Pi/N*p;
W.real=cos(ps);
W.imag=-sin(ps);
Val[n0]=EE(Val[n0],W);
//计算X2(k0,k0)的过程
p=r1*n0*k1;
ps=2*Pi/N*p;
W.real=cos(ps);
W.imag=-sin(ps);
Val[n0]=EE(Val[n0],W);
//将X2(k0,k0)通过X（k1*r1+k0)恢复为X(k)的过程
Xout[k1*r1+k0]=add(Xout[k1*r1+k0],Val[n0]);
} //end for(n0=0)
} //end for(k0=0)
} //end for(k1=0)
return ;
}
float ResultFFT[Num]; //定义FFT输出的幅值
struct Compx Source[Num]; //定义FFT的采样点存放数组
struct Compx Result[Num]; //定义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[icounter].real = sin(2*Pi*Fn*icounter/1000)+0.5*sin(2*Pi*(5*Fn)*icounter/1000); //采样信号的实部
Source[icounter].imag=0; //采样信号的虚部为0
printf("%.4f , ",Source[icounter].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[0] = ResultFFT[0]/Num; //换成实际幅值时直流量除以Num
printf("%.4f ",ResultFFT[0]);
for(i=1;i<Num;i++)
{
ResultFFT = ResultFFT/(Num/2); //换成实际幅值时交流量除以Num/2
printf("%.4f ",ResultFFT);
}
printf("\n运行结束n");
}
}
}
}
运行结果如下：
|
楼主
标题置顶
标题高亮
