FFT的计算—（使用FFTW库）

作者: 侠之大者_7d3f | 来源:发表于2019-07-06 11:32 被阅读0次

FFT的计算—（使用FFTW库）
在QT中使用FFTW库
FFTW使用指南
code fftw
Mac上交叉编译iOS静态库
Centos7编译安装fftw及其代码测试
FFT详解
快速傅里叶变换FFT(Fast Fourier Transfor
NumPy API（十四）——离散傅里叶变换
离散傅里叶的通用快速计算（fft的一般化）

前言

为了验证DSP(Digital Singal Processor 数字信号处理器)上的目标跟踪算法(KCF, Kernalized Correlation Filter 核化相关滤波)，但是DSP端调试相对麻烦，尤其是输入输出不像PC方便，因此打算将DSP端的算法在PC上测试一下(都是C语言的，移植也不难)。在DSP上算法使用了TI提供的dsplib, imagelib, mathlib库，这些库中的部分函数PC上没有直接对应的，比如dsplib中的FFT、IFFT计算等。

研一期间选过一门课，当时使用了Intel MKL计算FFT, MKL功能强大，专门为Intel处理器进行了优化，包含大量的数学计算函数，体积也相对大，下载还需要注册账户之类，比较麻烦。简单快速起见，本人决定使用FFTW库进行FFT/IFFT计算。

开发，测试环境

windows 10 64bit
pycharm, python 3.6
FFTW 64bit

FFTW安装配置

FFTW官网真心不错，提供源码，编译好的dll,方便多了，最重要的是还有在线文档。

下载FFTW

image.png

选择对应的平台，这里选择windwos, 64bit, precomplied

image.png

注意FFTW提供了windows预编译好的DLL, 如果要用VS调用FFTW, 还需要生成.lib
上面已经有详细的步骤生成lib, lib.exe是Visual Studio安装目录下自带的工具。

整理一下库：

image.png

将bin目录加入到系统环境变量，程序运行时候才可以加载到DLL。

测试代码


#include<iostream>
#include <fftw3.h>

using namespace std;

int main()
{
    int N = 32;
    fftw_complex *in, *out;
    fftw_plan p;
    
    in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
    out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);

    for (int i = 0; i < N; i++)
    {
        in[i][0] = i;
        in[i][1] = 0;
    }

    p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    
    fftw_execute(p); /* repeat as needed */

    for (int i = 0; i < N; i++)
    {
        cout << out[i][0] << " " << out[i][1] << endl;
    }
    
    fftw_destroy_plan(p);
    fftw_free(in); 
    fftw_free(out);

为了验证结果的正确性，笔者使用python scipy进行同样的fft

python测试代码

import scipy.fftpack as fftpack
import numpy as np


x = np.arange(0, 32)
y = fftpack.fft(x, 32)

y_real = np.real(y)
y_imag = np.imag(y)

print(y_real)
print(y_imag)

结果

C++ 调用FFTW的计算结果：

image.png

python scipy的计算结果：

image.png

FFT封装

直接调用FFTW库的函数不方便，因此对函数进行封装，使得API接口简单，这样在C++中也可以像Python一样直接fft()

一维FFT


/*
 计算1维FFT
 http://www.fftw.org/fftw3_doc/Complex-One_002dDimensional-DFTs.html#Complex-One_002dDimensional-DFTs
*/
void fft1d(const vector<complex<float>>& input, vector<complex<float>>& output, int n) {
    fftw_complex *in, *out;
    fftw_plan p;

    in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * input.size());
    out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);

    for (int i = 0; i < input.size(); i++)
    {
        in[i][0] = input.at(i).real();
        in[i][1] = input.at(i).imag();
    }

    p = fftw_plan_dft_1d(n, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

    fftw_execute(p); /* repeat as needed */

    for (int i = 0; i < n; i++)
    {
        output.push_back(complex<float>(out[i][0], out[i][1]));
    }

    fftw_destroy_plan(p);
    fftw_free(in);
    fftw_free(out);
}

二维FFT

/*
2维 FFT
@param input [in] 输入数据，复数类型， rows*cols
@param output [out] fft2d计算结果， same size as input
@param rows [in] 矩阵行数
@param cols [in] 矩阵列数

http://www.fftw.org/fftw3_doc/Complex-Multi_002dDimensional-DFTs.html#Complex-Multi_002dDimensional-DFTs
*/
void fft2d(const vector<complex<float>>& input, vector<complex<float>>& output, int rows, int cols) {

    fftw_complex *in, *out;
    fftw_plan p;

    in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * input.size());
    out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * rows*cols);

    for (int i = 0; i < input.size(); i++)
    {
        in[i][0] = input.at(i).real();
        in[i][1] = input.at(i).imag();
    }

    p = fftw_plan_dft_2d(rows, cols, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(p); /* repeat as needed */

    int size = rows * cols;
    for (int i = 0; i < size; i++)
    {
        output.push_back({ static_cast<float>(out[i][0]), static_cast<float>(out[i][1]) });
    }

    fftw_destroy_plan(p);
    fftw_free(in);
    fftw_free(out);

}

图像的二维FFT变换

利用上面封装好的C++ API函数 fft2d(), 图像的读取采用opencv, 为了验证结果，将fft的结果保存到文本，最后使用python加载，可视化。

C++测试函数

void test_fft2d() {

    int M = 512;
    int N = 512;

    vector<complex<float>> input, output;

    //使用opencv读取一张512x512图像
    Mat image = imread("./lena_gray.jpg", ImreadModes::IMREAD_GRAYSCALE);

    for (int i = 0; i < image.rows; i++)
    {
        for (int j = 0; j < image.cols; j++)
        {
            input.push_back({ static_cast<float>(image.at<unsigned char>(i,j)), 0 });
        }
    }
    
    // 调用FFT2d
    fft2d(input, output, 512, 512);

    //将FFT的结果保存到文本
    ofstream fileOut("./fft2d.dat", ios::out);
    for (auto& i:output)
    {
        fileOut << i.real() << " " << i.imag() << endl;
    }

    fileOut.close();

    cout << "FFT2d test finish!" << endl;

}

C++ fft2d可视化 (python代码)

from scipy.fftpack import fft, fftshift
import numpy as np
import matplotlib.pyplot as plt
import cv2


image = cv2.imread('./lena_gray.jpg',cv2.IMREAD_GRAYSCALE)

fft_data = np.loadtxt('./fft2d.dat')
print(fft_data.shape)

fft_result = np.empty(shape=[fft_data.shape[0]], dtype=np.complex)
for i in range(fft_data.shape[0]):
    data_real, data_imag = fft_data[i]
    fft_result[i] = np.complex(data_real, data_imag)

fft_result = fft_result.reshape([512, 512])

#
plt.figure()
plt.subplot(2, 3, 1)
plt.imshow(image, cmap=plt.gray())
plt.title('image')

plt.subplot(2, 3, 2)
plt.imshow(np.abs(fft_result), cmap=plt.gray())
plt.title('fft')

plt.subplot(2, 3, 3)
plt.imshow(np.abs(fftshift(fft_result)), cmap=plt.gray())
plt.title('fft-shift')

plt.subplot(2, 3, 4)
plt.imshow(np.log(1 + np.abs(fft_result)), cmap=plt.gray())
plt.title('fft-log')

plt.subplot(2, 3, 5)
plt.imshow(np.log(1 + np.abs(fftshift(fft_result))), cmap=plt.gray())
plt.title('fft-shiift-log')

plt.show()

image.png

为了对比，再次采用scipy中的fft2d直接对图像进行fft


from scipy.fftpack import fft, fftshift, fft2
import numpy as np
import cv2
import matplotlib.pyplot as plt


image = cv2.imread('./lena_gray.jpg', cv2.IMREAD_GRAYSCALE)

img_fft = fft2(image, shape=[512, 512])
print(img_fft.shape)

img_fftshift = fftshift(img_fft)

plt.figure()
plt.subplot(2, 3, 1)
plt.imshow(image, cmap=plt.gray())
plt.title('image')
plt.subplot(2, 3, 2)
plt.imshow(np.abs(img_fft), cmap=plt.gray())
plt.title('fft')
plt.subplot(2, 3, 3)
plt.imshow(np.abs(img_fftshift), cmap=plt.gray())
plt.subplot(2, 3, 4)
plt.imshow(np.log(1 + np.abs(img_fft)), cmap=plt.gray())
plt.title('fft-log')
plt.subplot(2, 3, 5)
plt.imshow(np.log(1 + np.abs(img_fftshift)), cmap=plt.gray())
plt.title('fftshift-log')
plt.show()

image.png

对比发现，C++ FFTW 与python scipy的FFT2d结果完全一致。

完整工程

#include<iostream>
#include <fftw3.h>
#include<vector>
#include<complex>
#include<fstream>
#include<opencv2/opencv.hpp>

using namespace std;
using namespace cv;

void fft1d(const vector<complex<float>>& input, vector<complex<float>>& output, int n);
void fft2d(const vector<complex<float>>& input, vector<complex<float>>& output, int rows, int cols);
void test_fft2d();
void test_fft1d();

int main()
{

    test_fft2d();

    system("pause");
}


void test_fft1d() {

    int N = 32;

    //fft1d测试
    vector<complex<float>> input, output;
    for (int i = 0; i < N; i++)
    {
        input.push_back({ static_cast<float>(i), 0 });
    }

    //fft1d
    fft1d(input, output, N);

    for (auto& i : output)
    {
        cout << i.real() << " " << i.imag() << endl;
    }
}


void test_fft2d() {

    int M = 512;
    int N = 512;

    vector<complex<float>> input, output;

    //使用opencv读取一张512x512图像
    Mat image = imread("./lena_gray.jpg", ImreadModes::IMREAD_GRAYSCALE);

    for (int i = 0; i < image.rows; i++)
    {
        for (int j = 0; j < image.cols; j++)
        {
            input.push_back({ static_cast<float>(image.at<unsigned char>(i,j)), 0 });
        }
    }
    
    // 调用FFT2d
    fft2d(input, output, 512, 512);

    //将FFT的结果保存到文本
    ofstream fileOut("./fft2d.dat", ios::out);
    for (auto& i:output)
    {
        fileOut << i.real() << " " << i.imag() << endl;
    }

    fileOut.close();

    cout << "FFT2d test finish!" << endl;

}

/*
 计算1维FFT
 http://www.fftw.org/fftw3_doc/Complex-One_002dDimensional-DFTs.html#Complex-One_002dDimensional-DFTs
*/
void fft1d(const vector<complex<float>>& input, vector<complex<float>>& output, int n) {
    fftw_complex *in, *out;
    fftw_plan p;

    in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * input.size());
    out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);

    for (int i = 0; i < input.size(); i++)
    {
        in[i][0] = input.at(i).real();
        in[i][1] = input.at(i).imag();
    }

    p = fftw_plan_dft_1d(n, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

    fftw_execute(p); /* repeat as needed */

    for (int i = 0; i < n; i++)
    {
        output.push_back(complex<float>(out[i][0], out[i][1]));
    }

    fftw_destroy_plan(p);
    fftw_free(in);
    fftw_free(out);
}

/*
2维 FFT
@param input [in] 输入数据，复数类型， rows*cols
@param output [out] fft2d计算结果， same size as input
@param rows [in] 矩阵行数
@param cols [in] 矩阵列数

http://www.fftw.org/fftw3_doc/Complex-Multi_002dDimensional-DFTs.html#Complex-Multi_002dDimensional-DFTs
*/
void fft2d(const vector<complex<float>>& input, vector<complex<float>>& output, int rows, int cols) {

    fftw_complex *in, *out;
    fftw_plan p;

    in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * input.size());
    out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * rows*cols);

    for (int i = 0; i < input.size(); i++)
    {
        in[i][0] = input.at(i).real();
        in[i][1] = input.at(i).imag();
    }

    p = fftw_plan_dft_2d(rows, cols, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(p); /* repeat as needed */

    int size = rows * cols;
    for (int i = 0; i < size; i++)
    {
        output.push_back({ static_cast<float>(out[i][0]), static_cast<float>(out[i][1]) });
    }

    fftw_destroy_plan(p);
    fftw_free(in);
    fftw_free(out);

}

python

from scipy.fftpack import fft, fftshift, fft2
import numpy as np
import cv2
import matplotlib.pyplot as plt


image = cv2.imread('./lena_gray.jpg', cv2.IMREAD_GRAYSCALE)

img_fft = fft2(image, shape=[512, 512])
print(img_fft.shape)

img_fftshift = fftshift(img_fft)

plt.figure()
plt.subplot(2, 3, 1)
plt.imshow(image, cmap=plt.gray())
plt.title('image')
plt.subplot(2, 3, 2)
plt.imshow(np.abs(img_fft), cmap=plt.gray())
plt.title('fft')
plt.subplot(2, 3, 3)
plt.imshow(np.abs(img_fftshift), cmap=plt.gray())
plt.subplot(2, 3, 4)
plt.imshow(np.log(1 + np.abs(img_fft)), cmap=plt.gray())
plt.title('fft-log')
plt.subplot(2, 3, 5)
plt.imshow(np.log(1 + np.abs(img_fftshift)), cmap=plt.gray())
plt.title('fftshift-log')
plt.show()


from scipy.fftpack import fft, fftshift
import numpy as np
import matplotlib.pyplot as plt
import cv2


image = cv2.imread('./lena_gray.jpg',cv2.IMREAD_GRAYSCALE)

fft_data = np.loadtxt('./fft2d.dat')
print(fft_data.shape)

fft_result = np.empty(shape=[fft_data.shape[0]], dtype=np.complex)
for i in range(fft_data.shape[0]):
    data_real, data_imag = fft_data[i]
    fft_result[i] = np.complex(data_real, data_imag)

fft_result = fft_result.reshape([512, 512])

#
plt.figure()
plt.subplot(2, 3, 1)
plt.imshow(image, cmap=plt.gray())
plt.title('image')

plt.subplot(2, 3, 2)
plt.imshow(np.abs(fft_result), cmap=plt.gray())
plt.title('fft')

plt.subplot(2, 3, 3)
plt.imshow(np.abs(fftshift(fft_result)), cmap=plt.gray())
plt.title('fft-shift')

plt.subplot(2, 3, 4)
plt.imshow(np.log(1 + np.abs(fft_result)), cmap=plt.gray())
plt.title('fft-log')

plt.subplot(2, 3, 5)
plt.imshow(np.log(1 + np.abs(fftshift(fft_result))), cmap=plt.gray())
plt.title('fft-shiift-log')

plt.show()

参考

https://www.cnblogs.com/youmuchen/p/8361713.html

FFT的计算—（使用FFTW库）
前言为了验证DSP(Digital Singal Processor 数字信号处理器)上的目标跟踪算法(KCF,...
在QT中使用FFTW库
本文主要介绍在windows 64位系统中，建立QT MSVC2015 x64项目时如何使用FFTW库。32位的请...
FFTW使用指南
FFTW使用指南编译说明引入fftw3.h头文件在Unix系统中需要加上"-lfftw3 -lm"编译选项。需...
code fftw
fftw在linux下的安装 FFTW（the Fastest Fourier Transform in the ...
Mac上交叉编译iOS静态库
本文主要以fftw为例子，进行iOS静态库的编译。环境要求：运行脚本如下：
Centos7编译安装fftw及其代码测试
FFTW 软件名称：the Faster Fourier Transform in the West FFTW 软...
FFT详解
FFT详解频域模值*2/N=时域幅度 fft/ifft前后的数据会在能量上差sqrt(n); 可以使用rms()...
快速傅里叶变换FFT(Fast Fourier Transfor
快速傅里叶变换(Fast Fourier Transform,FFT)用来计算离散傅里叶变换（Discrete F...
NumPy API（十四）——离散傅里叶变换
离散傅里叶变换标准的 FFTs fft(a[, n, axis, norm]) 计算一维离散傅立叶变换。 iff...
离散傅里叶的通用快速计算（fft的一般化）
之前有篇文章实现了合数的fft，和2的幂次方的fft，这里主要是关于更加一般形式的快速计算时间隔了很久，才出了第二...

FFT的计算—（使用FFTW库）

前言

开发，测试环境

FFTW安装配置

测试代码

结果

FFT封装

一维FFT

二维FFT

图像的二维FFT变换

完整工程

参考

相关文章

FFT的计算—（使用FFTW库）

在QT中使用FFTW库

FFTW使用指南

code fftw

Mac上交叉编译iOS静态库

Centos7编译安装fftw及其代码测试

FFT详解

快速傅里叶变换FFT(Fast Fourier Transfor

NumPy API（十四）——离散傅里叶变换

离散傅里叶的通用快速计算（fft的一般化）

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读