关于FFT定点算法

prince2010

要求：采集220V、50HZ交流电，32个采样点，计算直流分量、基波及各次谐波，速度要求较高，所以要求定点。

附采集数据样本(AD原始数据是12位的,这里已经做了转换,以下数据就是交流电压瞬时值)——

{0,43,109,175,224,262,288,309,315,316,293,262,225,172,117,59,6,-46,-120,-182,-232,-269,-294,-305,-309,-306,-289,-255,-216,-162,-105,-47};

目前移植了一个浮点运算的例程,但速度太慢,约1ms.

本人以前没有玩过FFT,定点FFT在论坛上貌似暂时搜不到相关帖子,所以特发悬赏帖于此,欢迎大家热烈讨论.

目前的例程——

1. 
   
2. 
   
3. #include "math.h"
   
4. #define PI 3.1416              
   
5. #define FFT_N 32
   
6. 
   
7. struct compx {float real,imag;};               
   
8. struct  compx  s[FFT_N];
   
9. float R[FFT_N]=       
   
10. {
    
11. 0,
    
12. 43,
    
13. 109,
    
14. 175,
    
15. 224,
    
16. 262,
    
17. 288,
    
18. 309,
    
19. 315,
    
20. 316,
    
21. 293,
    
22. 262,
    
23. 225,
    
24. 172,
    
25. 117,
    
26. 59,
    
27. 6,
    
28. -46,
    
29. -120,
    
30. -182,
    
31. -232,
    
32. -269,
    
33. -294,
    
34. -305,
    
35. -309,
    
36. -306,
    
37. -289,
    
38. -255,
    
39. -216,
    
40. -162,
    
41. -105,
    
42. -47
    
43. };
    
44. 
    
45. float        cos_tab[17]={0.0};                  
    
46. float        sin_tab[17]={0.0};
    
47. 
    
48. struct compx EE(struct compx a,struct compx b)      
    
49. {
    
50.         struct compx c;
    
51.         c.real=a.real*b.real-a.imag*b.imag;
    
52.         c.imag=a.real*b.imag+a.imag*b.real;
    
53.         return(c);
    
54. }
    
55. 
    
56. 
    
57. void fft(struct compx *s)
    
58. {
    
59.         int f,m,nv2,nm1,i,k,l,j=0;
    
60.         struct compx u,w,t;
    
61. 
    
62.         nv2=FFT_N/2;                  
    
63.         nm1=FFT_N-1;  
    
64.         for(i=0;i<nm1;i++)        
    
65.         {
    
66.                 if(i<j)                     
    
67.                 {
    
68.                         t=s[j];           
    
69.                         s[j]=s[i];
    
70.                         s[i]=t;
    
71.                 }
    
72.                 k=nv2;                  
    
73.                 while(k<=j)         
    
74.                 {           
    
75.                         j=j-k;            
    
76.                         k=k/2;               
    
77.                 }
    
78.                 j=j+k;                 
    
79.         }
    
80. 
    
81.         {
    
82.                 int le,lei,ip;                             
    
83.                 for(m=1;m<=5;m++)                        
    
84.                 {                                         
    
85.                         le=2<<(m-1);                           
    
86.                         lei=le/2;                              
    
87.                         u.real=1.0;                             
    
88.                         u.imag=0.0;
    
89. //                        w.real=cos(PI/lei);      
    
90. //                        w.imag=-sin(PI/lei);
    
91.                         w.real=cos_tab[lei];                    
    
92.                         w.imag=-sin_tab[lei];
    
93.                         for(j=0;j<=lei-1;j++)                  
    
94.                         {
    
95.                                 for(i=j;i<=FFT_N-1;i=i+le)            
    
96.                                 {
    
97.                                         ip=i+lei;                           
    
98.                                         t=EE(s[ip],u);                    
    
99.                                         s[ip].real=s[i].real-t.real;
    
100.                                         s[ip].imag=s[i].imag-t.imag;
     
101.                                         s[i].real=s[i].real+t.real;
     
102.                                         s[i].imag=s[i].imag+t.imag;
     
103.                                 }
     
104.                                 u=EE(u,w);                           
     
105.                         }
     
106.                 }
     
107.         }
     
108. }
     
109. 
     
110. 
     
111. int main(void)   
     
112. {  
     
113. 
     
114.         uint32_t i,k;
     
115. 
     
116.         USART2PC_Config();
     
117. 
     
118. 
     
119.         //************Table**************
     
120.         for(i=1;i<=16;i++)
     
121.         {
     
122.                 cos_tab[i]=cos(PI/i);
     
123.                 sin_tab[i]=sin(PI/i);
     
124.         }
     
125.        
     
126. 
     
127.         //************Initialize s[]**************
     
128.         for(i=0;i<FFT_N;i++)
     
129.         {
     
130.                 s[i].real=R[i];                                
     
131.                 s[i].imag=0;                                
     
132.         }
     
133. 
     
134. 
     
135.         fft(s);
     
136.        
     
137. 
     
138.         for(i=0;i<FFT_N;i++)                          
     
139.         {
     
140.                 s[i].real=sqrtf(s[i].real*s[i].real+s[i].imag*s[i].imag);
     
141.         }
     
142. 
     
143.         while(1)
     
144.         {
     
145.                 printf("\r\n%f\r\n",s[0].real/FFT_N);
     
146.                 printf("\r\n%f\r\n",s[1].real/(FFT_N/2));
     
147.                 printf("\r\n%f\r\n",s[2].real/(FFT_N/2));
     
148.                 printf("\r\n%f\r\n",s[3].real/(FFT_N/2));
     
149.                 printf("\r\n%f\r\n",s[4].real/(FFT_N/2));
     
150.                 printf("\r\n%f\r\n",s[5].real/(FFT_N/2));
     
151.                 printf("\r\n%f\r\n",s[6].real/(FFT_N/2));
     
152.                 printf("\r\n%f\r\n",s[7].real/(FFT_N/2));
     
153.                 printf("\r\n%f\r\n",s[8].real/(FFT_N/2));
     
154.                 printf("\r\n%f\r\n",s[9].real/(FFT_N/2));
     
155.                 printf("\r\n%f\r\n",s[10].real/(FFT_N/2));
     
156.                 printf("\r\n%f\r\n",s[11].real/(FFT_N/2));
     
157.                 printf("\r\n%f\r\n",s[12].real/(FFT_N/2));
     
158.                 printf("\r\n%f\r\n",s[13].real/(FFT_N/2));
     
159.                
     
160.                 printf("\r\n*****************\r\n");
     
161.                
     
162.         }
     
163. }

复制代码

prince2010

呵呵,是不是发错地方了,我用得是STM32,不是DSP.

foxpro2005

请使用官方库...

prince2010

foxpro2005 发表于 2014-1-15 12:32
请使用官方库...

我的芯片是F2的，不知道可以不可以用？

gaolf_2012

用库嘛，arm_cfft_radix4_q31

笑笑我笑了

32个点不如用DFT。

prince2010

笑笑我笑了 发表于 2014-1-15 14:16
32个点不如用DFT。

要求速度很快.

moen

1里面都是有符号的浮点数在运算，肯定慢，要改成定点的。
2乘法要确认是启动了CPU的硬件乘法器。如果编译后的汇编没有启用乘法器要手动地去把数据放到相关寄存器中。
3除法尽量用移位
4其它运算要用查表的方式。比如开跟号。
5实际物理量是实数，你的算法中把虚部填成0参与运算很多运算是冗余的，虚部为0的情况结果是对称的，利用这个原理可以优化。
6就32点，没必要玩那么多花样，相输入倒位序你就用代码一行一行地写出来，相互之间赋值，没必要搞个循环计算谁排哪里。

senjet

以下是用过的16bit定点FFT，带防溢出的，看看是否对你有用。
/* fix_fft.c - Fixed-point in-place Fast Fourier Transform  */
/*
  All data are fixed-point short integers, in which -32768
  to +32768 represent -1.0 to +1.0 respectively. Integer
  arithmetic is used for speed, instead of the more natural
  floating-point.

  For the forward FFT (time -> freq), fixed scaling is
  performed to prevent arithmetic overflow, and to map a 0dB
  sine/cosine wave (i.e. amplitude = 32767) to two -6dB freq
  coefficients. The return value is always 0.

  For the inverse FFT (freq -> time), fixed scaling cannot be
  done, as two 0dB coefficients would sum to a peak amplitude
  of 64K, overflowing the 32k range of the fixed-point integers.
  Thus, the fix_fft() routine performs variable scaling, and
  returns a value which is the number of bits LEFT by which
  the output must be shifted to get the actual amplitude
  (i.e. if fix_fft() returns 3, each value of fr[] and fi[]
  must be multiplied by 8 (2**3) for proper scaling.
  Clearly, this cannot be done within fixed-point short
  integers. In practice, if the result is to be used as a
  filter, the scale_shift can usually be ignored, as the
  result will be approximately correctly normalized as is.

  Written by:  Tom Roberts  11/8/89
  Made portable:  Malcolm Slaney 12/15/94 malcolm@interval.com
  Enhanced:  Dimitrios P. Bouras  14 Jun 2006 dbouras@ieee.org
*/
#include "fix_fft.h"
#define N_WAVE      1024    /* full length of Sinewave[] */
#define LOG2_N_WAVE 10      /* log2(N_WAVE) */

/*
  Henceforth "short" implies 16-bit word. If this is not
  the case in your architecture, please replace "short"
  with a type definition which *is* a 16-bit word.
*/

/*
  Since we only use 3/4 of N_WAVE, we define only
  this many samples, in order to conserve data space.
*/
short Sinewave[N_WAVE-N_WAVE/4] = {
      0,    201,    402,    603,    804,   1005,   1206,   1406,
   1607,   1808,   2009,   2209,   2410,   2610,   2811,   3011,
   3211,   3411,   3611,   3811,   4011,   4210,   4409,   4608,
   4807,   5006,   5205,   5403,   5601,   5799,   5997,   6195,
   6392,   6589,   6786,   6982,   7179,   7375,   7571,   7766,
   7961,   8156,   8351,   8545,   8739,   8932,   9126,   9319,
   9511,   9703,   9895,  10087,  10278,  10469,  10659,  10849,
  11038,  11227,  11416,  11604,  11792,  11980,  12166,  12353,
  12539,  12724,  12909,  13094,  13278,  13462,  13645,  13827,
  14009,  14191,  14372,  14552,  14732,  14911,  15090,  15268,
  15446,  15623,  15799,  15975,  16150,  16325,  16499,  16672,
  16845,  17017,  17189,  17360,  17530,  17699,  17868,  18036,
  18204,  18371,  18537,  18702,  18867,  19031,  19194,  19357,
  19519,  19680,  19840,  20000,  20159,  20317,  20474,  20631,
  20787,  20942,  21096,  21249,  21402,  21554,  21705,  21855,
  22004,  22153,  22301,  22448,  22594,  22739,  22883,  23027,
  23169,  23311,  23452,  23592,  23731,  23869,  24006,  24143,
  24278,  24413,  24546,  24679,  24811,  24942,  25072,  25201,
  25329,  25456,  25582,  25707,  25831,  25954,  26077,  26198,
  26318,  26437,  26556,  26673,  26789,  26905,  27019,  27132,
  27244,  27355,  27466,  27575,  27683,  27790,  27896,  28001,
  28105,  28208,  28309,  28410,  28510,  28608,  28706,  28802,
  28897,  28992,  29085,  29177,  29268,  29358,  29446,  29534,
  29621,  29706,  29790,  29873,  29955,  30036,  30116,  30195,
  30272,  30349,  30424,  30498,  30571,  30643,  30713,  30783,
  30851,  30918,  30984,  31049,  31113,  31175,  31236,  31297,
  31356,  31413,  31470,  31525,  31580,  31633,  31684,  31735,
  31785,  31833,  31880,  31926,  31970,  32014,  32056,  32097,
  32137,  32176,  32213,  32249,  32284,  32318,  32350,  32382,
  32412,  32441,  32468,  32495,  32520,  32544,  32567,  32588,
  32609,  32628,  32646,  32662,  32678,  32692,  32705,  32717,
  32727,  32736,  32744,  32751,  32757,  32761,  32764,  32766,
  32767,  32766,  32764,  32761,  32757,  32751,  32744,  32736,
  32727,  32717,  32705,  32692,  32678,  32662,  32646,  32628,
  32609,  32588,  32567,  32544,  32520,  32495,  32468,  32441,
  32412,  32382,  32350,  32318,  32284,  32249,  32213,  32176,
  32137,  32097,  32056,  32014,  31970,  31926,  31880,  31833,
  31785,  31735,  31684,  31633,  31580,  31525,  31470,  31413,
  31356,  31297,  31236,  31175,  31113,  31049,  30984,  30918,
  30851,  30783,  30713,  30643,  30571,  30498,  30424,  30349,
  30272,  30195,  30116,  30036,  29955,  29873,  29790,  29706,
  29621,  29534,  29446,  29358,  29268,  29177,  29085,  28992,
  28897,  28802,  28706,  28608,  28510,  28410,  28309,  28208,
  28105,  28001,  27896,  27790,  27683,  27575,  27466,  27355,
  27244,  27132,  27019,  26905,  26789,  26673,  26556,  26437,
  26318,  26198,  26077,  25954,  25831,  25707,  25582,  25456,
  25329,  25201,  25072,  24942,  24811,  24679,  24546,  24413,
  24278,  24143,  24006,  23869,  23731,  23592,  23452,  23311,
  23169,  23027,  22883,  22739,  22594,  22448,  22301,  22153,
  22004,  21855,  21705,  21554,  21402,  21249,  21096,  20942,
  20787,  20631,  20474,  20317,  20159,  20000,  19840,  19680,
  19519,  19357,  19194,  19031,  18867,  18702,  18537,  18371,
  18204,  18036,  17868,  17699,  17530,  17360,  17189,  17017,
  16845,  16672,  16499,  16325,  16150,  15975,  15799,  15623,
  15446,  15268,  15090,  14911,  14732,  14552,  14372,  14191,
  14009,  13827,  13645,  13462,  13278,  13094,  12909,  12724,
  12539,  12353,  12166,  11980,  11792,  11604,  11416,  11227,
  11038,  10849,  10659,  10469,  10278,  10087,   9895,   9703,
   9511,   9319,   9126,   8932,   8739,   8545,   8351,   8156,
   7961,   7766,   7571,   7375,   7179,   6982,   6786,   6589,
   6392,   6195,   5997,   5799,   5601,   5403,   5205,   5006,
   4807,   4608,   4409,   4210,   4011,   3811,   3611,   3411,
   3211,   3011,   2811,   2610,   2410,   2209,   2009,   1808,
   1607,   1406,   1206,   1005,    804,    603,    402,    201,
      0,   -201,   -402,   -603,   -804,  -1005,  -1206,  -1406,
  -1607,  -1808,  -2009,  -2209,  -2410,  -2610,  -2811,  -3011,
  -3211,  -3411,  -3611,  -3811,  -4011,  -4210,  -4409,  -4608,
  -4807,  -5006,  -5205,  -5403,  -5601,  -5799,  -5997,  -6195,
  -6392,  -6589,  -6786,  -6982,  -7179,  -7375,  -7571,  -7766,
  -7961,  -8156,  -8351,  -8545,  -8739,  -8932,  -9126,  -9319,
  -9511,  -9703,  -9895, -10087, -10278, -10469, -10659, -10849,
-11038, -11227, -11416, -11604, -11792, -11980, -12166, -12353,
-12539, -12724, -12909, -13094, -13278, -13462, -13645, -13827,
-14009, -14191, -14372, -14552, -14732, -14911, -15090, -15268,
-15446, -15623, -15799, -15975, -16150, -16325, -16499, -16672,
-16845, -17017, -17189, -17360, -17530, -17699, -17868, -18036,
-18204, -18371, -18537, -18702, -18867, -19031, -19194, -19357,
-19519, -19680, -19840, -20000, -20159, -20317, -20474, -20631,
-20787, -20942, -21096, -21249, -21402, -21554, -21705, -21855,
-22004, -22153, -22301, -22448, -22594, -22739, -22883, -23027,
-23169, -23311, -23452, -23592, -23731, -23869, -24006, -24143,
-24278, -24413, -24546, -24679, -24811, -24942, -25072, -25201,
-25329, -25456, -25582, -25707, -25831, -25954, -26077, -26198,
-26318, -26437, -26556, -26673, -26789, -26905, -27019, -27132,
-27244, -27355, -27466, -27575, -27683, -27790, -27896, -28001,
-28105, -28208, -28309, -28410, -28510, -28608, -28706, -28802,
-28897, -28992, -29085, -29177, -29268, -29358, -29446, -29534,
-29621, -29706, -29790, -29873, -29955, -30036, -30116, -30195,
-30272, -30349, -30424, -30498, -30571, -30643, -30713, -30783,
-30851, -30918, -30984, -31049, -31113, -31175, -31236, -31297,
-31356, -31413, -31470, -31525, -31580, -31633, -31684, -31735,
-31785, -31833, -31880, -31926, -31970, -32014, -32056, -32097,
-32137, -32176, -32213, -32249, -32284, -32318, -32350, -32382,
-32412, -32441, -32468, -32495, -32520, -32544, -32567, -32588,
-32609, -32628, -32646, -32662, -32678, -32692, -32705, -32717,
-32727, -32736, -32744, -32751, -32757, -32761, -32764, -32766,
};

/*
  FIX_MPY() - fixed-point multiplication & scaling.
  Substitute inline assembly for hardware-specific
  optimization suited to a particluar DSP processor.
  Scaling ensures that result remains 16-bit.
*/
short FIX_MPY(short a, short b)
{
        /* shift right one less bit (i.e. 15-1) */
        int c = ((int)a * (int)b) >> 14;
        /* last bit shifted out = rounding-bit */
        b = c & 0x01;
        /* last shift + rounding bit */
        a = (c >> 1) + b;
        return a;
}

/*
  fix_fft() - perform forward/inverse fast Fourier transform.
  fr[n],fi[n] are real and imaginary arrays, both INPUT AND
  RESULT (in-place FFT), with 0 <= n < 2**m; set inverse to
  0 for forward transform (FFT), or 1 for iFFT.
*/
int fix_fft(short fr[], short fi[], short m, short inverse)
{
        int mr, nn, i, j, l, k, istep, n, scale, shift;
        short qr, qi, tr, ti, wr, wi;

        n = 1 << m;

        /* max FFT size = N_WAVE */
        if (n > N_WAVE)
                return -1;

        mr = 0;
        nn = n - 1;
        scale = 0;

        /* decimation in time - re-order data */
        for (m=1; m<=nn; ++m) {
                l = n;
                do {
                        l >>= 1;
                } while (mr+l > nn);
                mr = (mr & (l-1)) + l;

                if (mr <= m)
                        continue;
                tr = fr[m];
                fr[m] = fr[mr];
                fr[mr] = tr;
                ti = fi[m];
                fi[m] = fi[mr];
                fi[mr] = ti;
        }

        l = 1;
        k = LOG2_N_WAVE-1;
        while (l < n) {
                if (inverse) {
                        /* variable scaling, depending upon data */
                        shift = 0;
                        for (i=0; i<n; ++i) {
                                j = fr;
                                if (j < 0)
                                        j = -j;
                                m = fi;
                                if (m < 0)
                                        m = -m;
                                if (j > 16383 || m > 16383) {
                                        shift = 1;
                                        break;
                                }
                        }
                        if (shift)
                                ++scale;
                } else {
                        /*
                          fixed scaling, for proper normalization --
                          there will be log2(n) passes, so this results
                          in an overall factor of 1/n, distributed to
                          maximize arithmetic accuracy.
                        */
                        shift = 1;
                }
                /*
                  it may not be obvious, but the shift will be
                  performed on each data point exactly once,
                  during this pass.
                */
                istep = l << 1;
                for (m=0; m<l; ++m) {
                        j = m << k;
                        /* 0 <= j < N_WAVE/2 */
                        wr =  Sinewave[j+N_WAVE/4];
                        wi = -Sinewave[j];
                        if (inverse)
                                wi = -wi;
                        if (shift) {
                                wr >>= 1;
                                wi >>= 1;
                        }
                        for (i=m; i<n; i+=istep) {
                                j = i + l;
                                tr = FIX_MPY(wr,fr[j]) - FIX_MPY(wi,fi[j]);
                                ti = FIX_MPY(wr,fi[j]) + FIX_MPY(wi,fr[j]);
                                qr = fr;
                                qi = fi;
                                if (shift) {
                                        qr >>= 1;
                                        qi >>= 1;
                                }
                                fr[j] = qr - tr;
                                fi[j] = qi - ti;
                                fr = qr + tr;
                                fi = qi + ti;
                        }
                }
                --k;
                l = istep;
        }
        return scale;
}

/*
  fix_fftr() - forward/inverse FFT on array of real numbers.
  Real FFT/iFFT using half-size complex FFT by distributing
  even/odd samples into real/imaginary arrays respectively.
  In order to save data space (i.e. to avoid two arrays, one
  for real, one for imaginary samples), we proceed in the
  following two steps: a) samples are rearranged in the real
  array so that all even samples are in places 0-(N/2-1) and
  all imaginary samples in places (N/2)-(N-1), and b) fix_fft
  is called with fr and fi pointing to index 0 and index N/2
  respectively in the original array. The above guarantees
  that fix_fft "sees" consecutive real samples as alternating
  real and imaginary samples in the complex array.
*/
int fix_fftr(short f[], int m, int inverse)
{
        int i, N = 1<<(m-1), scale = 0;
        short tt, *fr=f, *fi=&f[N];

        if (inverse)
                scale = fix_fft(fi, fr, m-1, inverse);
        for (i=1; i<N; i+=2) {
                tt = f[N+i-1];
                f[N+i-1] = f;
                f = tt;
        }
        if (! inverse)
                scale = fix_fft(fi, fr, m-1, inverse);
        return scale;
}

prince2010

senjet 发表于 2014-1-15 19:11
以下是用过的16bit定点FFT，带防溢出的，看看是否对你有用。
/* fix_fft.c - Fixed-point in-place Fast Fo ...

我看看

gaolf_2012

#include <arm_math.h>

#define LENGTH       64                 //定义做FFT的点数
#define HALF_LENGTH  (LENGTH / 2)       //FFT一半的点数

#define CFFT_FORWARD        0
#define CFFT_INVERSE        1
#define CFFT_NO_BITREVERSAL 0
#define CFFT_BITREVERSAL    1

// Global Sample Data
arm_cfft_radix4_instance_q31 S;
// Instance structure for Q31 CFFT/CIFFT function
q31_t q31SampleData[2 * LENGTH]; // Data format: real[0}, imag[0], real[1], imag[1], ...

INT32S dwpRe[LENGTH];          //存放结果的实部
INT32S dwpIm[LENGTH];          //存放结果的虚部
INT16S wpSampling[LENGTH];     //存放采样数据

void test(void)
{
        INT16U i;

        arm_cfft_radix4_init_q31(&S, LENGTH, CFFT_FORWARD, CFFT_BITREVERSAL); //蝶形因子
        for(i = 0; i < LENGTH; i++)
        {
                q31SampleData[2*i] = (q31_t)wpSampling[i];// real input data
                q31SampleData[2*i + 1] = 0;       // imaginary input data (clear to 0)
        }
        arm_cfft_radix4_q31(&S, q31SampleData);
        for(i=0; i<HALF_LENGTH; i++)
        {
                dwpReal[i] = q31SampleData[2*i];
                dwpImag[i] = q31SampleData[2*i + 1];
        }
}

prince2010

gaolf_2012 发表于 2014-1-15 19:54
#include

#define LENGTH       64                 //定义做FFT的点数

正在测试........

dianyuan

看下这个
http://www.amobbs.com/forum.php?mod=viewthread&tid=4472775

cjy3063

senjet 发表于 2014-1-15 19:11
以下是用过的16bit定点FFT，带防溢出的，看看是否对你有用。
/* fix_fft.c - Fixed-point in-place Fast Fo ...

你好！senjet

这段程序里有几个地方不是很清楚，请帮忙指点下，谢谢！

for (i=m; i<n; i+=istep) {
                                j = i + l;
                                tr = FIX_MPY(wr,fr[j]) - FIX_MPY(wi,fi[j]);
                                ti = FIX_MPY(wr,fi[j]) + FIX_MPY(wi,fr[j]);
                                qr = fr;        //[这里的fr 下标是什么？i,j,还是？]
                                qi = fi;        //[这里的fr 下标是什么？i,j,还是？]
                                if (shift) {
                                        qr >>= 1;
                                        qi >>= 1;
                                }
                                fr[j] = qr - tr;
                                fi[j] = qi - ti;
                                fr = qr + tr;        //[这里的fr 下标是什么？i,j,还是？]
                                fi = qi + ti;                //[这里的fr 下标是什么？i,j,还是？]

                        }

cjy3063

prince2010 发表于 2014-1-16 08:02
正在测试........

你好！
测试可以了？
我在测试过程中发现一点问题，不知道你怎么处理的？


for (i=m; i<n; i+=istep) {
                                j = i + l;
                                tr = FIX_MPY(wr,fr[j]) - FIX_MPY(wi,fi[j]);
                                ti = FIX_MPY(wr,fi[j]) + FIX_MPY(wi,fr[j]);
                                qr = fr;//[这里的fr 下标是什么？i,j,还是？]
                                qi = fi;//[这里的fr 下标是什么？i,j,还是？]
                                if (shift) {
                                        qr >>= 1;
                                        qi >>= 1;
                                }
                                fr[j] = qr - tr;
                                fi[j] = qi - ti;
                                fr = qr + tr;//[这里的fr 下标是什么？i,j,还是？]
                                fi = qi + ti;//[这里的fr 下标是什么？i,j,还是？]
                        }