离散傅里叶变换手工实现
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@W
F = 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()