做网站1g1核够吗,同泰公司网站公司查询,做网站代下,网站开发合同怎么写二维离散傅里叶变换的实现 1.使用Python包实现1.1 fftshift在numpy中的实现1.2 平移后的幅度谱 2.使用c实现之12.1 FFTW库安装2.2 结果比较 3.使用c实现之2参考文献 1.使用Python包实现
import numpy as np
import matplotlib.pyplot as plt
anp.array([0, 2, 4, 1,6, 1, 3, … 二维离散傅里叶变换的实现 1.使用Python包实现1.1 fftshift在numpy中的实现1.2 平移后的幅度谱 2.使用c实现之12.1 FFTW库安装2.2 结果比较 3.使用c实现之2参考文献 1.使用Python包实现
import numpy as np
import matplotlib.pyplot as plt
anp.array([0, 2, 4, 1,6, 1, 3, 2,5]).reshape(3,3)
fnp.fft.fft2(a)
fshiftnp.fft.fftshift(f)
mag20*np.log(np.abs(fshift))
plt.axis(off)
plt.imshow(mag)
plt.show()其中
f:
array([[24. 0.00000000e00j, -6. 4.44089210e-16j, -6. -4.44089210e-16j],[-3. 1.73205081e00j, -7.54.33012702e00j, 4.5-8.66025404e-01j],[-3. -1.73205081e00j, 4.58.66025404e-01j, -7.5-4.33012702e00j]])fshift:
array([[-7.5-4.33012702e00j, -3. -1.73205081e00j, 4.58.66025404e-01j],[-6. -4.44089210e-16j, 24. 0.00000000e00j, -6. 4.44089210e-16j],[ 4.5-8.66025404e-01j, -3. 1.73205081e00j, -7.54.33012702e00j]])1.1 fftshift在numpy中的实现
通过官方代码可以看出fftshift是通过np.roll实现的。
def fftshift(x, axesNone):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.73205081e00j, 4.58.66025404e-01j, -7.5-4.33012702e00j],[24. 0.00000000e00j, -6. 4.44089210e-16j, -6. -4.44089210e-16j],[-3. 1.73205081e00j, -7.54.33012702e00j, 4.5-8.66025404e-01j]])np.roll(np.roll(f, 1,0),1,1)
array([[-7.5-4.33012702e00j, -3. -1.73205081e00j, 4.58.66025404e-01j],[-6. -4.44089210e-16j, 24. 0.00000000e00j, -6. 4.44089210e-16j],[ 4.5-8.66025404e-01j, -3. 1.73205081e00j, -7.54.33012702e00j]])与np.fft.fftshift(f)结果一致。
1.2 平移后的幅度谱 2.使用c实现之1
#include iostream
using namespace std;
#includeEigen/Dense
using namespace Eigen;
#include fftw3.hint 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 52.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
#includearray
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 官方文档
