Windows编程

开发平台：
Visual C++

Pre_deconvolution.cpp：源码内容
							#include "iostream.h"
#include "math.h"
//#include <stdlib.h>
#define M 20
#define N 350
double ** pre_deconv(double **y,int m,int n,int ln,int da);
double ** zeros(int m,int n);
double *xxcorr(double *Y,int n);
double *div(double **A,double *b,int num);
double *conv(double *a,double *b,int num);
void main()
{
    double **p;
	int i,j;
   p= new double*[2];
	for(i=0;i<2;i++)
	   p[i]=new double[20];
  
  p=zeros(M,N);
  double **f=pre_deconv(p,M,N,20,1);
  for(i=0;i<2;i++)
	   {  cout<<endl;
	   for(j=0;j<10;j++)
		   cout<<f[i][j]<<"	";
	   }
	   
for(i=0;i<M;i++)
		   delete[] p[i];
	   delete[] p;
}
double ** pre_deconv(double **y,int m,int n,int ln,int da)
{
	
    int i,j,k;
	int yn=0;
	
	double **py=new double*[m];
	for(i=0;i<m;i++)
	 py[i]=new double[n];
	double **Dr;
	double *yy=new double[ln];
	while(yn<(n-n%ln))
	{
		yn=yn+ln;
		Dr=zeros(m,ln);
		for(i=0;i<m;i++)
		{
			for(k=0;k<ln;k++)
				yy[k]=y[i][yn+k];
			Dr[i]=xxcorr(yy,ln);
		}
		double *DR= new double[ln];
		for(i=0;i<ln;i++)
		{	
			DR[i]=0;
			for(j=0;j<m;j++)
			DR[i]+=Dr[j][i];
		}
	
		double **A;       //Toeplitz矩阵
		
		A=zeros(ln,ln);
		
		for(i=0;i<ln;i++)
		 for(j=0;j<ln;j++)
		{
			if(j==i)A[i][j]=DR[0];
		else if(j>i)A[i][j]=DR[j-i];	
			else A[i][j]=DR[i-j];	
		}
		
		 double *b=new double[ln];
		 for(i=0;i<ln;i++)
		 {
			 if((da+i)<ln)b[i]=DR[da+i];
			 else b[i]=0;
		 }
         double *x=div(A,b,ln);		//求Ax=b, x=Ab;
		 double *C=new double[ln+da];//求反滤波因子C[1,0...0,-c1,-c2,..]
		 C[0]=1;
		 for(i=0;i<ln+da;i++)
		 {
			 if(i<da)C[i]=0;
			else C[i]=-x[i-da];
         }
///////////////////////////////////做卷积...	
		 
		 for(i=0;i<m;i++)
		 {
			 double *s=conv(Dr[i],C,ln);
			 for(j=0;j<ln;j++)
				 py[i][yn-ln+j]=s[j];
		 }
	}
		
	  
		  return py;
}
double ** zeros(int m,int n)
{
	
	int i,j;
    double  **a= new double*[m];
	for(i=0;i<m;i++)
		a[i]=new double[n];
	
	for(i=0;i<m;i++)
		for(j=0;j<n;j++)
			a[i][j]=0;
		return a;
		
}
double *xxcorr(double *Y,int n)
{
	int i,j;
	double *X = new double[n];
	for(i=0;i<n;i++)
		{
			double xk=0;
			for(j=0;j+i<n;j++)
			xk+=Y[j]*Y[j+i];
		
			X[i]=xk;
		}
		return X;
}
double *div(double **A,double *b,int num)
{
	int i,j,k;
//共处理num-1行
	double *x0= new double[num];
	for(i=0;i<num-1;i++)
	{
		//每列选主元
		double lineMax=fabs(A[i][i]);
		int lineI=i;
		for(j=i;j<num;j++)
		if(fabs(A[j][i])>fabs(lineMax))lineI=j;
		//主元所在行和当前行交换
		for(j=i;j<num;j++)
		{	
			//
			lineMax=A[i][j];
			A[i][j]=A[lineI][j];
			A[lineI][j]=lineMax;
	
			lineMax=b[lineI];
			b[i]=b[lineI];
			b[lineI]=lineMax;
		}
		
		if(A[i][i]==0)continue;
		//将A变为上三角
		for(j=i+1;j<num;j++)
		{
			double temp=-A[j][i]/A[i][i];
			for(k=i;k<num;k++)
			{
				A[j][k]+=temp*A[i][k];
			}
			b[j]+=temp*b[i];
		}
		
	}
	double determinantA=1;
	for(i=0;i<num;i++)
		determinantA*=A[i][i];
	
	
	
//回代求解x初值
		x0[num-1]=0;
	for(i=num-1;i>=0;i--)
	{
		for(j=num-1;j>i;j--)
			b[i]-=A[i][j]*x0[j];
		x0[i]=b[i]/A[i][i];
		
	}
//////
	return x0;
	
}
double *conv(double *a,double *b,int num)
{
	int i,t;
	double *c1=new double[num];
    for(i=0;i<num;i++) 
	  
	 {
		 
	      c1[i]=0;
	  
	     for(t=num;t>=0;t--)                 //	做卷积
			if(t<=i)c1[i]+=*(a+t)*(*(b+i-t));
			else continue;
     }	
	
	return c1;
}

						