n = 1024;
name = ‘gaussiannoise’; % filtered gaussian noise
name = ‘pieceregular’; % piecewise regular
f = load_signal(name, n);
sigma = 0.03 * (max(f)min(f)); % noise level
fn = f + sigma*randn(1,n); % noisy signal
% plot signals
subplot(2,1,1);
plot(f); axis tight; title(‘Original’);
subplot(2,1,2);
plot(fn); axis tight; title(‘Noisy’);
% mu is width of the filter, approximat
mu = 5;
m = 101; % total size of the filter, odd number is better
h = build_gaussian_filter(m,mu/(4*n),n);
fnh = perform_convolution(fn,h);
% display
subplot(3,1,1);
plot(f); axis tight; title(‘Original’);
subplot(3,1,2);
plot(fn); axis tight; title(‘Noisy’);
subplot(3,1,3);
plot(fnh); axis tight; title(‘Estimated’);
Denoising with linear filtering.You can see that blurring tends to destroy singularities on the left. On the contrary, linear filtering is well adapted to stationary signals such as filtered noise (right). Linear denoising is bound to fail on complex signals (see the next lecture for a nonlinear denoising using wavelets).
mu_list = linspace( 0, 12, 15 );
err = [];
clf;
for i = 1:length(mu_list)
mu = mu_list(i);
h = build_gaussian_filter(m,mu/(4*n),n);
fnh = perform_convolution(fn,h);
% compute the error
e = sum( (f(:)fnh(:)).^2 );
err = [ err, e ];
end
save_image([rep ‘denoiseprogr‘ name]);
[tmp,i] = min(err); mu = mu_list(i);
% plot the denoising result here
…
Dectection of the optimal value for the variance of the filter. This value is quite low, which result in a conservative filter. This filter prefers not to remove too much noise.
ff = fft(f); % Fourier transform of the signal
ffn = fft(fn); % Fourier transform of the noisy input
pf = abs(ff).^2; % spectral power
% fourier transform of the wiener filter
hwf = pf./(pf+ n*sigma^2);
% filter in the spatial domain (just for display)
hw = fftshift( ifft(hwf) );
% perform the filtering over the fourier domain
fw = real( ifft(ffn .* hwf) );
% compute the error. Warning: remove the boundary.
sel = n/2256:n/2+256;
pwiener = psnr(f(sel),fw(sel));
pgauss = psnr(f(sel),fnh(sel));
% display
…
Comparison of Gaussian filtering vs. Wiener filtering.
Wiener is optimal among linear process.
But one could do better with nonlinear denoising for
non gaussian sources.
n = 256;
name = ‘gaussiannoise’;
name = ‘lena’;
I = load_image(name);
% reduce size to speed up
I = I(end/2n/2+1:end/2+n/2,end/2n/2+1:end/2+n/2);
% to avoid saturation
I = rescale(I,15,240);
% add noise
sigma = 0.18 * (max(I(:))min(I(:)));
In = I + sigma*randn(n);
mu = 5; m = 31;
h = build_gaussian_filter([m m],mu/(4*n),[n n]);
Inh = perform_convolution(In,h); % calcule la convolution
subplot(1,3,1);
image( I ); axis image; axis off; title(‘Original’);
subplot(1,3,2);
image( clamp(In,0,255) ); axis image; axis off; title(‘Noisy’);
subplot(1,3,3);
image( clamp(Inh,0,255) ); axis image; axis off; title(‘Smoothed’);
colormap gray(256);
Denoising with linear filtering.
Up: lena image.
Bottom: filtered gaussian noise.
mu_list = linspace( 0, 8, 12 );
err = [];
for mu = mu_list
…
end
% display the error function
plot(err, ‘.‘)
% selection of the best variance
[tmp,i] = min(err); mu = mu_list(i);
display the denoising result
…
Progression of the denoising for increasing filtering size.


Denoising result for the optimal value of the variance of the filters.
Once again, there is still a lot of noise after the blurring. 
fI = fft2(I);
fIn = fft2(In);
pf = abs(fI).^2; % spectral power
% fourier transform of the wiener filter
hwf = pf./(pf+ n^2*sigma^2);
% filter in the spatial domain
hw = fftshift( ifft2(hwf) );
% perform the filtering over the fourier
Iw = real( ifft2(fIn .* hwf) );
sel = n/2100:n/2+100;
pwiener = psnr(I(sel,sel),Iw(sel,sel));
pgauss = psnr(I(sel,sel),Inh(sel,sel));
% display
…
Comparison of best gaussian filter with Wiener filter.
As in 1D, Wiener filter is the best linear filtering denoising. Note however that nonlinear denoisings achieve much better results on natural images. 
Source:
http://www.ceremade.dauphine.fr/~peyre/teaching/wavelets/tp3.html
Virtual Fashion Education
"chúng tôi chỉ là tôi tớ của anh em, vì Đức Kitô" (2Cr 4,5b)
hienphap.net
News About Tech, Money and Innovation
Modern art using the GPU
Find the perfect theme for your blog.
Learn to Learn
Con tằm đến thác vẫn còn vương tơ
Khoa Vật lý, Đại học Sư phạm Tp.HCM  ĐT :(08)38352020  109
Blog Toán Cao Cấp (M4Ps)
Indulge Travel, Adventure, & New Experiences
"Behind every stack of books there is a flood of knowledge."
The latest news on WordPress.com and the WordPress community.