convlv.cpp

Visual C++ 常用数值算法集 源代码

int cint(double x)
{
	int temp;
	double iprt;
	if (x>0)
	{
	x=modf(x,&iprt);
		if(fabs(x)<0.5)
			temp=int(iprt);
		else
			temp=int(iprt+1);
	}
	else if(x==0)
		temp=0;
	else
	{
		x=modf(x,&iprt);
		if(fabs(x)<0.5)
			temp=int(iprt);
		else
			temp=int(iprt)-1;
	}
		return temp;
}

void convlv(double data[],int n,double respns[],int m,
			int isign,double ans[])
{
	int i,no2,temp;
    double fft[33],ans1,dum1,dum2,dum;
	temp=-1;
    for (i = 1; i<=cint(m - 1) / 2; i++)
	{
        respns[n + 1 - i] = respns[m + 1 - i];
    }
    for (i = cint(m + 3) / 2; i<=n - cint(m - 1) / 2; i++)
	{
        respns[i] = 0.0;
    }
    twofft(data, respns, fft, ans, n);
    no2 = cint(n / 2);
    for (i = 1; i<=no2 + 1; i++)
	{
        if( isign == 1)
		{
          dum = ans[2 * i - 1];
          dum1 = fft[2 * i - 1] * dum - fft[2 * i] * ans[2 * i];
          ans[2 * i - 1] = dum1 / no2;
          dum2 = fft[2 * i - 1] * ans[2 * i] + fft[2 * i] * dum;
          ans[2 * i] = dum2 / no2;
		}
        else if (isign == -1)
		{
          if( dum == 0.0 || ans[2 * i] == 0) 
		  {
            cout<<"deconvolving at a response zero"<<endl;
            _c_exit();
          }
          ans1 = fft[2 * i - 1] * dum + fft[2 * i] * ans[2 * i];
          dum1 = dum * dum + ans[2 * i] * ans[2 * i];
          ans[2 * i - 1] = ans1 / dum1 / no2;
          ans1 = fft[2 * i] * dum - fft[2 * i - 1] * ans[2 * i];
          dum2 = dum * dum + ans[2 * i] * ans[2 * i];
          ans[2 * i] = ans1 / dum2 / no2;
		}
        else
		{
          cout<<" no meaning for isign"<<endl;
		}
    }
    ans[2] = ans[2 * no2 + 1];
    realft(ans, no2, temp);
    for(i=1; i<=32; i++)
		fft[i]=0.0;
}