FFT的模版编程实现,麻烦大侠帮忙优化一下,编译时间太长了.

[zhang_db]

//对数计算
template<int N>struct Log{
	enum{value = Log<(N>>1)>::value + 1}; 
};
template<>struct Log<1>{ enum{value = 0};};
//数值反转
template<int N, int L> struct Inverse
{
	enum{value = ((N % 2) << (L - 1)) + Inverse<N / 2,L - 1>::value};
};
template<int N> struct Inverse<N,1>
{
	enum{value = (N % 2)};
};
//反转表
template<int N> struct InvPairs
{
	static int p[N];
};
template<int N>
int InvPairs<N>::p[N] = {0};

template<int n,int N> 
struct Initp
{
	//优化后,这是一个代价很小的内联函数
	void operator()(void)
	{
		InvPairs<N>::p[n] = Inverse<n,Log<N>::value>::value;
		Initp<n+1,N>()();
	}
};
template<int N>
struct Initp<N, N>
{
	void operator()(void)
	{
	}
};
//反转表生成函数
template<int N>
inline void Initps(void)
{
	Initp<0,N>()();
}

//系数表函数
template<int i, int Nn, int N>
struct IDX
{
	enum{value = i * N / Nn};
};

template<int idxParam, int idxFst,int idxSec, int N>
struct Meta_Fly
{
	template<typename T>
	static inline void Meta_Fly_Op(const complex<T> (¶m)[N/2], complex<T> (&data)[N])
	{
		complex<T> tmp(param[idxParam]);
		tmp *= data[idxSec];
		data[idxSec] = data[idxFst];
		data[idxSec]-=tmp;
		data[idxFst]+=tmp;
	}
};
///Loop 的初始值应该为 M / 2 - 1
template<int M, int N, int lvl, int Loop>
struct Mth_Fly
{
//	enum{logN = Log<N>::value};
//	enum{logM = Log<M>::value};
//	enum{logMd2 = logM / 2};
	template<typename T> 
	static inline void Mth_Fly_Op(const complex<T> (¶m)[N/2], complex<T> (&data)[N])
	{
		Meta_Fly<IDX<Loop, M, N>::value, lvl * M + Loop, lvl * M + M / 2 + Loop, N>::Meta_Fly_Op<T>(param,data);
		Mth_Fly<M,N,lvl,Loop -1>::Mth_Fly_Op<T>(param, data);
	}
};
template<int M,int N,int lvl>
struct Mth_Fly<M,N,lvl,0>
{
//	enum{logN = Log<N>::value};
//	enum{logM = Log<M>::value};
//	enum{logMd2 = logM / 2};
	template<typename T>
	static inline void Mth_Fly_Op(const complex<T> (¶m)[N/2], complex<T> (&data)[N])
	{
		Meta_Fly<IDX<0, M, N>::value, lvl * M, lvl * M + M / 2, N>::Meta_Fly_Op<T>(param,data);
	}
};
///一级蝶型运算 Loop 的初始值为N/M - 1, 终止值为 0

template<int M, int N, int Loop>
struct Mth_Fly_Full
{
	///第M级,一共存在N/M个蝶型运算.
	template<typename T>
	static inline void Mth_Fly_Full_Op(const complex<T> (¶m)[N/2], complex<T> (&data)[N])
	{
		Mth_Fly<M,N,Loop,M / 2 - 1>::Mth_Fly_Op<T>(param,data);
		Mth_Fly_Full<M,N,Loop - 1>::Mth_Fly_Full_Op(param,data);
	}
};
//终止递归的特化
template<int M,int N>
struct Mth_Fly_Full<M,N,0>
{
	template<typename T>
	static inline void Mth_Fly_Full_Op(const complex<T> (¶m)[N/2], complex<T> (&data)[N])
	{
		Mth_Fly<M,N,0, M / 2 - 1>::Mth_Fly_Op<T>(param,data);
	}
};
///FFT运算,Loop的初始值为2,即第一级的蝶型运算
template<int N, int Loop>
struct Fly_Full
{
	template<typename T>
	static inline void Fly_Full_Op(const complex<T> (¶m)[N/2], complex<T> (&data)[N])
	{
		Mth_Fly_Full<Loop,N,N/Loop-1>::Mth_Fly_Full_Op<T>(param,data);
		Fly_Full<N,Loop * 2>::Fly_Full_Op(param,data);
	}
};
///终止递归的特化
template<int N>
struct Fly_Full<N,N>
{
	template<typename T>
	static inline void Fly_Full_Op(const complex<T> (¶m)[N/2], complex<T>(&data)[N])
	{
		Mth_Fly_Full<N,N,0>::Mth_Fly_Full_Op<T>(param,data);
	}
};
template<int N,typename T>
void FlyFull(const complex<T> (¶m)[N/2], complex<T>(&data)[N])
{
	Fly_Full<N,2>::Fly_Full_Op(param,data);
}
template<int N, typename T = double>
class FFT
{
private:
	typedef complex<T> type;
	static const T pi;
	static const T _2pi;
	static const int Nd2 = N / 2;
	static type W[Nd2];
	int (&p)[N];
	type data_[N];
public:
	type result[N];
public:

	FFT() : p(InvPairs<N>::p)
	{
		Initps<N>();
		for(int i=0; i< Nd2; ++i)
		{
			W[i].real = cos( i * _2pi / N);
			W[i].image = -sin(i * _2pi / N);
		}
	}
	/*超出范围的数值补零处理.*/
	void SetData(T* data, unsigned int len)
	{
		if(len > N) len = N;
		int i;
		for(i = 0; i < len; ++i) data_[i].real = data[i];
		for(; i < N; ++i) data_[i] = 0;
	}
	void DoFFT()
	{
		for(int i=0; i< N; ++i)
		{
			result[i] = data_[p[i]];
		}
		Fly_Full<N,2>::Fly_Full_Op<T>(W,result);
	}
};
template<int N,typename T> typename FFT<N,T>::type FFT<N,T>::W[FFT<N,T>::Nd2];
template<int N,typename T> const typename T FFT<N,T>::pi = 3.1415926535897932384626433832795;
template<int N,typename T> const typename T FFT<N,T>::_2pi = 3.1415926535897932384626433832795 * 2;

那个complex<T> 是个复数的模版类,没什么特别的.需要有复数的*=, +=, -= 和赋值运算几个操作.
程序编译通过,可是FFT的数目上升到4096的时候编译器就会出错,无法编译.
那位大哥帮忙优化下,
示例测试代码:

#define SIZE 256
#include <fstream>
int _tmain(int argc, _TCHAR* argv[])
{
	double pi = 3.1415926;
	FFT<SIZE> fd;
	double d[SIZE];
	for(int i=0; i< SIZE; ++i)
	{
		d[i] = sin( i * pi / 64);
	}
	fd.SetData(d,SIZE);
	fd.DoFFT();

	std::ofstream of("result.txt",std::ios_base::out);
	for(int i=0; i< SIZE; ++i)
	{
		of << fd.result[i].real<<'\t'<<fd.result[i].image<<'\n';
	}
	of.flush();
        of.close();
}

分不够可以加.

[fallening]

如果真的想优化的话，去看fftw,这个是做fft最快的;
你的这个fft是残疾的，只能处理2^n量的数据，没有什么价值。

[K行天下]

看看FFTW库

[linglongyouzhi]

正解
[Quote=引用 2 楼 k2eats 的回复:]
看看FFTW库
[/Quote]

[taodm]

楼主，你到底要编译快还是运行快？

[hastings]

想都快，哈哈。

[xkyx_cn]

学习中