背景知识
- 线性系统
X 1 ( t ) − > Y 1 ( t ) X 2 ( t ) − > Y 2 ( t ) X_1(t) -> Y_1(t)\\ X_2(t)->Y_2(t) X1(t)−>Y1(t)X2(t)−>Y2(t)
则:
X
1
(
t
)
+
X
2
(
t
)
−
>
Y
1
(
t
)
+
Y
2
(
t
)
a
X
1
(
t
)
−
>
a
Y
1
(
t
)
X_1(t)+X_2(t)->Y_1(t)+Y_2(t)\\ aX_1(t)->aY_1(t)
X1(t)+X2(t)−>Y1(t)+Y2(t)aX1(t)−>aY1(t)
- 平移不变性
对于某系统,有 x ( t ) − > y ( t ) x(t)->y(t) x(t)−>y(t) 当输入信号平移 T T T 时, x ( t − T ) − > y ( t − T ) x(t-T)->y(t-T) x(t−T)−>y(t−T) 则称该系统具有平移不变性。
- 卷积
对于一个线性系统的输入f(x)和输出g(x),他们之间必然存在这样的关系:
g ( x ) = ∫ − ∞ + ∞ f ( τ ) h ( x − τ ) g(x)=\int_{-\infin}^{+\infin}f(\tau)h(x-\tau) g(x)=∫−∞+∞f(τ)h(x−τ)
其中 h ( x ) h(x) h(x)是这个系统的单位冲激响响应。该积分也称为卷积积分,记为: g ( x ) = f ( x ) ∗ h ( x ) g(x)=f(x)*h(x) g(x)=f(x)∗h(x)
离散一维卷积:
f
(
x
)
∗
h
(
x
)
=
1
M
∑
m
=
0
M
−
1
f
(
m
)
h
(
x
−
m
)
f(x)*h(x)=\frac{1}{M}\sum_{m=0}^{M-1}f(m)h(x-m)
f(x)∗h(x)=M1m=0∑M−1f(m)h(x−m)
二维卷积:
f
(
x
,
y
)
∗
g
(
x
,
y
)
=
∫
−
∞
+
∞
f
(
u
,
v
)
h
(
x
−
u
,
y
−
v
)
d
u
d
v
f(x,y)*g(x,y) = \int_{-\infin}^{+\infin}f(u,v)h(x-u,y-v)dudv
f(x,y)∗g(x,y)=∫−∞+∞f(u,v)h(x−u,y−v)dudv
离散二维卷积:
f
(
x
,
y
)
∗
g
(
x
,
y
)
=
1
M
N
∑
m
=
0
M
−
1
∑
n
=
0
N
−
1
f
(
m
,
n
)
h
(
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)h(x-m,y-n)
f(x,y)∗g(x,y)=MN1m=0∑M−1n=0∑N−1f(m,n)h(x−m,y−n)
傅里叶变换
一维傅里叶变换:
F
(
u
)
=
∫
−
∞
+
∞
f
(
x
)
e
−
j
2
π
u
x
d
x
F(u)=\int_{-\infin}^{+\infin}f(x)e^{-j2\pi ux}dx
F(u)=∫−∞+∞f(x)e−j2πuxdx
逆变换:
f
(
x
)
=
∫
−
∞
+
∞
F
(
u
)
e
j
2
π
u
x
d
u
f(x)=\int_{-\infin}^{+\infin}F(u)e^{j2\pi ux}du
f(x)=∫−∞+∞F(u)ej2πuxdu
将f(x)分解为不同频率的成分,F(u)代表频率为u的成分的强度
函数f(x)的傅立叶变换后一般是一个复量,它可以用下式表示:
F
(
u
)
=
R
(
u
)
+
j
I
(
u
)
F
(
u
)
=
∣
F
(
u
)
∣
e
j
ϕ
(
u
)
∣
F
(
u
)
∣
=
R
2
(
u
)
+
I
2
(
u
)
,
f
(
x
)
的
幅
度
谱
ϕ
(
u
)
=
a
r
t
a
n
[
I
(
u
)
R
(
u
)
]
,
f
(
x
)
的
相
位
谱
P
(
u
)
=
∣
F
(
u
)
∣
2
=
P
2
(
u
)
+
I
2
(
u
)
,
f
(
x
)
的
功
率
谱
F(u)=R(u)+jI(u)\\ F(u) = |F(u)|e^{j\phi(u)}\\ |F(u)|=\sqrt{R^2(u)+I^2(u)} , f(x)的幅度谱\\ \phi(u)=artan[\frac{I(u)}{R(u)}], f(x)的相位谱\\ P(u)=|F(u)|^2=P^2(u)+I^2(u),f(x)的功率谱
F(u)=R(u)+jI(u)F(u)=∣F(u)∣ejϕ(u)∣F(u)∣=R2(u)+I2(u),f(x)的幅度谱ϕ(u)=artan[R(u)I(u)],f(x)的相位谱P(u)=∣F(u)∣2=P2(u)+I2(u),f(x)的功率谱
二维连续傅里叶变换及其逆变换:
F
(
u
,
v
)
=
∫
−
∞
+
∞
∫
−
∞
+
∞
f
(
x
,
y
)
e
−
j
2
π
(
u
x
+
v
y
)
d
x
d
y
F(u,v)=\int_{-\infin}^{+\infin}\int_{-\infin}^{+\infin}f(x,y)e^{-j2\pi (ux+vy)}dxdy
F(u,v)=∫−∞+∞∫−∞+∞f(x,y)e−j2π(ux+vy)dxdy
f ( x , y ) = ∫ − ∞ + ∞ ∫ − ∞ + ∞ F ( u , v ) e j 2 π ( u x + v y ) d u d v f(x,y)=\int_{-\infin}^{+\infin}\int_{-\infin}^{+\infin}F(u,v)e^{j2\pi (ux+vy)}dudv f(x,y)=∫−∞+∞∫−∞+∞F(u,v)ej2π(ux+vy)dudv
二维傅里叶变换的幅度谱和相位谱为:
∣
F
(
u
,
v
)
∣
=
R
2
(
u
,
v
)
+
I
2
(
u
,
v
)
,
幅
度
谱
ϕ
(
u
,
v
)
=
a
r
t
a
n
[
I
(
u
,
v
)
R
(
u
,
v
)
]
,
相
位
谱
|F(u,v)|=\sqrt{R^2(u,v)+I^2(u,v)} , 幅度谱\\ \phi(u,v)=artan[\frac{I(u,v)}{R(u,v)}], 相位谱\\
∣F(u,v)∣=R2(u,v)+I2(u,v),幅度谱ϕ(u,v)=artan[R(u,v)I(u,v)],相位谱
离散傅里叶变换 D F T DFT DFT (Discrete Fourier Transform)
一维离散傅立叶变换:
F
(
u
)
=
∑
k
=
0
N
−
1
f
(
k
)
e
−
j
2
π
u
k
N
,
u
=
0
,
1
,
2
,
.
.
.
,
N
−
1
f
(
k
)
=
1
N
∑
u
=
0
N
−
1
F
(
u
)
e
j
2
π
u
k
N
,
k
=
0
,
1
,
2
,
.
.
.
,
N
−
1
F(u)=\sum_{k=0}^{N-1}f(k)e^{-\frac{j2\pi uk}{N}},u=0,1,2,...,N-1\\ f(k)=\frac{1}{N}\sum_{u=0}^{N-1}F(u)e^{\frac{j2\pi uk}{N}},k=0,1,2,...,N-1\\
F(u)=k=0∑N−1f(k)e−Nj2πuk,u=0,1,2,...,N−1f(k)=N1u=0∑N−1F(u)eNj2πuk,k=0,1,2,...,N−1
二维离散傅立叶变换:
F
(
u
,
v
)
=
∑
x
=
0
M
−
1
∑
y
=
0
N
−
1
f
(
x
,
y
)
e
−
j
2
π
(
u
x
M
+
v
y
N
)
u
=
0
,
1
,
2
,
.
.
.
,
M
−
1
v
=
0
,
1
,
2
,
.
.
.
,
N
−
1
f
(
x
,
y
)
=
1
M
N
∑
u
=
0
M
−
1
∑
v
=
0
N
−
1
F
(
u
,
v
)
e
j
2
π
(
u
x
M
+
v
y
N
)
x
=
0
,
1
,
2
,
.
.
.
,
M
−
1
y
=
0
,
1
,
2
,
.
.
.
,
N
−
1
F(u,v)=\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)e^{-j2\pi(\frac{ ux}{M}+\frac{vy}{N})}\\ u=0,1,2,...,M-1\\ v=0,1,2,...,N-1\\ f(x,y)=\frac{1}{MN}\sum_{u=0}^{M-1}\sum_{v=0}^{N-1}F(u,v)e^{j2\pi(\frac{ ux}{M}+\frac{vy}{N})}\\ x=0,1,2,...,M-1\\ y=0,1,2,...,N-1\\
F(u,v)=x=0∑M−1y=0∑N−1f(x,y)e−j2π(Mux+Nvy)u=0,1,2,...,M−1v=0,1,2,...,N−1f(x,y)=MN1u=0∑M−1v=0∑N−1F(u,v)ej2π(Mux+Nvy)x=0,1,2,...,M−1y=0,1,2,...,N−1
二维离散傅立叶变换特性:
- 可分离性
- 平移性
- 线性
- 旋转特性
- 均值性
- 微分特性
图像傅立叶变换的物理意义是得到图像的不同频率成分的分布情况。
频率域平滑滤波器
思想:边缘和其它尖锐变化(如噪声)在傅立叶变换中对于高频成分,因此平滑可通过衰减图像傅立叶变换的高频成分,保留低频成分来实现。
理想低通滤波器
一个二维的理想低通过滤器(
I
L
P
F
ILPF
ILPF)的系统函数满足(是一个分段函数):
H
(
u
,
v
)
=
{
1
D
(
u
,
v
)
≤
D
0
0
D
(
u
,
v
)
>
D
0
D
0
为
截
止
频
率
,
D
(
u
,
v
)
为
距
离
函
数
D
(
u
,
v
)
=
(
u
2
+
v
2
)
1
2
H(u,v) = \begin{cases} 1&D(u,v)\leq D_0\\ 0&D(u,v)>D_0\\ \end{cases}\\ D_0 为截止频率,D(u,v)为距离函数 D(u,v)=(u^2+v^2)^\frac{1}{2}
H(u,v)={10D(u,v)≤D0D(u,v)>D0D0为截止频率,D(u,v)为距离函数D(u,v)=(u2+v2)21
D
0
D_0
D0 半径内的频率分量无损通过,圆外的频率分量会被滤除。若滤除的高频分量中含有大量的边缘信息,会发生图像边缘模糊现象。钝化的图像被一种非常严重的振铃效果(理想低通滤波器的一种特性)所影响。
Code
import matplotlib.pyplot as plt
import numpy as np
import cv2 as cv
# 二维离散傅里叶变换
"""
def fft2(img):
cp = np.copy(img).astype(np.complex)
temp = np.copy(img)
M,N = cp.shape
def fun(u,v,f):
value = 0
for x in range(M):
for y in range(N):
exp = np.exp(np.complex(0,-2*np.pi*(u*x/M+v*y/N)))
value += f[x][y]*exp
return value
for u in range(M):
for v in range(N):
cp[u][v] = fun(u,v,temp)
print(cp[u][v])
return cp
"""
# 理想低通滤波
def ilpf(f,D0):
cp = np.copy(f)
r,c = cp.shape
H = np.zeros((r,c))
for i in range(r):
for j in range(c):
if (i-r/2)**2+(j-c/2)**2<= D0**2:
H[i][j] = 1
return H*cp
def show(img,title,sub):
plt.subplot(sub[0],sub[1],sub[2])
plt.imshow(img,cmap="gray")
plt.axis('off') # 去掉坐标轴
plt.title(title)
if __name__=='__main__':
img = cv.imread('10.png',0)
# 傅里叶变换
f = np.fft.fftn(img)
# 坐标原点平移到中心
fshift = np.fft.fftshift(f)
# 避免频谱范围过大影响视觉效果
fimg = np.log(1+np.abs(fshift))
D = [5,15,30,80,230]
i = 1
for D0 in D:
# 傅里叶逆变换
ishift = np.fft.ifftshift(ilpf(fshift,D0))
re = np.fft.ifftn(ishift)
reimg = np.abs(re)
show(reimg,'D0=%d'%D0,(2,3,i))
i+=1
show(img,'Original img',(2,3,6))
plt.show()
B u t t e r w o r t h Butterworth Butterworth 低通滤波器
一个截止频率在与原点距离为 D 0 D0 D0 的 n n n 阶 B u t t e r w o r t h Butterworth Butterworth 低通过滤器( B L P F BLPF BLPF)的变换函数如下:
H
(
u
,
v
)
=
1
1
+
[
D
(
u
,
v
)
D
0
]
2
n
H(u,v)=\frac{1}{1+[\frac{D(u,v)}{D_0}]^{2n}}
H(u,v)=1+[D0D(u,v)]2n1
Code
if __name__=='__main__':
x = np.linspace(0,10,200)
H = []
for i in range(1,5):
H.append(1/(1+np.power(x/5.,2*i)))
plt.plot(x, H[0])
plt.plot(x, H[1])
plt.plot(x, H[2])
plt.plot(x, H[3])
plt.legend(['n=1','n=2','n=3','n=4'])
plt.grid()
plt.ylabel('H(u,v)')
plt.xlabel('D(u,v)')
plt.show()
# 巴特沃斯低通滤波器
def blpf(f,D0,n):
"""
:param f: 傅里叶变换后的幅度谱
:param D0: 截止距离
:param n: 阶数
:return: 新的幅度谱
"""
cp = np.copy(f)
r,c = cp.shape
H = np.zeros((r, c))
for i in range(r):
for j in range(c):
D = np.sqrt((i-r/2)**2+(j-c/2)**2)
H[i][j] = 1/(1+np.power(D/D0,2*n))
return H*cp
- 在任何经二阶 B L P F BLPF BLPF 处理过的图像中都没有明显的振铃效果,这是过滤器在低频和高频之间的平滑过渡的结果。
- 同样的截止频率下,比理想低通保留相对较多高频成分,所以模糊程度减少,但是去噪效果不如理想低通。
高斯低通滤波器
H ( u , v ) = e x p ( − D 2 ( u , v ) 2 D 0 2 ) H(u,v)=exp(-\frac{D^2(u,v)}{2D_0^2}) H(u,v)=exp(−2D02D2(u,v))
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-3WI5B2RL-1584523520945)(image-20200318154429995.png)]
高斯低通滤波器 ( G L P F GLPF GLPF) 有更加平滑的过渡带,比 B L P F BLPF BLPF 相比衰减更快。无振铃,比 B L F P BLFP BLFP 去噪效果好。
Code
def glpf(f,D0):
cp = np.copy(f)
r,c = cp.shape
H = np.zeros((r, c))
for i in range(r):
for j in range(c):
D = (i-r/2)**2+(j-c/2)**2
H[i][j] = np.exp(-D/(2*D0**2))
return H*cp
频率域锐化滤波器
基本思想:低通滤波器的反操作,即 H h p ( u , v ) = 1 − H l p ( u , v ) H_{hp}(u,v)=1-H_{lp}(u,v) Hhp(u,v)=1−Hlp(u,v)。与低通滤波器相对应,频率域内常用的高通滤波器有3种:
- 理想高通滤波器
D 0 = 15 D0=15 D0=15 和 D 0 = 30 D0=30 D0=30 时,振铃问题十分明显。
- 巴特沃斯高通滤波器
- 高斯高通滤波器
高通滤波器的问题:
- 低频成分被严重地消弱了,图像失去主体成分,只剩下细节。
- 改进措施
- 加一个常数到变换函数 H ( u , v ) + A H(u,v) + A H(u,v)+A 这种方法被称为 高频提升
- 为了解决变暗的趋势,在变换结果图像上再进行一次直方图均衡化。这种方法被称为 后过滤处理
高频提升滤波:
f
h
b
=
A
f
(
x
,
y
)
−
f
l
p
(
x
,
y
)
=
(
A
−
1
)
f
(
x
,
y
)
+
f
h
p
(
x
,
y
)
H
h
p
=
(
A
−
1
)
+
H
h
p
(
u
,
v
)
f_{hb}=Af(x,y)-f_{lp}(x,y)\\ =(A-1)f(x,y)+f_{hp}(x,y)\\ H_{hp}=(A-1)+H_{hp}(u,v)
fhb=Af(x,y)−flp(x,y)=(A−1)f(x,y)+fhp(x,y)Hhp=(A−1)+Hhp(u,v)
高频加强滤波:在高通滤波器函数前乘以系数
b
b
b,再加上
a
a
a 倍的原始图像。
H
h
p
e
(
u
,
v
)
=
a
+
b
H
h
p
(
u
,
v
)
,
a
≥
0
,
b
>
a
当
b
=
1
时
,
高
频
加
强
转
换
为
高
频
提
升
滤
波
H_{hpe}(u,v)=a+bH_{hp}(u,v), a\geq0,b>a\\ 当b=1时,高频加强转换为高频提升滤波
Hhpe(u,v)=a+bHhp(u,v),a≥0,b>a当b=1时,高频加强转换为高频提升滤波
- 改进措施
- 加一个常数到变换函数 H ( u , v ) + A H(u,v) + A H(u,v)+A 这种方法被称为 高频提升
- 为了解决变暗的趋势,在变换结果图像上再进行一次直方图均衡化。这种方法被称为 后过滤处理
高频提升滤波:
f
h
b
=
A
f
(
x
,
y
)
−
f
l
p
(
x
,
y
)
=
(
A
−
1
)
f
(
x
,
y
)
+
f
h
p
(
x
,
y
)
H
h
p
=
(
A
−
1
)
+
H
h
p
(
u
,
v
)
f_{hb}=Af(x,y)-f_{lp}(x,y)\\ =(A-1)f(x,y)+f_{hp}(x,y)\\ H_{hp}=(A-1)+H_{hp}(u,v)
fhb=Af(x,y)−flp(x,y)=(A−1)f(x,y)+fhp(x,y)Hhp=(A−1)+Hhp(u,v)
高频加强滤波:在高通滤波器函数前乘以系数
b
b
b,再加上
a
a
a 倍的原始图像。
H
h
p
e
(
u
,
v
)
=
a
+
b
H
h
p
(
u
,
v
)
,
a
≥
0
,
b
>
a
当
b
=
1
时
,
高
频
加
强
转
换
为
高
频
提
升
滤
波
H_{hpe}(u,v)=a+bH_{hp}(u,v), a\geq0,b>a\\ 当b=1时,高频加强转换为高频提升滤波
Hhpe(u,v)=a+bHhp(u,v),a≥0,b>a当b=1时,高频加强转换为高频提升滤波