📄 fft_ifft.c
字号:
case 480 : sinw= -0.013090; break;
case 2400 : sinw= -0.002618; break;
case 288 : sinw= -0.021815; break;
case 1440 : sinw= -0.004363; break;
case 864 : sinw= -0.007272; break;
case 2592 : sinw= -0.002424; break;
case 192 : sinw= -0.032719; break;
case 960 : sinw= -0.006545; break;
case 576 : sinw= -0.010908; break;
case 2880 : sinw= -0.002182; break;
case 1728 : sinw= -0.003636; break;
case 384 : sinw= -0.016362; break;
case 1920 : sinw= -0.003272; break;
case 1152 : sinw= -0.005454; break;
case 768 : sinw= -0.008181; break;
case 2304 : sinw= -0.002727; break;
case 1536 : sinw= -0.004091; break;
case 3072 : sinw= -0.002045; break;
case 8 : sinw= -0.707107; break;
case 32 : sinw= -0.195090; break;
case 128 : sinw= -0.049068; break;
case 4 : sinw= -1.000000; break;
case 16 : sinw= -0.382683; break;
case 64 : sinw= -0.098017; break;
case 256 : sinw= -0.024541; break;
case 512 : sinw= -0.012272; break;
case 1024 : sinw= -0.006136; break;
default : sinw= 0; break;
}
}
void twiddleTransf(int sofarRadix, int sofar_mul_Radix, int radix, int remainRadix, int fft_length,
double yRe[], double yIm[])
{ /* twiddleTransf */
double gem;
double t1_re,t1_im, t2_re,t2_im, t3_re,t3_im;
double t4_re,t4_im, t5_re,t5_im;
double m2_re,m2_im, m3_re,m3_im, m4_re,m4_im;
double m1_re,m1_im, m5_re,m5_im;
double s1_re,s1_im, s2_re,s2_im, s3_re,s3_im;
double s4_re,s4_im, s5_re,s5_im;
/*Table
cosw=cos(2*pi/(double)(sofar_mul_Radix))
sinw-sin(2*pi/(double)(sofar_mul_Radix)) */
scin_table (sofar_mul_Radix);
tw_re = 1.0;
tw_im = 0;
dataOffset=0;
groupOffset=dataOffset;
adr=groupOffset;
//sofarRadix means the group Numbers in particular stage
for (dataNo=0; dataNo<sofarRadix; dataNo++)
{
/* sofarRadix =1 means the first stage */
if (sofarRadix>1)
{
twiddleRe[0] = 1.0;
twiddleIm[0] = 0.0;
twiddleRe[1] = tw_re;
twiddleIm[1] = tw_im;
for (twNo=2; twNo<radix; twNo++)
{
twiddleRe[twNo]=tw_re*twiddleRe[twNo-1]
- tw_im*twiddleIm[twNo-1];
twiddleIm[twNo]=tw_im*twiddleRe[twNo-1]
+ tw_re*twiddleIm[twNo-1];
}
gem = cosw*tw_re - sinw*tw_im;
tw_im = sinw*tw_re + cosw*tw_im;
tw_re = gem;
}
//remainRadix means the number of butterfly-units in one group
for (groupNo=0; groupNo<remainRadix; groupNo++)
{
/* dataNo=0 means group number equals zero,
where the twiddle factors are useless */
if ((sofarRadix>1) && (dataNo > 0))
{
zRe[0]=yRe[adr];
zIm[0]=yIm[adr];
blockNo=1;
do {
adr = adr + sofarRadix;
zRe[blockNo]= twiddleRe[blockNo] * yRe[adr]
- twiddleIm[blockNo] * yIm[adr];
zIm[blockNo]= twiddleRe[blockNo] * yIm[adr]
+ twiddleIm[blockNo] * yRe[adr];
blockNo++;
} while (blockNo < radix);
}
else
for (blockNo=0; blockNo<radix; blockNo++)
{
zRe[blockNo]=yRe[adr];
zIm[blockNo]=yIm[adr];
adr=adr+sofarRadix;
}
/* details of different butterfly-units*/
switch(radix) {
case 2 : gem=zRe[0] + zRe[1];
zRe[1]=zRe[0] - zRe[1]; zRe[0]=gem;
gem=zIm[0] + zIm[1];
zIm[1]=zIm[0] - zIm[1]; zIm[0]=gem;
break;
case 3 : t1_re=zRe[1] + zRe[2]; t1_im=zIm[1] + zIm[2];
/* restore to zRe[0] to save space */
zRe[0]=zRe[0] + t1_re; zIm[0]=zIm[0] + t1_im;
/* cos(2*pi/3)-1 is needed */
m1_re=c3_1*t1_re; m1_im=c3_1*t1_im;
m2_re=c3_2*(zIm[1] - zIm[2]);
m2_im=c3_2*(zRe[2] - zRe[1]);
s1_re=zRe[0] + m1_re; s1_im=zIm[0] + m1_im;
zRe[1]=s1_re + m2_re; zIm[1]=s1_im + m2_im;
zRe[2]=s1_re - m2_re; zIm[2]=s1_im - m2_im;
break;
case 4 : t1_re=zRe[0] + zRe[2]; t1_im=zIm[0] + zIm[2];
t2_re=zRe[1] + zRe[3]; t2_im=zIm[1] + zIm[3];
m2_re=zRe[0] - zRe[2]; m2_im=zIm[0] - zIm[2];
m3_re=zIm[1] - zIm[3]; m3_im=zRe[3] - zRe[1];
zRe[0]=t1_re + t2_re; zIm[0]=t1_im + t2_im;
zRe[2]=t1_re - t2_re; zIm[2]=t1_im - t2_im;
zRe[1]=m2_re + m3_re; zIm[1]=m2_im + m3_im;
zRe[3]=m2_re - m3_re; zIm[3]=m2_im - m3_im;
break;
case 5 : t1_re=zRe[1] + zRe[4]; t1_im=zIm[1] + zIm[4];
t2_re=zRe[2] + zRe[3]; t2_im=zIm[2] + zIm[3];
t3_re=zRe[1] - zRe[4]; t3_im=zIm[1] - zIm[4];
t4_re=zRe[3] - zRe[2]; t4_im=zIm[3] - zIm[2];
t5_re=t1_re + t2_re; t5_im=t1_im + t2_im;
/* restore to zRe[0] to save space */
zRe[0]=zRe[0] + t5_re; zIm[0]=zIm[0] + t5_im;
/* (cos(u5)+cos(2*u5))/2-1 is needed */
m1_re=c5_1*t5_re; m1_im=c5_1*t5_im;
m2_re=c5_2*(t1_re - t2_re);
m2_im=c5_2*(t1_im - t2_im);
m3_re=-c5_3*(t3_im + t4_im);
m3_im=c5_3*(t3_re + t4_re);
m4_re=-c5_4*t4_im; m4_im=c5_4*t4_re;
m5_re=-c5_5*t3_im; m5_im=c5_5*t3_re;
s3_re=m3_re - m4_re; s3_im=m3_im - m4_im;
s5_re=m3_re + m5_re; s5_im=m3_im + m5_im;
s1_re=zRe[0] + m1_re; s1_im=zIm[0] + m1_im;
s2_re=s1_re + m2_re; s2_im=s1_im + m2_im;
s4_re=s1_re - m2_re; s4_im=s1_im - m2_im;
zRe[1]=s2_re + s3_re; zIm[1]=s2_im + s3_im;
zRe[2]=s4_re + s5_re; zIm[2]=s4_im + s5_im;
zRe[3]=s4_re - s5_re; zIm[3]=s4_im - s5_im;
zRe[4]=s2_re - s3_re; zIm[4]=s2_im - s3_im;
break;
default : break;
} /* butterfly-units */
/* process of data restore to the formal array : yRe[] and yIm[]*/
adr=groupOffset;
for (blockNo=0; blockNo<radix; blockNo++)
{
yRe[adr]=zRe[blockNo]; yIm[adr]=zIm[blockNo];
adr=adr+sofarRadix;
}
/* generate new address within the same group */
groupOffset=groupOffset+sofar_mul_Radix;//sofarRadix*radix;
adr=groupOffset;
}
/* generate new address between closed groups */
dataOffset=dataOffset+1;
groupOffset=dataOffset;
adr=groupOffset;
}
} /* twiddleTransf */
/*suitable for any kinds of length of fft composed by the radix: 2,3,4,5
*/
void fft(int n, double xRe[], double xIm[],
double yRe[], double yIm[])
{
int sofarRadix[maxFactorCount],
actualRadix[maxFactorCount],
remainRadix[maxFactorCount];
int nFactor;
int count;
int sofar_mul_Radix;
pi = 4*atan(1);
/* set up tables using for control the times of buffterfly-unit in each stage
*/
transTableSetup(sofarRadix, actualRadix, remainRadix, &nFactor, &n);
/*the sequence needs to be resorted before fft
*/
permute(n, nFactor, actualRadix, remainRadix, xRe, xIm, yRe, yIm);
/* Every fft process contains "nFactor" stages
*/
for (count=1; count<=nFactor; count++)
{
if (count != nFactor) sofar_mul_Radix = sofarRadix[count+1];
else sofar_mul_Radix = remainRadix[0];
twiddleTransf(sofarRadix[count], sofar_mul_Radix, actualRadix[count],
remainRadix[count], remainRadix[0], yRe, yIm);
}
} /* fft */
/*interface defined
*/
int fftio(long n,double are[],double aim[],double bre[],double bim[])
{
fft((int)n, are, aim, bre, bim);
return(0);
}
/*interface defined
*/
int ifftio(long n,double are[],double aim[],double bre[],double bim[])
{
int i;
for (i=0; i<n; i++)
aim[i]= - aim[i];
fft((int)n, are, aim, bre, bim);
for (i=0; i<n; i++)
bim[i]= - bim[i];
return(0);
}
/*test program with input of x[i]=cos(pi*(i*i)/(n*n))+sin(pi*(i*i)/(n*n)) i=0~fft_length-1
*/
int test(long fftLength, int fft_ifft,
double xre[], double xim[], double yre[], double yim[])
{
int i, k, n;
double w, pi;
n = fftLength;
fprintf(fp, "\nwhen fftLength=%d : \n", n);
pi = 4 * atan(1);
k = 1;
w = k*pi/n/n;
/* input signals*/
for (i=0; i<n; i++) {
xre[i] = cos(i*i*w);
xim[i] = sin(i*i*w);
}
/* one fft caculate process */
if (fft_ifft) fftio(n, xre, xim, yre, yim);
else ifftio(n, xre, xim, yre, yim);
/* output the result to a file in order to check the results
*/
k=0;
for (i=0; i<n; i++)
if (k<4) {
k++;
fprintf(fp,"%f + %fj, ", yre[i], yim[i]);
}else {
k=0;
fprintf(fp,"%f + %fj\n", yre[i], yim[i]);
}
/**/
return(0);
} /* fft_test */
/* FFT benchmark main program.
*/
void main()
{
double *are, *aim, *bre, *bim;
int x2, x3, x5;
//int fft_totolNumber;
int nfft;
int fft_ifft=1;
//int i;
are = (double *) calloc(maxIndex, sizeof(double));
aim = (double *) calloc(maxIndex, sizeof(double));
bre = (double *) calloc(maxIndex, sizeof(double));
bim = (double *) calloc(maxIndex, sizeof(double));
if( bim == NULL )
{
printf( "Insufficient memory available\n" );
exit(2);
}
printf("\nFFT Benchmark Test");
printf("input fft or ifft: \n ");
printf("FFT : 1\n");
printf("IFFT : 0\n");
printf("Please choose: ");
scanf("%d", &fft_ifft);
fp = fopen("fft_result.txt", "w");
// fft_totolNumber = 0;
/* simulating all possible probability*/
for (x2=2; x2<=10; x2++)
for (x3=1; x3<=6; x3++)
for (x5=0; x5<=3; x5++)
{
nfft= (int)pow(2,x2)* (int)pow(3,x3)* (int)pow(5,x5);
if (pow(2,x2)* pow(3,x3)*pow(5,x5)<3120)
{
//fft_totolNumber++;
fprintf(fp, "\n2(%d)*3(%d)*5(%d)=%d: \n", x2,x3,x5, nfft);
test( nfft, fft_ifft, are, aim, bre, bim);
}
}
test( 128, fft_ifft, are, aim, bre, bim);
fprintf(fp, "\n\n");
test( 256, fft_ifft, are, aim, bre, bim);
fprintf(fp, "\n\n");
test( 512, fft_ifft, are, aim, bre, bim);
fprintf(fp, "\n\n");
test( 1024, fft_ifft, are, aim, bre, bim);
fprintf(fp,"\n");
fclose(fp);
free(bim);
free(bre);
free(aim);
free(are);
printf("\n\nFFT test terminated.\n");
} /* main */⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -