1.FFT:
// data为输入和输出的数据,N为长度
bool CFFT::Forward(complex *const Data, const unsigned int N)
{
if (!Data || N < 1 || N & (N - 1))
return false;
// 排序
Rearrange(Data, N);
// FFT计算:const bool Inverse = false
Perform(Data, N);
return true;
}
2.IFFT:
// Scale 为是否缩放
bool CFFT::Inverse(complex *const Data, const unsigned int N,
const bool Scale /* = true */)
{
if (!Data || N < 1 || N & (N - 1))
return false;
// 排序
Rearrange(Data, N);
// FFT计算,ture表示是逆运算
Perform(Data, N, true);
// 对结果进行缩放
if (Scale)
CFFT::Scale(Data, N);
return true;
}
3.排序:
void CFFT::Rearrange(complex *const Data, const unsigned int N)
{
// Swap position
unsigned int Target = 0;
// Process all positions of input signal
for (unsigned int Position = 0; Position < N; ++Position)
{
// Only for not yet swapped entries
if (Target > Position)
{
// Swap entries
const complex Temp(Data[Target]);
Data[Target] = Data[Position];
Data[Position] = Temp;
}
// Bit mask
unsigned int Mask = N;
// While bit is set
while (Target & (Mask >>= 1))
// Drop bit
Target &= ~Mask;
// The current bit is 0 - set it
Target |= Mask;
}
}
4.FFT计算:
void CFFT::Perform(complex *const Data, const unsigned int N,
const bool Inverse /* = false */)
{
const double pi = Inverse ? 3.14159265358979323846 : -3.14159265358979323846;
// Iteration through dyads, quadruples, octads and so on...
for (unsigned int Step = 1; Step < N; Step <<= 1)
{
// Jump to the next entry of the same transform factor
const unsigned int Jump = Step << 1;
// Angle increment
const double delta = pi / double(Step);
// Auxiliary sin(delta / 2)
const double Sine = sin(delta * .5);
// Multiplier for trigonometric recurrence
const complex Multiplier(-2. * Sine * Sine, sin(delta));
// Start value for transform factor, fi = 0
complex Factor(1.);
// Iteration through groups of different transform factor
for (unsigned int Group = 0; Group < Step; ++Group)
{
// Iteration within group
for (unsigned int Pair = Group; Pair < N; Pair += Jump)
{
// Match position
const unsigned int Match = Pair + Step;
// Second term of two-point transform
const complex Product(Factor * Data[Match]);
// Transform for fi + pi
Data[Match] = Data[Pair] - Product;
// Transform for fi
Data[Pair] += Product;
}
// Successive transform factor via trigonometric recurrence
Factor = Multiplier * Factor + Factor;
}
}
}
5.缩放:
void CFFT::Scale(complex *const Data, const unsigned int N)
{
const double Factor = 1. / double(N);
// Scale all data entries
for (unsigned int Position = 0; Position < N; ++Position)
Data[Position] *= Factor;
}
#include
#include
#include
typedef struct{
double r;
double i;
}my_complex
;
//检查a是否为2的整数次方数
#define NOT2POW(a) (((a)-1)&(a)||(a)<=0)
//pi
#define MYPI 3.14159265358979323846
my_complex* fft(const my_complex* x, unsigned int len){
unsigned int ex=0,t=len;
unsigned int i,j,k;
my_complex *y;
double tr,ti,rr,ri,yr,yi;
if(NOT2POW(len)) return NULL; //如果失败,返回空指针
for(;!(t&1);t>>=1) ex++; //len应该等于2的ex次方
y=(my_complex*)malloc(len*sizeof(my_complex));
if(!y) return NULL;
//变址计算,库里-图基算法
for(i=0;i
j=0;
t=ex;
while((t--)>0){
j<<=1;
j|=k&1;
k>>=1;
}
if(j>=i){
y[i]=x[j];
y[j]=x[i];
}
}
//用变址后的y向量进行计算
for(i=0;i
rr=cos(ti);
ri=sin(ti);
tr=y[j+k+t].r;
ti=y[j+k+t].i;
yr=rr*tr-ri*ti;
yi=rr*ti+ri*tr;
tr=y[j+k].r;
ti=y[j+k].i;
y[j+k].r=tr+yr;
y[j+k].i=ti+yi;
y[j+k+t].r=tr-yr;
y[j+k+t].i=ti-yi;
}
}
}
return y;
}
//以下为测试
int main()
{
int i,DATA_LEN;
my_complex *x,*y;
printf("基二FFT测试\n输入生成序列长度:");
scanf("%d",&DATA_LEN);
x=(my_complex*)malloc(DATA_LEN*sizeof(my_complex));
for(i=0;i
x[i].i=i-1;
}
printf("处理前...\n实部\t\t虚部\n");
for(i=0;i
y=fft(x,DATA_LEN);
if(!y){
printf("序列长度不为2的整数次方!\n");
return 0;
}
printf("处理后...\n实部\t\t虚部\n");
for(i=0;i
free(y);
free(x);
return 0;
}