📄 fft.c
字号:
#include <math.h> */
*/
#ifndef PI */
# define PI (3.14159265358979323846) */
#endif */
*/
short d2s(double d) */
{ */
d = floor(0.5 + d); // Explicit rounding to integer // */
if (d >= 32767.0) return 32767; */
if (d <= -32768.0) return -32768; */
return (short)d; */
} */
*/
void gen_twiddle(short *w, int n) */
{ */
double M = 32767.5; */
int i, j, k; */
*/
for (j = 1, k = 0; j < n >> 2; j = j << 2) */
{ */
for (i = 0; i < n >> 2; i += j << 1) */
{ */
w[k + 11] = d2s(M * cos(6.0 * PI * (i + j) / n)); */
w[k + 10] = d2s(M * sin(6.0 * PI * (i + j) / n)); */
w[k + 9] = d2s(M * cos(6.0 * PI * (i ) / n)); */
w[k + 8] = d2s(M * sin(6.0 * PI * (i ) / n)); */
*/
w[k + 7] = d2s(M * cos(4.0 * PI * (i + j) / n)); */
w[k + 6] = d2s(M * sin(4.0 * PI * (i + j) / n)); */
w[k + 5] = d2s(M * cos(4.0 * PI * (i ) / n)); */
w[k + 4] = d2s(M * sin(4.0 * PI * (i ) / n)); */
*/
w[k + 3] = d2s(M * cos(2.0 * PI * (i + j) / n)); */
w[k + 2] = d2s(M * sin(2.0 * PI * (i + j) / n)); */
w[k + 1] = d2s(M * cos(2.0 * PI * (i ) / n)); */
w[k + 0] = d2s(M * sin(2.0 * PI * (i ) / n)); */
*/
k += 12; */
} */
} */
w[2*n - 1] = w[2*n - 3] = w[2*n - 5] = 32767; */
w[2*n - 2] = w[2*n - 4] = w[2*n - 6] = 0; */
} */
/* */
/* ASSUMPTIONS */
/* n must be a power of 4 and n >= 16 & n < 32768. */
/* FFT data x are aligned on a double word boundary, in real/imag */
/* pairs, FFT twiddle factors w are also aligned on a double word */
/* boundary in real/imaginary pairs. */
/* */
/* Input FFT coeffs. are in signed Q.15 format. */
/* The memory Configuration is LITTLE ENDIAN. */
/* The complex data will be returned in natural order. This code is */
/* uninteruptable, interupts are disabled on entry to the function and */
/* re-enabled on exit. */
/* */
/* MEMORY NOTE */
/* No bank conflict stalls occur in this code. */
/* */
/* TECHNIQUES */
/* A special sequence of coefficients. are used (as generated above) */
/* to produce the DSP_fft. This collapses the inner 2 loops in the */
/* taditional Burrus and Parks implementation Fortran Code. */
/* */
/* The following C code represents an implementation of the Cooley */
/* Tukey radix 4 DIF FFT. It accepts the inputs in normal order and */
/* produces the outputs in digit reversed order. The natural C code */
/* shown in this file on the other hand, accepts the inputs in nor- */
/* mal order and produces the outputs in normal order. */
/* */
/* Several transformations have been applied to the original Cooley */
/* Tukey code to produce the natural C code description shown here. */
/* In order to understand these it would first be educational to */
/* understand some of the issues involved in the conventional Cooley */
/* Tukey FFT code. */
/* */
void radix4(int n, short x[], short wn[]) */
{ */
int n1, n2, ie, ia1, ia2, ia3; */
int i0, i1, i2, i3, i, j, k; */
short co1, co2, co3, si1, si2, si3; */
short xt0, yt0, xt1, yt1, xt2, yt2; */
short xh0, xh1, xh20, xh21, xl0, xl1,xl20,xl21; */
*/
n2 = n; */
ie = 1; */
for (k = n; k > 1; k >>= 2) */
{ */
n1 = n2; */
n2 >>= 2; */
ia1 = 0; */
*/
for (j = 0; j < n2; j++) */
{ */
ia2 = ia1 + ia1; */
ia3 = ia2 + ia1; */
*/
co1 = wn[2 * ia1 ]; */
si1 = wn[2 * ia1 + 1]; */
co2 = wn[2 * ia2 ]; */
si2 = wn[2 * ia2 + 1]; */
co3 = wn[2 * ia3 ]; */
si3 = wn[2 * ia3 + 1]; */
ia1 = ia1 + ie; */
*/
for (i0 = j; i0< n; i0 += n1) */
{ */
i1 = i0 + n2; */
i2 = i1 + n2; */
i3 = i2 + n2; */
*/
*/
xh0 = x[2 * i0 ] + x[2 * i2 ]; */
xh1 = x[2 * i0 + 1] + x[2 * i2 + 1]; */
xl0 = x[2 * i0 ] - x[2 * i2 ]; */
xl1 = x[2 * i0 + 1] - x[2 * i2 + 1]; */
*/
xh20 = x[2 * i1 ] + x[2 * i3 ]; */
xh21 = x[2 * i1 + 1] + x[2 * i3 + 1]; */
xl20 = x[2 * i1 ] - x[2 * i3 ]; */
xl21 = x[2 * i1 + 1] - x[2 * i3 + 1]; */
*/
x[2 * i0 ] = xh0 + xh20; */
x[2 * i0 + 1] = xh1 + xh21; */
*/
xt0 = xh0 - xh20; */
yt0 = xh1 - xh21; */
xt1 = xl0 + xl21; */
yt2 = xl1 + xl20; */
xt2 = xl0 - xl21; */
yt1 = xl1 - xl20; */
*/
x[2 * i1 ] = (xt1 * co1 + yt1 * si1) >> 15; */
x[2 * i1 + 1] = (yt1 * co1 - xt1 * si1) >> 15; */
x[2 * i2 ] = (xt0 * co2 + yt0 * si2) >> 15; */
x[2 * i2 + 1] = (yt0 * co2 - xt0 * si2) >> 15; */
x[2 * i3 ] = (xt2 * co3 + yt2 * si3) >> 15; */
x[2 * i3 + 1] = (yt2 * co3 - xt2 * si3) >> 15; */
} */
} */
*/
ie <<= 2; */
} */
} */
/* */
/* The conventional Cooley Tukey FFT, is written using three loops. */
/* The outermost loop "k" cycles through the stages. There are log */
/* N to the base 4 stages in all. The loop "j" cycles through the */
/* groups of butterflies with different twiddle factors, loop "i" */
/* reuses the twiddle factors for the different butterflies within */
/* a stage. It is interesting to note the following: */
/* */
/*-------------------------------------------------------------------------S*/
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -