相信内行一看标题就知道我要说谁了!没错,就是 FFTW 。这是一个求解快速傅里叶变换的开源库,就连大佬级别的Matlab也是用FFTW家的库。当然了,也有人知道这个库的存在,可能苦于这个库还不足够傻瓜,而长时间没有亲自体验。本文的目的就是一步一步道来FFTW的使用过程。linux和mac系统下的使用太简单了,问后稍微提一下就行。主要讲解windows系统下的使用,我这里使用的环境是:
windows 10
+Visual Studio 2017 Community版本
Windows
下载
获取FFTW库有两种方式,一种是从 源码编译 ;一种是直接使用 官方发布的dll,有 x63 和 win32 两种版本可供选择。我推荐 第二种方式 。我测试的是x64版本,直接 下载 后解压即可。解压后的文件列表如下图所示:
绿色方框内的 .dll 文件和 .def 文件是两个关键文件。还有类似的两组,分别后面带有 f
和 l
应该表示数值精度的区别,同样的用法。
生成lib
一般在VS中编写C/C++程序调用dll时候,同时也需要lib文件作为连接库文件。虽然FFTW下载文件里面没有对应的lib文件,但是没关系,可以使用一个工具生成lib文件。
这个工具就是vs2017自带的一个程序 lib.exe
,位于 C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Tools\MSVC\14.15.26726\bin\Hostx64\x64
目录(自己的系统下也可以按照类似的路径去寻找)。
然后在命令行中切换到此目录下,运行命令。
命令格式为:lib /out:保存lib的路径\fftw3-3.lib /machine:x64或者win32 /def:.def文件的路径\fftw3-3.def
运行命令: lib /out:D:\SoftwareProject\fftw-3.3.5-dll64\fftw3-3.lib /machine:x64 /def:D:\SoftwareProject\fftw-3.3.5-dll64\libfftw3-3.def
。
或者直接在随便一个目录下新建一个.bat文件(比如createlib.bat),里面写入如下代码,然后双击运行即可
C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Tools\MSVC\14.15.26726\bin\Hostx64\x64\lib /out:D:\SoftwareProject\fftw-3.3.5-dll64\fftw3-3.lib /machine:x64 /def:D:\SoftwareProject\fftw-3.3.5-dll64\libfftw3-3.def
然后在代码中填入的lib生成路径(D:\SoftwareProject\fftw-3.3.5-dll64)下看到新生成的两个文件: fftw3-3.lib和fftw3-3.exp。现在就可以正式编写代码使用fftw进行快速傅里叶变换啦。
使用fftw
这里举个一维傅里叶变换的例子,二维和三维的用法与一维的类似,细节可以参考 官方文档 。
问题描述和模拟数据生成
(1) 问题: 提取一个叠加正弦信号的频率。
(2)模拟信号: 用代码生成频率为30 Hz和90Hz的两个正弦信号的叠加作为模拟信号,可用下面这段matlab代码生成,并保存到文件 signal_30_90.txt。
% create a 1-D signal with frequancy of 30Hz and 90Hz
f1=30; f2=90;
x=linspace(-3*pi,3*pi,2000);
y=sin(2*pi*f1*x)+sin(2*pi*f2*x);
figure(1);
plot(x,y);
fpout=fopen('signal_30_90.txt','w');
for i=1:length(x)
fprintf(fpout,'%.5E %.5E\n',x(i),y(i));
end
fclose(fpout);
编写C++代码使用fftw进行频谱分析
前奏
(1)添加fftw头文件
有两种方式添加fftw的头文件到自己的项目中:第一种方式是直接将解压文件中的fftw3.h拷贝到项目目录下;第二种方式是在vs2017的项目属性中添加包含目录中加入fftw3.h所在的路径。然后就可以在自己的代码中用 #include "fftw3.h
包含头文件了,不然编译器会提示错误:找不到头文件。推荐第一种方式,因为毕竟文件很小很少,而且这种方式方便移植程序给别人 。
(2)添加链接库
与(1)类似,将libfftw3-3.lib拷贝到项目目录下,然后在vs2017项目属性中的连接库下面添加libfftw3-3.lib,操作过程如下图所示:
在代码中使用fftw
首先将fftw通过 #include "fftw3.h"
将fftw的头文件加入代码中,然后就可以用fftw中的函数进行正、逆傅里叶变换。
下面代码功能是从文件 signal_30_90.txt
中读取信号,然后对其进行傅里叶变化求出频谱并保存到文件 Spectrum_Singal.txt
中,结果见下图
// TestFFT.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
//
#define PI 3.141592653
#include <iostream>
#include <fstream>
using namespace std;
// 将fftw的使用进行一定的封装,🙆使用过程更简便
void FFT1d(vector<double> invector_r, fftw_complex* out_c);
//测试一维傅里叶变换
int testFFT1d();
int main()
{
//测试一维傅里叶变换:求解一个正弦叠加信号的频谱
testFFT1d();
return 0;
}
void FFT1d(vector<double> invector_r, fftw_complex* out_c)
{
//1. 申明数据输入输出数组
fftw_complex *in;
//2. 申明fft执行变量
fftw_plan p;//,q;
//3. 给出数据个数
int N = invector_r.size();
//4. 根据数据个数开辟动态数组空间
in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)* N);
//out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)* N);
//5. 给输入数据数组赋值
for (int i = 0; i < N; i++)
{
in[i][0] = invector_r[i];
in[i][1] = 0.0;
//printf("%6.2f \n", in[i][0]);
}
//6. 调用fft执行函数
p = fftw_plan_dft_1d(N, in, out_c, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(p); /* repeat as needed*/
//8. 最后销毁fft相关的变量
fftw_destroy_plan(p);
//fftw_destroy_plan(q);
fftw_free(in);
}
int testFFT1d()
{
string filename = "data\\signal_30_90.txt";
ifstream fin = ifstream(filename);
vector<double> t;
vector<double> amplitude;
vector<double> fvector;
vector<double> SpectrumVector;
double t0, amp0;
if (!fin)
{
cout << "打开文件失败: " << endl;
return 0;
}
while (!fin.eof())
{
fin >> t0 >> amp0;
t.push_back(t0);
amplitude.push_back(amp0);
}
cout << t[0] << " " << amplitude[0] << " " << t.size() << endl;
//2. Calculate FFT
fftw_complex* fftout=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*t.size());
FFT1d(amplitude, fftout);
//3. calculate spectrum
GetSpectrum1d(t.size(), t[1]-t[0],fftout, fvector, SpectrumVector);
//write spectrum to file
ofstream fout = ofstream("data\\Spectrum_Singal.txt");
for (int i = 0; i < SpectrumVector.size(); i++)
{
fout<<fvector[i]<<"\t" << SpectrumVector[i] << endl;
}
fout.close();
//free pointer
fftw_free(fftout);
return 0;
}