fft.cpp

在vc上做的802.16d ofdm phy的仿真

#include "math.h"
//#include "stdio.h"
#include "fft.h"
//#include "qam64.h"


void ChangeOrder (Int16 *changeDataIn, Int16 length, Int16 pow_4)
{
	Int16 i, j, k;
	Int16 temp1, temp2;
	Int16 temp3;
	Uint16 Index[10];
	for (i = 1; i < length - 1; i++)
	{
		temp1 = i;
		temp2 = 0;
		for (j = pow_4 - 1; j > 0; j--)
		{
			Index[j] = temp1 / int (pow (4, j) ) ;
			temp1 = temp1 % int (pow (4, j) );
		}
		Index[0] = temp1;
		for (k = 0; k < pow_4; k++)
		{
			temp2 = temp2 + Index[k] * int (pow (4, pow_4 - k - 1));
		}
		if (i < temp2)
		{
			temp3 = changeDataIn[2*i];
			changeDataIn[2*i] = changeDataIn[2*temp2];
			changeDataIn[2*temp2] = temp3;
			temp3 = changeDataIn[2*i+1];
			changeDataIn[2*i+1] = changeDataIn[2*temp2+1];
			changeDataIn[2*temp2+1] = temp3;

		}
	}
}
void fft (Int16 *pFftIn , 
		  Uint16  fftSize,
		  Int16 *pWeight,
		  Uint8  weightBitWidth,
		  Int16  *pFftOut)    
{                                                                       
	Int32	n1, n2, ie, ia1, ia2, ia3, i0, i1, i2, i3, j, k;    
    Int16   t, r1, r2, s1, s2, co1, co2, co3, si1, si2, si3;    
    Int16	pow_4;
	Int16	*x;
	Int16	iMaxValue;
	iMaxValue = 1 << (weightBitWidth - 1);
	
/*	for (j = 0; j < 2*fftSize; j++)
	{
		if (pWeight[j] >= iMaxValue)
		{
			pWeight[j] = iMaxValue;
		}
	}
	*/
	x = new Int16[2*fftSize];
	for (j = 0; j < fftSize; j++)
	{
		x[2*j] = pFftIn[2*j];
		x[2*j+1] = pFftIn[2*j+1];
	}
	pow_4 = (Int16)(log(fftSize) / log(4) + 0.5);                                                                    
    n2 = fftSize;                                                             
    ie = 1;                                                             
    for (k = fftSize; k > 4; k >>= 2)                                         
    {                                                                   
        n1 = n2;                                                        
        n2 >>= 2;                                                       
        ia1 = 0;                                                        
        for (j = 0; j < n2; j++)                                        
        {                                                               
            ia2 = ia1 + ia1;                                            
            ia3 = ia2 + ia1;                                            
            co1 = pWeight[ia1 * 2 + 1];                                       
            si1 = pWeight[ia1 * 2];                                           
            co2 = pWeight[ia2 * 2 + 1];                                       
            si2 = pWeight[ia2 * 2];                                           
            co3 = pWeight[ia3 * 2 + 1];                                       
            si3 = pWeight[ia3 * 2];                                           
            ia1 = ia1 + ie;                                             
            for (i0 = j; i0 < fftSize; i0 += n1)                              
            {                                                           
                i1 = i0 + n2;                                           
                i2 = i1 + n2;                                           
                i3 = i2 + n2;                                           
                r1 = x[2 * i0] + x[2 * i2];                             
                r2 = x[2 * i0] - x[2 * i2];                             
                t = x[2 * i1] + x[2 * i3];                              
                x[2 * i0] = (r1 + t + 1) >> 1;                                     
                r1 = r1 - t;                                            
                s1 = x[2 * i0 + 1] + x[2 * i2 + 1];                     
                s2 = x[2 * i0 + 1] - x[2 * i2 + 1];                     
                t = x[2 * i1 + 1] + x[2 * i3 + 1];                      
                x[2 * i0 + 1] = ( s1 + t + 1) >> 1;                                 
                s1 = s1 - t;                                            
                x[2 * i2] = (r1 * co2 + s1 * si2 + 0x00008000) >> 16;                
                x[2 * i2 + 1] = (s1 * co2-r1 * si2 + 0x00008000) >>16;                
                t = x[2 * i1 + 1] - x[2 * i3 + 1];                      
                r1 = r2 + t;                                            
                r2 = r2 - t;                                            
                t = x[2 * i1] - x[2 * i3];                              
                s1 = s2 - t;                                            
                s2 = s2 + t;                                            
                x[2 * i1] = (r1 * co1 + s1 * si1 + 0x00008000)  >>16;                
                x[2 * i1 + 1] = (s1 * co1-r1 * si1 + 0x00008000)>>16;                
                x[2 * i3] = (r2 * co3 + s2 * si3 + 0x00008000)  >>16;                
                x[2 * i3 + 1] = (s2 * co3-r2 * si3 + 0x00008000)>>16;                
            }                                                           
        }                                                               
        ie <<= 2;   
    }            
	for (i0 = 0; i0 < fftSize; i0 += 4)                              
	{                                                           
        i1 = i0 + 1;                                           
        i2 = i1 + 1;                                           
        i3 = i2 + 1;                                           
        r1 = x[2 * i0] + x[2 * i2];                             
        r2 = x[2 * i0] - x[2 * i2];                             
        t = x[2 * i1] + x[2 * i3];                              
        x[2 * i0] = r1 + t;                                     
        r1 = r1 - t;                                            
        s1 = x[2 * i0 + 1] + x[2 * i2 + 1];                     
        s2 = x[2 * i0 + 1] - x[2 * i2 + 1];                     
        t = x[2 * i1 + 1] + x[2 * i3 + 1];                      
        x[2 * i0 + 1] = s1 + t;                                 
        s1 = s1 - t;                                            
        x[2 * i2] = r1;                
        x[2 * i2 + 1] = s1;                
        t = x[2 * i1 + 1] - x[2 * i3 + 1];                      
        r1 = r2 + t;                                            
        r2 = r2 - t;                                            
        t = x[2 * i1] - x[2 * i3];                              
        s1 = s2 - t;                                            
        s2 = s2 + t;                                            
        x[2 * i1] = r1 ;                
        x[2 * i1 + 1] = s1 ;                
        x[2 * i3] = r2 ;                
        x[2 * i3 + 1] = s2 ;                
    }                                     
	  
	ChangeOrder (x, fftSize, pow_4);

	for (j = 0; j < fftSize; j++)
	{
		pFftOut[2*j] = x[2*j];
		pFftOut[2*j+1] = x[2*j+1];
	}
	delete[]x;
}    
void iFft (Int16 *pIfftIn, 
		   Int16  iFftSize,
		   Int16 *pWeight,
		   Uint8 weightBitWidth,
		   Int16  *pIfftOut)
{
	int i;
	Int16 *x;
	x = new Int16[2*iFftSize];
	for (i = 0; i < iFftSize; i++)
	{
		x[2*i] = pIfftIn[2*i];
		x[2*i+1] = -pIfftIn[2*i+1];
	}
    fft (x, iFftSize, pWeight, weightBitWidth, pIfftOut);
	for (i = 0; i < iFftSize; i++)
	{
		pIfftOut[2*i+1] = -pIfftOut[2*i+1];
	}
	delete[]x; 
}