利用Intel IPP函数库实现信号频谱计算

    Intel IPP（Intel Integrated Performance Primitives）函数库是一套跨平台的软件函数库，它为用户提供了一套高效、实用的函数集，可用于实现通信、图像、语音等多个数字信号处理领域，关于该函数库的介绍，可查询其官方使用手册或以下几篇博客

Intel IPP的基本使用方法

Intel IPP函数

    本文介绍如何利用Intel函数库实现对信号进行频谱计算。

一、算法原理

    本文中，信号频谱的计算主要采用Bartlett法进行信号频谱的估计，其基本思想是对信号进行分段FFT计算，然后将分段FFT的计算结果进行累加，从而到达平滑信号频谱的效果；

二、函数介绍

    计算信号的频谱，主要用到IPP函数库中的快速傅里叶变换（FFT）相关函数和一些常见的乘、加等基本运算。

    计算FFT，主要用到ippsFFTInitAlloc_C（）和ippsFFTFwd_CToC（）两个函数，其中ippsFFTInitAlloc_C（）函数用于初始化变量空间和开辟相关缓存，ippsFFTFwd_CToC（）函数用于实现FFT的计算。相关函数的介绍和输入输出接口定义可在官方文档ippsman.pdf中找到。

三、调用示例

    下面的代码展示了如何使用IPP相关函数实现对一实信号频谱的计算：

int   CalculateSpectrumByReal(int* rData, int samplePoint, float* specturmData,	int &outSpecturmDataLength, 
	int winOpt, int fftOrder)/*实数计算频谱*/
{
	IppStatus status;
	int segement,fftPoint;
	IppsFFTSpec_R_32f* spec;
	
	fftPoint = (int)pow(2.0,fftOrder);        //FFT 点数
	outSpecturmDataLength = fftPoint/2;       //输出数据长度
	if (fftOrder < 10)
	{
		return -2; //FFT阶数设置过小,FFT阶数至少要求10阶
	}
	else if (samplePoint < fftPoint)
	{
		return -1; //数据点数不够
	}

	Ipp32f *fData = ippsMalloc_32f(samplePoint);//转为浮点数
	Ipp32f *x = ippsMalloc_32f(fftPoint);      
	Ipp32f *X = ippsMalloc_32f(fftPoint+2);    
	Ipp32f *mag = ippsMalloc_32f(outSpecturmDataLength);    //求abs()后的幅度数据
	memset(fData,0,sizeof(Ipp32f)*samplePoint);
	memset(x,0,sizeof(Ipp32f)*fftPoint);
	memset(X,0,sizeof(Ipp32f)*(fftPoint+2));
	memset(mag,0,sizeof(Ipp32f)*outSpecturmDataLength);

	memset(specturmData,0,sizeof(Ipp32f)*outSpecturmDataLength);
	
	ippsConvert_32s32f(rData,fData,samplePoint);
	status = ippsFFTInitAlloc_R_32f(&spec,fftOrder,IPP_FFT_DIV_INV_BY_N,ippAlgHintNone);

	segement = (int)(samplePoint/fftPoint);  //做FFT时的分段数
	for (int i=0;i<segement;i++)
	{		
		memcpy(x,fData+i*fftPoint,fftPoint*sizeof(Ipp32f));
		status = ippsFFTFwd_RToCCS_32f(x,X,spec,NULL);  //do fft...
		ippsMagnitude_32fc((Ipp32fc*)X,mag,outSpecturmDataLength); //abs()...
		ippsAdd_32f_I(mag,specturmData,outSpecturmDataLength);		
	}
	ippsDivC_32f_I(segement,specturmData,outSpecturmDataLength);//求平均
	ippsLn_32f_I(specturmData,outSpecturmDataLength);   //求Ln，因为IPP函数不支持浮点的10Log10计算这里通过Ln进行转换计算
	ippsDivC_32f_I(2.3026,specturmData,outSpecturmDataLength);// log(e,10)=2.3026
	ippsMulC_32f_I(10.0,specturmData,outSpecturmDataLength);  //10log10

	ippsFFTFree_R_32f(spec);
	ippsFree(fData);
	ippsFree(x);
	ippsFree(X);
	ippsFree(mag); 

	return status;
}