rfft.cpp

信号处理的几种算法

#include <math.h>

void rfft(double x[64], int n)
{
	int i,j,k,m,i1,i2,i3,i4,n1,n2,n4;
	double a,e,cc,ss,xt,t1,t2;
	for(j=1,i=1;i<16;i++)
	{
		m=i;
		j=2*j;
		if(j==n) break;
	}
	n1=n-1;
	for(j=0,i=0;i<n1;i++)
	{
		if(i<j)
		{
			xt=x[j];
			x[j]=x[i];
			x[i]=xt;
		}
		k=n/2;
		while(k<(j+1))
		{
			j=j-k;
			k=k/2;
		}
		j=j+k;
	}
	for(i=0;i<n;i+=2)
	{
		xt=x[i];
		x[i]=xt+x[i+1];
		x[i+1]=xt-x[i+1];
	}
	n2=1;
	for(k=2;k<=m;k++)
	{
		n4=n2;
		n2=2*n4;
		n1=2*n2;
		e=6.28318530718/n1;
		for(i=0;i<n;i+=n1)
		{
			xt=x[i];
			x[i]=xt+x[i+n2];
			x[i+n2]=xt-x[i+n2];
			x[i+n2+n4]=-x[i+n2+n4];
			a=e;
			for(j=1;j<=(n4-1);j++)
			{
				i1=i+j;
				i2=i-j+n2;
				i3=i+j+n2;
				i4=i-j+n1;
				cc=cos(a);
				ss=sin(a);
				a=a+e;
				t1=cc*x[i3]+ss*x[i4];
				t2=ss*x[i3]-cc*x[i4];
				x[i4]=x[i2]-t2;
				x[i3]=-x[i2]-t2;
				x[i2]=x[i1]-t1;
				x[i1]=x[i1]+t1;
			}
		}
	}
 }