| 网站首页 | 业界新闻 | 小组 | 威客 | 人才 | 下载频道 | 博客 | 代码贴 | 在线编程 | 编程论坛
欢迎加入我们,一同切磋技术
用户名:   
 
密 码:  
共有 1915 人关注过本帖
标题:快速傅立叶变换FFT C程序的疑问
只看楼主 加入收藏
hczsea
Rank: 2
等 级:论坛游民
帖 子:129
专家分:68
注 册:2007-10-23
结帖率:100%
收藏
已结贴  问题点数:20 回复次数:21 
快速傅立叶变换FFT C程序的疑问
最近在看《数字信号处理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");
    }
}
搜索更多相关主题的帖子: include double C语言 
2013-07-24 10:46
TonyDeng
Rank: 20Rank: 20Rank: 20Rank: 20Rank: 20
等 级:贵宾
威 望:304
帖 子:25859
专家分:48889
注 册:2011-6-22
收藏
得分:0 
你已经发现了问题,麻烦你尊重一下别人,把发现的现象描述一遍,不要让人浪费时间重复你的发现,基本的礼貌。
此帖锁定!

授人以渔,不授人以鱼。
2013-07-24 10:50
rjsp
Rank: 20Rank: 20Rank: 20Rank: 20Rank: 20
等 级:版主
威 望:528
帖 子:9025
专家分:54030
注 册:2011-1-18
收藏
得分:0 
运行后的结果与书上的不一致
------ 运行后的结果是什么?书上的结果是什么?尤其是后者,你不说别人知道个球呀!
2013-07-24 16:10
beyondyf
Rank: 19Rank: 19Rank: 19Rank: 19Rank: 19Rank: 19
等 级:贵宾
威 望:103
帖 子:3282
专家分:12654
注 册:2008-1-21
收藏
得分:20 
这贴没什么问题,锁贴的理由很牵强。如果这贴都要锁,那这论坛里没几个值得打开的贴子了。

这段代码没有语法错误,所以也没什么现象可以描述。楼主说的也很清楚,程序的输出与书中的结果不一致。如果非要个现象,这就是现象。

输出结果与标准结果不一致这一定是代码中存在逻辑错误,这种错误往往很难发现。这题有标准结果可对比还好,而实际应用中多数没有标准结果可对比,或者是只有部分样本而程序的逻辑错误又只出现在部分问题域内,也就是说大部分情况下它是对的,只在少数点上会出错。说所有的软件都存在这样的错误有点绝对,但这样的错误绝不在少数,这也是我们时常要给操作系统打补丁的原因。

书归正传,仔细阅读了一下楼主的代码,其中FFT函数中存在一处逻辑错误

程序代码:
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]; //应为ti = y[j];
            x[j] = x[i];
            y[j] = y[i];
            x[i] = tr;
            y[i] = ti;
        }
    ...


另外想说说这段代码,太老了,从函数参数的定义格式以及上面FFT中第一个循环的迭代语句“i < 16”可以判断,这段代码写就与早期的16位编译器,估计就是TC2.0。我不知道现在的编译器是否还支持这样的参数定义格式,该废止了。

再者,#include "fft.c",这样的用法太不专业了,存在大量隐患。

原因解释起来话就长了,总之新手记住只用来包含.h头文件。

在自己写头文件时,里面只包含函数声明、结构定义、常量/参数宏,不要把函数定义、全局变量定义放进去。

重剑无锋,大巧不工
2013-07-24 16:24
beyondyf
Rank: 19Rank: 19Rank: 19Rank: 19Rank: 19Rank: 19
等 级:贵宾
威 望:103
帖 子:3282
专家分:12654
注 册:2008-1-21
收藏
得分:0 
回复 3楼 rjsp
你介不介意我删掉你3楼的发言?

重剑无锋,大巧不工
2013-07-24 16:29
岁月如刀
Rank: 6Rank: 6
来 自:冰冻星球
等 级:侠之大者
威 望:7
帖 子:165
专家分:477
注 册:2013-7-21
收藏
得分:0 
以下是引用beyondyf在2013-7-24 16:29:42的发言:

你介不介意我删掉你3楼的发言?

我靠。确实牛逼啊。我想去采访一下以上两位是什么心情~


Now. Game over.
2013-07-24 17:08
beyondyf
Rank: 19Rank: 19Rank: 19Rank: 19Rank: 19Rank: 19
等 级:贵宾
威 望:103
帖 子:3282
专家分:12654
注 册:2008-1-21
收藏
得分:0 
回复 6楼 岁月如刀
你介不介意我删掉你6楼的发言?

重剑无锋,大巧不工
2013-07-24 17:17
love云彩
Rank: 19Rank: 19Rank: 19Rank: 19Rank: 19Rank: 19
来 自:青藏高原
等 级:贵宾
威 望:53
帖 子:3663
专家分:11416
注 册:2012-11-17
收藏
得分:0 
以下是引用beyondyf在2013-7-24 16:24:07的发言:

这贴没什么问题,锁贴的理由很牵强。如果这贴都要锁,那这论坛里没几个值得打开的贴子了。
 
这段代码没有语法错误,所以也没什么现象可以描述。楼主说的也很清楚,程序的输出与书中的结果不一致。如果非要个现象,这就是现象。
 
输出结果与标准结果不一致这一定是代码中存在逻辑错误,这种错误往往很难发现。这题有标准结果可对比还好,而实际应用中多数没有标准结果可对比,或者是只有部分样本而程序的逻辑错误又只出现在部分问题域内,也就是说大部分情况下它是对的,只在少数点上会出错。说所有的软件都存在这样的错误有点绝对,但这样的错误绝不在少数,这也是我们时常要给操作系统打补丁的原因。
 
书归正传,仔细阅读了一下楼主的代码,其中FFT函数中存在一处逻辑错误
 
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
帖子还好吧,就是问题描述不太好,至少让新手可以借鉴一下

思考赐予新生,时间在于定义
2013-07-25 02:08
beyondyf
Rank: 19Rank: 19Rank: 19Rank: 19Rank: 19Rank: 19
等 级:贵宾
威 望:103
帖 子:3282
专家分:12654
注 册:2008-1-21
收藏
得分:0 
回复 8楼 love云彩
你想借鉴什么?这份代码写的很一般,当然这不是楼主的问题,他也是从书上抄的代码想试一下。原书作者是信号处理方面的专家,但不是编程的高手。

所以这贴没什么可借鉴的编程技巧。想看懂这份代码光有C语言的知识是不够的,更重要的是懂傅里叶变换。

这个论坛里的众人大部分是大学生吧,关于傅里叶级数的知识在高数下册中有讲解。你们需要的是重新复习一下这部分知识,然后找本数值分析或信号处理的书再看看快速傅里叶变换的算法原理。

重剑无锋,大巧不工
2013-07-25 08:57
rjsp
Rank: 20Rank: 20Rank: 20Rank: 20Rank: 20
等 级:版主
威 望:528
帖 子:9025
专家分:54030
注 册:2011-1-18
收藏
得分:0 
以下是引用beyondyf在2013-7-24 16:29:42的发言:

你介不介意我删掉你3楼的发言?
我觉得“运行后的结果与书上的不一致”并不是个现象的描述,因为文中没有定义“运行后的结果”和“书上的”,起码是个不完整的现象描述。

我模仿问一题:
int main() { return 0; }
运行后的结果与我猜想的不一致。请高手帮忙看看是咋回事。
2013-07-25 09:00
快速回复:快速傅立叶变换FFT C程序的疑问
数据加载中...
 
   



关于我们 | 广告合作 | 编程中国 | 清除Cookies | TOP | 手机版

编程中国 版权所有,并保留所有权利。
Powered by Discuz, Processed in 0.038565 second(s), 7 queries.
Copyright©2004-2024, BCCN.NET, All Rights Reserved