% Filtering noisy signal

np = 2^12;
a = 0;
b = 4*pi;
x = linspace(a, b, np);
c = sin(x);               % clear signal

%plot(x,c)
%xlim([a b])

cf = myradix2dft(c);
cf2 = conj(cf).*cf;

%plot(1:np,cf2)
%xlim([-100 np+100])

noise = std(c)/2;
r = noise*randn(np,1)';   % noise
s = c + r;                % noisy signal

%plot(x, s)
%xlim([a b])

sf = myradix2dft(s);
sf2 = conj(sf).*sf;

%plot(1:np,sf2)
%xlim([-100 np+100])
%ylim([0 10000])

phi = max(sf2-10000*noise,0) ./ sf2;   % filter
sout = real(myradix2invdft(phi .* sf));          % filtered signal

figure(1)
plot(x, c)
title('Original signal')
xlim([a b])

figure(2)
plot(x, sout)
title('Filtered signal')
xlim([a b])

figure(3)
plot(x, s)
hold on
plot(x, sout, 'r')
title('Noisy and filtered signals')
xlim([a b])

