ttxs_2013 发表于 2014-5-13 21:25

FFT的实现不明白,希望指点下

在学FFT,有两个地方没有弄明白
1 蝶形计算没看明白,希望能够详细解释下
dataR=dataR+dataR*cos_tab+dataI*sin_tab;
dataI=dataI-dataR*sin_tab+dataI*cos_tab;
dataR=TR-dataR*cos_tab-dataI*sin_tab;
dataI=TI+temp*sin_tab-dataI*cos_tab;
2 makewave函数中
        INPUT=sin(PI*2*i/SAMPLENUMBER*3)*1024;
乘以3再乘1024 是什么意思?
***********************************************************************
#include "myapp.h"
#include "csedu.h"
#include "scancode.h"
#include <math.h>

#define PI 3.1415926
#define SAMPLENUMBER 128

void InitForFFT();
void MakeWave();

int INPUT,DATA;
float fWaveR,fWaveI,w;
float sin_tab,cos_tab;

main()
{
        int i;
       
        InitForFFT();
        MakeWave();
        for ( i=0;i<SAMPLENUMBER;i++ )
        {
                fWaveR=INPUT;
                fWaveI=0.0f;
                w=0.0f;
        }
        FFT(fWaveR,fWaveI);//地址传递,数组首地址传递给dataR,dataI
        for ( i=0;i<SAMPLENUMBER;i++ )
        {
                DATA=w;
        }
        while ( 1 );        // break point
}
      
void FFT(float dataR,float dataI)
{
        int x0,x1,x2,x3,x4,x5,x6,xx;
        int i,j,k,b,p,L;
        float TR,TI,temp;
       
        /********** following code invert sequence ************/
        /*位倒序算法:128点FFT可用7位二进制表示顺序数,这里依次为x6-x0(高到低)
          x0到x6依次与0x01相与,得到顺序数i用7位二进制表示,xx则是计算倒序二进制的十进制数
          dataR数组是存放顺序数的(输入),dataI数组则是存放倒序数的
          注意 i/(2,4,8,16等)相当于i二进制右移1,2,3,4位等*/
        for ( i=0;i<SAMPLENUMBER;i++ )
        {
                x0=x1=x2=x3=x4=x5=x6=0;
                x0=i&0x01; x1=(i/2)&0x01; x2=(i/4)&0x01; x3=(i/8)&0x01;x4=(i/16)&0x01; x5=(i/32)&0x01; x6=(i/64)&0x01;
                xx=x0*64+x1*32+x2*16+x3*8+x4*4+x5*2+x6;
                dataI=dataR;   //dataIa 数组中存放的就是位倒序后的值   
        }
        for ( i=0;i<SAMPLENUMBER;i++ )
        {
                dataR=dataI; dataI=0; //把位倒序数组中的值逐位的赋值给原输入数组,这样就少用一个数组
        }

        /************** following code FFT *******************/
        for ( L=1;L<=7;L++ )      //第一层循环:蝶形的级数
        { /* for(1) */         
                b=1; i=L-1;
                while ( i>0 )
                {
                        b=b*2; i--;
                } /* b= 2^(L-1) */
       //第二层循环:获得旋转因子 P,第L级有2的(L-1)个不同的旋转因子,同时要计算每一个旋转因子对应的蝶形
                for ( j=0;j<=b-1;j++ ) /* for (2) */
                {
                        p=1; i=7-L;
                        while ( i>0 ) /* p=pow(2,7-L)*j; */
                        {
                                p=p*2; i--;
                        }
                        p=p*j;
         /*第三层循环:计算蝶形 ,第L级旋转因子的个数同每一级中每个蝶形的两个输入数据相差的点数 相同*/
               /*实际上这里b的值每一级中都是固定的,与循环的次数并没有关系*/
                        for ( k=j;k<128;k=k+2*b ) /* for (3) */
                        {
                                TR=dataR; //已经是位倒序过了
                                TI=dataI; // dataI中数据都赋值为0了,这里还有什么作用?
                                temp=dataR; // temp是蝶形中另一个输入
                       /****************************************************************/
                                dataR=dataR+dataR*cos_tab+dataI*sin_tab;
                                dataI=dataI-dataR*sin_tab+dataI*cos_tab;
                                dataR=TR-dataR*cos_tab-dataI*sin_tab;
                                dataI=TI+temp*sin_tab-dataI*cos_tab;
                        } /* END for (3) */
                } /* END for (2) */
        } /* END for (1) */
        for ( i=0;i<SAMPLENUMBER/2;i++ )
        {
                w=sqrt(dataR*dataR+dataI*dataI);
        }
} /* END FFT */
               

void InitForFFT()
{
        int i;
       
        for ( i=0;i<SAMPLENUMBER;i++ )
        {
                sin_tab=sin(PI*2*i/SAMPLENUMBER);
                cos_tab=cos(PI*2*i/SAMPLENUMBER);
        }
}

void MakeWave()
{
        int i;
       
        for ( i=0;i<SAMPLENUMBER;i++ )
        {
                INPUT=sin(PI*2*i/SAMPLENUMBER*3)*1024;
        }
}

zhangmangui 发表于 2014-5-13 21:37

没深入研究    向你学习啦

aresc 发表于 2014-5-13 22:34

本帖最后由 aresc 于 2014-5-13 22:39 编辑

1. INPUT=sin(PI*2*i/SAMPLENUMBER*3)*1024;
*1024实际没啥特别的含义,只是把一个范围为[-1.0,1.0]的正弦信号放大到了[-1024.0,1024.0]的范围.
那个3是正弦信号的频率,采样率为SAMPLENUMBER.

2. 关于蝶形运算建议你找本信号处理的书好好看看,你用的是时域抽取的算法,和频域抽取的蝶形运算是有差别的。几句话写不清楚,呵呵。关键就是旋转因子然后用欧拉恒等式展开成(cosA+jsinA)做复数乘、加。

3. DataI为0的为什么:因为你用的是复数运算的FFT,但你给的输入信号是一个实数的正弦信号,相当于复数信号只有实部,虚部为0。

NWPU_CHEN 发表于 2014-5-14 09:20

楼上说的不错,想要了解找本数字信号处理好好看看,就懂了,这东西一两句说不清

zhangmangui 发表于 2014-5-14 14:52

aresc 发表于 2014-5-13 22:34
1. INPUT=sin(PI*2*i/SAMPLENUMBER*3)*1024;
*1024实际没啥特别的含义,只是把一个范围为[-1.0,1.0]的正弦 ...

很给力啊

ttxs_2013 发表于 2014-5-14 19:37

NWPU_CHEN 发表于 2014-5-14 09:20 static/image/common/back.gif
楼上说的不错,想要了解找本数字信号处理好好看看,就懂了,这东西一两句说不清 ...

嗯,就是看书来理解的,但是理论和具体的代码实现还是有区别的。。。。。

aresc 发表于 2014-5-14 22:23

FFT算法比较难理解的是最终算法实现是和算法的推导是反过来的。比如N点FFT,按奇偶分解成两个N/2点,...最后拆成N/2个2点。最后算法实现的时候是反过来的,即先做最里层的N/2个两点的FFT,最后是两个N/2点的FFT。如果按算法推导的过程写程序只能用递归函数来实现,用非递归的方法就是现在大家都在采用的程序。

ttxs_2013 发表于 2014-5-14 22:56

aresc 发表于 2014-5-13 22:34 static/image/common/back.gif
1. INPUT=sin(PI*2*i/SAMPLENUMBER*3)*1024;
*1024实际没啥特别的含义,只是把一个范围为[-1.0,1.0]的正弦 ...

又分析了下,但还是要追问几个问题:
(1) INPUT中的3好像不是频率,应该是周期,因为按照对模拟输入信号进行采样,得到离散序列,进而进行DFT变换。Xa(t)=sin(2π f t) 令t=nT,T=1/fs,那么有Xa(t)=sin(2π f n/ fs),这样看来,3就是周期,而不是频率了。
(2) 我没有明白的地方应该是这里的 旋转因子用欧拉公式表示,后代入到蝶形运算中(之前不疑惑的地方同样是这个地方),比如
dataR=dataR+dataR*cos_tab+dataI*sin_tab; // 这个可以对应起来
dataR=TR-dataR*cos_tab-dataI*sin_tab;   //这个按照欧拉展开,sin-tab 不应该是加吗?

dataI=dataI-dataR*sin_tab+dataI*cos_tab;
dataI=TI+temp*sin_tab-dataI*cos_tab;               // 这两句不明白,或者说dataI 这里表示的是什么?
简单来说,这里的复数运算 没有看明白,希望前辈有针对性的解释下。

aresc 发表于 2014-5-15 12:33

1. sin(PI*2*i/SAMPLENUMBER*3) = sin(2*PI*3*i/SAMPLENUMBER) = sin(2*PI*f*i/fs), 所以f/fs=3/SAMPLENUMBER.
2. 正向FFT时,欧拉公式展开后是cos(-A)+j*sin(-A) = cosA-jsinA.
时域抽取的一个蝶形包含运算: a' = a+b*(cosA-jsinA), b' = a - b*(cosA-jsinA). 其中a、b都是复数。展开之后就是分别求a',b'的实部、虚部,即4个运算式。

ttxs_2013 发表于 2014-5-15 19:41

aresc 发表于 2014-5-14 22:23 static/image/common/back.gif
FFT算法比较难理解的是最终算法实现是和算法的推导是反过来的。比如N点FFT,按奇偶分解成两个N/2点,...最 ...

是的,不结合具体的FFT流图来理解,很难看懂,我是结合8点的FFT流图来分析的,和代码是一致的。

icekoor 发表于 2014-8-8 15:39

最近也在研究FFT,这个基数为2的算法真心不错,网上很多人都在用,花费了半天时间才看懂。

zhangmangui 发表于 2014-8-9 15:43

icekoor 发表于 2014-8-8 15:39 static/image/common/back.gif
最近也在研究FFT,这个基数为2的算法真心不错,网上很多人都在用,花费了半天时间才看懂。 ...

看懂了抽时间总结分享一下啊

王亚威 发表于 2021-12-19 14:07

你好
请问你现在弄明白了吗
页: [1]
查看完整版本: FFT的实现不明白,希望指点下