fft.cpp

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

		return(p);
    b0=p;
    if (t<8.0)
	{ 
		y=t*t; p=c[5]; q=d[5];
        for (i=4; i>=0; i--)
		{ 
			p=p*y+c[i]; q=q*y+d[i]; 
		}
        p=x*p/q;
	}
    else
	{ 
		z=8.0/t; y=z*z;
        p=g[4]; q=h[4];
        for (i=3; i>=0; i--)
		{ 
			p=p*y+g[i]; q=q*y+h[i];
		}
        s=t-2.356194491;
        p=p*cos(s)-z*q*sin(s);
        p=p*x*sqrt(0.636619772/t)/t;
	}
    if (n==1) 
		return(p);
    b1=p;
    if (x==0.0) 
		return(0.0);
    s=2.0/t;
    if (t>1.0*n)
	{ 
		if (x<0.0) b1=-b1;
        for (i=1; i<=n-1; i++)
		{ 
			p=s*i*b1-b0; b0=b1; b1=p;
		}
	}
    else
	{ 
		m=(n+(int)sqrt(40.0*n))/2;
        m=2*m;
        p=0.0; q=0.0; b0=1.0; b1=0.0;
        for (i=m-1; i>=0; i--)
		{ 
			t=s*(i+1)*b0-b1;
            b1=b0; b0=t;
            if (fabs(b0)>1.0e+10)
			{ 
				b0=b0*1.0e-10; b1=b1*1.0e-10;
                p=p*1.0e-10; q=q*1.0e-10;
			}
            if ((i+2)%2==0) 
				q=q+b0;
            if ((i+1)==n) 
				p=b1;
		}
        q=2.0*q-b0; p=p/q;
	}
    if ((x<0.0)&&(n%2==1)) 
		p=-p;
    return(p);
}


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

	n----整型变量，第一类贝塞尔函数的阶数，n>=0时
		 当n<0时，本函数按|n|计算
	x----双精度实型变量，自变量值
**************************************************************************************/
double WINAPI Jim_Bessel0_R(int n,double x)
{ 
	int i,m;
    double t,y,p,b0,b1,q;
    static double a[7]={ 1.0,3.5156229,3.0899424,1.2067492,
                         0.2659732,0.0360768,0.0045813};
    static double b[7]={ 0.5,0.87890594,0.51498869,
              0.15084934,0.02658773,0.00301532,0.00032411};
    static double c[9]={ 0.39894228,0.01328592,0.00225319,
                        -0.00157565,0.00916281,-0.02057706,
                         0.02635537,-0.01647633,0.00392377};
    static double d[9]={ 0.39894228,-0.03988024,-0.00362018,
                        0.00163801,-0.01031555,0.02282967,
                        -0.02895312,0.01787654,-0.00420059};
    if (n<0) 
		n=-n;
    t=fabs(x);
    if (n!=1)
	{ 
		if (t<3.75)
		{ 
			y=(x/3.75)*(x/3.75); p=a[6];
            for (i=5; i>=0; i--)
              p=p*y+a[i];
		}
        else
		{ 
			y=3.75/t; p=c[8];
            for (i=7; i>=0; i--)
              p=p*y+c[i];
            p=p*exp(t)/sqrt(t);
		}
	}
    if (n==0)
		return(p);
    q=p;
    if (t<3.75)
	{ 
		y=(x/3.75)*(x/3.75); p=b[6];
        for (i=5; i>=0; i--) 
			p=p*y+b[i];
        p=p*t;
	}
    else
	{ 
		y=3.75/t; p=d[8];
        for (i=7; i>=0; i--) 
			p=p*y+d[i];
        p=p*exp(t)/sqrt(t);
	}
    if (x<0.0) 
		p=-p;
    if (n==1)
		return(p);
    if (x==0.0) 
		return(0.0);
    y=2.0/t; t=0.0; b1=1.0; b0=0.0;
    m=n+(int)sqrt(40.0*n);
    m=2*m;
    for (i=m; i>0; i--)
	{
		p=b0+i*y*b1; b0=b1; b1=p;
        if (fabs(b1)>1.0e+10)
		{ 
			t=t*1.0e-10; b0=b0*1.0e-10;
            b1=b1*1.0e-10;
		}
        if (i==n) 
			t=b0;
	}
    p=t*q/b1;
    if ((x<0.0)&&(n%2==1))
		p=-p;
    return(p);
}

/**************************************************************************************
             五点三次平滑  by Jim Fang at 2007

	n----整型变量，给定等距的采样点个数n>5
	temp-双精度实型一维数组，长度为n，n个等距的采样电数据
	yy---双精度实型一维数组，长度为n，返回n个等距采样数据的平滑结果，结果不允许出现负数
**************************************************************************************/
void WINAPI Jim_Smooth(int n,unsigned char *temp,double *yy)
{ 
	int i;
	
    if (n<5)
    { 
		for (i=0; i<=n-1; i++) 
			yy[i]=temp[i];
	}
    else
    { 
		yy[0]=69.0*temp[0]+4.0*temp[1]-6.0*temp[2]+4.0*temp[3]-temp[4];
        yy[0]=yy[0]/70.0;
        yy[1]=2.0*temp[0]+27.0*temp[1]+12.0*temp[2]-8.0*temp[3];
        yy[1]=(yy[1]+2.0*temp[4])/35.0;
        for (i=2; i<=n-3; i++)
        { 
			yy[i]=-3.0*temp[i-2]+12.0*temp[i-1]+17.0*temp[i];
            yy[i]=(yy[i]+12.0*temp[i+1]-3.0*temp[i+2])/35.0;
			if(yy[i]-0.0<-0.000001)
				yy[i]=0.0;
        }
        yy[n-2]=2.0*temp[n-5]-8.0*temp[n-4]+12.0*temp[n-3];
        yy[n-2]=(yy[n-2]+27.0*temp[n-2]+2.0*temp[n-1])/35.0;
        yy[n-1]=-temp[n-5]+4.0*temp[n-4]-6.0*temp[n-3];
        yy[n-1]=(yy[n-1]+4.0*temp[n-2]+69.0*temp[n-1])/70.0;
	}
}

/**************************************************************************************
             最小二乘曲线拟合  by Jim Fang at 2007

    x----双精度实型一维数组，长度为n，存放给定n个数据点的X坐标
	y----双精度实型一维数组，长度为n，存放给定n个数据点的Y坐标
	n----整型变量，给定数据点的个数
	a----双精度实型一维数组，长度为m，返回m-1次拟合多项式的m个系数
	m----整型变量，拟合多项数的项数，即拟合多项式的最高次数为m-1。这里要求m<20
	dt---双精度实型一维数组，长度为3，其中
	     dt[0]返回拟合多项式与数据点误差的平方和
		 dt[1]返回拟合多项式与数据点误差的绝对值之和
		 dt[2]返回拟合多项式与数据点误差绝对值的最大值
**************************************************************************************/
void WINAPI Jim_Least(double *x,double *y,int n,double *a,int m,double *dt)
{ 
	int i,j,k;
    double z,p,c,g,q,d1,d2,s[20],t[20],b[20];
    for (i=0; i<=m-1; i++) 
		a[i]=0.0;
    if (m>n) 
		m=n;
    if (m>20) 
		m=20;
    z=0.0;
    for (i=0; i<=n-1; i++) 
		z=z+x[i]/(1.0*n);
    b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;
    for (i=0; i<=n-1; i++)
    { 
		p=p+(x[i]-z); 
		c=c+y[i];
	}
    c=c/d1; p=p/d1;
    a[0]=c*b[0];
    if (m>1)
    { 
		t[1]=1.0; t[0]=-p;
        d2=0.0; c=0.0; g=0.0;
        for (i=0; i<=n-1; i++)
        {
			q=x[i]-z-p; d2=d2+q*q;
            c=c+y[i]*q;
            g=g+(x[i]-z)*q*q;
		}
        c=c/d2; p=g/d2; q=d2/d1;
        d1=d2;
        a[1]=c*t[1]; a[0]=c*t[0]+a[0];
	}
    for (j=2; j<=m-1; j++)
	{ 
		s[j]=t[j-1];
        s[j-1]=-p*t[j-1]+t[j-2];
        if (j>=3)
          for (k=j-2; k>=1; k--)
            s[k]=-p*t[k]+t[k-1]-q*b[k];
        s[0]=-p*t[0]-q*b[0];
        d2=0.0; c=0.0; g=0.0;
        for (i=0; i<=n-1; i++)
		{ 
			q=s[j];
            for (k=j-1; k>=0; k--)
              q=q*(x[i]-z)+s[k];
            d2=d2+q*q; c=c+y[i]*q;
            g=g+(x[i]-z)*q*q;
		}
        c=c/d2; p=g/d2; q=d2/d1;
        d1=d2;
        a[j]=c*s[j]; t[j]=s[j];
        for (k=j-1; k>=0; k--)
		{ 
			a[k]=c*s[k]+a[k];
            b[k]=t[k]; t[k]=s[k];
		}
	}
    dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
    for (i=0; i<=n-1; i++)
	{ 
		q=a[m-1];
        for (k=m-2; k>=0; k--)
          q=a[k]+q*(x[i]-z);
        p=q-y[i];
        if (fabs(p)>dt[2]) 
			dt[2]=fabs(p);
        dt[0]=dt[0]+p*p;
        dt[1]=dt[1]+fabs(p);
	}
}

/**************************************************************************************
             图像膨胀算法  by Jim Fang at 2007

	为保证图256点频谱在屏幕上的显示效果，采用膨胀差值算法计算出过渡点
	暂用以下简单算法代替
	if(yPlusDis==0)
		rgb=imageBuffer[i];
	else
		rgb=imageBuffer[i-1]+((imageBuffer[i]-imageBuffer[i-1])/yPlusDis)*j;
**************************************************************************************/
void WINAPI Jim_Swell(int src,int dst,int *data,int no)
{
}

/**************************************************************************************
             窗函数法设计数字滤波器FIR  by Jim Fang at 2007

    l-------FIR滤波器的数据序列长度
	iband---滤波器类型
			1:低通，2:高通，3:带通，4:带阻
	fl------低频截止频率Hz
	fh------高频截止频率Hz(对于带通和带阻而言)
	fs------采样频率Hz
            0.<=fln,fhn<=.5
          |---        |   ---       |   ---           |--      --
          |   |       |  |          |  |   |          |  |    |
          |   |       |  |          |  |   |          |  |    |
        --|------    -|--------    -|-----------    --|--------------
          0  fln      0  fln        0 fln  fhn        0  fln  fhn
             fhn=*       fhn=*
        iband=1 LP    iband=2 HP      iband=3 BP        iband=4 BS

    iwindow-窗函数类型
            1: rectangular window ,   2: triangular window ,
            3: cosin window ,         4: Hanning window ,
            5: Hamming window ,       6: Blackman window ,
            7: Papoulis window .
	b-------滤波器系数数组h(z)=b(0)+b(1)z^(-1)+ ... +b(l-1)z^(-l+1)
	w------- l dimensioned real work array.
    ierror--函数运算结果
	        0:运行正确
            1:无效的数据长度(l<=0)
            2:无效的窗函数类型
            3:无效的滤波器类型
            4:无效的截止频率
    根据窗函数系数w[]，计算滤波器系数b[]
**************************************************************************************/
int WINAPI Jim_FirFilter(int l,int iband,double fl,double fh,double fs,int iwindow,float *b)
{							  
	//float *b=(float *)coe;
	int ierror;
	//double w[256];
	double pi,fln,fhn,wc1,wc2,s;
	int lim,i;
	double dly;
	pi=4.0*atan(1.0);
	fln=fl/fs;
	fhn=fh/fs;
	for(i=0;i<l;i++) 
		b[i]=0.0;
	ierror=0;
	dly=l/2.0;
	lim=l/2;
	//判断输入参数的正确性
	//if(dly==(float)lim) 
	//	ierror=1;
	if(iwindow<1 || iwindow>7) 
		ierror=2;
	if(iband<1 || iband>4) 
		ierror=3;
	if(fln<0.001 || fln>0.5) 
		ierror=4;
	if(iband>=3 && fln>=fhn) 
		ierror=4;
	//if(iband>=3 && fhn>=0.5) 
	//	ierror=4;
	if(ierror!=0) 
		return ierror;
	dly-=0.5;
	lim-=1;
	wc1=2.0*pi*fln;
	//带通/带阻类型的处理
	if(iband>=3) 
		wc2=2.0*pi*fhn;
	//Jim_Window(w,l,iwindow,ierror);
	switch (iband)
    {
        //低通
        case 1:
		{
			for(i=0;i<=lim;i++)
			{
				s=i-dly;
				b[i]=float((sin(wc1*s))/(pi*s));
				b[i]=float(b[i]*1.0);//w[i];
				b[l-1-i]=b[i];
			}
			b[(l-1)/2]=float(wc1/pi);
			return ierror;
		}
		//高通
		case 2:
		{
			
			for(i=0;i<=lim;i++)
			{
				s=i-dly;//dly数据的一半长度
				b[i]=float((sin(pi*s)-sin(wc1*s))/(pi*s));
				b[i]=float(b[i]*1.0);//w[i];
				//对称
				b[l-1-i]=b[i];
			}
			b[(l-1)/2]=float(1.0-wc1/pi);
			fprintf(fp,"*********************************************************************\n");
			return ierror;
		}
		//带通
		case 3:
		{
			
			for(i=0;i<=lim;i++)
			{
				s=i-dly;
				b[i]=float((sin(wc2*s)-sin(wc1*s))/(pi*s));
				b[i]=float(b[i]*1.0);//w[i];
				b[l-1-i]=b[i];
			}
			b[(l-1)/2]=float((wc2-wc1)/pi);
			return ierror;
		}
		//带阻
		case 4:
		{
			for(i=0;i<=lim;i++)
			{
				s=i-dly;
				b[i]=float((sin(wc1*s)+sin(pi*s)-sin(wc2*s))/(pi*s));
				b[i]=float(b[i]*1.0);//w[i];
				b[l-1-i]=b[i];
			}
			b[(l-1)/2]=float((wc1+pi-wc2)/pi);
			return ierror;
		}
		default:break;
     }
	return ierror;
}

/**************************************************************************************
             窗函数  by Jim Fang at 2007

    w-------n dimension real array.the result is in w(0) to w(n-1)
	n-------数据序列长度
    iwindow-窗函数类型
            1: rectangular window ,   2: triangular window ,
            3: cosin window ,         4: Hanning window ,
            5: Hamming window ,       6: Blackman window ,
            7: Papoulis window .
    ierror--出错类型
	        0: no error,    1: Iwindow out of range.
**************************************************************************************/
/*void WINAPI Jim_Window(double *w,int n,int iwindow,int ierror)
{
	double pi,pn,a,b,c;
	int i;
	ierror=1;
	if(iwindow<1) 
		return;
	ierror=0;
	pi=4.*atan(1.);
	pn=2.*pi/(double)(n);
    switch (iwindow)
	{
		case 1:
			{
				for(i=0;i<n;i++)
				w[i]=1.0;
				break;
			}
		case 2:
			{
				for(i=0;i<n;i++)
				w[i]=1.-fabs(1.-2.*i/(double)(n));
				break;
			}
		case 3:
			{
				for(i=0;i<n;i++)
				w[i]=sin(pn*i/2.);
				break;
			}
		case 4:
			{
				for(i=0;i<n;i++)
				w[i]=0.5*(1.0-cos(pn*i));
				break;
			}
		case 5:
			{
				for(i=0;i<n;i++)
				w[i]=0.54-0.46*cos(pn*i);
				break;
			}
		case 6:
			{
				for(i=0;i<n;i++)
				w[i]=0.42-0.5*cos(pn*i)+0.08*cos(2.*pn*i);
				break;
			}
		case 7:
			{
				for(i=0;i<n;i++)
				{
						a=fabs(sin(pn*i))/pi;
					b=1.-2.*(fabs(i-n/2.))/(double)(n);
					c=cos(pn*i);
					w[i]=a-b*c;
				}
    		}
	}
}*/