fcomplex.h

用Visual C++编写的地震记录显示程序

#include "math.h"

class fcomplex
{ 
public:
	float r , i ;
	fcomplex operator + (fcomplex a ) ;
	fcomplex operator - ( fcomplex a ) ;
	fcomplex operator *( fcomplex a ) ;
	fcomplex operator / (fcomplex b) ;
} ;

fcomplex fcomplex::operator + ( fcomplex a ){   /////两复数相加
	fcomplex c ;
	c.r = a.r + r ;
	c.i = a.i + i ;
	return c ;
}

fcomplex fcomplex::operator - ( fcomplex a ){   //////两复数相减
	fcomplex c ;
	c.r = r - a.r ;
	c.i = i - a.i ;
	return c ;
}

fcomplex fcomplex::operator *( fcomplex a ){ ////////两复数相乘
	fcomplex c ;
	c.r = a.r * r - a.i * i ;
	c.i = a.i * r + a.r * i ;
	return c ;
}

fcomplex Complex( float re , float im )  //////生成复数
{
	fcomplex c ;
	c.r = re ;
	c.i = im ;
	return c ;
}
 
fcomplex Complex( float re )  //////生成复数
{
	fcomplex c ;
	c.r = re ;
	c.i = 0 ;
	return c ;
}
 
fcomplex Complex(  )  //////生成复数
{
	fcomplex c ;
	c.r = 0 ;
	c.i = 0 ;
	return c ;
}

fcomplex Conjg( fcomplex z ){ //////////计算复数的共轭
	fcomplex c ;
	c.r = z.r ;
	c.i = - z.i ;
	return c ;
}

fcomplex fcomplex::operator / (fcomplex b){ /////////计算两复数相除
	fcomplex c ;
	float r , den ;
	if( fabs(b.r) >= fabs(b.i)){
		r = b.i / b.r ;
		den = b.r + r * b.i ;
		c.r = ( r + r * i ) / den ;
		c.i = ( i - r * r ) / den ;
	}
	else{
		r = b.r / b.i ;
		den = b.i + r * b.r ;
		c.r = ( r * r + i ) / den ;
		c.i = ( i * r - r ) / den ;
	}
	return c ;
}

float real(fcomplex z)
{
	return z.r;
}

float imag(fcomplex z)
{
	return z.i;
}

float Cabs(fcomplex z){ ////////计算复数的模
	float x , y , ans , temp ;
	x = fabs( z.r ) ;
	y = fabs( z.i ) ;
	if( x == 0 )
		ans = y ;
	else if( y == 0 )
		ans = x ;
	else if( x > y ){
		temp = y / x ;
		ans = x * sqrt( 1.0 + temp * temp ) ;
	}
	else{
		temp = x / y ;
		ans = y * sqrt( 1.0 + temp * temp) ;
	}
	return ans ;
}

fcomplex Csqrt( fcomplex z ){ ////////计算复数的二次方根
	fcomplex c ;
	float x , y , w , r ;
	if((z.r == 0.0) && (z.i == 0.0)){
		c.r = 0.0 ;
		c.i = 0.0 ;
		return c ;
	}
	else{
		x = fabs( z.r ) ;
		y = fabs( z.i ) ;
		if( x >= y ){
			r = y / x ;
			w = sqrt(x) * sqrt( 0.5 *(1.0 + sqrt(1.0 + r * r))) ;
		}
		else{
			r = x / y ;
			w = sqrt( y ) * sqrt( 0.5 * (r + sqrt( 1.0 + r * r ))) ;
		}
		if( z.r >= 0.0 ){
			c.r = w ;
			c.i = z.i / ( 2.0 * w ) ;
		}
		else{
			c.i = ( z.i >= 0 ) ? w : -w ;
			c.r = z.i / ( 2.0 * c.i ) ;
		}
		return c ;
	}
}

fcomplex RCmul( float x , fcomplex a ){      ////////计算实数与复数相乘
	fcomplex c ;
	c.r = x * a.r ;
	c.i = x * a.i ;
	return c ;
}

void FFT( fcomplex *a , fcomplex* b , int m)
{  ////////计算数组a[]的富氏变换存入数组b[]中， m为以2为底的次数
	int q,r1,m1,k,k1,k2,i,j,n;   //////////////
	float pi = 3.1415926535 ;
	n=pow(2,m) ;
	fcomplex w,W[10];
	w.r = cos( 2 * pi / n ) ; /////////// w = exp( - i * 2 * pi / n )  = cos( 2 * pi / n ) - i * sin( 2 * pi / n ) 
	w.i = -sin( 2 * pi / n ) ;
	for(i=0;i<n;i++)
	{
	//	cout<<a[i].r<<","<<a[i].i<<endl;
	}
	W[0].r=1;
	W[0].i=0;
	for(i=1;i<n/2;i++)
	{
		W[i]=W[i-1]*w;
	}
	for(q=1;q<=m;q++)
	{
		r1 = pow( 2 , q - 1 ) ;
		m1 = pow( 2 , m - 1 ) ;
		for(k=0;k<=pow(2,m-q)-1;k++)
		{
			k1 = k * r1 ;
				
			for(j=0;j<=pow(2,q-1)-1;j++)
			{
				if(q%2!=0)
				{
					b[2*k1+j]=a[k1+j]+a[k1+j+m1];
					b[2*k1+j+r1]=(a[k1+j]-a[k1+j+m1])*W[k1];
				}
				else
				{
					a[2*k1+j]=b[k1+j]+b[k1+j+m1];
					a[2*k1+j+r1]=(b[k1+j]-b[k1+j+m1])*W[k1];
				}
			}
		}
	}
	if(m%2==0 )
	{
		for(i=1;i<n;i++)
			b[i]=a[i] ;
	}
	for(i=0;i<n;i++)
	{
	//	cout<<b[i].r<<","<<b[i].i<<endl;
	}

}

void IFFT( fcomplex* a , fcomplex* b , int m){ ////////计算数组a[]的富氏反变换存入数组b[]中， m为以2为底的次数
	int r , r1 , m1 , k , k1 , k2 ,  i , j , ii ,n ;
	float pi = 3.1415926535 ;
	n = pow( 2 , m ) ;
	fcomplex w , ans;
	w.r = cos( 2 * pi / n ) ;
	w.i = sin( 2 * pi / n ) ;
	for( r = 1 ; r <= m ; r ++ )
	{
		r1 = pow( 2 , m - r ) ;
		m1 = pow( 2 , m - 1 ) ;
		for( k = 0 ; k < (float)m1 / r1 ; k ++)
		{
			k1 = k * r1 ;
			k2 = 2 * k1 + r1 ;
			if( r % 2 != 0 )
			{
				for( j = 0 ; j < r1 ; j ++ )
				{
					ans.i = 0 ;
					ans.r = 1 ;
					for( ii = 0 ; ii < k1 ; ii ++)
						ans = ans * w ;
					b[k1 + j] = a[ 2 * k1 + j] + a[k2 + j] * ans ;
					b[ k1 + j + m1 ] = a[ 2* k1 + j ] - a[ k2 + j ] * ans ;
				}
			}
			else
			{
				for( j = 0 ; j < r1 ; j ++ )
				{
					ans.i = 0 ;
					ans.r = 1 ;
					for( ii = 0 ; ii < k1 ; ii ++)
						ans = ans * w ;
					a[ k1 + j ] = b[ 2*k1 + j ] + b[ k2 + j ] * ans ;
					a[ k1 + j + m1 ] = b[ 2 * k1 + j ] - b[ k2 + j ] * ans ;
				}
			}
		}
	}
	if( m % 2 == 0 )
	{
		for( i = 0 ; i < n ; i ++)
			b[ i ] = a[ i ] ;
	}
	for( i = 0 ; i < n ; i ++ )
	{
		b[i].r /= n ;
		b[i].i /= n ;
	}
}