# 利用傅里叶变换去除下面图像中的条纹

lena_horiz_lines.bmp

# 实验思路 (含分析和结果):

# 图像数据输入并显示

close all;

image = imread('lena_corrupt.bmp');% 图像写入

subplot(1,4,2); imshow(image);

# 图像数据空域转频域，观察噪声类型设计滤波器

temp = im2gray(image);% 灰度转换

gray = double(temp);% 转为 double

F = fft2(gray);频谱转换

F1 = fftshift(F);% 频谱中心化处理

spectrum = log(1+abs(F1)); % 将复数变成实数并压缩数据范围

spectrum1 = uint8(mat2gray(spectrum)*255);% 标准化使范围在 [0.255]

subplot(1,4,1); imshow(spectrum1);

含噪声图像频谱图

原始图像无噪声频谱图是只有原点中心一个亮点的，除中心外其余亮点均为噪声频谱，且为周期噪声并分布在 y 轴上，所以需要设计一个垂直方向的陷波滤波器

# 设计并实现滤波器

[M,N] = size(image);% 图像大小

AV = 0; AH = 1;% 水平轴和垂直轴设置

SV = 30;SH = 30;% 水平轴和垂直轴长度设置

W = 61;% 带宽

H = ones(M,N,'single');

UC = floor(M/2)+1;

VC = floor(N/2)+1;

WL = (W-1)/2;% 折半

H(UC-WL:UC+WL,1:VC-SH) = AH;% 左边 X 轴

H(UC-WL:UC+WL,VC+SH:N) = AH;% 右边 X 轴

H(1:UC-SV,VC-WL:VC+WL) = AV;% 上边 Y 轴

H(UC+SV:M,VC-WL:VC+WL) = AV;% 下边 Y 轴

subplot(1,4,3);imshow(H);

陷波滤波器的一般形式
$H_{NR}(u,v) = \prod_{k=1}^{Q}H_k(u,v)H_{-k}(u,v)$,
$H_k(u,v)$ 和$H_{-k}(u,v)$ 是以中心为$(u_k,v_k)$ 和$(u_{-k},v_{-k})$ 的高通滤波器传递函数，它可以是理想高通滤波器、高斯高通滤波器和巴特沃斯高通滤波器。
其中$(u_k,v_k)$ 和$(u_{-k},v_{-k})$ 由相对频率矩形中心$\left [ floor(M/2),floor(N/2) \right ]$ 决定的
以理想高通滤波器举例
$H_{NR}(u,v)=\begin{cases} 0 & D(u,v)<=D0 \\ 1 & D(u,v)>D0 \end{cases}$
其中$D_k(u,v) = \left [(u-M/2-u_k)^2+(u-N/2-v_k)^2\right ]^{1/2}$,
$D_{-k}(u,v) = \left [(u-M/2+u_k)^2+(u-N/2+v_k)^2\right ]^{1/2}$,
$D_k(u,v)$ 和$D_{-k}(u,v)$ 是传递函数中心$(u,v)$ 到频率矩形中心的距离

理想带阻陷波滤波器频谱图

# 图像频域滤波后还原图像

result = ifftshift(H.*F1);

G = ifft2(result);

out=abs(G);

out=out/max(out(:));% 归一化 [0,1]

subplot(1,4,4);imshow(out);

一开始读取到的图像数据是整数型的灰度图像，值域在$[0,255]$, 经过滤波后得到的图像数据是浮点数型，但值域仍在$[0,255]$, 不属于灰度图像，需要经过归一化把值域压缩至$[0,1]$ 转换成二值图像 (浮点型) 才能显示图像，否则显示的只是一张白底，归一化后图像数据并没有发生变化。想了解不同图像类型的区别可戳此

滤波后图像频谱图

滤波后图像

实际最终得出的图像并没有完全去噪，放大后仍有条纹，由于直接$Y$ 轴放置滤波器，所以失去了很多图像细节

# 滤除下图的干扰信号

lena corrupt.bmp

# 实验思路 (含分析和结果):

# 图像数据输入并显示

close all;

image = imread('lena_corrupt.bmp');% 图像写入

subplot(1,4,2); imshow(image);

# 图像数据空域转频域，观察噪声类型设计滤波器

temp = im2gray(image);% 若是彩色图，返回灰度图，否则原样返回

gray = double(temp);% 转为 double

F = fft2(gray);

F1 = fftshift(F);% 频谱中心化处理

spectrum = log(1+abs(F1)); % 将复数变成实数并压缩数据范围

spectrum1 = uint8(mat2gray(spectrum)*255);% 标准化使范围在 [0.255]

subplot(1,4,1); imshow(spectrum1);

含噪声图像频谱图

同样是周期噪声，只是分布在中点周围，且很小，所以高斯带阻滤波器$(GBRF)$ 比较适合

# 设计并实现滤波器

[M,N] = size(image);% 图像大小

fbrf=ones(M,N);

W = 39;% 带宽

C0 = 170;% 带阻中心

for i=1:M

    for j=1:N

        D2 = (i-M/2)^2+(j-N/2)^2;

        D = sqrt((i-M/2)^2+(j-N/2)^2);

        fbrf(i,j)=1-exp(-((D2-C0^2)/(D*W))^2);

    end

end

 H=fbrf;

subplot(1,4,3);imshow(H);

高斯带阻滤波器$(GBRF)$ 的一般形式
$H(u,v) = 1-exp[\frac{D^2(u,v)-C_0^2}{D(u,v)W} ]$, 其中$D(u,v) = [(u-M/2)^2+(v-N/2)^2]^{1/2}$,
$C_0^2$ 是频带中心，$W$ 是带宽，$D(u,v)$ 是传递函数中心$(u,v)$ 到频率矩形中心的距离

高斯带阻滤波器频谱图

# 图像频域滤波后还原图像

result = ifftshift(H.*F1);

G = ifft2(result);

out=abs(G);

out=out/max(out(:));% 归一化 [0,1]

subplot(1,4,4);imshow(out);

滤波后图像频谱图

滤波后图像

因为这次噪声比较小，去噪效果比较好

源代码在 BigJod_1 文件夹内