fftmy.cpp

fft and convolution

// fftmy.cpp : implementation file
//

#include "stdafx.h"
#include "fft2.h"
#include "fftmy.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

/////////////////////////////////////////////////////////////////////////////
// fftmy

fftmy::fftmy()
{
}

fftmy::~fftmy()
{
}


BEGIN_MESSAGE_MAP(fftmy, CStatic)
	//{{AFX_MSG_MAP(fftmy)
		// NOTE - the ClassWizard will add and remove mapping macros here.
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// fftmy message handlers

#include "stdafx.h"
#include "math.h"


fushu fftmy:: fucheng(fushu a,fushu b)
{
	double shibu;
	double xubu;
	fushu result;
	shibu=a.real*b.real-a.image*b.image;
	xubu=a.real*b.image+a.image*b.real;
    result.real=shibu;
	result.image=xubu;
	return result;
}
fushu fftmy:: fujia(fushu a,fushu b)
{
	double shibu;
	double xubu;
	fushu result;
	shibu=a.real+b.real;
	xubu=a.image+b.image;
    result.real=shibu;
	result.image=xubu;
	return result;

} 
fushu fftmy:: fujian(fushu a,fushu b)
{
	b.real=-b.real;
	b.image=-b.image;
	return fujia(a,b);

}

void fftmy:: floorll(int L,int N,fushu A[])
{
	int LE=1;
	int l=0;
	double pi=acos(-1);
	while(l<L)
	{
		LE=2*LE;
		l++;
	}
	int LE1=LE/2;
	fushu U;
	fushu T;
	U.real=1.0;
	U.image=0.0;
	fushu W;
	W.real=cos(pi/(double)LE1);
	W.image=-sin(pi/(double)LE1);
	int I=0;
	int J=0;
	int IP;

	while(J<LE1){
	I=J;
	while((I<N)&&(I+LE1<N)){
	IP=I+LE1;
	T=fucheng(A[IP],U);
    A[IP]=fujian(A[I],T);
	A[I]=fujia(A[I],T);
	I=I+LE;
	}
	fushu temp;
	temp=fucheng(U,W);
	U=temp;
	J++;
	}

}


int fftmy:: nextinvorder(int J,int N)
{
	int NV2;
	NV2=N/2;
	int K=NV2;
	int J1=J;
    while(J1>=K)
	{
		J1=J1-K;
		K=K/2;

	}	
	J1=J1+K;
	return J1;

}///////N 为总数，J 为当前的逆序数，返回下一个逆序数（J<N-1不能取N－1）

void  fftmy:: reorder(fushu  A[],int N)
{

fushu T;
 int I=0;
int J=0;
while(I<N-1){/////////注意x(0)和x(N-1)不需要调动，雷道算法之精妙处
if(I<J)
{
	T=A[J];
	A[J]=A[I];
	A[I]=T;
}

J=nextinvorder(J,N);
I++;
}

}

void fftmy:: FFTmy(fushu A[],int N)
{
double  m=log(N)/log(2);
int M=(int)(m+0.1);
/*
CString str;

for(int i=1;i<100;i++)
{
	int test=1;
	for(int j=0;j<i;j++)
		test*=2;
str.Format(str+" %d ",(int)((log(test)/log(2))+0.1));
}
AfxMessageBox(str);
*/
int L=1;
reorder(A,N);
while(L<=M)
{
	floorll(L,N,A);
	L++;
}
	
}

void fftmy:: fuabs(fushu A[] , int N , double B[])
{
	for(int i=0;i<N;i++){
	B[i]=sqrt(A[i].real*A[i].real+A[i].image*A[i].image);
	}
}


void fftmy:: inver(fushu A[],int N){

for(int i=0;i<N;i++)
{
	A[i].image=-A[i].image;
}
}

void fftmy:: convolution(fushu A[],fushu B[],fushu Y[],int N){


	
	FFTmy(A,N);
	FFTmy(B,N);

	for(int i=0;i<N;i++)
	{
		Y[i]=fucheng(A[i],B[i]);
	
	}

	inver(Y,N);

	for(int j=0;j<N;j++)
	{
		Y[j].real=Y[j].real/N;
			Y[j].image=Y[j].image/N;
	}

	FFTmy(Y,N);

	inver(Y,N);


}

void  fftmy::test()
{
	fushu a[2];
	fushu b[2];
	a[0].real=a[1].real=1;
	b[0].real=b[1].real=1;
	a[0].image=a[1].image=0;
	b[0].image=b[1].image=0;

	fushu c[2];
	convolution(a,b,c,2);

}

void fftmy::gettimecost(){


}

void fftmy:: convolutioncommon(fushu A[],fushu B[],fushu Y[],int N,int M){

int N1=N+1;
int N2=M+1;
fushu *X= new fushu[N1];
fushu *H = new fushu[N2];
fushu *y= new fushu[N1+N2];
for(int i=0;i<N;i++)
    X[i+1]=A[i];
for(i=0;i<M;i++)
    H[i+1]=B[i];


				i = 1;
                while(i <= N1+N2) 
                {
                    y[i].real = 0.0;
                    for(int j = 1; j <= N1+N2; j++)
                    {
                        double mm;
                        if(i - j > N2 || i - j < 1)
                            mm = 0.0;
                        else
                            mm = H[i - j].real;
                        double mmm;
                        if(j > N1)
                            mmm = 0.0;
                        else
                            mmm = X[j].real;
                        y[i].real += mm * mmm;
                    }

                    i++;
                }

for(i=0;i<N+M;i++)
  Y[i]=y[i+1];
	

}


//////////////////////////////////////