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*!pi*v*y/N),-1*sin(2*!pi*v*y/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*!pi*u*x/M),-1*sin(2*!pi*u*x/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*!pi*v*y/N),sin(2*!pi*v*y/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*!pi*u*x/M),sin(2*!pi*u*x/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没有找到。如果你找到了相关的方法，欢迎留言，或者给本站投稿。