实验要求：

Objective:
To understand the principle of the notch filter and its periodic noise reducing ability.
Main requirements:
Ability of programming with C, C++, or Matlab.
Instruction manual:
(a) Write a program that implements sinusoidal noise of the form given in Problem 5.14. The inputs to the program must be the amplitude, A, and the two frequency components u0 and v0 shown in the problem equation.
(b) Download image 5.26(a) and add sinusoidal noise to it, with u0 = M/2 (the image is square) and v0 = 0. The value of A must be high enough for the noise to be quite visible in the image.
(c) Compute and display the spectrum of the image. If the FFT program you developed in Project 4.01 can only handle images of size equal to an integer power of 2, reduce the size of the image to 512 x 512 or 256 x 256 using the program from Project 02-04. Resize the image before adding noise to it.
(d) Notch-filter the image using a notch filter of the form shown in Fig. 5.19(c).

这个实验的要求是使用陷波滤波器消除周期噪声。
对于某种周期噪声常常是以一对共轭的点的形式，出现在图像的频谱中，所以为我们使用一个陷波滤波器，使其正好对应那对共轭点以消除周期噪声。

实验代码：

%
close all;
clc;
clear all;

%
img = imread('Fig5.26(a).jpg');
[M, N] = size(img);
figure;
subplot(2,3,1);
imshow(img);
title('原图像');

% 产生周期噪声
C = [343, 0];
[noise, R, S] = imnoise3(M, N, C, 40000000);
subplot(2,3,2);
imshow(noise, []);
title('噪声');

% 计算傅里叶变换
img1 = double(img);
img2 = img1 + noise;
img3 = fft2(img2);
subplot(2,3,3);
imshow(img3);
title('傅里叶变换');

% 生成滤波器
sig = 100;
H = lpfilter('gaussian', M, N,sig);
subplot(2,3,4);
imshow(H);
title('滤波器');

% 频率域滤波
img_G = H .* img3;
subplot(2,3,5);
imshow(img_G);
title('频域陷波滤波');

% 傅里叶逆变换还原图像到时域
img_g = real(ifft2(img_G));
img_g = mat2gray(img_g);
subplot(2,3,6);
imshow(img_g);
title('陷波滤波结果');

实验结果：