fft.cpp

fft傅立叶快速变换在图象处理方面的应用

	bPlus=true;bNeg=true;
	//startPlus=base+offset+nCenterPlus;
	//startNeg=base+offset-nCenterNeg;

	if(meanPlus-theshold<0.00001 )
	{
		//peakPlus=0;
		startPlus=base+offset+20;
		//bPlus=false;
	}
	else
	{
		startPlus=base+offset+nCenterPlus;
		//bPlus=true;
		
	}
	if(meanNeg-theshold<0.00001)
	{
		//peakNeg=0;
		startNeg=base+offset;
		//bNeg=false;
	}
	else
	{
		startNeg=base+offset-nCenterNeg;
		//bNeg=true;
		
	}

	//3.峰值分析，点包络，从质心位置开始搜索，当前仅求正向包络，注意：算法要避免死循环
	cntPlus=0,cntNeg=0;
	if(bPlus)
	{
		//peakPlus=nCenterPlus;
		for(i=startPlus;i<255;i++)
		{
			if(cntPlus>5)
			{
				peakPlus=i-(base+offset)-cntPlus;//包络点到基线的距离
				break;
			}
			else
			{
				peakPlus=i-(base+offset);
			}
			
			if(smooth[i]-meanPlus*1.3<0.000001)//阈值如何选取
			{
				cntPlus++;
			}
			else//保证检测的是连续点
			{
				if(--cntPlus<0)
					cntPlus=0;
			}
			
		}

	}
	if(bNeg)
	{
		//for(i=base+offset;i>=/*base+offset-nCenterNeg*/0;i--)
		for(i=startNeg;i>=0;i--)
		{
			if(cntNeg>5)
			{
				peakNeg=(base+offset)-(i+cntNeg);
				break;
			}
			else
			{
				peakNeg=(base+offset)-i;
			}
			if(smooth[i]-meanNeg*1.3<0.000001)
			{
				cntNeg++;
			}
			else
			{
				if(--cntNeg<0)
					cntNeg=0;
			}
		}
	}
	//cover=((nCenterPlus&0xffff)<<8)+nCenterNeg;
	cover=((peakPlus&0xffff)<<8)+peakNeg;
	return cover;
}

/**************************************************************************************
                    包络中值滤波算法  by Jim Fang at 2007

  对前后的七个包络点和前一运算的中值点，共8点进行中值滤波
**************************************************************************************/
int WINAPI Jim_Mean(int *pTemp,int no)
{
	int tmp[8];

	int m,k,j,i,d;
	k=0; 
	m=no-1;
	for(i=0;i<8;i++)
	{
		tmp[i]=pTemp[i];
	}
	while (k<m)
	{ 
		j=m-1; m=0;
		for (i=k; i<=j; i++)
			if (tmp[i]>tmp[i+1])
			{ 
				d=tmp[i]; tmp[i]=tmp[i+1]; tmp[i+1]=d; m=i;
			}
		j=k+1; k=0;
		for (i=m; i>=j; i--)
			if (tmp[i-1]>tmp[i])
			{ 
				d=tmp[i]; tmp[i]=tmp[i-1]; tmp[i-1]=d; k=i;
			}
	}
	//取中值输出，用于滤除包络点的奇异点
	//return tmp[4];
	return int((tmp[2]+tmp[3]+tmp[4]+tmp[5])/4.0);
}

/**************************************************************************************
                    图像增强算法  by Jim Fang at 2007

	根据包络的位置，增强包络内的图像，消弱包络外的图像
	实际检测中会出现双向血流的情况，为避免误操作，需要对基线上下的能量进行对比
**************************************************************************************/
void WINAPI	Jim_Boost(void *pSrcData,void *pDstData,int coverP,int coverN,int base,int splus,int image,double dnr,int method)
{
	unsigned char *srcPtr =(unsigned char *)pSrcData;
	unsigned char *dstPtr =(unsigned char *)pDstData;
	int i,nBoostDir;
	double nZoom,temp;
	double powerPlus=0.0,powerNeg=0.0;

	for(i=base-coverN;i<base;i++)
		powerNeg+=srcPtr[i];
	for(i=base;i<base+coverP;i++)
		powerPlus+=srcPtr[i];

	//判断需要增强的是正向还是负向
	if(powerPlus/powerNeg>dnr)
		nBoostDir=1;
	else if(powerPlus/powerNeg>(1.0/dnr))
		nBoostDir=0;
	else
		nBoostDir=2;


	if(nBoostDir==1)
	{
		for(i=0;i<fft_point;i++)
		{
			if(i>base+coverP)
			{
				dstPtr[i]=srcPtr[i];
			}
			
			else if(i>=base)
			{
				nZoom=1+splus*0.1;
				//注意数制转换存在的bug
				temp=srcPtr[i]*nZoom;
				if(temp-255>0.000001)
					dstPtr[i]=255;
				else
					dstPtr[i]=unsigned char(temp);
				
			}
			else//图象增强，将包络外的能量衰减，包络内的能量不变
			{
				nZoom=1-image*0.1;
				//注意数制转换存在的bug
				temp=srcPtr[i]*nZoom;
				if(temp-255>0.000001)
					dstPtr[i]=255;
				else
					dstPtr[i]=unsigned char(temp);
				
			}
		}
	}
	else if(nBoostDir==2)
	{
		for(i=0;i<fft_point;i++)
		{
			if(i<base-coverN)
			{
				dstPtr[i]=srcPtr[i];
			}
			
			else if(i<=base)
			{
				nZoom=1+splus*0.1;
				//注意数制转换存在的bug
				temp=srcPtr[i]*nZoom;
				if(temp-255>0.000001)
					dstPtr[i]=255;
				else
					dstPtr[i]=unsigned char(temp);
				
			}
			else//图象增强，将包络外的能量衰减，包络内的能量不变
			{
				nZoom=1-image*0.1;
				//注意数制转换存在的bug
				temp=srcPtr[i]*nZoom;
				if(temp-255>0.000001)
					dstPtr[i]=255;
				else
					dstPtr[i]=unsigned char(temp);
				
			}
		}
	}
	else
	{
		for(i=0;i<fft_point;i++)
		{
			dstPtr[i]=srcPtr[i];
		}
	}
	
}

/**************************************************************************************
                    趋势包络算法  by Jim Fang at 2007

	判断目标点dst和源点src的关系，控制包络的变化率
**************************************************************************************/
int WINAPI Jim_TrendCover(int src,int dst,double scale,int nTrend)
{
	int k;
	//scale=1.0;
	int varRise,varFall,nTrendLimit=5;
	varRise=int(50*scale);
	varFall=int(5.0*scale);//int(3*scale);

	if(nTrend>=0)
	{
		k=src+(nTrend+1)*varRise;
		//限制变化率
		if(dst>k)
		{
			dst=k;
		}
		else if(dst<src-varRise)
		{
			dst=src-varRise;
		}
	}
	else 
	{
		k = src + (nTrend-1) * varFall ;
		if ( dst < k ) 
		{
			dst = k ;
		} 
		else if ( dst > src + varFall )
		{
			dst = src + varFall ;
		}
	}
/*
	if ( src > dst ) 
	{
		if ( --nTrend < -1*nTrendLimit ) 
		{
			nTrend = -1*nTrendLimit ;
		}
	} 
	else if (src < dst) 
	{
		if ( ++nTrend > nTrendLimit )
		{
			nTrend = nTrendLimit ;
		}
	}
*/
	return dst;
}

/**************************************************************************************
                    FFT算法(Fourier Transform)  by Jim Fang at 2007
							    (Based on the priciple of Coolkey-Tukey)
	pr---双精度实型一维数组，长度为n
	     当l=0时，存放n个采样输入的实部，返回时存放离散傅立叶变换的模
		 当l=1时，存放傅立叶变换的n个实部，返回时存放逆傅立叶变换的模
	pi---双精度实型一维数组，长度为n
		 当l=0时，存放n个采样输入的虚部，返回时存放离散傅立叶变换的幅角
		 当l=1时，存放傅立叶变换的n个虚部，返回时存放逆傅立叶变换的幅角(幅角单位为度)
	n----整型变量，输入的点数
	k----整型变量，满足n=2^k
	fr---双精度实型一维数组，长度为n
		 当l=0时，返回傅立叶变换的实部
		 当l=1时，返回逆傅立叶变换的实部
	fi---双精度实型一维数组，长度为n
		 当l=0时，返回傅立叶变换的虚部
		 当l=1时，返回逆傅立叶变换的虚部
	l----整型变量
		 当l=0时，表示要求本函数计算傅立叶变换
		 当l=1时，表示要求本函数计算逆傅立叶变换
	il---整型变量
		 若il=0，表示不要求本函数计算傅立叶变换或逆变换的模与幅角
		 若il=1，表示要求本函数计算傅立叶变换或逆变换的模与幅角
**************************************************************************************/
void WINAPI Jim_FFT(double* pr,double* pi,int n,int k,int l,int il)
{ 

	int it,m,is,i,j,nv,l0;
	double p,q,s,vr,vi,poddr,poddi;
	double fr[256],fi[256];
	for (it=0; it<=n-1; it++)
	{ 
		m=it; is=0;
		for (i=0; i<=k-1; i++)
		{ 
			j=m/2; is=2*is+(m-2*j); m=j;
		}
		fr[it]=pr[is]; fi[it]=pi[is];
	}
	pr[0]=1.0; 
	pi[0]=0.0;
	
	p=(2*PI)/(1.0*n);
	pr[1]=cos(p); 
	pi[1]=-sin(p);
	
	//计算逆傅立叶变换
	if (l!=0) 
		pi[1]=-pi[1];
	for (i=2; i<=n-1; i++)
	{ 
		p=pr[i-1]*pr[1]; q=pi[i-1]*pi[1];
		s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
		pr[i]=p-q; pi[i]=s-p-q;
	}
	for (it=0; it<=n-2; it=it+2)
	{ 
		vr=fr[it]; vi=fi[it];
		fr[it]=vr+fr[it+1]; fi[it]=vi+fi[it+1];
		fr[it+1]=vr-fr[it+1]; fi[it+1]=vi-fi[it+1];
	}
	m=n/2; nv=2;
	for (l0=k-2; l0>=0; l0--)
	{ 
		m=m/2; nv=2*nv;
		for (it=0; it<=(m-1)*nv; it=it+nv)
			for (j=0; j<=(nv/2)-1; j++)
			{ 
				p=pr[m*j]*fr[it+j+nv/2];
				q=pi[m*j]*fi[it+j+nv/2];
				s=pr[m*j]+pi[m*j];
				s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
				poddr=p-q; poddi=s-p-q;
				fr[it+j+nv/2]=fr[it+j]-poddr;
				fi[it+j+nv/2]=fi[it+j]-poddi;
				fr[it+j]=fr[it+j]+poddr;
				fi[it+j]=fi[it+j]+poddi;
			}
	}
	//计算逆傅立叶变换
	if (l!=0)
	{
		for (i=0; i<=n-1; i++)
		{ 
			fr[i]=fr[i]/(1.0*n);
			fi[i]=fi[i]/(1.0*n);
		}
	}
	//计算傅立叶变换或逆变换的模pr[i]与幅角pi[i]
	if (il!=0)
	{
		for (i=0; i<=n-1; i++)
		{ 
			pr[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i]);
			if (fabs(fr[i])<0.000001*fabs(fi[i]))
			{ 
				if ((fi[i]*fr[i])>0) 
					pi[i]=90.0;
				else 
					pi[i]=-90.0;
			}
			else
				pi[i]=atan(fi[i]/fr[i])*360.0/(2*PI);
		}
	}
	return;
}

//计算400个点里的峰值、谷值、心率
void WINAPI Jim_CaculateParam(int *data,int *param,int point,int nParam)
{
	int		i,tmp=0;
	int		*result=(int *)param;
	int		peak=0,valley=128,middle=0,hr=0;
	int		offsetRise=0,offsetFall1=0,offsetFall2=0,offsetRiseLimit=20,offsetFallLimit=20; //控制下降幅度的检测
	bool    bFindPeak=false,bFindValley=false;
	int     pFirstPeak=0,pSecondPeak=0,pFirstValley=0;

	for(i=0;i<point;i++)
	{
		//寻找第一个峰值
		if(data[i]>peak)
			peak=data[i];
		if(data[i]<valley)
			valley=data[i];
	}
	middle=int((peak-valley)/3+valley);
	//找到起始搜索点
	for(i=0;i<point;i++)
	{
		if(data[i]==middle || data[i]-middle<3)
		{
			tmp=i;
			break;
		}
	}
	for(i=tmp;i<point;i++)
	{
	}
	result[0]=peak;
	result[1]=valley;
		//result[5]=pFirstValley;

}
/**************************************************************************************
             第一类整数阶贝塞尔函数算法(Bessel Function)  by Jim Fang at 2007

	n----整型变量，第一类贝塞尔函数的阶数，n>=0时
		 当n<0时，本函数按|n|计算
	x----双精度实型变量，自变量值

**************************************************************************************/
double WINAPI Jim_Bessel0(int n,double x)
{ 
	int i,m;
    double t,y,z,p,q,s,b0,b1;
    static double a[6]={ 57568490574.0,-13362590354.0,
             651619640.7,-11214424.18,77392.33017,
            -184.9052456};
    static double b[6]={ 57568490411.0,1029532985.0,
             9494680.718,59272.64853,267.8532712,1.0};
    static double c[6]={ 72362614232.0,-7895059235.0,
             242396853.1,-2972611.439,15704.4826,
             -30.16036606};
    static double d[6]={ 144725228443.0,2300535178.0,
             18583304.74,99447.43394,376.9991397,1.0};
    static double e[5]={ 1.0,-0.1098628627e-02,
             0.2734510407e-04,-0.2073370639e-05,
             0.2093887211e-06};
    static double f[5]={ -0.1562499995e-01,
             0.1430488765e-03,-0.6911147651e-05,
             0.7621095161e-06,-0.934935152e-07};
    static double g[5]={ 1.0,0.183105e-02,
             -0.3516396496e-04,0.2457520174e-05,
             -0.240337019e-06};
    static double h[5]={ 0.4687499995e-01,
             -0.2002690873e-03,0.8449199096e-05,
             -0.88228987e-06,0.105787412e-06};
    t=fabs(x);
    if (n<0) 
		n=-n;
    if (n!=1)
	{ 
		if (t<8.0)
		{ 
			y=t*t; p=a[5]; q=b[5];
			for (i=4; i>=0; i--)
			{ 
				p=p*y+a[i]; q=q*y+b[i];
			}
            p=p/q;
		}
        else
		{ 
			z=8.0/t; y=z*z;
            p=e[4]; q=f[4];
            for (i=3; i>=0; i--)
			{ 
				p=p*y+e[i]; q=q*y+f[i];
			}
            s=t-0.785398164;
            p=p*cos(s)-z*q*sin(s);
            p=p*sqrt(0.636619772/t);
		}
	}
    if (n==0)