matrixoper.cpp

矩阵运算处理,包括FFT,IFFT,FFTSHIFT等

#include "stdafx.h"
#include "math.h"
COMPLEX CMatrixOper::Add(COMPLEX c1,COMPLEX c2)
{
	COMPLEX   c;   
	c.re=c1.re+c2.re;   
	c.im=c1.im+c2.im;   
	return   c; 	
}

COMPLEX CMatrixOper::Sub(COMPLEX c1,COMPLEX c2)
{
	COMPLEX   c;   
	c.re=c1.re-c2.re;   
	c.im=c1.im-c2.im;   
	return   c; 
}

COMPLEX CMatrixOper::Mul(COMPLEX c1,COMPLEX c2)
{
	COMPLEX   c;   
	c.re=c1.re*c2.re-c1.im*c2.im;   
	c.im=c1.re*c2.im+c2.re*c1.im;   
	return   c;  	
}

void CMatrixOper::Convert(COMPLEX*c1)
{
	(*c1).im = -(*c1).im;	
}   
/*  参数：   
    
  FD为频域值   
  TD为时域值   
  power为2的幂数   
    
  返回值：   
    
  无   
    
  说明：   
    
  本函数利用快速傅立叶变换实现傅立叶变换   
  ****************************************************/    
void CMatrixOper::FFT(COMPLEX*TD,COMPLEX*FD,int power)
{
	int   count;   
	int   i,j,k,bfsize,p;   
	double   angle;   
	COMPLEX   *W,*X1,*X2,*X;   
    
	/*计算傅立叶变换点数*/   
	count=1<<power;   
    
	/*分配运算所需存储器*/   
	W = (COMPLEX *)malloc(sizeof(COMPLEX)*count/2);   
	X1 = (COMPLEX *)malloc(sizeof(COMPLEX)*count);   
	X2 = (COMPLEX *)malloc(sizeof(COMPLEX)*count);   
    
	/*计算加权系数*/   
	for(i=0;i<count/2;i++)   
	{   
		angle = -i*PI*2/count;   
		W[i].re = cos(angle);   
		W[i].im = sin(angle);   
	}   
    
	/*将时域点写入存储器*/   
	memcpy(X1,TD,sizeof(COMPLEX)*count);   
    
	/*蝶形运算*/   
	for(k = 0;k<power;k++)   
	{   
		for(j = 0;j<1<<k;j++)   
		{   
			bfsize = 1<<(power-k);   
			for(i = 0;i<bfsize/2;i++)   
			{   
				p = j*bfsize;   
				X2[i+p] = Add(X1[i+p],X1[i+p+bfsize/2]);   
				X2[i+p+bfsize/2] = Mul(Sub(X1[i+p],X1[i+p+bfsize/2]),W[i*(1<<k)]);   
			}   
		}   
		X = X1;   
		X1 = X2;   
		X2 = X;   
	}   
    
	/*重新排序*/   
	for(j = 0;j < count;j++)   
	{   
		p = 0;   
		for(i = 0;i<power;i++)   
		{   
			if(j&(1<<i))   
			p+=1<<(power-i-1);   
		}   
		FD[j] = X1[p];   
	}   
    
	/*释放存储器*/   
	free(W);   
	free(X1);   
	free(X2);   
}
/*  参数：   
    
  FD为频域值   
  TD为时域值   
  power为2的幂数   
    
  返回值：   
    
  无   
    
  说明：   
    
  本函数利用快速傅立叶变换实现傅立叶反变换   
  ****************************************************/  
void CMatrixOper::IFFT(COMPLEX*FD,COMPLEX*TD,int power)
{
	int   i,count;   
	COMPLEX *x;   
    
	/*计算傅立叶反变换点数*/   
	count = 1<<power;   
    
	/*分配运算所需存储器*/   
	x = (COMPLEX *)malloc(sizeof(COMPLEX)*count);   
    
	/*将频域点写入存储器*/   
	memcpy(x,FD,sizeof(COMPLEX)*count);   
    
	/*求频域点的共轭*/   
	for(i = 0 ;i < count; i++)   
		x[i].im = -x[i].im;   
    
	/*调用FFT*/   
	FFT(x,TD,power);   
    
	/*求时域点的共轭*/   
	for(i = 0;i<count;i++)   
	{   
		TD[i].re /= count;   
		TD[i].im = -TD[i].im / count;   
	}   
    
	/*释放存储器*/   
	free(x);     
}

void CMatrixOper::MatrCOMRePointMulMatr(double*Matr1,int row,int colum,COMPLEX*Matr2,COMPLEX*Matr)
{
	for(int i = 0;i < row*colum;i++)
	{
		(*(Matr + i)).re = *(Matr1 + i) * (*(Matr2 + i)).re;
	}
}

void CMatrixOper::MatrCOMImPointMulMatr(double*Matr1,int row,int colum,COMPLEX*Matr2,COMPLEX*Matr)
{
	for(int i = 0;i < row*colum;i++)
	{
		(*(Matr + i)).im = *(Matr1 + i) * (*(Matr2 + i)).im;
	}
}

void CMatrixOper::MatrPointMul(double*Matr1,int row,int colum,double*Matr2,double*Matr)
{
	for(int i = 0;i < row*colum;i++)
	{
		*(Matr + i) = *(Matr1 + i) * *(Matr2 + i);
	}
}

void CMatrixOper::MatrPointDiv(double*Matr1,int row,int colum,double*Matr2,double*Matr)
{
	for(int i = 0;i < row*colum;i++)
	{
		*(Matr + i) = *(Matr1 + i) / *(Matr2 + i);
	}
}

void CMatrixOper::MatrPointAdd(double*Matr1,int row,int colum,double*Matr2,double*Matr)
{
	for(int i = 0;i < row*colum;i++)
	{
		*(Matr + i) = *(Matr1 + i) + *(Matr2 + i);
	}
}

void CMatrixOper::MatrPointSub(double*Matr1,int row,int colum,double*Matr2,double*Matr)
{
	for(int i = 0;i < row*colum;i++)
	{
		*(Matr + i) = *(Matr1 + i) - *(Matr2 + i);
	}
}

void CMatrixOper::MatrreAddim(COMPLEX*Matr,int row,int colum)
{
	for(int i = 0;i < row*colum;i++)
	{
		(*(Matr + i)).re = (*(Matr + i)).re + (*(Matr + i)).im;
	}
}
void CMatrixOper::MatrPointAddNum(double*Matr,int row,int colum,double a)
{
	for(int i = 0;i < row*colum;i++)
	{
		*(Matr + i) = *(Matr + i) + a;
	}
}
void CMatrixOper::MatrPointMulNum(double*Matr,int row,int colum,double a)
{
	for(int i = 0;i < row*colum;i++)
	{
		*(Matr + i) = *(Matr + i) * a;
	}
}
void CMatrixOper::MatrPointMulNumre(COMPLEX*Matr,int row,int colum,double a)
{
	for(int i = 0;i < row*colum;i++)
	{
		(*(Matr + i)).re = (*(Matr + i)).re * a;
	}
}

void CMatrixOper::MatrPointMulNumim(COMPLEX*Matr,int row,int colum,double a)
{
	for(int i = 0;i < row*colum;i++)
	{
		(*(Matr + i)).im = (*(Matr + i)).im * a;
	}
}
void CMatrixOper::MatrPointMulNumretoim(COMPLEX*Matr1,int row,int colum,double a,COMPLEX*Matr)
{
	for(int i = 0;i < row*colum;i++)
	{
		(*(Matr + i)).im = (*(Matr1 + i)).re * a;
	}
}

void CMatrixOper::MatrGetre(double*Matr1,int row,int colum,COMPLEX*Matr)
{
	for(int i = 0;i < row*colum;i++)
	{
		(*(Matr + i)).re = cos(*(Matr1 + i));
	}
}
void CMatrixOper::MatrGetim(double*Matr1,int row,int colum,COMPLEX*Matr)
{
	for(int i = 0;i < row*colum;i++)
	{
		(*(Matr + i)).im = sin(*(Matr1 + i));
	}
}
void CMatrixOper::MatrPointMulNumimtore(COMPLEX*Matr1,int row,int colum,double a,COMPLEX*Matr)
{
	for(int i = 0;i < row*colum;i++)
	{
		(*(Matr + i)).re = (*(Matr1 + i)).im * a;
	}
}
void CMatrixOper::MatrPointPow2(double*Matr1,int row,int colum,double*Matr)
{
	for(int i = 0;i < row*colum;i++)
	{
		*(Matr + i) = *(Matr1 + i) * *(Matr1 + i);
	}
}
void CMatrixOper::MatrPointPow2re(COMPLEX*Matr1,int row,int colum,COMPLEX*Matr)
{
	for(int i = 0;i < row*colum;i++)
	{
		(*(Matr + i)).re = (*(Matr1 + i)).re * (*(Matr1 + i)).re;
	}
}

void CMatrixOper::MatrCOMreEquMatr(COMPLEX*Matr1,int row,int colum,double*Matr)
{
	for(int i = 0;i < row*colum;i++)
	{
		*(Matr + i) = (*(Matr1 + i)).re;
	}
}

void CMatrixOper::MatrCOMimEquMatr(COMPLEX*Matr1,int row,int colum,double*Matr)
{
	for(int i = 0;i < row*colum;i++)
	{
		*(Matr + i) = (*(Matr1 + i)).im;
	}
}
void CMatrixOper::MatrEqu(double*Matr1,int row,int colum,double*Matr)
{
	for(int i = 0;i < row*colum;i++)
	{
		*(Matr + i) = *(Matr1 + i);
	}
}
void CMatrixOper::MatrPointPow2im(COMPLEX*Matr1,int row,int colum,COMPLEX*Matr)
{
	for(int i = 0;i < row*colum;i++)
	{
		(*(Matr + i)).im = (*(Matr1 + i)).im * (*(Matr1 + i)).im;
	}
}
void CMatrixOper::MatrPointSqrt(double*Matr1,int row,int colum,double*Matr)
{
	for(int i = 0;i < row*colum;i++)
	{
		*(Matr + i) = sqrt(*(Matr1 + i));
	}
}

void CMatrixOper::MatrPointcos(double*Matr,int row,int colum)
{
	for(int i = 0;i < row*colum;i++)
	{
		*(Matr + i) = cos(*(Matr + i));
	}
}

void CMatrixOper::MatrPoinJudge(double*Matr,int row,int colum,double range1,double range2)
{
	for(int i = 0;i < row*colum;i++)
	{	
		if ((*(Matr + i) >= range1)&&(*(Matr + i)<=range2)) 
		{
			*(Matr + i) = 1;
		}
		else
		{
			*(Matr + i) = 0;
		}
	}
}

void CMatrixOper::MatrJudge(double*Matr1,int row,int colum,double*Matr)
{
	for(int i = 0;i < row*colum;i++)
	{	
		if ((*(Matr1 + i)<=*(Matr + i))) 
		{
			*(Matr + i) = 1;
		}
		else
		{
			*(Matr + i) = 0;
		}
	}
}

void CMatrixOper::MatrPoinJudgeSmal(double*Matr,int row,int colum,double range)
{
	for(int i = 0;i < row*colum;i++)
	{	
		if (*(Matr + i)<range) 
		{
			*(Matr + i) = 1;
		}
		else
		{