function [xs, g, rxy, ryy] = wiener_td(x, y, sigma_noise)
%rxy = ifft(fft(x).*conj(fft(y)));
x = x - mean(x);
y = y - mean(y);
rxy = circulant(y, 'step', -1).'*conj(x);
%ryy = ifft(fft(y).*conj(fft(y)));
ryy = circulant(y, 'step', -1).'*conj(y);

sum_rxy = sum(abs(rxy))
sum_ryy = sum(abs(ryy))

R = toeplitz(ryy + sigma_noise.^2*(1-j));
g = inv(R)*(rxy)./sum_ryy;

sum_g = mean(abs(g))
g = g./sum_g;

xs = (circulant((g)))*(y);