t1.cpp

同样是浮点型的fft算法一样为vc开发的

//信道探测

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

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


struct xi
     {
       short real;
       short imag;
     };


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

	//计算sinad发送连续正弦波，持续时间为Long*200us
    short carry[M/2*Len];//正弦波
	short *carry1=(short *)malloc(sizeof(short)*M/2*Len*Num);
    short *sigsinad=(short *)malloc(sizeof(short)*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);
	int 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]=(short)(sin(2*pi*fr[i]*j*0.000000390625)*pow(2,15));
            *(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,k,n1,n2;
	double e;
	short c,c1,s,s1,t,tr,ti;
	short *p=(short *)malloc(sizeof(short)*M/2*Len*Num/2);
    short *p1=(short *)malloc(sizeof(short)*M/2*Len*Num/2);
	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<=m;l++)
	{
		n1=2*n1;
		n2=n1/2;
		e=pi/n2;
		c=1;
		s=0;
		c1=(short)(cos(e)*pow(2,15));
		s1=-1*(short)(sin(e)*pow(2,15));
 	    *p=32767;
		*p1=0;


		for(i=0;i<Num*M/2*Len;i+=n1)
		{
				k=i+n2;
				tr=(fsig+k)->real;
				ti=(fsig+k)->imag;
				(fsig+k)->real=((long)(fsig+i)->real-tr)>>1;
				(fsig+k)->imag=((long)(fsig+i)->imag-ti)>>1;
				(fsig+i)->real=((long)(fsig+i)->real+tr)>>1;
				(fsig+i)->imag=((long)(fsig+i)->imag+ti)>>1;
		}
            
		c=c1;
		s=s1;
		*(p+1)=c;
		*(p1+1)=s;
		printf("%d ",c);

		for(j=1;j<n2/2;j++)
		{
			for(i=j;i<Num*M/2*Len;i+=n1)
			{
					k=i+n2;
					tr=((long)c*(fsig+k)->real-(long)s*(fsig+k)->imag)>>16;
					ti=((long)c*(fsig+k)->imag+(long)s*(fsig+k)->real)>>16;
					(fsig+k)->real=((fsig+i)->real>>1)-tr;
					(fsig+k)->imag=((fsig+i)->imag>>1)-ti;
					(fsig+i)->real=((fsig+i)->real>>1)+tr;
					(fsig+i)->imag=((fsig+i)->imag>>1)+ti;
			}
			
			t=c;
			c=((long)t*c1-(long)s*s1)>>15;
			s=((long)t*s1+(long)s*c1)>>15;
			*(p+j+1)=c;
			*(p1+j+1)=s;
			printf("%d ",c);
		}

		if(n2/2>0)
		{
		for(j=n2/2;j<n2;j++)
		{
			c=-*(p+n2-j);
			s=*(p1+n2-j);
			for(i=j;i<Num*M/2*Len;i+=n1)
			{
					k=i+n2;
					tr=((long)c*(fsig+k)->real-(long)s*(fsig+k)->imag)>>16;
					ti=((long)c*(fsig+k)->imag+(long)s*(fsig+k)->real)>>16;
					(fsig+k)->real=((fsig+i)->real>>1)-tr;
					(fsig+k)->imag=((fsig+i)->imag>>1)-ti;
					(fsig+i)->real=((fsig+i)->real>>1)+tr;
					(fsig+i)->imag=((fsig+i)->imag>>1)+ti;
			}
		
			printf("%d ",c);
		}
		}		
		printf("\n");
        //printf("%d, %d ", (fsig+k)->real,(fsig+k)->imag);
	    printf("\n");
	}
	for (i=0;i<Num*M/2*Len;i++)
	{
	    (fsig+i)->real=(fsig+i)->real<<6;
	    (fsig+i)->imag=(fsig+i)->imag<<6;
	}

    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]+(long)(sig1+j)->real*(sig1+j)->real+(long)(sig1+j)->imag*(sig1+j)->imag;
		}
		psig[i]=psig[i]/(Num*M/2*Len);
	
	}

    //for(i=0;i<Num;i++)
		//printf("%d ",psig[i]);

	free(carry1);
	free(sigsinad);
	free(fsig);
	free(sig1);
	free(p);
	free(p1);
}