当前位置：首页 > article >正文

FFTW基本概念与安装使用

article 2024/12/26 12:05:21

FFTW基本概念与安装使用

1 基本概念
2 编译安装
3 使用实例
- 3.1 单线程
- 3.2 多线程

本文主要介绍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_plan
  plan = fftw_plan_dft_1d(len, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
  //plan = fftw_plan_dft_1d(len, in, out, FFTW_FORWARD, FFTW_MEASURE);

  // 执行fft
  auto 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 64

int 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_plan
  plan = fftw_plan_dft_2d(ROWS, COLS, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
  //plan = fftw_plan_dft_1d(len, in, out, FFTW_FORWARD, FFTW_MEASURE);

  // 执行fft
  auto 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 64

int 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_plan
  plan = fftw_plan_dft_3d(HIGH, ROWS, COLS, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

  // 执行fft
  auto 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 64

int 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_plan

  auto 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;
}

查看全文

http://www.kler.cn/a/452114.html