离散傅里叶变换手工实现

离散傅里叶变换手工实现

import numpy as np
import matplotlib.pyplot as plt
import cv2
f = [[4,6],[2,9]]
F = np.complex64(f)
def dft(img,N):W = np.zeros((N,N),dtype=np.complex64)for x in range(N):for y in range(N):W[x][y] = ((np.cos(2*np.pi/N))-(np.sin(2*np.pi/N)*1j))**(x*y) #计算excp(-2πjxy/N) = (cos(2π/N)-jsin(2π/N))^(xy)return W@img@WF = dft(F,2)
i = F.imag
r = F.real
print(F)
F = np.abs(F) #求模即频谱
plt.subplot(121)
plt.title("spatial")
plt.imshow(f,"gray")
plt.subplot(122)
plt.title("frequency")
plt.imshow(F,cmap="gray")
plt.show()

[[21.+0.0000000e+00j -9.-1.8369703e-15j][-1.-1.3471115e-15j  5.+1.2246468e-15j]]

与opencv中的dft对比

在上一文中使用dft后并未进行取模操作，所以最终结果并不是频域，仅仅是取了结果的实部，现已纠正

f = [[4,6],[2,9]]
f = np.float64(f)
F =  cv2.dft(f,flags=cv2.DFT_COMPLEX_OUTPUT)
fr = np.zeros((2,2),np.complex64)
for i in range(2):for k in range(2):for l in range(2):fr[i][k] = F[i][k][0] + F[i][k][1]*1j
fr = np.abs(fr)
plt.subplot(121)
plt.title("spatial")
plt.imshow(f,"gray")
plt.subplot(122)
plt.title("frequency")
plt.imshow(fr,cmap="gray")
plt.show()