DSP

1024点的FFT算法

2019-07-13 19:49发布

站内文章 / DSP
13047 0 1082
牛人写的1024点FFT，结果直接输出到XIn里，需要最少的内存空间 void FFT1024(float *XIn,float *YIn) { //============================================= // XIn 采样数据 // YIn 虚部全为0 //============================================= float TWOPI = 6.283185307179586; int N,M,N2,iN,ID,I0,I1,I2,I3,J,K,N1; int N4; int I; float R1,R2,S1,S2,S3; float X1,E,A; float A3,CC1,SS1,CC3,SS3; float XT; N = 1024; M = 10; N2 = 2 * N; for(K = 1;K < M;K ++) { N2 = N2 / 2; N4 = N2 / 4; X1 = N2; E = TWOPI / X1; A = 0; for(J = 1;J <= N4;J ++) { A3 = 3 * A; CC1 = cos(A); SS1 = sin(A); CC3 = cos(A3); SS3 = sin(A3); X1 = J; A = X1 * E; iN = J; ID = 2 * N2; do{ for(I0 = iN;I0 < N;I0 += ID) { I1 = I0 + N4; I2 = I1 + N4; I3 = I2 + N4; R1 = XIn[I0 - 1] - XIn[I2 - 1]; XIn[I0 - 1] = XIn[I0 - 1] + XIn[I2 - 1]; R2 = XIn[I1 - 1] - XIn[I3 - 1]; XIn[I1 - 1] = XIn[I1 - 1] + XIn[I3 - 1]; S1 = YIn[I0 - 1] - YIn[I2 - 1]; YIn[I0 - 1] = YIn[I0 - 1] + YIn[I2 - 1]; S2 = YIn[I1 - 1] - YIn[I3 - 1]; YIn[I1 - 1] = YIn[I1 - 1] + YIn[I3 - 1]; S3 = R1 - S2; R1 = R1 + S2; S2 = R2 - S1; R2 = R2 + S1; XIn[I2 - 1] = R1 * CC1 - S2 * SS1; YIn[I2 - 1] = -S2 * CC1 - R1 * SS1; XIn[I3 - 1] = S3 * CC3 + R2 * SS3; YIn[I3 - 1] = R2 * CC3 - S3 * SS3; } iN = 2 * ID - N2 + J; ID = 4 * ID; }while(iN < N); } } //-------LAST STAGE, LENGTH-2 BUTTERFLY------ iN = 1; ID = 4; do{ for(I0 = iN;I0 <= N;I0 += ID) { I1 = I0 + 1; R1 = XIn[I0]; XIn[I0 - 1] = R1 + XIn[I1 - 1]; XIn[I1 - 1] = R1 - XIn[I1 - 1]; R1 = YIn[I0 - 1]; YIn[I0 - 1] = R1 + YIn[I1 - 1]; YIn[I1 - 1] = R1 - YIn[I1 - 1]; } iN = 2 * ID - 1; ID = 4 * ID; }while(iN < N); //------BIT REVERSE COUNTER------ J = 1; N1 = N - 1; for(I = 1;I <= N1;I ++) { if(I < J) { XT = XIn[J - 1]; XIn[J - 1] = XIn[I - 1]; XIn[I - 1] = XT; XT = YIn[J - 1]; YIn[J - 1] = YIn[I - 1]; YIn[I - 1] = XT; } K = N / 2; while(K < J) { J = J - K; K = K / 2; } J = J + K; } }

热门文章