fastfourier.cpp

这是一个非常经典的FourierTransform 算法比较有新意

// FastFourier.cpp: implementation of the FastFourier class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "FastFourier.h"

const double PI = 4*atan(1);

void initialize(Complex *Z, int n)
{
	int k;
	Z[0] = 1;
	if(n == 1)
		return;
	else if(n == 2)
	{
		Z[1].m_dReal = -1;
		Z[1].m_dImage = 0;
	}
	else
	{
		int Quarter = n/4;
		Complex w(cos(2*PI/n),-sin(2*PI/n));
		for(k = 1; k < Quarter; k++)
		{
			Z[k] = Z[k-1]*w; 
		}
	
		for(k = 0; k < Quarter; k++)
		{
			Z[k+Quarter].m_dImage = -Z[k].m_dReal; 
			Z[k+Quarter].m_dReal = Z[k].m_dImage;
		}
	}
}


//注意：传递参数时，必须为与形参C对应的实参分配N个Complex空间
void FastFourierTransform(int m, Complex(*pf)(double), Complex *C)
{
	if(m<0)
	{
		cout << "结点个数必须为2的非负幂" << endl;
		exit(1);
	}
	int n, j, k;
	int N = (int)pow(2, m);
	
	Complex *D = new Complex[N];
	Complex *Z = new Complex[N];
	Complex *pTemp;
	
	initialize(Z, N);

	for(k = 0; k <= N-1; k++)
		C[k] = (*pf)(2*PI*k/N);

	
	for(k = 0; k <= N-1; k++)
		Z[k] /= 2;


	Complex u, v;
	
	int halfN = N/2;
	int NextGroupNum = halfN;
	int HalfGroupSize = 1;
	int const &PairGap = HalfGroupSize;
	
	int CurrentIndex1, CurrentIndex2;
	int NextIndex1, NextIndex2;
	int CofficientIndex;
	for(n = 0; n<m; n++)
	{
		CofficientIndex = 0;
		CurrentIndex1 = 0; CurrentIndex2 = CurrentIndex1+halfN;
		NextIndex1 = 0; NextIndex2 = NextIndex1+PairGap;
		for(k = 0; k < NextGroupNum; k++)
		{
			for(j = 0; j < HalfGroupSize; j++)
			{
				u = C[CurrentIndex1++] / 2;
				v = Z[CofficientIndex]*C[CurrentIndex2++];
				D[NextIndex1] = u + v;
				D[NextIndex2] = u - v;
				if((NextIndex1+1)%HalfGroupSize)
				{
					NextIndex1++;
					CofficientIndex = CofficientIndex+ NextGroupNum;
				}
				else
				{
					NextIndex1 = NextIndex1+HalfGroupSize+1; 
					CofficientIndex = 0;
				}
				NextIndex2 = NextIndex1+PairGap;
			}
		}

		NextGroupNum /= 2;
		HalfGroupSize *= 2;

		pTemp = D;
		D = C;
		C = pTemp;
	}

	if(m%2)
	{
		pTemp = D;
		D = C;
		C = pTemp;
		for(k = 0; k < N; k++)
		{
			C[k] = D[k];
		}
	}

	delete []D;
	delete []Z;
}