d12r6.cpp

离散相关和自相关的快速算法

#include <iostream.h>
#include <math.h>
#include <iomanip.h>
#include <stdlib.h>
#include <fstream.h>
#include <string>
#include <process.h>
#include<stdio.h>

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

void realft(double data[], int &n, int &isign)
{
	int n2p3,i,i1,i2,i3,i4;
	double theta,c1,c2,wpr,wpi,wr,wi,wrs,wis,h1r,h1i,h2r,h2i,wtemp,wri;
	wri=0.0;
    theta = 6.28318530717959 / 2.0 / n;
    c1 = 0.5;
    if (isign == 1)
	{
        c2 = -0.5;
        four1(data, n, 1);
	}
    else
	{
        c2 = 0.5;
        theta = -theta;
	}
    wpr = -2.0 * sin(0.5 * theta)*sin(0.5 * theta);
    wpi = sin(theta);
    wr = 1.0 + wpr;
    wi = wpi;
    n2p3 = 2 * n + 3;
    for (i = 2;i<=n / 2 + 1;i++)
	{
        i1 = 2 * i - 1;
        i2 = i1 + 1;
        i3 = n2p3 - i2;
        i4 = i3 + 1;
        wrs = float(wr);
        wis = float(wi);
        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 + wrs * h2r - wis * h2i;
        data[i2] = h1i + wrs * h2i + wis * h2r;
        data[i3] = h1r - wrs * h2r + wis * h2i;
        data[i4] = -h1i + wrs * h2i + wis * h2r;
        wtemp = wr;
        wr = wr * wpr - wi * wri + wr;
        wi = wi * wpr + wtemp * wpi + wi;
    }
    if( isign == 1 )
	{
        h1r = data[1];
        data[1] = h1r + data[2];
        data[2] = h1r - data[2];
	}
    else
	{
        h1r = data[1];
        data[1] = c1 * (h1r + data[2]);
        data[2] = c1 * (h1r - data[2]);
        four1(data, n, -1);
    }
}

void twofft(double data1[], double data2[], double fft1[], double fft2[], int &n)
{
	int j,n2,j2;
	double c1r,c1i,c2r,c2i,conjr,conji,h1r,h1i,h2r,h2i;
    c1r = 0.5;
    c1i = 0.0;
    c2r = 0.0;
    c2i = -0.5;
    for (j = 1;j<=n;j++)
	{
        fft1[2 * j - 1] = data1[j];
        fft1[2 * j] = data2[j];
    }
    four1(fft1, n, 1);
    fft2[1] = fft1[2];
    fft2[2] = 0.0;
    fft1[2] = 0.0;
    n2 = 2 * (n + 2);
    for (j = 2;j<=n / 2 + 1;j++)
	{
        j2 = 2 * j;
        conjr = fft1[n2 - j2 - 1];
        conji = -fft1[n2 - j2];
        h1r = c1r * (fft1[j2 - 1] + conjr) - c1i * (fft1[j2] + conji);
        h1i = c1i * (fft1[j2 - 1] + conjr) + c1r * (fft1[j2] + conji);
        h2r = c2r * (fft1[j2 - 1] - conjr) - c2i * (fft1[j2] - conji);
        h2i = c2i * (fft1[j2 - 1] - conjr) + c2r * (fft1[j2] - conji);
        fft1[j2 - 1] = h1r;
        fft1[j2] = h1i;
        fft1[n2 - j2 - 1] = h1r;
        fft1[n2 - j2] = -h1i;
        fft2[j2 - 1] = h2r;
        fft2[j2] = h2i;
        fft2[n2 - j2 - 1] = h2r;
        fft2[n2 - j2] = -h2i;
    }
}

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;
}

void main()
{
    //program d12r6
    //driver for routine convlv
	int n,n2,m,i,j,isign;
	double data[17], respns[10], resp[17], ans[34],cmp;
    n = 16;
    n2 = 34;
    m = 9;
    const double pi = 3.1415926;   
    for (i = 1; i<=n; i++)
	{
        data[i] = 0.0;
        if (i >= (n / 2 - n / 8) && i <= (n / 2 + n / 8))
		{
            data[i] = 1.0;
        }
    }
    cout<<endl;
    for( i = 1; i<=m; i++)
	{
        respns[i] = 0.0;
        if (i > 2 && i < 7) 
			respns[i] = 1.0;
        resp[i] = respns[i];
    }
    isign = 1;
    convlv(data, n, resp, m, isign, ans);
    //compare with a direct convolution
    cout<<"   i       convlv       expected"<<endl;
    for( i = 1; i<=n; i++)
	{
        cmp = 0.0;
        for (j = 1; j<=m / 2; j++)
		{
          cmp = cmp + data[((i - j - 1 + n) % n) + 1] * respns[j + 1];
          cmp = cmp + data[((i + j - 1) % n) + 1] * respns[m - j + 1];
        }
        cmp = cmp + data[i] * respns[1];
        cout<<setw(4)<<i;
        cout<<setprecision(6)<<setw(14)<<ans[i];
        cout<<setprecision(6)<<setiosflags(ios::fixed)<<setw(14)<<cmp<<endl;
    }
}