wib_fft_f.c

802.16e物理层FFT算法的实现。严格按照802.16e实现。

            fft_tbl fft_tbl_inst
            )
{
    //complex16 data_reorder[256];
    unsigned short i,j,k,p;
    

    signed short exp_statge[3] = {0,0,0},exponent = 0;

    unsigned char bit_width;   //width of fft compute
    unsigned short N;                //point of fft
            

    unsigned char COEF_W ;
    unsigned short BOUND1;
    unsigned short BOUND2,BOUND3;
    unsigned char log2n;
   
   
    
    complex16 *p_dat;
    
    unsigned short base,index,index_wp,base_index_wp;
    complex16 x1,x2,x3,x4;

    
    unsigned short temp1;

    complex16 *p_w1p,*p_w2p,*p_w3p;
    unsigned short *p_reorder_idx;
    unsigned char degree;
	FILE *fp_outdbg;
//	signed short test;
    
    p_w1p = fft_tbl_inst.p_w1p;
    p_w2p = fft_tbl_inst.p_w2p;
    p_w3p = fft_tbl_inst.p_w3p;
    p_reorder_idx = fft_tbl_inst.p_reorder_idx;
    p_dat = p_dat_out;
    N = fft_tbl_inst.N;
    bit_width = fft_tbl_inst.bit_width;

    BOUND1 = (unsigned short)(1 << (bit_width-1));
    BOUND2 = (unsigned short)(1 << (bit_width));
	BOUND3 = (unsigned short)(1 << (bit_width+1));
    COEF_W = bit_width-1;
    degree = fft_tbl_inst.degree;
    log2n = fft_tbl_inst.log2n;
    
   
    //data reorder
    for (i = 0; i < N; i++)
    {
        p_dat[i].r = (p_dat_in + *(p_reorder_idx+i)) -> r;
        if (mode) {
            p_dat[i].i = (signed short)(-((p_dat_in + *(p_reorder_idx+i)) -> i));
        } else {
            p_dat[i].i = (p_dat_in + *(p_reorder_idx+i)) -> i;
        }

	}
	for (i = 0; i < N; i++) {
		if (((p_dat[i].r) >= BOUND1) || ((p_dat[i].r) < -BOUND1)){
            exp_statge[0] = 1;
		}
		if (((p_dat[i].i) >= BOUND1) || ((p_dat[i].i) < -BOUND1)){
            exp_statge[0] = 1;
		}
	}

	if(exp_statge[0]) {
		for (i = 0; i < N; i++) {
			p_dat[i].r = p_dat[i].r >> 1;
			p_dat[i].i = p_dat[i].i >> 1;
		}
	}
   

 //printf("N:%d\n",N >> 6);  
   
	exp_statge[0] = 0;
    k = 0;
    base_index_wp = 0;
    for (i = 0; i < degree; i++) {
        base_index_wp = (unsigned short)(base_index_wp + k);
        index_wp = base_index_wp;
        //t2  = N;
        temp1 = (unsigned short)(N >> (unsigned short)((i + 1) * 2));
        //temp1 = (N / pow(4,i+1));
        // printf("tem:%d\n%d\n",temp1,i);  

        for (j = 0; j < temp1; j++) {      //for k=1:N/(4^degree) 
            base = (unsigned short)(j * POW4I);
            
            k = 0;
            do {
                index = (unsigned short)(base + k);
                //if result is too large then shorten bits
				if (exp_statge[2]) {
					x1.r = (signed short)(p_dat[index].r >> 3);
                    x1.i = (signed short)(p_dat[index].i >> 3);
                    x2.r = (signed short)(p_dat[index+POW4I_1].r >> 3);
                    x2.i = (signed short)(p_dat[index+POW4I_1].i >> 3);
                    x3.r = (signed short)(p_dat[index+2*POW4I_1].r >> 3);
                    x3.i = (signed short)(p_dat[index+2*POW4I_1].i >> 3);
                    x4.r = (signed short)(p_dat[index+3*POW4I_1].r >> 3);
                    x4.i = (signed short)(p_dat[index+3*POW4I_1].i >> 3);
				}
                else if (exp_statge[1]) {
                    x1.r = (signed short)(p_dat[index].r >> 2);
                    x1.i = (signed short)(p_dat[index].i >> 2);
                    x2.r = (signed short)(p_dat[index+POW4I_1].r >> 2);
                    x2.i = (signed short)(p_dat[index+POW4I_1].i >> 2);
                    x3.r = (signed short)(p_dat[index+2*POW4I_1].r >> 2);
                    x3.i = (signed short)(p_dat[index+2*POW4I_1].i >> 2);
                    x4.r = (signed short)(p_dat[index+3*POW4I_1].r >> 2);
                    x4.i = (signed short)(p_dat[index+3*POW4I_1].i >> 2);
                } else if (exp_statge[0]) {
                    x1.r = (signed short)(p_dat[index].r >> 1);
                    x1.i = (signed short)(p_dat[index].i >> 1);
                    x2.r = (signed short)(p_dat[index+POW4I_1].r >> 1);
                    x2.i = (signed short)(p_dat[index+POW4I_1].i >> 1);
                    x3.r = (signed short)(p_dat[index+2*POW4I_1].r >> 1);
                    x3.i = (signed short)(p_dat[index+2*POW4I_1].i >> 1);
                    x4.r = (signed short)(p_dat[index+3*POW4I_1].r >> 1);
                    x4.i = (signed short)(p_dat[index+3*POW4I_1].i >> 1);
                } else {
                    x1.r = p_dat[index].r;
                    x1.i = p_dat[index].i;
                    x2.r = p_dat[index+POW4I_1].r;
                    x2.i = p_dat[index+POW4I_1].i;
                    x3.r = p_dat[index+2*POW4I_1].r;
                    x3.i = p_dat[index+2*POW4I_1].i;
                    x4.r = p_dat[index+3*POW4I_1].r;
                    x4.i = p_dat[index+3*POW4I_1].i;
                }
                index_wp = (unsigned short)(base_index_wp + k);
                // base4 butterfly              
                dft4p(x1,x2,x3,x4,p_w1p[index_wp],p_w2p[index_wp],p_w3p[index_wp],
                    COEF_W,&p_dat[index],&p_dat[index+POW4I_1],
                    &p_dat[index+2*POW4I_1],&p_dat[index+3*POW4I_1]);
                
                
             }while(++k < POW4I_1);
        } 
        
       	if (exp_statge[2]) {
			exponent = (signed short)(exponent + 3);
		} else if (exp_statge[1]) {
			exponent = (signed short)(exponent + 2);
		} else if(exp_statge[0]) {
			exponent = (signed short)(exponent + 1);
		}

        exp_statge[0] = 0;
        exp_statge[1] = 0;
		exp_statge[2] = 0;
		
        for (p = 0; p < N; p++) {
			if (((p_dat[p].r) >= BOUND3) || ((p_dat[p].r) < -BOUND3)){
                exp_statge[2] = 1;
            } else if (((p_dat[p].r) >= BOUND2) || ((p_dat[p].r) < -BOUND2)){
                exp_statge[1] = 1;
            } else if (((p_dat[p].r) >= BOUND1) || ((p_dat[p].r) < -BOUND1)){
                exp_statge[0] = 1;
            }

			if (((p_dat[p].i) >= BOUND3) || ((p_dat[p].i) < -BOUND3)){
                exp_statge[2] = 1;
            } else if (((p_dat[p].i) >= BOUND2) || ((p_dat[p].i) < -BOUND2)){
                exp_statge[1] = 1;
            } else if (((p_dat[p].i) >= BOUND1) || ((p_dat[p].i) < -BOUND1)) {
                exp_statge[0] = 1;
            }
			
        }

    }
   // }

   // xx = clock()-xx;
   // printf("FFT time is %d",xx);
    exp_statge[0] = 0;
    exp_statge[1] = 0;
	exp_statge[2] = 0;
		
    for (p = 0; p < N; p++) {
		if (((p_dat[p].r) >= BOUND3) || ((p_dat[p].r) < -BOUND3)){
            exp_statge[2] = 1;
        } else if (((p_dat[p].r) >= BOUND2) || ((p_dat[p].r) < -BOUND2)){
            exp_statge[1] = 1;
        } else if (((p_dat[p].r) >= BOUND1) || ((p_dat[p].r) < -BOUND1)){
            exp_statge[0] = 1;
        }

		if (((p_dat[p].i) >= BOUND3) || ((p_dat[p].i) < -BOUND3)){
            exp_statge[2] = 1;
        } else if (((p_dat[p].i) >= BOUND2) || ((p_dat[p].i) < -BOUND2)){
            exp_statge[1] = 1;
        } else if (((p_dat[p].i) >= BOUND1) || ((p_dat[p].i) < -BOUND1)) {
            exp_statge[0] = 1;
        }
			
    }

	if (exp_statge[2]) {
        for (i = 0; i < N; i++) {
            (p_dat+i)->r = (signed short)(((p_dat+i)->r) >> 3);
            (p_dat+i)->i = (signed short)(((p_dat+i)->i) >> 3);
        }
    } else if (exp_statge[1]) {
        for (i = 0; i < N; i++) {
            (p_dat+i)->r = (signed short)(((p_dat+i)->r) >> 2);
            (p_dat+i)->i = (signed short)(((p_dat+i)->i) >> 2);
        }
    } else if (exp_statge[0]) {
        for (i = 0; i < N; i++) {
            (p_dat+i)->r = (signed short)(((p_dat+i)->r) >> 1);
            (p_dat+i)->i = (signed short)(((p_dat+i)->i) >> 1);
        }
    }

	if (exp_statge[2]) {
        exponent = (signed short)(exponent + 3);
    } else if (exp_statge[1]) {
        exponent = (signed short)(exponent + 2);
    } else if(exp_statge[0]) {
        exponent = (signed short)(exponent + 1);
    }

    if (mode) {
        for (i = 0; i < N; i++) {
			if(((p_dat+i)->i) == -BOUND1) {
				(p_dat+i)->i = BOUND1-1;
			} else {
				(p_dat+i)->i = (signed short)(-((p_dat+i)->i));
			}
        }
    }

    if (mode) {
        exponent = (signed short)(exponent-log2n);
    }

	if(exponent==-8) {
		exponent=-8;
	}

	#ifdef  DEBUG_OUT_FFT_EXPONENT
	if (mode) {
		fp_outdbg   = fopen("fft_exponent.txt","at");   //append, txt
		fprintf(fp_outdbg,"%d\n",exponent);
		fclose(fp_outdbg);
	}
    #endif
    
    return exponent;
}