fft3.cpp

关于fft的混合基运算

#include <stdio.h>
#include <math.h>

#define PI 3.1415926

typedef struct{
	double real;
	double imag;
}COMPLEX;

COMPLEX Add(COMPLEX a, COMPLEX b);	//复数加法
COMPLEX Sub(COMPLEX a, COMPLEX b);	//复数减法
COMPLEX Mul(COMPLEX a, COMPLEX b);	//复数乘法
COMPLEX GetW(int order, int size);	//计算旋转因子
void L_Atom(COMPLEX* TD, int size);	//分裂基算子
void FFT4(COMPLEX* TD, int size);	//基4FFT算子

//复数加法
COMPLEX Add(COMPLEX a, COMPLEX b)
{
	COMPLEX temp;
	temp.real = a.real + b.real;
	temp.imag = a.imag + b.imag;
	return temp;
}

//复数减法
COMPLEX Sub(COMPLEX a, COMPLEX b)
{
	COMPLEX temp;
	temp.real = a.real + b.real;
	temp.imag = a.imag + b.imag;
	return temp;
}

//复数乘法
COMPLEX Mul(COMPLEX a, COMPLEX b)
{
	COMPLEX temp;
	temp.real = a.real * b.real - a.imag * b.imag;
	temp.imag = a.real * b.imag + a.imag * b.real;
	return temp;
}

//计算旋转因子
COMPLEX GetW(int order, int size)	
{
	COMPLEX temp;
	temp.real = cos(2*PI*order/size);
	temp.imag = sin(2*PI*order/size);
	return temp;
}

//分裂基FFT的递归算法
void L_Atom(COMPLEX* TD, int size)
{
    int i, half1, half2;
	half1 = size/2;
	half2 = size/4;

	for(i = 0; i < half1; i++)
	{
		TD[i] = Add(TD[i], TD[i+half1]);
		TD[i+half1] = Sub(TD[i], TD[i+half1]);
	}
	if(size > 2)
	{
		for(i = 0; i < half2; i++)
		{
			TD[i+half1] = Mul(Sub(TD[i+half1], TD[i+3*half2]), GetW(i, size));
			TD[i+3*half2] = Mul(Add(TD[i+half1], TD[i+3*half2]), GetW(3*i, size));
		}
		if(size > 4)
		{
			L_Atom(&TD[half1], size/4);
			L_Atom(&TD[3*half2], size/4);
		}
		L_Atom(TD, size/2);
    }
}

//基4FFT算子
void FFT4(COMPLEX* TD, int size)
{
	int i, half1, half2;
	half1 = size/2;
	half2 = size/4;

	for(i = 0; i < half1; i++)
	{
		TD[i] = Add(TD[i], TD[i+half1]);
		TD[i+half1] = Sub(TD[i], TD[i+half1]);
	}
	if(size > 2)
	{
		for(i = 0; i < half2; i++)
		{
			TD[i] = Add(TD[i], TD[i+half2]);
			TD[i+half1] = Mul(Sub(TD[i+half1], TD[i+3*half2]), GetW(i, size));
			TD[i+half2] = Mul(Sub(TD[i], TD[i+half2]), GetW(2*i, size));
			TD[i+3*half2] = Mul(Add(TD[i+half1], TD[i+3*half2]), GetW(3*i, size));
		}
		if(size > 4)
		{
			FFT4(TD, size/4);
			FFT4(&TD[half2], size/4);
			FFT4(&TD[half1], size/4);
			FFT4(&TD[3*half2], size/4);
		}
		if(size == 4)
		{
			FFT4(TD, size/2);
			FFT4(&TD[half1], size/2);
		}
    }
}