[music-DSP] fixed point FFT

[yogurthu at arnet.com.ar]

I\'ve been searching the subject for implementations on microcontrollers. Have found a suitable algorithm but can\'t remember where. I\'ll paste just the core here, email me if you need further info.
BTW, if you feel adventurous, there is a great implementation of an integer (block floating point) FFT for the ADSP21xx on the Analog Devices DSP books, available at their web site. It is in 21xx assembler, though. There are old TMS320 sources also around.
Unless you resort to block floating point, you\'ll have to make sure there are enough bits on the left to accomodate room for bit growth inside the butterflies, integer FFT algos can easily overflow

int fix_fft(fixed fr[], fixed fi[], int m)
{
        int mr,nn,i,j,l,k,istep, n, scale, shift;
        fixed qr,qi,tr,ti,wr,wi,t;

        n = 1<<m;            

        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) {
                        /* 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. */

                /* 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];

                        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[i];
                                qi = fi[i];

                                qr >>= 1;
                                qi >>= 1;

                                fr[j] = qr - tr;
                                fi[j] = qi - ti;
                                fr[i] = qr + tr;
                                fi[i] = qi + ti;
                        }
                }
                --k;
                l = istep;
        }

        return scale;
}

--
Sergio R. Caprile, Electronics Engineer
http://www.geocities.com/scaprile

__________________________________
Registrate desde http://servicios.arnet.com.ar/registracion/registracion.asp?origenid=9 y participá de todos los beneficios del Portal Arnet.