已结贴√
问题点数:20 回复次数:21
快速傅立叶变换FFT C程序的疑问
最近在看《数字信号处理C语言程序集》,其中1.2 快速傅立叶变换的程序fft.c,按照书上的程序和例题,运行后的结果与书上的不一致。请高手帮忙看看是咋回事。
fft.c如下:
#include
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
{
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
{
for (i=j; i
{
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
{
x[i] /= n;
y[i] /= n;
}
}
}
主程序如下:
#include
#include
#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
{
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
{
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
{
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
{
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
{
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
{
for (j=0; j<2; j++)
{
printf(" %10.7f + J %10.7f",x[2*i+j],y[2*i+j]);
}
printf("\n");
}
}