rceps.cpp

基于TMS320C6713B的最小相位法复倒谱算法实现

#include <stdio.h>
#include <math.h>

#define		PI					3.14159265358979f

//信号点数
#define		FFT_NUM					8

//FFT级数
#define		FFT_M					3


//复数结构
typedef struct
	{
		float real ;
		float image ;
	} COMPLEX ;								



/*=====================================================
函数名称：cbrev
函数功能：码位反转
输入参数: x 复数点序列,
		  m FFT的级数
输出参数：无,原址操作
====================================================*/
void cbrev(COMPLEX *input, int x);

/*=====================================================
函数名称：FFT
函数功能：FFT变换，基2DIT
输入参数: x 复数点序列,
		  m FFT的级数
输出参数：无,原址操作
====================================================*/
void FFT(COMPLEX *input, int x);

/*=====================================================
函数名称：com_abs_log
函数功能：log(abs(input))
输入参数: x 复数点序列,		  
输出参数：
====================================================*/
void com_abs_log(COMPLEX *input,int x);

/*=====================================================
函数名称：com_gonge
函数功能：共扼
输入参数: x 复数点序列,		  
输出参数：
====================================================*/
void com_gonge(COMPLEX *input,int x);

/*=====================================================
函数名称：com_real
函数功能：实部
输入参数: x 复数点序列,		  
输出参数：
====================================================*/
void com_real(COMPLEX *input,float *output,int x);

float cSOUECEData[FFT_NUM]={
62.97,
8.14,
14.5,
13.2,
0.1,
12.3,
5.6,
48.2};			//用来保存原始数据
COMPLEX cFFTData[FFT_NUM];			//用来保存复数的原始数据
float cCCEPSData[FFT_NUM];			//用来保存复倒谱

/*
 *  ======== main ========
 */
int main()
{
	int       i;	
	
	for (i = 0; i < FFT_NUM; i++) {
        cFFTData[i].real = cSOUECEData[i];  
		cFFTData[i].image = 0;
    }

	//cbrev(cFFTData,FFT_M);
	FFT(cFFTData,FFT_M);

	com_abs_log(cFFTData,FFT_NUM);

	com_gonge(cFFTData,FFT_NUM);
	FFT(cFFTData,FFT_M);
	com_gonge(cFFTData,FFT_NUM);
	for(i=0;i<FFT_NUM;i++)
	{
		cFFTData[i].real=cFFTData[i].real/FFT_NUM;
		cFFTData[i].image=cFFTData[i].image/FFT_NUM;
	}
	
	/*
	for(i=1;i<FFT_NUM;i++)
	{
		cFFTData[i].real=cFFTData[i].real*2;
		cFFTData[i].image=cFFTData[i].image*2;
	}
	*/
	
	//com_real(cFFTData,cCCEPSData,FFT_NUM);

	//输出复倒谱
	printf("输出复倒谱;\n");
	for (i = 0; i < FFT_NUM; i++) {
		printf("cCCEPSData[%d]=%f+j%f;\n",i,cFFTData[i].real,cFFTData[i].image);
	}
	
	return 0;
}

void cbrev(COMPLEX *input, int x)
{
	int n , i , nv2 , j , k , nm1 ;
	COMPLEX t;

	//2^M=N
	n = 1;
	for(i=0; i<x; i++)
		n = n*2 ;

	nv2 = n / 2 ;
	nm1 = n - 1 ;
  	j = 1 ;


	for (i = 1 ; i <= nm1 ; i ++)
	{
		if (i < j)
		{
			//input[i - 1]<->input[j - 1]
			t.real = input[i - 1].real ;
			t.image = input[i - 1].image ;
			input[i - 1].real = input[j - 1].real ;
			input[i - 1].image = input[j - 1].image ;
			input[j - 1].real = t.real ;
			input[j - 1].image = t.image ;
		}

		k = nv2 ;

		while (k < j)
		{
			j -= k ;
			k /= 2 ;
		}
		j += k ;
	}
}

void FFT(COMPLEX *input, int x)
{
	int n , i ,  j , le , l , le1 , ip  ;
	COMPLEX t , u , w ;

	n = 1;
	for(i=0; i<x; i++)
		n = n*2 ;


  	j = 1 ;

	cbrev(input,x);

	le = 1 ;
	for (l= 1 ; l <= x ; l ++)
	{
		le *= 2 ;
   		le1 = le / 2 ;
		u.real = 1.0f ;
		u.image = 0.0f ;
		w.real = (float) cos(PI / le1) ;
		w.image =(float) -sin(PI / le1) ;

		for (j = 1 ; j <= le1 ; j ++)
		{
			for (i = j ; i <= n ; i += le)
			{
				ip = i + le1 ;
				t.real = input[ip - 1].real * u.real - input[ip - 1].image * u.image ;
				t.image = input[ip - 1].real * u.image + input[ip - 1].image * u.real ;
				input[ip - 1].real = input[i - 1].real - t.real ;
				input[ip - 1].image = input[i - 1].image - t.image ;
				input[i - 1].real = t.real + input[i - 1].real ;
				input[i - 1].image = t.image + input[i - 1].image ;
			}

			t.real = u.real * w.real - u.image * w.image ;
			t.image = u.image * w.real + u.real * w.image ;
			u.real = t.real ;
      		u.image = t.image ;
		}
	}
}


void com_abs_log(COMPLEX *input,int x)
{
	int i;
	for(i=0;i<x;i++)
	{
		input[i].real=input[i].real*input[i].real+input[i].image*input[i].image;
		input[i].real=(float)log((double)input[i].real);
		input[i].real=input[i].real/2;
		input[i].image=0.0;
	}
}

void com_gonge(COMPLEX *input,int x)
{
	int i;
	for(i=0;i<x;i++)
	{
		input[i].image=-input[i].image;
	}
}

void com_real(COMPLEX *input,float *output,int x)
{
	int i;
	for(i=0;i<x;i++)
	{
		output[i]=input[i].real;
	}
}