美文网首页
FFT的计算—（使用FFTW库）

FFT的计算—（使用FFTW库）

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

    前言

    为了验证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

    • include目录


      image.png

    • lib 目录


      image.png

    • bin目录


      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结果完全一致。

    完整工程

    • C++
    #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库）

          本文链接：https://www.haomeiwen.com/subject/zrpjsqtx.html

          延伸阅读

          深度阅读

          • 您也可以注册成为美文阅读网的作者，发表您的原创作品、分享您的心情！

          栏目导航

          热点阅读

          关于我们|服务条款|联系我们|FFT的计算—（使用FFTW库）|投稿指南|网站地图|RSS订阅|排版工具|手机版

          提供经典美文摘抄,优美散文欣赏,现代诗歌精选,短篇小说,心情随笔,表白情书范文,故事会在线阅读欣赏

          Copyright © 2014-2023 Haomeiwen.com All Rights Reserved. 好美文阅读网 版权所有

          备案信息：桂公网安备 45052102000051号 · 桂ICP备13007215号-3

          本站所收录作品、热点评论等信息部分来源互联网，目的只是为了系统归纳学习和传递资讯

          所有作品版权归原创作者所有，与本站立场无关，如不慎侵犯了你的权益，请联系我们告知，我们将做删除处理！