卷積和快速傅立葉變換（FFT）的實現
阿新 • • 發佈：2019-01-01
卷積運算
卷積可以說是影象處理中最基本的操作。線性濾波通過不同的卷積核，可以產生很多不同的效果。假如有一個要處理的二維影象，通過二維的濾波矩陣（卷積核），對於影象的每一個畫素點，計算它的領域畫素和濾波器矩陣的對應元素的乘積，然後累加，作為該畫素位置的值。關於影象卷積和濾波的一些知識點可以參考這篇部落格。
下面是通過python模擬實現的影象卷積操作，模擬了sobel運算元，prewitt運算元和拉普拉斯運算元。python的np包中有convolve函式可以直接呼叫，scipy中也有scipy.signal.convolve函式可以直接呼叫。
import matplotlib.pyplot as 快速傅立葉變換（FFT）
通過上面的方法實現卷積操作之後，發現可以通過快速傅立葉變換提升卷積運算的效率。python提供了很多標準工具和封裝來計算它。NumPy 和 SciPy 都有經過充分測試的封裝好的FFT庫，分別位於子模組 numpy.fft 和 scipy.fftpack 。有關FFT演算法的原理和推導可以參見參考連結的部落格。
離散傅立葉變換
xn 到 Xk 的轉化就是空域到頻域的轉換，轉化為點值表示法。
def DFT_slow(x):
# Compute the discrete Fourier Transform of the 1D array x
x = np.asarray(x, dtype=float) # 轉化為ndarray
N = x.shape[0] # 維度
n = np.arange(N) # 0~N組成一個一維向量
k = n.reshape((N, 1)) # 轉換為一個N維向量
M = np.exp(-2j * np.pi * k * n / N) # 離散傅立葉公式 -2j複數表示
return np.dot(M, x)
快速傅立葉變換
離散傅立葉變換具有對稱性，利用這種對稱性，可以推出遞迴的快速傅立葉演算法。
def FFT(x):
# A recursive implementation of the 1D Cooley-Tukey FFT
x = np.asarray(x, dtype=float) # 淺拷貝
N = x.shape[0]
if N % 2 > 0:
raise ValueError("size of x must be a power of 2")
elif N <= 32:
return DFT_slow(x)
else:
X_even = FFT(x[::2]) # 從0開始，2為間隔
X_odd = FFT(x[1::2]) # 從1開始，2為間隔
factor = np.exp(-2j * np.pi * np.arange(N) / N)
'''
使用/會出現下面的錯誤，改為// 向下取整
TypeError: slice indices must be integers or None or have an __index__ method
'''
return np.concatenate([X_even + factor[:N // 2] * X_odd,
X_even + factor[N // 2:] * X_odd])向量化的FFT
使用numpy進行矩陣向量化。
def FFT_vectorized(x):
# A vectorized, non-recurisive version of the Cooley_Tukey FFT
x = np.asarray(x, dtype=float)
N = x.shape[0]
if np.log2(N) % 1 > 0:
raise ValueError("size of x must be a power of 2")
# N_min here is equivalent to the stopping condition above,
# and should be a power of 2
N_min = min(N, 32)
# Perform an O[N^2] DFT on all length-N_min sub-problems at once
n = np.arange(N_min)
k = n[:, None]
M = np.exp(-2j * np.pi * n * k / N_min)
X = np.dot(M, x.reshape((N_min, -1)))
# build-up each level of the recursive calculation all at once
while X.shape[0] < N:
X_even = X[:, :X.shape[1] // 2]
X_odd = X[:, X.shape[1] // 2:]
factor = np.exp(-1j * np.pi * np.arange(X.shape[0]) / X.shape[0])[:, None]
X = np.vstack([X_even + factor * X_odd,
X_even - factor * X_odd])
return X.ravel() # 降維用快速傅立葉變換實現卷積
根據卷積定理，時域上卷積運算對應於頻域上的傅立葉變換的乘積。
def fft_convolve(a, b):
n = len(a) + len(b) -1
N = 2 ** (int(np.log2(n))+1)
A = np.fft.fft(a, N)
B = np.fft.fft(b, N)
return np.fft.ifft(A*B)[:n]c++實現的遞迴版
void FFT(Complex* a,int len){
if(len==1) return;
Complex* a0=new Complex[len/2];
Complex* a1=new Complex[len/2];
for(int i=0;i<len;i+=2){
a0[i/2]=a[i];
a1[i/2]=a[i+1];
}
FFT(a0,len/2);FFT(a1,len/2);
Complex wn(cos(2*Pi/len),sin(2*Pi/len));
Complex w(1,0);
for(int i=0;i<(len/2);i++){
a[i]=a0[i]+w*a1[i];
a[i+len/2]=a0[i]-w*a1[i];
w=w*wn;
}
return;
}c++實現的非遞迴版
const double PI = acos(-1.0);
// 複數
struct complex
{
double r,i;
complex(double _r = 0.0, double _i = 0.0)
{
r = _r;
i = _i;
}
complex operator +(const complex &b)
{
return complex(r + b.r, i + b.i);
}
complex operator -(const complex &b)
{
return complex(r - b.r, i - b.i);
}
complex operator *(const complex &b)
{
return complex(r*b.r - i*b.i, r*b.i + i*b.r);
}
};
// 雷德演算法 -- 倒位序
void Rader(complex F[], int len)
{
int j = len >> 1;
for(int i=1; i<len-1; i++)
{
if(i < j)
swap(F[i], F[j]);
int k = len >> 1;
while(j >= k)
{
j -= k;
k >>= 1;
}
if(j < k)
j += k;
}
}
void FFT(complex F[], int len, int on)
{
Rader(F, len);
for(int h=2; h<=len; h<<=1) //計算長度為h的DFT
{
complex wn(cos(-on*2*PI/h), sin(-on*2*PI/h)); //單位復根e^(2*PI/m),用尤拉公式展開
for(int j=0; j<len; j+=h)
{
complex w(1,0); //旋轉因子
for(int k=j; k<j+h/2; k++)
{
complex u = F[k];
complex t = w * F[k + h/2];
F[k] = u + t; //蝴蝶合併操作
F[k + h/2] = u - t;
w = w * wn; //更新旋轉因子
}
}
}
//求逆傅立葉變換
if(on == -1)
{
for(int i=0; i<len; i++)
{
F[i].r /= len;
}
}
}
//求卷積
void Conv(complex a[], complex b[], int len)
{
FFT(a, len, 1);
FFT(b, len, 1);
for(int i=0; i<len; i++)
{
a[i] = a[i] * b[i]; //卷積定理
}
FFT(a, len, -1);
}