最近在看《数字信号处理C语言程序集》，其中1.2 快速傅立叶变换的程序fft.c，按照书上的程序和例题，运行后的结果与书上的不一致。请高手帮忙看看是咋回事。
fft.c如下：
#include <math.h>

void fft(x, y, n, sign)
int n, sign;
double x[], y[];
{
    int i, j, k, l, m, n1, n2;
    double c, c1, e, s, s1, t, tr, ti;

    for (j=1,i=1; i<16; i++)
    {
        m = i;
        j = 2 * j;
        if (j == n)
        {
            break;
        }
    }

    n1 = n - 1;
    for (j=0,i=0; i<n1; i++)
    {
        if (i < j)
        {
            tr = x[j];
            ti = y[i];
            x[j] = x[i];
            y[j] = y[i];
            x[i] = tr;
            y[i] = ti;
        }

        k = n / 2;
        while (k < (j+1))
        {
            j = j - k;
            k = k / 2;
        }

        j = j + k;
    }

    n1 = 1;
    for (l=1; l<=m; l++)
    {
        n1 = 2 * n1;
        n2 = n1 / 2;
        e = 3.14159265359 / n2;
        c = 1.0;
        s = 0.0;
        c1 = cos(e);
        s1 = -sign * sin(e);

        for (j=0; j<n2; j++)
        {
            for (i=j; i<n; i+=n1)
            {
                k = i + n2;
                tr = c * x[k] - s * y[k];
                ti = c * y[k] + s * x[k];
                x[k] = x[i] - tr;
                y[k] = y[i] - ti;
                x[i] = x[i] + tr;
                y[i] = y[i] + ti;
            }
        

            t = c;
            c = c * c1 - s * s1;
            s = t * s1 + s * c1;
        }
    }

    if (sign == -1)
    {
        for (i=0; i<n; i++)
        {
            x[i] /= n;
            y[i] /= n;
        }
    }
}

主程序如下：
#include <math.h>
#include <stdio.h>
#include "fft.c"

main()
{
    int i, j, n;
    double a1, a2, c, c1, c2, d1, d2, q1, q2, w, w1, w2;
    double x[32], y[32], a[32], b[32];

    n = 32;
    a1 = 0.9;
    a2 = 0.3;
    x[0] = 1.0;
    y[0] = 0.0;

    for (i=1; i<n; i++)
    {
        x[i] = a1 * x[i-1] - a2 * y[i-1];
        y[i] = a2 * x[i-1] + a1 * y[i-1];
    }

    printf("\nCOMPLEX INPUT SEQUENCE\n");
    for (i=0; i<n/2; i++)
    {
        for (j=0; j<2; j++)
        {
            printf("  %10.7f + J %10.7f",x[2*i+j],y[2*i+j]);
        }
        printf("\n");
    }

    q1 = x[n-1];
    q2 = y[n-1];
    for (i=0; i<n; i++)
    {
        w = 6.28318530718 / n * i;
        w1 = cos(w);
        w2 = -sin(w);
        c1 = 1.0 - a1 * w1 + a2 * w2;
        c2 = a1 * w2 + a2 * w1;
        c = c1 * c1 + c2 * c2;
        d1 = 1.0 - a1 * q1 + a2 * q2;
        d2 = a1 * q2 + a2 * q1;
        a[i] = (c1 * d1 + c2 * d2) / c;
        b[i] = (c2 * d1 - c1 * d2) / c;
    }   

    printf("\nTHEORETICAL DFT\n");
    for (i=0; i<n/2; i++)
    {
        for (j=0; j<2; j++)
        {
            printf("  %10.7f + J%10.7f",a[2*i+j],b[2*i+j]);
        }
        printf("\n");
    }

    fft(x, y, n, 1);

    printf("\nDISCRETE FOUTIER TRANSFORM\n");
    for (i=0; i<n/2; i++)
    {
        for (j=0; j<2; j++)
        {
            printf("  %10.7f + J %10.7f",x[2*i+j],y[2*i+j]);
        }
        printf("\n");
    }

    fft(x, y, n, -1);

    printf("\nINVERSE DISCRETE FOUTIER TRANSFORM\n");
    for (i=0; i<n/2; i++)
    {
        for (j=0; j<2; j++)
        {
            printf("  %10.7f + J %10.7f",x[2*i+j],y[2*i+j]);
        }
        printf("\n");
    }
}

好了，结帖了。自己找到问题了，跟4楼的相同，感谢！
顺便说一句，这个论坛的版主不需要这么多，充数的没意思，宁缺勿滥吧，不然会毁了一个好的论坛。