IDL实现傅里叶(傅立叶)变换

关于傅里叶变换

这个网上的描述实在太多,原理上就是把空域的图片转换成频域。可以看这张维基百科的示意图

实现公式

IDL实现

;Author:尹全超
;2014-01-17

PRO DFTSELF

;读取图像
file = DIALOG_PICKFILE(/READ, FILTER = '*.bmp')
binary_img = READ_BMP(file)
img01 = image(binary_img,LAYOUT=[2,2,1],title='Original')

;傅里叶变换
N_col=N_elements(binary_img[,0])
M_row=N_elements(binary_img[0,
])
M=M_row
N=N_col
n_dft=complexarr(N,M)
n_dft_tp=complexarr(N,M)
dft_img=complexarr(N,M)
;列(一维DFT)
FOR v=0,N-1,1 DO BEGIN
FOR x=0,M-1,1 DO BEGIN
fxv=complex(0.0d,0.0d)
FOR y=0,N-1,1 DO BEGIN
fxv=fxv+binary_img[y,x]complex(cos(2!pivy/N),-1sin(2!pivy/N))
ENDFOR
n_dft_tp[v,x]=fxv
ENDFOR
ENDFOR
;行(一维DFT)
FOR u=0,M-1,1 DO BEGIN
FOR y=0,N-1,1 DO BEGIN
fxv=complex(0.0d,0.0d)
FOR x=0,M-1,1 DO BEGIN
fxv=fxv+n_dft_tp[y,x]complex(cos(2!piux/M),-1sin(2!piux/M))
ENDFOR
n_dft[y,u]=fxv
ENDFOR
ENDFOR
n_dft=n_dft/(M*N)

powerSpectrum = ABS(n_dft)^2
scaledPowerSpect = ALOG10(powerSpectrum)
img02 = IMAGE(scaledPowerSpect,/CURRENT,LAYOUT=[2,2,2],title='MyFFT')

;自带FFT
ffTransform = FFT(binary_img)
powerSpectrum = ABS(ffTransform)^2
scaledPowerSpect = ALOG10(powerSpectrum)
img03 = IMAGE(scaledPowerSpect,/CURRENT,LAYOUT=[2,2,3],title='FFT')

;反傅立叶变换
;列(一维DFT)
FOR v=0,N-1,1 DO BEGIN
FOR x=0,M-1,1 DO BEGIN
fxv=complex(0.0d,0.0d)
FOR y=0,N-1,1 DO BEGIN
fxv=fxv+n_dft[y,x]complex(cos(2!pivy/N),sin(2!pivy/N))
ENDFOR
n_dft_tp[v,x]=fxv
ENDFOR
ENDFOR
;行(一维DFT)
FOR u=0,M-1,1 DO BEGIN
FOR y=0,N-1,1 DO BEGIN
fxv=complex(0.0d,0.0d)
FOR x=0,M-1,1 DO BEGIN
fxv=fxv+n_dft_tp[y,x]
complex(cos(2!piux/M),sin(2!piux/M))
ENDFOR
dft_img[y,u]=fxv
ENDFOR
ENDFOR
;dft_img=dft_img/(M*N)

img04 = IMAGE(dft_img,/CURRENT,LAYOUT=[2,2,4],title='DFT to IMAGE')

END

这个代码实现了DFT的变换方法,将图像做傅立叶变换,并和IDL自带FFT函数做了一个对比,然后再做了一个反变换。

结果如图

缺陷

目前反傅立叶变换的时候,会产生色差,具体的bug没有找到。如果你找到了相关的方法,欢迎留言,或者给本站投稿。

IDL实现马赫带效应 IDL实现巴特沃斯(butterworth)高通滤波

作者:,GIS爱好者。
分享本文,请您带上本文链接
分享到:

已有 2 条评论

  1. Cook
    1#
    Cook  · 2015-03-24 15:20

发表评论