1.使用Python包实现

import numpy as np 
import matplotlib.pyplot as plt
a=np.array([0, 2, 4, 1,6, 1, 3, 2,5]).reshape(3,3)
f=np.fft.fft2(a)
fshift=np.fft.fftshift(f)
mag=20*np.log(np.abs(fshift))
plt.axis("off")
plt.imshow(mag)
plt.show()

其中：

f:
array([[24. +0.00000000e+00j, -6. +4.44089210e-16j, -6. -4.44089210e-16j],[-3. +1.73205081e+00j, -7.5+4.33012702e+00j,  4.5-8.66025404e-01j],[-3. -1.73205081e+00j,  4.5+8.66025404e-01j, -7.5-4.33012702e+00j]])fshift:
array([[-7.5-4.33012702e+00j, -3. -1.73205081e+00j,  4.5+8.66025404e-01j],[-6. -4.44089210e-16j, 24. +0.00000000e+00j, -6. +4.44089210e-16j],[ 4.5-8.66025404e-01j, -3. +1.73205081e+00j, -7.5+4.33012702e+00j]])

1.1 fftshift在numpy中的实现

通过官方代码，可以看出fftshift是通过np.roll实现的。

def fftshift(x, axes=None):x = asarray(x)if axes is None:axes = tuple(range(x.ndim))shift = [dim // 2 for dim in x.shape]return roll(x, shift, axes)

对于3*3的二维矩阵，对应np.roll(x,1,(0,1)。

np.roll(f, 1,0)
array([[-3. -1.73205081e+00j,  4.5+8.66025404e-01j, -7.5-4.33012702e+00j],[24. +0.00000000e+00j, -6. +4.44089210e-16j, -6. -4.44089210e-16j],[-3. +1.73205081e+00j, -7.5+4.33012702e+00j,  4.5-8.66025404e-01j]])np.roll(np.roll(f, 1,0),1,1)
array([[-7.5-4.33012702e+00j, -3. -1.73205081e+00j,  4.5+8.66025404e-01j],[-6. -4.44089210e-16j, 24. +0.00000000e+00j, -6. +4.44089210e-16j],[ 4.5-8.66025404e-01j, -3. +1.73205081e+00j, -7.5+4.33012702e+00j]])

与np.fft.fftshift(f)结果一致。

1.2 平移后的幅度谱

2.使用c++实现之1

#include <iostream>
using namespace std;
#include"Eigen/Dense"
using namespace Eigen;
#include "fftw3.h"int main()
{MatrixXd a(3, 3), out(a.rows(), a.cols());MatrixXcd FTa(a.rows() / 2 + 1, a.cols());a << 0, 2, 4, 1,6, 1, 3, 2,5;fftw_plan P;P = fftw_plan_dft_r2c_2d(a.cols(), a.rows(), a.data(), (fftw_complex*)FTa.data(), FFTW_ESTIMATE);fftw_execute(P);cout << "dft" << endl;cout << FTa << endl;cout << endl;P = fftw_plan_dft_c2r_2d(a.cols(), a.rows(), (fftw_complex*)FTa.data(), out.data(), FFTW_ESTIMATE);fftw_execute(P);cout << "idft" << endl;out = out / (a.cols() * a.rows());cout << out << endl;return 0;
}

结果如下：

dft(24,0)          (-6,0)          (-6,0)(-3,1.73205)  (-7.5,4.33013) (4.5,-0.866025)idft
0 2 4
1 6 1
3 2 5

2.1 FFTW库安装

这里用到了FFTW c++库，具体编译及调用可参考Windows下FFTW_2.1.5的编译及使用。

这里仅列出生成lib文件，用到的vs中powetshell打开方式：

2.2 结果比较

从结果可以看出，与Python代码相比，FFTW的输出未进行shift，而且仅输出部分有用信息。

3.使用c++实现之2

第2节中的实现使用eigen MatrixXcd 来接收fftw_complex*类型。

本节实现使用double类型，通过(double(*)[2])数组指针来接收fftw_complex*类型。

#include <iostream>
using namespace std;#include "fftw3.h"
#include<array>
int main()
{double vecIn[9] = { 0, 2, 4, 1,6, 1, 3, 2,5};double* vecOut;vecOut = vecIn + 9*sizeof(double);fftw_plan P = fftw_plan_dft_r2c_2d(3, 3, vecIn, (double(*)[2])vecOut, FFTW_ESTIMATE);fftw_execute(P);cout << endl;for (int i = 0; i < 12; i++)cout << *(vecOut + i) << endl;return 0;
}

结果如下：

24
0
-6
4.44089e-16
-3
1.73205
-7.5
4.33013
-3
-1.73205
4.5
0.866025

可以看出，实部和虚部分别存放。

参考文献

[1] FFTW库官网
[2] Windows下FFTW_2.1.5的编译及使用
[3] FFTW 官方文档