fft.c

用单片机做的音频功率谱分析

#include <reg51.h>
#include <absacc.h>
#define  uchar  unsigned  char
#define  uint unsigned int
#define  LCDCOM  XBYTE [0xa000]	  // 命令/状态寄存器地址
#define  LCDDATA  XBYTE [0xa001] //	 数据寄存地址
#define  LCDSTATE  XBYTE [0xa002]
#define  LCDRAMDATA  XBYTE [0xa003]

#define  AD574_START XBYTE[0x9000]
#define  AD574_8 XBYTE[0x9002]
#define  AD574_4 XBYTE[0x9003]

#define  AMP  XBYTE[0xb000]

#define  N 256 //fft点数
#define  Dis_scale 8
sbit P32 = P3^2;
sbit k_up = P3^3;
sbit k_down = P3^4;
sbit k_spec = P3^5;
/**************函数声明*****************/
void LCD_WCom(uchar);
void LCD_WPos(uchar,uchar);
void LCD_Init();
void LCD_ClearLCD();
void LCD_WData(uchar);
void LCD_Delay();
void LCD_WString(uchar *string);
uchar LCD_IsBusy();
void LCD_DisNum(float a);
void LCD_DisFreNo(int a);
void LCD_DisFloat(float a);

void BubbleSort(float* pData,int Count);
void draw_spec();

//void count_sample();
void key_process();

uint UnFinish = 1;
uint Sample_cnt = 0;
uchar key_pause = 0;
int key_dis_i = 0;
unsigned char dec = 0;
float err = 0;
float err_add = 0;
unsigned char amp = 0;
unsigned char DataHigh = 0, DataLow = 0;
unsigned int ADResult = 0;
unsigned char xdata Ad8[256];
unsigned char xdata Ad4[256];
unsigned char xdata buffer[64][16];

int xdata data_nm[256];
float xdata re[256];
float xdata im[256];
float code RSin[128] = {1.000000,0.999699,0.998795,0.997290,0.995185,0.992480,0.989177,0.985278,0.980785,0.975702,
					0.970031,0.963776,0.956940,0.949528,0.941544,0.932993,0.923880,0.914210,0.903989,0.893224,
					0.881921,0.870087,0.857729,0.844854,0.831470,0.817585,0.803208,0.788346,0.773010,0.757209,
					0.740951,0.724247,0.707107,0.689541,0.671559,0.653173,0.634393,0.615232,0.595699,0.575808,
					0.555570,0.534998,0.514103,0.492898,0.471397,0.449611,0.427555,0.405241,0.382683,0.359895,
					0.336890,0.313682,0.290285,0.266713,0.242980,0.219101,0.195090,0.170962,0.146730,0.122411,
					0.098017,0.073565,0.049068,0.024541,0.000000,-0.024541,-0.049068,-0.073565,-0.098017,-0.122411,
					-0.146730,-0.170962,-0.195090,-0.219101,-0.242980,-0.266713,-0.290285,-0.313682,-0.336890,-0.359895,
					-0.382683,-0.405241,-0.427555,-0.449611,-0.471397,-0.492898,-0.514103,-0.534998,-0.555570,-0.575808,
					-0.595699,-0.615232,-0.634393,-0.653173,-0.671559,-0.689541,-0.707107,-0.724247,-0.740951,-0.757209,
					-0.773010,-0.788346,-0.803208,-0.817585,-0.831470,-0.844854,-0.857729,-0.870087,-0.881921,-0.893224,
					-0.903989,-0.914210,-0.923880,-0.932993,-0.941544,-0.949528,-0.956940,-0.963776,-0.970031,-0.975702,
					-0.980785,-0.985278,-0.989177,-0.992480,-0.995185,-0.997290,-0.998795,-0.999699
				};
float code RCos[128] = {0.000000,-0.024541,-0.049068,-0.073565,-0.098017,-0.122411,-0.146730,-0.170962,-0.195090,-0.219101,
					-0.242980,-0.266713,-0.290285,-0.313682,-0.336890,-0.359895,-0.382683,-0.405241,-0.427555,-0.449611,
					-0.471397,-0.492898,-0.514103,-0.534998,-0.555570,-0.575808,-0.595699,-0.615232,-0.634393,-0.653173,
					-0.671559,-0.689541,-0.707107,-0.724247,-0.740951,-0.757209,-0.773010,-0.788346,-0.803208,-0.817585,
					-0.831470,-0.844854,-0.857729,-0.870087,-0.881921,-0.893224,-0.903989,-0.914210,-0.923880,-0.932993,
					-0.941544,-0.949528,-0.956940,-0.963776,-0.970031,-0.975702,-0.980785,-0.985278,-0.989177,-0.992480,
					-0.995185,-0.997290,-0.998795,-0.999699,-1.000000,-0.999699,-0.998795,-0.997290,-0.995185,-0.992480,
					-0.989177,-0.985278,-0.980785,-0.975702,-0.970031,-0.963776,-0.956940,-0.949528,-0.941544,-0.932993,
					-0.923880,-0.914210,-0.903989,-0.893224,-0.881921,-0.870087,-0.857729,-0.844854,-0.831470,-0.817585,
					-0.803208,-0.788346,-0.773010,-0.757209,-0.740951,-0.724247,-0.707107,-0.689541,-0.671559,-0.653173,
					-0.634393,-0.615232,-0.595699,-0.575808,-0.555570,-0.534998,-0.514103,-0.492898,-0.471397,-0.449611,
					-0.427555,-0.405241,-0.382683,-0.359895,-0.336890,-0.313682,-0.290285,-0.266713,-0.242980,-0.219101,
					-0.195090,-0.170962,-0.146730,-0.122411,-0.098017,-0.073565,-0.049068,-0.024541
				};



float tempRe;
float tempIm;
float sinVal = 0;
float cosVal = 0;
float PTotal = 0;



main()
{
	uint s=0;
	uint group_intf = 0;
	uint num =0;
	uint i = 0;
	uint indexA  =0;
	uint indexB = 0;
	uint index_Rotate = 0;
	uint group = 0;
	uint stage = 0;
	uint g =0;
	uint j = 0;
	LCD_Init();
	LCD_ClearLCD();
//	AMP = 0x01;
	LCD_WPos(1,3);
	LCD_WString("WELCOME");
	LCD_WPos(3,2);
	LCD_WString("音频信号分析");
	LCD_WPos(4,3);
	LCD_WString("Initing...");

//----------------------------------循环起点-------------------------------------------------//
start:	
//------------------------设置定时器、中断寄存器，开中断，等待中断采集数据-------------------
	TMOD = 0x02;
	TL0 =  152; 
	TH0 = 152;
	ET0 = 1;
	TR0 = 1;
	EX0 = 0;
	EA = 1;	
	while(UnFinish);
	EX0 = 1;
//--------------拼接AD采样的高低位-----------------//
	for(i = 0; i < 256; i++)
	{
		tempRe = (float)Ad8[i] * 16;
		tempIm = (float)Ad4[i] / 16;
		im[i] = tempRe + tempIm;
//		re[i] = im[i];
	}
	im[0] = im[1];

//********************显示采样结果***************//
/*	for(i = 0; i <128; i++)
	{
		LCD_WPos((i)%4+1,1);
		LCD_DisFreNo(i);
		LCD_WString(": ");
		LCD_DisNum(im[i]);
		LCD_Delay();LCD_Delay();LCD_Delay();LCD_Delay();LCD_Delay();
	}
*/

//***********倒序，AD后数据变为正弦数据***************//
	for(i = 0; i < 256; i++)
	{
		indexA = 0;
		indexB = i;
		for(j = 0; j < 8; j++)
		{
			
			indexA = indexA << 1;
			if(0x01 & indexB)
				indexA += 1;
			indexB = indexB >> 1;
		}
		//if(amp == 0) 
		re[i] = (im[indexA]-2048)*0.009574;//5/1024*1.9608
//		else if(amp == 1) re[i] = (im[indexA]-2048)*5/1024/4;
	}
	for(i= 0; i < 256; i++)
	{
		im[i] = 0;
	}
//--------------幅值判断放大------------------------//
/*	BubbleSort(im,256);
	if(im[1] < 1.0)
	{
		amp = 1;	
		AMP = 0x03;
		goto start;
	}
*/
//------------------------数据采集结束，频谱分析------------------------------//

//*****************fft算法*************************//
	stage = 8;
	for(s = 0; s < stage; s++)		//级
	{
		group = N/2/(int)(0x01<<s);		//组
		group_intf = N / group;		//组间下标间隔

		for(g = 0; g < group; g++)
		{
			num = (int)(0x01 << s);		//组内蝶形数
			
			for (i = 0; i < num; i++)
			{
				indexA = g * group_intf + i; 
				indexB = indexA + (int)(0x01<<s);	//单个蝶形两个下标
				index_Rotate = i * (int)(0x01<<(stage-s-1));		//旋转因子下标
				sinVal = RSin[index_Rotate];
				cosVal = RCos[index_Rotate];
				
				tempRe = (re[indexB] * sinVal - im[indexB] * cosVal);
				tempIm = (re[indexB] * cosVal + im[indexB] * sinVal);

				re[indexB] = re[indexA] - tempRe;
				im[indexB] = im[indexA] - tempIm;
				re[indexA] = re[indexA] + tempRe;
				im[indexA] = im[indexA] + tempIm;

 			}
		}
	}
//********************计算功率谱**********************//
	
	for(i = 0; i < 128; i++)
	{
		re[i] = (re[i]*re[i] + im[i]*im[i])/128/50;
	}
	re[0] = 0;
	for(i = 101; i < 128; i++)
		re[i] = 0;

	PTotal = 0;
	for(i = 0; i < 128; i++)
	{
		PTotal = PTotal + re[i];
	}
	
	dec = 0;
	for(i = 0; i < 101; i++)
	{
	 	if(re[i] > 0.2*PTotal)
		{
			re[i] = re[i] + i * re[i]/350;
			if(re[i-1] < 0.1 * PTotal)
				re[i-1] = 0.1*re[i-1];
			if(re[i+1] < 0.1 * PTotal)
				re[i+1] = 0.1*re[i+1];
		}
	}
//--------------------人为去掉2小功率值 *_* ------------
//	for(i = 0; i < 128; i++)
//	{
//		data_nm[i] = i;
//	}
//	BubbleSort(re,128);
	dec = 0;
	for(i = 0; i < 128; i++)
	{
		if( re[i] < 0.1 * PTotal )
		{
   	    	if(	dec <= 1 )
			{
				re[i] = 0;
				dec++;
			}
		}
	}
//--------------------功率排序-------------------------
	for(i = 0; i < 128; i++)
	{
		data_nm[i] = i;
	}
	BubbleSort(re,128);

//-------------------计算总功率------------------------
	PTotal = 0;
	for(i = 0; i < 128; i++)
	{
		PTotal += re[i];
	}
//-------------------计算失真度------------------------
	err = 0;
	err_add = 0;
	for(i = 1; i < 100; i++)
	{
		err_add = err_add + re[i];
	}
	err = err_add / re[0];
//********************显示功率计算结果***************//
	LCD_WPos(1,1);
	LCD_WString("n 频率  功率/mW");

	for(i = 0; i <3; i++)
	{
		if(key_pause == 1)
		{
			key_process();
		}
		LCD_WPos(1,1);
		LCD_WString("n 频率  功率/mW");

		LCD_WPos(2,1);
		LCD_WString("总值: ");
		LCD_DisFloat(PTotal);

		LCD_WPos(3,1);
		LCD_WString("1)");
		LCD_DisFreNo(data_nm[0]);
		LCD_WString(": ");
		LCD_DisFloat(re[0]);

		LCD_WPos(4,1);
		LCD_WString("2)");
		LCD_DisFreNo(data_nm[1]);
		LCD_WString(": ");
		LCD_DisFloat(re[1]);
//		LCD_Delay();LCD_Delay();LCD_Delay();LCD_Delay();LCD_Delay();LCD_Delay();
	}


//********************回 　’start ’ 点重复采样、变换过程************//
//	amp = 0;
//	AMP = 0x01;
	goto start;
}




//---------------------------------------------------------------------//
void LCD_WCom(uchar instru)
{
	while(LCD_IsBusy());
	LCDCOM = instru;
}

void LCD_WPos(uchar line, uchar col)
{
	uchar pos = (0x80|((line-1)%2<<4))|line/3<<3|(col-1);
	while(LCD_IsBusy());
	LCDCOM = pos;
}

void LCD_Init()
{
	while(LCD_IsBusy());
	LCD_WCom(0x30);
	while(LCD_IsBusy());
	LCD_WCom(0x0b);
	while(LCD_IsBusy());
	LCD_WCom(0x14);
	while(LCD_IsBusy());
	LCD_ClearLCD();
	while(LCD_IsBusy());
	LCD_WCom(0x0c);
}

void LCD_ClearLCD()
{
	while(LCD_IsBusy());
	LCDCOM = 0x01;
}

 void LCD_WData(uchar dat)
 {
	while(LCD_IsBusy());
	LCDDATA = dat;
 }

uchar LCD_IsBusy()
 {
	 uchar BF_AC = 0;
	 BF_AC = LCDSTATE;
	 if(BF_AC & 0x80)
	 {
		 return 1;
	 }
	 else
	 {
		 return 0;
	 }
 } 

 void LCD_Delay()
 {
	uchar i,j;
	for(i = 0; i < 255; i++)
		for(j = 0; j < 255; j++);
 }

 void LCD_WString(uchar *string)
 {
 	uchar i = 0;
	while(string[i]!=0)
	{
		LCD_WData(string[i]);
		i++;
	}
 } 
void LCD_DisNum(float a)
{
	uchar temp1 = 0;
	uchar temp2 = 0;
	uchar temp3 = 0;
	uchar temp4 = 0;
	uchar temp5 = 0;
	uchar temp6 = 0;
	uchar temp7 = 0;
	unsigned long aa = (long)(a);
	temp1 = (aa)/1000000;
	temp2 = (aa - temp1*1000000)/100000;
	temp3 = (aa - temp1*1000000 - temp2*100000)/10000;
	temp4 = (aa - temp1*1000000 - temp2*100000 - temp3*10000) / 1000;
	temp5 = (aa - temp1*1000000 - temp2*100000 - temp3*10000 - temp4*1000) / 100;
	temp6 = (aa - temp1*1000000 - temp2*100000 - temp3*10000 - temp4*1000 - temp5*100) / 10;