數字圖像處理作業_二維離散快速傅里葉變換_2dfft

笨比來學fft了

先不加推理得給出一維得離散傅里葉變換公式

dft代碼實現

def dft(pic_line):#一維離散傅里葉變換暴力n方
    pic_line=np.squeeze(pic_line)
    n=pic_line.size
    w=np.array([[np.exp(-1j*2*np.pi*i*k/n) for k in range(0,n) ] for i in range(0,n)]).T#系數矩陣，生成之后需要轉置
    return pic_line.dot(w)

直接這樣處理一行像素復雜度為n^2.所以需要快速傅里葉變換

旋轉因子具有的性質

推論可以理解為的復數在復平面轉過半個周期

 

base=2 FFT算法

由上面算法我們可以知道一個N點dft可以分治為兩個N/2的DFT

過程是如下圖的蝶形運算式

代碼實現如下

def Butterfly_operation(seq1,seq2):#蝶形運算 seq1為偶序列，seq2為奇序列
    n=seq1.size
    N=n*2
    seq1=np.squeeze(seq1)#降維
    seq2=np.squeeze(seq2)
    res=np.zeros(2*n,dtype=seq1.dtype)
    for i in range(0,n):
        res[i]+=seq1[i]+np.exp(-1j*2*np.pi*i/N)*seq2[i]
        res[n+i]=seq1[i]-np.exp(-1j*2*np.pi*i/N)*seq2[i]
    return  res

之后我們就可以通過遞歸分治，或者二進制倒序數實現重排，來實現一維的fft ,時間復雜度為o(nlogn)

從一維得到二維其實很簡單，就是先對輸入的圖像（m×n）的每一行先做一維fft，把結果存到一個（m×n）的矩陣A中，再對矩陣A的每一列做fft。這時得到的結果就是二維圖像（m×n）的fft。

實現圖像空間域變到頻率域通過二維fft算法大媽如下

 

import numpy as np
import math
import matplotlib.pyplot as plt
from numpy.fft import *
from PIL import Image as IMG
import cv2 as cv
def dft(pic_line):#一維離散傅里葉變換暴力n方
    pic_line=np.squeeze(pic_line)
    n=pic_line.size
    w=np.array([[np.exp(-1j*2*np.pi*i*k/n) for k in range(0,n) ] for i in range(0,n)]).T#系數矩陣，生成之后需要轉置
    return pic_line.dot(w)

def Butterfly_operation(seq1,seq2):#蝶形運算 seq1為偶序列，seq2為奇序列
    n=seq1.size
    N=n*2
    seq1=np.squeeze(seq1)#降維
    seq2=np.squeeze(seq2)
    res=np.zeros(2*n,dtype=seq1.dtype)
    for i in range(0,n):
        res[i]+=seq1[i]+np.exp(-1j*2*np.pi*i/N)*seq2[i]
        res[n+i]=seq1[i]-np.exp(-1j*2*np.pi*i/N)*seq2[i]
    return  res

def fft(img):#遞歸分治實現fft,只對圖像的每一行做fft
    n=img.shape[1]
    if(n%2!=0):
        return ValueError("圖片大小不是2的次冪")#因為fft分治時每次都分為等大的奇數列和偶數列
    if(n<=8):#n足夠小，直接n方對每行求dft
        return np.array([dft(img[i, :]) for i in range(img.shape[0])])
    res_odd=fft(img[:,1::2])
    res_even=fft(img[:,0::2])
    return  np.array([Butterfly_operation(res_even[i:i+1,:],res_odd[i:i+1,:]) for i in range(0,img.shape[0])])

def TwoD_fft(img):
    return fft(fft(img).T).T
def FFT_SHIFT(img):
    M, N = img.shape
    M = int(M / 2)
    N = int(N / 2)
    return np.vstack((np.hstack((img[M:, N:], img[M:, :N])), np.hstack((img[:M, N:], img[:M, :N]))))

if __name__ == '__main__':
    img = cv.imread(r"C:\Users\USER\Pictures\Saved Pictures\test.jpg", cv.IMREAD_COLOR)
    gray = cv.cvtColor(img, cv.COLOR_BGR2GRAY)
    gray=gray[0:1024,0:1024]
    fig, ax = plt.subplots(1, 2)
    my_fft = abs(FFT_SHIFT(TwoD_fft(gray)))
    target = abs(fftshift(fft2(gray)))
    plt.subplot(2, 2, 1)
    plt.imshow(gray)
    plt.subplot(2, 2, 2)
    plt.title('FFT2D')
    plt.imshow(np.log(1 + my_fft))
    plt.subplot(2,2,3)
    plt.title('numpy.fft2')
    plt.imshow(np.log(1+target))
    plt.show()
    plt.title('Lenna')

 

運行結果

參考文獻：

數字信號處理--FFT與蝶形算法--學習筆記

十分簡明易懂的FFT（快速傅里葉變換）