卷積與濾波 -- Convolution and Filtering
在《
卷積與線性時間不變系統--Convolution and LTI Systems
》中，我們展示了用矩形脈衝對訊號進行卷積相當於對訊號應用
movmean
運算。
現在觀察當正弦波的組合與矩形脈衝訊號進行卷積時，頻域的結果。
% Time domain -> convolution
fs
=
10000
；
% Sampling frequency
n
=
-
10
：
1
/
fs
：
10
；
% Discretized time
x
=
sin
（
n
）
+
sin
（
10
*
n
）
+
cos
（
5
*
n
）；
% Input signal is a combination of sinusoids at different frequencies
h
=
rectpuls
（
n
）；
% Impulse response of an LTI system/operation
y
=
conv
（
h
，
x
，
“
same
”
）
。/
fs
；
% Convolution of h（n） and x（n）
% Frequency domain -> multiplication
L
=
2
^
nextpow2
（
length
（
n
））；
% Make length of fft to be the nearest next power of 2 for computation efficiency
f
=
0
：
fs
/
L
：
fs
-
（
fs
/
L
）；
% Discretized frequencies
X
=
abs
（
fft
（
x
，
L
）
。/
fs
）；
% Magnitude of FFT of x
H
=
abs
（
fft
（
h
，
L
）
。/
fs
）；
% Magnitude of FFT of the impulse response h
Y
=
X
。*
H
；
% Convolution theorem
plotSignals
（
n
，
h
，
x
，
y
，
f
，
H
，
X
，
Y
）；
plotSignals定義，在時域中視覺化卷積，在頻域中視覺化相應訊號的功能。
function
plotSignals
（
n，hn，xn，yn，f，hf，xf，yf
）
figure
（
“
Position
”
，［
0
0
1000
350
］，
“
Color
”
，［
0。9
0。9
0。9
］）
tiledlayout
（
1
，
2
，
“
TileSpacing
”
，
“
compact
”
）；
nexttile
plot
（
n
，
hn
，
“
b
——
”
，
“
LineWidth
”
，
2
）；
hold
on
plot
（
n
，
xn
，
“
r
”
，
“
LineWidth
”
，
1。5
）；
plot
（
n
，
yn
，
“
g
”
，
“
LineWidth
”
，
1。8
）；
hold
off
xlim
（［
-
3
3
］）
ylim
（［
-
3
3
］）
xlabel
（
“
Time
”
）
ylabel
（
“
Amplitude
”
）
title
（
“
Time
domain
”
）
legend
（
“
h
（
n
）
”
，
“
x
（
n
）
”
，
“
y
（
n
）
”
，
“
Location
”
，
“
northwest
”
）
nexttile
；
plot
（
f
，
hf
，
“
b
——
”
，
“
LineWidth
”
，
2
）；
hold
on
plot
（
f
，
xf
，
“
r
”
，
“
LineWidth
”
，
1。5
）；
plot
（
f
，
yf
，
“
g
”
，
“
LineWidth
”
，
1。8
）
hold
off
xlim
（［
0
，
4
］）
ylim
（［
0
，
2
］）
xlabel
（
“
Frequency
”
）
ylabel
（
“
Magnitude
”
）
title
（
“
Frequency
domain
”
）
legend
（
“
H
（
f
）
”
，
“
X
（
f
）
”
，
“
Y
（
f
）
”
，
“
Location
”
，
“
northeast
”
）
end
觀察與思考
在這個例子中，你可以認為
是一個濾波器的脈衝響應。簡單的說，濾波就是去除不理想的成分/特徵。movmean操作可以被認為是應用了一個“低通濾波器”，其脈衝響應是一個矩形脈衝。這可以去除較高頻率的成分。
加權平均法 #FormatImgID_4# 低通濾波法
將
作為濾波後的輸入訊號，
作為輸出。當是
是矩形脈衝時，在頻域中看到一個低通濾波器。如果我們希望在輸出中只保留低頻分量，那麼（藍色虛線）的理想頻率響應應該小於0。5。
二維卷積--空間濾波 Convolution in 2-D - Spatial Filtering
設計、實現或分析濾波器有兩種思路。
由時域濾波器卷積，即所需的脈衝響應。
乘以頻域濾波器，即所需的頻率響應。
我們重點介紹第一種實現濾波器的方式。這被廣泛用於處理數字影象和數字音訊。由於影象是二維（2-D）離散訊號，我們首先定義一個離散的二維卷積，它只是將一維求和擴充套件到二維。
對於影象來說，2維是空間的（空間的），而不是一維訊號中通常遇到的時間的（時間的）。因此，我們將使用2-D卷積，透過2-D平均濾波卷積來模糊這張花的影象。
加權平均
空間模糊
低通濾波
使用imread讀取影象，使用imshow顯示影象。
% Read image， visualize it， and display its size and class。
I
=
imread
（
“
flower
。
jpg
”
）；
% Convert to double precision
I
=
im2double
（
I
）；
% Display
figure
，
imshow
（
I
）
whos
I
注：在MATLAB中，表示影象的double型矩陣預計在
［0，1］
範圍內有值。 像這樣的彩色影象是三個二維訊號的組合——分別為R，G，B。因此，我們將分別進行三次二維卷積，以模糊三個顏色強度矩陣。
% Split the image into the 3 component color intensity planes or matrices。
［
R
，
G
，
B
］
=
imsplit
（
I
）；
% Create a 2-D convolution filter or kernel for averaging or smoothing
W
=
5
；
h
=
ones
（
W
）
。/
W
。^
2
% Convolve R with h to blur the image
Ravg
=
conv2
（
R
，
h
，
“
same
”
）；
% Convolve G with h to blur the image
Gavg
=
conv2
（
G
，
h
，
“
same
”
）；
% Convolve B with h to blur the image
Bavg
=
conv2
（
B
，
h
，
“
same
”
）；
% Concatenate the three matrices to put together the color image
Iavg
=
cat
（
3
，
Ravg
，
Gavg
，
Bavg
）；
whos
Iavg
imfilter
這個功能可以讓過濾圖片的過程變得更加簡單。你可以透過使用
imfilter
達到和我們上面卷積中一樣的結果。
Ifilt
=
imfilter
（
I
，
h
，
“
conv
”
）；
whos
Ifilt
imshow
（
Ifilt
，［］）；
觀察與思考
噪聲的產生是因為影象中的某些畫素的灰度值發生了突變，使得和周圍區域不和諧。除噪其實去除高頻噪聲，使得影象中的噪聲畫素的灰度值不那麼突兀。
噪聲去除有高斯濾波，均值濾波，中值濾波等。從頻率域觀點來看這些濾波器是一種低通濾波器，高頻訊號將會去掉，因此可以幫助消除影象尖銳噪聲，實現影象平滑，模糊等功能。
上文中的 W 是一個均值濾波器。增加濾波器的大小， 會使得圖片更模糊
使用進行均值濾波操作來模糊影象。輸出影象的每一個畫素灰度值是卷積核在輸入影象中對應的畫素的平均值（ 所有畫素加權係數相等）。 均值濾波卷積核所覆蓋的九個畫素點具有同樣權重，該卷積核的作用在於取九個值的平均值代替中間畫素值，所以起到的平滑的效果。
相比於高斯濾波，它不能很好地保護影象細節，在影象去噪的同時也破壞了影象的細節部分，丟失了影象本身的一些屬性，從而使影象變得模糊，不能很好地去除噪聲點。
更多空間濾波器 -- More Spatial Filters
kernel 或 mask（h） 中的值決定了要保留的影象的空間成分/特徵。kernel 可以包含任何值。然而，有一些kernel 或 mask，是經常使用的。讓我們生成一些常用的 mask，並觀察它們對影象的影響。
常用的 mask 有 Prewitt 運算元，Sobel 運算元，Laplacian 運算元
Prewitt 運算元和Sobel 運算元都是透過求一階導數來計算梯度的，用於線的檢測，通常用於邊緣檢測。在影象處理過程中，除了檢測線，有時候也需要檢測特殊點，這就需要用二階導數進行檢測，著名的就是拉普拉斯（Laplacian）運算元。
對影象求兩次導數，公式如下：
下面程式碼可以測試不同運算元偏遠的
% Read image from file
I
=
imread
（
“
flower
。
jpg
”
）；
% Define 3 different filter kernels
h1
=
［
1
2
1
；
% Sobel operator， calculate vertical gradients and detect horizontal edges
0
0
0
；
-
1
-
2
-
1
］；
h2
=
［
1
0
-
1
；
% Sobel operator calculate horizontal gradients and detect vertical edges
2
0
-
2
；
1
0
-
1
］；
h3
=
［
0
-
1
0
；
-
1
5
-
1
；
% Laplace Operator detect edges
0
-
1
0
］；
h4
=
［
-
1
-
1
-
1
；
-
1
9
-
1
；
-
1
-
1
-
1
］；
%% Laplace Operator detect edges
% Apply filters
Ifilt
=
imfilter
（
I
，
h3
，
“
conv
”
）；
% Display the original image and the filtered image side-by-side
figure
；
imshowpair
（
I
，
Ifilt
，
“
montage
”
）；
拉普拉斯運算元在邊緣檢測的應用中並不侷限於水平方向或垂直方向， 這是Laplacian 與Soble 的區別。
因為一階二階導數都能放大孤立點和孤立線（噪聲）的影響，所以如果存在噪聲，那麼一階二階導數處理過後的影象將會有更多更大的噪聲。所以對影象進行一階二階導數運算之前需要先對影象做平滑去噪處理。
觀察與思考
卷積 -- Convolution
影象的線性濾波是透過一種叫做卷積的操作來完成的。卷積是一種鄰域操作，其中每個輸出畫素是相鄰輸入畫素的加權和。權重矩陣稱為卷積核，也稱為濾波器。卷積核是一個旋轉了180度的相關核。
例如，假設影象是
![A=\left[\begin{array}{rrrrr} 17 & 24 & 1 & 8 & 15 \\ 23 & 5 & 7 & 14 & 16 \\ 4 & 6 & 13 & 20 & 22 \\ 10 & 12 & 19 & 21 & 3 \\ 11 & 18 & 25 & 2 & 9\end{array}\right] \\](https://img.heatask.com/upload/website_attach/yglfuZgUgfWIzd36PaJAhzAoGSaIbQElhMMf4Fy7+Rvqgn_Mh=cW-nOL1_OVjbtbvIA+pclr3Qo6jv4mFN=ezMN9aQ=K0__mgD-2u19Jqb6ALmUqku0CcJ=Zs8TaU+owRoTIEGJXDqj+QGWfi08_=cAJ6IpLcfJ7ZKzXeuxrh7NWLmUqku0Ccn=Zs8TaU+ocSdA+pm9+pf6eQmr1r4wjUcA4a7ptcfJ7ZKzXh3aXDqj+QGs1i08_=cAJ6IptcbJ7ZKzXhv=XDqj+QGWKi08e0+p4PrpKc7BwZdpretgre8YWXvWWh=tjUdL7tIAL0KJmmnprh31XhK6qX3HW8t8_=iOoHdAb=gpmgyM6Dvu9DqAryXkwBcnk0p=Zs8O6-qEn0.jpeg)
相關kernel是
![\mathrm{h}=\left[\begin{array}{rrr}8 & 1 & 6 \\ 3 & 5 & 7 \\ 4 & 9 & 2\end{array}\right] \\](https://img.heatask.com/upload/website_attach/CWnc_aN=kbyfjf0S_yIBh_E6CiYwOf8=PlXchABSLipoQpzxGqbxz0txyPekaduSi44q-sYVi7Gpia8yiZrj=6Ef5iSkHfUfi44LzBudALmo0cJxGhtf_ARqLgsRt0oM8GUEUGod8fxeUHoisdrj=+n2QcEea4Sa_9rH02Sj_b9d=HoisdU_U+E1L+nSt0oM8GUH-s8kB9-F-t1ZPmtfz0BqLgskaDYMCl4Uj9oi4Xxo0lNxGqykH.jpeg)
您可以使用以下步驟來計算位置
（2，4）
的輸出畫素。
將相關kernel圍繞其中心元素旋轉180度，建立一個卷積kernel。
滑動卷積kernel的中心元素，使其位於A的（2，4）元素之上。
將旋轉後的卷積核中的每個權重乘以下面A的畫素。
將步驟3中的各個乘積相加。
因此，
（2，4）
的輸出畫素為
如下圖所示。
透過
fft2
計算二維線性卷積，可以直觀地看到頻域中的濾波過程，理解影象中的頻率概念。這與透過
fft
計算一維訊號的卷積輸出類似。
% Import an image and display it。
figure
；
I
=
imread
（
“
flower
。
jpg
”
）；
% For ease of analysis， convert the image to grayscale before further processing。
I
=
rgb2gray
（
I
）；
imshow
（
I
）
% Create some filter kernels to test。 You can add your own to the list。
W
=
15
；
h1
=
ones
（
W
）
。/
W
。^
2
；
%% Mean Filter
h2
=
［
1
0
-
1
；
2
0
-
2
；
1
0
-
1
］；
% Sobel Operator
h3
=
［
1
2
1
；
0
0
0
；
-
1
-
2
-
1
］；
%%Sobel Operator
h4
=
［
-
1
-
1
-
1
；
-
1
9
-
1
；
-
1
-
1
-
1
］；
%% Laplace Operator detect edges
h
=
h2
根據影象的大小和 kernel 計算全線性卷積的大小。
% Get the sizes of the image and filter kernel。
［
MI
，
NI
］
=
size
（
I
）；
［
Mh
，
Nh
］
=
size
（
h
）；
% Compute size of full convolution output for each dimension。
m
=
MI
+
Mh
+
1
；
% rows
n
=
NI
+
Nh
+
1
；
% columns
我們使用
fft2
來計算影象 I 的二維傅立葉變換，大小為 m 行 n 列。
% Compute FFT of grayscale image。
F1
=
fft2
（
I
，
m
，
n
）；
我們使用
imagesc
來視覺化 2-D FFT的輸出，用縮放的顏色來表示FFT的幅度。
% Display FFT output using a colormap and log scale for better visualization。
imagesc
（
fftshift
（
log
（
1
+
abs
（
F1
））））
colorbar
使用
fft2
計算濾波器 h 的二維傅立葉變換，大小為 m 行 n 列。
% Compute FFT of filter kernel。
F2
=
fft2
（
h
，
m
，
n
）；
% Display FFT output using a colormap and log scale for better visualization。
imagesc
（
fftshift
（
log
（
1
+
abs
（
F2
））））
colorbar
將 F1 和 F2 元素相乘。
% Convolution in spatial domain is multiplication in frequency domain。
F
=
F1
。*
F2
；
% Display FFT output using a colormap and log scale for better
% visualization。
imagesc
（
fftshift
（
log
（
1
+
abs
（
F
））））
colorbar
使用
ifft2
計算乘積F的逆傅立葉變換。使用
ifft2
。這將使我們得到濾波後的影象。
% Compute inverse FFT of the product。
Ifilt
=
ifft2
（
F
，
m
，
n
）；
將影象修剪到合適的大小，以實現相當於線性卷積的效果。將影象轉換為無符號8位整數，並檢視輸出。
Ifilt
=
Ifilt
（
ceil
（
Mh
/
2
）：
end
-
（
ceil
（
Mh
/
2
）
+
1
），
ceil
（
Nh
/
2
）：
end
-
（
ceil
（
Nh
/
2
）
+
1
））；
Ifilt
=
uint8
（
real
（
Ifilt
））；
imshow
（
Ifilt
）
使用
imfilter
對 I 和 h 進行等效線性空間卷積。驗證發現2者一致。
Ifilt2
=
imfilter
（
I
，
h
，
“
conv
”
）；
imshow
（
Ifilt2
）
我們先將 I 和 h 透過
fft2
轉為頻域，將頻域結果點乘，再轉為時域。與直接將 I 與 h 做卷積的結果一致。
進一步探索
影象的空間過濾通常透過應用相關操作來進行
模板匹配--Template Matching
在這個例子中，我們將計算字母 “t ”在這段文字中出現的次數。
文字圖片
% Read an input image and a template image to match。
figure
；
I
=
imread
（
“
text
。
png
”
）；
imshow
（
I
）
title
（
“
Input
text
image
”
）
template
=
imread
（
“
lettert
。
png
”
）；
imshow
（
template
）
title
（
“
Template
to
match
”
）
% Convert the image and the template to double precision floating point to
% perform further computations。
I
=
im2double
（
I
）；
template
=
im2double
（
template
）；
result
=
imfilter
（
I
，
template
，
“
corr
”
）；
% Rescale image between 0 and 1， visualize the output。
result
=
rescale
（
result
）；
imshow
（
result
）
title
（
“
Result
of
applying
filtering
operation
”
）
% Locate the peak responses indicating a match and count the number of
% times it occurs。
t
=
0。95
；
imshow
（
result
>
t
）
disp
（
“
Number
of
standing
t
‘
s
in
the
text
=
”
+
nnz
（
result
>
t
））
結果顯示圖片中有5個匹配的t。
結論
對訊號進行濾波時，等價於與時域濾波器卷積，得到時域輸出。
也可以將訊號用
fft
轉為頻域，乘以頻域濾波器，再用
ifft
轉為時域。
圖片除噪其實去除高頻噪聲，保留低頻訊號，使得影象中的噪聲畫素的灰度值不那麼突兀。
上一篇:全職媽媽，能做什麼副業？
下一篇：關於腎的小知識
