图像处理中, 为了方便处理,便于抽取特征,数据压缩等目的,常常要将图像进行变换。
一般有如下变换方法
- 傅立叶变换Fourier Transform
- 离散余弦变换Discrete Cosine Transform
- 沃尔希-哈德玛变换Walsh-Hadamard Transform
- 斜变换Slant Transform
- 哈尔变换Haar Transform
- 离散K-L变换Discrete Karhunen-Leave Transform
- 奇异值分解SVD变换Singular-Value Decomposition
- 离散小波变换Discrete Wavelet Transform
这篇文章介绍一下傅里叶变换
0.1. 定义
0.1.1. 连续
积分形式
如果一个函数的绝对值的积分存在,即
∫
−
∞
∞
∣
h
(
t
)
∣
d
t
<
∞
\int_{-\infty} ^\infty |h(t)|dt<\infty
∫−∞∞∣h(t)∣dt<∞
并且函数是连续的或者只有有限个不连续点,则对于 x 的任何值, 函数的傅里叶变换存在
- 一维傅里叶变换
H ( f ) = ∫ − ∞ ∞ h ( t ) e − j 2 π f t d t H(f)=\int_{-\infty} ^\infty h(t)e^{-j2\pi ft}dt H(f)=∫−∞∞h(t)e−j2πftdt - 一维傅里叶逆变换
H ( f ) = ∫ − ∞ ∞ h ( t ) e j 2 π f t d t H(f)=\int_{-\infty} ^\infty h(t)e^{j2\pi ft}dt H(f)=∫−∞∞h(t)ej2πftdt
同理多重积分
0.1.2. 离散
实际应用中,多用离散傅里叶变换 DFT.
- 一维傅里叶变换
F ( u ) = ∑ x = 0 N − 1 f ( x ) e − 2 π j N u x F(u)=\sum_{x=0} ^{N-1} f(x)e^{\frac{-2\pi j}{N} ux} F(u)=x=0∑N−1f(x)eN−2πjux - 一维傅里叶逆变换
f ( x ) = 1 N ∑ u = 0 N − 1 F ( u ) e 2 π j N u x f(x)=\frac{1}{N}\sum_{u=0} ^{N-1} F(u)e^{\frac{2\pi j}{N} ux} f(x)=N1u=0∑N−1F(u)eN2πjux
需要注意的是, 逆变换乘以 1 N \frac{1}{N} N1 是为了归一化,这个系数可以随意改变, 即可以正变换乘以 1 N \frac{1}{N} N1, 逆变换就不乘,或者两者都乘以 1 N \frac{1}{\sqrt{N}} N1等系数。 - 二维傅里叶变换
F ( u , v ) = 1 N ∑ x = 0 N − 1 ∑ y = 0 N − 1 f ( x , y ) e − 2 π j N ( u x + v y ) F(u,v)=\frac{1}{N}\sum_{x=0}^{N-1}\sum_{y=0} ^{N-1} f(x,y)e^{\frac{-2\pi j}{N} (ux+vy)} F(u,v)=N1x=0∑N−1y=0∑N−1f(x,y)eN−2πj(ux+vy) - 二维傅里叶逆变换
f ( x , y ) = 1 N ∑ u = 0 N − 1 ∑ v = 0 N − 1 F ( u , v ) e 2 π j N ( u x + v y ) f(x,y)=\frac{1}{N}\sum_{u=0}^{N-1}\sum_{v=0} ^{N-1} F(u,v)e^{\frac{2\pi j}{N} (ux+vy)} f(x,y)=N1u=0∑N−1v=0∑N−1F(u,v)eN2πj(ux+vy)
幅度
∣
F
(
u
,
v
)
∣
=
r
e
a
l
(
F
)
2
+
i
m
a
g
(
F
)
2
|F(u,v)| = \sqrt{real(F)^2+imag(F)^2}
∣F(u,v)∣=real(F)2+imag(F)2
相位
a
r
c
t
a
n
i
m
a
g
(
F
)
r
e
a
l
(
F
)
arctan{\frac{imag(F)}{real(F)}}
arctanreal(F)imag(F)
对于图像的幅度谱显示,由于 |F(u,v)| 变换范围太大,一般显示
D
=
l
o
g
(
∣
F
(
u
,
v
)
+
1
)
D= log(|F(u,v)+1)
D=log(∣F(u,v)+1)
用 <=>
表示傅里叶变换对
f
(
x
)
<
=
>
F
(
u
)
f
(
x
,
y
)
<
=
>
F
(
u
,
v
)
f(x)<=>F(u)\\ f(x,y)<=>F(u,v)
f(x)<=>F(u)f(x,y)<=>F(u,v)
f,g,h 对应的傅里叶变换 F,G,H
F ∗ F^* F∗ 表示 F F F 的共轭
0.2. 性质
0.2.1. 分离性
F
(
x
,
v
)
=
∑
y
=
0
N
−
1
f
(
x
,
y
)
e
−
2
π
j
N
v
y
F
(
u
,
v
)
=
1
N
∑
x
=
0
N
−
1
F
(
x
,
v
)
e
−
2
π
j
N
u
x
\begin{aligned} &F(x,v)=\sum_{y=0} ^{N-1} f(x,y)e^{\frac{-2\pi j}{N} vy}\\ &F(u,v)=\frac{1}{N}\sum_{x=0}^{N-1}F(x,v)e^{\frac{-2\pi j}{N}ux} \end{aligned}
F(x,v)=y=0∑N−1f(x,y)eN−2πjvyF(u,v)=N1x=0∑N−1F(x,v)eN−2πjux
进行多维变换时,可以依次对每一维进行变换。 下面在代码中就是这样实现的。
0.2.2. 位移定理
f
(
x
,
y
)
e
2
π
j
N
(
u
0
x
+
v
0
y
)
<
=
>
F
(
u
−
u
0
,
v
−
v
0
)
f(x,y)e^{\frac{2\pi j}{N}(u_0x+v_0y)} <=>F(u-u_0,v-v_0)
f(x,y)eN2πj(u0x+v0y)<=>F(u−u0,v−v0)
f
(
x
−
x
0
,
y
−
y
0
)
<
=
>
F
(
u
,
v
)
e
−
2
π
j
N
(
u
x
0
+
v
y
0
)
f(x-x_0,y-y_0)<=>F(u,v)e^{\frac{-2\pi j}{N}(ux_0+vy_0)}
f(x−x0,y−y0)<=>F(u,v)eN−2πj(ux0+vy0)
0.2.3. 周期性
F ( u , v ) = F ( u + N , v + N ) F(u,v) = F(u+N,v+N) F(u,v)=F(u+N,v+N)
0.2.4. 共轭对称性
F
(
u
,
v
)
=
F
∗
(
−
u
,
−
v
)
F(u,v) = F^*(-u,-v)
F(u,v)=F∗(−u,−v)
a)偶分量函数在变换中产生偶分量函数;
b)奇分量函数在变换中产生奇分量函数;
c)奇分量函数在变换中引入系数-j;
d)偶分量函数在变换中不引入系数.
0.2.5. 旋转性
if
f
(
r
,
θ
)
<
=
>
F
(
ω
,
ϕ
)
f(r,\theta)<=>F(\omega,\phi)
f(r,θ)<=>F(ω,ϕ)
then
f
(
r
,
θ
+
t
)
<
=
>
F
(
ω
,
ϕ
+
t
)
f(r,\theta+t)<=>F(\omega,\phi+t)
f(r,θ+t)<=>F(ω,ϕ+t)
0.2.6. 加法定理
F
o
u
r
i
e
r
[
f
+
g
]
=
F
o
u
r
i
e
r
[
f
]
+
F
o
u
r
i
e
r
[
g
]
Fourier[f+g]=Fourier[f]+Fourier[g]
Fourier[f+g]=Fourier[f]+Fourier[g]
2.
a
f
(
x
,
y
)
<
=
>
a
F
[
u
,
v
]
af(x,y)<=>aF[u,v]
af(x,y)<=>aF[u,v]
0.2.7. 平均值
1 N 2 ∑ x = 0 N − 1 ∑ y = 0 N − 1 f ( x , y ) = 1 N F ( 0 , 0 ) \frac{1}{N^2}\sum_{x=0}^{N-1}\sum_{y=0} ^{N-1} f(x,y) = \frac{1}{N}F(0,0) N21x=0∑N−1y=0∑N−1f(x,y)=N1F(0,0)
0.2.8. 相似性定理
尺度变换
f
(
a
x
,
b
y
)
<
=
>
F
(
u
a
,
v
b
)
a
b
f(ax,by)<=>\frac{F(\frac{u}{a},\frac{v}{b})}{ab}
f(ax,by)<=>abF(au,bv)
0.2.9. 卷积定理
卷积定义
1d
f
∗
g
=
1
M
∑
m
=
0
M
−
1
f
(
m
)
g
(
x
−
m
)
f*g = \frac{1}{M}\sum_{m=0}^{M-1}f(m)g(x-m)
f∗g=M1m=0∑M−1f(m)g(x−m)
2d
f
(
x
,
y
)
∗
g
(
x
,
y
)
=
1
M
N
∑
m
=
0
M
−
1
∑
n
=
0
N
−
1
f
(
m
,
n
)
g
(
x
−
m
,
y
−
n
)
f(x,y)*g(x,y) = \frac{1}{MN}\sum_{m=0}^{M-1}\sum_{n=0}^{N-1}f(m,n)g(x-m,y-n)
f(x,y)∗g(x,y)=MN1m=0∑M−1n=0∑N−1f(m,n)g(x−m,y−n)
卷积定理
f
(
x
,
y
)
∗
g
(
x
,
y
)
<
=
>
F
(
u
,
v
)
G
(
u
,
v
)
f(x,y)*g(x,y) <=> F(u,v)G(u,v)
f(x,y)∗g(x,y)<=>F(u,v)G(u,v)
f
(
x
,
y
)
g
(
x
,
y
)
<
=
>
F
(
u
,
v
)
∗
G
(
u
,
v
)
f(x,y)g(x,y)<=>F(u,v)*G(u,v)
f(x,y)g(x,y)<=>F(u,v)∗G(u,v)
离散卷积
用
∑
i
=
0
N
−
1
x
(
i
T
)
h
[
(
k
−
i
)
T
]
<
=
>
X
(
n
N
T
)
H
(
n
N
T
)
\sum_{i=0}^{N-1}x(iT)h[(k-i)T] <=> X(\frac{n}{NT})H(\frac{n}{NT})
i=0∑N−1x(iT)h[(k−i)T]<=>X(NTn)H(NTn)
即两个周期为 N 的抽样函数, 他们的卷积的离散傅里叶变换等于他们的离散傅里叶变换的卷积
卷积的应用:
去除噪声, 特征增强
两个不同周期的信号卷积需要周期扩展的原因:如果直接进行傅里叶变换和乘积,会产生折叠误差(卷绕)。
0.2.10. 相关定理
下面用$ \infty$ 表示相关。
相关函数描述了两个信号之间的相似性,其相关性大小有相关系数衡量
- 相关函数的定义
离散
f ( x , y ) ∞ g ( x , y ) = 1 M N ∑ m = 0 M − 1 ∑ n = 0 N − 1 f ∗ ( m , n ) g ( x + m , y + n ) f(x,y)\quad \infty \quad g(x,y) = \frac{1}{MN}\sum_{m=0}^{M-1}\sum_{n=0}^{N-1}f^*(m,n)g(x+m,y+n) f(x,y)∞g(x,y)=MN1m=0∑M−1n=0∑N−1f∗(m,n)g(x+m,y+n)
连续
z ( t ) = ∫ − ∞ ∞ x ∗ ( τ ) h ( t + τ ) d τ z(t) = \int_{-\infty}^{\infty}x^*(\tau) h(t+\tau)d\tau z(t)=∫−∞∞x∗(τ)h(t+τ)dτ - 定理
f ( x , y ) ∞ g ( x , y ) < = > F ∗ ( u , v ) G ( u , v ) f(x,y)\quad \infty \quad g(x,y)<=>F^*(u,v)G(u,v) f(x,y)∞g(x,y)<=>F∗(u,v)G(u,v)
0.2.11. Rayleigh 定理
能量变换
对于有限区间非零函数 f(t), 其能量为
E
=
∫
−
∞
∞
∣
f
(
t
)
∣
2
d
t
E = \int_{-\infty}^{\infty}|f(t)|^2dt
E=∫−∞∞∣f(t)∣2dt
其变换函数与原函数有相同的能量
∫
−
∞
∞
∣
f
(
t
)
∣
2
d
t
=
∫
−
∞
∞
∣
F
(
u
)
∣
2
d
t
\int_{-\infty}^{\infty}|f(t)|^2dt = \int_{-\infty}^{\infty}|F(u)|^2dt
∫−∞∞∣f(t)∣2dt=∫−∞∞∣F(u)∣2dt
0.3. 快速傅里叶变换
由上面离散傅里叶变换的性质易知,直接计算 1维 dft 的时间复杂度维 O ( N 2 ) O(N^2) O(N2)。
利用到单位根的对称性,快速傅里叶变换可以达到 O ( n l o g n ) O(nlogn) O(nlogn)的时间复杂度。
0.3.1. 复数中的单位根
我们知道, 在复平面,复数
c
o
s
θ
+
i
s
i
n
θ
cos\theta +i\ sin\theta
cosθ+i sinθk可以表示成
e
i
θ
e^{i\theta}
eiθ, 可以对应一个向量。
θ
\theta
θ即为幅角。
在单位圆中 ,单位圆被分成
2
π
θ
\frac{2\pi}{\theta}
θ2π 份, 由单位圆的对称性
e
i
θ
=
e
i
(
θ
+
2
π
)
e^{i\theta} = e^{i(\theta+2\pi)}
eiθ=ei(θ+2π)
现在记 $ n =\frac{ 2\pi }{\theta}$ , 即被分成 n 份,幅度角为正且最小的向量称为 n 次单位向量, 记为
ω
n
\omega _n
ωn,
其余的 n-1 个向量分别为
ω
n
2
,
ω
n
3
,
…
,
ω
n
n
\omega_{n}^{2},\omega_{n}^{3},\ldots,\omega_{n}^{n}
ωn2,ωn3,…,ωnn ,它们可以由复数之间的乘法得来 $w_{n}{k}=w_{n}{k-1}\cdot w_{n}^{1}\ (2 \leq k \leq n) $。
单位根的性质
- 这个可以用 e 表示出来证明
ω 2 n 2 k = ω n k \omega_{2n}^{2k}=\omega_{n}^{k} ω2n2k=ωnk - 可以写成三角函数证明
ω n k + n 2 = − ω n k \omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^{k} ωnk+2n=−ωnk
容易看出 $w_{n}{n}=w_{n}{0}=1 $。
对于$ w_{n}^{k}$ , 它事实上就是 e 2 π i n k e^{\frac{2\pi i}{n}k} en2πik 。
0.3.2. 快速傅里叶变换的计算
下面的推导假设 n = 2 k n=2^k n=2k,以及代码实现 FFT 部分也是 如此。
利用上面的对称性,
将傅里叶计算进行奇偶分组
F
(
u
)
=
∑
i
=
0
n
−
1
ω
n
i
u
a
i
=
∑
i
=
0
n
2
−
1
ω
n
2
i
u
a
2
i
+
∑
i
=
0
n
2
−
1
ω
n
(
2
i
+
1
)
u
a
2
i
+
1
=
∑
i
=
0
n
2
−
1
ω
n
2
i
u
a
2
i
+
ω
n
u
∑
i
=
0
n
2
−
1
ω
n
2
i
u
a
2
i
+
1
=
F
e
v
e
n
(
u
)
+
ω
n
u
F
o
d
d
(
u
)
\begin{aligned} F(u)&=\sum_{i=0}^{n-1}\omega_n ^{iu} a^i\\ &= \sum_{i=0}^{\frac{n}{2}-1}\omega_n ^{2iu} a^{2i}+\sum_{i=0}^{\frac{n}{2}-1}\omega_n ^{(2i+1)u} a^{2i+1}\\ &=\sum_{i=0}^{\frac{n}{2}-1}\omega_{\frac{n}{2}} ^{iu} a^{2i}+\omega_n^u\sum_{i=0}^{\frac{n}{2}-1}\omega_{\frac{n}{2}} ^{iu} a^{2i+1}\\ & = F_{even}(u)+\omega_n^u F_{odd}(u) \end{aligned}
F(u)=i=0∑n−1ωniuai=i=0∑2n−1ωn2iua2i+i=0∑2n−1ωn(2i+1)ua2i+1=i=0∑2n−1ω2niua2i+ωnui=0∑2n−1ω2niua2i+1=Feven(u)+ωnuFodd(u)
F
e
v
e
n
F_{even}
Feven表示将 输入的次序中偶数点进行 Fourier 变换,
F
o
d
d
F_{odd}
Fodd 同理,这样就形成递推公式。
现在还没有减少计算量,下面通过将分别计算的 奇项,偶项利用起来,只计算 前
n
2
−
1
\frac{n}{2}-1
2n−1项,后面的一半可以利用此结果马上算出来。每一次可以减少一半的计算量。
对于
n
2
≤
i
+
n
2
≤
n
−
1
\frac{n}{2}\leq i+\frac{n}{2}\leq n-1
2n≤i+2n≤n−1
F
(
ω
n
i
+
n
2
)
=
F
e
v
e
n
(
ω
n
2
i
+
n
)
+
ω
n
i
+
n
2
⋅
F
o
d
d
(
ω
n
2
i
+
n
)
=
F
e
v
e
n
(
ω
n
2
i
+
n
2
)
+
ω
n
2
i
+
n
2
⋅
F
o
d
d
(
ω
n
2
i
+
n
2
)
=
F
e
v
e
n
(
ω
n
2
i
)
−
ω
n
2
i
⋅
F
o
d
d
(
ω
n
2
i
)
\begin{aligned} F(\omega_{n}^{i+\frac{n}{2}})&=F_{even}(\omega_{n}^{2i+n})+\omega_{n}^{i+\frac{n}{2}}\cdot F_{odd}(\omega_{n}^{2i+n})\\ &=F_{even}(\omega_{\frac{n}{2}}^{i+\frac{n}{2}})+\omega_{\frac{n}{2}}^{i+\frac{n}{2}}\cdot F_{odd}(\omega_{\frac{n}{2}}^{i+\frac{n}{2}})\\ & =F_{even}(\omega_{\frac{n}{2}}^{i})-\omega_{\frac{n}{2}}^{i}\cdot F_{odd}(\omega_{\frac{n}{2}}^{i}) \end{aligned}
F(ωni+2n)=Feven(ωn2i+n)+ωni+2n⋅Fodd(ωn2i+n)=Feven(ω2ni+2n)+ω2ni+2n⋅Fodd(ω2ni+2n)=Feven(ω2ni)−ω2ni⋅Fodd(ω2ni)
现在很清楚了,在每次计算 a[0…n-1] 的傅里叶变换F[0…n-1],分别计算出奇 odd[0…n/2-1],偶even[0…n/2-1](可以递归地进行),
那么傅里叶变换为:
F
[
i
]
=
{
e
v
e
n
[
i
]
+
ω
i
⋅
o
d
d
[
i
]
,
i
<
n
2
e
v
e
n
[
i
]
−
ω
i
⋅
o
d
d
[
i
]
,
e
l
s
e
F[i] = \begin{cases} even[i]+ \omega^i \cdot odd[i], \quad i<\frac{n}{2}\\ even[i]- \omega^i \cdot odd[i], \quad else \end{cases}
F[i]={even[i]+ωi⋅odd[i],i<2neven[i]−ωi⋅odd[i],else
0.4. 代码
下面是 python 实现
一维用 FFT 实现, 不过 只实现了 2 的幂。/ 对于非 2 的幂,用 FFT 实现有点困难,还需要插值,所以我 用
O
(
n
2
)
O(n^2)
O(n2) 直接实现。
二维的 DFT利用 分离性,直接调用 一维 FFT。
GitHub
import numpy as np
def _fft(a, invert=False):
N = len(a)
if N == 1:
return [a[0]]
elif N & (N - 1) == 0: # O(nlogn), 2^k
even = _fft(a[::2], invert)
odd = _fft(a[1::2], invert)
i = 2j if invert else -2j
factor = np.exp(i * np.pi * np.arange(N // 2) / N)
prod = factor * odd
return np.concatenate([even + prod, even - prod])
else: # O(n^2)
w = np.arange(N)
i = 2j if invert else -2j
m = w.reshape((N, 1)) * w
W = np.exp(m * i * np.pi / N)
return np.concatenate(np.dot(W, a.reshape(
(N, 1)))) # important, cannot use *
def fft(a):
'''fourier[a]'''
n = len(a)
if n == 0:
raise Exception("[Error]: Invalid length: 0")
return _fft(a)
def ifft(a):
'''invert fourier[a]'''
n = len(a)
if n == 0:
raise Exception("[Error]: Invalid length: 0")
return _fft(a, True) / n
def fft2(arr):
return np.apply_along_axis(fft, 0,
np.apply_along_axis(fft, 1, np.asarray(arr)))
def ifft2(arr):
return np.apply_along_axis(ifft, 0,
np.apply_along_axis(ifft, 1, np.asarray(arr)))
def test(n=128):
print('\nsequence length:', n)
print('fft')
li = np.random.random(n)
print(np.allclose(fft(li), np.fft.fft(li)))
print('ifft')
li = np.random.random(n)
print(np.allclose(ifft(li), np.fft.ifft(li)))
print('fft2')
li = np.random.random(n * n).reshape((n, n))
print(np.allclose(fft2(li), np.fft.fft2(li)))
print('ifft2')
li = np.random.random(n * n).reshape((n, n))
print(np.allclose(ifft2(li), np.fft.ifft2(li)))
if __name__ == '__main__':
for i in range(1, 3):
test(i * 16)
0.5. 参考
图像处理中, 为了方便处理,便于抽取特征,数据压缩等目的,常常要将图像进行变换。
一般有如下变换方法
- 傅立叶变换Fourier Transform
- 离散余弦变换Discrete Cosine Transform
- 沃尔希-哈德玛变换Walsh-Hadamard Transform
- 斜变换Slant Transform
- 哈尔变换Haar Transform
- 离散K-L变换Discrete Karhunen-Leave Transform
- 奇异值分解SVD变换Singular-Value Decomposition
- 离散小波变换Discrete Wavelet Transform
这篇文章介绍一下傅里叶变换
0.1. 定义
0.1.1. 连续
积分形式
如果一个函数的绝对值的积分存在,即
∫
−
∞
∞
∣
h
(
t
)
∣
d
t
<
∞
\int_{-\infty} ^\infty |h(t)|dt<\infty
∫−∞∞∣h(t)∣dt<∞
并且函数是连续的或者只有有限个不连续点,则对于 x 的任何值, 函数的傅里叶变换存在
- 一维傅里叶变换
H ( f ) = ∫ − ∞ ∞ h ( t ) e − j 2 π f t d t H(f)=\int_{-\infty} ^\infty h(t)e^{-j2\pi ft}dt H(f)=∫−∞∞h(t)e−j2πftdt - 一维傅里叶逆变换
H ( f ) = ∫ − ∞ ∞ h ( t ) e j 2 π f t d t H(f)=\int_{-\infty} ^\infty h(t)e^{j2\pi ft}dt H(f)=∫−∞∞h(t)ej2πftdt
同理多重积分
0.1.2. 离散
实际应用中,多用离散傅里叶变换 DFT.
- 一维傅里叶变换
F ( u ) = ∑ x = 0 N − 1 f ( x ) e − 2 π j N u x F(u)=\sum_{x=0} ^{N-1} f(x)e^{\frac{-2\pi j}{N} ux} F(u)=x=0∑N−1f(x)eN−2πjux - 一维傅里叶逆变换
f ( x ) = 1 N ∑ u = 0 N − 1 F ( u ) e 2 π j N u x f(x)=\frac{1}{N}\sum_{u=0} ^{N-1} F(u)e^{\frac{2\pi j}{N} ux} f(x)=N1u=0∑N−1F(u)eN2πjux
需要注意的是, 逆变换乘以 1 N \frac{1}{N} N1 是为了归一化,这个系数可以随意改变, 即可以正变换乘以 1 N \frac{1}{N} N1, 逆变换就不乘,或者两者都乘以 1 N \frac{1}{\sqrt{N}} N1等系数。 - 二维傅里叶变换
F ( u , v ) = 1 N ∑ x = 0 N − 1 ∑ y = 0 N − 1 f ( x , y ) e − 2 π j N ( u x + v y ) F(u,v)=\frac{1}{N}\sum_{x=0}^{N-1}\sum_{y=0} ^{N-1} f(x,y)e^{\frac{-2\pi j}{N} (ux+vy)} F(u,v)=N1x=0∑N−1y=0∑N−1f(x,y)eN−2πj(ux+vy) - 二维傅里叶逆变换
f ( x , y ) = 1 N ∑ u = 0 N − 1 ∑ v = 0 N − 1 F ( u , v ) e 2 π j N ( u x + v y ) f(x,y)=\frac{1}{N}\sum_{u=0}^{N-1}\sum_{v=0} ^{N-1} F(u,v)e^{\frac{2\pi j}{N} (ux+vy)} f(x,y)=N1u=0∑N−1v=0∑N−1F(u,v)eN2πj(ux+vy)
幅度
∣
F
(
u
,
v
)
∣
=
r
e
a
l
(
F
)
2
+
i
m
a
g
(
F
)
2
|F(u,v)| = \sqrt{real(F)^2+imag(F)^2}
∣F(u,v)∣=real(F)2+imag(F)2
相位
a
r
c
t
a
n
i
m
a
g
(
F
)
r
e
a
l
(
F
)
arctan{\frac{imag(F)}{real(F)}}
arctanreal(F)imag(F)
对于图像的幅度谱显示,由于 |F(u,v)| 变换范围太大,一般显示
D
=
l
o
g
(
∣
F
(
u
,
v
)
+
1
)
D= log(|F(u,v)+1)
D=log(∣F(u,v)+1)
用 <=>
表示傅里叶变换对
f
(
x
)
<
=
>
F
(
u
)
f
(
x
,
y
)
<
=
>
F
(
u
,
v
)
f(x)<=>F(u)\\ f(x,y)<=>F(u,v)
f(x)<=>F(u)f(x,y)<=>F(u,v)
f,g,h 对应的傅里叶变换 F,G,H
F ∗ F^* F∗ 表示 F F F 的共轭
0.2. 性质
0.2.1. 分离性
F
(
x
,
v
)
=
∑
y
=
0
N
−
1
f
(
x
,
y
)
e
−
2
π
j
N
v
y
F
(
u
,
v
)
=
1
N
∑
x
=
0
N
−
1
F
(
x
,
v
)
e
−
2
π
j
N
u
x
\begin{aligned} &F(x,v)=\sum_{y=0} ^{N-1} f(x,y)e^{\frac{-2\pi j}{N} vy}\\ &F(u,v)=\frac{1}{N}\sum_{x=0}^{N-1}F(x,v)e^{\frac{-2\pi j}{N}ux} \end{aligned}
F(x,v)=y=0∑N−1f(x,y)eN−2πjvyF(u,v)=N1x=0∑N−1F(x,v)eN−2πjux
进行多维变换时,可以依次对每一维进行变换。 下面在代码中就是这样实现的。
0.2.2. 位移定理
f
(
x
,
y
)
e
2
π
j
N
(
u
0
x
+
v
0
y
)
<
=
>
F
(
u
−
u
0
,
v
−
v
0
)
f(x,y)e^{\frac{2\pi j}{N}(u_0x+v_0y)} <=>F(u-u_0,v-v_0)
f(x,y)eN2πj(u0x+v0y)<=>F(u−u0,v−v0)
f
(
x
−
x
0
,
y
−
y
0
)
<
=
>
F
(
u
,
v
)
e
−
2
π
j
N
(
u
x
0
+
v
y
0
)
f(x-x_0,y-y_0)<=>F(u,v)e^{\frac{-2\pi j}{N}(ux_0+vy_0)}
f(x−x0,y−y0)<=>F(u,v)eN−2πj(ux0+vy0)
0.2.3. 周期性
F ( u , v ) = F ( u + N , v + N ) F(u,v) = F(u+N,v+N) F(u,v)=F(u+N,v+N)
0.2.4. 共轭对称性
F
(
u
,
v
)
=
F
∗
(
−
u
,
−
v
)
F(u,v) = F^*(-u,-v)
F(u,v)=F∗(−u,−v)
a)偶分量函数在变换中产生偶分量函数;
b)奇分量函数在变换中产生奇分量函数;
c)奇分量函数在变换中引入系数-j;
d)偶分量函数在变换中不引入系数.
0.2.5. 旋转性
if
f
(
r
,
θ
)
<
=
>
F
(
ω
,
ϕ
)
f(r,\theta)<=>F(\omega,\phi)
f(r,θ)<=>F(ω,ϕ)
then
f
(
r
,
θ
+
t
)
<
=
>
F
(
ω
,
ϕ
+
t
)
f(r,\theta+t)<=>F(\omega,\phi+t)
f(r,θ+t)<=>F(ω,ϕ+t)
0.2.6. 加法定理
F
o
u
r
i
e
r
[
f
+
g
]
=
F
o
u
r
i
e
r
[
f
]
+
F
o
u
r
i
e
r
[
g
]
Fourier[f+g]=Fourier[f]+Fourier[g]
Fourier[f+g]=Fourier[f]+Fourier[g]
2.
a
f
(
x
,
y
)
<
=
>
a
F
[
u
,
v
]
af(x,y)<=>aF[u,v]
af(x,y)<=>aF[u,v]
0.2.7. 平均值
1 N 2 ∑ x = 0 N − 1 ∑ y = 0 N − 1 f ( x , y ) = 1 N F ( 0 , 0 ) \frac{1}{N^2}\sum_{x=0}^{N-1}\sum_{y=0} ^{N-1} f(x,y) = \frac{1}{N}F(0,0) N21x=0∑N−1y=0∑N−1f(x,y)=N1F(0,0)
0.2.8. 相似性定理
尺度变换
f
(
a
x
,
b
y
)
<
=
>
F
(
u
a
,
v
b
)
a
b
f(ax,by)<=>\frac{F(\frac{u}{a},\frac{v}{b})}{ab}
f(ax,by)<=>abF(au,bv)
0.2.9. 卷积定理
卷积定义
1d
f
∗
g
=
1
M
∑
m
=
0
M
−
1
f
(
m
)
g
(
x
−
m
)
f*g = \frac{1}{M}\sum_{m=0}^{M-1}f(m)g(x-m)
f∗g=M1m=0∑M−1f(m)g(x−m)
2d
f
(
x
,
y
)
∗
g
(
x
,
y
)
=
1
M
N
∑
m
=
0
M
−
1
∑
n
=
0
N
−
1
f
(
m
,
n
)
g
(
x
−
m
,
y
−
n
)
f(x,y)*g(x,y) = \frac{1}{MN}\sum_{m=0}^{M-1}\sum_{n=0}^{N-1}f(m,n)g(x-m,y-n)
f(x,y)∗g(x,y)=MN1m=0∑M−1n=0∑N−1f(m,n)g(x−m,y−n)
卷积定理
f
(
x
,
y
)
∗
g
(
x
,
y
)
<
=
>
F
(
u
,
v
)
G
(
u
,
v
)
f(x,y)*g(x,y) <=> F(u,v)G(u,v)
f(x,y)∗g(x,y)<=>F(u,v)G(u,v)
f
(
x
,
y
)
g
(
x
,
y
)
<
=
>
F
(
u
,
v
)
∗
G
(
u
,
v
)
f(x,y)g(x,y)<=>F(u,v)*G(u,v)
f(x,y)g(x,y)<=>F(u,v)∗G(u,v)
离散卷积
用
∑
i
=
0
N
−
1
x
(
i
T
)
h
[
(
k
−
i
)
T
]
<
=
>
X
(
n
N
T
)
H
(
n
N
T
)
\sum_{i=0}^{N-1}x(iT)h[(k-i)T] <=> X(\frac{n}{NT})H(\frac{n}{NT})
i=0∑N−1x(iT)h[(k−i)T]<=>X(NTn)H(NTn)
即两个周期为 N 的抽样函数, 他们的卷积的离散傅里叶变换等于他们的离散傅里叶变换的卷积
卷积的应用:
去除噪声, 特征增强
两个不同周期的信号卷积需要周期扩展的原因:如果直接进行傅里叶变换和乘积,会产生折叠误差(卷绕)。
0.2.10. 相关定理
下面用$ \infty$ 表示相关。
相关函数描述了两个信号之间的相似性,其相关性大小有相关系数衡量
- 相关函数的定义
离散
f ( x , y ) ∞ g ( x , y ) = 1 M N ∑ m = 0 M − 1 ∑ n = 0 N − 1 f ∗ ( m , n ) g ( x + m , y + n ) f(x,y)\quad \infty \quad g(x,y) = \frac{1}{MN}\sum_{m=0}^{M-1}\sum_{n=0}^{N-1}f^*(m,n)g(x+m,y+n) f(x,y)∞g(x,y)=MN1m=0∑M−1n=0∑N−1f∗(m,n)g(x+m,y+n)
连续
z ( t ) = ∫ − ∞ ∞ x ∗ ( τ ) h ( t + τ ) d τ z(t) = \int_{-\infty}^{\infty}x^*(\tau) h(t+\tau)d\tau z(t)=∫−∞∞x∗(τ)h(t+τ)dτ - 定理
f ( x , y ) ∞ g ( x , y ) < = > F ∗ ( u , v ) G ( u , v ) f(x,y)\quad \infty \quad g(x,y)<=>F^*(u,v)G(u,v) f(x,y)∞g(x,y)<=>F∗(u,v)G(u,v)
0.2.11. Rayleigh 定理
能量变换
对于有限区间非零函数 f(t), 其能量为
E
=
∫
−
∞
∞
∣
f
(
t
)
∣
2
d
t
E = \int_{-\infty}^{\infty}|f(t)|^2dt
E=∫−∞∞∣f(t)∣2dt
其变换函数与原函数有相同的能量
∫
−
∞
∞
∣
f
(
t
)
∣
2
d
t
=
∫
−
∞
∞
∣
F
(
u
)
∣
2
d
t
\int_{-\infty}^{\infty}|f(t)|^2dt = \int_{-\infty}^{\infty}|F(u)|^2dt
∫−∞∞∣f(t)∣2dt=∫−∞∞∣F(u)∣2dt
0.3. 快速傅里叶变换
由上面离散傅里叶变换的性质易知,直接计算 1维 dft 的时间复杂度维 O ( N 2 ) O(N^2) O(N2)。
利用到单位根的对称性,快速傅里叶变换可以达到 O ( n l o g n ) O(nlogn) O(nlogn)的时间复杂度。
0.3.1. 复数中的单位根
我们知道, 在复平面,复数
c
o
s
θ
+
i
s
i
n
θ
cos\theta +i\ sin\theta
cosθ+i sinθk可以表示成
e
i
θ
e^{i\theta}
eiθ, 可以对应一个向量。
θ
\theta
θ即为幅角。
在单位圆中 ,单位圆被分成
2
π
θ
\frac{2\pi}{\theta}
θ2π 份, 由单位圆的对称性
e
i
θ
=
e
i
(
θ
+
2
π
)
e^{i\theta} = e^{i(\theta+2\pi)}
eiθ=ei(θ+2π)
现在记 $ n =\frac{ 2\pi }{\theta}$ , 即被分成 n 份,幅度角为正且最小的向量称为 n 次单位向量, 记为
ω
n
\omega _n
ωn,
其余的 n-1 个向量分别为
ω
n
2
,
ω
n
3
,
…
,
ω
n
n
\omega_{n}^{2},\omega_{n}^{3},\ldots,\omega_{n}^{n}
ωn2,ωn3,…,ωnn ,它们可以由复数之间的乘法得来 $w_{n}{k}=w_{n}{k-1}\cdot w_{n}^{1}\ (2 \leq k \leq n) $。
单位根的性质
- 这个可以用 e 表示出来证明
ω 2 n 2 k = ω n k \omega_{2n}^{2k}=\omega_{n}^{k} ω2n2k=ωnk - 可以写成三角函数证明
ω n k + n 2 = − ω n k \omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^{k} ωnk+2n=−ωnk
容易看出 $w_{n}{n}=w_{n}{0}=1 $。
对于$ w_{n}^{k}$ , 它事实上就是 e 2 π i n k e^{\frac{2\pi i}{n}k} en2πik 。
0.3.2. 快速傅里叶变换的计算
下面的推导假设 n = 2 k n=2^k n=2k,以及代码实现 FFT 部分也是 如此。
利用上面的对称性,
将傅里叶计算进行奇偶分组
F
(
u
)
=
∑
i
=
0
n
−
1
ω
n
i
u
a
i
=
∑
i
=
0
n
2
−
1
ω
n
2
i
u
a
2
i
+
∑
i
=
0
n
2
−
1
ω
n
(
2
i
+
1
)
u
a
2
i
+
1
=
∑
i
=
0
n
2
−
1
ω
n
2
i
u
a
2
i
+
ω
n
u
∑
i
=
0
n
2
−
1
ω
n
2
i
u
a
2
i
+
1
=
F
e
v
e
n
(
u
)
+
ω
n
u
F
o
d
d
(
u
)
\begin{aligned} F(u)&=\sum_{i=0}^{n-1}\omega_n ^{iu} a^i\\ &= \sum_{i=0}^{\frac{n}{2}-1}\omega_n ^{2iu} a^{2i}+\sum_{i=0}^{\frac{n}{2}-1}\omega_n ^{(2i+1)u} a^{2i+1}\\ &=\sum_{i=0}^{\frac{n}{2}-1}\omega_{\frac{n}{2}} ^{iu} a^{2i}+\omega_n^u\sum_{i=0}^{\frac{n}{2}-1}\omega_{\frac{n}{2}} ^{iu} a^{2i+1}\\ & = F_{even}(u)+\omega_n^u F_{odd}(u) \end{aligned}
F(u)=i=0∑n−1ωniuai=i=0∑2n−1ωn2iua2i+i=0∑2n−1ωn(2i+1)ua2i+1=i=0∑2n−1ω2niua2i+ωnui=0∑2n−1ω2niua2i+1=Feven(u)+ωnuFodd(u)
F
e
v
e
n
F_{even}
Feven表示将 输入的次序中偶数点进行 Fourier 变换,
F
o
d
d
F_{odd}
Fodd 同理,这样就形成递推公式。
现在还没有减少计算量,下面通过将分别计算的 奇项,偶项利用起来,只计算 前
n
2
−
1
\frac{n}{2}-1
2n−1项,后面的一半可以利用此结果马上算出来。每一次可以减少一半的计算量。
对于
n
2
≤
i
+
n
2
≤
n
−
1
\frac{n}{2}\leq i+\frac{n}{2}\leq n-1
2n≤i+2n≤n−1
F
(
ω
n
i
+
n
2
)
=
F
e
v
e
n
(
ω
n
2
i
+
n
)
+
ω
n
i
+
n
2
⋅
F
o
d
d
(
ω
n
2
i
+
n
)
=
F
e
v
e
n
(
ω
n
2
i
+
n
2
)
+
ω
n
2
i
+
n
2
⋅
F
o
d
d
(
ω
n
2
i
+
n
2
)
=
F
e
v
e
n
(
ω
n
2
i
)
−
ω
n
2
i
⋅
F
o
d
d
(
ω
n
2
i
)
\begin{aligned} F(\omega_{n}^{i+\frac{n}{2}})&=F_{even}(\omega_{n}^{2i+n})+\omega_{n}^{i+\frac{n}{2}}\cdot F_{odd}(\omega_{n}^{2i+n})\\ &=F_{even}(\omega_{\frac{n}{2}}^{i+\frac{n}{2}})+\omega_{\frac{n}{2}}^{i+\frac{n}{2}}\cdot F_{odd}(\omega_{\frac{n}{2}}^{i+\frac{n}{2}})\\ & =F_{even}(\omega_{\frac{n}{2}}^{i})-\omega_{\frac{n}{2}}^{i}\cdot F_{odd}(\omega_{\frac{n}{2}}^{i}) \end{aligned}
F(ωni+2n)=Feven(ωn2i+n)+ωni+2n⋅Fodd(ωn2i+n)=Feven(ω2ni+2n)+ω2ni+2n⋅Fodd(ω2ni+2n)=Feven(ω2ni)−ω2ni⋅Fodd(ω2ni)
现在很清楚了,在每次计算 a[0…n-1] 的傅里叶变换F[0…n-1],分别计算出奇 odd[0…n/2-1],偶even[0…n/2-1](可以递归地进行),
那么傅里叶变换为:
F
[
i
]
=
{
e
v
e
n
[
i
]
+
ω
i
⋅
o
d
d
[
i
]
,
i
<
n
2
e
v
e
n
[
i
]
−
ω
i
⋅
o
d
d
[
i
]
,
e
l
s
e
F[i] = \begin{cases} even[i]+ \omega^i \cdot odd[i], \quad i<\frac{n}{2}\\ even[i]- \omega^i \cdot odd[i], \quad else \end{cases}
F[i]={even[i]+ωi⋅odd[i],i<2neven[i]−ωi⋅odd[i],else
0.4. 代码
下面是 python 实现
一维用 FFT 实现, 不过 只实现了 2 的幂。/ 对于非 2 的幂,用 FFT 实现有点困难,还需要插值,所以我 用
O
(
n
2
)
O(n^2)
O(n2) 直接实现。
二维的 DFT利用 分离性,直接调用 一维 FFT。
GitHub
import numpy as np
def _fft(a, invert=False):
N = len(a)
if N == 1:
return [a[0]]
elif N & (N - 1) == 0: # O(nlogn), 2^k
even = _fft(a[::2], invert)
odd = _fft(a[1::2], invert)
i = 2j if invert else -2j
factor = np.exp(i * np.pi * np.arange(N // 2) / N)
prod = factor * odd
return np.concatenate([even + prod, even - prod])
else: # O(n^2)
w = np.arange(N)
i = 2j if invert else -2j
m = w.reshape((N, 1)) * w
W = np.exp(m * i * np.pi / N)
return np.concatenate(np.dot(W, a.reshape(
(N, 1)))) # important, cannot use *
def fft(a):
'''fourier[a]'''
n = len(a)
if n == 0:
raise Exception("[Error]: Invalid length: 0")
return _fft(a)
def ifft(a):
'''invert fourier[a]'''
n = len(a)
if n == 0:
raise Exception("[Error]: Invalid length: 0")
return _fft(a, True) / n
def fft2(arr):
return np.apply_along_axis(fft, 0,
np.apply_along_axis(fft, 1, np.asarray(arr)))
def ifft2(arr):
return np.apply_along_axis(ifft, 0,
np.apply_along_axis(ifft, 1, np.asarray(arr)))
def test(n=128):
print('\nsequence length:', n)
print('fft')
li = np.random.random(n)
print(np.allclose(fft(li), np.fft.fft(li)))
print('ifft')
li = np.random.random(n)
print(np.allclose(ifft(li), np.fft.ifft(li)))
print('fft2')
li = np.random.random(n * n).reshape((n, n))
print(np.allclose(fft2(li), np.fft.fft2(li)))
print('ifft2')
li = np.random.random(n * n).reshape((n, n))
print(np.allclose(ifft2(li), np.fft.ifft2(li)))
if __name__ == '__main__':
for i in range(1, 3):
test(i * 16)