打印
[DSP编程]

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

[复制链接]
6387|12
手机看帖
扫描二维码
随时随地手机跟帖
跳转到指定楼层
楼主
ttxs_2013|  楼主 | 2014-5-13 21:25 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
在学FFT,有两个地方没有弄明白
1 蝶形计算没看明白,希望能够详细解释下
dataR[k]=dataR[k]+dataR[k+b]*cos_tab[p]+dataI[k+b]*sin_tab[p];
dataI[k]=dataI[k]-dataR[k+b]*sin_tab[p]+dataI[k+b]*cos_tab[p];
dataR[k+b]=TR-dataR[k+b]*cos_tab[p]-dataI[k+b]*sin_tab[p];
dataI[k+b]=TI+temp*sin_tab[p]-dataI[k+b]*cos_tab[p];
2 makewave函数中
        INPUT[i]=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[SAMPLENUMBER],DATA[SAMPLENUMBER];
float fWaveR[SAMPLENUMBER],fWaveI[SAMPLENUMBER],w[SAMPLENUMBER];
float sin_tab[SAMPLENUMBER],cos_tab[SAMPLENUMBER];

main()
{
        int i;
       
        InitForFFT();
        MakeWave();
        for ( i=0;i<SAMPLENUMBER;i++ )
        {
                fWaveR[i]=INPUT[i];
                fWaveI[i]=0.0f;
                w[i]=0.0f;
        }
        FFT(fWaveR,fWaveI);//地址传递,数组首地址传递给dataR[128],dataI[128]
        for ( i=0;i<SAMPLENUMBER;i++ )
        {
                DATA[i]=w[i];
        }
        while ( 1 );        // break point
}
      
void FFT(float dataR[SAMPLENUMBER],float dataI[SAMPLENUMBER])
{
        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[xx]=dataR[i];   //dataIa 数组中存放的就是位倒序后的值     
        }
        for ( i=0;i<SAMPLENUMBER;i++ )
        {
                dataR[i]=dataI[i]; dataI[i]=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[k]; //已经是位倒序过了
                                TI=dataI[k]; // dataI中数据都赋值为0了,这里还有什么作用?
                                temp=dataR[k+b]; // temp是蝶形中另一个输入
                         /****************************************************************/
                                dataR[k]=dataR[k]+dataR[k+b]*cos_tab[p]+dataI[k+b]*sin_tab[p];
                                dataI[k]=dataI[k]-dataR[k+b]*sin_tab[p]+dataI[k+b]*cos_tab[p];
                                dataR[k+b]=TR-dataR[k+b]*cos_tab[p]-dataI[k+b]*sin_tab[p];
                                dataI[k+b]=TI+temp*sin_tab[p]-dataI[k+b]*cos_tab[p];
                        } /* END for (3) */
                } /* END for (2) */
        } /* END for (1) */
        for ( i=0;i<SAMPLENUMBER/2;i++ )
        {
                w[i]=sqrt(dataR[i]*dataR[i]+dataI[i]*dataI[i]);
        }
} /* END FFT */
               

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

void MakeWave()
{
        int i;
       
        for ( i=0;i<SAMPLENUMBER;i++ )
        {
                INPUT[i]=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 | 只看该作者
楼上说的不错,想要了解找本数字信号处理好好看看,就懂了,这东西一两句说不清

使用特权

评论回复
5
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]的正弦 ...

很给力啊

使用特权

评论回复
6
ttxs_2013|  楼主 | 2014-5-14 19:37 | 只看该作者
NWPU_CHEN 发表于 2014-5-14 09:20
楼上说的不错,想要了解找本数字信号处理好好看看,就懂了,这东西一两句说不清 ...

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

使用特权

评论回复
7
aresc| | 2014-5-14 22:23 | 只看该作者
FFT算法比较难理解的是最终算法实现是和算法的推导是反过来的。比如N点FFT,按奇偶分解成两个N/2点,...最后拆成N/2个2点。最后算法实现的时候是反过来的,即先做最里层的N/2个两点的FFT,最后是两个N/2点的FFT。如果按算法推导的过程写程序只能用递归函数来实现,用非递归的方法就是现在大家都在采用的程序。

使用特权

评论回复
8
ttxs_2013|  楼主 | 2014-5-14 22:56 | 只看该作者
aresc 发表于 2014-5-13 22:34
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[k]=dataR[k]+dataR[k+b]*cos_tab[p]+dataI[k+b]*sin_tab[p]; // 这个可以对应起来
dataR[k+b]=TR-dataR[k+b]*cos_tab[p]-dataI[k+b]*sin_tab[p];   //  这个按照欧拉展开,sin-tab 不应该是加吗?

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

使用特权

评论回复
9
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个运算式。

使用特权

评论回复
10
ttxs_2013|  楼主 | 2014-5-15 19:41 | 只看该作者
aresc 发表于 2014-5-14 22:23
FFT算法比较难理解的是最终算法实现是和算法的推导是反过来的。比如N点FFT,按奇偶分解成两个N/2点,...最 ...

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

使用特权

评论回复
11
icekoor| | 2014-8-8 15:39 | 只看该作者
最近也在研究FFT,这个基数为2的算法真心不错,网上很多人都在用,花费了半天时间才看懂。

使用特权

评论回复
12
zhangmangui| | 2014-8-9 15:43 | 只看该作者
icekoor 发表于 2014-8-8 15:39
最近也在研究FFT,这个基数为2的算法真心不错,网上很多人都在用,花费了半天时间才看懂。 ...

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

使用特权

评论回复
13
王亚威| | 2021-12-19 14:07 | 只看该作者
你好
请问你现在弄明白了吗

使用特权

评论回复
发新帖 我要提问
您需要登录后才可以回帖 登录 | 注册

本版积分规则

31

主题

125

帖子

4

粉丝