快速傅里叶变换很难实现吗?其实有它就足够了

相信内行一看标题就知道我要说谁了!没错,就是 FFTW 。这是一个求解快速傅里叶变换的开源库,就连大佬级别的Matlab也是用FFTW家的库。当然了,也有人知道这个库的存在,可能苦于这个库还不足够傻瓜,而长时间没有亲自体验。本文的目的就是一步一步道来FFTW的使用过程。linux和mac系统下的使用太简单了,问后稍微提一下就行。主要讲解windows系统下的使用,我这里使用的环境是:windows 10 + Visual Studio 2017 Community版本

Windows

下载

获取FFTW库有两种方式,一种是从 源码编译 ;一种是直接使用 官方发布的dll,有 x63win32 两种版本可供选择。我推荐 第二种方式 。我测试的是x64版本,直接 下载 后解压即可。解压后的文件列表如下图所示:

FFTW-3.3.5-dll64.zip解压后的文件列表

绿色方框内的 .dll 文件和 .def 文件是两个关键文件。还有类似的两组,分别后面带有 fl 应该表示数值精度的区别,同样的用法。

生成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

首先将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;
}

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 205,386评论 6 479
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 87,939评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 151,851评论 0 341
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,953评论 1 278
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,971评论 5 369
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,784评论 1 283
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,126评论 3 399
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,765评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 43,148评论 1 300
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,744评论 2 323
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,858评论 1 333
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,479评论 4 322
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,080评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,053评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,278评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,245评论 2 352
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,590评论 2 343

推荐阅读更多精彩内容