attenuating high frequency by Fourier transform

Fourier平面距$(0,0)$较远的点是高噪声，这Fourier变换的形式可以看出
$$F(k_1,k_2)=\int f(x,y)exp(-i{k}_{1}x-i{k}_{2}y)$$较大的$(k_1,k_2)$是$f$向较大的$sinkx和coskx$投影。
而选取$(0,0)$的邻域进行Fourier逆变换自然是达到去噪的目的。

额外需要指出的是，由Riemann-Lebesgue引理$\underset{n\to \infty }{\lim }\int f(x)e^{inx}dx=0$，高噪声的部分系数是偏小的。

f = imread('pic.jpg');
f = double(f);
s = size(f); N = s(1); M = s(2);

K = 75;
re = zeros(2*K+1,2*K+1);
im = zeros(2*K+1,2*K+1);
F = complex(re, im);

% Fourier transform
debug = 0;
for k1 = -K:K
	for k2 = -K:K
		for x = 1:N
			for y = 1:M
				F(k1+K+1,k2+K+1) = F(k1+K+1,k2+K+1) + ...
				f(x,y) * exp(-1i*pi*x*k1/N-1i*pi*y*k2/M);
			end
		end
		debug = debug + 1;
		disp(debug);
	end
end

% Fourier inverse transform
re = zeros(N,M);
im = zeros(N,M);
F_inv = complex(re,im);
debug = 0;
for x = 1:N
	for y = 1:M
		for k1 = -K:K
			for k2 = -K:K
				F_inv(x,y) = F_inv(x,y) + ...
				F(k1+K+1,k2+K+1) * exp(1i*pi*(x/N)*k1+1i*pi*(y/M)*k2);
			end
		end
		debug = debug + 1;
		disp(debug);
	end
end

A = max(max(F_inv));
B = min(min(F_inv));
imwrite(uint8((F_inv-B)/(A-B)*255),'pic2.jpg');

效果示下：
原图pic.jpg

K = 75，大概运行了三四个小时

K = 50