首页 > 编程知识 正文

fft算法的matlab实现(FFT算法的完整DSP实现)

时间:2023-05-05 05:55:43 阅读:121691 作者:564

傅立叶变换或FFT的理论参考:

[1] http://www.DSP guide.com/ch12/2.htm

thescientistandengineer ' sguidetodigitalsignalprocessing,By Steven W. Smith,Ph.D。

[2] http://blog.csdn.net/v _ July _ v/article/details/6196862,可供[1]的中文参考

[3]所有数字信号处理教材都详细描述了将DCT求解转换为FFT求解的过程

[4] TI文档:基于TMS320C64x DSP的FFT实现。 可以在百度/谷歌上搜索。

1 .关于FFT理论的一点说明关于FFT,这里只提两点:

)1) DFT转换对的表达式(必须记住) )。

——称为旋转因子

)2) FFT用途——只有一个目标,提高DFT的计算效率。

DFT计算x(k )需要N^2次复数乘法和n ) n-1 )次复数加法; ft会降低N^2的计算量。

“FFT其实很难,多年来在这个领域工作的科学家也未必能很好地写FFT的算法。 ”——参考上面提供的参考文献[1]摘录

因此,不需要太在意细节。 明白了FFT理论后,沿用现有的算法就可以了。 对于在关闭教材的状态下写不了FFT这件事,没有必要太挑剔。

FFT的BASIC程序伪代码如下。

IFFT的BASIC程序伪代码如下。 IFFT是通过调用FFT计算的。

FFT算法的流程图如下图所示,归纳为三个流程三个周期。

(1) 3过程;点时域分解;反序过程;点时域计算;点谱频域合成

)2) 3个循环:外循环——分解次数、中循环——sub-DFT运算、内循环——2点蝶形算法

分解过程或反序获取参考下图理解:

2. FFT的DSP实现

以下是本人使用c语言实现的FFT及IFFT算法的例子,可以计算以任意2为对数底的采样点数的FFT。 算法参考上面列举的流程图。

/* * zx _ FFT.h * * created on :2013-8-5 * author : monkey zx */# ifndef zx _ FFT _ h _ # define zx _ FFT _ h FFT_TYPE img; } complex; intFFT(complex*x,int N ); intIFFT(complex*x,int N ); voidzx_FFT(void; #endif /* ZX_FFT_H_ */

/* * zx _ FFT.c * *实施offastfouriertransform (FFT ) * andreversalfastfouriertransform (IFFT ) * created on 330 include math.h # include stdlib.h/* * bit reverse *==input==n : nodesoffft.@ nshouldbepowerof 2,thatis2^ * * @L=ceil{log2(n ) } *==output==* r : resultsafterreversed.* note 3360 iusealocalvariable @ tempthatresult @ to @ int j=0; 短行程=0; 静态完成*时间=0; temp=(complex* ) malloc(sizeof ) complex ) * n ); if (! temp({return; (for ) I=0; in; I ) {stk=0; j=0; o{STK|=(I ) j ) )0x01; if(JL ) {stk=

1;}}while(j<l);if(stk < n) { /* 满足倒位序输出 */temp[stk] = x[i];}}/* copy @temp to @r */for (i=0; i<n; i++) {r[i] = temp[i];}free(temp);}/* * FFT Algorithm * === Inputs === * x : complex numbers * N : nodes of FFT. @N should be power of 2, that is 2^(*) * === Output === * the @x contains the result of FFT algorithm, so the original data * in @x is destroyed, please store them before using FFT. */int fft(complex *x, int N){int i,j,l,ip;static int M = 0;static int le,le2;static FFT_TYPE sR,sI,tR,tI,uR,uI;M = (int)(log(N) / log(2));/* * bit reversal sorting */BitReverse(x,x,N,M);/* * For Loops */for (l=1; l<=M; l++) { /* loop for ceil{log2(N)} */le = (int)pow(2,l);le2 = (int)(le / 2);uR = 1;uI = 0;sR = cos(PI / le2);sI = -sin(PI / le2);for (j=1; j<=le2; j++) { /* loop for each sub DFT *///jm1 = j - 1;for (i=j-1; i<=N-1; i+=le) { /* loop for each butterfly */ip = i + le2;tR = x[ip].real * uR - x[ip].img * uI;tI = x[ip].real * uI + x[ip].img * uR;x[ip].real = x[i].real - tR;x[ip].img = x[i].img - tI;x[i].real += tR;x[i].img += tI;} /* Next i */tR = uR;uR = tR * sR - uI * sI;uI = tR * sI + uI *sR;} /* Next j */} /* Next l */return 0;}/* * Inverse FFT Algorithm * === Inputs === * x : complex numbers * N : nodes of FFT. @N should be power of 2, that is 2^(*) * === Output === * the @x contains the result of FFT algorithm, so the original data * in @x is destroyed, please store them before using FFT. */int ifft(complex *x, int N){int k = 0;for (k=0; k<=N-1; k++) {x[k].img = -x[k].img;}fft(x, N); /* using FFT */for (k=0; k<=N-1; k++) {x[k].real = x[k].real / N;x[k].img = -x[k].img / N;}return 0;}/* * Code below is an example of using FFT and IFFT. */#define SAMPLE_NODES (128)complex x[SAMPLE_NODES];int INPUT[SAMPLE_NODES];int OUTPUT[SAMPLE_NODES];static void MakeInput(){int i;for ( i=0;i<SAMPLE_NODES;i++ ){x[i].real = sin(PI*2*i/SAMPLE_NODES);x[i].img = 0.0f;INPUT[i]=sin(PI*2*i/SAMPLE_NODES)*1024;}}static void MakeOutput(){int i;for ( i=0;i<SAMPLE_NODES;i++ ){OUTPUT[i] = sqrt(x[i].real*x[i].real + x[i].img*x[i].img)*1024;}}void zx_fft(void){MakeInput();fft(x,128);MakeOutput();ifft(x,128);MakeOutput();}
程序在TMS320C6713上实验,主函数中调用zx_fft()函数即可。

FFT的采样点数为128,输入信号的实数域为正弦信号,虚数域为0,数据精度定义FFT_TYPE为float类型,MakeInput和MakeOutput函数分别用于产生输入数据INPUT和输出数据OUTPUT的函数,便于使用CCS 的Graph功能绘制波形图。这里调试时使用CCS v5中的Tools -> Graph功能得到下面的波形图(怎么用自己琢磨,不会的使用CCS 的Help)。

输入波形


输入信号的频域幅值表示


FFT运算结果


对FFT运算结果逆变换(IFFT)



如何检验运算结果是否正确呢?有几种方法:

(1)使用matlab验证,下面为相同情况的matlab图形验证代码

SAMPLE_NODES = 128;i = 1:SAMPLE_NODES;x = sin(pi*2*i / SAMPLE_NODES);subplot(2,2,1); plot(x);title('Inputs');axis([0 128 -1 1]);y = fft(x, SAMPLE_NODES);subplot(2,2,2); plot(abs(y));title('FFT');axis([0 128 0 80]);z = ifft(y, SAMPLE_NODES);subplot(2,2,3); plot(abs(z));title('IFFT');axis([0 128 0 1]);

(2)使用IFFT验证:输入信号的FFT获得的信号再IFFT,则的到的信号与原信号相同

可能大家发现输入信号上面的最后IFFT的信号似乎不同,这是因为FFT和IFFT存在精度截断误差(也叫数据截断噪声,意思就是说,我们使用的float数据类型数据位数有限,没法完全保留原始信号的信息)。因此,IFFT之后是复数(数据截断噪声引入了虚数域,只不过值很小),所以在绘图时使用了计算幅值的方法,

C代码中:

OUTPUT[i] = sqrt(x[i].real*x[i].real + x[i].img*x[i].img)*1024;

matlab代码中:

subplot(2,2,3); plot(abs(z));title('IFFT');

所以IFFT的结果将sin函数的负y轴数据翻到了正y轴。另外,在CCS v5的图形中我们将显示信号的幅度放大了1024倍便于观察,而matlab中没有放大。


版权声明:该文观点仅代表作者本人。处理文章:请发送邮件至 三1五14八八95#扣扣.com 举报,一经查实,本站将立刻删除。