wavedlg.cpp

自己设计的新型FFT算法

  struct compx v,t,w;
  nv2=N/2;
  f=N;
  for(m=1;(f=f/2)!=1;m++){;}
  nm1=N-1;

  /*位反转程序*/
  for (i=1;i<=nm1;i++)   
  { 
	  if(i<j){t=xin[j];xin[j]=xin[i];xin[i]=t;}
	  k=nv2;
	  while(k<j){j=j-k;k=k/2;}
	  j=j+k;
  }   

  /*蝶形运算*/
  
  for(l=1;l<=m;l++)        
  {le=pow(2,l);
   lei=le/2;
   v.real=1.0;
   v.imag=0.0;
   w.real=cos(PI/lei);
   w.imag=-sin(PI/lei);

   for(j=1;j<=lei;j++)    
   {for(i=j;i<=N;i=i+le)
   {ip=i+lei;
    t=EE(xin[ip],v);
	       
    xin[ip].real=xin[i].real-t.real;
	xin[ip].imag=xin[i].imag-t.imag;
    xin[i].real=xin[i].real+t.real;
	xin[i].imag=xin[i].imag+t.imag;
   }
   v=EE(v,w);
  }
  }
}

void CWaveDlg::rfft(int m0, double *x)
{
    int nn, nn2, i, j;
    double d, ti0, tr0, ti1, tr1, ac, as;
    double sstep, cstep, s, c, ww;

    nn = 1 << m0;
    nn2 = nn/2;

    cfftall(m0-1, x, 1.0);
    d   = PI * 2.0 / (double)nn;
    cstep = cos(d);
    sstep = sin(d);

    c = cstep;
    s = sstep;

    for (i = 2; i < nn2; i += 2) {
    	j = nn - i;
        tr0 = (x[i]   + x[j])   * 0.5;
        ti0 = (x[i+1] - x[j+1]) * 0.5;
        tr1 = (x[i+1] + x[j+1]) * 0.5;
        ti1 = (x[j]   - x[i])   * 0.5;

        ac = tr1 * c - ti1 * s;
        as = ti1 * c + tr1 * s;

        x[j]   =  tr0 - ac;
        x[j+1] = -ti0 + as;
        x[i]   =  tr0 + ac;
        x[i+1] =  ti0 + as;

        ww = c * cstep - s * sstep;
        s  = s * cstep + c * sstep;
        c  = ww;
    }
    tr0     = x[0];
    tr1     = x[1];
    x[0]    = tr0 + tr1;
    x[1]   = tr0 - tr1;
    //x[1]    = 0.0;
    //x[nn]   = tr0 - tr1;
    //x[nn+1] = 0.0;
}
void CWaveDlg::irfft(int m0, double *x)
{
    int nn, i, j;
    double d, ti0, tr0, ti1, tr1, ac, as;
    double sstep, cstep, s, c, ww;

    nn = 1 << m0;

    d   = PI * 2.0 / (double)nn;
    cstep = cos(d);
    sstep = sin(d);

    c = cstep;
    s = sstep;

    for (i = 2; i < (j = nn - i); i += 2) {
        tr0 = (x[i]   + x[j])* 0.5;
        ti0 = (x[i+1] - x[j+1]) * 0.5;
        as  = (x[i+1] + x[j+1]) * 0.5;
        ac  = (x[i]   - x[j]) * 0.5;

        tr1 = ac * c + as * s;
        ti1 = as * c - ac * s;

        x[i]   =  tr0 - ti1;
        x[i+1] =  ti0 + tr1;
        x[j]   =  tr0 + ti1;
        x[j+1] =  tr1 - ti0;

        ww = c * cstep - s * sstep;
        s  = s * cstep + c * sstep;
        c  = ww;
    }
    //tr0  = x[0] + x[nn];
    //tr1  = x[0] - x[nn];
	tr0  = x[0] + x[1];
    tr1  = x[0] - x[1];
    x[0] = tr0 * 0.5;
    x[1] = tr1 * 0.5;
    cfftall(m0-1, x, -1.0);
}


void CWaveDlg::cfftall(int m0, double *x, double ainv)
{
    int    i, j, lm, li, k, lmx, lmx2, np, lix;
    double  temp1, temp2;
    double  c, s, csave, sstep, cstep;
    double  c0, s0, c1, s1;

    lmx = 1 << m0;

    csave = PI * 2.0 / (double)lmx;
    cstep = cos(csave);
    sstep = sin(csave);

    lmx += lmx;
    np   = lmx;

/*----- fft butterfly numeration */
    for (i = 3; i <= m0; ++i) {
	lix = lmx;
	lmx >>= 1;
	lmx2 = lmx >> 1;
	c = cstep;
	s = sstep;
	s0 = ainv * s;
	c1 = -s;
	s1 = ainv * c;
	for (li = 0; li < np; li += lix ) {
	    j = li;
	    k = j + lmx;
	    temp1  = x[j] - x[k];
	    x[j]  += x[k];
	    x[k]   = temp1;
	    temp2  = x[++j] - x[++k];
	    x[j]  += x[k];
	    x[k]   = temp2;

	    temp1  = x[++j] - x[++k];
	    x[j]  += x[k];
	    temp2  = x[++j] - x[++k];
	    x[j]  += x[k];
	    x[k-1] = c * temp1 - s0 * temp2;
	    x[k]   = s0 * temp1 + c * temp2;

	    j = li + lmx2;
	    k = j + lmx;
	    temp1  = x[j] - x[k];
	    x[j]  += x[k];
	    temp2  = x[++j] - x[++k];
	    x[j]  += x[k];
	    x[k-1] = -ainv * temp2;
	    x[k]   =  ainv * temp1;

	    temp1  = x[++j] - x[++k];
	    x[j]  += x[k];
	    temp2  = x[++j] - x[++k];
	    x[j]  += x[k];
	    x[k-1] = c1 * temp1 - s1 * temp2;
	    x[k]   = s1 * temp1 + c1 * temp2;

	}
	for (lm = 4; lm < lmx2; lm += 2) {
	    csave = c;
	    c = cstep * c - sstep * s;
            s = sstep * csave + cstep * s;

	    s0 = ainv * s;
	    c1 = -s;
	    s1 = ainv * c;

	    for (li = 0; li < np; li += lix ) {
		j = li + lm;
		k = j + lmx;
		temp1  = x[j] - x[k];
	        x[j]  += x[k];
		temp2  = x[++j] - x[++k];
		x[j]  += x[k];
	        x[k-1] = c * temp1 - s0 * temp2;
	        x[k]   = s0 * temp1 + c * temp2;

		j = li + lm + lmx2;
		k = j + lmx;
		temp1  = x[j] - x[k];
	        x[j]  += x[k];
		temp2  = x[++j] - x[++k];
		x[j]  += x[k];
	        x[k-1] = c1 * temp1 - s1 * temp2;
	        x[k]   = s1 * temp1 + c1 * temp2;
	    }
	}
	csave = cstep;
	cstep = 2.0 * cstep * cstep - 1.0;
	sstep = 2.0 * sstep * csave;
    }
    if (m0 >= 2)
	for (li = 0; li < np; li += 8) {
	    j = li;
	    k = j + 4;
	    temp1 = x[j] - x[k];
	    x[j] += x[k];
	    temp2 = x[++j] - x[++k];
	    x[j] += x[k];
	    x[k-1] = temp1;
	    x[k]   = temp2;
	    temp1  = x[++j] - x[++k];
	    x[j]  += x[k];
	    temp2  = x[++j] - x[++k];
	    x[j]  += x[k];
	    x[k-1] = -ainv * temp2;
	    x[k]   =  ainv * temp1;
	}
    for (li = 0; li < np; li += 4) {
	j = li;
	k = j + 2;
	temp1 = x[j] - x[k];
	x[j]  += x[k];
	x[k]   = temp1;
	temp2  = x[++j] - x[++k];
	x[j]  += x[k];
	x[k]   = temp2;
    }

/*----- fft bit reversal */
    lmx = np / 2;
    j = 0;
    for (i = 2; i < np - 2; i += 2) {
	k = lmx;
	while(k <= j) {
	    j -= k;
	    k >>= 1;
	}
	j += k;
        if ( i < j ) {
	    temp1 = x[j];
	    x[j] = x[i];
	    x[i] = temp1;
	    lm = j + 1;
	    li = i + 1;
	    temp1 = x[lm];
	    x[lm] = x[li];
	    x[li] = temp1;
	}
    }
    if (ainv == 1.0) return;

    temp1 = 1.0 / (double)lmx;
    for (i = 0; i < np; ++i) *x++ *= temp1;
    return;
}

/*
int CWaveDlg::rfft(int dpts, double *yin, double *yout, int dir)
{
	 int i;
   double fact;

   if( FT_BACK == dir ) {
      fact = 2. / ((double) dpts) ;
      for( i = 0 ; i < dpts ; i++ )
         *(yout+i) = fact * *(yin+i);
      realft(yout-1, dpts/2, -1);  /* yout: unit offset in Num. Rec.*/
/*	 for(i=3;i<dpts;i+=2)
        *(yout+i)=-*(yout+i);
   } else {
      for( i = 0 ; i < dpts ; i++ )
         *(yout+i) = *(yin+i);
      realft(yout-1, dpts/2, 1);
	  for(i=3;i<dpts;i+=2)
	 	  *(yout+i)=-*(yout+i);
   }

   return 1;

}*/


 

int CWaveDlg::rfft_mply(int dpts, double *x, double *y, double *z)
{
	int i;
   double real, imag;

   if( dpts >= 2 ) {                                /* these are not complex */
      *z     = *x     * *y     ;
      *(z+1) = *(x+1) * *(y+1) ;
   }

   for( i = 2 ; i < dpts ; i += 2 ) {                           /* other are */
      real     = *(x+i) * *(y+i)   - *(x+i+1) * *(y+i+1) ;
      imag     = *(x+i) * *(y+i+1) + *(x+i+1) * *(y+i)   ;
      *(z+i)   = real ;
      *(z+i+1) = imag ;
   }

   return 1;

}

int CWaveDlg::rfft_dvde(int dpts, double *x, double *y, double *z)
{
	int i;
   double real, imag, denom;

   if( dpts >= 2 ) {                                /* these are not complex */
      if( ( 0. == *y ) || ( 0. == *(y+1) ) )
         return EXIT_FAILURE;
      *z     = *x     / *y     ;
      *(z+1) = *(x+1) / *(y+1) ;
   }

   for( i = 2 ; i < dpts ; i += 2 ) {                          /* others are */
      denom    = *(y+i) * *(y+i) + *(y+i+1) * *(y+i+1) ;
      if( 0.0 == denom )
         return EXIT_FAILURE;                         /* avoid division by 0 */
      real     = *(x+i) * *(y+i) + *(x+i+1) * *(y+i+1) ;
      imag     = *(x+i+1) * *(y+i) - *(x+i) * *(y+i+1) ;
      *(z+i)   = real / denom;
      *(z+i+1) = imag / denom;
	   }

   return 1;

}

int CWaveDlg::rfft_corr(int dpts, double *x, double *y, double *z)
{
	int i;
   double real, imag;

   if( dpts >= 2 ) {                                /* these are not complex */
      *z     = *x     * *y     ;
      *(z+1) = *(x+1) * *(y+1) ;
   }

   for( i = 2 ; i < dpts ; i += 2 ) {                           /* other are */
      real     = *(x+i) * *(y+i) + *(x+i+1) * *(y+i+1) ;
      imag     = - *(x+i) * *(y+i+1) + *(x+i+1) * *(y+i)   ;
      *(z+i)   = real ;
      *(z+i+1) = imag ;
   }

   return 1;

}

void CWaveDlg::realft(double *data, int n, int isign)
{
	 int i,i1,i2,i3,i4,n2p3;
    double c1=0.5,c2,h1r,h1i,h2r,h2i;
    double wr,wi,wpr,wpi,wtemp,theta;

    theta=3.141592653589793/(double) n;
    if (isign == 1) {
        c2 = -0.5;
        four1(data,n,1);
    } else {
        c2=0.5;
        theta = -theta;
    }
    wtemp=sin(0.5*theta);
    wpr = -2.0*wtemp*wtemp;
    wpi=sin(theta);
    wr=1.0+wpr;
    wi=wpi;
    n2p3=2*n+3;
    for (i=2;i<=n/2;i++) {
        i4=1+(i3=n2p3-(i2=1+(i1=i+i-1)));
        h1r=c1*(data[i1]+data[i3]);
        h1i=c1*(data[i2]-data[i4]);
        h2r = -c2*(data[i2]+data[i4]);
        h2i=c2*(data[i1]-data[i3]);
        data[i1]=h1r+wr*h2r-wi*h2i;
        data[i2]=h1i+wr*h2i+wi*h2r;
        data[i3]=h1r-wr*h2r+wi*h2i;
        data[i4] = -h1i+wr*h2i+wi*h2r;
        wr=(wtemp=wr)*wpr-wi*wpi+wr;
        wi=wi*wpr+wtemp*wpi+wi;
    }
    if (isign == 1) {
        data[1] = (h1r=data[1])+data[2];
        data[2] = h1r-data[2];
    } else {
        data[1]=c1*((h1r=data[1])+data[2]);
        data[2]=c1*(h1r-data[2]);
        four1(data,n,-1);
    }

}

void CWaveDlg::four1(double *data, int nn, int isign)
{
	int n,mmax,m,j,istep,i;
    double wtemp,wr,wpr,wpi,wi,theta;
    double tempr,tempi;

    n=nn << 1;
    j=1;
    for (i=1;i<n;i+=2) {
        if (j > i) {
            SWAP(data[j],data[i]);
            SWAP(data[j+1],data[i+1]);
        }
        m=n >> 1;
        while (m >= 2 && j > m) {
            j -= m;
            m >>= 1;
        }
        j += m;
    }
    mmax=2;
    while (n > mmax) {
        istep=2*mmax;
        theta=6.28318530717959/(isign*mmax);
        wtemp=sin(0.5*theta);
        wpr = -2.0*wtemp*wtemp;
        wpi=sin(theta);
        wr=1.0;
        wi=0.0;
        for (m=1;m<mmax;m+=2) {
            for (i=m;i<=n;i+=istep) {
                j=i+mmax;
                tempr=wr*data[j]-wi*data[j+1];
                tempi=wr*data[j+1]+wi*data[j];
                data[j]=data[i]-tempr;
                data[j+1]=data[i+1]-tempi;
                data[i] += tempr;
                data[i+1] += tempi;
            }
            wr=(wtemp=wr)*wpr-wi*wpi+wr;
            wi=wi*wpr+wtemp*wpi+wi;
        }
        mmax=istep;
    }

}


double CWaveDlg::maxNum(double x, double y)
{
  double mx;
  if (x>y)
      mx=x;
  else mx=y;
  return (mx);
}