📄 fft.c
字号:
#include<c8051f120.h>
#include<fft.h>
#include"math.h"
#include"math.h"
/*完成日期:5.14
/*修改日期:
/*入口参数:
X[]:需要变换的数据.变换后放回源数据区域,必须为符号
FFT_N:FFT点数.
RealFlag:实数变换的标志.如果为实数,由于处理安复数格式处理.需要变换,此标志控制变换与否.
//Int16Flag:输入数据为16位满刻度的标志.如果是满刻度的十六位数据.需要先右移一位,否则将溢出.
SignFlag:输入数据为有符号数据的标志,如果不是有符号数据,需要转换为有符号数据.
*/
/*一定要注意输出数据为有符号数据,所以如果显示处理后得数据要小心,
*/
void Int_FFT(int X[],unsigned int FFT_N,char windowtype,char RealFlag,char SignFlag)
{idata int DotDiv=2048/FFT_N; //由于查的表是的角度步进为2*pi/2048;故如果FFT点数少于2048时应查的表的步进要加倍.
idata int FFT_HFNum; //FFT点数
idata int FFT_Level=0; //FFT的级数
idata int LCnt; //级数的计数器
idata int FFT_Seg; //当前级FFT中的段数
idata int SCnt; //段数计数器
idata int SeGBTF; //每段包含的碟形算子的个数
idata int ButterFlySCnt; //每段中碟形算子序号
idata int ButterFlyCnt; //碟形算子总的序号地址
idata int ButterFlyCpl; //碟形算子对的另外一个数据的地址
idata int Ele_Inteval; //每个碟形算子涉及的两个元素的间隔
// idata int RealDirCrt,ImagDirCrt;
#if MAC120
idata int Rx1buff,Rx2buff,Ix1buff,Ix2buff;
idata int RSxbuff,ISxbuff;
idata char SFRPAGBUFF;
#else
idata long Rx1buff,Rx2buff,Ix1buff,Ix2buff;
idata long RSxbuff,ISxbuff;
idata long MACBUF,MACBUF1,MACBUF2,MACBUF3;
#endif
idata int Wnk;
idata int WnkInterval;//
idata int EXP2[16]={0x0001,0x0002,0x0004,0x0008,0x0010,0x0020,0x0040,0x0080,
0x0100,0x0200,0x0400,0x0800,0x1000,0x2000,0x4000,0x8000,
};//2的N次幂表.
/*120MCU时的页寄存器处理:
*/
#if MAC120
SFRPAGBUFF=SFRPAGE;
SFRPAGE=MAC0_PAGE;
#endif
/*FFT前期数据准备工作.
*/
if(SignFlag!=1)//如果为不是有符合数据,需要转化为有符合数据.
for(LCnt=0;LCnt<FFT_N;LCnt++)
{FFT_Seg=X[LCnt];
SCnt=FFT_Seg-0x8000;
X[LCnt]=SCnt;
}
FFT_N=FFT_N/2; //用N/2点求,是标准FFT数据量的1/2
FFT_HFNum=FFT_N;
FFT_Level=0;
for(;FFT_N!=0x0001;FFT_N=FFT_N>>1)FFT_Level++; //得到FFT级数
/*FFT加窗;
*/
#if MAC120
MAC0CF=0x01; //工作为乘法方式;
#endif
switch(windowtype)
{
case 'h':if(RealFlag)
{
for(LCnt=0,SCnt=0;SCnt<FFT_HFNum*2;LCnt=LCnt+DotDiv)
{RSxbuff=X[SCnt];
#if MAC120
MAC0A=hammingwindow[LCnt];
MAC0B=RSxbuff;
X[SCnt]=MACW1;
#else
// MACBUF1=hammingwindow[LCnt];
MACBUF=MACBUF1*RSxbuff;
MACBUF=MACBUF>>16;
X[SCnt]=MACBUF;
#endif
SCnt++;
}
}
break;
case 't':break;
}
//if(Int16FullFlag)for(LCnt=0;LCnt<FFT_N;LCnt++)X[LCnt]=X[LCnt]>>1;//如果输入为16位满刻度数据,需要转换为15位以防溢出.
/*FFT碟形运算:
*/ /*
#if MAC120
#else
Rx1buff=0;
Rx2buff=0;
RSxbuff=0;
for(LCnt=0;LCnt<FFT_HFNum*2;LCnt++)
{Rx1buff=X[LCnt];
RSxbuff= RSxbuff+Rx1buff;
}
#endif*/
/*****遍历所有级***********/
//DIF算法
for(LCnt=0;LCnt<FFT_Level;LCnt++) //一共要计算的级数
{FFT_Seg=EXP2[LCnt]; //得到本级FFT需要计算的段数
SeGBTF=EXP2[FFT_Level-LCnt-1];//得到本级FFT中每段的算子个数
Ele_Inteval=EXP2[FFT_Level-LCnt];//得到每个碟形算子涉及的两个元素的距离*2
ButterFlyCnt=0; //碟形算子序号复位
#if MAC120
MAC0CF=0x01; //工作于乘法方式,不是乘加方式
MAC0A=FFT_Seg;
MAC0B=2;
MAC0A=DotDiv;
MAC0B=MACW0;
WnkInterval=MACW0;//得到旋转因子的旋转间隔.WnkInterval=2*FFT_Seg*DotDiv;2是由于实数FFT时将N个点的FFT减少为N/2个点的FFT
//当FFT点数少于2048时步进为2048时的DotDiv(2048/FFT_N)倍.
MAC0CF=0x00; //工作于乘加方式,
#else
WnkInterval=FFT_Seg*DotDiv;
WnkInterval=WnkInterval<<1;
#endif
/*****遍历本级的所有段***********/
for(SCnt=0;SCnt<FFT_Seg;SCnt++)//本级FFT一共要计算的段数
{
Wnk=0;//复位旋转因子;
/*****遍历本段的所有碟形算子***********/
for(ButterFlySCnt=0;ButterFlySCnt<SeGBTF;ButterFlyCnt=ButterFlyCnt+2)//ButterFlySCnt++在循环体中.
{ButterFlyCpl=ButterFlyCnt+Ele_Inteval;//
Rx1buff=X[ButterFlyCnt];
Rx1buff=Rx1buff>>1;
Ix1buff=X[ButterFlyCnt+1];
Ix1buff=Ix1buff>>1;
Rx2buff=X[ButterFlyCpl];
Rx2buff=Rx2buff>>1;
Ix2buff=X[ButterFlyCpl+1];
Ix2buff=Ix2buff>>1; //以上为缓存源数据.每一级除以2防止溢出;
X[ButterFlyCnt]=Rx1buff+Rx2buff; //得到实部,这句可以一方面延时以提供MAC需要的时钟周期;
X[ButterFlyCnt+1]=Ix1buff+Ix2buff;//得到虚部,碟形算子的 A+B
RSxbuff=Rx1buff-Rx2buff; //碟形算子的(A-B)WnK;
RSxbuff=RSxbuff<<1;
ISxbuff=Ix1buff-Ix2buff;
ISxbuff=ISxbuff<<1;
#if MAC120 //如果MCU有硬件MAC则相应用MAC处理
MACW1=0;
MACW0=0;
MAC0A=RSxbuff;
MAC0B=WnC[Wnk]; //Re(A-B)*Cos(a);
MAC0A=ISxbuff; //Im(A-B)*Sin(a);
MAC0B=WnS[Wnk];
X[ButterFlyCpl]=MACW1; //由于正弦余弦表扩大了2^15倍.直接读入高16为.
MACW1=0; //得到Im(A-B)*Cos(a)-Re(A-B)*Sin(a)
MACW0=0;
MAC0A=~RSxbuff+1;
MAC0B=WnS[Wnk];
MAC0A=ISxbuff;
MAC0B=WnC[Wnk];
X[ButterFlyCpl+1]=MACW1;
#else
MACBUF=WnC[Wnk]; //得到Re(A-B)*Cos(a)+Im(A-B)*Sin(a)
MACBUF1=MACBUF*RSxbuff;
MACBUF=WnS[Wnk];
MACBUF2=MACBUF*ISxbuff;
MACBUF=MACBUF1+MACBUF2;
MACBUF=MACBUF>>16;
X[ButterFlyCpl]=MACBUF;
MACBUF2=~RSxbuff+1; //得到Im(A-B)*Cos(a)-Re(A-B)*Sin(a)
MACBUF=WnS[Wnk];
MACBUF1=MACBUF*MACBUF2;
MACBUF=WnC[Wnk];
MACBUF2=ISxbuff*MACBUF;
MACBUF=MACBUF1+MACBUF2;
MACBUF=MACBUF>>16;
X[ButterFlyCpl+1]=MACBUF;
#endif
ButterFlySCnt++;
Wnk=Wnk+WnkInterval; //旋转因子旋转相应角度.
}
ButterFlyCnt=ButterFlyCnt+SeGBTF*2;
}
}
/*********以下程序将倒位序的数据顺序排列*********/
//借用了变量LCnt,ButterFlySCnt,SCnt
LCnt=10-FFT_Level;
for(ButterFlySCnt=0;ButterFlySCnt<FFT_HFNum;ButterFlySCnt++)
{FFT_Seg=BitReversedTable1[ButterFlySCnt];
FFT_Seg=FFT_Seg>>LCnt; //比1024点小的倒位序表可由词表直接右移LCnt位得到.
if(FFT_Seg>ButterFlySCnt) //保证只交换一次
{
SCnt=2*FFT_Seg;//得到要与他相交换的元素序号的地址.
Rx1buff=X[SCnt];
Ix1buff=X[SCnt+1];
X[SCnt]=X[2*ButterFlySCnt];
X[SCnt+1]=X[2*ButterFlySCnt+1];
X[2*ButterFlySCnt]=Rx1buff;
X[2*ButterFlySCnt+1]=Ix1buff;
}
}
/*由于用N/2点计算N点FFT,这里由X[K]=X1[K]+W2nK*X2[k]恢复前N/2个点的值.*/
/*由于用N/2点计算N点FFT,这里由X[K]=X1[K]+W2nK*X2[k]恢复前N/2个点的值.*/
/*X1[K]=0.5*[Y(k)+Y~(N-K)]*/
/*X2[K]=-0.5*j*[Y(K)-Y~(N -K)]*/
/*得到X[K]=0.5*[Y(K)+Y~(N-K)]-0.5*j*W2nK*[Y(K)-Y~(N-K)]*/
/*得到X[N-K]=0.5*[Y(N-K)+Y~(K)]-0.5*j*W2n(N-K)*[Y(N-K)-Y~(K)]*/
/*只计算前N/2个FFT的值,后N/2是镜频.*/
if(RealFlag)//如果是实数运算则要作以下变换
{ButterFlyCpl=FFT_HFNum*2;
#if MAC120
MAC0CF=0x00;
#endif
Rx1buff=X[0];
Ix1buff=X[1];
Rx1buff=Rx1buff>>1;
Ix1buff=Ix1buff>>1;
X[0]=Rx1buff+Ix1buff;
X[1]=0;
WnkInterval=1024/FFT_HFNum;
Wnk=WnkInterval;
for(ButterFlySCnt=2;ButterFlySCnt<FFT_HFNum;ButterFlySCnt=ButterFlySCnt+2,Wnk=Wnk+WnkInterval)/*由于X[K]和X[N-k]配对处理,故计算FFT_HFNum/2次.*/
{Rx1buff=X[ButterFlySCnt];
Ix1buff=X[ButterFlySCnt+1];
Rx2buff=X[ButterFlyCpl-ButterFlySCnt];
Ix2buff=X[ButterFlyCpl-ButterFlySCnt+1];
RSxbuff=Rx1buff+Rx2buff;
ISxbuff=Ix1buff+Ix2buff;
#if MAC120
MACW1=0;
MACW0=0;
MAC0A=WnC[Wnk];
MAC0B=ISxbuff;
MAC0A=WnS[Wnk];
MAC0B=Rx2buff-Rx1buff;
FFT_Seg= RSxbuff;
FFT_Seg=FFT_Seg>>1;
X[ButterFlySCnt]=FFT_Seg+MACW1;//化简为:Rx1+Rx2+W2nC[k]*(Ix1+Ix2)+W2nS[k]*(Rx2-Rx1)得到实部.
MACW1=0;
MACW0=0;
MAC0A=WnC[Wnk];
MAC0B=Rx2buff-Rx1buff;
MAC0A=WnS[Wnk];
MAC0B=~(ISxbuff)+1;
FFT_Seg=Ix1buff-Ix2buff;
FFT_Seg=FFT_Seg>>1;
X[ButterFlySCnt+1]=FFT_Seg+MACW1; //化简为:Ix1-Ix2+W2nC[k]*(Rx2-Rx1)-W2nS[k]*(Ix1+Ix2)得到虚部.
/*以上为得到X[K]的值.*/
MACW1=0;
MACW0=0;
MAC0A=WnC[1024-Wnk];
MAC0B=ISxbuff;
MAC0A=WnS[1024-Wnk];
MAC0B=Rx1buff-Rx2buff;
FFT_Seg= RSxbuff;
FFT_Seg=FFT_Seg>>1;
X[ButterFlyCpl-ButterFlySCnt]=FFT_Seg+MACW1;//化简为:Rx1+Rx2+W2n[N-k]*(Ix1+Ix2)+W2nS[k]*(Rx1-Rx2)得到实部.
MACW1=0;
MACW0=0;
MAC0A=WnC[1024-Wnk];
MAC0B=Rx1buff-Rx2buff;
MAC0A=WnS[1024-Wnk];
MAC0B=~(ISxbuff)+1;
FFT_Seg= Ix2buff-Ix1buff;
FFT_Seg=FFT_Seg>>1;
X[ButterFlyCpl-ButterFlySCnt+1]=FFT_Seg+MACW1;//化简为:Ix2-Ix1+W2nC[N-k]*(Rx1-Rx2)-W2nS[N-k]*(Ix1+Ix2)得到虚部.
#else
MACBUF1=WnC[Wnk];
MACBUF2=ISxbuff;
MACBUF=MACBUF1*MACBUF2;
MACBUF1=WnS[Wnk];
MACBUF2=Rx2buff-Rx1buff;
MACBUF=MACBUF+MACBUF1*MACBUF2;
MACBUF=MACBUF>>16;
MACBUF3=RSxbuff;
MACBUF3=MACBUF3>>1;
MACBUF=MACBUF+MACBUF3;
X[ButterFlySCnt]=MACBUF; //
MACBUF1=WnC[Wnk];
MACBUF2=Rx2buff-Rx1buff;
MACBUF=MACBUF1*MACBUF2;
MACBUF1=WnS[Wnk];
MACBUF2=~(ISxbuff)+1;
MACBUF=MACBUF+MACBUF1*MACBUF2;
MACBUF=MACBUF>>16;
MACBUF3=Ix1buff-Ix2buff;
MACBUF3=MACBUF3>>1;
MACBUF=MACBUF+MACBUF3;
X[ButterFlySCnt+1]=MACBUF; //
MACBUF1=WnC[1024-Wnk];
MACBUF2=ISxbuff;
MACBUF=MACBUF1*MACBUF2;
MACBUF1=WnS[Wnk];
MACBUF2=Rx1buff-Rx2buff;
MACBUF=MACBUF+MACBUF1*MACBUF2;
MACBUF=MACBUF>>16;
MACBUF3=RSxbuff;
MACBUF3=MACBUF3>>1;
MACBUF=MACBUF+MACBUF3;
X[ButterFlyCpl-ButterFlySCnt]=MACBUF; //
MACBUF1=WnC[1024-Wnk];
MACBUF2=Rx1buff-Rx2buff;
MACBUF=MACBUF1*MACBUF2;
MACBUF1=WnS[Wnk];
MACBUF2=~(ISxbuff)+1;
MACBUF=MACBUF+MACBUF1*MACBUF2;
MACBUF=MACBUF>>16;
MACBUF3=Ix2buff-Ix1buff;
MACBUF3=MACBUF3>>1;
MACBUF=MACBUF+MACBUF3;
X[ButterFlyCpl-ButterFlySCnt+1]=MACBUF; //
#endif
}
}
#if MAC120
SFRPAGE=SFRPAGBUFF;
#endif
}
/*void shizhd(void)
{
float val1=0,val2=0;
unsigned int j;
long int tempval=0,temp1=0,temp2=0;
SFRPAGE=MAC0_PAGE;
MAC0CF=0x00; //工作于乘加方式,
for(j=8;j<256;j=j+16)
{
MACW1=0;
MACW0=0;
MAC0A=source[j];
MAC0B=source[j];
MAC0A=source[j+1];
MAC0B=source[j+1];
tempval=MACW1;
tempval=tempval<<16;
tempval=tempval+MACW0;
if(j>8)
{
temp1=temp1+tempval;
}
else
{
temp2=tempval;
}
}
r1=(float)temp1/temp2;
r1=sqrt(r1);
} */⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -