#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 //采樣點數(shù)
typedef struct Compex //復數(shù)結(jié)構(gòu)體
{
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}}; //輸入輸出數(shù)據(jù)存儲數(shù)組
void Init_forward(void); //倒位序，采用所謂的雷德算法
COMPLEX MUL_EE(COMPLEX X1,COMPLEX X2); //復數(shù)乘法公式
COMPLEX ADD_EE(COMPLEX X1,COMPLEX X2); //復數(shù)加法公式
COMPLEX SUB_EE(COMPLEX X1,COMPLEX X2); //復數(shù)減法公式
/**************************雷德算法思想*************************
* 自然順序 倒位序
* 000 000
* 001 100
* 010 010
* 011 110
* 100 001
* 101 101
* 110 011
* 111 000
由此可見到位序?qū)嶋H上就是鏡像運算，然而我們沒有采用鏡像算法，（據(jù)說可以用匯編來實現(xiàn)比較容易）
我們所要做的工作是：
1.如果我們已知自然順序的一個數(shù)想要知道下一個數(shù)只需要將當前數(shù)加1即可
2.再觀察倒位序之后數(shù)據(jù)的規(guī)律，，，
3.如果我們知道倒位序后的其中一個數(shù)，要想求出下一個數(shù)，同樣也可以采用加法，
然而，此處的加法跟我們從小學習的不同，我們要實現(xiàn)的是向低位進位的加法（這里看不懂腦袋砸兩下）
話不多說，此函數(shù)就是實現(xiàn)這個功能仔細分析，如果能看懂說明邏輯思維不錯
*/
/**********************************************************
*函數(shù)名稱：COMPLEX MUL_EE(COMPLEX X1,COMPLEX X2)
*函數(shù)功能：實現(xiàn)復數(shù)乘法
*輸入?yún)?shù)：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;
}
/**********************************************************
*函數(shù)名稱：COMPLEX ADD_EE(COMPLEX X1,COMPLEX X2)
*函數(shù)功能：進行復數(shù)加法
*輸入?yún)?shù)：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;
}
/**********************************************************
*函數(shù)名稱：COMPLEX Dcc_EE(COMPLEX X1,COMPLEX X2)
*函數(shù)功能：進行復數(shù)減法
*輸入?yún)?shù)：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;
}
/**********************************************************
*函數(shù)名稱：void Init_forward(void)
*函數(shù)功能：倒位序
*輸入函數(shù)：void
*返回值：void
***********************************************************/
void Init_forward(void)
{
u8 I,J,LH,N1,K ;
COMPLEX T; //替換結(jié)構(gòu)體
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時不交換因為之前已經(jīng)交換了
T = FFT_result[I];
FFT_result[I] = FFT_result[J];
FFT_result[J] = T;
}
K = LH; //將給K賦值N/2
while(J>=K) //循環(huán)，，判斷所需判斷的位是否為1
{
J = J-K;
K = K/2;
}
J = J+K;
}
}
/**************************************************************
*函數(shù)名稱：void FFT_Run(void)
*函數(shù)功能：進行快速傅里葉運算
*輸入?yún)?shù)：void
*返回值：void
***************************************************************/
void FFT_Run(void)
{
u8 B,P,K;
u8 L = 0; //蝶形變換級數(shù)
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++); //計算蝶形級數(shù)
for(L=1;L<=M;L++)
{
B = pow(2,L-1); // 旋轉(zhuǎn)因子個數(shù)
for(J=0;J<=B-1;J++)
{
P = pow(2,M-L)*J; //旋轉(zhuǎn)因子系數(shù)
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]); //復數(shù)乘法
Result_ADD = ADD_EE(FFT_result[K],Result_MUL); //復數(shù)加法
Result_SUB = SUB_EE(FFT_result[K],Result_MUL); //復數(shù)減法
FFT_result[K] = Result_ADD; //把加法后的結(jié)果放到 FFT_result[K]
FFT_result[K+B] = Result_SUB; //把減法之后的結(jié)果放到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");
}
經(jīng)過一星期已搞定，學弟學妹可以看看。。
|
