tfloat.cpp

fft算法 在vc上的实现 是以基2dit算法实现的

//信道探测

#include "stdio.h" 
#include "string.h"
#include "math.h"
#include "stdlib.h"

const double pi=3.14159;
const  Num=512;
const  M=1024;
const  Len=8;


struct xi
     {
       double  real;
       double  imag;
     };

//****************************************均匀分布的随机数******************************

double  uniform(double a,double  b,long int *seed)
{
	double  t;
	*seed=2045*(*seed)+1;
	*seed=*seed-(*seed/1048576)*1048576;
	t=(*seed)/1048576;
	t=a+(b-a)*t;
	return(t);
}

//************************************产生高斯白噪声随机数*******************************

double gauss(double  mean,double sigma,long int *s)
{
	double  i;
	double  x,y;
	double  uniform(double,double,long int *);
	for(x=0,i=0;i<12;i++)
		x+=uniform(0,1,s);
	x=x-6;
	y=mean+x*sigma;
	return (y);
}

//************************************fir数字滤波器**************************************

void filtery (double b[],double a[],int m, int n, double *x,double *y,double lenth,double px[],double py[])
{
	int i,k;
	double sum;
	for (k=0;k<lenth;k++)
	{
		px[0]=*(x+k);
		sum=0;
		for (i=0;i<=m;i++)
		{
			sum=sum+b[i]*px[i];
		}
		for (i=1;i<=n;i++)
		{
			sum=sum-a[i]*py[i];
		}
		if (fabs(*(x+k))>1.0e10)
		{
			printf("This is an unstable filter!\n");
			exit(0);
		}
		for(i=m;i>=1;i--)
		{
			px[i]=px[i-1];
		}
		for(i=n;i>=2;i--)
		{
			py[i]=py[i-1];
		}
		py[1]=sum;
		y[k]+=sum;
	}
}

//***********************************产生有色噪声****************************************
void colnoise(double *noise, double *noi, long int nn)
{
    int m=1,n=0; 
	int i;
	
    double bb[2]={1,-0.9};
	double aa[1]={1};
	double px[2]={0,0};
	double py[1]={0};
	int lenth=nn;
    long int s=13579;

	for(i=0;i<nn;i++)
	{
		*(noise+i)=0;
	}
   
    for (i=0;i<nn;i++)
	{
		*(noi+i)=gauss(0,1,&s);
		
	}

	filtery(bb,aa,m,n,noi,noise,lenth,px,py);

}



 /*Function of Fast Fourier Transform*/
//void fft(struct xi *x,double n,double sign)
//{
	
//}

void main() 
{ 
    int i,j;
	int fr[Num];
	for (i=0;i<Num;i++)
	{
		fr[i]=2500+5000*i;
	}

	//计算sinad发送连续正弦波，持续时间为Long*200us
    double carry[M/2*Len];//正弦波
	double *carry1=(double *)malloc(sizeof(double)*M/2*Len*Num);
    double *sigsinad=(double *)malloc(sizeof(double)*M/2*Len*Num);//检测sinad时的发送信号
    struct xi *fsig=(struct xi *)malloc(sizeof(struct xi)*M/2*Len*Num);//信道探测，发送信号的频谱
    struct xi *sig1=(struct xi *)malloc(sizeof(struct xi)*Len*Num);
	double psig[Num];//计算信号功率



  
	//计算sinad时的发送信号
    for (i=0;i<M/2*Len*Num;i++)
	{
        *(sigsinad+i)=0;
		*(carry1+i)=0;
	}
	for (i=0;i<Num;i++)
	{
		psig[i]=0;

	}

    //连续发送正弦波
	for (i=0;i<Num;i++)
	{
		for (j=0;j<M/2*Len;j++)
		{
			carry[j]=sin(2*pi*fr[i]*j*0.000000390625);
            *(carry1+M/2*Len*i+j)=carry[j];
		}
   
        for (j=0;j<M/2*Len*Num;j++)
		{
			*(sigsinad+j)=*(sigsinad+j)+*(carry1+j);
			*(carry1+j)=0;
		}

	}
    
    //计算各频段内混合信号的功率
	for (i=0;i<M/2*Len*Num;i++)
	{
		(fsig+i)->real=*(sigsinad+i);
		(fsig+i)->imag=0;
	}


    
    int l,m,n2;
	int k,n1;
	double c,c1,s,s1,t,tr,ti;
	double e;
	for(j=1,i=1;i<25;i++)
	{
		m=i;
		j=2*j;
		if(j==Num*M/2*Len)
			break;
	}
	n1=Num*M/2*Len-1;
	for(j=0,i=0;i<n1;i++)
	{
		if(i<j)
		{
			tr=(fsig+j)->real;
			ti=(fsig+j)->imag;
			(fsig+j)->real=(fsig+i)->real;
			(fsig+j)->imag=(fsig+i)->imag;
			(fsig+i)->real=tr;
			(fsig+i)->imag=ti;
		}
		k=Num*M/2*Len/2;
		while(k<(j+1))
		{
			j=j-k;
			k=k/2;
		}
		j=j+k;
	}

	n1=1;
	for (l=1;l<=13;l++)
	{
		n1=2*n1;
		n2=n1/2;
		e=pi/n2;
		c=1;
		s=0;
		c1=cos(e);
		s1=-1*sin(e);
        

		for(j=0;j<n2;j++)
		{
			for(i=j;i<Num*M/2*Len;i+=n1)
			{
				k=i+n2;
				tr=c*(fsig+k)->real-s*(fsig+k)->imag;
				ti=c*(fsig+k)->imag+s*(fsig+k)->real;
				(fsig+k)->real=(fsig+i)->real-tr;
				(fsig+k)->imag=(fsig+i)->imag-ti;
				(fsig+i)->real=(fsig+i)->real+tr;
				(fsig+i)->imag=(fsig+i)->imag+ti;
			}

			t=c;
			c=t*c1-s*s1;
			s=t*s1+s*c1;
            printf("%f ",c);
		}

		printf("\n");
		//printf("%f %f ",(fsig+k)->real,(fsig+k)->imag);
				printf("\n");

	}



   /* for (i=0;i<Num;i++)
	{
		for (j=0;j<Len*Num;j++)
		{
			(sig1+j)->real=(fsig+j+i*Len*Num)->real;
			(sig1+j)->imag=(fsig+j+i*Len*Num)->imag;
		}
		for (j=0;j<Len*Num;j++)
		{
            psig[i]=psig[i]+((sig1+j)->real*(sig1+j)->real+(sig1+j)->imag*(sig1+j)->imag)/(Num*M/2*Len);
		}
	
	}*/

		

	free(carry1);
	free(sigsinad);
	free(fsig);
	free(sig1);

}