f_fft.cpp

神经网络中的声速测定的源代码 学神经网络的必用程序

#include <math.h>
#include <stdlib.h>

#define PI (4*atan(1))

void fft(double *rx,double *ix,int n,int ctrl)
{
	int i,j,k,m,num,t;
	double temp1,temp2,temp3,re,im;
	if (ctrl!=0)
		ctrl=1;
//码位倒置
	temp1=n;	m=0;
	while(temp1!=1)
	{
		temp1/=2;
		m++;
	}
	for (i=0;i<n;i++)
	{
		num=0;
		t=i;
		for (j=0;j<m;j++)
		{
			num+=(t%2)*int(pow(2,(m-1-j)));
			t=(t-t%2)/2;
		}
		if(num>i)
		{
			temp1=*(rx+i);
			temp2=*(ix+i);
			*(rx+i)=*(rx+num);
			*(ix+i)=*(ix+num);
			*(rx+num)=temp1;
			*(ix+num)=temp2;
		}
	}
//蝶形运算
	for (i=0;i<m;i++)
	{
		for (j=0;j<(n/int(pow(2,(i+1))));j++)
		{
			for (k=0;k<int(pow(2,i));k++)
			{
				t=k+j*int(pow(2,(i+1)));
				re=rx[t+int(pow(2,i))]*cos(2*PI*k/pow(2,(i+1)))+pow(-1,ctrl)*ix[t+int(pow(2,i))]*sin(2*PI*k/pow(2,(i+1)));
				im=ix[t+int(pow(2,i))]*cos(2*PI*k/pow(2,(i+1)))-pow(-1,ctrl)*rx[t+int(pow(2,i))]*sin(2*PI*k/pow(2,(i+1)));
				temp1=rx[t]+re;
				temp2=rx[t]-re;
				temp3=ix[t]+im;
				ix[t+int(pow(2,i))]=ix[t]-im;
				rx[t]=temp1;
				rx[t+int(pow(2,i))]=temp2;
				ix[t]=temp3;
			}
		}
	}
	if (ctrl==1)
	{
		for (i=0;i<n;i++)
		{
			*(rx+i)/=n;
			*(ix+i)/=n;
		}
	}
	return;
}