telafft.c

用于Pattern Recongnition的一个图像处理程序－－Trace Transform。

/*
*This file is part of tela the Tensor Language.
*Copyright(c)1994-1995 Pekka Janhunen
*/

/*
  fftpack.C : A set of FFT routines in C.
  Algorithmically based on Fortran-77 FFTPACK by Paul N. Swarztrauber(Version 4, 1985).

  This file is also C++compatible.

  Pekka Janhunen 23.2.1995

 (reformatted by joerg arndt)
*/

/*isign is +1 for backward and -1 for forward transforms*/


#include <math.h>






/*----------------------------------------------------------------------
   passf2, passf3, passf4, passf5, passf. Complex FFT passes fwd and bwd.
  ----------------------------------------------------------------------*/

static void 
passf2(int ido, int l1, const double cc[], double ch[], const double wa1[], int isign)
     /*isign==+1 for backward transform*/
{
    int     i, k, ah, ac;
    double  ti2, tr2;

    if(ido<=2)
    {
	for(k=0; k<l1; k++)
	{
	    ah=k*ido;
	    ac=2*k*ido;
	    ch[ah]=cc[ac]+cc[ac+ido];
	    ch[ah+ido*l1]=cc[ac]-cc[ac+ido];
	    ch[ah+1]=cc[ac+1]+cc[ac+ido+1];
	    ch[ah+ido*l1+1]=cc[ac+1]-cc[ac+ido+1];
	}
    }
    else
    {
	for(k=0; k<l1; k++)
	{
	    for(i=0; i<ido-1; i+=2)
	    {
		ah=i+k*ido;
		ac=i+2*k*ido;
		ch[ah]=cc[ac]+cc[ac+ido];
		tr2=cc[ac]-cc[ac+ido];
		ch[ah+1]=cc[ac+1]+cc[ac+1+ido];
		ti2=cc[ac+1]-cc[ac+1+ido];
		ch[ah+l1*ido+1]=wa1[i]*ti2+isign*wa1[i+1]*tr2;
		ch[ah+l1*ido]=wa1[i]*tr2-isign*wa1[i+1]*ti2;
	    }
	}
    }
}				/*passf2*/


static void 
passf3(int ido, int l1, const double cc[], double ch[],
	const double wa1[], const double wa2[], int isign)
     /*isign==+1 for backward transform*/
{
    static const double taur=-0.5;
    static const double taui=0.866025403784439;
    int     i, k, ac, ah;
    double  ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;

    if(ido==2)
    {
	for(k=1; k<=l1; k++)
	{
	    ac=(3*k-2)*ido;
	    tr2=cc[ac]+cc[ac+ido];
	    cr2=cc[ac-ido]+taur*tr2;
	    ah=(k-1)*ido;
	    ch[ah]=cc[ac-ido]+tr2;

	    ti2=cc[ac+1]+cc[ac+ido+1];
	    ci2=cc[ac-ido+1]+taur*ti2;
	    ch[ah+1]=cc[ac-ido+1]+ti2;

	    cr3=isign*taui*(cc[ac]-cc[ac+ido]);
	    ci3=isign*taui*(cc[ac+1]-cc[ac+ido+1]);
	    ch[ah+l1*ido]=cr2-ci3;
	    ch[ah+2*l1*ido]=cr2+ci3;
	    ch[ah+l1*ido+1]=ci2+cr3;
	    ch[ah+2*l1*ido+1]=ci2-cr3;
	}
    }
    else
    {
	for(k=1; k<=l1; k++)
	{
	    for(i=0; i<ido-1; i+=2)
	    {
		ac=i+(3*k-2)*ido;
		tr2=cc[ac]+cc[ac+ido];
		cr2=cc[ac-ido]+taur*tr2;
		ah=i+(k-1)*ido;
		ch[ah]=cc[ac-ido]+tr2;
		ti2=cc[ac+1]+cc[ac+ido+1];
		ci2=cc[ac-ido+1]+taur*ti2;
		ch[ah+1]=cc[ac-ido+1]+ti2;
		cr3=isign*taui*(cc[ac]-cc[ac+ido]);
		ci3=isign*taui*(cc[ac+1]-cc[ac+ido+1]);
		dr2=cr2-ci3;
		dr3=cr2+ci3;
		di2=ci2+cr3;
		di3=ci2-cr3;
		ch[ah+l1*ido+1]=wa1[i]*di2+isign*wa1[i+1]*dr2;
		ch[ah+l1*ido]=wa1[i]*dr2-isign*wa1[i+1]*di2;
		ch[ah+2*l1*ido+1]=wa2[i]*di3+isign*wa2[i+1]*dr3;
		ch[ah+2*l1*ido]=wa2[i]*dr3-isign*wa2[i+1]*di3;
	    }
	}
    }
}				/*passf3*/


static void 
passf4(int ido, int l1, const double cc[], double ch[],
      const double wa1[], const double wa2[], const double wa3[], int isign)
     /*isign==-1 for forward transform and+1 for backward transform*/
{
    int     i, k, ac, ah;
    double  ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2,
            tr3, tr4;

    if(ido==2)
    {
	for(k=0; k<l1; k++)
	{
	    ac=4*k*ido+1;
	    ti1=cc[ac]-cc[ac+2*ido];
	    ti2=cc[ac]+cc[ac+2*ido];
	    tr4=cc[ac+3*ido]-cc[ac+ido];
	    ti3=cc[ac+ido]+cc[ac+3*ido];
	    tr1=cc[ac-1]-cc[ac+2*ido-1];
	    tr2=cc[ac-1]+cc[ac+2*ido-1];
	    ti4=cc[ac+ido-1]-cc[ac+3*ido-1];
	    tr3=cc[ac+ido-1]+cc[ac+3*ido-1];
	    ah=k*ido;
	    ch[ah]=tr2+tr3;
	    ch[ah+2*l1*ido]=tr2-tr3;
	    ch[ah+1]=ti2+ti3;
	    ch[ah+2*l1*ido+1]=ti2-ti3;
	    ch[ah+l1*ido]=tr1+isign*tr4;
	    ch[ah+3*l1*ido]=tr1-isign*tr4;
	    ch[ah+l1*ido+1]=ti1+isign*ti4;
	    ch[ah+3*l1*ido+1]=ti1-isign*ti4;
	}
    }
    else
    {
	for(k=0; k<l1; k++)
	{
	    for(i=0; i<ido-1; i+=2)
	    {
		ac=i+1+4*k*ido;
		ti1=cc[ac]-cc[ac+2*ido];
		ti2=cc[ac]+cc[ac+2*ido];
		ti3=cc[ac+ido]+cc[ac+3*ido];
		tr4=cc[ac+3*ido]-cc[ac+ido];
		tr1=cc[ac-1]-cc[ac+2*ido-1];
		tr2=cc[ac-1]+cc[ac+2*ido-1];
		ti4=cc[ac+ido-1]-cc[ac+3*ido-1];
		tr3=cc[ac+ido-1]+cc[ac+3*ido-1];
		ah=i+k*ido;
		ch[ah]=tr2+tr3;
		cr3=tr2-tr3;
		ch[ah+1]=ti2+ti3;
		ci3=ti2-ti3;
		cr2=tr1+isign*tr4;
		cr4=tr1-isign*tr4;
		ci2=ti1+isign*ti4;
		ci4=ti1-isign*ti4;
		ch[ah+l1*ido]=wa1[i]*cr2-isign*wa1[i+1]*ci2;
		ch[ah+l1*ido+1]=wa1[i]*ci2+isign*wa1[i+1]*cr2;
		ch[ah+2*l1*ido]=wa2[i]*cr3-isign*wa2[i+1]*ci3;
		ch[ah+2*l1*ido+1]=wa2[i]*ci3+isign*wa2[i+1]*cr3;
		ch[ah+3*l1*ido]=wa3[i]*cr4-isign*wa3[i+1]*ci4;
		ch[ah+3*l1*ido+1]=wa3[i]*ci4+isign*wa3[i+1]*cr4;
	    }
	}
    }
}				/*passf4*/


static void 
passf5(int ido, int l1, const double cc[], double ch[],
	const double wa1[], const double wa2[], const double wa3[], const double wa4[], int isign)
     /*isign==-1 for forward transform and+1 for backward transform*/
{
    static const double tr11=0.309016994374947;
    static const double ti11=0.951056516295154;
    static const double tr12=-0.809016994374947;
    static const double ti12=0.587785252292473;
    int     i, k, ac, ah;
    double  ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4,
            ti2, ti3, ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;

    if(ido==2)
    {
	for(k=1; k<=l1;++k)
	{
	    ac=(5*k-4)*ido+1;
	    ti5=cc[ac]-cc[ac+3*ido];
	    ti2=cc[ac]+cc[ac+3*ido];
	    ti4=cc[ac+ido]-cc[ac+2*ido];
	    ti3=cc[ac+ido]+cc[ac+2*ido];
	    tr5=cc[ac-1]-cc[ac+3*ido-1];
	    tr2=cc[ac-1]+cc[ac+3*ido-1];
	    tr4=cc[ac+ido-1]-cc[ac+2*ido-1];
	    tr3=cc[ac+ido-1]+cc[ac+2*ido-1];
	    ah=(k-1)*ido;
	    ch[ah]=cc[ac-ido-1]+tr2+tr3;
	    ch[ah+1]=cc[ac-ido]+ti2+ti3;
	    cr2=cc[ac-ido-1]+tr11*tr2+tr12*tr3;
	    ci2=cc[ac-ido]+tr11*ti2+tr12*ti3;
	    cr3=cc[ac-ido-1]+tr12*tr2+tr11*tr3;
	    ci3=cc[ac-ido]+tr12*ti2+tr11*ti3;
	    cr5=isign*(ti11*tr5+ti12*tr4);
	    ci5=isign*(ti11*ti5+ti12*ti4);
	    cr4=isign*(ti12*tr5-ti11*tr4);
	    ci4=isign*(ti12*ti5-ti11*ti4);
	    ch[ah+l1*ido]=cr2-ci5;
	    ch[ah+4*l1*ido]=cr2+ci5;
	    ch[ah+l1*ido+1]=ci2+cr5;
	    ch[ah+2*l1*ido+1]=ci3+cr4;
	    ch[ah+2*l1*ido]=cr3-ci4;
	    ch[ah+3*l1*ido]=cr3+ci4;
	    ch[ah+3*l1*ido+1]=ci3-cr4;
	    ch[ah+4*l1*ido+1]=ci2-cr5;
	}
    }
    else
    {
	for(k=1; k<=l1; k++)
	{
	    for(i=0; i<ido-1; i+=2)
	    {
		ac=i+1+(k*5-4)*ido;
		ti5=cc[ac]-cc[ac+3*ido];
		ti2=cc[ac]+cc[ac+3*ido];
		ti4=cc[ac+ido]-cc[ac+2*ido];
		ti3=cc[ac+ido]+cc[ac+2*ido];
		tr5=cc[ac-1]-cc[ac+3*ido-1];
		tr2=cc[ac-1]+cc[ac+3*ido-1];
		tr4=cc[ac+ido-1]-cc[ac+2*ido-1];
		tr3=cc[ac+ido-1]+cc[ac+2*ido-1];
		ah=i+(k-1)*ido;
		ch[ah]=cc[ac-ido-1]+tr2+tr3;
		ch[ah+1]=cc[ac-ido]+ti2+ti3;
		cr2=cc[ac-ido-1]+tr11*tr2+tr12*tr3;

		ci2=cc[ac-ido]+tr11*ti2+tr12*ti3;
		cr3=cc[ac-ido-1]+tr12*tr2+tr11*tr3;

		ci3=cc[ac-ido]+tr12*ti2+tr11*ti3;
		cr5=isign*(ti11*tr5+ti12*tr4);
		ci5=isign*(ti11*ti5+ti12*ti4);
		cr4=isign*(ti12*tr5-ti11*tr4);
		ci4=isign*(ti12*ti5-ti11*ti4);
		dr3=cr3-ci4;
		dr4=cr3+ci4;
		di3=ci3+cr4;
		di4=ci3-cr4;
		dr5=cr2+ci5;
		dr2=cr2-ci5;
		di5=ci2-cr5;
		di2=ci2+cr5;
		ch[ah+l1*ido]=wa1[i]*dr2-isign*wa1[i+1]*di2;
		ch[ah+l1*ido+1]=wa1[i]*di2+isign*wa1[i+1]*dr2;
		ch[ah+2*l1*ido]=wa2[i]*dr3-isign*wa2[i+1]*di3;
		ch[ah+2*l1*ido+1]=wa2[i]*di3+isign*wa2[i+1]*dr3;
		ch[ah+3*l1*ido]=wa3[i]*dr4-isign*wa3[i+1]*di4;
		ch[ah+3*l1*ido+1]=wa3[i]*di4+isign*wa3[i+1]*dr4;
		ch[ah+4*l1*ido]=wa4[i]*dr5-isign*wa4[i+1]*di5;
		ch[ah+4*l1*ido+1]=wa4[i]*di5+isign*wa4[i+1]*dr5;
	    }
	}
    }
}				/*passf5*/


static void 
passf(int*nac, int ido, int ip, int l1, int idl1,
     const double cc[], double c1[], double c2[], double ch[], double ch2[],
       const double wa[], int isign)
     /*isign is-1 for forward transform and+1 for backward transform*/
{
    int     idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, nt, idj, idl,
            inc, idp;
    double  wai, war;

    idot=ido / 2;
    nt=ip*idl1;
    ipph=(ip+1)/ 2;
    idp=ip*ido;
    if(ido>=l1)
    {
	for(j=1; j<ipph; j++)
	{
	    jc=ip-j;
	    for(k=0; k<l1; k++)
	    {
		for(i=0; i<ido; i++)
		{
		    ch[i+(k+j*l1)*ido]=
			cc[i+(j+k*ip)*ido]+cc[i+(jc+k*ip)*ido];
		    ch[i+(k+jc*l1)*ido]=
			cc[i+(j+k*ip)*ido]-cc[i+(jc+k*ip)*ido];
		}
	    }
	}
	for(k=0; k<l1; k++)
	    for(i=0; i<ido; i++)
		ch[i+k*ido]=cc[i+k*ip*ido];
    }
    else
    {
	for(j=1; j<ipph; j++)
	{
	    jc=ip-j;
	    for(i=0; i<ido; i++)
	    {
		for(k=0; k<l1; k++)
		{
		    ch[i+(k+j*l1)*ido]=cc[i+(j+k*ip)*ido]+cc[i+(jc+k*ip)*ido];
		    ch[i+(k+jc*l1)*ido]=cc[i+(j+k*ip)*ido]-cc[i+(jc+k*ip)*ido];
		}
	    }
	}
	for(i=0; i<ido; i++)
	    for(k=0; k<l1; k++)
		ch[i+k*ido]=cc[i+k*ip*ido];
    }

    idl=2-ido;
    inc=0;
    for(l=1; l<ipph; l++)
    {
	lc=ip-l;
	idl+=ido;
	for(ik=0; ik<idl1; ik++)
	{
	    c2[ik+l*idl1]=ch2[ik]+wa[idl-2]*ch2[ik+idl1];
	    c2[ik+lc*idl1]=isign*wa[idl-1]*ch2[ik+(ip-1)*idl1];
	}
	idlj=idl;
	inc+=ido;
	for(j=2; j<ipph; j++)
	{
	    jc=ip-j;
	    idlj+=inc;
	    if(idlj>idp)
		idlj-=idp;
	    war=wa[idlj-2];
	    wai=wa[idlj-1];
	    for(ik=0; ik<idl1; ik++)
	    {
		c2[ik+l*idl1]+=war*ch2[ik+j*idl1];
		c2[ik+lc*idl1]+=isign*wai*ch2[ik+jc*idl1];
	    }
	}
    }
    for(j=1; j<ipph; j++)
	for(ik=0; ik<idl1; ik++)
	    ch2[ik]+=ch2[ik+j*idl1];
    for(j=1; j<ipph; j++)
    {
	jc=ip-j;
	for(ik=1; ik<idl1; ik+=2)
	{
	    ch2[ik-1+j*idl1]=c2[ik-1+j*idl1]-c2[ik+jc*idl1];
	    ch2[ik-1+jc*idl1]=c2[ik-1+j*idl1]+c2[ik+jc*idl1];
	    ch2[ik+j*idl1]=c2[ik+j*idl1]+c2[ik-1+jc*idl1];
	    ch2[ik+jc*idl1]=c2[ik+j*idl1]-c2[ik-1+jc*idl1];
	}
    }
   *nac=1;
    if(ido==2)
	return;
   *nac=0;
    for(ik=0; ik<idl1; ik++)
	c2[ik]=ch2[ik];
    for(j=1; j<ip; j++)
    {
	for(k=0; k<l1; k++)
	{
	    c1[(k+j*l1)*ido+0]=ch[(k+j*l1)*ido+0];
	    c1[(k+j*l1)*ido+1]=ch[(k+j*l1)*ido+1];
	}
    }
    if(idot<=l1)
    {
	idij=0;
	for(j=1; j<ip; j++)
	{
	    idij+=2;
	    for(i=3; i<ido; i+=2)
	    {
		idij+=2;
		for(k=0; k<l1; k++)
		{
		    c1[i-1+(k+j*l1)*ido]=
			wa[idij-2]*ch[i-1+(k+j*l1)*ido]-
			isign*wa[idij-1]*ch[i+(k+j*l1)*ido];
		    c1[i+(k+j*l1)*ido]=
			wa[idij-2]*ch[i+(k+j*l1)*ido]+
			isign*wa[idij-1]*ch[i-1+(k+j*l1)*ido];
		}
	    }
	}
    }
    else
    {
	idj=2-ido;
	for(j=1; j<ip; j++)
	{
	    idj+=ido;
	    for(k=0; k<l1; k++)
	    {
		idij=idj;
		for(i=3; i<ido; i+=2)
		{
		    idij+=2;
		    c1[i-1+(k+j*l1)*ido]=
			wa[idij-2]*ch[i-1+(k+j*l1)*ido]-
			isign*wa[idij-1]*ch[i+(k+j*l1)*ido];
		    c1[i+(k+j*l1)*ido]=
			wa[idij-2]*ch[i+(k+j*l1)*ido]+
			isign*wa[idij-1]*ch[i-1+(k+j*l1)*ido];
		}
	    }
	}
    }
}				/*passf*/


/*----------------------------------------------------------------------
   radf2,radb2, radf3,radb3, radf4,radb4, radf5,radb5, radfg,radbg.
   Real FFT passes fwd and bwd.
  ----------------------------------------------------------------------*/

static void 
radf2(int ido, int l1, const double cc[], double ch[], const double wa1[])
{
    int     i, k, ic;
    double  ti2, tr2;

    for(k=0; k<l1; k++)
    {
	ch[2*k*ido]=
	    cc[k*ido]+cc[(k+l1)*ido];
	ch[(2*k+1)*ido+ido-1]=
	    cc[k*ido]-cc[(k+l1)*ido];
    }
    if(ido<2)
	return;
    if(ido !=2)
    {
	for(k=0; k<l1; k++)
	{
	    for(i=2; i<ido; i+=2)
	    {
		ic=ido-i;
		tr2=wa1[i-2]*cc[i-1+(k+l1)*ido]+wa1[i-1]*cc[i+(k+l1)*ido];
		ti2=wa1[i-2]*cc[i+(k+l1)*ido]-wa1[i-1]*cc[i-1+(k+l1)*ido];
		ch[i+2*k*ido]=cc[i+k*ido]+ti2;
		ch[ic+(2*k+1)*ido]=ti2-cc[i+k*ido];
		ch[i-1+2*k*ido]=cc[i-1+k*ido]+tr2;
		ch[ic-1+(2*k+1)*ido]=cc[i-1+k*ido]-tr2;
	    }
	}
	if(ido%2==1)
	    return;
    }
    for(k=0; k<l1; k++)
    {
	ch[(2*k+1)*ido]=-cc[ido-1+(k+l1)*ido];
	ch[ido-1+2*k*ido]=cc[ido-1+k*ido];
    }
}				/*radf2*/


static void 
radb2(int ido, int l1, const double cc[], double ch[], const double wa1[])
{
    int     i, k, ic;
    double  ti2, tr2;

    for(k=0; k<l1; k++)
    {
	ch[k*ido]=
	    cc[2*k*ido]+cc[ido-1+(2*k+1)*ido];
	ch[(k+l1)*ido]=
	    cc[2*k*ido]-cc[ido-1+(2*k+1)*ido];
    }
    if(ido<2)
	return;
    if(ido !=2)
    {
	for(k=0; k<l1;++k)
	{
	    for(i=2; i<ido; i+=2)
	    {
		ic=ido-i;
		ch[i-1+k*ido]=
		    cc[i-1+2*k*ido]+cc[ic-1+(2*k+1)*ido];
		tr2=cc[i-1+2*k*ido]-cc[ic-1+(2*k+1)*ido];
		ch[i+k*ido]=
		    cc[i+2*k*ido]-cc[ic+(2*k+1)*ido];
		ti2=cc[i+(2*k)*ido]+cc[ic+(2*k+1)*ido];
		ch[i-1+(k+l1)*ido]=
		    wa1[i-2]*tr2-wa1[i-1]*ti2;
		ch[i+(k+l1)*ido]=
		    wa1[i-2]*ti2+wa1[i-1]*tr2;
	    }
	}
	if(ido%2==1)
	    return;
    }
    for(k=0; k<l1; k++)
    {
	ch[ido-1+k*ido]=2*cc[ido-1+2*k*ido];
	ch[ido-1+(k+l1)*ido]=-2*cc[(2*k+1)*ido];
    }
}				/*radb2*/


static void 
radf3(int ido, int l1, const double cc[], double ch[],