fft.c

2007年全国赛题A题音频信号分析仪的源码

#include "msp430x22x4.h"

#include "math.h"

//////////////////////////////////////////////////////////////
//快速傅利叶变换的源码
/*
This computes an in-place complex-to-complex FFT 
x and y are the real and imaginary arrays of 2^m points.
dir = -1 gives forward transform
dir = 1 gives reverse transform 
*/
//////////////////////////////////////////////////////////////

void FFT(signed char dir,int m,float*z1,float*z2)//m=num
{
	unsigned char   n,i,i1,j,k,i2,l,l1,l2;
	float  c1,c2,tx,ty,t1,t2,u1,u2,z;
	/* Calculate the number of points */
	n = 1;
	for (i=0;i<m;i++) 
		n *= 2;
	/* Do the bit reversal */
	i2 = n >> 1;
	j = 0;
	//64=n;
	for (i=0;i<n-1;i++) 
	{
		if (i < j)
		{
			tx = z1[i];
			ty = z2[i];
			z1[i] = z1[j];
			z2[i] = z2[j];
			z1[j] = tx;
			z2[j] = ty;
		}
		k = i2;
		while (k <= j) 
		{
			j -= k;
			k >>= 1;
		}
		j += k;
	}

	/* Compute the FFT */
	c1 = -1.0; 
	c2 = 0.0;
	l2 = 1;
	for (l=0;l<m;l++) {
		l1 = l2;
		l2 <<= 1;
		u1 = 1.0; 
		u2 = 0.0;
		for (j=0;j<l1;j++)
		{
			for (i=j;i<n;i+=l2) 
			{
				i1 = i + l1;
				t1 = u1 * z1[i1] - u2 * z2[i1];
				t2 = u1 * z2[i1] + u2 * z1[i1];
				z1[i1] = z1[i] - t1; 
				z2[i1] = z2[i] - t2;
				z1[i] += t1;
				z2[i] += t2;
			}
			z = u1 * c1 - u2 * c2;
			u2 = u1 * c2 + u2 * c1;
			u1 = z;
		}
		c2 = sqrt((1.0 - c1) / 2.0);
		//if (dir == 1) 
			//c2 = -c2;
		c1 = sqrt((1.0 + c1) / 2.0);
	}
	/* Scaling for forward transform */
	/*if (dir == 1)
	{
		for (i=0;i<n;i++)
		{
			z1[i] /= n;
			z2[i] /= n;
		}
	}*/
	for(i=0;i<n;i++)
			z2[i]=-z2[i];
	
}

float Wpower(float*m1,float*m2)//定点功率计算和功率分贝计算,m1,m2是所得数据的实部和虚部
{
	float tatal;//功率和	
	unsigned char i;
	tatal=0;
    for( i=0;i<64;i++)
	{
		m1[i]=(m1[i]*m1[i]+m2[i]*m2[i])/64/64*34;//定点功率
		m2[i]=0;
	}
	for(i=1;i<32;i++)
		tatal+=m1[i];
	//tatal_t=tatal/64*2;//功率
    return tatal;

}