windows下使⽤fftw进⾏傅⾥叶变换及其编程实例傅⾥叶变换应该是上⼤⼆的时候《信号与系统》课上学过,上研后在《数字信号处理》课上⼜学了⼀遍。当初⼀直在想傅⾥叶变换到底有什么⽤呢?什么时候能⽤上呢?时间如梭,没想到毕业四年后,⼀个⼩项⽬要⽤到傅⾥叶变换,喜⼤普奔啊,当初晦涩的概念、眼晕的公式,终于没⽩学。是的,其实任何知识都不是⽩学的,即使⼯作中⼀直⽤不到傅⾥叶变换,⾄少思维得到了锻炼,都是有益的。在应⽤傅⾥叶变换过程中,可以按照公式⾃⼰编程实现,也可以使⽤已有的开源傅⾥叶变换函数库。经过查阅资料,本⼈选⽤fftw函数库来进⾏傅⾥叶变换,因为 fftw 官⽹上提到说该函数库性能⾮常出众,甚⾄可以和同类型的商业软件相竞争。fftw官⽹地址为:。
⼀、fftw下载
由于本⼈是在windows平台开发,于是到下⾯链接:下载已编译版本。如下图所⽰,该图信息量很⼤,既有32/64两个版本的下载链接,同时还有⽣成lib⽂件命令。即32位系统命令为: lib /def:libfftw3-3.def,⽽64位系统为: lib /machine:x64 /def:libfftw3-3.def。
下载解压后,会发现只有dll和头⽂件,在没有lib⽂件时,使⽤fftw库时会⽐较繁琐。于是在应⽤fftw编程前还是动⼿先⽣成lib⽂件。⽅法很简单,⾸先打开VS命令提⽰⼯具。本⼈使⽤的开发⼯具为VS2005,如下图所⽰:
打开后,cd到已下载的fftw⽂件⽬录,然后执⾏前⾯提到的lib⽣成命令,如下图所⽰:
命令执⾏完毕后,即可在当前⽬录内看到新⽣成的lib⽂件,我这边分别⽣成了double、float、long double三个版本,实际使⽤时可根据项⽬需要,选择相应精度的版本。⽣成的lib库⽂件如下图所⽰:
⼆、简单编程实例
有了头⽂件、lib/dll⽂件,那就什么都有了。新建⼀个⼯程,马上⽤起来吧。下⾯为⼀个⾮常简单的编程实例,主要完成实数傅⾥叶变换,输⼊为正弦曲线,输出为傅⾥叶变换的幅度谱。
#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
vs编程软件#include "include/fftw3.h"
#pragma comment(lib, "libfftw3-3.lib") // double版本
// #pragma comment(lib, "libfftw3f-3.lib")// float版本
// #pragma comment(lib, "libfftw3l-3.lib")// long double版本
#define PI 3.1415926
int main()
{
int len = 8;
double *in = NULL;
// 如果要使⽤float版本,需先引⽤float版本的lib库,然后在fftw后⾯加上f后缀即可.
fftw_complex *out = NULL;// fftwf_complex --> 即为float版本
fftw_plan p;
in = (double *)fftw_malloc(sizeof(double) * len);
out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * len);
double dx = 1.0 / len;
// 输⼊纯实数
for (int i = 0; i < len; i++)
{
in[i] = sin(2*PI * dx*i) + sin(4*PI * dx*i);
printf("%.2f ", in[i]);
}
printf("\n\n");
/
/ 傅⾥叶变换
p = fftw_plan_dft_r2c_1d(len, in, out, FFTW_ESTIMATE);
fftw_execute(p);
// 输出幅度谱
for (int i = 0; i < len; i++)
{
float len = sqrt(out[i][0]*out[i][0] + out[i][1]*out[i][1]);
printf("%.2f ", len);
}
printf("\n");
// 释放资源
fftw_destroy_plan(p);
fftw_free(in);
fftw_free(out);
system("pause");
return 0;
}
运⾏结果如下图所⽰,从结果来看,在频率为1及2两个位置有幅度输出,这与输⼊的频率为1及2的两个叠加正弦波吻合,结果正确。(8⽉23⽇补充:由于实数DFT变换的输出序列具有共轭对称性,因此fftw_plan_dft_r2c_1d()函数返回的序列实际长度为len /2 + 1,因此上述代码在为out分配内存时仅需分配len / 2 + 1即可。)
版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。
发表评论