FFT.rar_dsp_fft_fft

/************************************************************ Copyright (C), 2007 by SEED Electronic Technology LTD. FileName: FFT_function.c Author: Marine Version : V1.0 Date:Dec25,2006 Description: FFT函数 *************************************************************/ #include "i_cmplx.h" /* definition of the complex type */ #include "twiddle1024.h" /* quantised and scaled Twiddle factors */ #define LL 1024 /* Maximum length of FFT */ /*********************************************************** Function: fft1024 Description: 1024点的FFT变换函数 Calls: No Called By: main Input: *Y：输入采样点数组 N：FFT变换的点数 Output: No Return: No Others: No ************************************************************/ void fft1024(COMPLEX *Y, int N) /* FFT(input sample array, # of points) */ { int temp1R, temp1I, temp2R,temp2I; /* 32 bits temporary storage for */ /* intermediate results */ short tempR, tempI, c, s; /* 16 bits temporary storages */ /* variables */ int TwFStep, /* Step between twiddle factors */ TwFIndex, /* Index of twiddle factors */ BLStep, /* Step for incrementing butterfly index */ BLdiff, /* Difference between upper and lower butterfly legs */ upperIdx, lowerIdx, /* upper and lower indexes of buterfly leg */ i, j, k; /* loop control variables */ BLdiff=N; TwFStep=1; for(k=N;k>1;k=(k>>1)) /* Do Log(base 2)(N) Stages */ { BLStep=BLdiff; BLdiff=BLdiff>>1; TwFIndex=0; for(j=0;j<BLdiff;j++)/* Nbr of twiddle factors to use=BLDiff */ { c=w[TwFIndex].real; s=w[TwFIndex].imag; TwFIndex=TwFIndex+TwFStep; /* Now do N/BLStep butterflies */ for(upperIdx=j; upperIdx<N; upperIdx+=BLStep) { /* Calculations inside this loop avoid overflow by shifting left once the result of every adittion/substration and by shifting left 15 places the result of every multiplication. Double precision temporary results (32-bit) are used in order to avoid losing information because of overflow. Final DFT result is scaled by N (number of points), i.e., 2^(Nbr of stages) =2^(log(base 2) N) = N */ lowerIdx=upperIdx+BLdiff; temp1R = (Y[upperIdx].real - Y[lowerIdx].real)>>1; temp2R = (Y[upperIdx].real + Y[lowerIdx].real)>>1; Y[upperIdx].real = (short) temp2R; temp1I = (Y[upperIdx].imag - Y[lowerIdx].imag)>>1; temp2I = (Y[upperIdx].imag + Y[lowerIdx].imag)>>1; Y[upperIdx].imag = (short) temp2I; temp2R = (c*temp1R - s*temp1I)>>15; Y[lowerIdx].real = (short) temp2R; temp2I = (c*temp1I + s*temp1R)>>15; Y[lowerIdx].imag = (short) temp2I; } } TwFStep = TwFStep<<1; /* update separation of twiddle factors)*/ } /* bit reversal for resequencing data */ j=0; for (i=1;i<(N-1);i++) { k=N/2; while (k<=j) { j = j-k; k = k/2; } j=j+k; if ( i < j ) { tempR = Y[j].real; tempI = Y[j].imag; Y[j].real = Y[i].real; Y[j].imag = Y[i].imag; Y[i].real = tempR; Y[i].imag = tempI; } } return; } /*********************************************************** Function: fft512 Description: 512点的FFT变换函数 Calls: No Called By: main Input: *Y：输入采样点数组 N：FFT变换的点数 Output: No Return: No Others: No ************************************************************/ void fft512(COMPLEX *Y, int N) /* FFT(input sample array, # of points) */ { int temp1R, temp1I, temp2R,temp2I; /* 32 bits temporary storage for */ /* intermediate results */ short tempR, tempI, c, s; /* 16 bits temporary storages */ /* variables */ int TwFStep, /* Step between twiddle factors */ TwFIndex, /* Index of twiddle factors */ BLStep, /* Step for incrementing butterfly index */ BLdiff, /* Difference between upper and lower butterfly legs */ upperIdx, lowerIdx, /* upper and lower indexes of buterfly leg */ i, j, k; /* loop control variables */ BLdiff=N; TwFStep=1; for(k=N; k>1; k=(k>>1)) /* Do Log(base 2)(N) Stages */ { BLStep=BLdiff; BLdiff=BLdiff>>1; TwFIndex=0; for(j=0; j<BLdiff; j++)/* Nbr of twiddle factors to use=BLDiff */ { c=w512[TwFIndex].real; s=w512[TwFIndex].imag; TwFIndex=TwFIndex + TwFStep; /* Now do N/BLStep butterflies */ for(upperIdx=j;upperIdx<N;upperIdx+=BLStep) { lowerIdx=upperIdx+BLdiff; temp1R = (Y[upperIdx].real - Y[lowerIdx].real)>>1; temp2R = (Y[upperIdx].real + Y[lowerIdx].real)>>1; Y[upperIdx].real = (short) temp2R; temp1I = (Y[upperIdx].imag - Y[lowerIdx].imag)>>1; temp2I = (Y[upperIdx].imag + Y[lowerIdx].imag)>>1; Y[upperIdx].imag = (short) temp2I; temp2R = (c*temp1R - s*temp1I)>>15; Y[lowerIdx].real = (short) temp2R; temp2I = (c*temp1I + s*temp1R)>>15; Y[lowerIdx].imag = (short) temp2I; } } TwFStep = TwFStep<<1; /* update separation of twiddle factors)*/ } /* bit reversal for resequencing data */ j=0; for (i=1;i<(N-1);i++) { k=N/2; while (k<=j) { j = j-k; k = k/2; } j=j+k; if ( i < j ) { tempR = Y[j].real; tempI = Y[j].imag; Y[j].real = Y[i].real; Y[j].imag = Y[i].imag; Y[i].real = tempR; Y[i].imag = tempI; } } return; } /*********************************************************** Function: fft256 Description: 256点的FFT变换函数 Calls: No Called By: main Input: *Y：输入采样点数组 N：FFT变换的点数 Output: No Return: No Others: No ************************************************************/ void fft256(COMPLEX *Y, int N) /* FFT(input sample array, # of points) */ { int temp1R, temp1I, temp2R,temp2I; /* 32 bits temporary storage for */ /* intermediate results */ short tempR, tempI, c, s; /* 16 bits temporary storages */ /* variables */ int TwFStep, /* Step between twiddle factors */ TwFIndex, /* Index of twiddle factors */ BLStep, /* Step for incrementing butterfly index */ BLdiff, /* Difference between upper and lower butterfly legs */ upperIdx, lowerIdx, /* upper and lower indexes of buterfly leg */ i, j, k; /* loop control variables */ BLdiff=N; TwFStep=1; for(k=N;k>1;k=(k>>1)) /* Do Log(base 2)(N) Stages */ { BLStep=BLdiff; BLdiff=BLdiff>>1; TwFIndex=0; for(j=0;j<BLdiff;j++)/* Nbr of twiddle factors to use=BLDiff */ { c=w256[TwFIndex].real; s=w256[TwFIndex].imag; TwFIndex=TwFIndex+TwFStep; /* Now do N/BLStep butterflies */ for(upperIdx=j; upperIdx<N; upperIdx+=BLStep) {