FFTW基本概念与安装使用

本文主要介绍FFTW库的基本概念、编译安装和使用方法。

1 基本概念

FFTW (Fastest Fourier Transform in the West)是C语言的一个子程序库，用于计算一维或多维离散傅里叶变换(Discrete Fourier transform, DFT)，输入可以是任意长度的实数或复数数组。最新的版本为3.3.10，所以也记为FFTW3。
FFTW的计算流程包括规划(planning)和执行(execution)两个阶段。在规划阶段，规划器(planner)会根据用户输入的复杂问题递归的分解为更简单的子问题，当子问题足够简单时，便调用预定义的代码直接解决。由此形成诸多解决问题的方案(plans)，并从其中选择速度最快的一种，方案确定后便可进入执行阶段。在执行阶段，用户按照数据结构输入待变换数组，然后按照方案逐步执行，最后输入变换后的数组。
对于给定傅里叶变换问题，FFTW规划器能够给出解决问题的一系列方案，并从中选择最快的一种。FFTW能够解决多维的傅里叶变换问题，定义二维傅里叶变换的输入数据数组为n0×n1，三维傅里叶变换的输入数据数组为n0×n1×n2。数组中数据的顺序为row-major顺序，即最后一个维度的索引变化最快(n0×n1表示n0个n1点的序列，n1连续变化)。FFTW有4中工作模式：

• 单线程单内存：使用一个线程、且数据存储到一个内存中。
• 多线程共享内存：为提升计算速度，可以使用多线程并行进行一维或多维的傅里叶变换，所有线程使用同一共享内存。
• 多核分布式内存：当变换的矩阵非常大以至于单处理器的内存容量不足时，分布式内存并行处理架构则变得非常实用。在该架构下，每个处理器具备独立的内存，成千上万个处理器构成一个集群。每个处理器对应的内存仅存储待变换数组的一部分数据，从内存中读取相应的数据进行傅里叶变换。因此在进行多维傅里叶变换时，需要在分布式内存之间进行数据交互，数据交互方式使用MPI (Message-Passing Interface)。
• 多处理器多线程：FFTW也支持多处理器和多线程同时使用，例如具备4个支持共享内存的处理器节点，可以在每个节点内使用多线程进行并行计算，在节点间使用MPI进行数据交互。由于计算流程和数据交互更加复杂，需要进行详细的设计。

2 编译安装

从FFTW官网下载FFTW3.3.10，依次执行下述命令对软件进行配置、编译和安装。

• 配置：./configure --enable-threads
• 编译：make
• 安装：make install
• 安装依赖
  • Ubuntu：sudo apt-get install libfftw3-dev libfftw3-doc libfftw3-double3
  • CentOS：yum install fftw-devel
• 卸载：make distclean

根据第三节说明编写程序，程序命名为fftw.cpp，执行下述编译语句进行编译和运行，-lfftw3表示使用fftw3的库进行编译，-lfftw3_threads表示使用多线程的库。

• 编译
  • 单线程编译：g++ -o fftw fftw.cpp -lfftw3
  • 多线程编译：g++ -o fftw fftw.cpp -lfftw3 -lfftw3_threads -lm
• 运行：./fftw

3 使用实例

3.1 单线程

• 一维傅里叶变换

#include <complex>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <ctime>
#include <chrono>
#include <unistd.h>
#include "fftw3.h"using namespace std;int main() {const int len = 8; //序列长度fftw_complex *in, *out; //定义输入输出变量fftw_plan plan; //定义方案in  = (fftw_complex*) fftw_malloc (sizeof(fftw_complex) * len); //分配输入空间out = (fftw_complex*) fftw_malloc (sizeof(fftw_complex) * len); //分配输出空间// 初始化输入信号for (int i = 0; i < len; i++) {in[i][0] = 35.24 * i; //实部in[i][1] = -27.62 * i; //虚部std::cout << in[i][0] << " " << in[i][1] << std::endl;}// 定义fftw_planplan = fftw_plan_dft_1d(len, in, out, FFTW_FORWARD, FFTW_ESTIMATE);//plan = fftw_plan_dft_1d(len, in, out, FFTW_FORWARD, FFTW_MEASURE);// 执行fftauto start = std::chrono::high_resolution_clock::now();fftw_execute(plan);//usleep(1);auto end = std::chrono::high_resolution_clock::now();auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();std::cout << "程序运行时间: " << duration << " 微秒" << std::endl;// 输出结果ofstream fout("out.txt");int w = 36;for (int i = 0; i < len; i++) {if (out[i][1] < 0) {fout << out[i][0] << out[i][1] << ' ';}else {fout << out[i][0] << '+' << out[i][1] << ' ';}}// 释放内存fftw_cleanup();fftw_free(in);fftw_free(out);return 0;
}

• 二维傅里叶变换

#include <complex>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <ctime>
#include <chrono>
#include <unistd.h>
#include "fftw3.h"using namespace std;#define ROWS 64
#define COLS 64int main() {fftw_complex *in, *out; //定义输入输出变量fftw_plan plan; //定义方案in  = (fftw_complex*) fftw_malloc (sizeof(fftw_complex)*ROWS*COLS); //分配输入空间out = (fftw_complex*) fftw_malloc (sizeof(fftw_complex)*ROWS*COLS); //分配输出空间// 初始化输入信号for (int i = 0; i < ROWS; i++) {for (int j = 0; j < COLS; j++) {in[i*COLS+j][0] = 35.24 * i; //实部in[i*COLS+j][1] = -27.62 * j; //虚部std::cout << in[i*COLS+j][0] << " " << in[i*COLS+j][1] << std::endl;}}// 定义fftw_planplan = fftw_plan_dft_2d(ROWS, COLS, in, out, FFTW_FORWARD, FFTW_ESTIMATE);//plan = fftw_plan_dft_1d(len, in, out, FFTW_FORWARD, FFTW_MEASURE);// 执行fftauto start = std::chrono::high_resolution_clock::now();fftw_execute(plan);//usleep(1);auto end = std::chrono::high_resolution_clock::now();auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();std::cout << "程序运行时间: " << duration << " 微秒" << std::endl;// 输出结果ofstream fout("results_cpp");int w = 36;for (int i = 0; i < ROWS; i++) {for (int j = 0; j < COLS; j++) {if (out[i*COLS+j][1] < 0) {fout << out[i*COLS+j][0] << out[i*COLS+j][1] << ' ';}else {fout << out[i*COLS+j][0] << '+' << out[i*COLS+j][1] << ' ';}if (j == COLS-1) {fout << endl;}}}// 释放内存fftw_cleanup();fftw_free(in);fftw_free(out);return 0;
}

• 三维傅里叶变换

#include <complex>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <ctime>    //计时
#include <chrono>   //计时
#include <unistd.h> //usleep函数
#include <string.h>
#include <stdio.h>
#include "fftw3.h"using namespace std;#define HIGH 64
#define ROWS 64
#define COLS 64int main() {fftw_complex *in, *out; //定义输入输出变量fftw_plan plan; //定义方案int i, j, k;in  = (fftw_complex*) fftw_malloc (sizeof(fftw_complex)*HIGH*ROWS*COLS); //分配输入空间out = (fftw_complex*) fftw_malloc (sizeof(fftw_complex)*HIGH*ROWS*COLS); //分配输出空间// 初始化输入信号for (k = 0; k < HIGH; k++) {for (j = 0; j < ROWS; j++) {for (i = 0; i < COLS; i++) {in[k*ROWS*COLS+j*COLS+i][0] = k*215.42 + 35.24*i; //实部in[k*ROWS*COLS+j*COLS+i][1] = k*215.42 - 27.62*j; //虚部//std::cout << in[k*ROWS*COLS+j*COLS+i][0] << " " << in[k*ROWS*COLS+j*COLS+i][1] << std::endl;}}}/*ifstream fin("data");string line;  //读取文件的每一行string value; //每一行的数据int cnt;int pos; //分割符位置int len; //字符串长度k = j = cnt = 0;while (getline(fin, line)) {j = cnt / COLS;if (cnt >= COLS * ROWS) {k++;cnt = 0;} else {cnt = cnt + COLS;}len = line.size(); //字符串长度if (len == 0) {continue; //跳过空行}for (i = 0; i < COLS; i++) {len = line.size();pos = line.find(' '); //查找字符在字符串中的位置value = line.substr(0,pos);//std::cout << value << ' ';line = line.substr(pos+1,len); //剩余字符串//std::cout << line << endl;//if (i == COLS-1) {//  std::cout << endl;//}in[k*ROWS*COLS+j*COLS+i][0] = stof(value); //实部in[k*ROWS*COLS+j*COLS+i][1] = 0; //虚部//std::cout << in[k*ROWS*COLS+j*COLS+i][0] << " " << in[k*ROWS*COLS+j*COLS+i][1] << std::endl;}}*/// 定义fftw_planplan = fftw_plan_dft_3d(HIGH, ROWS, COLS, in, out, FFTW_FORWARD, FFTW_ESTIMATE);// 执行fftauto start = std::chrono::high_resolution_clock::now();fftw_execute(plan);//usleep(1);auto end = std::chrono::high_resolution_clock::now();auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();std::cout << "程序运行时间: " << duration << " 微秒" << std::endl;// 输出结果ofstream fout("results_cpp");for (k = 0; k < HIGH; k++) {for (j = 0; j < ROWS; j++) {for (i = 0; i < COLS; i++) {if (out[k*ROWS*COLS+j*COLS+i][1] < 0) {fout << out[k*ROWS*COLS+j*COLS+i][0] << out[k*ROWS*COLS+j*COLS+i][1] << ' ';} else {fout << out[k*ROWS*COLS+j*COLS+i][0] << '+' << out[k*ROWS*COLS+j*COLS+i][1] << ' ';}if (i == COLS-1) {fout << endl;}if (j==ROWS-1 && i==COLS-1) {fout << endl;}}}}// 释放内存fftw_cleanup();fftw_free(in);fftw_free(out);return 0;
}

3.2 多线程

• 三维傅里叶变换

#include <complex>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <ctime>    //计时
#include <chrono>   //计时
#include <unistd.h> //usleep函数
#include <string.h>
#include <stdio.h>
#include "fftw3.h"using namespace std;#define HIGH 64
#define ROWS 64
#define COLS 64int main() {fftw_complex *in, *out; //定义输入输出变量fftw_plan plan; //定义方案int nthreads; //使用的线程数int i, j, k;in  = (fftw_complex*) fftw_malloc (sizeof(fftw_complex)*HIGH*ROWS*COLS); //分配输入空间out = (fftw_complex*) fftw_malloc (sizeof(fftw_complex)*HIGH*ROWS*COLS); //分配输出空间nthreads = 12;// 初始化输入信号for (k = 0; k < HIGH; k++) {for (j = 0; j < ROWS; j++) {for (i = 0; i < COLS; i++) {in[k*ROWS*COLS+j*COLS+i][0] = k*215.42 + 35.24*i; //实部in[k*ROWS*COLS+j*COLS+i][1] = k*215.42 - 27.62*j; //虚部//std::cout << in[k*ROWS*COLS+j*COLS+i][0] << " " << in[k*ROWS*COLS+j*COLS+i][1] << std::endl;}}}/*// 从文件输入数据 ifstream fin("data");string line;  //读取文件的每一行string value; //每一行的数据int cnt;int pos; //分割符位置int len; //字符串长度k = j = cnt = 0;while (getline(fin, line)) {j = cnt / COLS;if (cnt >= COLS * ROWS) {k++;cnt = 0;} else {cnt = cnt + COLS;}len = line.size(); //字符串长度if (len == 0) {continue; //跳过空行}for (i = 0; i < COLS; i++) {len = line.size();pos = line.find(' '); //查找字符在字符串中的位置value = line.substr(0,pos);//std::cout << value << ' ';line = line.substr(pos+1,len); //剩余字符串//std::cout << line << endl;//if (i == COLS-1) {//  std::cout << endl;//}in[k*ROWS*COLS+j*COLS+i][0] = stof(value); //实部in[k*ROWS*COLS+j*COLS+i][1] = 0; //虚部//std::cout << in[k*ROWS*COLS+j*COLS+i][0] << " " << in[k*ROWS*COLS+j*COLS+i][1] << std::endl;}}*/fftw_init_threads();fftw_plan_with_nthreads(nthreads);plan = fftw_plan_dft_3d(HIGH, ROWS, COLS, in, out, FFTW_FORWARD, FFTW_ESTIMATE); // 定义fftw_planauto start = std::chrono::high_resolution_clock::now();fftw_execute(plan); // 执行fft//usleep(1);auto end = std::chrono::high_resolution_clock::now();auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();std::cout << "程序运行时间: " << duration << " 微秒" << std::endl;// 输出结果ofstream fout("results_cpp");for (k = 0; k < HIGH; k++) {for (j = 0; j < ROWS; j++) {for (i = 0; i < COLS; i++) {if (out[k*ROWS*COLS+j*COLS+i][1] < 0) {fout << out[k*ROWS*COLS+j*COLS+i][0] << out[k*ROWS*COLS+j*COLS+i][1] << ' ';} else {fout << out[k*ROWS*COLS+j*COLS+i][0] << '+' << out[k*ROWS*COLS+j*COLS+i][1] << ' ';}if (i == COLS-1) {fout << endl;}if (j==ROWS-1 && i==COLS-1) {fout << endl;}}}}// 释放内存fftw_cleanup();fftw_free(in);fftw_free(out);return 0;
}