#include"stdio.h"
#include "math.h"
typedef unsigned char u8;
typedef unsigned int u16;
typedef unsigned long u32;
#define PI 3.141592653589793238462643 //定义圆周率
#define FFT_N 64 //采样点数
typedef struct Compex //复数结构体
{
double real;
double image;
}COMPLEX;
COMPLEX FFT_result[FFT_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}}; //输入输出数据存储数组
void Init_forward(void); //倒位序，采用所谓的雷德算法
COMPLEX MUL_EE(COMPLEX X1,COMPLEX X2); //复数乘法公式
COMPLEX ADD_EE(COMPLEX X1,COMPLEX X2); //复数加法公式
COMPLEX SUB_EE(COMPLEX X1,COMPLEX X2); //复数减法公式
/**************************雷德算法思想*************************
* 自然顺序 倒位序
* 000 000
* 001 100
* 010 010
* 011 110
* 100 001
* 101 101
* 110 011
* 111 000
由此可见到位序实际上就是镜像运算，然而我们没有采用镜像算法，（据说可以用汇编来实现比较容易）
我们所要做的工作是：
1.如果我们已知自然顺序的一个数想要知道下一个数只需要将当前数加1即可
2.再观察倒位序之后数据的规律，，，
3.如果我们知道倒位序后的其中一个数，要想求出下一个数，同样也可以采用加法，
然而，此处的加法跟我们从小学习的不同，我们要实现的是向低位进位的加法（这里看不懂脑袋砸两下）
话不多说，此函数就是实现这个功能仔细分析，如果能看懂说明逻辑思维不错
*/
/**********************************************************
*函数名称：COMPLEX MUL_EE(COMPLEX X1,COMPLEX X2)
*函数功能：实现复数乘法
*输入参数：COMPLEX X1,COMPLEX X2
*返回值：COMPLEX c
***********************************************************/
COMPLEX MUL_EE(COMPLEX X1,COMPLEX X2)
{
COMPLEX c;
c.real = X1.real*X2.real - X1.image*X2.image;
c.image = X1.real*X2.image + X1.image*X2.real;
return c;
}
/**********************************************************
*函数名称：COMPLEX ADD_EE(COMPLEX X1,COMPLEX X2)
*函数功能：进行复数加法
*输入参数：COMPLEX X1,COMPLEX X2
*返回值：
***********************************************************/
COMPLEX ADD_EE(COMPLEX X1,COMPLEX X2)
{
COMPLEX c;
c.real = X1.real + X2.real;
c.image = X1.image + X2.image;
return c;
}
/**********************************************************
*函数名称：COMPLEX Dcc_EE(COMPLEX X1,COMPLEX X2)
*函数功能：进行复数减法
*输入参数：COMPLEX X1,COMPLEX X2
*返回值：COMPLEX c
***********************************************************/
COMPLEX SUB_EE(COMPLEX X1,COMPLEX X2)
{
COMPLEX c;
c.real = X1.real - X2.real;
c.image = X1.image - X2.image;
return c;
}
/**********************************************************
*函数名称：void Init_forward(void)
*函数功能：倒位序
*输入函数：void
*返回值：void
***********************************************************/
void Init_forward(void)
{
u8 I,J,LH,N1,K ;
COMPLEX T; //替换结构体
LH = FFT_N/2; //N/2
J = LH;
N1 = FFT_N - 2;
for(I=1;I<N1;I++) //从1到N-2开始倒位序
{
if(I<J) //此处的意思是当I不等于J时交换位置
{ //然而I>J时不交换因为之前已经交换了
T = FFT_result[I];
FFT_result[I] = FFT_result[J];
FFT_result[J] = T;
}
K = LH; //将给K赋值N/2
while(J>=K) //循环，，判断所需判断的位是否为1
{
J = J-K;
K = K/2;
}
J = J+K;
}
}
/**************************************************************
*函数名称：void FFT_Run(void)
*函数功能：进行快速傅里叶运算
*输入参数：void
*返回值：void
***************************************************************/
void FFT_Run(void)
{
u8 B,P,K;
u8 L = 0; //蝶形变换级数
u8 M = 0; //N = 2^M
u8 J;
u8 FFT_N1 = FFT_N;
COMPLEX Result_Wn,Result_MUL,Result_ADD,Result_SUB;
Init_forward(); //进行倒位序运算
for(M=1; (FFT_N1=FFT_N1/2)!=1;M++); //计算蝶形级数
for(L=1;L<=M;L++)
{
B = pow(2,L-1); // 旋转因子个数
for(J=0;J<=B-1;J++)
{
P = pow(2,M-L)*J; //旋转因子系数
for(K=J;K<FFT_N;K=K+pow(2,L))
{
Result_Wn.real = cos((2*PI/FFT_N)*P);
Result_Wn.image = -sin((2*PI/FFT_N)*P);
Result_MUL = MUL_EE(Result_Wn,FFT_result[K+B]); //复数乘法
Result_ADD = ADD_EE(FFT_result[K],Result_MUL); //复数加法
Result_SUB = SUB_EE(FFT_result[K],Result_MUL); //复数减法
FFT_result[K] = Result_ADD; //把加法后的结果放到 FFT_result[K]
FFT_result[K+B] = Result_SUB; //把减法之后的结果放到FFT_result[K+B]
}
}
}
}
void main(void)
{
u8 a;
u8 M;
u8 FFT_N1 = FFT_N;
FFT_Run();
for(a=0;a<FFT_N;a++)
{
printf("%f",FFT_result[a].real/100);
printf(" ");
printf("%f",FFT_result[a].image/100);
printf("\n");
}
a= pow(2,3);
printf("%d",a);
printf("\n");
}
经过一星期已搞定，学弟学妹可以看看。。
|
QQ好友和群
QQ空间
腾讯微博
腾讯朋友
收藏
淘帖
顶
踩
